diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e79b3eab7..ac6825d88 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -53,7 +53,7 @@ jobs: - name: Setup QuantEcon shell: bash -l {0} run: | - pip install -U pip + python -m pip install -U pip pip install .[testing] - name: flake8 Tests diff --git a/quantecon/markov/__init__.py b/quantecon/markov/__init__.py index b4cca45df..ce906da51 100644 --- a/quantecon/markov/__init__.py +++ b/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 diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index e00b55e98..a2c7184d9 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -1,12 +1,15 @@ """ -tauchen -------- -Discretizes Gaussian linear AR(1) processes via Tauchen's method +Collection of functions to approximate a continuous random process by a +discrete Markov chain. """ 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 +245,198 @@ 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""" + Generate an `MarkovChain` instance that discretizes a multivariate + autorregressive process by a simulation of the process. + + This function 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; and :math:`{C}` is a volatility matrix. + Internally, from this process a sample time series of length + `sim_length` is produced, and with a cartesian grid as specified by + `grid_sizes` and `std_devs` a Markov chain is estimated that fits + the time series, where the states that are never visited under the + simulation are removed. + + 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. In particular, we follow Schmitt-Grohé + and Uribe's method in contructing the grid for approximation. + + Parameters + ---------- + A : array_like(float, ndim=2) + An m x m matrix containing the process' autocorrelation + parameters. Its eigenvalues are assumed to have moduli + bounded by unity. + C : array_like(float, ndim=2) + An m x r volatility matrix + grid_sizes : array_like(int, ndim=1), optional(default=None) + An m-vector containing the number of grid points in the + discretization of each dimension of x_t. If None, then 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 : optional(default=None) + Object that represents the disturbance term u_t. If None, then + standard normal distribution from numpy.random is used. + Alternatively, one can pass a "frozen" object of a multivariate + distribution from `scipy.stats`. It must have a zero mean and + unit standard deviation (of dimension m). + order : str, optional(default='C') + ('C' or 'F') order in which the states in the cartesian grid are + enumerated. + random_state : int or np.random.RandomState/Generator, optional + Random seed (integer) or np.random.RandomState or Generator + instance to set the initial state of the random number generator + for reproducibility. If None, a randomly initialized RandomState + is used. + + Returns + ------- + mc : MarkovChain + An instance of the MarkovChain class that stores the transition + matrix and state values returned by the discretization method, + in the following attributes: + + `P` : ndarray(float, ndim=2) + A 2-dim array containing the transition probability + matrix over the discretized states. + + `state_values` : ndarray(float, ndim=2) + A 2-dim array containing the state vectors (of dimension + m) as rows, which are ordered according to the `order` + option. + + Examples + -------- + 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. + + >>> rng = np.random.default_rng(12345) + >>> A = [[0.7901, -1.3570], + ... [-0.0104, 0.8638]] + >>> Omega = [[0.0012346, -0.0000776], + ... [-0.0000776, 0.0000401]] + >>> C = scipy.linalg.sqrtm(Omega) + >>> grid_sizes = [21, 11] + >>> mc = discrete_var(A, C, grid_sizes, random_state=rng) + >>> mc.P.shape + (145, 145) + >>> mc.state_values.shape + (145, 2) + >>> mc.state_values[:10] # First 10 states + array([[-0.38556417, 0.02155098], + [-0.38556417, 0.03232648], + [-0.38556417, 0.04310197], + [-0.38556417, 0.05387746], + [-0.34700776, 0.01077549], + [-0.34700776, 0.02155098], + [-0.34700776, 0.03232648], + [-0.34700776, 0.04310197], + [-0.34700776, 0.05387746], + [-0.30845134, 0. ]]) + >>> mc.simulate(10, random_state=rng) + array([[ 0.11566925, -0.01077549], + [ 0.11566925, -0.01077549], + [ 0.15422567, 0. ], + [ 0.15422567, 0. ], + [ 0.15422567, -0.01077549], + [ 0.11566925, -0.02155098], + [ 0.11566925, -0.03232648], + [ 0.15422567, -0.03232648], + [ 0.15422567, -0.03232648], + [ 0.19278209, -0.03232648]]) + + The simulation below uses the same parameters with :math:`{u_t}` + drawn from a multivariate t-distribution + + >>> df = 100 + >>> Sigma = np.diag(np.tile((df-2)/df, 2)) + >>> mc = discrete_var(A, C, grid_sizes, + ... rv=scipy.stats.multivariate_t(shape=Sigma, df=df), + ... random_state=rng) + >>> mc.P.shape + (146, 146) + >>> mc.state_values.shape + (146, 2) + >>> mc.state_values[:10] + array([[-0.38556417, 0.02155098], + [-0.38556417, 0.03232648], + [-0.38556417, 0.04310197], + [-0.38556417, 0.05387746], + [-0.34700776, 0.01077549], + [-0.34700776, 0.02155098], + [-0.34700776, 0.03232648], + [-0.34700776, 0.04310197], + [-0.34700776, 0.05387746], + [-0.30845134, 0. ]]) + >>> mc.simulate(10, random_state=rng) + array([[-0.03855642, -0.02155098], + [ 0.03855642, -0.03232648], + [ 0.07711283, -0.03232648], + [ 0.15422567, -0.03232648], + [ 0.15422567, -0.04310197], + [ 0.15422567, -0.03232648], + [ 0.15422567, -0.03232648], + [ 0.2313385 , -0.04310197], + [ 0.2313385 , -0.03232648], + [ 0.26989492, -0.03232648]]) + """ + 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, 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) + + return mc diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 2bd6ef24c..ed3268683 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/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,124 @@ 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.P_out_non_square = [ + [3.70370370e-02, 2.22222222e-01, 0.00000000e+00, 4.25925926e-01, + 3.14814815e-01, 0.00000000e+00], + [6.92865939e-05, 8.50870664e-01, 3.83177215e-02, 4.17954615e-04, + 1.10275202e-01, 4.91711312e-05], + [0.00000000e+00, 3.08160000e-01, 6.91360000e-01, 0.00000000e+00, + 3.91111111e-04, 8.88888889e-05], + [1.08405001e-04, 5.05890005e-04, 0.00000000e+00, 6.95707162e-01, + 3.03678543e-01, 0.00000000e+00], + [3.40242525e-05, 1.11864937e-01, 4.44583566e-04, 3.77260912e-02, + 8.49844169e-01, 8.61947730e-05], + [0.00000000e+00, 4.70588235e-01, 3.08823529e-01, 0.00000000e+00, + 1.76470588e-01, 4.41176471e-02]] + + self.A, self.C, self.S_out, self.P_out, self.S_out_orderF,\ + self.P_out_orderF, self.P_out_non_square \ + = map(np.array,(self.A, self.C, self.S_out, self.P_out, + self.S_out_orderF, self.P_out_orderF, + self.P_out_non_square)) + + 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(cov=np.identity(2)), + 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) + + def test_order_non_squareC(self): + new_col = np.array([0, 0]) + self.C = np.insert(self.C, 2, new_col, axis=1) + mc = discrete_var( + self.A, self.C, + grid_sizes=self.sizes, + sim_length=self.T, + order='C', + random_state=self.random_state) + assert_allclose(mc.state_values, self.S_out) + assert_allclose(mc.P, self.P_out_non_square)