From 4ae4ee5bd7bd037d0ef43e92cb17efa914ac99cd Mon Sep 17 00:00:00 2001 From: Smit-create Date: Thu, 3 Nov 2022 19:48:35 +0530 Subject: [PATCH 01/25] Add discrete_var MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: John Stachurski Co-Authored-By: Carlos Rondón-Moreno <31420139+crondonm@users.noreply.github.com> --- quantecon/markov/approximation.py | 200 ++++++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 2256fadba..9a8af5be6 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -7,8 +7,10 @@ from math import erfc, sqrt from .core import MarkovChain +from quantecon import matrix_eqn as qme import numpy as np +import scipy as sp from numba import njit @@ -212,3 +214,201 @@ 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 cartesian_product(arrays, row_major=True): + """ + Create a Cartesian product from the list of grids in `arrays` and then + flatten it so that S[i, :] is an array of grid points corresponding to + state i. + + If row_major is True, then the Cartesian product of the grids is + enumerated in row major order. Otherwise it is enumerated in column major. + + Currently the column major operation is copied from MATLAB and could + surely be made more efficient. + """ + + m = len(arrays) + if row_major: + s = np.stack(np.meshgrid(*arrays, indexing='ij'), axis=-1) + S = np.reshape(s, (-1, m)) + else: + V = arrays + gs = [len(v) for v in V] + n = np.prod(gs) + S = np.zeros((n, m)) # Discrete state space + + for i in range(m): + if i == 0: + S0 = np.ravel(V[i]) + S[:, i] = np.ravel(np.tile(S0, [np.prod(gs[i+1:]), 1])) + else: + S0 = np.sort(np.ravel(np.tile(V[i], [np.prod(gs[0:i]), 1]))) + S[:, i] = np.ravel(np.tile(S0, [np.prod(gs[i+1:]), 1])) + return S + + +def discrete_var(A, + Omega, + grid_sizes=None, + std_devs=np.sqrt(10), + seed=1234, + sim_length=1_000_000, + burn_in=100_000, + row_major=True, + return_sim=False): + r""" + This code discretizes a VAR(1) process of the form: + + .. math:: + + x_t = A x_{t-1} + u_t + + where :math:`{u_t}` is zero-mean Gaussian with variance-covariance + matrix Omega. + + 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 by Carlos Rondón-Moreno from Schmitt-Grohé and + Uribe's code for MATLAB. + + Parameters + ---------- + A : array_like(float) + An m x m matrix containing the process' autocorrelation parameters + Omega : array_like(float) + An m x m variance-covariance matrix + grid_sizes : array_like(int) or None + An m-vector containing the number of grid points in the discretization + of each dimension of x_t. If grid_sizes is None, then grid_sizes is + set to (10, ..., 10). + std_devs : float + The number of standard deviations the grid should stretch in each + dimension, where standard deviations are measured under the stationary + distribution. + sim_length : int + The the length of the simulated time series (default, 1_000_000). + burn_in : int + The number of burn-in draws from the simulated series (default, + 100_000). + row_major : bool + If True then return form the cartesian product state grid using row + major ordering. Otherwise use column major. + return_sim : bool + If True then return the state space simulation generated by the model. + + Returns + ------- + Pi : array_like(float) + A square matrix containing the transition probability + matrix of the discretized state. + S : array_like(float) + 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. + Xvec : array_like(float) + A matrix of size m x sim_length containing the simulated time + series of the m discretized states. + + + Notes + ----- + The code presently assumes normal shocks but normality is not required + for the algorithm to work. The draws from the multivariate standard + normal generator can be replaced by any other random number generator + with mean 0 and unit standard deviation. + + + 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]]) + grid_sizes = np.array([21, 11]) + Pi, Xvec, S = discrete_var(A, Omega, grid_sizes, + sim_length=1_000_000, burn_in = 100_000) + """ + + m = len(A) # The number of dimensions of the original state x_t + default_grid_size = 10 + + if grid_sizes is None: + # Set the size of every grid to default_grid_size + grid_sizes = np.full(m, default_grid_size) + + n = grid_sizes.prod() # Size of the discretized state + + # Compute stationary variance-covariance matrix of AR process and use + # it to obtain grid bounds. + Sigma = qme.solve_discrete_lyapunov(A, Omega) + sigma_vector = np.sqrt(np.diagonal(Sigma)) # Stationary std dev + upper_bounds = std_devs * sigma_vector + + # Build the individual grids along each dimension + V = [] + for i in range(m): + b = np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i]) + V.append(b) + + S = cartesian_product(V, row_major=row_major) + + Pi = np.zeros((n, n)) + Xvec = np.zeros((m, sim_length)) + C = sp.linalg.sqrtm(Omega) + + # Run simulation to compute transition probabilities + _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed) + + # Cut states where the column sum of Pi is zero (i.e., inaccesible states + # according to the simulation) + indx = np.where(np.sum(Pi, axis=0) > 0) + Pi = Pi[indx[0], :] + Pi = Pi[:, indx[0]] + S = S[indx[0], :] + + # Normalize + sum_row = np.sum(Pi, axis=1) + for i in range(len(Pi)): + Pi[i, :] = Pi[i, :] / sum_row[i] + + if return_sim: + return Pi, S, Xvec + return Pi, S + + +@njit +def _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed): + m = len(A) + np.random.seed(seed) + x0 = np.zeros((m, 1)) + d = np.sum(S**2, axis=1) + ind_i = np.argmin(d) + + for t in range(sim_length + burn_in): + # Update state + drw = C @ np.random.randn(m, 1) + x = A @ x0 + drw + # Find the index of the state closest to x + xx = np.reshape(x, (1, m)) + d = np.sum((S - xx)**2, axis=1) + ind_j = np.argmin(d) + + if t > burn_in: + Pi[ind_i, ind_j] += 1 + Xvec[:, t-burn_in] = x.T + x0 = x + ind_i = ind_j From 4709925c62d1cbbb7ed78f41d36e491df2f0d81a Mon Sep 17 00:00:00 2001 From: Smit-create Date: Thu, 3 Nov 2022 19:49:07 +0530 Subject: [PATCH 02/25] Add suggestions by @crondonm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Carlos Rondón-Moreno <31420139+crondonm@users.noreply.github.com> --- quantecon/markov/approximation.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 9a8af5be6..22ea606d3 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -275,8 +275,7 @@ def discrete_var(A, VAR Processes: A Simulation Approach* by Stephanie Schmitt-Grohé and Martín Uribe, July 11, 2010. - This code was adapted by Carlos Rondón-Moreno from Schmitt-Grohé and - Uribe's code for MATLAB. + This code was adapted from Schmitt-Grohé and Uribe's original MATLAB code. Parameters ---------- From 8295961f7347d5e9badfbb2c9ea20703e8b1049c Mon Sep 17 00:00:00 2001 From: Smit-create Date: Thu, 3 Nov 2022 19:49:41 +0530 Subject: [PATCH 03/25] Add discrete_var in approximation --- quantecon/markov/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quantecon/markov/__init__.py b/quantecon/markov/__init__.py index 6e2a07df0..c8e60102c 100644 --- a/quantecon/markov/__init__.py +++ b/quantecon/markov/__init__.py @@ -8,6 +8,6 @@ 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 2e7526d554534e0f28b125eca1dea4cd26669ee3 Mon Sep 17 00:00:00 2001 From: Smit-create Date: Thu, 3 Nov 2022 19:51:51 +0530 Subject: [PATCH 04/25] Add tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: John Stachurski Co-Authored-By: Carlos Rondón-Moreno <31420139+crondonm@users.noreply.github.com> --- quantecon/markov/tests/test_approximation.py | 71 +++++++++++++++++++- 1 file changed, 68 insertions(+), 3 deletions(-) diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index b1d33d4f6..7a33b6ce6 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -3,11 +3,9 @@ """ import numpy as np -from quantecon.markov import tauchen, rouwenhorst +from quantecon.markov import tauchen, rouwenhorst, discrete_var from numpy.testing import assert_, assert_allclose -#from quantecon.markov.approximation import rouwenhorst - class TestTauchen: @@ -98,3 +96,70 @@ 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.T, self.burn_in = 1_000_000, 100_000 + + self.A = [[0.7901, -1.3570], + [-0.0104, 0.8638]] + self.Omega = [[0.0012346, -0.0000776], + [-0.0000776, 0.0000401]] + + self.sizes = np.array((2, 3)) + + # Expected outputs (column major) + self.S_out = [[-0.38556417, -0.05387746], + [ 0.38556417, -0.05387746], + [-0.38556417, 0. ], + [ 0.38556417, 0. ], + [-0.38556417, 0.05387746], + [ 0.38556417, 0.05387746]] + + self.Pi_out = [[8.06451613e-02, 3.54838710e-01, 1.93548387e-01, 3.70967742e-01, + 0.00000000e+00, 0.00000000e+00], + [1.03646634e-04, 6.99407487e-01, 5.52782048e-04, 2.99936085e-01, + 0.00000000e+00, 0.00000000e+00], + [8.39514352e-05, 4.96901738e-04, 8.50806955e-01, 1.10647992e-01, + 3.79097454e-02, 5.44549850e-05], + [3.14483775e-05, 3.85467256e-02, 1.09581871e-01, 8.51289608e-01, + 4.53755161e-04, 9.65914451e-05], + [0.00000000e+00, 0.00000000e+00, 3.01394785e-01, 5.70755895e-04, + 6.97927443e-01, 1.07016730e-04], + [0.00000000e+00, 0.00000000e+00, 3.63636364e-01, 2.46753247e-01, + 3.37662338e-01, 5.19480519e-02]] + + self.A, self.Omega, self.S_out, self.Pi_out = map(np.array, + (self.A, self.Omega, self.S_out, self.Pi_out)) + + def tearDown(self): + del self.A, self.Omega, self.S_out, self.Pi_out + + def test_column_major_discretization(self): + Pi, S = discrete_var( + self.A, self.Omega, grid_sizes=self.sizes, + sim_length=self.T, burn_in=self.burn_in, + row_major=False) + assert_allclose(S, self.S_out) + assert_allclose(Pi, self.Pi_out) + + def test_column_row_parity(self): + """ + Test that column major and row major discretization produces the same + output, up to a column-major / row-major reordering. + """ + Pi, S = discrete_var( + self.A, self.Omega, grid_sizes=self.sizes, + sim_length=self.T, burn_in=self.burn_in, + row_major=False) + Pi_r, S_r = discrete_var( + self.A, self.Omega, grid_sizes=self.sizes, + sim_length=self.T, burn_in=self.burn_in, + row_major=True) + + # State 1 under column major becomes state 3 under row major, given + # the size of the state space (multigrid state index is (1, 0)). + assert_allclose(S[1, :], S_r[3, :]) + assert_(Pi[1, 1] == Pi_r[3, 3]) From 84e751d8becdc4a066dc0370e18d04b257826588 Mon Sep 17 00:00:00 2001 From: Shu Date: Sun, 6 Nov 2022 01:04:01 +1100 Subject: [PATCH 05/25] add suggestions 1, 2, 5 by @jstac --- quantecon/markov/approximation.py | 26 ++++++-------------- quantecon/markov/tests/test_approximation.py | 14 +++++------ 2 files changed, 14 insertions(+), 26 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 22ea606d3..789ffd2fa 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -256,8 +256,7 @@ def discrete_var(A, seed=1234, sim_length=1_000_000, burn_in=100_000, - row_major=True, - return_sim=False): + row_major=True): r""" This code discretizes a VAR(1) process of the form: @@ -299,21 +298,12 @@ def discrete_var(A, row_major : bool If True then return form the cartesian product state grid using row major ordering. Otherwise use column major. - return_sim : bool - If True then return the state space simulation generated by the model. Returns ------- - Pi : array_like(float) - A square matrix containing the transition probability - matrix of the discretized state. - S : array_like(float) - 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. - Xvec : array_like(float) - A matrix of size m x sim_length containing the simulated time - series of the m discretized states. + mc : MarkovChain + An instance of the MarkovChain class that stores the transition + matrix and state values returned by the discretization method Notes @@ -338,7 +328,7 @@ def discrete_var(A, Omega = np.array([[0.0012346, -0.0000776], [-0.0000776, 0.0000401]]) grid_sizes = np.array([21, 11]) - Pi, Xvec, S = discrete_var(A, Omega, grid_sizes, + mc = discrete_var(A, Omega, grid_sizes, sim_length=1_000_000, burn_in = 100_000) """ @@ -384,9 +374,8 @@ def discrete_var(A, for i in range(len(Pi)): Pi[i, :] = Pi[i, :] / sum_row[i] - if return_sim: - return Pi, S, Xvec - return Pi, S + mc = MarkovChain(Pi, state_values=S) + return mc @njit @@ -408,6 +397,5 @@ def _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed): if t > burn_in: Pi[ind_i, ind_j] += 1 - Xvec[:, t-burn_in] = x.T x0 = x ind_i = ind_j diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 7a33b6ce6..2569285e4 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -138,28 +138,28 @@ def tearDown(self): del self.A, self.Omega, self.S_out, self.Pi_out def test_column_major_discretization(self): - Pi, S = discrete_var( + mc = discrete_var( self.A, self.Omega, grid_sizes=self.sizes, sim_length=self.T, burn_in=self.burn_in, row_major=False) - assert_allclose(S, self.S_out) - assert_allclose(Pi, self.Pi_out) + assert_allclose(mc.state_values, self.S_out) + assert_allclose(mc.P, self.Pi_out) def test_column_row_parity(self): """ Test that column major and row major discretization produces the same output, up to a column-major / row-major reordering. """ - Pi, S = discrete_var( + mc = discrete_var( self.A, self.Omega, grid_sizes=self.sizes, sim_length=self.T, burn_in=self.burn_in, row_major=False) - Pi_r, S_r = discrete_var( + mc_r = discrete_var( self.A, self.Omega, grid_sizes=self.sizes, sim_length=self.T, burn_in=self.burn_in, row_major=True) # State 1 under column major becomes state 3 under row major, given # the size of the state space (multigrid state index is (1, 0)). - assert_allclose(S[1, :], S_r[3, :]) - assert_(Pi[1, 1] == Pi_r[3, 3]) + assert_allclose(mc.state_values[1, :], mc_r.state_values[3, :]) + assert_(mc.P[1, 1] == mc_r.P[3, 3]) From ba099268edaae192d0d6a2c7b37a91c95d7e4caf Mon Sep 17 00:00:00 2001 From: Smit-create Date: Sun, 6 Nov 2022 11:45:12 +0530 Subject: [PATCH 06/25] Add docs of Pi and S in output. --- quantecon/markov/approximation.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 789ffd2fa..965f02593 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -301,9 +301,16 @@ def discrete_var(A, Returns ------- + mc : MarkovChain An instance of the MarkovChain class that stores the transition - matrix and state values returned by the discretization method + matrix and state values returned by the discretization method. + The MarkovChain instance contains: + Pi : 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. Notes From cc5217c6b155427a1a727d19a3579a15a71e1655 Mon Sep 17 00:00:00 2001 From: Smit-create Date: Sun, 6 Nov 2022 11:49:14 +0530 Subject: [PATCH 07/25] Remove row_major argument and use cartesian from gridtools --- quantecon/markov/approximation.py | 45 ++++--------------------------- 1 file changed, 5 insertions(+), 40 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 965f02593..fe09e5809 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -8,6 +8,7 @@ from math import erfc, sqrt from .core import MarkovChain from quantecon import matrix_eqn as qme +from quantecon.gridtools import cartesian import numpy as np import scipy as sp @@ -216,47 +217,13 @@ def _fill_tauchen(x, P, n, rho, sigma, half_step): std_norm_cdf((z - half_step) / sigma)) -def cartesian_product(arrays, row_major=True): - """ - Create a Cartesian product from the list of grids in `arrays` and then - flatten it so that S[i, :] is an array of grid points corresponding to - state i. - - If row_major is True, then the Cartesian product of the grids is - enumerated in row major order. Otherwise it is enumerated in column major. - - Currently the column major operation is copied from MATLAB and could - surely be made more efficient. - """ - - m = len(arrays) - if row_major: - s = np.stack(np.meshgrid(*arrays, indexing='ij'), axis=-1) - S = np.reshape(s, (-1, m)) - else: - V = arrays - gs = [len(v) for v in V] - n = np.prod(gs) - S = np.zeros((n, m)) # Discrete state space - - for i in range(m): - if i == 0: - S0 = np.ravel(V[i]) - S[:, i] = np.ravel(np.tile(S0, [np.prod(gs[i+1:]), 1])) - else: - S0 = np.sort(np.ravel(np.tile(V[i], [np.prod(gs[0:i]), 1]))) - S[:, i] = np.ravel(np.tile(S0, [np.prod(gs[i+1:]), 1])) - return S - - def discrete_var(A, Omega, grid_sizes=None, std_devs=np.sqrt(10), seed=1234, sim_length=1_000_000, - burn_in=100_000, - row_major=True): + burn_in=100_000): r""" This code discretizes a VAR(1) process of the form: @@ -295,9 +262,6 @@ def discrete_var(A, burn_in : int The number of burn-in draws from the simulated series (default, 100_000). - row_major : bool - If True then return form the cartesian product state grid using row - major ordering. Otherwise use column major. Returns ------- @@ -310,7 +274,8 @@ def discrete_var(A, 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. + unique values yields the grid values. The cartesian product + state grid uses a row major ordering. Notes @@ -360,7 +325,7 @@ def discrete_var(A, b = np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i]) V.append(b) - S = cartesian_product(V, row_major=row_major) + S = cartesian(V) Pi = np.zeros((n, n)) Xvec = np.zeros((m, sim_length)) From c26d265323a5088bd04423b43a5edc383d5b1e77 Mon Sep 17 00:00:00 2001 From: Smit-create Date: Sun, 6 Nov 2022 11:56:18 +0530 Subject: [PATCH 08/25] Modify tests to drop row_major --- quantecon/markov/tests/test_approximation.py | 62 +++++++------------- 1 file changed, 21 insertions(+), 41 deletions(-) diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 2569285e4..6225f5165 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -110,26 +110,26 @@ def setup_method(self): self.sizes = np.array((2, 3)) - # Expected outputs (column major) + # Expected outputs self.S_out = [[-0.38556417, -0.05387746], - [ 0.38556417, -0.05387746], - [-0.38556417, 0. ], - [ 0.38556417, 0. ], - [-0.38556417, 0.05387746], - [ 0.38556417, 0.05387746]] - - self.Pi_out = [[8.06451613e-02, 3.54838710e-01, 1.93548387e-01, 3.70967742e-01, - 0.00000000e+00, 0.00000000e+00], - [1.03646634e-04, 6.99407487e-01, 5.52782048e-04, 2.99936085e-01, - 0.00000000e+00, 0.00000000e+00], - [8.39514352e-05, 4.96901738e-04, 8.50806955e-01, 1.10647992e-01, - 3.79097454e-02, 5.44549850e-05], - [3.14483775e-05, 3.85467256e-02, 1.09581871e-01, 8.51289608e-01, - 4.53755161e-04, 9.65914451e-05], - [0.00000000e+00, 0.00000000e+00, 3.01394785e-01, 5.70755895e-04, - 6.97927443e-01, 1.07016730e-04], - [0.00000000e+00, 0.00000000e+00, 3.63636364e-01, 2.46753247e-01, - 3.37662338e-01, 5.19480519e-02]] + [-0.38556417, 0. ], + [-0.38556417, 0.05387746], + [ 0.38556417, -0.05387746], + [ 0.38556417, 0. ], + [ 0.38556417, 0.05387746]] + + self.Pi_out = [[8.06451613e-02, 1.93548387e-01, 0.00000000e+00, 3.54838710e-01, + 3.70967742e-01, 0.00000000e+00], + [8.39514352e-05, 8.50806955e-01, 3.79097454e-02, 4.96901738e-04, + 1.10647992e-01, 5.44549850e-05], + [0.00000000e+00, 3.01394785e-01, 6.97927443e-01, 0.00000000e+00, + 5.70755895e-04, 1.07016730e-04], + [1.03646634e-04, 5.52782048e-04, 0.00000000e+00, 6.99407487e-01, + 2.99936085e-01, 0.00000000e+00], + [3.14483775e-05, 1.09581871e-01, 4.53755161e-04, 3.85467256e-02, + 8.51289608e-01, 9.65914451e-05], + [0.00000000e+00, 3.63636364e-01, 3.37662338e-01, 0.00000000e+00, + 2.46753247e-01, 5.19480519e-02]] self.A, self.Omega, self.S_out, self.Pi_out = map(np.array, (self.A, self.Omega, self.S_out, self.Pi_out)) @@ -137,29 +137,9 @@ def setup_method(self): def tearDown(self): del self.A, self.Omega, self.S_out, self.Pi_out - def test_column_major_discretization(self): + def test_discretization(self): mc = discrete_var( self.A, self.Omega, grid_sizes=self.sizes, - sim_length=self.T, burn_in=self.burn_in, - row_major=False) + sim_length=self.T, burn_in=self.burn_in) assert_allclose(mc.state_values, self.S_out) assert_allclose(mc.P, self.Pi_out) - - def test_column_row_parity(self): - """ - Test that column major and row major discretization produces the same - output, up to a column-major / row-major reordering. - """ - mc = discrete_var( - self.A, self.Omega, grid_sizes=self.sizes, - sim_length=self.T, burn_in=self.burn_in, - row_major=False) - mc_r = discrete_var( - self.A, self.Omega, grid_sizes=self.sizes, - sim_length=self.T, burn_in=self.burn_in, - row_major=True) - - # State 1 under column major becomes state 3 under row major, given - # the size of the state space (multigrid state index is (1, 0)). - assert_allclose(mc.state_values[1, :], mc_r.state_values[3, :]) - assert_(mc.P[1, 1] == mc_r.P[3, 3]) From 00f457eae7f7449daf6c7971084737372d642473 Mon Sep 17 00:00:00 2001 From: Smit-create Date: Sun, 6 Nov 2022 11:59:32 +0530 Subject: [PATCH 09/25] Rename Pi to P --- quantecon/markov/approximation.py | 26 ++++++++++---------- quantecon/markov/tests/test_approximation.py | 12 ++++----- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index fe09e5809..e2fff7abb 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -270,7 +270,7 @@ def discrete_var(A, An instance of the MarkovChain class that stores the transition matrix and state values returned by the discretization method. The MarkovChain instance contains: - Pi : A square matrix containing the transition probability + 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 @@ -327,31 +327,31 @@ def discrete_var(A, S = cartesian(V) - Pi = np.zeros((n, n)) + P = np.zeros((n, n)) Xvec = np.zeros((m, sim_length)) C = sp.linalg.sqrtm(Omega) # Run simulation to compute transition probabilities - _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed) + _run_sim(A, C, P, Xvec, S, sim_length, burn_in, seed) - # Cut states where the column sum of Pi is zero (i.e., inaccesible states + # Cut states where the column sum of P is zero (i.e., inaccesible states # according to the simulation) - indx = np.where(np.sum(Pi, axis=0) > 0) - Pi = Pi[indx[0], :] - Pi = Pi[:, indx[0]] + indx = np.where(np.sum(P, axis=0) > 0) + P = P[indx[0], :] + P = P[:, indx[0]] S = S[indx[0], :] # Normalize - sum_row = np.sum(Pi, axis=1) - for i in range(len(Pi)): - Pi[i, :] = Pi[i, :] / sum_row[i] + sum_row = np.sum(P, axis=1) + for i in range(len(P)): + P[i, :] = P[i, :] / sum_row[i] - mc = MarkovChain(Pi, state_values=S) + mc = MarkovChain(P, state_values=S) return mc @njit -def _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed): +def _run_sim(A, C, P, Xvec, S, sim_length, burn_in, seed): m = len(A) np.random.seed(seed) x0 = np.zeros((m, 1)) @@ -368,6 +368,6 @@ def _run_sim(A, C, Pi, Xvec, S, sim_length, burn_in, seed): ind_j = np.argmin(d) if t > burn_in: - Pi[ind_i, ind_j] += 1 + P[ind_i, ind_j] += 1 x0 = x ind_i = ind_j diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 6225f5165..7dd6bb121 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -118,7 +118,7 @@ def setup_method(self): [ 0.38556417, 0. ], [ 0.38556417, 0.05387746]] - self.Pi_out = [[8.06451613e-02, 1.93548387e-01, 0.00000000e+00, 3.54838710e-01, + self.P_out = [[8.06451613e-02, 1.93548387e-01, 0.00000000e+00, 3.54838710e-01, 3.70967742e-01, 0.00000000e+00], [8.39514352e-05, 8.50806955e-01, 3.79097454e-02, 4.96901738e-04, 1.10647992e-01, 5.44549850e-05], @@ -131,15 +131,15 @@ def setup_method(self): [0.00000000e+00, 3.63636364e-01, 3.37662338e-01, 0.00000000e+00, 2.46753247e-01, 5.19480519e-02]] - self.A, self.Omega, self.S_out, self.Pi_out = map(np.array, - (self.A, self.Omega, self.S_out, self.Pi_out)) + self.A, self.Omega, self.S_out, self.P_out = map(np.array, + (self.A, self.Omega, self.S_out, self.P_out)) - def tearDown(self): - del self.A, self.Omega, self.S_out, self.Pi_out + def teardown_method(self): + del self.A, self.Omega, self.S_out, self.P_out def test_discretization(self): mc = discrete_var( self.A, self.Omega, grid_sizes=self.sizes, sim_length=self.T, burn_in=self.burn_in) assert_allclose(mc.state_values, self.S_out) - assert_allclose(mc.P, self.Pi_out) + assert_allclose(mc.P, self.P_out) From 740cf739b884124afd182cf8f792fb7d00ac2801 Mon Sep 17 00:00:00 2001 From: Shu Date: Sun, 6 Nov 2022 21:06:32 +1100 Subject: [PATCH 10/25] use mlinspace --- quantecon/markov/approximation.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index e2fff7abb..86bed0fa4 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -8,7 +8,7 @@ from math import erfc, sqrt from .core import MarkovChain from quantecon import matrix_eqn as qme -from quantecon.gridtools import cartesian +from quantecon.gridtools import mlinspace import numpy as np import scipy as sp @@ -320,12 +320,7 @@ def discrete_var(A, upper_bounds = std_devs * sigma_vector # Build the individual grids along each dimension - V = [] - for i in range(m): - b = np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i]) - V.append(b) - - S = cartesian(V) + S = mlinspace(-upper_bounds, upper_bounds, grid_sizes) P = np.zeros((n, n)) Xvec = np.zeros((m, sim_length)) From 745c0d447a95277aadbbc66c50a6dc9656bafe36 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sat, 17 Dec 2022 13:04:56 +0900 Subject: [PATCH 11/25] Implement discrete_var using quantecon functions --- quantecon/markov/approximation.py | 87 ++++++++++++------------------- 1 file changed, 34 insertions(+), 53 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 5ab31e0a1..55b5c25cb 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -7,8 +7,11 @@ from math import erfc, sqrt from .core import MarkovChain -from quantecon import matrix_eqn as qme -from quantecon.gridtools import mlinspace +from .estimate import estimate_mc +from ..gridtools import cartesian, cartesian_nearest_index +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 @@ -244,12 +247,13 @@ def _fill_tauchen(x, P, n, rho, sigma, half_step): def discrete_var(A, - Omega, + C, grid_sizes=None, std_devs=np.sqrt(10), - seed=1234, sim_length=1_000_000, - burn_in=100_000): + rv=None, + order='C', + random_state=None): r""" This code discretizes a VAR(1) process of the form: @@ -329,66 +333,43 @@ def discrete_var(A, mc = discrete_var(A, Omega, grid_sizes, sim_length=1_000_000, burn_in = 100_000) """ + A = np.asarray(A) + C = np.asarray(C) + m, r = C.shape - m = len(A) # The number of dimensions of the original state x_t - default_grid_size = 10 + # Run simulation to compute transition probabilities + random_state = check_random_state(random_state) - if grid_sizes is None: - # Set the size of every grid to default_grid_size - grid_sizes = np.full(m, default_grid_size) + 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) - n = grid_sizes.prod() # Size of the discretized 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 = qme.solve_discrete_lyapunov(A, Omega) + 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 - S = mlinspace(-upper_bounds, upper_bounds, grid_sizes) - - P = np.zeros((n, n)) - Xvec = np.zeros((m, sim_length)) - C = sp.linalg.sqrtm(Omega) + 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) - # Run simulation to compute transition probabilities - _run_sim(A, C, P, Xvec, S, sim_length, burn_in, seed) + V = [np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i]) + for i in range(m)] - # Cut states where the column sum of P is zero (i.e., inaccesible states - # according to the simulation) - indx = np.where(np.sum(P, axis=0) > 0) - P = P[indx[0], :] - P = P[:, indx[0]] - S = S[indx[0], :] + # Estimate the Markov chain + X_indices = cartesian_nearest_index(X.T, V, order=order) + mc = estimate_mc(X_indices) - # Normalize - sum_row = np.sum(P, axis=1) - for i in range(len(P)): - P[i, :] = P[i, :] / sum_row[i] + # Assign the visited states in the cartesian product as the state values + prod = cartesian(V, order=order) + mc.state_values = prod[mc.state_values] - mc = MarkovChain(P, state_values=S) return mc - - -@njit -def _run_sim(A, C, P, Xvec, S, sim_length, burn_in, seed): - m = len(A) - np.random.seed(seed) - x0 = np.zeros((m, 1)) - d = np.sum(S**2, axis=1) - ind_i = np.argmin(d) - - for t in range(sim_length + burn_in): - # Update state - drw = C @ np.random.randn(m, 1) - x = A @ x0 + drw - # Find the index of the state closest to x - xx = np.reshape(x, (1, m)) - d = np.sum((S - xx)**2, axis=1) - ind_j = np.argmin(d) - - if t > burn_in: - P[ind_i, ind_j] += 1 - x0 = x - ind_i = ind_j From f4ebd6f89b949a3c917204c58fc2ede5d7a781a6 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Fri, 6 Jan 2023 14:40:25 +0900 Subject: [PATCH 12/25] Use `fit_discrete_mc` --- quantecon/markov/approximation.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 796a5036e..36d360c20 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -7,7 +7,7 @@ from math import erfc, sqrt from .core import MarkovChain -from .estimate import estimate_mc +from .estimate import fit_discrete_mc from ..gridtools import cartesian, cartesian_nearest_index from ..lss import simulate_linear_model from ..matrix_eqn import solve_discrete_lyapunov @@ -367,12 +367,7 @@ def discrete_var(A, V = [np.linspace(-upper_bounds[i], upper_bounds[i], grid_sizes[i]) for i in range(m)] - # Estimate the Markov chain - X_indices = cartesian_nearest_index(X.T, V, order=order) - mc = estimate_mc(X_indices) - - # Assign the visited states in the cartesian product as the state values - prod = cartesian(V, order=order) - mc.state_values = prod[mc.state_values] + # Fit the Markov chain + mc = fit_discrete_mc(X.T, V, order='C') return mc From 7f7eec46c97947d40b1298093214f96c42fbf9e7 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sun, 8 Jan 2023 20:49:56 +0900 Subject: [PATCH 13/25] Remove unused import --- quantecon/markov/approximation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 36d360c20..8da576cf0 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -8,7 +8,7 @@ from math import erfc, sqrt from .core import MarkovChain from .estimate import fit_discrete_mc -from ..gridtools import cartesian, cartesian_nearest_index +from ..gridtools import cartesian from ..lss import simulate_linear_model from ..matrix_eqn import solve_discrete_lyapunov from ..util import check_random_state From 951595168fdb52b2396d22866c9292fc10eca3d3 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 9 Jan 2023 12:40:59 +0900 Subject: [PATCH 14/25] Remove unused import --- quantecon/markov/approximation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 8da576cf0..4fc122ae1 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -8,7 +8,6 @@ from math import erfc, sqrt from .core import MarkovChain from .estimate import fit_discrete_mc -from ..gridtools import cartesian from ..lss import simulate_linear_model from ..matrix_eqn import solve_discrete_lyapunov from ..util import check_random_state From bb918643f5d9e681ff63795b8886da82f8eb4227 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Fri, 13 Jan 2023 00:40:09 +1100 Subject: [PATCH 15/25] adjust test and docstring --- quantecon/markov/approximation.py | 63 ++++++++------------ quantecon/markov/tests/test_approximation.py | 51 ++++++++-------- 2 files changed, 51 insertions(+), 63 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 4fc122ae1..6f6f0b8f8 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -258,46 +258,44 @@ def discrete_var(A, random_state=None): r""" This code discretizes a VAR(1) process of the form: - .. math:: - - x_t = A x_{t-1} + u_t - - where :math:`{u_t}` is zero-mean Gaussian with variance-covariance - matrix Omega. - + 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 - Omega : array_like(float) - An m x m variance-covariance matrix - grid_sizes : array_like(int) or None + 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 grid_sizes is None, then grid_sizes is + of each dimension of x_t. If None, then grid_sizes is set to (10, ..., 10). - std_devs : float + 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 - The the length of the simulated time series (default, 1_000_000). - burn_in : int - The number of burn-in draws from the simulated series (default, - 100_000). - + 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. + 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. @@ -308,32 +306,21 @@ def discrete_var(A, 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. - - - Notes - ----- - The code presently assumes normal shocks but normality is not required - for the algorithm to work. The draws from the multivariate standard - normal generator can be replaced by any other random number generator - with mean 0 and unit standard deviation. - - + 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, Omega, grid_sizes, - sim_length=1_000_000, burn_in = 100_000) + mc = discrete_var(A, C, grid_sizes, sim_length=1_000_000) """ A = np.asarray(A) C = np.asarray(C) @@ -367,6 +354,6 @@ def discrete_var(A, for i in range(m)] # Fit the Markov chain - mc = fit_discrete_mc(X.T, V, order='C') + mc = fit_discrete_mc(X.T, V, order=order) - return mc + return mc \ No newline at end of file diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index ee10c5c35..453e0bb1c 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -6,7 +6,7 @@ import pytest from quantecon.markov import tauchen, rouwenhorst, discrete_var from numpy.testing import assert_, assert_allclose, assert_raises - +import scipy as sp class TestTauchen: @@ -110,13 +110,14 @@ def test_control_case(self): class TestDiscreteVar: def setup_method(self): - self.T, self.burn_in = 1_000_000, 100_000 + self.random_state = np.random.RandomState(0) - self.A = [[0.7901, -1.3570], + Omega = np.array([[0.0012346, -0.0000776], + [-0.0000776, 0.0000401]]) + self.A = [[0.7901, -1.3570], [-0.0104, 0.8638]] - self.Omega = [[0.0012346, -0.0000776], - [-0.0000776, 0.0000401]] - + self.C = sp.linalg.sqrtm(Omega) + self.T= 1_000_000 self.sizes = np.array((2, 3)) # Expected outputs @@ -127,28 +128,28 @@ def setup_method(self): [ 0.38556417, 0. ], [ 0.38556417, 0.05387746]] - self.P_out = [[8.06451613e-02, 1.93548387e-01, 0.00000000e+00, 3.54838710e-01, - 3.70967742e-01, 0.00000000e+00], - [8.39514352e-05, 8.50806955e-01, 3.79097454e-02, 4.96901738e-04, - 1.10647992e-01, 5.44549850e-05], - [0.00000000e+00, 3.01394785e-01, 6.97927443e-01, 0.00000000e+00, - 5.70755895e-04, 1.07016730e-04], - [1.03646634e-04, 5.52782048e-04, 0.00000000e+00, 6.99407487e-01, - 2.99936085e-01, 0.00000000e+00], - [3.14483775e-05, 1.09581871e-01, 4.53755161e-04, 3.85467256e-02, - 8.51289608e-01, 9.65914451e-05], - [0.00000000e+00, 3.63636364e-01, 3.37662338e-01, 0.00000000e+00, - 2.46753247e-01, 5.19480519e-02]] - - self.A, self.Omega, self.S_out, self.P_out = map(np.array, - (self.A, self.Omega, self.S_out, self.P_out)) + 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.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.Omega, self.S_out, self.P_out + del self.A, self.C, self.S_out, self.P_out def test_discretization(self): mc = discrete_var( - self.A, self.Omega, grid_sizes=self.sizes, - sim_length=self.T, burn_in=self.burn_in) + 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) + assert_allclose(mc.P, self.P_out) \ No newline at end of file From 03e1a8ddcedcfeabe1df9ee31c651942404a9d23 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sun, 22 Jan 2023 22:19:19 +0900 Subject: [PATCH 16/25] Edit docstring --- quantecon/markov/approximation.py | 147 +++++++++++++++++++----------- 1 file changed, 96 insertions(+), 51 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 6f6f0b8f8..7f5eca242 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -1,7 +1,6 @@ """ -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. """ @@ -257,42 +256,59 @@ def discrete_var(A, order='C', random_state=None): r""" - This code discretizes a VAR(1) process of the form: + 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; - :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. + + 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. - This code was adapted from Schmitt-Grohé and Uribe's original MATLAB code. + 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) - 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). + A : array_like(float, ndim=2) + An m x m matrix containing the process' autocorrelation + parameters. Its eigenvalues must 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. + 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. - 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. + 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 ------- @@ -300,27 +316,56 @@ def discrete_var(A, 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. + + * `mc.P`: A 2-dim array containing the transition probability + matrix over the discretized states. + * `mc.state_values`: A 2-dim array containing the state vectors + (of dimension m) as rows, which are ordered according to the + `order` option. - 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) + 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]]) + """ A = np.asarray(A) C = np.asarray(C) @@ -356,4 +401,4 @@ def discrete_var(A, # Fit the Markov chain mc = fit_discrete_mc(X.T, V, order=order) - return mc \ No newline at end of file + return mc From ac83c52f0081be896a9633cd94d3e9862f5c60b5 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Wed, 25 Jan 2023 18:09:55 +0900 Subject: [PATCH 17/25] Import _lss and _matrix_eqn --- quantecon/markov/approximation.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 7f5eca242..7e497c237 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -7,8 +7,8 @@ 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 .._lss import simulate_linear_model +from .._matrix_eqn import solve_discrete_lyapunov from ..util import check_random_state import warnings @@ -277,7 +277,7 @@ def discrete_var(A, 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) @@ -309,7 +309,7 @@ def discrete_var(A, instance to set the initial state of the random number generator for reproducibility. If None, a randomly initialized RandomState is used. - + Returns ------- mc : MarkovChain @@ -322,7 +322,7 @@ def discrete_var(A, * `mc.state_values`: 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 From f8e9551c91eff13530fdd714969ae4fa05cc3cc3 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Thu, 16 Feb 2023 01:07:43 +1100 Subject: [PATCH 18/25] add new tests and bug fixing --- quantecon/markov/approximation.py | 7 +- quantecon/markov/tests/test_approximation.py | 69 ++++++++++++++++---- 2 files changed, 60 insertions(+), 16 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 6f6f0b8f8..723d0d8e2 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -291,6 +291,8 @@ def discrete_var(A, 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. @@ -332,7 +334,7 @@ def discrete_var(A, 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) + u = rv.rvs(size=(sim_length-1, r), random_state=random_state) v = C @ u.T x0 = np.zeros(m) @@ -352,8 +354,7 @@ def discrete_var(A, 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 \ No newline at end of file diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 453e0bb1c..7f4a59a60 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -128,18 +128,41 @@ def setup_method(self): [ 0.38556417, 0. ], [ 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.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)) @@ -152,4 +175,24 @@ def test_discretization(self): 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) \ No newline at end of file + 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) From b7f4497860c1b6f18871ef568066cd5d13e802b2 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Thu, 16 Feb 2023 12:00:08 +1100 Subject: [PATCH 19/25] add tests for non-square C --- quantecon/markov/tests/test_approximation.py | 54 +++++++++++++++----- 1 file changed, 42 insertions(+), 12 deletions(-) diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 7f4a59a60..7368e9660 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -113,15 +113,16 @@ def setup_method(self): self.random_state = np.random.RandomState(0) Omega = np.array([[0.0012346, -0.0000776], - [-0.0000776, 0.0000401]]) + [-0.0000776, 0.0000401]]) self.A = [[0.7901, -1.3570], - [-0.0104, 0.8638]] + [-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], + self.S_out = [ + [-0.38556417, -0.05387746], [-0.38556417, 0. ], [-0.38556417, 0.05387746], [ 0.38556417, -0.05387746], @@ -129,12 +130,12 @@ def setup_method(self): [ 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]] + [-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, @@ -163,9 +164,26 @@ def setup_method(self): 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)) + + 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 @@ -196,3 +214,15 @@ def test_order_F(self): 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) From 15c4f78206e8c7e9630202af6b88e714818ceec7 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Thu, 16 Feb 2023 12:07:14 +1100 Subject: [PATCH 20/25] adjust docstring --- quantecon/markov/approximation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 230837706..719f43c8e 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -282,7 +282,8 @@ def discrete_var(A, ---------- A : array_like(float, ndim=2) An m x m matrix containing the process' autocorrelation - parameters. Its eigenvalues must have moduli bounded by unity. + 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) From 496942dff0d9c670d382fff8183c94588c707055 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 16 Feb 2023 21:59:45 +0900 Subject: [PATCH 21/25] Minor edits in docstring --- quantecon/markov/approximation.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 719f43c8e..876d0fc6a 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -315,14 +315,17 @@ def discrete_var(A, ------- 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: - - * `mc.P`: A 2-dim array containing the transition probability - matrix over the discretized states. - * `mc.state_values`: A 2-dim array containing the state vectors - (of dimension m) as rows, which are ordered according to the - `order` option. + 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 -------- From 4392ef944aa2fe049800a6f372a1e248f6a49b64 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Fri, 17 Feb 2023 00:58:07 +1100 Subject: [PATCH 22/25] adjust based on review --- quantecon/markov/approximation.py | 9 +++++---- quantecon/markov/tests/test_approximation.py | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 719f43c8e..8c80f8561 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -339,7 +339,9 @@ def discrete_var(A, ... [-0.0000776, 0.0000401]] >>> C = scipy.linalg.sqrtm(Omega) >>> grid_sizes = [21, 11] - >>> mc = discrete_var(A, C, grid_sizes, random_state=rng) + >>> mc = discrete_var(A, C, grid_sizes, + ... rv=scipy.stats.multivariate_normal(cov=np.identity(2)), + ... random_state=rng) >>> mc.P.shape (145, 145) >>> mc.state_values.shape @@ -366,7 +368,6 @@ def discrete_var(A, [ 0.15422567, -0.03232648], [ 0.15422567, -0.03232648], [ 0.19278209, -0.03232648]]) - """ A = np.asarray(A) C = np.asarray(C) @@ -378,8 +379,8 @@ def discrete_var(A, 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) - + 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) diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py index 7368e9660..ed3268683 100644 --- a/quantecon/markov/tests/test_approximation.py +++ b/quantecon/markov/tests/test_approximation.py @@ -129,7 +129,7 @@ def setup_method(self): [ 0.38556417, 0. ], [ 0.38556417, 0.05387746]] - self.S_out_orderF =[ + self.S_out_orderF = [ [-0.38556417, -0.05387746], [ 0.38556417, -0.05387746], [-0.38556417, 0. ], @@ -200,7 +200,7 @@ def test_sp_distributions(self): self.A, self.C, grid_sizes=self.sizes, sim_length=self.T, - rv=sp.stats.multivariate_normal, + 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) From 0a8a7f698bf5832b207fa92c577ae28a9d0455f2 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Fri, 17 Feb 2023 10:35:17 +1100 Subject: [PATCH 23/25] update docstring --- quantecon/markov/approximation.py | 40 +++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 8c80f8561..2af95d59a 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -339,9 +339,7 @@ def discrete_var(A, ... [-0.0000776, 0.0000401]] >>> C = scipy.linalg.sqrtm(Omega) >>> grid_sizes = [21, 11] - >>> mc = discrete_var(A, C, grid_sizes, - ... rv=scipy.stats.multivariate_normal(cov=np.identity(2)), - ... random_state=rng) + >>> mc = discrete_var(A, C, grid_sizes, random_state=rng) >>> mc.P.shape (145, 145) >>> mc.state_values.shape @@ -368,6 +366,39 @@ def discrete_var(A, [ 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 + + >>> mc = discrete_var(A, C, grid_sizes, + ... rv=scipy.stats.multivariate_t(shape=np.identity(2)), + ... random_state=rng) + >>> mc.P.shape + (231, 231) + >>> mc.state_values.shape + (231, 2) + >>> mc.state_values[:10] + array([[-0.38556417, -0.05387746], + [-0.38556417, -0.04310197], + [-0.38556417, -0.03232648], + [-0.38556417, -0.02155098], + [-0.38556417, -0.01077549], + [-0.38556417, 0. ], + [-0.38556417, 0.01077549], + [-0.38556417, 0.02155098], + [-0.38556417, 0.03232648], + [-0.38556417, 0.04310197]]) + >>> mc.simulate(10, random_state=rng) + array([[ 0.2313385 , -0.02155098], + [ 0.2313385 , -0.02155098], + [ 0.11566925, -0.02155098], + [ 0.15422567, -0.02155098], + [ 0.11566925, 0. ], + [ 0. , 0. ], + [ 0.03855642, 0.01077549], + [ 0. , 0.02155098], + [ 0.03855642, 0.03232648], + [-0.03855642, 0.02155098]]) """ A = np.asarray(A) C = np.asarray(C) @@ -380,7 +411,7 @@ def discrete_var(A, 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) @@ -399,6 +430,7 @@ def discrete_var(A, 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) From a85a6cd626082bc1c36e7f42f050139e592ece36 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Fri, 17 Feb 2023 12:17:58 +1100 Subject: [PATCH 24/25] set unit Cov for t-distribution --- quantecon/markov/approximation.py | 48 ++++++++++++++++--------------- 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py index 2af95d59a..ce7f19967 100644 --- a/quantecon/markov/approximation.py +++ b/quantecon/markov/approximation.py @@ -370,35 +370,37 @@ def discrete_var(A, 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=np.identity(2)), - ... random_state=rng) + ... rv=scipy.stats.multivariate_t(shape=Sigma, df=df), + ... random_state=rng) >>> mc.P.shape - (231, 231) + (146, 146) >>> mc.state_values.shape - (231, 2) + (146, 2) >>> mc.state_values[:10] - array([[-0.38556417, -0.05387746], - [-0.38556417, -0.04310197], - [-0.38556417, -0.03232648], - [-0.38556417, -0.02155098], - [-0.38556417, -0.01077549], - [-0.38556417, 0. ], - [-0.38556417, 0.01077549], - [-0.38556417, 0.02155098], + array([[-0.38556417, 0.02155098], [-0.38556417, 0.03232648], - [-0.38556417, 0.04310197]]) + [-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.2313385 , -0.02155098], - [ 0.2313385 , -0.02155098], - [ 0.11566925, -0.02155098], - [ 0.15422567, -0.02155098], - [ 0.11566925, 0. ], - [ 0. , 0. ], - [ 0.03855642, 0.01077549], - [ 0. , 0.02155098], - [ 0.03855642, 0.03232648], - [-0.03855642, 0.02155098]]) + 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) From 08dd3214fb845bf2580c6d8c33a0aee2a600b38f Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Wed, 22 Feb 2023 11:12:57 +0900 Subject: [PATCH 25/25] CI: Workaround for windows (#694) --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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