Skip to content

Commit

Permalink
Unify tauchen and rouwenhorst API (#664)
Browse files Browse the repository at this point in the history 
* [Breaking Change] Modify `tauchen` API to be same as julia

* Fix tests

* Fix docs

* [Breaking Change] Modify rouwenhorst API to be same as julia

* Fix tests

* Add a note on p & q values

* Add warnings

* Test warnings
  • Loading branch information
@Smit-create
Smit-create authored Dec 17, 2022
1 parent b0c4bcf commit 0e21854
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 37 deletions.
74 changes: 50 additions & 24 deletions quantecon/markov/approximation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,16 @@
from math import erfc, sqrt
from .core import MarkovChain

import warnings
import numpy as np
from numba import njit


def rouwenhorst(n, ybar, sigma, rho):
def rouwenhorst(n, rho, sigma, mu=0.):
r"""
Takes as inputs n, p, q, psi. It will then construct a markov chain
Takes as inputs n, mu, sigma, rho. It will then construct a markov chain
that estimates an AR(1) process of:
:math:`y_t = \bar{y} + \rho y_{t-1} + \varepsilon_t`
:math:`y_t = \mu + \rho y_{t-1} + \varepsilon_t`
where :math:`\varepsilon_t` is i.i.d. normal of mean 0, std dev of sigma
The Rouwenhorst approximation uses the following recursive defintion
Expand Down Expand Up @@ -54,22 +55,23 @@ def rouwenhorst(n, ybar, sigma, rho):
0 & \theta_n \\
\end{bmatrix}
where :math:`{p = q = \frac{(1 + \rho)}{2}}`
Parameters
----------
n : int
The number of points to approximate the distribution
ybar : float
The value :math:`\bar{y}` in the process. Note that the mean of this
AR(1) process, :math:`y`, is simply :math:`\bar{y}/(1 - \rho)`
rho : float
Persistence parameter in AR(1) process, if you are approximating
an AR(1) process then this is the autocorrelation across periods.
sigma : float
The value of the standard deviation of the :math:`\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
mu : float, optional(default=0.0)
The value :math:`\mu` in the process. Note that the mean of this
AR(1) process, :math:`y`, is simply :math:`\mu/(1 - \rho)`
Returns
-------
Expand All @@ -78,8 +80,19 @@ def rouwenhorst(n, ybar, sigma, rho):
An instance of the MarkovChain class that stores the transition
matrix and state values returned by the discretization method
Note
----
UserWarning: The API of `rouwenhorst` was changed from
`rouwenhorst(n, ybar, sigma, rho)` to
`rouwenhorst(n, rho, sigma, mu=0.)` in version 0.6.0.
"""

warnings.warn("The API of rouwenhorst has changed from `rouwenhorst(n, ybar, sigma, rho)`"
" to `rouwenhorst(n, rho, sigma, mu=0.)`. To find more details please visit:"
" https://github.com/QuantEcon/QuantEcon.py/issues/663.",
UserWarning, stacklevel=2)
# Get the standard deviation of y
y_sd = sqrt(sigma**2 / (1 - rho**2))

Expand Down Expand Up @@ -130,35 +143,37 @@ def row_build_mat(n, p, q):

theta = row_build_mat(n, p, q)

bar += ybar / (1 - rho)
bar += mu / (1 - rho)

return MarkovChain(theta, bar)


def tauchen(rho, sigma_u, b=0., m=3, n=7):
def tauchen(n, rho, sigma, mu=0., n_std=3):
r"""
Computes a Markov chain associated with a discretized version of
the linear Gaussian AR(1) process
.. math::
y_{t+1} = b + \rho y_t + u_{t+1}
y_t = \mu + \rho y_{t-1} + \epsilon_t
using Tauchen's method. Here :math:`{u_t}` is an i.i.d. Gaussian process
using Tauchen's method. Here :math:`{\epsilon_t}` is an i.i.d. Gaussian process
with zero mean.
Parameters
----------
b : scalar(float)
The constant term of {y_t}
n : scalar(int)
The number of states to use in the approximation
rho : scalar(float)
The autocorrelation coefficient
sigma_u : scalar(float)
The autocorrelation coefficient, Persistence parameter in AR(1) process
sigma : scalar(float)
The standard deviation of the random process
m : scalar(int), optional(default=3)
mu : scalar(float), optional(default=0.0)
The value :math:`\mu` in the process. Note that the mean of this
AR(1) process, :math:`y`, is simply :math:`\mu/(1 - \rho)`
n_std : scalar(int), optional(default=3)
The number of standard deviations to approximate out to
n : scalar(int), optional(default=7)
The number of states to use in the approximation
Returns
-------
Expand All @@ -167,13 +182,24 @@ def tauchen(rho, sigma_u, b=0., m=3, n=7):
An instance of the MarkovChain class that stores the transition
matrix and state values returned by the discretization method
Note
----
UserWarning: The API of `tauchen` was changed from
`tauchen(rho, sigma_u, b=0., m=3, n=7)` to
`tauchen(n, rho, sigma, mu=0., n_std=3)` in version 0.6.0.
"""
warnings.warn("The API of tauchen has changed from `tauchen(rho, sigma_u, b=0., m=3, n=7)`"
" to `tauchen(n, rho, sigma, mu=0., n_std=3)`. To find more details please visit:"
" https://github.com/QuantEcon/QuantEcon.py/issues/663.",
UserWarning, stacklevel=2)

# standard deviation of demeaned y_t
std_y = np.sqrt(sigma_u**2 / (1 - rho**2))
std_y = np.sqrt(sigma**2 / (1 - rho**2))

# top of discrete state space for demeaned y_t
x_max = m * std_y
x_max = n_std * std_y

# bottom of discrete state space for demeaned y_t
x_min = -x_max
Expand All @@ -187,10 +213,10 @@ def tauchen(rho, sigma_u, b=0., m=3, n=7):

# approximate Markov transition matrix for
# demeaned y_t
_fill_tauchen(x, P, n, rho, sigma_u, half_step)
_fill_tauchen(x, P, n, rho, sigma, half_step)

# shifts the state values by the long run mean of y_t
mu = b / (1 - rho)
mu = mu / (1 - rho)

mc = MarkovChain(P, state_values=x+mu)

Expand Down
31 changes: 18 additions & 13 deletions quantecon/markov/tests/test_approximation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""
import numpy as np
import pytest
from quantecon.markov import tauchen, rouwenhorst
from numpy.testing import assert_, assert_allclose

Expand All @@ -12,24 +13,26 @@
class TestTauchen:

def setup_method(self):
self.rho, self.sigma_u = np.random.rand(2)
self.rho, self.sigma = np.random.rand(2)
self.n = np.random.randint(3, 25)
self.m = np.random.randint(5)
self.n_std = np.random.randint(5)
self.tol = 1e-12
self.b = 0.
self.mu = 0.

mc = tauchen(self.rho, self.sigma_u, self.b, self.m, self.n)
with pytest.warns(UserWarning):
mc = tauchen(self.n, self.rho, self.sigma, self.mu, self.n_std)
self.x, self.P = mc.state_values, mc.P

def teardown_method(self):
del self.x
del self.P

def testStateCenter(self):
for b in [0., 1., -1.]:
mu = b / (1 - self.rho)
mc = tauchen(self.rho, self.sigma_u, b, self.m, self.n)
assert_allclose(mu, np.mean(mc.state_values), atol=self.tol)
for mu in [0., 1., -1.]:
mu_expect = mu / (1 - self.rho)
with pytest.warns(UserWarning):
mc = tauchen(self.n, self.rho, self.sigma, mu, self.n_std)
assert_allclose(mu_expect, np.mean(mc.state_values), atol=self.tol)

def testShape(self):
i, j = self.P.shape
Expand Down Expand Up @@ -57,10 +60,11 @@ class TestRouwenhorst:
def setup_method(self):
self.rho, self.sigma = np.random.uniform(0, 1, size=2)
self.n = np.random.randint(3, 26)
self.ybar = np.random.randint(0, 11)
self.mu = np.random.randint(0, 11)
self.tol = 1e-10

mc = rouwenhorst(self.n, self.ybar, self.sigma, self.rho)
with pytest.warns(UserWarning):
mc = rouwenhorst(self.n, self.rho, self.sigma, self.mu)
self.x, self.P = mc.state_values, mc.P

def teardown_method(self):
Expand All @@ -84,12 +88,13 @@ def test_positive_probs(self):
assert_(np.all(self.P > -self.tol))

def test_states_sum_0(self):
tol = self.tol + self.n*(self.ybar/(1 - self.rho))
tol = self.tol + self.n*(self.mu/(1 - self.rho))
assert_(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)
n = 3; mu = 1; sigma = 0.5; rho = 0.8;
with pytest.warns(UserWarning):
mc_rouwenhorst = rouwenhorst(n, rho, sigma, mu)
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)
Expand Down

0 comments on commit 0e21854

Please sign in to comment.