Merge pull request #1521 from QingyuanFang/master

Add an alternative way of implementing Tauchen's method to make_tauchen_ar1()

Author: mnwhite

HARK/distributions/__init__.py

@@ -22,6 +22,7 @@
     "Uniform",
     "MarkovProcess",
     "add_discrete_outcome_constant_mean",
+    "make_tauchen_ar1",
 ]
 
 from HARK.distributions.base import (
@@ -53,4 +54,5 @@
     combine_indep_dstns,
     distr_of_function,
     expected,
+    make_tauchen_ar1,
 )

HARK/distributions/utils.py

@@ -180,10 +180,10 @@ def make_markov_approx_to_normal_by_monte_carlo(x_grid, mu, sigma, N_draws=10000
     return p_vec
 
 
-def make_tauchen_ar1(N, sigma=1.0, ar_1=0.9, bound=3.0):
+def make_tauchen_ar1(N, sigma=1.0, ar_1=0.9, bound=3.0, inflendpoint=True):
     """
     Function to return a discretized version of an AR1 process.
-    See https://www.fperri.net/TEACHING/macrotheory08/numerical.pdf for details
+    See http://www.fperri.net/TEACHING/macrotheory08/numerical.pdf for details
 
     Parameters
     ----------
@@ -196,6 +196,12 @@ def make_tauchen_ar1(N, sigma=1.0, ar_1=0.9, bound=3.0):
     bound: float
         The highest (lowest) grid point will be bound (-bound) multiplied by the unconditional
         standard deviation of the process
+    inflendpoint: Bool
+        If True: implement the standard method as in Tauchen (1986):
+            assign the probability of jumping to a point outside the grid to the closest endpoint
+        If False: implement an alternative method:
+            discard the probability of jumping to a point outside the grid, effectively
+            reassigning it to the remaining points in proportion to their probability of being reached
 
     Returns
     -------
@@ -208,16 +214,25 @@ def make_tauchen_ar1(N, sigma=1.0, ar_1=0.9, bound=3.0):
     y = np.linspace(-yN, yN, N)
     d = y[1] - y[0]
     trans_matrix = np.ones((N, N))
-    for j in range(N):
-        for k_1 in range(N - 2):
-            k = k_1 + 1
-            trans_matrix[j, k] = stats.norm.cdf(
-                (y[k] + d / 2.0 - ar_1 * y[j]) / sigma
-            ) - stats.norm.cdf((y[k] - d / 2.0 - ar_1 * y[j]) / sigma)
-        trans_matrix[j, 0] = stats.norm.cdf((y[0] + d / 2.0 - ar_1 * y[j]) / sigma)
-        trans_matrix[j, N - 1] = 1.0 - stats.norm.cdf(
-            (y[N - 1] - d / 2.0 - ar_1 * y[j]) / sigma
-        )
+    if inflendpoint:
+        for j in range(N):
+            for k_1 in range(N - 2):
+                k = k_1 + 1
+                trans_matrix[j, k] = stats.norm.cdf(
+                    (y[k] + d / 2.0 - ar_1 * y[j]) / sigma
+                ) - stats.norm.cdf((y[k] - d / 2.0 - ar_1 * y[j]) / sigma)
+            trans_matrix[j, 0] = stats.norm.cdf((y[0] + d / 2.0 - ar_1 * y[j]) / sigma)
+            trans_matrix[j, N - 1] = 1.0 - stats.norm.cdf(
+                (y[N - 1] - d / 2.0 - ar_1 * y[j]) / sigma
+            )
+    else:
+        for j in range(N):
+            for k in range(N):
+                trans_matrix[j, k] = stats.norm.cdf(
+                    (y[k] + d / 2.0 - ar_1 * y[j]) / sigma
+                ) - stats.norm.cdf((y[k] - d / 2.0 - ar_1 * y[j]) / sigma)
+        ## normalize: each row sums to 1
+        trans_matrix = trans_matrix / trans_matrix.sum(axis=1)[:, np.newaxis]
 
     return y, trans_matrix
 

HARK/tests/test_distribution.py

@@ -22,6 +22,7 @@
     combine_indep_dstns,
     distr_of_function,
     expected,
+    make_tauchen_ar1,
 )
 from HARK.tests import HARK_PRECISION
 
@@ -657,3 +658,38 @@ def transition(shocks, state):
 
         assert np.all(exp1["m"] == exp2["m"]).item()
         assert np.all(exp1["n"] == exp2["n"]).item()
+
+
+class TestTauchenAR1(unittest.TestCase):
+    def test_tauchen(self):
+        # Test with a simple AR(1) process
+        N = 5
+        sigma = 1.0
+        ar_1 = 0.9
+        bound = 3.0
+
+        # By default, inflendpoint = True
+        standard = make_tauchen_ar1(N, sigma, ar_1, bound)
+        alternative = make_tauchen_ar1(N, sigma, ar_1, bound, inflendpoint=False)
+
+        # Check that the grid points of the two methods are identical
+        self.assertTrue(np.all(np.equal(standard[0], alternative[0])))
+
+        # Check the shape of the transition matrix
+        self.assertEqual(standard[1].shape, (N, N))
+        self.assertEqual(alternative[1].shape, (N, N))
+
+        # Check that the sum of each row in the transition matrix is 1
+        self.assertTrue(np.allclose(np.sum(standard[1], axis=1), np.ones(N)))
+        self.assertTrue(np.allclose(np.sum(alternative[1], axis=1), np.ones(N)))
+
+        # Check that [k]-th column ./ [k-1]-th column are identical (k = 3, ..., N-1)
+        # Note: the first and the last column of the 'standard' transition matrix are inflated
+        if N > 3:
+            for i in range(2, N - 1):
+                self.assertTrue(
+                    np.allclose(
+                        standard[1][:, i] * alternative[1][:, i - 1],
+                        alternative[1][:, i] * standard[1][:, i - 1],
+                    )
+                )