ENH: Adding Discrete Approximation of VAR Methods

quantecon/markov/__init__.py

@@ -8,7 +8,7 @@
 from .gth_solve import gth_solve
 from .random import random_markov_chain, random_stochastic_matrix, \
     random_discrete_dp
-from .approximation import tauchen, rouwenhorst
+from .approximation import tauchen, rouwenhorst, discrete_var
 from .ddp import DiscreteDP, backward_induction
 from .utilities import sa_indices
 from .estimate import estimate_mc, fit_discrete_mc

quantecon/markov/approximation.py

@@ -7,6 +7,10 @@
 
 from math import erfc, sqrt
 from .core import MarkovChain
+from .estimate import fit_discrete_mc
+from ..lss import simulate_linear_model
+from ..matrix_eqn import solve_discrete_lyapunov
+from ..util import check_random_state
 
 import warnings
 import numpy as np
@@ -242,3 +246,115 @@ def _fill_tauchen(x, P, n, rho, sigma, half_step):
             z = x[j] - rho * x[i]
             P[i, j] = (std_norm_cdf((z + half_step) / sigma) -
                        std_norm_cdf((z - half_step) / sigma))
+
+
+def discrete_var(A,
+                 C,
+                 grid_sizes=None,
+                 std_devs=np.sqrt(10),
+                 sim_length=1_000_000,
+                 rv=None,
+                 order='C',
+                 random_state=None):
+    r"""
+    This code discretizes a VAR(1) process of the form:
+    .. math::
+        x_t = A x_{t-1} + C u_t
+    where 
+    :math:`{u_t}` is drawn iid from a distribution with mean 0 
+    and unit standard deviaiton;
+    :math:`{C}` is a volatility matrix in linear state space model.
+    By default, the code removes the states that are never visited under the
+    simulation that computes the transition probabilities.
+    For a mathematical derivation check *Finite-State Approximation Of
+    VAR Processes:  A Simulation Approach* by Stephanie Schmitt-Grohé and
+    Martín Uribe, July 11, 2010.
+    This code was adapted from Schmitt-Grohé and Uribe's original MATLAB code.
+
+    Parameters
+    ----------
+    A : array_like(float)
+        An m x m matrix containing the process' autocorrelation parameters
+    C : array_like(float)
+        An m x m volatility matrix
+    grid_sizes : array_like(int) or None, optional(default=None)
+        An m-vector containing the number of grid points in the discretization
+        of each dimension of x_t. If None, then grid_sizes is
+        set to (10, ..., 10).
+    std_devs : float, optional(default=np.sqrt(10))
+        The number of standard deviations the grid should stretch in each
+        dimension, where standard deviations are measured under the stationary
+        distribution.
+    sim_length : int, optional(default=1_000_000)
+        The length of the simulated time series.
+    rv: a `scipy.stats` instance or None, optional(default=None)
+        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.
+
+    Returns
+    -------
+    mc : MarkovChain
+        An instance of the MarkovChain class that stores the transition
+        matrix and state values returned by the discretization method.
+        The MarkovChain instance contains:
+        P : A square matrix containing the transition probability
+            matrix of the discretized state.
+        S : An array where element (i,j) of S is the discretized
+            value of the j-th element of x_t in state i. Reducing S to its
+            unique values yields the grid values. The cartesian product
+            state grid uses a row major ordering.
+
+    Example
+    -------
+        This example discretizes the stochastic process used to calibrate
+        the economic model included in ``Downward Nominal Wage Rigidity,
+        Currency Pegs, and Involuntary Unemployment'' by Stephanie
+        Schmitt-Grohé and Martín Uribe, Journal of Political Economy 124,
+        October 2016, 1466-1514.
+            A     = np.array([[0.7901, -1.3570],
+                              [-0.0104, 0.8638]])
+            Omega = np.array([[0.0012346, -0.0000776],
+                              [-0.0000776, 0.0000401]])
+            C = sp.linalg.sqrtm(Omega)
+            grid_sizes = np.array([21, 11])
+            mc = discrete_var(A, C, grid_sizes, sim_length=1_000_000)
+    """
+    A = np.asarray(A)
+    C = np.asarray(C)
+    m, r = C.shape
+
+    # Run simulation to compute transition probabilities
+    random_state = check_random_state(random_state)
+
+    if rv is None:
+        u = random_state.standard_normal(size=(sim_length-1, r))
+    else:
+        u = rv.rvs(size=(sim_length-1, r), random_state=random_state)
+
+    v = C @ u.T
+    x0 = np.zeros(m)
+    X = simulate_linear_model(A, x0, v, ts_length=sim_length)
+
+    # Compute stationary variance-covariance matrix of AR process and use
+    # it to obtain grid bounds.
+    Sigma = solve_discrete_lyapunov(A, C @ C.T)
+    sigma_vector = np.sqrt(np.diagonal(Sigma))    # Stationary std dev
+    upper_bounds = std_devs * sigma_vector
+
+    # Build the individual grids along each dimension
+    if grid_sizes is None:
+        # Set the size of every grid to default_grid_size
+        default_grid_size = 10
+        grid_sizes = np.full(m, default_grid_size)
+
+    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

@@ -4,11 +4,9 @@
 """
 import numpy as np
 import pytest
-from quantecon.markov import tauchen, rouwenhorst
+from quantecon.markov import tauchen, rouwenhorst, discrete_var
 from numpy.testing import assert_, assert_allclose, assert_raises
-
-#from quantecon.markov.approximation import rouwenhorst
-
+import scipy as sp
 
 class TestTauchen:
 
@@ -107,3 +105,94 @@ def test_control_case(self):
             [[0.81, 0.18, 0.01], [0.09, 0.82, 0.09], [0.01, 0.18, 0.81]])
         assert_(np.sum(mc_rouwenhorst.x - known_x) < self.tol and
                 np.sum(mc_rouwenhorst.P - known_P) < self.tol)
+
+
+class TestDiscreteVar:
+
+    def setup_method(self):
+        self.random_state = np.random.RandomState(0)
+
+        Omega = np.array([[0.0012346, -0.0000776],
+                              [-0.0000776, 0.0000401]])
+        self.A = [[0.7901, -1.3570], 
+             [-0.0104, 0.8638]]
+        self.C = sp.linalg.sqrtm(Omega)
+        self.T= 1_000_000
+        self.sizes = np.array((2, 3))
+
+        # Expected outputs
+        self.S_out = [[-0.38556417, -0.05387746],
+                     [-0.38556417,  0.        ],
+                     [-0.38556417,  0.05387746],
+                     [ 0.38556417, -0.05387746],
+                     [ 0.38556417,  0.        ],
+                     [ 0.38556417,  0.05387746]]
+
+        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))
+
+    def teardown_method(self):
+        del self.A, self.C, self.S_out, self.P_out
+
+    def test_discretization(self):
+        mc = discrete_var(
+                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)
+
+    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)