From 5825dc9229f219f6f025596619dd35df7e5be8c6 Mon Sep 17 00:00:00 2001 From: Spencer Lyon Date: Mon, 20 Feb 2017 20:18:37 -0500 Subject: [PATCH 1/2] Add rouwenhorst method for approx AR(1) with MC NOTE: This needs tests... --- quantecon/markov/approximation.py | 97 +++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index b7a475631..cf842ab87 100644 --- a/quantecon/markov/approximation.py +++ b/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 From 4ebb7f5f04d2eeef09d5e6c0502e33a33fb6a26b Mon Sep 17 00:00:00 2001 From: lbui01 Date: Wed, 1 Mar 2017 13:20:23 +1100 Subject: [PATCH 2/2] Adding test for rouwenhorst --- quantecon/__init__.py | 2 +- quantecon/markov/__init__.py | 2 +- quantecon/markov/tests/test_approximation.py | 50 +++++++++++++++++++- 3 files changed, 50 insertions(+), 4 deletions(-) diff --git a/quantecon/__init__.py b/quantecon/__init__.py index 28ec60f13..9350f8897 100644 --- a/quantecon/__init__.py +++ b/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 diff --git a/quantecon/markov/__init__.py b/quantecon/markov/__init__.py index 22c8e9f10..4eb8aeb80 100644 --- a/quantecon/markov/__init__.py +++ b/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 diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 92d44d576..dbb4c0f36 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/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)