Add rouwenhorst method for approx AR(1) with MC

quantecon/__init__.py

@@ -30,7 +30,7 @@
 from .matrix_eqn import solve_discrete_lyapunov, solve_discrete_riccati
 from .quadsums import var_quadratic_sum, m_quadratic_sum
 #->Propose Delete From Top Level
-from .markov import MarkovChain, random_markov_chain, random_stochastic_matrix, gth_solve, tauchen 	 	#Promote to keep current examples working
+from .markov import MarkovChain, random_markov_chain, random_stochastic_matrix, gth_solve, tauchen, rouwenhorst 	 	#Promote to keep current examples working
 from .markov import mc_compute_stationary, mc_sample_path 												#Imports that Should be Deprecated with markov package
 #<-
 from .rank_nullspace import rank_est, nullspace

quantecon/markov/__init__.py

@@ -7,6 +7,6 @@
 from .gth_solve import gth_solve
 from .random import random_markov_chain, random_stochastic_matrix, \
     random_discrete_dp
-from .approximation import tauchen
+from .approximation import tauchen, rouwenhorst
 from .ddp import DiscreteDP, backward_induction
 from .utilities import sa_indices

quantecon/markov/approximation.py

@@ -16,6 +16,103 @@
 from numba import njit
 
 
+def rouwenhorst(n, ybar, sigma, rho):
+    """
+    Takes as inputs n, p, q, psi.  It will then construct a markov chain
+    that estimates an AR(1) process of:
+    y_t = \bar{y} + \rho y_{t-1} + \varepsilon_t
+    where \varepsilon_t is i.i.d. normal of mean 0, std dev of sigma
+
+    The Rouwenhorst approximation uses the following recursive defintion
+    for approximating a distribution:
+
+    theta_2 = [p    , 1 - p]
+              [1 - q, q    ]
+
+    theta_{n+1} = p [theta_n, 0] + (1 - p) [0, theta_n]
+                    [0      , 0]           [0,       0]
+                + q [0       , 0] + (1 - q) [0,        ]
+                    [theta_n , 0]           [0, theta_n]
+
+    Parameters
+    ----------
+    n : int
+        The number of points to approximate the distribution
+
+    ybar : float
+        The value \bar{y} in the process.  Note that the mean of this
+        AR(1) process, y, is simply ybar/(1 - rho)
+
+    sigma : float
+        The value of the standard deviation of the \varepsilon process
+
+    rho : float
+        By default this will be 0, but if you are approximating an AR(1)
+        process then this is the autocorrelation across periods
+
+    Returns
+    -------
+
+    mc : MarkovChain
+        An instance of the MarkovChain class that stores the transition 
+        matrix and state values returned by the discretization method
+
+    """
+
+    # Get the standard deviation of y
+    y_sd = sqrt(sigma**2 / (1 - rho**2))
+
+    # Given the moments of our process we can find the right values
+    # for p, q, psi because there are analytical solutions as shown in
+    # Gianluca Violante's notes on computational methods
+    p = (1 + rho) / 2
+    q = p
+    psi = y_sd * np.sqrt(n - 1)
+
+    # Find the states
+    ubar = psi
+    lbar = -ubar
+
+    bar = np.linspace(lbar, ubar, n)
+
+    def row_build_mat(n, p, q):
+        """
+        This method uses the values of p and q to build the transition
+        matrix for the rouwenhorst method
+        """
+
+        if n == 2:
+            theta = np.array([[p, 1 - p], [1 - q, q]])
+
+        elif n > 2:
+            p1 = np.zeros((n, n))
+            p2 = np.zeros((n, n))
+            p3 = np.zeros((n, n))
+            p4 = np.zeros((n, n))
+
+            new_mat = row_build_mat(n - 1, p, q)
+
+            p1[:n - 1, :n - 1] = p * new_mat
+            p2[:n - 1, 1:] = (1 - p) * new_mat
+            p3[1:, :-1] = (1 - q) * new_mat
+            p4[1:, 1:] = q * new_mat
+
+            theta = p1 + p2 + p3 + p4
+            theta[1:n - 1, :] = theta[1:n - 1, :] / 2
+
+        else:
+            raise ValueError("The number of states must be positive " +
+                             "and greater than or equal to 2")
+
+        return theta
+
+    theta = row_build_mat(n, p, q)
+
+    bar += ybar / (1 - rho)
+
+    return MarkovChain(theta, bar)
+
+
 def tauchen(rho, sigma_u, m=3, n=7):
     """
     Computes a Markov chain associated with a discretized version of

quantecon/markov/tests/test_approximation.py

@@ -11,7 +11,8 @@
 import unittest
 import numpy as np
 from numpy.testing import assert_allclose
-from quantecon.markov import tauchen
+from quantecon.markov import tauchen, rouwenhorst
+#from quantecon.markov.approximation import rouwenhorst
 
 
 class TestTauchen(unittest.TestCase):
@@ -50,7 +51,52 @@ def test_states_sum_0(self):
         self.assertTrue(abs(np.sum(self.x)) < self.tol)
 
 
+class TestRouwenhorst(unittest.TestCase):
+
+    def setUp(self):
+        self.rho, self.sigma = np.random.uniform(0,1, size=2)
+        self.n = np.random.random_integers(3,25)
+        self.ybar = np.random.random_integers(0,10)
+        self.tol = 1e-12
+
+        mc = rouwenhorst(self.n, self.ybar, self.sigma, self.rho)
+        self.x, self.P = mc.state_values, mc.P
+
+    def tearDown(self):
+        del self.x
+        del self.P
+
+    def testShape(self):
+        i, j = self.P.shape
+
+        self.assertTrue(i == j)
+
+    def testDim(self):
+        dim_x = self.x.ndim
+        dim_P = self.P.ndim
+        self.assertTrue(dim_x == 1 and dim_P == 2)
+
+    def test_transition_mat_row_sum_1(self):
+        self.assertTrue(np.allclose(np.sum(self.P, axis=1), 1, atol=self.tol))
+
+    def test_positive_probs(self):
+        self.assertTrue(np.all(self.P) > -self.tol)
+
+    def test_states_sum_0(self):
+        tol = self.tol + self.n*(self.ybar/(1 - self.rho))
+        self.assertTrue(abs(np.sum(self.x)) < tol)
+
+    def test_control_case(self):
+        n = 3; ybar = 1; sigma = 0.5; rho = 0.8;
+        mc_rouwenhorst = rouwenhorst(n, ybar, sigma, rho)
+        mc_rouwenhorst.x, mc_rouwenhorst.P = mc_rouwenhorst.state_values, mc_rouwenhorst.P 
+        sigma_y = np.sqrt(sigma**2 / (1-rho**2))
+        psi = sigma_y * np.sqrt(n-1)
+        known_x = np.array([-psi+5.0, 5., psi+5.0])
+        known_P = np.array([[0.81, 0.18, 0.01], [0.09, 0.82, 0.09], [0.01, 0.18, 0.81]])
+        self.assertTrue(sum(mc_rouwenhorst.x - known_x) < self.tol and sum(sum(mc_rouwenhorst.P - known_P)) < self.tol)
+
 
 if __name__ == '__main__':
-    suite = unittest.TestLoader().loadTestsFromTestCase(TestTauchen)
+    suite = unittest.TestLoader().loadTestsFromTestCase([TestTauchen, TestRouwenhorst])
     unittest.TextTestRunner(verbosity=2, stream=sys.stderr).run(suite)