Add method option in robustlq.py

Author: unknown

quantecon/robustlq.py

@@ -152,7 +152,7 @@ def b_operator(self, P):
 
         return F, new_P
 
-    def robust_rule(self):
+    def robust_rule(self, method='doubling'):
         """
         This method solves the robust control problem by tricking it
         into a stacked LQ problem, as described in chapter 2 of Hansen-
@@ -165,6 +165,12 @@ def robust_rule(self):
 
         And the value function is :math:`-x'Px`
 
+        Parameters
+        ----------
+        method : str, optional(default='doubling')
+            Solution method used in solving the associated Riccati
+            equation, str in {'doubling', 'qz'}.
+            
         Returns
         -------
         F : array_like(float, ndim=2)
@@ -189,7 +195,7 @@ def robust_rule(self):
             lq = LQ(-beta*I*theta, R, A, C, beta=beta)
 
             # == Solve and convert back to robust problem == #
-            P, f, d = lq.stationary_values()
+            P, f, d = lq.stationary_values(method=method)
             F = np.zeros((self.k, self.n))
             K = -f[:k, :]
 
@@ -199,7 +205,7 @@ def robust_rule(self):
             lq = LQ(Qa, R, A, Ba, beta=beta)
 
             # == Solve and convert back to robust problem == #
-            P, f, d = lq.stationary_values()
+            P, f, d = lq.stationary_values(method=method)
             F = f[:k, :]
             K = -f[k:f.shape[0], :]
 
@@ -254,14 +260,17 @@ def robust_rule_simple(self, P_init=None, max_iter=80, tol=1e-8):
 
         return F, K, P
 
-    def F_to_K(self, F):
+    def F_to_K(self, F, method='doubling'):
         """
         Compute agent 2's best cost-minimizing response K, given F.
 
         Parameters
         ----------
         F : array_like(float, ndim=2)
             A k x n array
+        method : str, optional(default='doubling')
+            Solution method used in solving the associated Riccati
+            equation, str in {'doubling', 'qz'}.
 
         Returns
         -------
@@ -276,18 +285,21 @@ def F_to_K(self, F):
         A2 = self.A - dot(self.B, F)
         B2 = self.C
         lq = LQ(Q2, R2, A2, B2, beta=self.beta)
-        neg_P, neg_K, d = lq.stationary_values()
+        neg_P, neg_K, d = lq.stationary_values(method=method)
 
         return -neg_K, -neg_P
 
-    def K_to_F(self, K):
+    def K_to_F(self, K, method='doubling'):
         """
         Compute agent 1's best value-maximizing response F, given K.
 
         Parameters
         ----------
         K : array_like(float, ndim=2)
             A j x n array
+        method : str, optional(default='doubling')
+            Solution method used in solving the associated Riccati
+            equation, str in {'doubling', 'qz'}.
 
         Returns
         -------
@@ -302,7 +314,7 @@ def K_to_F(self, K):
         Q1 = self.Q
         R1 = self.R - self.beta * self.theta * dot(K.T, K)
         lq = LQ(Q1, R1, A1, B1, beta=self.beta)
-        P, F, d = lq.stationary_values()
+        P, F, d = lq.stationary_values(method=method)
 
         return F, P
 
@@ -388,3 +400,4 @@ def evaluate_F(self, F):
         o_F = (ho + beta * tr) / (1 - beta)
 
         return K_F, P_F, d_F, O_F, o_F
+        

quantecon/tests/test_robustlq.py

@@ -51,9 +51,7 @@ def setUp(self):
         self.rblq_test = RBLQ(Q, R, A, B, C, beta, theta)
         self.rblq_test_pf = RBLQ(Q_pf, R, A, B_pf, C, beta, theta)
         self.lq_test = LQ(Q, R, A, B, C, beta)
-
-        self.Fr, self.Kr, self.Pr = self.rblq_test.robust_rule()
-        self.Fr_pf, self.Kr_pf, self.Pr_pf = self.rblq_test_pf.robust_rule()
+        self.methods = ['doubling', 'qz']
 
     def tearDown(self):
         del self.rblq_test
@@ -65,52 +63,48 @@ def test_pure_forecasting(self):
     def test_robust_rule_vs_simple(self):
         rblq = self.rblq_test
         rblq_pf = self.rblq_test_pf
-        Fr, Kr, Pr = self.Fr, self.Kr, self.Pr
-        Fr_pf, Kr_pf, Pr_pf = self.Fr_pf, self.Kr_pf, self.Pr_pf
 
-        Fs, Ks, Ps = rblq.robust_rule_simple(P_init=Pr, tol=1e-12)
-        Fs_pf, Ks_pf, Ps_pf = rblq_pf.robust_rule_simple(P_init=Pr_pf, tol=1e-12)
+        for method in self.methods:
+            Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)
+            Fr_pf, Kr_pf, Pr_pf = self.rblq_test_pf.robust_rule(method=method)
+
+            Fs, Ks, Ps = rblq.robust_rule_simple(P_init=Pr, tol=1e-12)
+            Fs_pf, Ks_pf, Ps_pf = rblq_pf.robust_rule_simple(P_init=Pr_pf, tol=1e-12)
 
-        assert_allclose(Fr, Fs, rtol=1e-4)
-        assert_allclose(Kr, Ks, rtol=1e-4)
-        assert_allclose(Pr, Ps, rtol=1e-4)
+            assert_allclose(Fr, Fs, rtol=1e-4)
+            assert_allclose(Kr, Ks, rtol=1e-4)
+            assert_allclose(Pr, Ps, rtol=1e-4)
 
-        atol = 1e-10
-        assert_allclose(Fr_pf, Fs_pf, rtol=1e-4)
-        assert_allclose(Kr_pf, Ks_pf, rtol=1e-4, atol=atol)
-        assert_allclose(Pr_pf, Ps_pf, rtol=1e-4, atol=atol)
+            atol = 1e-10
+            assert_allclose(Fr_pf, Fs_pf, rtol=1e-4)
+            assert_allclose(Kr_pf, Ks_pf, rtol=1e-4, atol=atol)
+            assert_allclose(Pr_pf, Ps_pf, rtol=1e-4, atol=atol)
 
 
     def test_f2k_and_k2f(self):
         rblq = self.rblq_test
-        Fr, Kr, Pr = self.Fr, self.Kr, self.Pr
-
-        K_f2k, P_f2k = rblq.F_to_K(Fr)
-        F_k2f, P_k2f = rblq.K_to_F(Kr)
 
-        assert_allclose(K_f2k, Kr, rtol=1e-4)
-        assert_allclose(F_k2f, Fr, rtol=1e-4)
-        assert_allclose(P_f2k, P_k2f, rtol=1e-4)
+        for method in self.methods:
+            Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)
+            K_f2k, P_f2k = rblq.F_to_K(Fr, method=method)
+            F_k2f, P_k2f = rblq.K_to_F(Kr, method=method)
+            assert_allclose(K_f2k, Kr, rtol=1e-4)
+            assert_allclose(F_k2f, Fr, rtol=1e-4)
+            assert_allclose(P_f2k, P_k2f, rtol=1e-4)
 
     def test_evaluate_F(self):
         rblq = self.rblq_test
-        Fr, Kr, Pr = self.Fr, self.Kr, self.Pr
-
-        Kf, Pf, df, Of, of = rblq.evaluate_F(Fr)
-
-        # In the future if we wanted, we could check more things, but I
-        # think the other pieces are basically just plugging these into
-        # equations so if these hold then the others should be correct
-        # as well.
-        assert_allclose(Pf, Pr)
-        assert_allclose(Kf, Kr)
-
-
-
-
-
+        for method in self.methods:
+            Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)
 
+            Kf, Pf, df, Of, of = rblq.evaluate_F(Fr)
 
+            # In the future if we wanted, we could check more things, but I
+            # think the other pieces are basically just plugging these into
+            # equations so if these hold then the others should be correct
+            # as well.
+            assert_allclose(Pf, Pr)
+            assert_allclose(Kf, Kr)
 
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestRBLQControl)