[WIP] t-indexed inputs for lifecycle models

Documentation/CHANGELOG.md

@@ -14,6 +14,7 @@ Release Data: TBD
 
 #### Major Changes
 
+* Input parameters for cyclical models now indexed by t [#1039](https://github.com/econ-ark/HARK/pull/1039)
 * A IndexDistribution class for representing time-indexed probability distributions [#1018](https://github.com/econ-ark/pull/1018/).
 * Adds new consumption-savings-portfolio model `RiskyContrib`, which represents an agent who can save in risky and risk-free assets but faces
 frictions to moving funds between them. To circumvent these frictions, he has access to an income-deduction scheme to accumulate risky assets.

HARK/Calibration/Income/IncomeTools.py

@@ -629,26 +629,20 @@ def parse_income_spec(
         raise NotImplementedError()
 
     # Volatilities
-    # In this section, it is important to keep in mind that IncomeDstn[t]
-    # is the income distribution from period t to t+1, as perceived in period
-    # t.
-    # Therefore (assuming an annual model with agents entering at age 0),
-    # IncomeDstn[3] would contain the distribution of income shocks that occur
-    # at the start of age 4.
     if SabelhausSong:
 
         if age_ret is None:
 
             IncShkStds = sabelhaus_song_var_profile(
-                cohort=1950, age_min=age_min + 1, age_max=age_max
+                cohort=1950, age_min=age_min, age_max=age_max
             )
             PermShkStd = IncShkStds["PermShkStd"]
             TranShkStd = IncShkStds["TranShkStd"]
 
         else:
 
             IncShkStds = sabelhaus_song_var_profile(
-                cohort=1950, age_min=age_min + 1, age_max=age_ret
+                cohort=1950, age_min=age_min, age_max=age_ret
             )
             PermShkStd = IncShkStds["PermShkStd"] + [0.0] * (N_ret_periods + 1)
             TranShkStd = IncShkStds["TranShkStd"] + [0.0] * (N_ret_periods + 1)
@@ -664,10 +658,10 @@ def parse_income_spec(
 
             else:
 
-                PermShkStd = [PermShkStd] * (N_work_periods - 1) + [0.0] * (
+                PermShkStd = [PermShkStd] * (N_work_periods) + [0.0] * (
                     N_ret_periods + 1
                 )
-                TranShkStd = [TranShkStd] * (N_work_periods - 1) + [0.0] * (
+                TranShkStd = [TranShkStd] * (N_work_periods) + [0.0] * (
                     N_ret_periods + 1
                 )
 

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -60,6 +60,7 @@
     "IndShockConsumerType",
     "KinkedRconsumerType",
     "init_perfect_foresight",
+    "init_perfect_foresight_infinite",
     "init_idiosyncratic_shocks",
     "init_kinked_R",
     "init_lifecycle",
@@ -1536,8 +1537,8 @@ def prepare_to_calc_EndOfPrdvP(self):
     'CRRA': 2.0,          # Coefficient of relative risk aversion,
     'Rfree': 1.03,        # Interest factor on assets
     'DiscFac': 0.96,      # Intertemporal discount factor
-    'LivPrb': [0.98],     # Survival probability
-    'PermGroFac': [1.01],  # Permanent income growth factor
+    'LivPrb': [0.98, 0.98],     # Survival probability
+    'PermGroFac': [1.01, 1.01],  # Permanent income growth factor
     'BoroCnstArt': None,  # Artificial borrowing constraint
     'MaxKinks': 400,      # Maximum number of grid points to allow in cFunc (should be large)
     'AgentCount': 10000,  # Number of agents of this type (only matters for simulation)
@@ -1549,9 +1550,16 @@ def prepare_to_calc_EndOfPrdvP(self):
     # Aggregate permanent income growth factor: portion of PermGroFac attributable to aggregate productivity growth (only matters for simulation)
     'PermGroFacAgg': 1.0,
     'T_age': None,       # Age after which simulated agents are automatically killed
-    'T_cycle': 1         # Number of periods in the cycle for this agent type
+    'T_cycle': 2         # Number of periods in the cycle for this agent type
 }
 
+init_perfect_foresight_infinite = init_perfect_foresight.copy()
+init_perfect_foresight_infinite.update({
+    'cycles' : 0,         # Finite, non-cyclic model
+    'LivPrb': [0.98],     # Survival probability
+    'PermGroFac': [1.01],
+    'T_cycle': 1  
+})
 
 class PerfForesightConsumerType(AgentType):
     """
@@ -1729,7 +1737,7 @@ def sim_death(self):
         # Determine who dies
         DiePrb_by_t_cycle = 1.0 - np.asarray(self.LivPrb)
         DiePrb = DiePrb_by_t_cycle[
-            self.t_cycle - 1
+            self.t_cycle - 1 if self.cycles == 1 else self.t_cycle
         ]  # Time has already advanced, so look back one
 
         # In finite-horizon problems the previous line gives newborns the
@@ -2019,9 +2027,9 @@ def check_conditions(self, verbose=None):
             None
         ],  # Some other value of "assets above minimum" to add to the grid, not used
         # Income process variables
-        "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+        "PermShkStd": [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": [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
         "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
@@ -2173,12 +2181,13 @@ def get_shocks(self):
         newborn = self.t_age == 0
         for t in range(self.T_cycle):
             these = t == self.t_cycle
+
             N = np.sum(these)
             if N > 0:
                 IncShkDstnNow = self.IncShkDstn[
-                    t - 1
+                    t
                 ]  # set current income distribution
-                PermGroFacNow = self.PermGroFac[t - 1]  # and permanent growth factor
+                PermGroFacNow = self.PermGroFac[t]  # and permanent growth factor
                 # Get random draws of income shocks from the discrete distribution
                 IncShks = IncShkDstnNow.draw(N)
 
@@ -2187,8 +2196,7 @@ def get_shocks(self):
                 )  # permanent "shock" includes expected growth
                 TranShkNow[these] = IncShks[1, :]
 
-        # That procedure used the *last* period in the sequence for newborns, but that's not right
-        # Redraw shocks for newborns, using the *first* period in the sequence.  Approximation.
+        # This is now redundant and can be safely removed. #1022
         N = np.sum(newborn)
         if N > 0:
             these = newborn
@@ -3047,7 +3055,7 @@ def construct_assets_grid(parameters):
 
 # Make a dictionary to specify an infinite consumer with a four period cycle
 init_cyclical = copy(init_idiosyncratic_shocks)
-init_cyclical['PermGroFac'] = [1.082251, 2.8, 0.3, 1.1]
+init_cyclical['PermGroFac'] = [1.1, 1.082251, 2.8, 0.3]
 init_cyclical['PermShkStd'] = [0.1, 0.1, 0.1, 0.1]
 init_cyclical['TranShkStd'] = [0.1, 0.1, 0.1, 0.1]
 init_cyclical['LivPrb'] = 4*[0.98]

HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py

@@ -23,9 +23,9 @@ def test_get_shocks(self):
 
         self.agent.get_shocks()
 
-        self.assertEqual(self.agent.shocks['PermShk'][0], 1.0427376294215103)
-        self.assertAlmostEqual(self.agent.shocks['PermShk'][1], 0.9278094171517413)
-        self.assertAlmostEqual(self.agent.shocks['TranShk'][0], 0.881761797501595)
+        self.assertEqual(self.agent.shocks['PermShk'][0], 1.0050166461586711)
+        self.assertAlmostEqual(self.agent.shocks['PermShk'][1], 1.1780702264015421)
+        self.assertAlmostEqual(self.agent.shocks['TranShk'][0], 1.0704497811153597)
 
     def test_ConsIndShockSolverBasic(self):
         LifecycleExample = IndShockConsumerType(**init_lifecycle)
@@ -35,18 +35,18 @@ def test_ConsIndShockSolverBasic(self):
         # test the solution_terminal
         self.assertAlmostEqual(LifecycleExample.solution[-1].cFunc(2).tolist(), 2)
 
-        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.79429538)
-        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.79391692)
-        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.79253095)
+        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.78943688)
+        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.78934515)
+        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.78821644)
 
         self.assertAlmostEqual(
-            LifecycleExample.solution[0].cFunc(1).tolist(), 0.7506184692092213
+            LifecycleExample.solution[0].cFunc(1).tolist(), 0.7472898578766399
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[1].cFunc(1).tolist(), 0.7586358637239385
+            LifecycleExample.solution[1].cFunc(1).tolist(), 0.7551914217221962
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[2].cFunc(1).tolist(), 0.7681247572911291
+            LifecycleExample.solution[2].cFunc(1).tolist(), 0.7645965714276972
         )
 
         solver = ConsIndShockSolverBasic(
@@ -73,20 +73,21 @@ def test_ConsIndShockSolverBasic(self):
 
         EndOfPrdvP = solver.calc_EndOfPrdvP()
 
-        self.assertAlmostEqual(EndOfPrdvP[0], 6657.839372100613)
-        self.assertAlmostEqual(EndOfPrdvP[-1], 0.2606075215645896)
+        self.assertAlmostEqual(EndOfPrdvP[0], 6657.849481857674)
+        self.assertAlmostEqual(EndOfPrdvP[-1], 0.2635560321220346)
 
         solution = solver.make_basic_solution(
             EndOfPrdvP, solver.aNrmNow, solver.make_linear_cFunc
         )
         solver.add_MPC_and_human_wealth(solution)
 
-        self.assertAlmostEqual(solution.cFunc(4).tolist(), 1.0028005137373956)
+        self.assertAlmostEqual(solution.cFunc(4).tolist(), 0.9935251977248167)
 
     def test_simulated_values(self):
         self.agent.initialize_sim()
         self.agent.simulate()
 
+        ## uses simulated values -- needs simulation code update.
         self.assertAlmostEqual(self.agent.MPCnow[1], 0.5711503906043797)
 
         self.assertAlmostEqual(self.agent.state_now['aLvl'][1], 0.18438326264597635)
@@ -251,12 +252,12 @@ def test_infinite_horizon(self):
     "CRRA": 2.0,  # Coefficient of relative risk aversion
     "Rfree": 1.03,  # Interest factor on assets
     "DiscFac": 0.96,  # Intertemporal discount factor
-    "LivPrb": [0.99, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1],
-    "PermGroFac": [1.01, 1.01, 1.01, 1.02, 1.02, 1.02, 0.7, 1.0, 1.0, 1.0],
+    "LivPrb": [1.0, 0.99, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1],
+    "PermGroFac": [1.0, 1.01, 1.01, 1.01, 1.02, 1.02, 1.02, 0.7, 1.0, 1.0, 1.0],
     # Parameters that specify the income distribution over the lifecycle
-    "PermShkStd": [0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0, 0, 0],
+    "PermShkStd": [0.1, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0, 0, 0],
     "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
-    "TranShkStd": [0.3, 0.2, 0.1, 0.3, 0.2, 0.1, 0.3, 0, 0, 0],
+    "TranShkStd": [0.3, 0.3, 0.2, 0.1, 0.3, 0.2, 0.1, 0.3, 0, 0, 0],
     "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
@@ -274,7 +275,7 @@ def test_infinite_horizon(self):
     "BoroCnstArt": 0.0,  # Artificial borrowing constraint; imposed minimum level of end-of period assets
     "vFuncBool": True,  # Whether to calculate the value function during solution
     "CubicBool": False,  # Preference shocks currently only compatible with linear cFunc
-    "T_cycle": 10,  # Number of periods in the cycle for this agent type
+    "T_cycle": 11,  # Number of periods in the cycle for this agent type
     # Parameters only used in simulation
     "AgentCount": 10000,  # Number of agents of this type
     "T_sim": 120,  # Number of periods to simulate
@@ -284,26 +285,19 @@ def test_infinite_horizon(self):
     "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
     "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
     "T_age": 11,  # Age after which simulated agents are automatically killed
+    "cycles" : 1
 }
 
 
 class testIndShockConsumerTypeLifecycle(unittest.TestCase):
     def test_lifecyle(self):
         LifecycleExample = IndShockConsumerType(**LifecycleDict)
-        LifecycleExample.cycles = 1
         LifecycleExample.solve()
 
         self.assertEqual(len(LifecycleExample.solution), 11)
 
-        mMin = np.min(
-            [
-                LifecycleExample.solution[t].mNrmMin
-                for t in range(LifecycleExample.T_cycle)
-            ]
-        )
-
         self.assertAlmostEqual(
-            LifecycleExample.solution[5].cFunc(3).tolist(), 2.129983771775666
+            LifecycleExample.solution[5].cFunc(3).tolist(), 2.16741812
         )
 
 
@@ -313,7 +307,7 @@ def test_lifecyle(self):
     "Rfree": 1.03,  # Interest factor on assets
     "DiscFac": 0.96,  # Intertemporal discount factor
     "LivPrb": 4 * [0.98],  # Survival probability
-    "PermGroFac": [1.082251, 2.8, 0.3, 1.1],
+    "PermGroFac": [1.1, 1.082251, 2.8, 0.3],
     # Parameters that specify the income distribution over the lifecycle
     "PermShkStd": [0.1, 0.1, 0.1, 0.1],
     "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
@@ -358,6 +352,11 @@ def test_cyclical(self):
             CyclicalExample.solution[3].cFunc(3).tolist(), 1.5958390056965004
         )
 
+        CyclicalExample.initialize_sim()
+        CyclicalExample.simulate()
+
+        self.assertAlmostEqual(CyclicalExample.state_now['aLvl'][1], 0.41839957)
+
 # %% Tests of 'stable points'
 
 

HARK/ConsumptionSaving/tests/test_PerfForesightConsumerType.py

@@ -1,12 +1,15 @@
-from HARK.ConsumptionSaving.ConsIndShockModel import PerfForesightConsumerType
+from HARK.ConsumptionSaving.ConsIndShockModel import (
+    PerfForesightConsumerType,
+    init_perfect_foresight_infinite
+)
 import numpy as np
 import unittest
 
 
 class testPerfForesightConsumerType(unittest.TestCase):
     def setUp(self):
         self.agent = PerfForesightConsumerType()
-        self.agent_infinite = PerfForesightConsumerType(cycles=0)
+        self.agent_infinite = PerfForesightConsumerType(**init_perfect_foresight_infinite)
 
         PF_dictionary = {
             "CRRA": 2.5,

HARK/ConsumptionSaving/tests/test_PerfForesightFastConsumerType.py

@@ -1,5 +1,8 @@
 from HARK.ConsumptionSaving.ConsIndShockModelFast import PerfForesightConsumerTypeFast
-from HARK.ConsumptionSaving.ConsIndShockModel import PerfForesightConsumerType
+from HARK.ConsumptionSaving.ConsIndShockModel import (
+    PerfForesightConsumerType,
+    init_perfect_foresight_infinite
+)
 from HARK.ConsumptionSaving.tests.test_PerfForesightConsumerType import (
     testPerfForesightConsumerType,
 )
@@ -9,8 +12,8 @@ class testPerfForesightFastConsumerType(testPerfForesightConsumerType):
     def setUp(self):
         self.agent = PerfForesightConsumerTypeFast()
         self.agent_slow = PerfForesightConsumerType()
-        self.agent_infinite = PerfForesightConsumerTypeFast(cycles=0)
-        self.agent_infinite_slow = PerfForesightConsumerType(cycles=0)
+        self.agent_infinite = PerfForesightConsumerTypeFast(**init_perfect_foresight_infinite)
+        self.agent_infinite_slow = PerfForesightConsumerType(**init_perfect_foresight_infinite)
 
         PF_dictionary = {
             "CRRA": 2.5,

HARK/core.py

@@ -1004,10 +1004,13 @@ def solve_one_cycle(agent, solution_last):
     # Initialize the solution for this cycle, then iterate on periods
     solution_cycle = []
     solution_next = solution_last
-    for t in range(T):
+
+    lifecycle_range = range(T-1, 0, -1) # All but the first
+    cycle_range = [0] + list(range(T - 1, 0, -1))
+    for k in (lifecycle_range if agent.cycles == 1 else cycle_range):
         # Update which single period solver to use (if it depends on time)
         if hasattr(agent.solve_one_period, "__getitem__"):
-            solve_one_period = agent.solve_one_period[T - 1 - t]
+            solve_one_period = agent.solve_one_period[k]
         else:
             solve_one_period = agent.solve_one_period
 
@@ -1019,8 +1022,8 @@ def solve_one_cycle(agent, solution_last):
         # Update time-varying single period inputs
         for name in agent.time_vary:
             if name in these_args:
-                # solve_dict[name] = eval('agent.' + name + '[t]')
-                solve_dict[name] = agent.__dict__[name][T - 1 - t]
+                solve_dict[name] = agent.__dict__[name][k]
+
         solve_dict["solution_next"] = solution_next
 
         # Make a temporary dictionary for this period