Coverage for HARK/ConsumptionSaving/ConsMarkovModel.py: 95%
354 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +0000
1"""
2Classes to solve and simulate consumption-savings model with a discrete, exogenous,
3stochastic Markov state. The only solver here extends ConsIndShockModel to
4include a Markov state; the interest factor, permanent growth factor, and income
5distribution can vary with the discrete state.
6"""
8import numpy as np
10from HARK import AgentType, NullFunc
11from HARK.Calibration.Income.IncomeProcesses import (
12 construct_markov_lognormal_income_process_unemployment,
13 get_PermShkDstn_from_IncShkDstn_markov,
14 get_TranShkDstn_from_IncShkDstn_markov,
15)
16from HARK.ConsumptionSaving.ConsIndShockModel import (
17 ConsumerSolution,
18 IndShockConsumerType,
19 make_basic_CRRA_solution_terminal,
20 make_lognormal_kNrm_init_dstn,
21 make_lognormal_pLvl_init_dstn,
22)
23from HARK.distributions import MarkovProcess, Uniform, expected, DiscreteDistribution
24from HARK.interpolation import (
25 CubicInterp,
26 LinearInterp,
27 LowerEnvelope,
28 IndexedInterp,
29 MargMargValueFuncCRRA,
30 MargValueFuncCRRA,
31 ValueFuncCRRA,
32)
33from HARK.rewards import (
34 UtilityFuncCRRA,
35 CRRAutility,
36 CRRAutility_inv,
37 CRRAutility_invP,
38 CRRAutilityP,
39 CRRAutilityP_inv,
40 CRRAutilityP_invP,
41 CRRAutilityPP,
42)
43from HARK.utilities import make_assets_grid
45__all__ = ["MarkovConsumerType"]
47utility = CRRAutility
48utilityP = CRRAutilityP
49utilityPP = CRRAutilityPP
50utilityP_inv = CRRAutilityP_inv
51utility_invP = CRRAutility_invP
52utility_inv = CRRAutility_inv
53utilityP_invP = CRRAutilityP_invP
56###############################################################################
58# Define some functions that can be used as constructors for MrkvArray
61def make_simple_binary_markov(T_cycle, Mrkv_p11, Mrkv_p22):
62 """
63 Make a list of very simple Markov arrays between two binary states by specifying
64 diagonal elements in each period (probability of remaining in that state).
66 Parameters
67 ----------
68 T_cycle : int
69 Number of non-terminal periods in this instance's sequential problem.
70 Mrkv_p11 : [float]
71 List or array of probabilities of remaining in the first state between periods.
72 Mrkv_p22 : [float]
73 List or array of probabilities of remaining in the second state between periods.
75 Returns
76 -------
77 MrkvArray : [np.array]
78 List of 2x2 Markov transition arrays, one for each non-terminal period.
79 """
80 p11 = np.array(Mrkv_p11)
81 p22 = np.array(Mrkv_p22)
83 if len(p11) != T_cycle or len(p22) != T_cycle:
84 raise ValueError("Length of p11 and p22 probabilities must equal T_cycle!")
85 if np.any(p11 > 1.0) or np.any(p22 > 1.0):
86 raise ValueError("The p11 and p22 probabilities must not exceed 1!")
87 if np.any(p11 < 0.0) or np.any(p22 < 0.0):
88 raise ValueError("The p11 and p22 probabilities must not be less than 0!")
90 MrkvArray = [
91 np.array([[p11[t], 1.0 - p11[t]], [1.0 - p22[t], p22[t]]])
92 for t in range(T_cycle)
93 ]
94 return MrkvArray
97def make_ratchet_markov(T_cycle, Mrkv_ratchet_probs):
98 """
99 Make a list of "ratchet-style" Markov transition arrays, in which transitions
100 are strictly *one way* and only by one step. Each element of the ratchet_probs
101 list is a size-N vector giving the probability of progressing from state i to
102 state to state i+1 in that period; progress from the topmost state reverts the
103 agent to the 0th state. Set ratchet_probs[t][-1] to zero to make absorbing state.
105 Parameters
106 ----------
107 T_cycle : int
108 Number of non-terminal periods in this instance's sequential problem.
109 Mrkv_ratchet_probs : [np.array]
110 List of vectors of "ratchet probabilities" for each period.
112 Returns
113 -------
114 MrkvArray : [np.array]
115 List of NxN Markov transition arrays, one for each non-terminal period.
116 """
117 if len(Mrkv_ratchet_probs) != T_cycle:
118 raise ValueError("Length of Mrkv_ratchet_probs must equal T_cycle!")
120 N = Mrkv_ratchet_probs[0].size # number of discrete states
121 StateCount = np.array([Mrkv_ratchet_probs[t].size for t in range(T_cycle)])
122 if np.any(StateCount != N):
123 raise ValueError(
124 "All periods of the problem must have the same number of discrete states!"
125 )
127 MrkvArray = []
128 for t in range(T_cycle):
129 if np.any(Mrkv_ratchet_probs[t] > 1.0):
130 raise ValueError("Ratchet probabilities cannot exceed 1!")
131 if np.any(Mrkv_ratchet_probs[t] < 0.0):
132 raise ValueError("Ratchet probabilities cannot be below 0!")
134 MrkvArray_t = np.zeros((N, N))
135 for i in range(N):
136 p_go = Mrkv_ratchet_probs[t][i]
137 p_stay = 1.0 - p_go
138 if i < (N - 1):
139 i_next = i + 1
140 else:
141 i_next = 0
142 MrkvArray_t[i, i] = p_stay
143 MrkvArray_t[i, i_next] = p_go
145 MrkvArray.append(MrkvArray_t)
147 return MrkvArray
150def make_MrkvInitDstn(MrkvPrbsInit, RNG):
151 """
152 The constructor function for MrkvInitDstn, the distribution of Markov states
153 at model birth.
155 Parameters
156 ----------
157 MrkvPrbsInit : np.array
158 Stochastic vector specifying the distribution of initial discrete states.
159 RNG : np.random.RandomState
160 Agent's internal random number generator.
162 Returns
163 -------
164 MrkvInitDstn : DiscreteDistribution
165 Distribution from which discrete states at birth can be drawn.
166 """
167 seed = RNG.integers(0, 2**31 - 1)
168 vals = np.arange(MrkvPrbsInit.size, dtype=int)
169 MrkvInitDstn = DiscreteDistribution(pmv=MrkvPrbsInit, atoms=vals, seed=seed)
170 return MrkvInitDstn
173###############################################################################
176def make_markov_solution_terminal(CRRA, MrkvArray):
177 """
178 Make the terminal period solution for a consumption-saving model with a discrete
179 Markov state. Simply makes a basic terminal solution for IndShockConsumerType
180 and then replicates the attributes N times for the N states in the terminal period.
182 Parameters
183 ----------
184 CRRA : float
185 Coefficient of relative risk aversion.
186 MrkvArray : [np.array]
187 List of Markov transition probabilities arrays. Only used to find the
188 number of discrete states in the terminal period.
190 Returns
191 -------
192 solution_terminal : ConsumerSolution
193 Terminal period solution to the Markov consumption-saving problem.
194 """
195 solution_terminal_basic = make_basic_CRRA_solution_terminal(CRRA)
196 StateCount_T = MrkvArray[-1].shape[1]
197 N = StateCount_T # for shorter typing
199 # Make replicated terminal period solution: consume all resources, no human wealth, minimum m is 0
200 solution_terminal = ConsumerSolution(
201 cFunc=N * [solution_terminal_basic.cFunc],
202 vFunc=N * [solution_terminal_basic.vFunc],
203 vPfunc=N * [solution_terminal_basic.vPfunc],
204 vPPfunc=N * [solution_terminal_basic.vPPfunc],
205 mNrmMin=np.zeros(N),
206 hNrm=np.zeros(N),
207 MPCmin=np.ones(N),
208 MPCmax=np.ones(N),
209 )
210 solution_terminal.cFuncX = IndexedInterp(solution_terminal.cFunc)
211 return solution_terminal
214def solve_one_period_ConsMarkov(
215 solution_next,
216 IncShkDstn,
217 LivPrb,
218 DiscFac,
219 CRRA,
220 Rfree,
221 PermGroFac,
222 MrkvArray,
223 BoroCnstArt,
224 aXtraGrid,
225 vFuncBool,
226 CubicBool,
227):
228 """
229 Solves a single period consumption-saving problem with risky income and
230 stochastic transitions between discrete states, in a Markov fashion. Has
231 identical inputs as the ConsIndShock, except for a discrete Markov transition
232 rule MrkvArray. Markov states can differ in their interest factor, permanent
233 growth factor, and income distribution, so the inputs Rfree, PermGroFac, and
234 IncShkDstn are lists specifying those values in each (succeeding) Markov state.
236 Parameters
237 ----------
238 solution_next : ConsumerSolution
239 The solution to next period's one period problem.
240 IncShkDstn : [distribution.Distribution]
241 A length N list of income distributions in each succeeding Markov state.
242 Each income distribution is a discrete approximation to the income process
243 at the beginning of the succeeding period.
244 LivPrb : [float]
245 Survival probability; likelihood of being alive at the beginning of the
246 succeeding period conditional on the current state.
247 DiscFac : float
248 Intertemporal discount factor for future utility.
249 CRRA : float
250 Coefficient of relative risk aversion.
251 Rfree : [float]
252 Risk free interest factor on end-of-period assets for each Markov
253 state in the succeeding period.
254 PermGroGac : [float]
255 Expected permanent income growth factor at the end of this period
256 for each Markov state in the succeeding period.
257 MrkvArray : numpy.array
258 An NxN array representing a Markov transition matrix between discrete
259 states. The i,j-th element of MrkvArray is the probability of
260 moving from state i in period t to state j in period t+1.
261 BoroCnstArt: float or None
262 Borrowing constraint for the minimum allowable assets to end the
263 period with. If it is less than the natural borrowing constraint,
264 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
265 rowing constraint.
266 aXtraGrid: np.array
267 Array of "extra" end-of-period asset values-- assets above the
268 absolute minimum acceptable level.
269 vFuncBool: boolean
270 An indicator for whether the value function should be computed and
271 included in the reported solution.
272 CubicBool: boolean
273 An indicator for whether the solver should use cubic or linear inter-
274 polation.
276 Returns
277 -------
278 solution : ConsumerSolution
279 The solution to the single period consumption-saving problem. Includes
280 a consumption function cFunc (using cubic or linear splines), a marg-
281 inal value function vPfunc, a minimum acceptable level of normalized
282 market resources mNrmMin, normalized human wealth hNrm, and bounding
283 MPCs MPCmin and MPCmax. It might also have a value function vFunc
284 and marginal marginal value function vPPfunc. All of these attributes
285 are lists or arrays, with elements corresponding to the current
286 Markov state. E.g. solution.cFunc[0] is the consumption function
287 when in the i=0 Markov state this period.
288 """
289 # Relabel the inputs that vary across Markov states
290 IncShkDstn_list = IncShkDstn
291 Rfree_list = np.array(Rfree)
292 LivPrb_list = np.array(LivPrb)
293 PermGroFac_list = np.array(PermGroFac)
294 StateCountNow = MrkvArray.shape[0]
295 StateCountNext = MrkvArray.shape[1]
297 # Define the utility function
298 uFunc = UtilityFuncCRRA(CRRA)
300 # Initialize the natural borrowing constraint when entering each succeeding state
301 BoroCnstNat_temp = np.zeros(StateCountNext) + np.nan
303 # Find the natural borrowing constraint conditional on next period's state
304 for j in range(StateCountNext):
305 PermShkMinNext = np.min(IncShkDstn_list[j].atoms[0])
306 TranShkMinNext = np.min(IncShkDstn_list[j].atoms[1])
307 BoroCnstNat_temp[j] = (
308 (solution_next.mNrmMin[j] - TranShkMinNext)
309 * (PermGroFac_list[j] * PermShkMinNext)
310 / Rfree_list[j]
311 )
313 # Initialize the natural borrowing constraint and minimum value of mNrm for
314 # *this* period's Markov states, as well as a "dependency table"
315 BoroCnstNat_list = np.zeros(StateCountNow) + np.nan
316 mNrmMin_list = np.zeros(StateCountNow) + np.nan
317 BoroCnstDependency = np.zeros((StateCountNow, StateCountNext)) + np.nan
319 # The natural borrowing constraint in each current state is the *highest*
320 # among next-state-conditional natural borrowing constraints that could
321 # occur from this current state.
322 for i in range(StateCountNow):
323 possible_next_states = MrkvArray[i, :] > 0
324 BoroCnstNat_list[i] = np.max(BoroCnstNat_temp[possible_next_states])
326 # Explicitly handle the "None" case:
327 if BoroCnstArt is None:
328 mNrmMin_list[i] = BoroCnstNat_list[i]
329 else:
330 mNrmMin_list[i] = np.max([BoroCnstNat_list[i], BoroCnstArt])
331 BoroCnstDependency[i, :] = BoroCnstNat_list[i] == BoroCnstNat_temp
332 # Also creates a Boolean array indicating whether the natural borrowing
333 # constraint *could* be hit when transitioning from i to j.
335 # Initialize end-of-period (marginal) value functions, expected income conditional
336 # on the next state, and probability of getting the worst income shock in each
337 # succeeding period state
338 BegOfPrd_vFunc_list = []
339 BegOfPrd_vPfunc_list = []
340 Ex_IncNextAll = np.zeros(StateCountNext) + np.nan
341 WorstIncPrbAll = np.zeros(StateCountNext) + np.nan
343 # Loop through each next-period-state and calculate the beginning-of-period
344 # (marginal) value function
345 for j in range(StateCountNext):
346 # Condition values on next period's state (and record a couple for later use)
347 Rfree = Rfree_list[j]
348 PermGroFac = PermGroFac_list[j]
349 LivPrb = LivPrb_list[j]
350 # mNrmMinNow = self.mNrmMin_list[state_index]
351 BoroCnstNat = BoroCnstNat_temp[j]
353 # Unpack the income distribution in next period's Markov state
354 IncShkDstn = IncShkDstn_list[j]
355 ShkPrbsNext = IncShkDstn.pmv
356 PermShkValsNext = IncShkDstn.atoms[0]
357 TranShkValsNext = IncShkDstn.atoms[1]
358 PermShkMinNext = np.min(PermShkValsNext)
359 TranShkMinNext = np.min(TranShkValsNext)
361 # Calculate the probability that we get the worst possible income draw
362 IncNext = PermShkValsNext * TranShkValsNext
363 WorstIncNext = PermShkMinNext * TranShkMinNext
364 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
365 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
367 DiscFacEff = DiscFac # survival probability LivPrb represents probability
368 # from *current* state, so DiscFacEff is just DiscFac for now
370 # Unpack next period's (marginal) value function
371 vFuncNext = solution_next.vFunc[j] # This is None when vFuncBool is False
372 vPfuncNext = solution_next.vPfunc[j]
373 vPPfuncNext = solution_next.vPPfunc[j] # This is None when CubicBool is False
375 # Compute expected income next period and record worst income probability
376 Ex_IncNextAll[j] = np.dot(ShkPrbsNext, PermShkValsNext * TranShkValsNext)
377 WorstIncPrbAll[j] = WorstIncPrb
379 # Construct the BEGINNING-of-period marginal value function conditional
380 # on next period's state and add it to the list of value functions
382 # Get data to construct the end-of-period marginal value function (conditional on next state)
383 aNrmNext = np.asarray(aXtraGrid) + BoroCnstNat
385 # Define local functions for taking future expectations
386 def calc_mNrmNext(S, a, R):
387 return R / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
389 def calc_vNext(S, a, R):
390 return (
391 S["PermShk"] ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)
392 ) * vFuncNext(calc_mNrmNext(S, a, R))
394 def calc_vPnext(S, a, R):
395 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a, R))
397 def calc_vPPnext(S, a, R):
398 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, a, R))
400 # Calculate beginning-of-period marginal value of assets at each gridpoint
401 vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA)
402 BegOfPrd_vPnext = vPfacEff * expected(
403 calc_vPnext, IncShkDstn, args=(aNrmNext, Rfree)
404 )
406 # "Decurved" marginal value
407 BegOfPrd_vPnvrsNext = uFunc.derinv(BegOfPrd_vPnext, order=(1, 0))
409 # Make the beginning-of-period pseudo-inverse marginal value of assets
410 # function conditionalon next period's state
411 if CubicBool:
412 # Calculate end-of-period marginal marginal value of assets at each gridpoint
413 vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0)
414 BegOfPrd_vPPnext = vPPfacEff * expected(
415 calc_vPPnext, IncShkDstn, args=(aNrmNext, Rfree)
416 )
417 # "Decurved" marginal marginal value
418 BegOfPrd_vPnvrsPnext = BegOfPrd_vPPnext * uFunc.derinv(
419 BegOfPrd_vPnext, order=(1, 1)
420 )
422 # Construct the end-of-period marginal value function conditional on the next state.
423 BegOfPrd_vPnvrsFunc = CubicInterp(
424 aNrmNext,
425 BegOfPrd_vPnvrsNext,
426 BegOfPrd_vPnvrsPnext,
427 lower_extrap=True,
428 )
429 # TODO: Should not be lower extrap, add point at BoroCnstNat
430 else:
431 BegOfPrd_vPnvrsFunc = LinearInterp(
432 aNrmNext, BegOfPrd_vPnvrsNext, lower_extrap=True
433 )
434 # TODO: Should not be lower extrap, add point at BoroCnstNat
436 # "Recurve" the pseudo-inverse marginal value function
437 BegOfPrd_vPfunc = MargValueFuncCRRA(BegOfPrd_vPnvrsFunc, CRRA)
438 BegOfPrd_vPfunc_list.append(BegOfPrd_vPfunc)
440 # Construct the beginning-of-period value functional conditional on next
441 # period's state and add it to the list of value functions
442 if vFuncBool:
443 # Calculate end-of-period value, its derivative, and their pseudo-inverse
444 BegOfPrd_vNext = DiscFacEff * expected(
445 calc_vNext, IncShkDstn, args=(aNrmNext, Rfree)
446 )
447 # value transformed through inverse utility
448 BegOfPrd_vNvrsNext = uFunc.inv(BegOfPrd_vNext)
449 BegOfPrd_vNvrsPnext = BegOfPrd_vPnext * uFunc.derinv(
450 BegOfPrd_vNext, order=(0, 1)
451 )
452 BegOfPrd_vNvrsNext = np.insert(BegOfPrd_vNvrsNext, 0, 0.0)
453 BegOfPrd_vNvrsPnext = np.insert(
454 BegOfPrd_vNvrsPnext, 0, BegOfPrd_vNvrsPnext[0]
455 )
456 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
458 # Construct the end-of-period value function
459 aNrm_temp = np.insert(aNrmNext, 0, BoroCnstNat)
460 BegOfPrd_vNvrsFunc = LinearInterp(
461 aNrm_temp,
462 BegOfPrd_vNvrsNext,
463 )
464 BegOfPrd_vFunc = ValueFuncCRRA(BegOfPrd_vNvrsFunc, CRRA)
465 BegOfPrd_vFunc_list.append(BegOfPrd_vFunc)
467 # BegOfPrdvP is marginal value conditional on *next* period's state.
468 # Take expectations over Markov transitions to get EndOfPrdvP conditional on
469 # *this* period's Markov state.
471 # Find unique values of minimum acceptable end-of-period assets (and the
472 # current period states for which they apply).
473 aNrmMin_unique, Mrkv_inverse = np.unique(BoroCnstNat_list, return_inverse=True)
474 possible_transitions = MrkvArray > 0
476 # Initialize end-of-period marginal value (and marg marg value) at each
477 # asset gridpoint for each current period state
478 EndOfPrd_vP = np.zeros((StateCountNow, aXtraGrid.size))
479 EndOfPrd_vPP = np.zeros((StateCountNow, aXtraGrid.size))
481 # Calculate end-of-period marginal value (and marg marg value) at each
482 # asset gridpoint for each current period state, grouping current states
483 # by their natural borrowing constraint
484 for k in range(aNrmMin_unique.size):
485 # Get the states for which this minimum applies amd the aNrm grid for
486 # this set of current states
487 aNrmMin = aNrmMin_unique[k] # minimum assets for this pass
488 which_states = Mrkv_inverse == k
489 aNrmNow = aNrmMin + aXtraGrid # assets grid for this pass
491 # Make arrays to hold successor period's beginning-of-period (marginal)
492 # marginal value if we transition to it
493 BegOfPrd_vPnext = np.zeros((StateCountNext, aXtraGrid.size))
494 BegOfPrd_vPPnext = np.zeros((StateCountNext, aXtraGrid.size))
496 # Loop through future Markov states and fill in those values, but only
497 # look at future states that can actually be reached from our current
498 # set of states (for this natural borrowing constraint value)
499 for j in range(StateCountNext):
500 if not np.any(np.logical_and(possible_transitions[:, j], which_states)):
501 continue
503 BegOfPrd_vPnext[j, :] = BegOfPrd_vPfunc_list[j](aNrmNow)
504 # Add conditional end-of-period (marginal) marginal value to the arrays
505 if CubicBool:
506 BegOfPrd_vPPnext[j, :] = BegOfPrd_vPfunc_list[j].derivativeX(aNrmNow)
508 # Weight conditional marginal values by transition probabilities
509 # to get unconditional marginal (marginal) value at each gridpoint.
510 EndOfPrd_vP_temp = np.dot(MrkvArray, BegOfPrd_vPnext)
512 # Only take the states for which this asset minimum applies
513 EndOfPrd_vP[which_states, :] = EndOfPrd_vP_temp[which_states, :]
515 # Do the same thing for marginal marginal value
516 if CubicBool:
517 EndOfPrd_vPP_temp = np.dot(MrkvArray, BegOfPrd_vPPnext)
518 EndOfPrd_vPP[which_states, :] = EndOfPrd_vPP_temp[which_states, :]
520 # Store the results as attributes of self, scaling end of period marginal value by survival probability from each current state
521 LivPrb_tiled = np.tile(
522 np.reshape(LivPrb_list, (StateCountNow, 1)), (1, aXtraGrid.size)
523 )
524 EndOfPrd_vP = LivPrb_tiled * EndOfPrd_vP
525 if CubicBool:
526 EndOfPrd_vPP = LivPrb_tiled * EndOfPrd_vPP
528 # Calculate the bounding MPCs and PDV of human wealth for each state
530 # Calculate probability of getting the "worst" income shock and transition
531 # from each current state
532 WorstIncPrb_array = BoroCnstDependency * np.tile(
533 np.reshape(WorstIncPrbAll, (1, StateCountNext)), (StateCountNow, 1)
534 )
535 temp_array = MrkvArray * WorstIncPrb_array
536 WorstIncPrbNow = np.sum(temp_array, axis=1)
538 # Calculate expectation of upper bound of next period's MPC
539 Ex_MPCmaxNext = (
540 np.dot(temp_array, Rfree_list ** (1.0 - CRRA) * solution_next.MPCmax ** (-CRRA))
541 / WorstIncPrbNow
542 ) ** (-1.0 / CRRA)
544 # Calculate limiting upper bound of MPC this period for each Markov state
545 DiscFacEff_temp = DiscFac * LivPrb_list
546 MPCmaxNow = 1.0 / (
547 1.0 + ((DiscFacEff_temp * WorstIncPrbNow) ** (1.0 / CRRA)) / Ex_MPCmaxNext
548 )
549 MPCmaxEff = MPCmaxNow
550 MPCmaxEff[BoroCnstNat_list < mNrmMin_list] = 1.0
552 # Calculate the current Markov-state-conditional PDV of human wealth, correctly
553 # accounting for risky returns and risk aversion
554 hNrmPlusIncNext = Ex_IncNextAll + solution_next.hNrm
555 R_adj = np.dot(MrkvArray, Rfree_list ** (1.0 - CRRA))
556 hNrmNow = (
557 np.dot(MrkvArray, (PermGroFac_list / Rfree_list**CRRA) * hNrmPlusIncNext)
558 / R_adj
559 )
561 # Calculate the lower bound on MPC as m gets arbitrarily large
562 temp = (
563 DiscFacEff_temp
564 * np.dot(
565 MrkvArray, solution_next.MPCmin ** (-CRRA) * Rfree_list ** (1.0 - CRRA)
566 )
567 ) ** (1.0 / CRRA)
568 MPCminNow = 1.0 / (1.0 + temp)
570 # Find consumption and market resources corresponding to each end-of-period
571 # assets point for each state (and add an additional point at the lower bound)
572 aNrmNow = (aXtraGrid)[np.newaxis, :] + np.array(BoroCnstNat_list)[:, np.newaxis]
573 cNrmNow = uFunc.derinv(EndOfPrd_vP, order=(1, 0))
574 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
575 cNrmNow = np.hstack((np.zeros((StateCountNow, 1)), cNrmNow))
576 mNrmNow = np.hstack((np.reshape(mNrmMin_list, (StateCountNow, 1)), mNrmNow))
578 # Calculate the MPC at each market resource gridpoint in each state (if desired)
579 if CubicBool:
580 dcda = EndOfPrd_vPP / uFunc.der(cNrmNow[:, 1:], order=2) # drop first
581 MPCnow = dcda / (dcda + 1.0)
582 MPCnow = np.hstack((np.reshape(MPCmaxNow, (StateCountNow, 1)), MPCnow))
584 # Initialize an empty solution to which we'll add state-conditional solutions
585 solution = ConsumerSolution()
587 # Loop through each current period state and add its solution to the overall solution
588 for i in range(StateCountNow):
589 # Set current-Markov-state-conditional human wealth and MPC bounds
590 hNrmNow_i = hNrmNow[i]
591 MPCminNow_i = MPCminNow[i]
592 mNrmMin_i = mNrmMin_list[i]
594 # Construct the consumption function by combining the constrained and unconstrained portions
595 cFuncNowCnst = LinearInterp(
596 np.array([mNrmMin_list[i], mNrmMin_list[i] + 1.0]), np.array([0.0, 1.0])
597 )
598 if CubicBool:
599 cFuncNowUnc = CubicInterp(
600 mNrmNow[i, :],
601 cNrmNow[i, :],
602 MPCnow[i, :],
603 MPCminNow_i * hNrmNow_i,
604 MPCminNow_i,
605 )
606 else:
607 cFuncNowUnc = LinearInterp(
608 mNrmNow[i, :], cNrmNow[i, :], MPCminNow_i * hNrmNow_i, MPCminNow_i
609 )
610 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst)
612 # Make the marginal (marginal) value function
613 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
614 if CubicBool:
615 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
616 else:
617 vPPfuncNow = NullFunc()
619 # Make the value function for this state if requested
620 if vFuncBool:
621 # Make state-conditional grids of market resources and consumption
622 mNrm_for_vFunc = mNrmMin_i + aXtraGrid
623 cNrm_for_vFunc = cFuncNow(mNrm_for_vFunc)
624 aNrm_for_vFunc = mNrm_for_vFunc - cNrm_for_vFunc
626 # Calculate end-of-period value at each gridpoint
627 BegOfPrd_v_temp = np.zeros((StateCountNow, aXtraGrid.size))
628 for j in range(StateCountNext):
629 if possible_transitions[i, j]:
630 BegOfPrd_v_temp[j, :] = BegOfPrd_vFunc_list[j](aNrm_for_vFunc)
631 EndOfPrd_v = np.dot(MrkvArray[i, :], BegOfPrd_v_temp)
633 # Calculate (normalized) value and marginal value at each gridpoint
634 v_now = uFunc(cNrm_for_vFunc) + EndOfPrd_v
635 vP_now = uFunc.der(cNrm_for_vFunc)
637 # Make a "decurved" value function with the inverse utility function
638 # value transformed through inverse utility
639 vNvrs_now = uFunc.inv(v_now)
640 vNvrsP_now = vP_now * uFunc.derinv(v_now, order=(0, 1))
641 mNrm_temp = np.insert(mNrm_for_vFunc, 0, mNrmMin_i) # add the lower bound
642 vNvrs_now = np.insert(vNvrs_now, 0, 0.0)
643 vNvrsP_now = np.insert(
644 vNvrsP_now, 0, MPCmaxEff[i] ** (-CRRA / (1.0 - CRRA))
645 )
646 # MPCminNvrs = MPCminNow[i] ** (-CRRA / (1.0 - CRRA))
647 vNvrsFuncNow = LinearInterp(
648 mNrm_temp,
649 vNvrs_now,
650 # vNvrsP_now,
651 ) # MPCminNvrs * hNrmNow_i, MPCminNvrs)
652 # The bounding function for the pseudo-inverse value function is incorrect.
653 # TODO: Resolve this strange issue; extrapolation is suppressed for now.
655 # "Recurve" the decurved value function and add it to the list
656 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
658 else:
659 vFuncNow = NullFunc()
661 # Make the current-Markov-state-conditional solution
662 solution_cond = ConsumerSolution(
663 cFunc=cFuncNow,
664 vFunc=vFuncNow,
665 vPfunc=vPfuncNow,
666 vPPfunc=vPPfuncNow,
667 mNrmMin=mNrmMin_i,
668 )
670 # Add the current-state-conditional solution to the overall period solution
671 solution.append_solution(solution_cond)
673 # Add the lower bounds of market resources, MPC limits, human resources,
674 # and the value functions to the overall solution, then return it
675 solution.mNrmMin = mNrmMin_list
676 solution.hNrm = hNrmNow
677 solution.MPCmin = MPCminNow
678 solution.MPCmax = MPCmaxNow
679 solution.cFuncX = IndexedInterp(solution.cFunc)
680 return solution
683####################################################################################################
684####################################################################################################
686# Make a dictionary of constructors for the markov consumption-saving model
687markov_constructor_dict = {
688 "IncShkDstn": construct_markov_lognormal_income_process_unemployment,
689 "PermShkDstn": get_PermShkDstn_from_IncShkDstn_markov,
690 "TranShkDstn": get_TranShkDstn_from_IncShkDstn_markov,
691 "aXtraGrid": make_assets_grid,
692 "MrkvArray": make_simple_binary_markov,
693 "solution_terminal": make_markov_solution_terminal,
694 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
695 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
696 "MrkvInitDstn": make_MrkvInitDstn,
697}
699# Make a dictionary with parameters for the default constructor for kNrmInitDstn
700default_kNrmInitDstn_params = {
701 "kLogInitMean": -12.0, # Mean of log initial capital
702 "kLogInitStd": 0.0, # Stdev of log initial capital
703 "kNrmInitCount": 15, # Number of points in initial capital discretization
704}
706# Make a dictionary with parameters for the default constructor for pLvlInitDstn
707default_pLvlInitDstn_params = {
708 "pLogInitMean": 0.0, # Mean of log permanent income
709 "pLogInitStd": 0.0, # Stdev of log permanent income
710 "pLvlInitCount": 15, # Number of points in initial capital discretization
711}
713# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
714default_IncShkDstn_params = {
715 "PermShkStd": np.array(
716 [[0.1, 0.1]]
717 ), # Standard deviation of log permanent income shocks
718 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
719 "TranShkStd": np.array(
720 [[0.1, 0.1]]
721 ), # Standard deviation of log transitory income shocks
722 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
723 "UnempPrb": np.array([0.05, 0.05]), # Probability of unemployment while working
724 "IncUnemp": np.array(
725 [0.3, 0.3]
726 ), # Unemployment benefits replacement rate while working
727 "T_retire": 0, # Period of retirement (0 --> no retirement)
728 "UnempPrbRet": None, # Probability of "unemployment" while retired
729 "IncUnempRet": None, # "Unemployment" benefits when retired
730}
732# Default parameters to make aXtraGrid using make_assets_grid
733default_aXtraGrid_params = {
734 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
735 "aXtraMax": 20, # Maximum end-of-period "assets above minimum" value
736 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
737 "aXtraCount": 48, # Number of points in the grid of "assets above minimum"
738 "aXtraExtra": None, # Additional other values to add in grid (optional)
739}
741# Default parameters to make MrkvArray using make_simple_binary_markov
742default_MrkvArray_params = {
743 "Mrkv_p11": [0.9], # Probability of remaining in binary state 1
744 "Mrkv_p22": [0.4], # Probability of remaining in binary state 2
745}
747# Make a dictionary to specify an idiosyncratic income shocks consumer type
748init_indshk_markov = {
749 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
750 "cycles": 1, # Finite, non-cyclic model
751 "T_cycle": 1, # Number of periods in the cycle for this agent type
752 "constructors": markov_constructor_dict, # See dictionary above
753 "pseudo_terminal": False, # Terminal period really does exist
754 "global_markov": False, # Whether the Markov state is shared across agents
755 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
756 "CRRA": 2.0, # Coefficient of relative risk aversion
757 "Rfree": [np.array([1.03, 1.03])], # Interest factor on retained assets
758 "DiscFac": 0.96, # Intertemporal discount factor
759 "LivPrb": [np.array([0.98, 0.98])], # Survival probability after each period
760 "PermGroFac": [np.array([0.99, 1.03])], # Permanent income growth factor
761 "BoroCnstArt": 0.0, # Artificial borrowing constraint
762 "vFuncBool": False, # Whether to calculate the value function during solution
763 "CubicBool": False, # Whether to use cubic spline interpolation when True
764 # (Uses linear spline interpolation for cFunc when False)
765 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
766 "AgentCount": 10000, # Number of agents of this type
767 "T_age": None, # Age after which simulated agents are automatically killed
768 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
769 "MrkvPrbsInit": np.array([1.0, 0.0]), # Initial distribution of discrete state
770 # (The portion of PermGroFac attributable to aggregate productivity growth)
771 "NewbornTransShk": False, # Whether Newborns have transitory shock
772 # ADDITIONAL OPTIONAL PARAMETERS
773 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
774 # (Forces Newborns to follow solution path of the agent they replaced if True)
775 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
776}
777init_indshk_markov.update(default_IncShkDstn_params)
778init_indshk_markov.update(default_aXtraGrid_params)
779init_indshk_markov.update(default_MrkvArray_params)
780init_indshk_markov.update(default_kNrmInitDstn_params)
781init_indshk_markov.update(default_pLvlInitDstn_params)
784class MarkovConsumerType(IndShockConsumerType):
785 """
786 An agent in the Markov consumption-saving model. His problem is defined by a sequence
787 of income distributions, survival probabilities, discount factors, and permanent
788 income growth rates, as well as time invariant values for risk aversion, the
789 interest rate, the grid of end-of-period assets, and how he is borrowing constrained.
790 """
792 time_vary_ = IndShockConsumerType.time_vary_ + ["MrkvArray"]
794 # Mrkv is both a shock and a state
795 shock_vars_ = IndShockConsumerType.shock_vars_ + ["Mrkv"]
796 state_vars = IndShockConsumerType.state_vars + ["Mrkv"]
797 default_ = {
798 "params": init_indshk_markov,
799 "solver": solve_one_period_ConsMarkov,
800 "model": "ConsMarkov.yaml",
801 }
802 distributions = [
803 "IncShkDstn",
804 "PermShkDstn",
805 "TranShkDstn",
806 "kNrmInitDstn",
807 "pLvlInitDstn",
808 "MrkvInitDstn",
809 ]
811 def check_markov_inputs(self):
812 """
813 Many parameters used by MarkovConsumerType are arrays. Make sure those arrays are the
814 right shape.
816 Parameters
817 ----------
818 None
820 Returns
821 -------
822 None
823 """
824 StateCount = self.MrkvArray[0].shape[0]
826 # Check that arrays are the right shape
827 for t in range(self.T_cycle):
828 if not isinstance(self.Rfree[t], np.ndarray) or self.Rfree[t].shape != (
829 StateCount,
830 ):
831 raise ValueError(
832 "Rfree[t] not the right shape, it should be an array of Rfree of all the states."
833 )
835 # Check that arrays in lists are the right shape
836 for MrkvArray_t in self.MrkvArray:
837 if not isinstance(MrkvArray_t, np.ndarray) or MrkvArray_t.shape != (
838 StateCount,
839 StateCount,
840 ):
841 raise ValueError(
842 "MrkvArray not the right shape, it should be of the size states*states."
843 )
844 for LivPrb_t in self.LivPrb:
845 if not isinstance(LivPrb_t, np.ndarray) or LivPrb_t.shape != (StateCount,):
846 raise ValueError(
847 "Array in LivPrb is not the right shape, it should be an array of length equal to number of states"
848 )
849 for PermGroFac_t in self.PermGroFac:
850 if not isinstance(PermGroFac_t, np.ndarray) or PermGroFac_t.shape != (
851 StateCount,
852 ):
853 raise ValueError(
854 "Array in PermGroFac is not the right shape, it should be an array of length equal to number of states"
855 )
857 # Now check the income distribution.
858 # Note IncShkDstn is (potentially) time-varying, so it is in time_vary.
859 # Therefore it is a list, and each element of that list responds to the income distribution
860 # at a particular point in time.
861 for IncShkDstn_t in self.IncShkDstn:
862 if not isinstance(IncShkDstn_t, list):
863 raise ValueError(
864 "self.IncShkDstn is time varying and so must be a list"
865 + "of lists of Distributions, one per Markov State. Found "
866 + f"{self.IncShkDstn} instead"
867 )
868 elif len(IncShkDstn_t) != StateCount:
869 raise ValueError(
870 "List in IncShkDstn is not the right length, it should be length equal to number of states"
871 )
873 def pre_solve(self):
874 """
875 Check to make sure that the inputs that are specific to MarkovConsumerType
876 are of the right shape (if arrays) or length (if lists).
878 Parameters
879 ----------
880 None
882 Returns
883 -------
884 None
885 """
886 AgentType.pre_solve(self)
887 self.check_markov_inputs()
888 self.construct("solution_terminal")
890 def initialize_sim(self):
891 self.shocks["Mrkv"] = np.zeros(self.AgentCount, dtype=int)
892 IndShockConsumerType.initialize_sim(self)
893 if (
894 self.global_markov
895 ): # Need to initialize markov state to be the same for all agents
896 base_draw = Uniform(seed=self.RNG.integers(0, 2**31 - 1)).draw(1)
897 Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit))
898 self.shocks["Mrkv"] = np.ones(self.AgentCount) * np.searchsorted(
899 Cutoffs, base_draw
900 ).astype(int)
901 self.shocks["Mrkv"] = self.shocks["Mrkv"].astype(int)
903 def sim_death(self):
904 """
905 Determines which agents die this period and must be replaced. Uses the sequence in LivPrb
906 to determine survival probabilities for each agent.
908 Parameters
909 ----------
910 None
912 Returns
913 -------
914 which_agents : np.array(bool)
915 Boolean array of size AgentCount indicating which agents die.
916 """
917 # Determine who dies
918 LivPrb = np.array(self.LivPrb)[
919 self.t_cycle - 1, self.shocks["Mrkv"]
920 ] # Time has already advanced, so look back one
921 DiePrb = 1.0 - LivPrb
922 DeathShks = Uniform(seed=self.RNG.integers(0, 2**31 - 1)).draw(
923 N=self.AgentCount
924 )
925 which_agents = DeathShks < DiePrb
926 if self.T_age is not None: # Kill agents that have lived for too many periods
927 too_old = self.t_age >= self.T_age
928 which_agents = np.logical_or(which_agents, too_old)
929 return which_agents
931 def sim_birth(self, which_agents):
932 """
933 Makes new Markov consumer by drawing initial normalized assets, permanent income levels, and
934 discrete states. Calls IndShockConsumerType.sim_birth, then draws from initial Markov distribution.
936 Parameters
937 ----------
938 which_agents : np.array(Bool)
939 Boolean array of size self.AgentCount indicating which agents should be "born".
941 Returns
942 -------
943 None
944 """
945 # Get initial assets and permanent income
946 IndShockConsumerType.sim_birth(self, which_agents)
948 # Markov state is not changed if it is set at the global level
949 if not self.global_markov:
950 N = np.sum(which_agents)
951 self.state_now["Mrkv"][which_agents] = self.MrkvInitDstn.draw(N)
953 def get_markov_states(self):
954 """
955 Draw new Markov states for each agent in the simulated population, using
956 the attribute MrkvArray to determine transition probabilities.
958 Parameters
959 ----------
960 None
962 Returns
963 -------
964 None
965 """
966 dont_change = (
967 self.t_age == 0
968 ) # Don't change Markov state for those who were just born (unless global_markov)
969 if self.t_sim == 0: # Respect initial distribution of Markov states
970 dont_change[:] = True
972 # Determine which agents are in which states right now
973 MrkvPrev = self.shocks["Mrkv"]
974 MrkvNow = np.zeros(self.AgentCount, dtype=int)
976 # Draw new Markov states for each agent
977 for t in range(self.T_cycle):
978 markov_process = MarkovProcess(
979 self.MrkvArray[t], seed=self.RNG.integers(0, 2**31 - 1)
980 )
981 right_age = self.t_cycle == t
982 MrkvNow[right_age] = markov_process.draw(MrkvPrev[right_age])
983 if not self.global_markov:
984 MrkvNow[dont_change] = MrkvPrev[dont_change]
986 self.shocks["Mrkv"] = MrkvNow.astype(int)
988 def get_shocks(self):
989 """
990 Gets new Markov states and permanent and transitory income shocks for this period. Samples
991 from IncShkDstn for each period-state in the cycle.
993 Parameters
994 ----------
995 None
997 Returns
998 -------
999 None
1000 """
1001 self.get_markov_states()
1002 MrkvNow = self.shocks["Mrkv"]
1004 # Now get income shocks for each consumer, by cycle-time and discrete state
1005 PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
1006 TranShkNow = np.zeros(self.AgentCount)
1007 for t in range(self.T_cycle):
1008 for j in range(self.MrkvArray[t].shape[0]):
1009 these = np.logical_and(t == self.t_cycle, j == MrkvNow)
1010 N = np.sum(these)
1011 if N > 0:
1012 IncShkDstnNow = self.IncShkDstn[t - 1][
1013 j
1014 ] # set current income distribution
1015 PermGroFacNow = self.PermGroFac[t - 1][
1016 j
1017 ] # and permanent growth factor
1019 # Get random draws of income shocks from the discrete distribution
1020 EventDraws = IncShkDstnNow.draw_events(N)
1021 PermShkNow[these] = (
1022 IncShkDstnNow.atoms[0][EventDraws] * PermGroFacNow
1023 ) # permanent "shock" includes expected growth
1024 TranShkNow[these] = IncShkDstnNow.atoms[1][EventDraws]
1025 newborn = self.t_age == 0
1026 PermShkNow[newborn] = 1.0
1027 TranShkNow[newborn] = 1.0
1028 self.shocks["PermShk"] = PermShkNow
1029 self.shocks["TranShk"] = TranShkNow
1031 def read_shocks_from_history(self):
1032 """
1033 A slight modification of AgentType.read_shocks that makes sure that MrkvNow is int, not float.
1035 Parameters
1036 ----------
1037 None
1039 Returns
1040 -------
1041 None
1042 """
1043 IndShockConsumerType.read_shocks_from_history(self)
1044 self.shocks["Mrkv"] = self.shocks["Mrkv"].astype(int)
1046 def get_Rfree(self):
1047 """
1048 Returns an array of size self.AgentCount with interest factor that varies with discrete state.
1050 Parameters
1051 ----------
1052 None
1054 Returns
1055 -------
1056 RfreeNow : np.array
1057 Array of size self.AgentCount with risk free interest rate for each agent.
1058 """
1059 RfreeNow = np.zeros(self.AgentCount)
1060 for t in range(self.T_cycle):
1061 these = self.t_cycle == t
1062 RfreeNow[these] = self.Rfree[t][self.shocks["Mrkv"][these]]
1063 return RfreeNow
1065 def get_controls(self):
1066 """
1067 Calculates consumption for each consumer of this type using the consumption functions.
1069 Parameters
1070 ----------
1071 None
1073 Returns
1074 -------
1075 None
1076 """
1077 cNrmNow = np.zeros(self.AgentCount) + np.nan
1078 MPCnow = np.zeros(self.AgentCount) + np.nan
1079 J = self.MrkvArray[0].shape[0]
1081 MrkvBoolArray = np.zeros((J, self.AgentCount), dtype=bool)
1082 for j in range(J):
1083 MrkvBoolArray[j, :] = j == self.shocks["Mrkv"]
1085 for t in range(self.T_cycle):
1086 right_t = t == self.t_cycle
1087 for j in range(J):
1088 these = np.logical_and(right_t, MrkvBoolArray[j, :])
1089 cNrmNow[these], MPCnow[these] = (
1090 self.solution[t]
1091 .cFunc[j]
1092 .eval_with_derivative(self.state_now["mNrm"][these])
1093 )
1094 self.controls["cNrm"] = cNrmNow
1095 self.MPCnow = MPCnow
1097 def get_poststates(self):
1098 super().get_poststates()
1099 self.state_now["Mrkv"] = self.shocks["Mrkv"].copy()
1101 def calc_bounding_values(self):
1102 """
1103 Calculate human wealth plus minimum and maximum MPC in an infinite
1104 horizon model with only one period repeated indefinitely. Store results
1105 as attributes of self. Human wealth is the present discounted value of
1106 expected future income after receiving income this period, ignoring mort-
1107 ality. The maximum MPC is the limit of the MPC as m --> mNrmMin. The
1108 minimum MPC is the limit of the MPC as m --> infty. Results are all
1109 np.array with elements corresponding to each Markov state.
1111 NOT YET IMPLEMENTED FOR THIS CLASS
1113 Parameters
1114 ----------
1115 None
1117 Returns
1118 -------
1119 None
1120 """
1121 raise NotImplementedError()
1123 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
1124 """
1125 Creates a "normalized Euler error" function for this instance, mapping
1126 from market resources to "consumption error per dollar of consumption."
1127 Stores result in attribute eulerErrorFunc as an interpolated function.
1128 Has option to use approximate income distribution stored in self.IncShkDstn
1129 or to use a (temporary) very dense approximation.
1131 NOT YET IMPLEMENTED FOR THIS CLASS
1133 Parameters
1134 ----------
1135 mMax : float
1136 Maximum normalized market resources for the Euler error function.
1137 approx_inc_dstn : Boolean
1138 Indicator for whether to use the approximate discrete income distri-
1139 bution stored in self.IncShkDstn[0], or to use a very accurate
1140 discrete approximation instead. When True, uses approximation in
1141 IncShkDstn; when False, makes and uses a very dense approximation.
1143 Returns
1144 -------
1145 None
1147 Notes
1148 -----
1149 This method is not used by any other code in the library. Rather, it is here
1150 for expository and benchmarking purposes.
1151 """
1152 raise NotImplementedError()