ENH: Adding Discrete Approximation of VAR Methods

quantecon/markov/approximation.py

@@ -291,6 +291,8 @@ def discrete_var(A,
         An instance of `scipy.stats` for a distribution that :math:`{u_t}`
         is drawn. It must have a zero mean and unit standard deviation. 
         If None, then standard normal distribution is used.
+    order : str, optional(default='C')
+            ('C' or 'F') order in which the cartesian product is enumerated 
     random_state: a `np.random.RandomState` or `np.random.Generator` instance, or None, optional(default=None)
         If None, the `np.random.RandomState` singleton is returned.
 
@@ -332,7 +334,7 @@ def discrete_var(A,
     if rv is None:
         u = random_state.standard_normal(size=(sim_length-1, r))
     else:
-        u = rv.rvs(size=sim_length-1, random_state=random_state)
+        u = rv.rvs(size=(sim_length-1, r), random_state=random_state)
 
     v = C @ u.T
     x0 = np.zeros(m)
@@ -352,8 +354,7 @@ def discrete_var(A,
 
     V = [np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i])
          for i in range(m)]
-
     # Fit the Markov chain
     mc = fit_discrete_mc(X.T, V, order=order)
-
+    print(mc)
     return mc

quantecon/markov/tests/test_approximation.py

@@ -128,18 +128,41 @@ def setup_method(self):
                      [ 0.38556417,  0.        ],
                      [ 0.38556417,  0.05387746]]
 
-        self.P_out = [[1.61290323e-02, 1.12903226e-01, 0.00000000e+00, 
-        3.70967742e-01, 5.00000000e-01, 0.00000000e+00],
-        [1.00964548e-04, 8.51048124e-01, 3.82857566e-02, 4.03858192e-04,
-         1.10111936e-01, 4.93604457e-05],
-        [0.00000000e+00, 3.02295449e-01, 6.97266822e-01, 0.00000000e+00,
-         3.85201268e-04, 5.25274456e-05],
-        [3.60600761e-05, 4.86811027e-04, 0.00000000e+00, 6.97473992e-01,
-         3.02003137e-01, 0.00000000e+00],
-        [3.17039037e-05, 1.11090478e-01, 4.55177474e-04, 3.75374219e-02,
-         8.50778784e-01, 1.06434534e-04],
-        [0.00000000e+00, 4.45945946e-01, 3.37837838e-01, 0.00000000e+00,
-         1.89189189e-01, 2.70270270e-02]]
+        self.S_out_orderF =[
+                  [-0.38556417, -0.05387746],
+                  [ 0.38556417, -0.05387746],
+                  [-0.38556417,  0.        ],
+                  [ 0.38556417,  0.        ],
+                  [-0.38556417,  0.05387746],
+                  [ 0.38556417,  0.05387746]]
+
+        self.P_out = [
+            [1.61290323e-02, 1.12903226e-01, 0.00000000e+00, 3.70967742e-01, 
+            5.00000000e-01, 0.00000000e+00],
+            [1.00964548e-04, 8.51048124e-01, 3.82857566e-02, 4.03858192e-04,
+            1.10111936e-01, 4.93604457e-05],
+            [0.00000000e+00, 3.02295449e-01, 6.97266822e-01, 0.00000000e+00,
+            3.85201268e-04, 5.25274456e-05],
+            [3.60600761e-05, 4.86811027e-04, 0.00000000e+00, 6.97473992e-01,
+            3.02003137e-01, 0.00000000e+00],
+            [3.17039037e-05, 1.11090478e-01, 4.55177474e-04, 3.75374219e-02,
+            8.50778784e-01, 1.06434534e-04],
+            [0.00000000e+00, 4.45945946e-01, 3.37837838e-01, 0.00000000e+00,
+            1.89189189e-01, 2.70270270e-02]]
+
+        self.P_out_orderF =[
+            [1.61290323e-02, 3.70967742e-01, 1.12903226e-01, 5.00000000e-01, 
+            0.00000000e+00, 0.00000000e+00],
+            [3.60600761e-05, 6.97473992e-01, 4.86811027e-04, 3.02003137e-01, 
+            0.00000000e+00, 0.00000000e+00],
+            [1.00964548e-04, 4.03858192e-04, 8.51048124e-01, 1.10111936e-01, 
+            3.82857566e-02, 4.93604457e-05],
+            [3.17039037e-05, 3.75374219e-02, 1.11090478e-01, 8.50778784e-01, 
+            4.55177474e-04, 1.06434534e-04],
+            [0.00000000e+00, 0.00000000e+00, 3.02295449e-01, 3.85201268e-04, 
+            6.97266822e-01, 5.25274456e-05],
+            [0.00000000e+00, 0.00000000e+00, 4.45945946e-01, 1.89189189e-01, 
+            3.37837838e-01, 2.70270270e-02]]
 
         self.A, self.C, self.S_out, self.P_out = map(np.array,
                                 (self.A, self.C, self.S_out, self.P_out))
@@ -152,4 +175,24 @@ def test_discretization(self):
                 self.A, self.C, grid_sizes=self.sizes,
                 sim_length=self.T, random_state=self.random_state)
         assert_allclose(mc.state_values, self.S_out)
-        assert_allclose(mc.P, self.P_out)
+        assert_allclose(mc.P, self.P_out)
+
+    def test_sp_distributions(self):
+        mc = discrete_var(
+                self.A, self.C, 
+                grid_sizes=self.sizes,
+                sim_length=self.T, 
+                rv=sp.stats.multivariate_normal,
+                random_state=self.random_state)
+        assert_allclose(mc.state_values, self.S_out)
+        assert_allclose(mc.P, self.P_out)
+
+    def test_order_F(self):
+        mc = discrete_var(
+                self.A, self.C, 
+                grid_sizes=self.sizes,
+                sim_length=self.T, 
+                order='F',
+                random_state=self.random_state)
+        assert_allclose(mc.state_values, self.S_out_orderF)
+        assert_allclose(mc.P, self.P_out_orderF)