Merge pull request #1013 from wdu9/WDU

Perfect Foresight MIT Shock Option Added

Documentation/CHANGELOG.md

@@ -25,7 +25,7 @@ PR: [#832](https://github.com/econ-ark/HARK/pull/832). See [this forthcoming REM
 * Using Lognormal.from_mean_std in the forward simulation of the RiskyAsset model [#1019](https://github.com/econ-ark/HARK/pull/1019)
 * Fix bug in DCEGM's primary kink finder due to numpy no longer accepting NaN in integer arrays [#990](https://github.com/econ-ark/HARK/pull/990).
 * Add a general class for consumers who can save using a risky asset [#1012](https://github.com/econ-ark/HARK/pull/1012/).
-
+Add Boolean attribute 'PerfMITShk' to consumption models. When true, allows perfect foresight MIT shocks to be simulated. [#1013](https://github.com/econ-ark/HARK/pull/1013).
 ### 0.11.0
 
 Release Data: March 4, 2021

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -1549,7 +1549,8 @@ 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': 1,        # Number of periods in the cycle for this agent type
+    "PerfMITShk": False    # Do Perfect Foresight MIT Shock: Forces Newborns to follow solution path of the agent he/she replaced when True
 }
 
 
@@ -1577,7 +1578,7 @@ class PerfForesightConsumerType(AgentType):
         MPCmax=1.0,
     )
     time_vary_ = ["LivPrb", "PermGroFac"]
-    time_inv_ = ["CRRA", "Rfree", "DiscFac", "MaxKinks", "BoroCnstArt"]
+    time_inv_ = ["CRRA", "Rfree", "DiscFac", "MaxKinks", "BoroCnstArt" ]
     state_vars = ['pLvl', 'PlvlAgg', 'bNrm', 'mNrm', "aNrm"]
     shock_vars_ = []
 
@@ -1707,9 +1708,11 @@ def sim_birth(self, which_agents):
             seed=self.RNG.randint(0, 2 ** 31 - 1)
         ).draw(N)
         self.t_age[which_agents] = 0  # How many periods since each agent was born
-        self.t_cycle[
-            which_agents
-        ] = 0  # Which period of the cycle each agent is currently in
+        if self.PerfMITShk == False:  # If True, Newborns inherit t_cycle of agent they replaced (i.e. t_cycles are not reset). 
+            self.t_cycle[
+                which_agents
+            ] = 0  # Which period of the cycle each agent is currently in
+
         return None
 
     def sim_death(self):
@@ -3051,4 +3054,4 @@ def construct_assets_grid(parameters):
 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]
-init_cyclical['T_cycle'] = 4
+init_cyclical['T_cycle'] = 4

HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py

@@ -6,7 +6,7 @@
 )
 import numpy as np
 import unittest
-from copy import copy
+from copy import copy, deepcopy
 
 
 class testIndShockConsumerType(unittest.TestCase):
@@ -399,3 +399,181 @@ def test_IndShock_stable_points(self):
         decimalPlacesTo = 10
         self.assertAlmostEqual(mNrmStE, 1.37731133865, decimalPlacesTo)
         self.assertAlmostEqual(mNrmTrg, 1.39101653806, decimalPlacesTo)
+
+
+JACDict={
+    # Parameters shared with the perfect foresight model
+    "CRRA":2,                             # Coefficient of relative risk aversion
+    "Rfree": 1.05**.25,                  # Interest factor on assets
+    "DiscFac": 0.972,                    # Intertemporal discount factor
+    "LivPrb" : [.99375],                    # Survival probability
+    "PermGroFac" :[1.00],                 # Permanent income growth factor
+
+    # Parameters that specify the income distribution over the lifecycle
+
+    "PermShkStd" :  [(0.01*4/11)**0.5],    # Standard deviation of log permanent shocks to income
+    "PermShkCount" : 5,                    # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd" : [.2],                   # Standard deviation of log transitory shocks to income
+    "TranShkCount" : 5,                    # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb" : 0.05,                     # Probability of unemployment while working
+    "IncUnemp" : 0.1,      # Unemployment benefits replacement rate
+    "UnempPrbRet" : 0.0005,                # Probability of "unemployment" while retired
+    "IncUnempRet" : 0.0,                   # "Unemployment" benefits when retired
+    "T_retire" : 0,                        # Period of retirement (0 --> no retirement)
+    "tax_rate" : 0.2,      # Flat income tax rate (legacy parameter, will be removed in future)
+
+    # Parameters for constructing the "assets above minimum" grid
+    "aXtraMin" : 0.001,                    # Minimum end-of-period "assets above minimum" value
+    "aXtraMax" : 15,                       # Maximum end-of-period "assets above minimum" value
+    "aXtraCount" : 48,                     # Number of points in the base grid of "assets above minimum"
+    "aXtraNestFac" : 3,                    # Exponential nesting factor when constructing "assets above minimum" grid
+    "aXtraExtra" : [None],                 # Additional values to add to aXtraGrid
+
+    # A few other parameters
+    "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" : 1,                         # Number of periods in the cycle for this agent type
+
+    # Parameters only used in simulation
+    "AgentCount" : 5000,                 # Number of agents of this type
+    "T_sim" : 100,                         # Number of periods to simulate
+    "aNrmInitMean" : np.log(2)-(.5**2)/2,# Mean of log initial assets
+    "aNrmInitStd"  : .5,                   # Standard deviation of log initial assets
+    "pLvlInitMean" : 0,                    # Mean of log initial permanent income
+    "pLvlInitStd"  : 0,                    # Standard deviation of log initial permanent income
+    "PermGroFacAgg" : 1.0,                 # Aggregate permanent income growth factor
+    "T_age" : None,                        # Age after which simulated agents are automatically killed
+
+}
+
+
+class testPerfMITShk(unittest.TestCase):
+
+    def jacobian(self):
+
+        class Test_agent(IndShockConsumerType):
+
+            def __init__(self, cycles= 0, **kwds):
+
+                IndShockConsumerType.__init__(self, cycles = 0, **kwds)
+
+            def get_Rfree(self):
+                """
+                Returns an array of size self.AgentCount with self.Rfree in every entry.
+                Parameters
+                ----------
+                None
+                Returns
+                -------
+                RfreeNow : np.array
+                     Array of size self.AgentCount with risk free interest rate for each agent.
+                """
+
+                if type(self.Rfree) == list:
+                    RfreeNow = self.Rfree[self.t_sim]* np.ones(self.AgentCount)
+                else:
+                    RfreeNow = ss.Rfree * np.ones(self.AgentCount)
+
+                return RfreeNow
+
+        ss = Test_agent(**JACDict )
+        ss.cycles = 0
+        ss.T_sim= 1200
+        ss.solve()
+        ss.initialize_sim()
+        ss.simulate()
+
+        class Test_agent2(Test_agent):
+
+             def transition(self):
+
+                pLvlPrev = self.state_prev['pLvl']
+                aNrmPrev = self.state_prev['aNrm']
+                RfreeNow = self.get_Rfree()
+
+                # Calculate new states: normalized market resources and permanent income level
+                pLvlNow = pLvlPrev*self.shocks['PermShk']  # Updated permanent income level
+                # Updated aggregate permanent productivity level
+                PlvlAggNow = self.state_prev['PlvlAgg']*self.PermShkAggNow
+                # "Effective" interest factor on normalized assets
+                ReffNow = RfreeNow/self.shocks['PermShk']
+                bNrmNow = ReffNow*aNrmPrev         # Bank balances before labor income
+                mNrmNow = bNrmNow + self.shocks['TranShk']  # Market resources after income
+
+
+                if self.t_sim == 0:
+
+                        mNrmNow = ss.state_now['mNrm']
+                        pLvlNow = ss.state_now['pLvl']
+
+                return pLvlNow, PlvlAggNow, bNrmNow, mNrmNow, None        
+
+
+
+        listA_g = []
+        params = deepcopy(JACDict)
+        params['T_cycle']= 200
+        params['LivPrb']= params['T_cycle']*[ss.LivPrb[0]]
+        params['PermGroFac']=params['T_cycle']*[1]
+        params['PermShkStd'] = params['T_cycle']*[(0.01*4/11)**0.5]
+        params['TranShkStd']= params['T_cycle']*[.2]
+        params['Rfree'] = params['T_cycle']*[ss.Rfree]
+
+        ss_dx = Test_agent2(**params )
+        ss_dx.pseudo_terminal = False
+        ss_dx.PerfMITShk = True
+        ss_dx.track_vars = ['aNrm','mNrm','cNrm','pLvl','aLvl']
+        ss_dx.cFunc_terminal_ = deepcopy(ss.solution[0].cFunc)
+        ss_dx.T_sim = params['T_cycle']
+        ss_dx.cycles= 1
+        ss_dx.IncShkDstn = params['T_cycle']*ss_dx.IncShkDstn
+        ss_dx.del_from_time_inv('Rfree')
+        ss_dx.add_to_time_vary('Rfree')
+
+        ss_dx.solve()
+        ss_dx.initialize_sim()
+        ss_dx.simulate()
+
+
+        for j in range(ss_dx.T_sim):
+
+            Ag = np.mean(ss_dx.history['aLvl'][j,:])
+            listA_g.append(Ag)
+
+        A_dx0 = np.array(listA_g)
+
+
+        ##############################################################################
+
+        example = Test_agent2(**params )
+        example.pseudo_terminal=False 
+        example.cFunc_terminal_ = deepcopy(ss.solution[0].cFunc)
+        example.T_sim = params['T_cycle']
+        example.cycles = 1
+        example.PerfMITShk = True
+        example.track_vars = ['aNrm','mNrm','cNrm','pLvl','aLvl']
+        example.del_from_time_inv('Rfree')
+        example.add_to_time_vary('Rfree')
+        example.IncShkDstn = params['T_cycle']*example.IncShkDstn
+
+        AHist =[]
+        listA = []
+        dx = .001
+        i = 50
+
+        example.Rfree = i *[ss.Rfree] + [ss.Rfree + dx] + (params['T_cycle']  - i - 1)*[ss.Rfree]
+
+        example.solve()
+        example.initialize_sim()
+        example.simulate()
+
+        for j in range(example.T_sim):
+
+            a = np.mean(example.history['aLvl'][j,:])
+            listA.append(a)
+
+        AHist.append(np.array(listA))
+        JACA = (AHist[0]-A_dx0)/(dx)
+
+        self.assertAlmostEqual(JACA[175], 6.441930322509393e-06)