Merge pull request #1484 from econ-ark/ImproveConsMarkovConstructors

Add IncShkDstn constructor for ConsMarkov

Documentation/CHANGELOG.md

@@ -26,6 +26,7 @@ Release Date: TBD
 - Fixes bug in `AgentPopulation` that caused discretization of distributions to not work. [1275](https://github.com/econ-ark/HARK/pull/1275)
 - Adds support for distributions, booleans, and callables as parameters in the `Parameters` class. [1387](https://github.com/econ-ark/HARK/pull/1387)
 - Removes a specific way of accounting for ``employment'' in the idiosyncratic-shocks income process. [1473](https://github.com/econ-ark/HARK/pull/1473)
+- Adds income process constructor for the discrete Markov state consumption-saving model. [1484](https://github.com/econ-ark/HARK/pull/1484)
 - Changes the behavior of make_lognormal_RiskyDstn so that the standard deviation represents the standard deviation of log(returns)
 - Adds detailed parameter and latex documentation to most models.
 - Add PermGroFac constructor that explicitly combines idiosyncratic and aggregate sources of growth. [1489](https://github.com/econ-ark/HARK/pull/1489)

HARK/Calibration/Income/IncomeProcesses.py

@@ -462,6 +462,189 @@ def construct_lognormal_income_process_unemployment(
     return IncShkDstn
 
 
+def construct_markov_lognormal_income_process_unemployment(
+    T_cycle,
+    PermShkStd,
+    PermShkCount,
+    TranShkStd,
+    TranShkCount,
+    T_retire,
+    UnempPrb,
+    IncUnemp,
+    UnempPrbRet,
+    IncUnempRet,
+    RNG,
+    neutral_measure=False,
+):
+    """
+    Generates a nested list of discrete approximations to the income process for each
+    life period, from end of life to beginning of life, for each discrete Markov
+    state in the problem. This function calls construct_lognormal_income_process_unemployment
+    for each Markov state, then rearranges the output. Permanent shocks are mean
+    one lognormally distributed with standard deviation PermShkStd[t] during the
+    working life, and degenerate at 1 in the retirement period.  Transitory shocks
+    are mean one lognormally distributed with a point mass at IncUnemp with
+    probability UnempPrb while working; they are mean one with a point mass at
+    IncUnempRet with probability UnempPrbRet.  Retirement occurs after t=T_retire
+    periods of working. The problem is specified as having T=T_cycle periods and
+    K discrete Markov states.
+
+    Parameters
+    ----------
+    PermShkStd : np.array
+        2D array of shape (T,K) of standard deviations of log permanent income
+        uncertainty during the agent's life.
+    PermShkCount : int
+        The number of approximation points to be used in the discrete approxima-
+        tion to the permanent income shock distribution.
+    TranShkStd : np.array
+        2D array of shape (T,K) standard deviations in log transitory income
+        uncertainty during the agent's life.
+    TranShkCount : int
+        The number of approximation points to be used in the discrete approxima-
+        tion to the permanent income shock distribution.
+    UnempPrb : np.array
+        1D or 2D array of the probability of becoming unemployed during the working
+        portion of the agent's life. If 1D, it should be size K, providing one
+        unemployment probability per discrete Markov state. If 2D, it should be
+        shape (T,K).
+    UnempPrbRet : np.array or None
+        The probability of not receiving typical retirement income when retired.
+        If not None, should be size K. Can be set to None when T_retire is 0.
+    T_retire : int
+        The index value for the final working period in the agent's life.
+        If T_retire <= 0 then there is no retirement.
+    IncUnemp : np.array
+        1D or 2D array of transitory income received when unemployed during the
+        working portion of the agent's life. If 1D, it should be size K, providing
+        one unemployment probability per discrete Markov state. If 2D, it should be
+        shape (T,K).
+    IncUnempRet : np.array or None
+        Transitory income received while "unemployed" when retired. If provided,
+        should be size K. Can be None when T_retire is 0.
+    T_cycle :  int
+        Total number of non-terminal periods in the consumer's sequence of periods.
+    RNG : np.random.RandomState
+        Random number generator for this type.
+    neutral_measure : bool
+        Indicator for whether the permanent-income-neutral measure should be used.
+
+    Returns
+    -------
+    IncShkDstn : [[distribution.Distribution]]
+        A list with T_cycle elements, each of which is a discrete approximation
+        to the income process in a period.
+    """
+    if T_retire > 0:
+        normal_length = T_retire
+        retire_length = T_cycle - T_retire
+    else:
+        normal_length = T_cycle
+        retire_length = 0
+
+    # Check dimensions of inputs
+    try:
+        PermShkStd_K = PermShkStd.shape[1]
+        TranShkStd_K = TranShkStd.shape[1]
+        if UnempPrb.ndim == 2:
+            UnempPrb_K = UnempPrb.shape[1]
+        else:
+            UnempPrb_K = UnempPrb.shape[0]
+        if IncUnemp.ndim == 2:
+            IncUnemp_K = IncUnemp.shape[1]
+        else:
+            IncUnemp_K = IncUnemp.shape[0]
+        K = PermShkStd_K
+        assert K == TranShkStd_K
+        assert K == TranShkStd_K
+        assert K == UnempPrb_K
+        assert K == IncUnemp_K
+    except:
+        raise Exception(
+            "The last dimension of PermShkStd, TranShkStd, IncUnemp,"
+            + " and UnempPrb must all be K, the number of discrete states."
+        )
+    try:
+        if T_retire > 0:
+            assert K == UnempPrbRet.size
+            assert K == IncUnempRet.size
+    except:
+        raise Exception(
+            "When T_retire is not zero, UnempPrbRet and IncUnempRet"
+            + " must be specified as arrays of size K, the number of "
+            + "discrete Markov states."
+        )
+    try:
+        D = UnempPrb.ndim
+        assert D == IncUnemp.ndim
+        if T_retire > 0:
+            assert D == UnempPrbRet.ndim
+            assert D == IncUnempRet.ndim
+    except:
+        raise Exception(
+            "If any of UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet "
+            + "are 2D arrays, then they must *all* be 2D arrays."
+        )
+    try:
+        assert D == 1 or D == 2
+    except:
+        raise Exception(
+            "UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet must "
+            + "all be 1D or 2D arrays."
+        )
+
+    # Prepare lists that don't vary by Markov state
+    PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
+    TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
+    neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
+
+    # Loop through the Markov states, constructing the lifecycle income process for each one
+    IncShkDstn_by_Mrkv = []
+    for k in range(K):
+        if D == 1:  # Unemployment parameters don't vary by age other than retirement
+            if T_retire > 0:
+                UnempPrb_list = [UnempPrb[k]] * normal_length + [
+                    UnempPrbRet[k]
+                ] * retire_length
+                IncUnemp_list = [IncUnemp[k]] * normal_length + [
+                    IncUnempRet[k]
+                ] * retire_length
+            else:
+                UnempPrb_list = [UnempPrb[k]] * normal_length
+                IncUnemp_list = [IncUnemp[k]] * normal_length
+        else:  # Unemployment parameters vary by age
+            UnempPrb_list = UnempPrb[:, k].tolist()
+            IncUnemp_list = IncUnemp[:, k].tolist()
+
+        PermShkStd_k = PermShkStd[:, k].tolist()
+        TranShkStd_k = TranShkStd[:, k].tolist()
+
+        IncShkDstn_k = IndexDistribution(
+            engine=BufferStockIncShkDstn,
+            conditional={
+                "sigma_Perm": PermShkStd_k,
+                "sigma_Tran": TranShkStd_k,
+                "n_approx_Perm": PermShkCount_list,
+                "n_approx_Tran": TranShkCount_list,
+                "neutral_measure": neutral_measure_list,
+                "UnempPrb": UnempPrb_list,
+                "IncUnemp": IncUnemp_list,
+            },
+            RNG=RNG,
+            seed=RNG.integers(0, 2**31 - 1),
+        )
+        IncShkDstn_by_Mrkv.append(IncShkDstn_k)
+
+    # Rearrange the list that was just constructed, so that element [t][k] represents
+    # the income shock distribution in period t, Markov state k.
+    IncShkDstn = []
+    for t in range(T_cycle):
+        IncShkDstn_t = [IncShkDstn_by_Mrkv[k][t] for k in range(K)]
+        IncShkDstn.append(IncShkDstn_t)
+
+    return IncShkDstn
+
+
 def construct_HANK_lognormal_income_process_unemployment(
     T_cycle,
     PermShkStd,
@@ -613,6 +796,28 @@ def get_TranShkDstn_from_IncShkDstn(IncShkDstn, RNG):
     return TimeVaryingDiscreteDistribution(TranShkDstn, seed=RNG.integers(0, 2**31 - 1))
 
 
+def get_PermShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
+    PermShkDstn = [
+        [
+            this.make_univariate(0, seed=RNG.integers(0, 2**31 - 1))
+            for this in IncShkDstn_t
+        ]
+        for IncShkDstn_t in IncShkDstn
+    ]
+    return PermShkDstn
+
+
+def get_TranShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
+    TranShkDstn = [
+        [
+            this.make_univariate(1, seed=RNG.integers(0, 2**31 - 1))
+            for this in IncShkDstn_t
+        ]
+        for IncShkDstn_t in IncShkDstn
+    ]
+    return TranShkDstn
+
+
 def get_TranShkGrid_from_TranShkDstn(T_cycle, TranShkDstn):
     TranShkGrid = [TranShkDstn[t].atoms.flatten() for t in range(T_cycle)]
     return TranShkGrid

HARK/ConsumptionSaving/ConsMarkovModel.py

@@ -9,9 +9,9 @@
 
 from HARK import AgentType, NullFunc
 from HARK.Calibration.Income.IncomeProcesses import (
-    construct_lognormal_income_process_unemployment,
-    get_PermShkDstn_from_IncShkDstn,
-    get_TranShkDstn_from_IncShkDstn,
+    construct_markov_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn_markov,
+    get_TranShkDstn_from_IncShkDstn_markov,
 )
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
@@ -657,25 +657,31 @@ def calc_vPPnext(S, a, R):
 
 # Make a dictionary of constructors for the markov consumption-saving model
 markov_constructor_dict = {
-    "IncShkDstn": construct_lognormal_income_process_unemployment,
-    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
-    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "IncShkDstn": construct_markov_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn_markov,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn_markov,
     "aXtraGrid": make_assets_grid,
     "MrkvArray": make_simple_binary_markov,
     "solution_terminal": make_markov_solution_terminal,
 }
 
 # Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
 default_IncShkDstn_params = {
-    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkStd": np.array(
+        [[0.1, 0.1]]
+    ),  # Standard deviation of log permanent income shocks
     "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
-    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkStd": np.array(
+        [[0.1, 0.1]]
+    ),  # Standard deviation of log transitory income shocks
     "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
-    "UnempPrb": 0.05,  # Probability of unemployment while working
-    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "UnempPrb": np.array([0.05, 0.05]),  # Probability of unemployment while working
+    "IncUnemp": np.array(
+        [0.3, 0.3]
+    ),  # Unemployment benefits replacement rate while working
     "T_retire": 0,  # Period of retirement (0 --> no retirement)
-    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
-    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+    "UnempPrbRet": None,  # Probability of "unemployment" while retired
+    "IncUnempRet": None,  # "Unemployment" benefits when retired
 }
 
 # Default parameters to make aXtraGrid using make_assets_grid

HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py

@@ -1,5 +1,5 @@
 import unittest
-from copy import copy
+from copy import deepcopy
 
 import numpy as np
 
@@ -56,19 +56,19 @@ def setUp(self):
             ]
         )
 
-        init_serial_unemployment = copy(init_indshk_markov)
+        init_serial_unemployment = deepcopy(init_indshk_markov)
         init_serial_unemployment["MrkvArray"] = [MrkvArray]
-        init_serial_unemployment["UnempPrb"] = 0.0
-        # to make income distribution when employed
+        init_serial_unemployment["UnempPrb"] = np.zeros(2)
+        # Income process is overwritten below to make income distribution when employed
         init_serial_unemployment["global_markov"] = False
         init_serial_unemployment["Rfree"] = np.array([1.03, 1.03, 1.03, 1.03])
         init_serial_unemployment["LivPrb"] = [np.array([0.98, 0.98, 0.98, 0.98])]
         init_serial_unemployment["PermGroFac"] = [np.array([1.01, 1.01, 1.01, 1.01])]
+        init_serial_unemployment["constructors"]["MrkvArray"] = None
 
         self.model = MarkovConsumerType(**init_serial_unemployment)
         self.model.cycles = 0
         self.model.vFuncBool = False  # for easy toggling here
-        self.model.MrkvArray = [MrkvArray]
 
         # Replace the default (lognormal) income distribution with a custom one
         employed_income_dist = DiscreteDistributionLabeled(
@@ -224,3 +224,8 @@ def main_test(self):
         self.assertAlmostEqual(
             Markov_vFuncBool_example.solution[0].vFunc[1](0.4), -4.12794
         )
+
+
+if __name__ == "__main__":
+    # Run all the tests
+    unittest.main()

HARK/ConsumptionSaving/tests/test_TractableBufferStockModel.py

@@ -24,3 +24,8 @@ def test_simulation(self):
             self.tct.history["mLvl"][15][0] - self.tct.history["cLvlNow"][15][0],
             self.tct.history["aLvl"][15][0],
         )
+
+
+if __name__ == "__main__":
+    # Run all the tests
+    unittest.main()

HARK/ConsumptionSaving/tests/test_modelInits.py

@@ -88,10 +88,14 @@ def test_MarkovConsumerType(self):
             # Make a consumer with serially correlated unemployment, subject to boom and bust cycles
             init_serial_unemployment = {}
             init_serial_unemployment["MrkvArray"] = [MrkvArray]
-            init_serial_unemployment["UnempPrb"] = (
-                0.0  # to make income distribution when employed
-            )
+            init_serial_unemployment["UnempPrb"] = np.zeros(2)
+            # Income process is overwritten below to make income distribution when employed
             init_serial_unemployment["global_markov"] = False
+            init_serial_unemployment["Rfree"] = np.array([1.03, 1.03, 1.03, 1.03])
+            init_serial_unemployment["LivPrb"] = [np.array([0.98, 0.98, 0.98, 0.98])]
+            init_serial_unemployment["PermGroFac"] = [
+                np.array([1.01, 1.01, 1.01, 1.01])
+            ]
             SerialUnemploymentExample = MarkovConsumerType(**init_serial_unemployment)
         except:
             self.fail(

HARK/ConsumptionSaving/tests/test_modelcomparisons.py

@@ -123,14 +123,14 @@ def setUp(self):
             "Mrkv_p11": [1.0 - base_primitives["UnempPrb"]],
             "Mrkv_p22": [1.0],
             "BoroCnstArt": None,
-            "PermShkStd": [0.0],
+            "PermShkStd": np.array([[0.0, 0.0]]),
             "PermShkCount": 1,
-            "TranShkStd": [0.0],
+            "TranShkStd": np.array([[0.0, 0.0]]),
             "TranShkCount": 1,
-            "UnempPrb": 0.0,
+            "UnempPrb": np.array([[0.0, 0.0]]),  # This will be overwritten
             "UnempPrbRet": 0.0,
             "T_retire": 0,
-            "IncUnemp": 0.0,
+            "IncUnemp": np.array([[0.0, 0.0]]),  # This will be overwritten
             "IncUnempRet": 0.0,
             "aXtraMin": 0.001,
             "aXtraMax": TBSType.mUpperBnd,
@@ -149,6 +149,8 @@ def setUp(self):
             "vFuncBool": False,
             "CubicBool": True,
             "T_cycle": 1,
+            "MrkvArray": [np.eye(2)],
+            # Will be overwritten, might prevent glitch in Ubuntu
         }
 
         MarkovType = MarkovConsumerType(**Markov_primitives)

HARK/core.py

@@ -601,6 +601,12 @@ def construct(self, *args, force=False):
                     else:
                         raise ValueError("No constructor found for " + key) from None
 
+                # If this constructor is None, do nothing and mark it as completed
+                if constructor is None:
+                    keys_complete[i] = True
+                    anything_accomplished_this_pass = True  # We did something!
+                    continue
+
                 # Get the names of arguments for this constructor and try to gather them
                 args_needed = get_arg_names(constructor)
                 has_no_default = {