Coverage for HARK / ConsumptionSaving / ConsIndShockModel.py: 99%
882 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-07 05:16 +0000
1"""
2Classes to solve canonical consumption-saving models with idiosyncratic shocks
3to income. All models here assume CRRA utility with geometric discounting, no
4bequest motive, and income shocks that are fully transitory or fully permanent.
6It currently solves three types of models:
7 1) A very basic "perfect foresight" consumption-savings model with no uncertainty.
8 2) A consumption-savings model with risk over transitory and permanent income shocks.
9 3) The model described in (2), with an interest rate for debt that differs
10 from the interest rate for savings.
12See NARK https://github.com/econ-ark/HARK/blob/master/docs/NARK/NARK.pdf for information on variable naming conventions.
13See HARK documentation for mathematical descriptions of the models being solved.
14"""
16from copy import copy
18import numpy as np
19from HARK.Calibration.Income.IncomeTools import (
20 Cagetti_income,
21 parse_income_spec,
22 parse_time_params,
23)
24from HARK.Calibration.Income.IncomeProcesses import (
25 construct_lognormal_income_process_unemployment,
26 get_PermShkDstn_from_IncShkDstn,
27 get_TranShkDstn_from_IncShkDstn,
28)
29from HARK.Calibration.life_tables.us_ssa.SSATools import parse_ssa_life_table
30from HARK.Calibration.SCF.WealthIncomeDist.SCFDistTools import (
31 income_wealth_dists_from_scf,
32)
33from HARK.distributions import (
34 Lognormal,
35 MeanOneLogNormal,
36 Uniform,
37 add_discrete_outcome_constant_mean,
38 combine_indep_dstns,
39 expected,
40)
41from HARK.interpolation import (
42 LinearInterp,
43 LowerEnvelope,
44 MargMargValueFuncCRRA,
45 MargValueFuncCRRA,
46 ValueFuncCRRA,
47)
48from HARK.interpolation import CubicHermiteInterp as CubicInterp
49from HARK.metric import MetricObject
50from HARK.rewards import (
51 CRRAutility,
52 CRRAutility_inv,
53 CRRAutility_invP,
54 CRRAutilityP,
55 CRRAutilityP_inv,
56 CRRAutilityP_invP,
57 CRRAutilityPP,
58 UtilityFuncCRRA,
59)
60from HARK.utilities import make_assets_grid
61from scipy.optimize import newton
63from HARK import (
64 AgentType,
65 NullFunc,
66 _log,
67 set_verbosity_level,
68)
70__all__ = [
71 "ConsumerSolution",
72 "PerfForesightConsumerType",
73 "IndShockConsumerType",
74 "KinkedRconsumerType",
75 "init_perfect_foresight",
76 "init_idiosyncratic_shocks",
77 "init_kinked_R",
78 "init_lifecycle",
79 "init_cyclical",
80]
82utility = CRRAutility
83utilityP = CRRAutilityP
84utilityPP = CRRAutilityPP
85utilityP_inv = CRRAutilityP_inv
86utility_invP = CRRAutility_invP
87utility_inv = CRRAutility_inv
88utilityP_invP = CRRAutilityP_invP
91# =====================================================================
92# === Classes that help solve consumption-saving models ===
93# =====================================================================
96class ConsumerSolution(MetricObject):
97 r"""
98 A class representing the solution of a single period of a consumption-saving
99 problem. The solution must include a consumption function and marginal
100 value function.
102 Here and elsewhere in the code, Nrm indicates that variables are normalized
103 by permanent income.
105 Parameters
106 ----------
107 cFunc : function
108 The consumption function for this period, defined over normalized market
109 resources: cNrm = cFunc(mNrm).
110 vFunc : function
111 The beginning-of-period value function for this period, defined over
112 normalized market resources: vNrm = vFunc(mNrm).
113 vPfunc : function
114 The beginning-of-period marginal value function for this period,
115 defined over normalized market resources: vNrmP = vPfunc(mNrm).
116 vPPfunc : function
117 The beginning-of-period marginal marginal value function for this
118 period, defined over normalized market resources: vNrmPP = vPPfunc(mNrm).
119 mNrmMin : float
120 The minimum allowable normalized market resources for this period; the consump-
121 tion function (etc) are undefined for m < mNrmMin.
122 hNrm : float
123 Normalized human wealth after receiving income this period: PDV of all future
124 income, ignoring mortality.
125 MPCmin : float
126 Infimum of the marginal propensity to consume this period.
127 MPC --> MPCmin as m --> infinity.
128 MPCmax : float
129 Supremum of the marginal propensity to consume this period.
130 MPC --> MPCmax as m --> mNrmMin.
132 """
134 distance_criteria = ["vPfunc"]
136 def __init__(
137 self,
138 cFunc=None,
139 vFunc=None,
140 vPfunc=None,
141 vPPfunc=None,
142 mNrmMin=None,
143 hNrm=None,
144 MPCmin=None,
145 MPCmax=None,
146 ):
147 # Change any missing function inputs to NullFunc
148 self.cFunc = cFunc if cFunc is not None else NullFunc()
149 self.vFunc = vFunc if vFunc is not None else NullFunc()
150 self.vPfunc = vPfunc if vPfunc is not None else NullFunc()
151 # vPFunc = NullFunc() if vPfunc is None else vPfunc
152 self.vPPfunc = vPPfunc if vPPfunc is not None else NullFunc()
153 self.mNrmMin = mNrmMin
154 self.hNrm = hNrm
155 self.MPCmin = MPCmin
156 self.MPCmax = MPCmax
158 def append_solution(self, new_solution):
159 """
160 Appends one solution to another to create a ConsumerSolution whose
161 attributes are lists. Used in ConsMarkovModel, where we append solutions
162 *conditional* on a particular value of a Markov state to each other in
163 order to get the entire solution.
165 Parameters
166 ----------
167 new_solution : ConsumerSolution
168 The solution to a consumption-saving problem; each attribute is a
169 list representing state-conditional values or functions.
171 Returns
172 -------
173 None
174 """
175 if type(self.cFunc) != list:
176 # Then we assume that self is an empty initialized solution instance.
177 # Begin by checking this is so.
178 assert NullFunc().distance(self.cFunc) == 0, (
179 "append_solution called incorrectly!"
180 )
182 # We will need the attributes of the solution instance to be lists. Do that here.
183 self.cFunc = [new_solution.cFunc]
184 self.vFunc = [new_solution.vFunc]
185 self.vPfunc = [new_solution.vPfunc]
186 self.vPPfunc = [new_solution.vPPfunc]
187 self.mNrmMin = [new_solution.mNrmMin]
188 else:
189 self.cFunc.append(new_solution.cFunc)
190 self.vFunc.append(new_solution.vFunc)
191 self.vPfunc.append(new_solution.vPfunc)
192 self.vPPfunc.append(new_solution.vPPfunc)
193 self.mNrmMin.append(new_solution.mNrmMin)
196# =====================================================================
197# == Functions for initializing newborns in consumption-saving models =
198# =====================================================================
201def make_lognormal_kNrm_init_dstn(kLogInitMean, kLogInitStd, kNrmInitCount, RNG):
202 """
203 Construct a lognormal distribution for (normalized) initial capital holdings
204 of newborns, kNrm. This is the default constructor for kNrmInitDstn.
206 Parameters
207 ----------
208 kLogInitMean : float
209 Mean of log capital holdings for newborns.
210 kLogInitStd : float
211 Stdev of log capital holdings for newborns.
212 kNrmInitCount : int
213 Number of points in the discretization.
214 RNG : np.random.RandomState
215 Agent's internal RNG.
217 Returns
218 -------
219 kNrmInitDstn : DiscreteDistribution
220 Discretized distribution of initial capital holdings for newborns.
221 """
222 dstn = Lognormal(
223 mu=kLogInitMean,
224 sigma=kLogInitStd,
225 seed=RNG.integers(0, 2**31 - 1),
226 )
227 kNrmInitDstn = dstn.discretize(kNrmInitCount)
228 return kNrmInitDstn
231def make_lognormal_pLvl_init_dstn(pLogInitMean, pLogInitStd, pLvlInitCount, RNG):
232 """
233 Construct a lognormal distribution for initial permanent income level of
234 newborns, pLvl. This is the default constructor for pLvlInitDstn.
236 Parameters
237 ----------
238 pLogInitMean : float
239 Mean of log permanent income for newborns.
240 pLogInitStd : float
241 Stdev of log capital holdings for newborns.
242 pLvlInitCount : int
243 Number of points in the discretization.
244 RNG : np.random.RandomState
245 Agent's internal RNG.
247 Returns
248 -------
249 pLvlInitDstn : DiscreteDistribution
250 Discretized distribution of initial permanent income for newborns.
251 """
252 dstn = Lognormal(
253 mu=pLogInitMean,
254 sigma=pLogInitStd,
255 seed=RNG.integers(0, 2**31 - 1),
256 )
257 pLvlInitDstn = dstn.discretize(pLvlInitCount)
258 return pLvlInitDstn
261# =====================================================================
262# === Classes and functions that solve consumption-saving models ===
263# =====================================================================
266def calc_human_wealth(h_nrm_next, perm_gro_fac, rfree, ex_inc_next):
267 """Calculate human wealth this period given human wealth next period.
269 Args:
270 h_nrm_next (float): Normalized human wealth next period.
271 perm_gro_fac (float): Permanent income growth factor.
272 rfree (float): Risk free interest factor.
273 ex_inc_next (float): Expected income next period.
274 """
275 return (perm_gro_fac / rfree) * (h_nrm_next + ex_inc_next)
278def calc_patience_factor(rfree, disc_fac_eff, crra):
279 """Calculate the patience factor for the agent.
281 Args:
282 rfree (float): Risk free interest factor.
283 disc_fac_eff (float): Effective discount factor.
284 crra (float): Coefficient of relative risk aversion.
286 """
287 return ((rfree * disc_fac_eff) ** (1.0 / crra)) / rfree
290def calc_mpc_min(mpc_min_next, pat_fac):
291 """Calculate the lower bound of the marginal propensity to consume.
293 Args:
294 mpc_min_next (float): Lower bound of the marginal propensity to
295 consume next period.
296 pat_fac (float): Patience factor.
297 """
298 return 1.0 / (1.0 + pat_fac / mpc_min_next)
301def solve_one_period_ConsPF(
302 solution_next,
303 DiscFac,
304 LivPrb,
305 CRRA,
306 Rfree,
307 PermGroFac,
308 BoroCnstArt,
309 MaxKinks,
310):
311 """Solves one period of a basic perfect foresight consumption-saving model with
312 a single risk free asset and permanent income growth.
314 Parameters
315 ----------
316 solution_next : ConsumerSolution
317 The solution to next period's one-period problem.
318 DiscFac : float
319 Intertemporal discount factor for future utility.
320 LivPrb : float
321 Survival probability; likelihood of being alive at the beginning of
322 the next period.
323 CRRA : float
324 Coefficient of relative risk aversion.
325 Rfree : float
326 Risk free interest factor on end-of-period assets.
327 PermGroFac : float
328 Expected permanent income growth factor at the end of this period.
329 BoroCnstArt : float or None
330 Artificial borrowing constraint, as a multiple of permanent income.
331 Can be None, indicating no artificial constraint.
332 MaxKinks : int
333 Maximum number of kink points to allow in the consumption function;
334 additional points will be thrown out. Only relevant in infinite
335 horizon model with artificial borrowing constraint.
337 Returns
338 -------
339 solution_now : ConsumerSolution
340 Solution to the current period of a perfect foresight consumption-saving
341 problem.
343 """
344 # Define the utility function and effective discount factor
345 uFunc = UtilityFuncCRRA(CRRA)
346 DiscFacEff = DiscFac * LivPrb # Effective = pure x LivPrb
348 # Prevent comparing None and float if there is no borrowing constraint
349 # Can borrow as much as we want
350 BoroCnstArt = -np.inf if BoroCnstArt is None else BoroCnstArt
352 # Calculate human wealth this period
353 hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, 1.0)
355 # Calculate the lower bound of the marginal propensity to consume
356 PatFac = calc_patience_factor(Rfree, DiscFacEff, CRRA)
357 MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac)
359 # Extract the discrete kink points in next period's consumption function;
360 # don't take the last one, as it only defines the extrapolation and is not a kink.
361 mNrmNext = solution_next.cFunc.x_list[:-1]
362 cNrmNext = solution_next.cFunc.y_list[:-1]
363 vFuncNvrsNext = solution_next.vFunc.vFuncNvrs.y_list[:-1]
364 EndOfPrdv = DiscFacEff * PermGroFac ** (1.0 - CRRA) * uFunc(vFuncNvrsNext)
366 # Calculate the end-of-period asset values that would reach those kink points
367 # next period, then invert the first order condition to get consumption. Then
368 # find the endogenous gridpoint (kink point) today that corresponds to each kink
369 aNrmNow = (PermGroFac / Rfree) * (mNrmNext - 1.0)
370 cNrmNow = (DiscFacEff * Rfree) ** (-1.0 / CRRA) * (PermGroFac * cNrmNext)
371 mNrmNow = aNrmNow + cNrmNow
373 # Calculate (pseudo-inverse) value at each consumption kink point
374 vNow = uFunc(cNrmNow) + EndOfPrdv
375 vNvrsNow = uFunc.inverse(vNow)
376 vNvrsSlopeMin = MPCminNow ** (-CRRA / (1.0 - CRRA))
378 # Add an additional point to the list of gridpoints for the extrapolation,
379 # using the new value of the lower bound of the MPC.
380 mNrmNow = np.append(mNrmNow, mNrmNow[-1] + 1.0)
381 cNrmNow = np.append(cNrmNow, cNrmNow[-1] + MPCminNow)
382 vNvrsNow = np.append(vNvrsNow, vNvrsNow[-1] + vNvrsSlopeMin)
384 # If the artificial borrowing constraint binds, combine the constrained and
385 # unconstrained consumption functions.
386 if BoroCnstArt > mNrmNow[0]:
387 # Find the highest index where constraint binds
388 cNrmCnst = mNrmNow - BoroCnstArt
389 CnstBinds = cNrmCnst < cNrmNow
390 idx = np.where(CnstBinds)[0][-1]
392 if idx < (mNrmNow.size - 1):
393 # If it is not the *very last* index, find the the critical level
394 # of mNrm where the artificial borrowing contraint begins to bind.
395 d0 = cNrmNow[idx] - cNrmCnst[idx]
396 d1 = cNrmCnst[idx + 1] - cNrmNow[idx + 1]
397 m0 = mNrmNow[idx]
398 m1 = mNrmNow[idx + 1]
399 alpha = d0 / (d0 + d1)
400 mCrit = m0 + alpha * (m1 - m0)
402 # Adjust the grids of mNrm and cNrm to account for the borrowing constraint.
403 cCrit = mCrit - BoroCnstArt
404 mNrmNow = np.concatenate(([BoroCnstArt, mCrit], mNrmNow[(idx + 1) :]))
405 cNrmNow = np.concatenate(([0.0, cCrit], cNrmNow[(idx + 1) :]))
407 # Adjust the vNvrs grid to account for the borrowing constraint
408 v0 = vNvrsNow[idx]
409 v1 = vNvrsNow[idx + 1]
410 vNvrsCrit = v0 + alpha * (v1 - v0)
411 vNvrsNow = np.concatenate(([0.0, vNvrsCrit], vNvrsNow[(idx + 1) :]))
413 else:
414 # If it *is* the very last index, then there are only three points
415 # that characterize the consumption function: the artificial borrowing
416 # constraint, the constraint kink, and the extrapolation point.
417 mXtra = (cNrmNow[-1] - cNrmCnst[-1]) / (1.0 - MPCminNow)
418 mCrit = mNrmNow[-1] + mXtra
419 cCrit = mCrit - BoroCnstArt
420 mNrmNow = np.array([BoroCnstArt, mCrit, mCrit + 1.0])
421 cNrmNow = np.array([0.0, cCrit, cCrit + MPCminNow])
423 # Adjust vNvrs grid for this three node structure
424 mNextCrit = BoroCnstArt * Rfree + 1.0
425 vNextCrit = PermGroFac ** (1.0 - CRRA) * solution_next.vFunc(mNextCrit)
426 vCrit = uFunc(cCrit) + DiscFacEff * vNextCrit
427 vNvrsCrit = uFunc.inverse(vCrit)
428 vNvrsNow = np.array([0.0, vNvrsCrit, vNvrsCrit + vNvrsSlopeMin])
430 # If the mNrm and cNrm grids have become too large, throw out the last
431 # kink point, being sure to adjust the extrapolation.
432 if mNrmNow.size > MaxKinks:
433 mNrmNow = np.concatenate((mNrmNow[:-2], [mNrmNow[-3] + 1.0]))
434 cNrmNow = np.concatenate((cNrmNow[:-2], [cNrmNow[-3] + MPCminNow]))
435 vNvrsNow = np.concatenate((vNvrsNow[:-2], [vNvrsNow[-3] + vNvrsSlopeMin]))
437 # Construct the consumption function as a linear interpolation.
438 cFuncNow = LinearInterp(mNrmNow, cNrmNow)
440 # Calculate the upper bound of the MPC as the slope of the bottom segment.
441 MPCmaxNow = (cNrmNow[1] - cNrmNow[0]) / (mNrmNow[1] - mNrmNow[0])
442 mNrmMinNow = mNrmNow[0]
444 # Construct the (marginal) value function for this period
445 # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations
446 vFuncNvrs = LinearInterp(mNrmNow, vNvrsNow)
447 vFuncNow = ValueFuncCRRA(vFuncNvrs, CRRA)
448 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
450 # Construct and return the solution
451 solution_now = ConsumerSolution(
452 cFunc=cFuncNow,
453 vFunc=vFuncNow,
454 vPfunc=vPfuncNow,
455 mNrmMin=mNrmMinNow,
456 hNrm=hNrmNow,
457 MPCmin=MPCminNow,
458 MPCmax=MPCmaxNow,
459 )
460 return solution_now
463def calc_worst_inc_prob(inc_shk_dstn, use_infimum=False):
464 """Calculate the probability of the worst income shock.
466 Args:
467 inc_shk_dstn (DiscreteDistribution): Distribution of shocks to income.
468 use_infimum (bool): Indicator for whether to try to use the infimum of the limiting (true) income distribution.
469 """
470 probs = inc_shk_dstn.pmv
471 perm, tran = inc_shk_dstn.atoms
472 income = perm * tran
473 if use_infimum:
474 worst_inc = np.prod(inc_shk_dstn.limit["infimum"])
475 else:
476 worst_inc = np.min(income)
477 return np.sum(probs[income == worst_inc])
480def calc_boro_const_nat(
481 m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac, use_infimum=False
482):
483 """Calculate the natural borrowing constraint.
485 Args:
486 m_nrm_min_next (float): Minimum normalized market resources next period.
487 inc_shk_dstn (DiscreteDstn): Distribution of shocks to income.
488 rfree (float): Risk free interest factor.
489 perm_gro_fac (float): Permanent income growth factor.
490 use_infimum (bool): Indicator for whether to use the infimum of the limiting (true) income distribution
491 """
492 if use_infimum:
493 perm_min, tran_min = inc_shk_dstn.limit["infimum"]
494 else:
495 perm, tran = inc_shk_dstn.atoms
496 perm_min = np.min(perm)
497 tran_min = np.min(tran)
499 temp_fac = (perm_gro_fac * perm_min) / rfree
500 boro_cnst_nat = (m_nrm_min_next - tran_min) * temp_fac
501 return boro_cnst_nat
504def calc_m_nrm_min(boro_const_art, boro_const_nat):
505 """Calculate the minimum normalized market resources this period.
507 Args:
508 boro_const_art (float): Artificial borrowing constraint.
509 boro_const_nat (float): Natural borrowing constraint.
510 """
511 return (
512 boro_const_nat
513 if boro_const_art is None
514 else max(boro_const_nat, boro_const_art)
515 )
518def calc_mpc_max(
519 mpc_max_next, worst_inc_prob, crra, pat_fac, boro_const_nat, boro_const_art
520):
521 """Calculate the upper bound of the marginal propensity to consume.
523 Args:
524 mpc_max_next (float): Upper bound of the marginal propensity to
525 consume next period.
526 worst_inc_prob (float): Probability of the worst income shock.
527 crra (float): Coefficient of relative risk aversion.
528 pat_fac (float): Patience factor.
529 boro_const_nat (float): Natural borrowing constraint.
530 boro_const_art (float): Artificial borrowing constraint.
531 """
532 temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac
533 return 1.0 / (1.0 + temp_fac / mpc_max_next)
536def calc_m_nrm_next(shock, a, rfree, perm_gro_fac):
537 """Calculate normalized market resources next period.
539 Args:
540 shock (float): Realization of shocks to income.
541 a (np.ndarray): Exogenous grid of end-of-period assets.
542 rfree (float): Risk free interest factor.
543 perm_gro_fac (float): Permanent income growth factor.
544 """
545 return rfree / (perm_gro_fac * shock["PermShk"]) * a + shock["TranShk"]
548def calc_v_next(shock, a, rfree, crra, perm_gro_fac, vfunc_next):
549 """Calculate continuation value function with respect to
550 end-of-period assets.
552 Args:
553 shock (float): Realization of shocks to income.
554 a (np.ndarray): Exogenous grid of end-of-period assets.
555 rfree (float): Risk free interest factor.
556 crra (float): Coefficient of relative risk aversion.
557 perm_gro_fac (float): Permanent income growth factor.
558 vfunc_next (Callable): Value function next period.
559 """
560 return (
561 shock["PermShk"] ** (1.0 - crra) * perm_gro_fac ** (1.0 - crra)
562 ) * vfunc_next(calc_m_nrm_next(shock, a, rfree, perm_gro_fac))
565def calc_vp_next(shock, a, rfree, crra, perm_gro_fac, vp_func_next):
566 """Calculate the continuation marginal value function with respect to
567 end-of-period assets.
569 Args:
570 shock (float): Realization of shocks to income.
571 a (np.ndarray): Exogenous grid of end-of-period assets.
572 rfree (float): Risk free interest factor.
573 crra (float): Coefficient of relative risk aversion.
574 perm_gro_fac (float): Permanent income growth factor.
575 vp_func_next (Callable): Marginal value function next period.
576 """
577 return shock["PermShk"] ** (-crra) * vp_func_next(
578 calc_m_nrm_next(shock, a, rfree, perm_gro_fac),
579 )
582def calc_vpp_next(shock, a, rfree, crra, perm_gro_fac, vppfunc_next):
583 """Calculate the continuation marginal marginal value function
584 with respect to end-of-period assets.
586 Args:
587 shock (float): Realization of shocks to income.
588 a (np.ndarray): Exogenous grid of end-of-period assets.
589 rfree (float): Risk free interest factor.
590 crra (float): Coefficient of relative risk aversion.
591 perm_gro_fac (float): Permanent income growth factor.
592 vppfunc_next (Callable): Marginal marginal value function next period.
593 """
594 return shock["PermShk"] ** (-crra - 1.0) * vppfunc_next(
595 calc_m_nrm_next(shock, a, rfree, perm_gro_fac),
596 )
599def solve_one_period_ConsIndShock(
600 solution_next,
601 IncShkDstn,
602 LivPrb,
603 DiscFac,
604 CRRA,
605 Rfree,
606 PermGroFac,
607 BoroCnstArt,
608 aXtraGrid,
609 vFuncBool,
610 CubicBool,
611):
612 """Solves one period of a consumption-saving model with idiosyncratic shocks to
613 permanent and transitory income, with one risk free asset and CRRA utility.
615 Parameters
616 ----------
617 solution_next : ConsumerSolution
618 The solution to next period's one period problem.
619 IncShkDstn : distribution.Distribution
620 A discrete approximation to the income process between the period being
621 solved and the one immediately following (in solution_next).
622 LivPrb : float
623 Survival probability; likelihood of being alive at the beginning of
624 the succeeding period.
625 DiscFac : float
626 Intertemporal discount factor for future utility.
627 CRRA : float
628 Coefficient of relative risk aversion.
629 Rfree : float
630 Risk free interest factor on end-of-period assets.
631 PermGroFac : float
632 Expected permanent income growth factor at the end of this period.
633 BoroCnstArt: float or None
634 Borrowing constraint for the minimum allowable assets to end the
635 period with. If it is less than the natural borrowing constraint,
636 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
637 rowing constraint.
638 aXtraGrid: np.array
639 Array of "extra" end-of-period asset values-- assets above the
640 absolute minimum acceptable level.
641 vFuncBool: boolean
642 An indicator for whether the value function should be computed and
643 included in the reported solution.
644 CubicBool: boolean
645 An indicator for whether the solver should use cubic or linear interpolation.
647 Returns
648 -------
649 solution_now : ConsumerSolution
650 Solution to this period's consumption-saving problem with income risk.
652 """
653 # Define the current period utility function and effective discount factor
654 uFunc = UtilityFuncCRRA(CRRA)
655 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
657 # Calculate the probability that we get the worst possible income draw
658 WorstIncPrb = calc_worst_inc_prob(IncShkDstn)
659 Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn)
660 hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, Ex_IncNext)
662 # Unpack next period's (marginal) value function
663 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
664 vPfuncNext = solution_next.vPfunc
665 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
667 # Calculate the minimum allowable value of money resources in this period
668 BoroCnstNat = calc_boro_const_nat(
669 solution_next.mNrmMin, IncShkDstn, Rfree, PermGroFac
670 )
671 # Set the minimum allowable (normalized) market resources based on the natural
672 # and artificial borrowing constraints
673 mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat)
675 # Update the bounding MPCs and PDV of human wealth:
676 PatFac = calc_patience_factor(Rfree, DiscFacEff, CRRA)
677 MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac)
678 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
679 # or artificial borrowing constraint actually binds
680 MPCmaxUnc = calc_mpc_max(
681 solution_next.MPCmax, WorstIncPrb, CRRA, PatFac, BoroCnstNat, BoroCnstArt
682 )
683 MPCmaxNow = 1.0 if BoroCnstNat < mNrmMinNow else MPCmaxUnc
685 cFuncLimitIntercept = MPCminNow * hNrmNow
686 cFuncLimitSlope = MPCminNow
688 # Define the borrowing-constrained consumption function
689 cFuncNowCnst = LinearInterp(
690 np.array([mNrmMinNow, mNrmMinNow + 1.0]),
691 np.array([0.0, 1.0]),
692 )
694 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
695 aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
697 # Calculate end-of-period marginal value of assets at each gridpoint
698 vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA)
699 EndOfPrdvP = vPfacEff * expected(
700 calc_vp_next,
701 IncShkDstn,
702 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPfuncNext),
703 )
705 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
706 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
707 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
709 # Limiting consumption is zero as m approaches mNrmMin
710 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
711 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
713 # Construct the consumption function as a cubic or linear spline interpolation
714 if CubicBool:
715 # Calculate end-of-period marginal marginal value of assets at each gridpoint
716 vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0)
717 EndOfPrdvPP = vPPfacEff * expected(
718 calc_vpp_next,
719 IncShkDstn,
720 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPPfuncNext),
721 )
722 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
723 MPC = dcda / (dcda + 1.0)
724 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxUnc)
726 # Construct the unconstrained consumption function as a cubic interpolation
727 cFuncNowUnc = CubicInterp(
728 m_for_interpolation,
729 c_for_interpolation,
730 MPC_for_interpolation,
731 cFuncLimitIntercept,
732 cFuncLimitSlope,
733 )
734 else:
735 # Construct the unconstrained consumption function as a linear interpolation
736 cFuncNowUnc = LinearInterp(
737 m_for_interpolation,
738 c_for_interpolation,
739 cFuncLimitIntercept,
740 cFuncLimitSlope,
741 )
743 # Combine the constrained and unconstrained functions into the true consumption function.
744 # LowerEnvelope should only be used when BoroCnstArt is True
745 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
747 # Make the marginal value function and the marginal marginal value function
748 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
750 # Define this period's marginal marginal value function
751 if CubicBool:
752 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
753 else:
754 vPPfuncNow = NullFunc() # Dummy object
756 # Construct this period's value function if requested
757 if vFuncBool:
758 # Calculate end-of-period value, its derivative, and their pseudo-inverse
759 EndOfPrdv = DiscFacEff * expected(
760 calc_v_next,
761 IncShkDstn,
762 args=(aNrmNow, Rfree, CRRA, PermGroFac, vFuncNext),
763 )
764 EndOfPrdvNvrs = uFunc.inv(
765 EndOfPrdv,
766 ) # value transformed through inverse utility
767 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
768 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
769 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
770 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
772 # Construct the end-of-period value function
773 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
774 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
775 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
777 # Compute expected value and marginal value on a grid of market resources
778 mNrm_temp = mNrmMinNow + aXtraGrid
779 cNrm_temp = cFuncNow(mNrm_temp)
780 aNrm_temp = mNrm_temp - cNrm_temp
781 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
782 vP_temp = uFunc.der(cNrm_temp)
784 # Construct the beginning-of-period value function
785 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
786 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
787 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
788 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
789 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
790 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
791 vNvrsFuncNow = CubicInterp(
792 mNrm_temp,
793 vNvrs_temp,
794 vNvrsP_temp,
795 MPCminNvrs * hNrmNow,
796 MPCminNvrs,
797 )
798 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
799 else:
800 vFuncNow = NullFunc() # Dummy object
802 # Create and return this period's solution
803 solution_now = ConsumerSolution(
804 cFunc=cFuncNow,
805 vFunc=vFuncNow,
806 vPfunc=vPfuncNow,
807 vPPfunc=vPPfuncNow,
808 mNrmMin=mNrmMinNow,
809 hNrm=hNrmNow,
810 MPCmin=MPCminNow,
811 MPCmax=MPCmaxNow,
812 )
813 return solution_now
816def solve_one_period_ConsKinkedR(
817 solution_next,
818 IncShkDstn,
819 LivPrb,
820 DiscFac,
821 CRRA,
822 Rboro,
823 Rsave,
824 PermGroFac,
825 BoroCnstArt,
826 aXtraGrid,
827 vFuncBool,
828 CubicBool,
829):
830 """Solves one period of a consumption-saving model with idiosyncratic shocks to
831 permanent and transitory income, with a risk free asset and CRRA utility.
832 In this variation, the interest rate on borrowing Rboro exceeds the interest
833 rate on saving Rsave.
835 Parameters
836 ----------
837 solution_next : ConsumerSolution
838 The solution to next period's one period problem.
839 IncShkDstn : distribution.Distribution
840 A discrete approximation to the income process between the period being
841 solved and the one immediately following (in solution_next).
842 LivPrb : float
843 Survival probability; likelihood of being alive at the beginning of
844 the succeeding period.
845 DiscFac : float
846 Intertemporal discount factor for future utility.
847 CRRA : float
848 Coefficient of relative risk aversion.
849 Rboro: float
850 Interest factor on assets between this period and the succeeding
851 period when assets are negative.
852 Rsave: float
853 Interest factor on assets between this period and the succeeding
854 period when assets are positive.
855 PermGroFac : float
856 Expected permanent income growth factor at the end of this period.
857 BoroCnstArt: float or None
858 Borrowing constraint for the minimum allowable assets to end the
859 period with. If it is less than the natural borrowing constraint,
860 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
861 rowing constraint.
862 aXtraGrid: np.array
863 Array of "extra" end-of-period asset values-- assets above the
864 absolute minimum acceptable level.
865 vFuncBool: boolean
866 An indicator for whether the value function should be computed and
867 included in the reported solution.
868 CubicBool: boolean
869 An indicator for whether the solver should use cubic or linear inter-
870 polation.
872 Returns
873 -------
874 solution_now : ConsumerSolution
875 Solution to this period's consumption-saving problem with income risk.
877 """
878 # Verifiy that there is actually a kink in the interest factor
879 assert Rboro >= Rsave, (
880 "Interest factor on debt less than interest factor on savings!"
881 )
882 # If the kink is in the wrong direction, code should break here. If there's
883 # no kink at all, then just use the ConsIndShockModel solver.
884 if Rboro == Rsave:
885 solution_now = solve_one_period_ConsIndShock(
886 solution_next,
887 IncShkDstn,
888 LivPrb,
889 DiscFac,
890 CRRA,
891 Rboro,
892 PermGroFac,
893 BoroCnstArt,
894 aXtraGrid,
895 vFuncBool,
896 CubicBool,
897 )
898 return solution_now
900 # Define the current period utility function and effective discount factor
901 uFunc = UtilityFuncCRRA(CRRA)
902 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
904 # Calculate the probability that we get the worst possible income draw
905 WorstIncPrb = calc_worst_inc_prob(IncShkDstn, use_infimum=False)
906 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
907 Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn)
908 hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext)
910 # Unpack next period's (marginal) value function
911 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
912 vPfuncNext = solution_next.vPfunc
913 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
915 # Calculate the minimum allowable value of money resources in this period
916 BoroCnstNat = calc_boro_const_nat(
917 solution_next.mNrmMin,
918 IncShkDstn,
919 Rboro,
920 PermGroFac,
921 use_infimum=False,
922 )
923 # Set the minimum allowable (normalized) market resources based on the natural
924 # and artificial borrowing constraints
925 mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat)
927 # Update the bounding MPCs and PDV of human wealth:
928 PatFacSave = calc_patience_factor(Rsave, DiscFacEff, CRRA)
929 PatFacBoro = calc_patience_factor(Rboro, DiscFacEff, CRRA)
930 MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave)
931 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
932 # or artificial borrowing constraint actually binds
933 MPCmaxUnc = calc_mpc_max(
934 solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro, BoroCnstNat, BoroCnstArt
935 )
936 MPCmaxNow = 1.0 if BoroCnstNat < mNrmMinNow else MPCmaxUnc
938 cFuncLimitIntercept = MPCminNow * hNrmNow
939 cFuncLimitSlope = MPCminNow
941 # Define the borrowing-constrained consumption function
942 cFuncNowCnst = LinearInterp(
943 np.array([mNrmMinNow, mNrmMinNow + 1.0]),
944 np.array([0.0, 1.0]),
945 )
947 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
948 aNrmNow = np.sort(
949 np.hstack((np.asarray(aXtraGrid) + mNrmMinNow, np.array([0.0, 1e-15]))),
950 )
952 # Make a 1D array of the interest factor at each asset gridpoint
953 Rfree = Rsave * np.ones_like(aNrmNow)
954 Rfree[aNrmNow <= 0] = Rboro
955 i_kink = np.argwhere(aNrmNow == 0.0)[0][0]
957 # Calculate end-of-period marginal value of assets at each gridpoint
958 vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA)
959 EndOfPrdvP = vPfacEff * expected(
960 calc_vp_next,
961 IncShkDstn,
962 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPfuncNext),
963 )
965 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
966 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
967 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
969 # Limiting consumption is zero as m approaches mNrmMin
970 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
971 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
973 # Construct the consumption function as a cubic or linear spline interpolation
974 if CubicBool:
975 # Calculate end-of-period marginal marginal value of assets at each gridpoint
976 vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0)
977 EndOfPrdvPP = vPPfacEff * expected(
978 calc_vpp_next,
979 IncShkDstn,
980 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPPfuncNext),
981 )
982 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
983 MPC = dcda / (dcda + 1.0)
984 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxUnc)
986 # Construct the unconstrained consumption function as a cubic interpolation
987 cFuncNowUnc = CubicInterp(
988 m_for_interpolation,
989 c_for_interpolation,
990 MPC_for_interpolation,
991 cFuncLimitIntercept,
992 cFuncLimitSlope,
993 )
994 # Adjust the coefficients on the kinked portion of the cFunc
995 cFuncNowUnc.coeffs[i_kink + 2] = [
996 c_for_interpolation[i_kink + 1],
997 m_for_interpolation[i_kink + 2] - m_for_interpolation[i_kink + 1],
998 0.0,
999 0.0,
1000 ]
1001 else:
1002 # Construct the unconstrained consumption function as a linear interpolation
1003 cFuncNowUnc = LinearInterp(
1004 m_for_interpolation,
1005 c_for_interpolation,
1006 cFuncLimitIntercept,
1007 cFuncLimitSlope,
1008 )
1010 # Combine the constrained and unconstrained functions into the true consumption function.
1011 # LowerEnvelope should only be used when BoroCnstArt is True
1012 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1014 # Make the marginal value function and the marginal marginal value function
1015 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1017 # Define this period's marginal marginal value function
1018 if CubicBool:
1019 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1020 else:
1021 vPPfuncNow = NullFunc() # Dummy object
1023 # Construct this period's value function if requested
1024 if vFuncBool:
1025 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1026 EndOfPrdv = DiscFacEff * expected(
1027 calc_v_next,
1028 IncShkDstn,
1029 args=(aNrmNow, Rfree, CRRA, PermGroFac, vFuncNext),
1030 )
1031 EndOfPrdvNvrs = uFunc.inv(
1032 EndOfPrdv,
1033 ) # value transformed through inverse utility
1034 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1035 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1036 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1037 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1039 # Construct the end-of-period value function
1040 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1041 EndOfPrdvNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1042 EndOfPrdvFunc = ValueFuncCRRA(EndOfPrdvNvrsFunc, CRRA)
1044 # Compute expected value and marginal value on a grid of market resources
1045 mNrm_temp = mNrmMinNow + aXtraGrid
1046 cNrm_temp = cFuncNow(mNrm_temp)
1047 aNrm_temp = mNrm_temp - cNrm_temp
1048 v_temp = uFunc(cNrm_temp) + EndOfPrdvFunc(aNrm_temp)
1049 vP_temp = uFunc.der(cNrm_temp)
1051 # Construct the beginning-of-period value function
1052 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1053 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1054 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1055 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1056 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
1057 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1058 vNvrsFuncNow = CubicInterp(
1059 mNrm_temp,
1060 vNvrs_temp,
1061 vNvrsP_temp,
1062 MPCminNvrs * hNrmNow,
1063 MPCminNvrs,
1064 )
1065 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1066 else:
1067 vFuncNow = NullFunc() # Dummy object
1069 # Create and return this period's solution
1070 solution_now = ConsumerSolution(
1071 cFunc=cFuncNow,
1072 vFunc=vFuncNow,
1073 vPfunc=vPfuncNow,
1074 vPPfunc=vPPfuncNow,
1075 mNrmMin=mNrmMinNow,
1076 hNrm=hNrmNow,
1077 MPCmin=MPCminNow,
1078 MPCmax=MPCmaxNow,
1079 )
1080 return solution_now
1083def make_basic_CRRA_solution_terminal(CRRA):
1084 """
1085 Construct the terminal period solution for a consumption-saving model with
1086 CRRA utility and only one state variable.
1088 Parameters
1089 ----------
1090 CRRA : float
1091 Coefficient of relative risk aversion. This is the only relevant parameter.
1093 Returns
1094 -------
1095 solution_terminal : ConsumerSolution
1096 Terminal period solution for someone with the given CRRA.
1097 """
1098 cFunc_terminal = LinearInterp([0.0, 1.0], [0.0, 1.0]) # c=m at t=T
1099 vFunc_terminal = ValueFuncCRRA(cFunc_terminal, CRRA)
1100 vPfunc_terminal = MargValueFuncCRRA(cFunc_terminal, CRRA)
1101 vPPfunc_terminal = MargMargValueFuncCRRA(cFunc_terminal, CRRA)
1102 solution_terminal = ConsumerSolution(
1103 cFunc=cFunc_terminal,
1104 vFunc=vFunc_terminal,
1105 vPfunc=vPfunc_terminal,
1106 vPPfunc=vPPfunc_terminal,
1107 mNrmMin=0.0,
1108 hNrm=0.0,
1109 MPCmin=1.0,
1110 MPCmax=1.0,
1111 )
1112 return solution_terminal
1115# ============================================================================
1116# == Classes for representing types of consumer agents (and things they do) ==
1117# ============================================================================
1119# Make a dictionary of constructors (very simply for perfect foresight model)
1120PerfForesightConsumerType_constructors_default = {
1121 "solution_terminal": make_basic_CRRA_solution_terminal,
1122 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1123 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1124}
1126# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1127PerfForesightConsumerType_kNrmInitDstn_default = {
1128 "kLogInitMean": -12.0, # Mean of log initial capital
1129 "kLogInitStd": 0.0, # Stdev of log initial capital
1130 "kNrmInitCount": 15, # Number of points in initial capital discretization
1131}
1133# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1134PerfForesightConsumerType_pLvlInitDstn_default = {
1135 "pLogInitMean": 0.0, # Mean of log permanent income
1136 "pLogInitStd": 0.0, # Stdev of log permanent income
1137 "pLvlInitCount": 15, # Number of points in initial capital discretization
1138}
1140# Make a dictionary to specify a perfect foresight consumer type
1141PerfForesightConsumerType_solving_defaults = {
1142 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1143 "cycles": 1, # Finite, non-cyclic model
1144 "T_cycle": 1, # Number of periods in the cycle for this agent type
1145 "pseudo_terminal": False, # Terminal period really does exist
1146 "constructors": PerfForesightConsumerType_constructors_default, # See dictionary above
1147 # PARAMETERS REQUIRED TO SOLVE THE MODEL
1148 "CRRA": 2.0, # Coefficient of relative risk aversion
1149 "Rfree": [1.03], # Interest factor on retained assets
1150 "DiscFac": 0.96, # Intertemporal discount factor
1151 "LivPrb": [0.98], # Survival probability after each period
1152 "PermGroFac": [1.01], # Permanent income growth factor
1153 "BoroCnstArt": None, # Artificial borrowing constraint
1154 "MaxKinks": 400, # Maximum number of grid points to allow in cFunc
1155}
1156PerfForesightConsumerType_simulation_defaults = {
1157 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1158 "AgentCount": 10000, # Number of agents of this type
1159 "T_age": None, # Age after which simulated agents are automatically killed
1160 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1161 # (The portion of PermGroFac attributable to aggregate productivity growth)
1162 # ADDITIONAL OPTIONAL PARAMETERS
1163 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1164 # (Forces Newborns to follow solution path of the agent they replaced if True)
1165}
1166PerfForesightConsumerType_defaults = {}
1167PerfForesightConsumerType_defaults.update(PerfForesightConsumerType_solving_defaults)
1168PerfForesightConsumerType_defaults.update(
1169 PerfForesightConsumerType_kNrmInitDstn_default
1170)
1171PerfForesightConsumerType_defaults.update(
1172 PerfForesightConsumerType_pLvlInitDstn_default
1173)
1174PerfForesightConsumerType_defaults.update(PerfForesightConsumerType_simulation_defaults)
1175init_perfect_foresight = PerfForesightConsumerType_defaults
1178class PerfForesightConsumerType(AgentType):
1179 r"""
1180 A perfect foresight consumer type who has no uncertainty other than mortality.
1181 Their problem is defined by a coefficient of relative risk aversion (:math:`\rho`), intertemporal
1182 discount factor (:math:`\beta`), interest factor (:math:`\mathsf{R}`), an optional artificial borrowing constraint (:math:`\underline{a}`)
1183 and time sequences of the permanent income growth rate (:math:`\Gamma`) and survival probability (:math:`1-\mathsf{D}`).
1184 Their assets and income are normalized by permanent income.
1186 .. math::
1187 \newcommand{\CRRA}{\rho}
1188 \newcommand{\DiePrb}{\mathsf{D}}
1189 \newcommand{\PermGroFac}{\Gamma}
1190 \newcommand{\Rfree}{\mathsf{R}}
1191 \newcommand{\DiscFac}{\beta}
1192 \begin{align*}
1193 v_t(m_t) &= \max_{c_t}u(c_t) + \DiscFac (1 - \DiePrb_{t+1}) \PermGroFac_{t+1}^{1-\CRRA} v_{t+1}(m_{t+1}), \\
1194 & \text{s.t.} \\
1195 a_t &= m_t - c_t, \\
1196 a_t &\geq \underline{a}, \\
1197 m_{t+1} &= \Rfree_{t+1} a_t/\PermGroFac_{t+1} + 1, \\
1198 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}
1199 \end{align*}
1202 Solving Parameters
1203 ------------------
1204 cycles: int
1205 0 specifies an infinite horizon model, 1 specifies a finite model.
1206 T_cycle: int
1207 Number of periods in the cycle for this agent type.
1208 CRRA: float, :math:`\rho`
1209 Coefficient of Relative Risk Aversion.
1210 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1211 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1212 DiscFac: float, :math:`\beta`
1213 Intertemporal discount factor.
1214 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1215 Survival probability after each period.
1216 PermGroFac: list[float], time varying, :math:`\Gamma`
1217 Permanent income growth factor.
1218 BoroCnstArt: float, :math:`\underline{a}`
1219 The minimum Asset/Perminant Income ratio, None to ignore.
1220 MaxKinks: int
1221 Maximum number of gridpoints to allow in cFunc.
1223 Simulation Parameters
1224 ---------------------
1225 AgentCount: int
1226 Number of agents of this kind that are created during simulations.
1227 T_age: int
1228 Age after which to automatically kill agents, None to ignore.
1229 T_sim: int, required for simulation
1230 Number of periods to simulate.
1231 track_vars: list[strings]
1232 List of variables that should be tracked when running the simulation.
1233 For this agent, the options are 'kNrm', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
1235 kNrm is beginning-of-period capital holdings (last period's assets)
1237 aLvl is the nominal asset level
1239 aNrm is the normalized assets
1241 bNrm is the normalized resources without this period's labor income
1243 cNrm is the normalized consumption
1245 mNrm is the normalized market resources
1247 pLvl is the permanent income level
1249 who_dies is the array of which agents died
1250 aNrmInitMean: float
1251 Mean of Log initial Normalized Assets.
1252 aNrmInitStd: float
1253 Std of Log initial Normalized Assets.
1254 pLvlInitMean: float
1255 Mean of Log initial permanent income.
1256 pLvlInitStd: float
1257 Std of Log initial permanent income.
1258 PermGroFacAgg: float
1259 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1260 PerfMITShk: boolean
1261 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1263 Attributes
1264 ----------
1265 solution: list[Consumer solution object]
1266 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1267 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1269 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1270 history: Dict[Array]
1271 Created by running the :func:`.simulate()` method.
1272 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1273 Visit :class:`HARK.core.AgentType.simulate` for more information.
1274 """
1276 solving_defaults = PerfForesightConsumerType_solving_defaults
1277 simulation_defaults = PerfForesightConsumerType_simulation_defaults
1279 default_ = {
1280 "params": PerfForesightConsumerType_defaults,
1281 "solver": solve_one_period_ConsPF,
1282 "model": "ConsPerfForesight.yaml",
1283 }
1285 time_vary_ = ["LivPrb", "PermGroFac", "Rfree"]
1286 time_inv_ = ["CRRA", "DiscFac", "MaxKinks", "BoroCnstArt"]
1287 state_vars = ["kNrm", "pLvl", "bNrm", "mNrm", "aNrm", "aLvl"]
1288 shock_vars_ = []
1289 distributions = ["kNrmInitDstn", "pLvlInitDstn"]
1291 def pre_solve(self):
1292 """
1293 Method that is run automatically just before solution by backward iteration.
1294 Solves the (trivial) terminal period and does a quick check on the borrowing
1295 constraint and MaxKinks attribute (only relevant in constrained, infinite
1296 horizon problems).
1297 """
1298 self.check_restrictions()
1299 self.construct("solution_terminal") # Solve the terminal period problem
1300 self.check_conditions(verbose=self.verbose)
1302 def post_solve(self):
1303 """
1304 Method that is run automatically at the end of a call to solve. Here, it
1305 simply calls calc_stable_points() if appropriate: an infinite horizon
1306 problem with a single repeated period in its cycle.
1308 Parameters
1309 ----------
1310 None
1312 Returns
1313 -------
1314 None
1315 """
1316 if (self.cycles == 0) and (self.T_cycle == 1):
1317 self.calc_stable_points()
1319 def check_restrictions(self):
1320 """
1321 A method to check that various restrictions are met for the model class.
1322 """
1323 if self.DiscFac < 0:
1324 raise ValueError("DiscFac is below zero with value: " + str(self.DiscFac))
1326 def initialize_sim(self):
1327 self.PermShkAggNow = self.PermGroFacAgg # This never changes during simulation
1328 self.state_now["PlvlAgg"] = 1.0
1329 super().initialize_sim()
1331 def sim_birth(self, which_agents):
1332 """
1333 Makes new consumers for the given indices. Initialized variables include aNrm and pLvl, as
1334 well as time variables t_age and t_cycle. Normalized assets and permanent income levels
1335 are drawn from lognormal distributions given by aNrmInitMean and aNrmInitStd (etc).
1337 Parameters
1338 ----------
1339 which_agents : np.array(Bool)
1340 Boolean array of size self.AgentCount indicating which agents should be "born".
1342 Returns
1343 -------
1344 None
1345 """
1346 # Get and store states for newly born agents
1347 N = np.sum(which_agents) # Number of new consumers to make
1348 self.state_now["aNrm"][which_agents] = self.kNrmInitDstn.draw(N)
1349 self.state_now["pLvl"][which_agents] = self.pLvlInitDstn.draw(N)
1350 self.state_now["pLvl"][which_agents] *= self.state_now["PlvlAgg"]
1351 self.t_age[which_agents] = 0 # How many periods since each agent was born
1353 # Because of the timing of the simulation system, kNrm gets written to
1354 # the *previous* period's aNrm after that aNrm has already been copied
1355 # to the history array (if it's being tracked). It will be loaded into
1356 # the simulation as kNrm, however, when the period is simulated.
1358 # If PerfMITShk not specified, let it be False
1359 if not hasattr(self, "PerfMITShk"):
1360 self.PerfMITShk = False
1361 if not self.PerfMITShk:
1362 # If True, Newborns inherit t_cycle of agent they replaced (i.e. t_cycles are not reset).
1363 self.t_cycle[which_agents] = 0
1364 # Which period of the cycle each agent is currently in
1366 def sim_death(self):
1367 """
1368 Determines which agents die this period and must be replaced. Uses the sequence in LivPrb
1369 to determine survival probabilities for each agent.
1371 Parameters
1372 ----------
1373 None
1375 Returns
1376 -------
1377 which_agents : np.array(bool)
1378 Boolean array of size AgentCount indicating which agents die.
1379 """
1380 # Determine who dies
1381 DiePrb_by_t_cycle = 1.0 - np.asarray(self.LivPrb)
1382 DiePrb = DiePrb_by_t_cycle[
1383 self.t_cycle - 1 if self.cycles == 1 else self.t_cycle
1384 ] # Time has already advanced, so look back one
1386 # In finite-horizon problems the previous line gives newborns the
1387 # survival probability of the last non-terminal period. This is okay,
1388 # however, since they will be instantly replaced by new newborns if
1389 # they die.
1390 # See: https://github.com/econ-ark/HARK/pull/981
1392 DeathShks = Uniform(seed=self.RNG.integers(0, 2**31 - 1)).draw(
1393 N=self.AgentCount
1394 )
1395 which_agents = DeathShks < DiePrb
1396 if self.T_age is not None: # Kill agents that have lived for too many periods
1397 too_old = self.t_age >= self.T_age
1398 which_agents = np.logical_or(which_agents, too_old)
1399 return which_agents
1401 def get_shocks(self):
1402 """
1403 Finds permanent and transitory income "shocks" for each agent this period. As this is a
1404 perfect foresight model, there are no stochastic shocks: PermShkNow = PermGroFac for each
1405 agent (according to their t_cycle) and TranShkNow = 1.0 for all agents.
1407 Parameters
1408 ----------
1409 None
1411 Returns
1412 -------
1413 None
1414 """
1415 PermGroFac = np.array(self.PermGroFac)
1416 # Cycle time has already been advanced
1417 self.shocks["PermShk"] = PermGroFac[self.t_cycle - 1]
1418 # self.shocks["PermShk"][self.t_cycle == 0] = 1. # Add this at some point
1419 self.shocks["TranShk"] = np.ones(self.AgentCount)
1421 def get_Rfree(self):
1422 """
1423 Returns an array of size self.AgentCount with Rfree in every entry.
1425 Parameters
1426 ----------
1427 None
1429 Returns
1430 -------
1431 RfreeNow : np.array
1432 Array of size self.AgentCount with risk free interest rate for each agent.
1433 """
1434 Rfree_array = np.array(self.Rfree)
1435 return Rfree_array[self.t_cycle]
1437 def transition(self):
1438 pLvlPrev = self.state_prev["pLvl"]
1439 kNrm = self.state_prev["aNrm"]
1440 RfreeNow = self.get_Rfree()
1442 # Calculate new states: normalized market resources and permanent income level
1443 # Updated permanent income level
1444 pLvlNow = pLvlPrev * self.shocks["PermShk"]
1445 # "Effective" interest factor on normalized assets
1446 ReffNow = RfreeNow / self.shocks["PermShk"]
1447 bNrmNow = ReffNow * kNrm # Bank balances before labor income
1448 # Market resources after income
1449 mNrmNow = bNrmNow + self.shocks["TranShk"]
1451 return kNrm, pLvlNow, bNrmNow, mNrmNow, None
1453 def get_controls(self):
1454 """
1455 Calculates consumption for each consumer of this type using the consumption functions.
1457 Parameters
1458 ----------
1459 None
1461 Returns
1462 -------
1463 None
1464 """
1465 cNrmNow = np.full(self.AgentCount, np.nan)
1466 MPCnow = np.full(self.AgentCount, np.nan)
1467 for t in np.unique(self.t_cycle):
1468 idx = self.t_cycle == t
1469 if np.any(idx):
1470 cNrmNow[idx], MPCnow[idx] = self.solution[t].cFunc.eval_with_derivative(
1471 self.state_now["mNrm"][idx]
1472 )
1473 self.controls["cNrm"] = cNrmNow
1475 # MPCnow is not really a control
1476 self.MPCnow = MPCnow
1478 def get_poststates(self):
1479 """
1480 Calculates end-of-period assets for each consumer of this type.
1482 Parameters
1483 ----------
1484 None
1486 Returns
1487 -------
1488 None
1489 """
1490 self.state_now["aNrm"] = self.state_now["mNrm"] - self.controls["cNrm"]
1491 self.state_now["aLvl"] = self.state_now["aNrm"] * self.state_now["pLvl"]
1492 # Update aggregate permanent productivity level
1493 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1495 def log_condition_result(self, name, result, message, verbose):
1496 """
1497 Records the result of one condition check in the attribute condition_report
1498 of the bilt dictionary, and in the message log.
1500 Parameters
1501 ----------
1502 name : string or None
1503 Name for the condition; if None, no test result is added to conditions.
1504 result : bool
1505 An indicator for whether the condition was passed.
1506 message : str
1507 The messages to record about the condition check.
1508 verbose : bool
1509 Indicator for whether verbose messages should be included in the report.
1510 """
1511 if name is not None:
1512 self.conditions[name] = result
1513 set_verbosity_level((4 - verbose) * 10)
1514 _log.info(message)
1515 self.bilt["conditions_report"] += message + "\n"
1517 def check_AIC(self, verbose=None):
1518 """
1519 Evaluate and report on the Absolute Impatience Condition.
1520 """
1521 name = "AIC"
1522 APFac = self.bilt["APFac"]
1523 result = APFac < 1.0
1525 messages = {
1526 True: f"APFac={APFac:.5f} : The Absolute Patience Factor satisfies the Absolute Impatience Condition (AIC) Þ < 1.",
1527 False: f"APFac={APFac:.5f} : The Absolute Patience Factor violates the Absolute Impatience Condition (AIC) Þ < 1.",
1528 }
1529 verbose = self.verbose if verbose is None else verbose
1530 self.log_condition_result(name, result, messages[result], verbose)
1532 def check_GICRaw(self, verbose=None):
1533 """
1534 Evaluate and report on the Growth Impatience Condition for the Perfect Foresight model.
1535 """
1536 name = "GICRaw"
1537 GPFacRaw = self.bilt["GPFacRaw"]
1538 result = GPFacRaw < 1.0
1540 messages = {
1541 True: f"GPFacRaw={GPFacRaw:.5f} : The Growth Patience Factor satisfies the Growth Impatience Condition (GICRaw) Þ/G < 1.",
1542 False: f"GPFacRaw={GPFacRaw:.5f} : The Growth Patience Factor violates the Growth Impatience Condition (GICRaw) Þ/G < 1.",
1543 }
1544 verbose = self.verbose if verbose is None else verbose
1545 self.log_condition_result(name, result, messages[result], verbose)
1547 def check_RIC(self, verbose=None):
1548 """
1549 Evaluate and report on the Return Impatience Condition.
1550 """
1551 name = "RIC"
1552 RPFac = self.bilt["RPFac"]
1553 result = RPFac < 1.0
1555 messages = {
1556 True: f"RPFac={RPFac:.5f} : The Return Patience Factor satisfies the Return Impatience Condition (RIC) Þ/R < 1.",
1557 False: f"RPFac={RPFac:.5f} : The Return Patience Factor violates the Return Impatience Condition (RIC) Þ/R < 1.",
1558 }
1559 verbose = self.verbose if verbose is None else verbose
1560 self.log_condition_result(name, result, messages[result], verbose)
1562 def check_FHWC(self, verbose=None):
1563 """
1564 Evaluate and report on the Finite Human Wealth Condition.
1565 """
1566 name = "FHWC"
1567 FHWFac = self.bilt["FHWFac"]
1568 result = FHWFac < 1.0
1570 messages = {
1571 True: f"FHWFac={FHWFac:.5f} : The Finite Human Wealth Factor satisfies the Finite Human Wealth Condition (FHWC) G/R < 1.",
1572 False: f"FHWFac={FHWFac:.5f} : The Finite Human Wealth Factor violates the Finite Human Wealth Condition (FHWC) G/R < 1.",
1573 }
1574 verbose = self.verbose if verbose is None else verbose
1575 self.log_condition_result(name, result, messages[result], verbose)
1577 def check_FVAC(self, verbose=None):
1578 """
1579 Evaluate and report on the Finite Value of Autarky Condition under perfect foresight.
1580 """
1581 name = "PFFVAC"
1582 PFVAFac = self.bilt["PFVAFac"]
1583 result = PFVAFac < 1.0
1585 messages = {
1586 True: f"PFVAFac={PFVAFac:.5f} : The Finite Value of Autarky Factor satisfies the Finite Value of Autarky Condition βG^(1-ρ) < 1.",
1587 False: f"PFVAFac={PFVAFac:.5f} : The Finite Value of Autarky Factor violates the Finite Value of Autarky Condition βG^(1-ρ) < 1.",
1588 }
1589 verbose = self.verbose if verbose is None else verbose
1590 self.log_condition_result(name, result, messages[result], verbose)
1592 def describe_parameters(self):
1593 """
1594 Make a string describing this instance's parameter values, including their
1595 representation in code and symbolically.
1597 Returns
1598 -------
1599 param_desc : str
1600 Description of parameters as a unicode string.
1601 """
1602 params_to_describe = [
1603 # [name, description, symbol, time varying]
1604 ["DiscFac", "intertemporal discount factor", "β", False],
1605 ["Rfree", "risk free interest factor", "R", True],
1606 ["PermGroFac", "permanent income growth factor", "G", True],
1607 ["CRRA", "coefficient of relative risk aversion", "ρ", False],
1608 ["LivPrb", "survival probability", "ℒ", True],
1609 ["APFac", "absolute patience factor", "Þ=(βℒR)^(1/ρ)", False],
1610 ]
1612 param_desc = ""
1613 for j in range(len(params_to_describe)):
1614 this_entry = params_to_describe[j]
1615 if this_entry[3]:
1616 val = getattr(self, this_entry[0])[0]
1617 else:
1618 try:
1619 val = getattr(self, this_entry[0])
1620 except:
1621 val = self.bilt[this_entry[0]]
1622 this_line = (
1623 this_entry[2]
1624 + f"={val:.5f} : "
1625 + this_entry[1]
1626 + " ("
1627 + this_entry[0]
1628 + ")\n"
1629 )
1630 param_desc += this_line
1632 return param_desc
1634 def calc_limiting_values(self):
1635 """
1636 Compute various scalar values that are relevant to characterizing the
1637 solution to an infinite horizon problem. This method should only be called
1638 when T_cycle=1 and cycles=0, otherwise the values generated are meaningless.
1639 This method adds the following values to the instance in the dictionary
1640 attribute called bilt.
1642 APFac : Absolute Patience Factor
1643 GPFacRaw : Growth Patience Factor
1644 FHWFac : Finite Human Wealth Factor
1645 RPFac : Return Patience Factor
1646 PFVAFac : Perfect Foresight Value of Autarky Factor
1647 cNrmPDV : Present Discounted Value of Autarky Consumption
1648 MPCmin : Limiting minimum MPC as market resources go to infinity
1649 MPCmax : Limiting maximum MPC as market resources approach minimum level.
1650 hNrm : Human wealth divided by permanent income.
1651 Delta_mNrm_ZeroFunc : Linear consumption function where expected change in market resource ratio is zero
1652 BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income
1654 Returns
1655 -------
1656 None
1657 """
1658 aux_dict = self.bilt
1659 aux_dict["APFac"] = (self.Rfree[0] * self.DiscFac * self.LivPrb[0]) ** (
1660 1 / self.CRRA
1661 )
1662 aux_dict["GPFacRaw"] = aux_dict["APFac"] / self.PermGroFac[0]
1663 aux_dict["FHWFac"] = self.PermGroFac[0] / self.Rfree[0]
1664 aux_dict["RPFac"] = aux_dict["APFac"] / self.Rfree[0]
1665 aux_dict["PFVAFac"] = (self.DiscFac * self.LivPrb[0]) * self.PermGroFac[0] ** (
1666 1.0 - self.CRRA
1667 )
1668 aux_dict["cNrmPDV"] = 1.0 / (1.0 - aux_dict["RPFac"])
1669 aux_dict["MPCmin"] = np.maximum(1.0 - aux_dict["RPFac"], 0.0)
1670 constrained = (
1671 hasattr(self, "BoroCnstArt")
1672 and (self.BoroCnstArt is not None)
1673 and (self.BoroCnstArt > -np.inf)
1674 )
1676 if constrained:
1677 aux_dict["MPCmax"] = 1.0
1678 else:
1679 aux_dict["MPCmax"] = aux_dict["MPCmin"]
1680 if aux_dict["FHWFac"] < 1.0:
1681 aux_dict["hNrm"] = 1.0 / (1.0 - aux_dict["FHWFac"])
1682 else:
1683 aux_dict["hNrm"] = np.inf
1685 # Generate the "Delta m = 0" function, which is used to find target market resources
1686 Ex_Rnrm = self.Rfree[0] / self.PermGroFac[0]
1687 aux_dict["Delta_mNrm_ZeroFunc"] = (
1688 lambda m: (1.0 - 1.0 / Ex_Rnrm) * m + 1.0 / Ex_Rnrm
1689 )
1691 # Generate the "E[M_tp1 / M_t] = G" function, which is used to find balanced growth market resources
1692 PF_Rnrm = self.Rfree[0] / self.PermGroFac[0]
1693 aux_dict["BalGroFunc"] = lambda m: (1.0 - 1.0 / PF_Rnrm) * m + 1.0 / PF_Rnrm
1695 self.bilt = aux_dict
1697 def check_conditions(self, verbose=None):
1698 """
1699 This method checks whether the instance's type satisfies the
1700 Absolute Impatience Condition (AIC), the Return Impatience Condition (RIC),
1701 the Finite Human Wealth Condition (FHWC), the perfect foresight model's
1702 Growth Impatience Condition (GICRaw) and Perfect Foresight Finite Value
1703 of Autarky Condition (FVACPF). Depending on the configuration of parameter
1704 values, somecombination of these conditions must be satisfied in order
1705 for the problem to have a nondegenerate solution. To check which conditions
1706 are required, in the verbose mode a reference to the relevant theoretical
1707 literature is made.
1709 Parameters
1710 ----------
1711 verbose : boolean
1712 Specifies different levels of verbosity of feedback. When False, it
1713 only reports whether the instance's type fails to satisfy a particular
1714 condition. When True, it reports all results, i.e. the factor values
1715 for all conditions.
1717 Returns
1718 -------
1719 None
1720 """
1721 self.conditions = {}
1722 self.bilt["conditions_report"] = ""
1723 self.degenerate = False
1724 verbose = self.verbose if verbose is None else verbose
1726 # This method only checks for the conditions for infinite horizon models
1727 # with a 1 period cycle. If these conditions are not met, we exit early.
1728 if self.cycles != 0 or self.T_cycle > 1:
1729 trivial_message = "No conditions report was produced because this functionality is only supported for infinite horizon models with a cycle length of 1."
1730 self.log_condition_result(None, None, trivial_message, verbose)
1731 if not self.quiet:
1732 _log.info(self.bilt["conditions_report"])
1733 return
1735 # Calculate some useful quantities that will be used in the condition checks
1736 self.calc_limiting_values()
1737 param_desc = self.describe_parameters()
1738 self.log_condition_result(None, None, param_desc, verbose)
1740 # Check individual conditions and add their results to the report
1741 self.check_AIC(verbose)
1742 self.check_RIC(verbose)
1743 self.check_GICRaw(verbose)
1744 self.check_FVAC(verbose)
1745 self.check_FHWC(verbose)
1746 constrained = (
1747 hasattr(self, "BoroCnstArt")
1748 and (self.BoroCnstArt is not None)
1749 and (self.BoroCnstArt > -np.inf)
1750 )
1752 # Exit now if verbose output was not requested.
1753 if not verbose:
1754 if not self.quiet:
1755 _log.info(self.bilt["conditions_report"])
1756 return
1758 # Report on the degeneracy of the consumption function solution
1759 if not constrained:
1760 if self.conditions["FHWC"]:
1761 RIC_message = "\nBecause the FHWC is satisfied, the solution is not c(m)=Infinity."
1762 if self.conditions["RIC"]:
1763 RIC_message += " Because the RIC is also satisfied, the solution is also not c(m)=0 for all m, so a non-degenerate linear solution exists."
1764 degenerate = False
1765 else:
1766 RIC_message += " However, because the RIC is violated, the solution is degenerate at c(m) = 0 for all m."
1767 degenerate = True
1768 else:
1769 RIC_message = "\nBecause the FHWC condition is violated and the consumer is not constrained, the solution is degenerate at c(m)=Infinity."
1770 degenerate = True
1771 else:
1772 if self.conditions["RIC"]:
1773 RIC_message = "\nBecause the RIC is satisfied and the consumer is constrained, the solution is not c(m)=0 for all m."
1774 if self.conditions["GICRaw"]:
1775 RIC_message += " Because the GICRaw is also satisfied, the solution is non-degenerate. It is piecewise linear with an infinite number of kinks, approaching the unconstrained solution as m goes to infinity."
1776 degenerate = False
1777 else:
1778 RIC_message += " Because the GICRaw is violated, the solution is non-degenerate. It is piecewise linear with a single kink at some 0 < m < 1; it equals the unconstrained solution above that kink point and has c(m) = m below it."
1779 degenerate = False
1780 else:
1781 if self.conditions["GICRaw"]:
1782 RIC_message = "\nBecause the RIC is violated but the GIC is satisfied, the FHWC is necessarily also violated. In this case, the consumer's pathological patience is offset by his infinite human wealth, against which he cannot borrow arbitrarily; a non-degenerate solution exists."
1783 degenerate = False
1784 else:
1785 RIC_message = "\nBecause the RIC is violated but the FHWC is satisfied, the solution is degenerate at c(m)=0 for all m."
1786 degenerate = True
1787 self.log_condition_result(None, None, RIC_message, verbose)
1789 if (
1790 degenerate
1791 ): # All of the other checks are meaningless if the solution is degenerate
1792 if not self.quiet:
1793 _log.info(self.bilt["conditions_report"])
1794 return
1796 # Report on the consequences of the Absolute Impatience Condition
1797 if self.conditions["AIC"]:
1798 AIC_message = "\nBecause the AIC is satisfied, the absolute amount of consumption is expected to fall over time."
1799 else:
1800 AIC_message = "\nBecause the AIC is violated, the absolute amount of consumption is expected to grow over time."
1801 self.log_condition_result(None, None, AIC_message, verbose)
1803 # Report on the consequences of the Growth Impatience Condition
1804 if self.conditions["GICRaw"]:
1805 GIC_message = "\nBecause the GICRaw is satisfed, the ratio of individual wealth to permanent income is expected to fall indefinitely."
1806 elif self.conditions["FHWC"]:
1807 GIC_message = "\nBecause the GICRaw is violated but the FHWC is satisfied, the ratio of individual wealth to permanent income is expected to rise toward infinity."
1808 else:
1809 pass # pragma: nocover
1810 # This can never be reached! If GICRaw and FHWC both fail, then the RIC also fails, and we would have exited by this point.
1811 self.log_condition_result(None, None, GIC_message, verbose)
1813 if not self.quiet:
1814 _log.info(self.bilt["conditions_report"])
1816 def calc_stable_points(self, force=False):
1817 """
1818 If the problem is one that satisfies the conditions required for target ratios of different
1819 variables to permanent income to exist, and has been solved to within the self-defined
1820 tolerance, this method calculates the target values of market resources.
1822 Parameters
1823 ----------
1824 force : bool
1825 Indicator for whether the method should be forced to be run even if
1826 the agent seems to be the wrong type. Default is False.
1828 Returns
1829 -------
1830 None
1831 """
1832 # Child classes should not run this method
1833 is_perf_foresight = type(self) is PerfForesightConsumerType
1834 is_ind_shock = type(self) is IndShockConsumerType
1835 if not (is_perf_foresight or is_ind_shock or force):
1836 return
1838 infinite_horizon = self.cycles == 0
1839 single_period = self.T_cycle == 1
1840 if not infinite_horizon:
1841 raise ValueError(
1842 "The calc_stable_points method works only for infinite horizon models."
1843 )
1844 if not single_period:
1845 raise ValueError(
1846 "The calc_stable_points method works only with a single infinitely repeated period."
1847 )
1848 if not hasattr(self, "conditions"):
1849 raise ValueError(
1850 "The check_conditions method must be run before the calc_stable_points method."
1851 )
1852 if not hasattr(self, "solution"):
1853 raise ValueError(
1854 "The solve method must be run before the calc_stable_points method."
1855 )
1857 # Extract balanced growth and delta m_t+1 = 0 functions
1858 BalGroFunc = self.bilt["BalGroFunc"]
1859 Delta_mNrm_ZeroFunc = self.bilt["Delta_mNrm_ZeroFunc"]
1861 # If the GICRaw holds, then there is a balanced growth market resources ratio
1862 if self.conditions["GICRaw"]:
1863 cFunc = self.solution[0].cFunc
1864 func_to_zero = lambda m: BalGroFunc(m) - cFunc(m)
1865 m0 = 1.0
1866 try:
1867 mNrmStE = newton(func_to_zero, m0)
1868 except:
1869 mNrmStE = np.nan
1871 # A target level of assets *might* exist even if the GICMod fails, so check no matter what
1872 func_to_zero = lambda m: Delta_mNrm_ZeroFunc(m) - cFunc(m)
1873 m0 = 1.0 if np.isnan(mNrmStE) else mNrmStE
1874 try:
1875 mNrmTrg = newton(func_to_zero, m0, maxiter=200)
1876 except:
1877 mNrmTrg = np.nan
1878 else:
1879 mNrmStE = np.nan
1880 mNrmTrg = np.nan
1882 self.solution[0].mNrmStE = mNrmStE
1883 self.solution[0].mNrmTrg = mNrmTrg
1884 self.bilt["mNrmStE"] = mNrmStE
1885 self.bilt["mNrmTrg"] = mNrmTrg
1888###############################################################################
1890# Make a dictionary of constructors for the idiosyncratic income shocks model
1891IndShockConsumerType_constructors_default = {
1892 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1893 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1894 "IncShkDstn": construct_lognormal_income_process_unemployment,
1895 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1896 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1897 "aXtraGrid": make_assets_grid,
1898 "solution_terminal": make_basic_CRRA_solution_terminal,
1899}
1901# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1902IndShockConsumerType_kNrmInitDstn_default = {
1903 "kLogInitMean": -12.0, # Mean of log initial capital
1904 "kLogInitStd": 0.0, # Stdev of log initial capital
1905 "kNrmInitCount": 15, # Number of points in initial capital discretization
1906}
1908# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1909IndShockConsumerType_pLvlInitDstn_default = {
1910 "pLogInitMean": 0.0, # Mean of log permanent income
1911 "pLogInitStd": 0.0, # Stdev of log permanent income
1912 "pLvlInitCount": 15, # Number of points in initial capital discretization
1913}
1915# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
1916IndShockConsumerType_IncShkDstn_default = {
1917 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1918 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1919 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1920 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1921 "UnempPrb": 0.05, # Probability of unemployment while working
1922 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1923 "T_retire": 0, # Period of retirement (0 --> no retirement)
1924 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1925 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1926}
1928# Default parameters to make aXtraGrid using make_assets_grid
1929IndShockConsumerType_aXtraGrid_default = {
1930 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1931 "aXtraMax": 20, # Maximum end-of-period "assets above minimum" value
1932 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
1933 "aXtraCount": 48, # Number of points in the grid of "assets above minimum"
1934 "aXtraExtra": None, # Additional other values to add in grid (optional)
1935}
1937# Make a dictionary to specify an idiosyncratic income shocks consumer type
1938IndShockConsumerType_solving_default = {
1939 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1940 "cycles": 1, # Finite, non-cyclic model
1941 "T_cycle": 1, # Number of periods in the cycle for this agent type
1942 "pseudo_terminal": False, # Terminal period really does exist
1943 "constructors": IndShockConsumerType_constructors_default, # See dictionary above
1944 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
1945 "CRRA": 2.0, # Coefficient of relative risk aversion
1946 "Rfree": [1.03], # Interest factor on retained assets
1947 "DiscFac": 0.96, # Intertemporal discount factor
1948 "LivPrb": [0.98], # Survival probability after each period
1949 "PermGroFac": [1.01], # Permanent income growth factor
1950 "BoroCnstArt": 0.0, # Artificial borrowing constraint
1951 "vFuncBool": False, # Whether to calculate the value function during solution
1952 "CubicBool": False, # Whether to use cubic spline interpolation when True
1953 # (Uses linear spline interpolation for cFunc when False)
1954}
1955IndShockConsumerType_simulation_default = {
1956 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1957 "AgentCount": 10000, # Number of agents of this type
1958 "T_age": None, # Age after which simulated agents are automatically killed
1959 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1960 # (The portion of PermGroFac attributable to aggregate productivity growth)
1961 "NewbornTransShk": False, # Whether Newborns have transitory shock
1962 # ADDITIONAL OPTIONAL PARAMETERS
1963 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1964 # (Forces Newborns to follow solution path of the agent they replaced if True)
1965 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1966}
1968IndShockConsumerType_defaults = {}
1969IndShockConsumerType_defaults.update(IndShockConsumerType_IncShkDstn_default)
1970IndShockConsumerType_defaults.update(IndShockConsumerType_kNrmInitDstn_default)
1971IndShockConsumerType_defaults.update(IndShockConsumerType_pLvlInitDstn_default)
1972IndShockConsumerType_defaults.update(IndShockConsumerType_aXtraGrid_default)
1973IndShockConsumerType_defaults.update(IndShockConsumerType_solving_default)
1974IndShockConsumerType_defaults.update(IndShockConsumerType_simulation_default)
1975init_idiosyncratic_shocks = IndShockConsumerType_defaults # Here so that other models which use the old convention don't break
1978class IndShockConsumerType(PerfForesightConsumerType):
1979 r"""
1980 A consumer type with idiosyncratic shocks to permanent and transitory income.
1981 Their problem is defined by a sequence of income distributions, survival probabilities
1982 (:math:`1-\mathsf{D}`), and permanent income growth rates (:math:`\Gamma`), as well
1983 as time invariant values for risk aversion (:math:`\rho`), discount factor (:math:`\beta`),
1984 the interest rate (:math:`\mathsf{R}`), the grid of end-of-period assets, and an artificial
1985 borrowing constraint (:math:`\underline{a}`).
1987 .. math::
1988 \newcommand{\CRRA}{\rho}
1989 \newcommand{\DiePrb}{\mathsf{D}}
1990 \newcommand{\PermGroFac}{\Gamma}
1991 \newcommand{\Rfree}{\mathsf{R}}
1992 \newcommand{\DiscFac}{\beta}
1993 \begin{align*}
1994 v_t(m_t) &= \max_{c_t}u(c_t) + \DiscFac (1 - \DiePrb_{t+1}) \mathbb{E}_{t} \left[ (\PermGroFac_{t+1} \psi_{t+1})^{1-\CRRA} v_{t+1}(m_{t+1}) \right], \\
1995 & \text{s.t.} \\
1996 a_t &= m_t - c_t, \\
1997 a_t &\geq \underline{a}, \\
1998 m_{t+1} &= a_t \Rfree_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
1999 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
2000 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1, \\
2001 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}
2002 \end{align*}
2005 Constructors
2006 ------------
2007 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2008 The agent's income shock distributions.
2010 It's default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2011 aXtraGrid: Constructor
2012 The agent's asset grid.
2014 It's default constructor is :func:`HARK.utilities.make_assets_grid`
2016 Solving Parameters
2017 ------------------
2018 cycles: int
2019 0 specifies an infinite horizon model, 1 specifies a finite model.
2020 T_cycle: int
2021 Number of periods in the cycle for this agent type.
2022 CRRA: float, :math:`\rho`
2023 Coefficient of Relative Risk Aversion.
2024 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
2025 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
2026 DiscFac: float, :math:`\beta`
2027 Intertemporal discount factor.
2028 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2029 Survival probability after each period.
2030 PermGroFac: list[float], time varying, :math:`\Gamma`
2031 Permanent income growth factor.
2032 BoroCnstArt: float, :math:`\underline{a}`
2033 The minimum Asset/Perminant Income ratio, None to ignore.
2034 vFuncBool: bool
2035 Whether to calculate the value function during solution.
2036 CubicBool: bool
2037 Whether to use cubic spline interpoliation.
2039 Simulation Parameters
2040 ---------------------
2041 AgentCount: int
2042 Number of agents of this kind that are created during simulations.
2043 T_age: int
2044 Age after which to automatically kill agents, None to ignore.
2045 T_sim: int, required for simulation
2046 Number of periods to simulate.
2047 track_vars: list[strings]
2048 List of variables that should be tracked when running the simulation.
2049 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2051 PermShk is the agent's permanent income shock
2053 TranShk is the agent's transitory income shock
2055 aLvl is the nominal asset level
2057 aNrm is the normalized assets
2059 bNrm is the normalized resources without this period's labor income
2061 cNrm is the normalized consumption
2063 mNrm is the normalized market resources
2065 pLvl is the permanent income level
2067 who_dies is the array of which agents died
2068 aNrmInitMean: float
2069 Mean of Log initial Normalized Assets.
2070 aNrmInitStd: float
2071 Std of Log initial Normalized Assets.
2072 pLvlInitMean: float
2073 Mean of Log initial permanent income.
2074 pLvlInitStd: float
2075 Std of Log initial permanent income.
2076 PermGroFacAgg: float
2077 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2078 PerfMITShk: boolean
2079 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2080 NewbornTransShk: boolean
2081 Whether Newborns have transitory shock.
2083 Attributes
2084 ----------
2085 solution: list[Consumer solution object]
2086 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2087 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2089 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2090 history: Dict[Array]
2091 Created by running the :func:`.simulate()` method.
2092 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2093 Visit :class:`HARK.core.AgentType.simulate` for more information.
2094 """
2096 IncShkDstn_defaults = IndShockConsumerType_IncShkDstn_default
2097 aXtraGrid_defaults = IndShockConsumerType_aXtraGrid_default
2098 solving_defaults = IndShockConsumerType_solving_default
2099 simulation_defaults = IndShockConsumerType_simulation_default
2100 default_ = {
2101 "params": IndShockConsumerType_defaults,
2102 "solver": solve_one_period_ConsIndShock,
2103 "model": "ConsIndShock.yaml",
2104 }
2106 time_inv_ = PerfForesightConsumerType.time_inv_ + [
2107 "vFuncBool",
2108 "CubicBool",
2109 "aXtraGrid",
2110 ]
2111 time_vary_ = PerfForesightConsumerType.time_vary_ + [
2112 "IncShkDstn",
2113 "PermShkDstn",
2114 "TranShkDstn",
2115 ]
2116 # This is in the PerfForesight model but not ConsIndShock
2117 time_inv_.remove("MaxKinks")
2118 shock_vars_ = ["PermShk", "TranShk"]
2119 distributions = [
2120 "IncShkDstn",
2121 "PermShkDstn",
2122 "TranShkDstn",
2123 "kNrmInitDstn",
2124 "pLvlInitDstn",
2125 ]
2127 def update_income_process(self):
2128 self.update("IncShkDstn", "PermShkDstn", "TranShkDstn")
2130 def get_shocks(self):
2131 """
2132 Gets permanent and transitory income shocks for this period. Samples from IncShkDstn for
2133 each period in the cycle.
2135 Parameters
2136 ----------
2137 NewbornTransShk : boolean, optional
2138 Whether Newborns have transitory shock. The default is False.
2140 Returns
2141 -------
2142 None
2143 """
2144 # Whether Newborns have transitory shock. The default is False.
2145 NewbornTransShk = self.NewbornTransShk
2147 PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
2148 TranShkNow = np.zeros(self.AgentCount)
2149 newborn = self.t_age == 0
2150 for t in np.unique(self.t_cycle):
2151 idx = self.t_cycle == t
2153 # temporary, see #1022
2154 if self.cycles == 1:
2155 t = t - 1
2157 N = np.sum(idx)
2158 if N > 0:
2159 # set current income distribution
2160 IncShkDstnNow = self.IncShkDstn[t]
2161 # and permanent growth factor
2162 PermGroFacNow = self.PermGroFac[t]
2163 # Get random draws of income shocks from the discrete distribution
2164 IncShks = IncShkDstnNow.draw(N)
2166 PermShkNow[idx] = (
2167 IncShks[0, :] * PermGroFacNow
2168 ) # permanent "shock" includes expected growth
2169 TranShkNow[idx] = IncShks[1, :]
2171 # That procedure used the *last* period in the sequence for newborns, but that's not right
2172 # Redraw shocks for newborns, using the *first* period in the sequence. Approximation.
2173 N = np.sum(newborn)
2174 if N > 0:
2175 idx = newborn
2176 # set current income distribution
2177 IncShkDstnNow = self.IncShkDstn[0]
2178 PermGroFacNow = self.PermGroFac[0] # and permanent growth factor
2180 # Get random draws of income shocks from the discrete distribution
2181 EventDraws = IncShkDstnNow.draw_events(N)
2182 PermShkNow[idx] = (
2183 IncShkDstnNow.atoms[0][EventDraws] * PermGroFacNow
2184 ) # permanent "shock" includes expected growth
2185 TranShkNow[idx] = IncShkDstnNow.atoms[1][EventDraws]
2187 # Whether Newborns have transitory shock. The default is False.
2188 if not NewbornTransShk:
2189 TranShkNow[newborn] = 1.0
2191 # Store the shocks in self
2192 self.shocks["PermShk"] = PermShkNow
2193 self.shocks["TranShk"] = TranShkNow
2195 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
2196 """
2197 Creates a "normalized Euler error" function for this instance, mapping
2198 from market resources to "consumption error per dollar of consumption."
2199 Stores result in attribute eulerErrorFunc as an interpolated function.
2200 Has option to use approximate income distribution stored in self.IncShkDstn
2201 or to use a (temporary) very dense approximation.
2203 Only works on (one period) infinite horizon models at this time, will
2204 be generalized later.
2206 Parameters
2207 ----------
2208 mMax : float
2209 Maximum normalized market resources for the Euler error function.
2210 approx_inc_dstn : Boolean
2211 Indicator for whether to use the approximate discrete income distri-
2212 bution stored in self.IncShkDstn[0], or to use a very accurate
2213 discrete approximation instead. When True, uses approximation in
2214 IncShkDstn; when False, makes and uses a very dense approximation.
2216 Returns
2217 -------
2218 None
2220 Notes
2221 -----
2222 This method is not used by any other code in the library. Rather, it is here
2223 for expository and benchmarking purposes.
2224 """
2225 # Get the income distribution (or make a very dense one)
2226 if approx_inc_dstn:
2227 IncShkDstn = self.IncShkDstn[0]
2228 else:
2229 TranShkDstn = MeanOneLogNormal(sigma=self.TranShkStd[0]).discretize(
2230 N=200,
2231 method="equiprobable",
2232 tail_N=50,
2233 tail_order=1.3,
2234 tail_bound=[0.05, 0.95],
2235 )
2236 TranShkDstn = add_discrete_outcome_constant_mean(
2237 TranShkDstn, p=self.UnempPrb, x=self.IncUnemp
2238 )
2239 PermShkDstn = MeanOneLogNormal(sigma=self.PermShkStd[0]).discretize(
2240 N=200,
2241 method="equiprobable",
2242 tail_N=50,
2243 tail_order=1.3,
2244 tail_bound=[0.05, 0.95],
2245 )
2246 IncShkDstn = combine_indep_dstns(PermShkDstn, TranShkDstn)
2248 # Make a grid of market resources
2249 mNowMin = self.solution[0].mNrmMin + 10 ** (
2250 -15
2251 ) # add tiny bit to get around 0/0 problem
2252 mNowMax = mMax
2253 mNowGrid = np.linspace(mNowMin, mNowMax, 1000)
2255 # Get the consumption function this period and the marginal value function
2256 # for next period. Note that this part assumes a one period cycle.
2257 cFuncNow = self.solution[0].cFunc
2258 vPfuncNext = self.solution[0].vPfunc
2260 # Calculate consumption this period at each gridpoint (and assets)
2261 cNowGrid = cFuncNow(mNowGrid)
2262 aNowGrid = mNowGrid - cNowGrid
2264 # Tile the grids for fast computation
2265 ShkCount = IncShkDstn.pmv.size
2266 aCount = aNowGrid.size
2267 aNowGrid_tiled = np.tile(aNowGrid, (ShkCount, 1))
2268 PermShkVals_tiled = (np.tile(IncShkDstn.atoms[0], (aCount, 1))).transpose()
2269 TranShkVals_tiled = (np.tile(IncShkDstn.atoms[1], (aCount, 1))).transpose()
2270 ShkPrbs_tiled = (np.tile(IncShkDstn.pmv, (aCount, 1))).transpose()
2272 # Calculate marginal value next period for each gridpoint and each shock
2273 mNextArray = (
2274 self.Rfree[0] / (self.PermGroFac[0] * PermShkVals_tiled) * aNowGrid_tiled
2275 + TranShkVals_tiled
2276 )
2277 vPnextArray = vPfuncNext(mNextArray)
2279 # Calculate expected marginal value and implied optimal consumption
2280 ExvPnextGrid = (
2281 self.DiscFac
2282 * self.Rfree[0]
2283 * self.LivPrb[0]
2284 * self.PermGroFac[0] ** (-self.CRRA)
2285 * np.sum(
2286 PermShkVals_tiled ** (-self.CRRA) * vPnextArray * ShkPrbs_tiled, axis=0
2287 )
2288 )
2289 cOptGrid = ExvPnextGrid ** (
2290 -1.0 / self.CRRA
2291 ) # This is the 'Endogenous Gridpoints' step
2293 # Calculate Euler error and store an interpolated function
2294 EulerErrorNrmGrid = (cNowGrid - cOptGrid) / cOptGrid
2295 eulerErrorFunc = LinearInterp(mNowGrid, EulerErrorNrmGrid)
2296 self.eulerErrorFunc = eulerErrorFunc
2298 def pre_solve(self):
2299 self.check_restrictions()
2300 self.construct("solution_terminal")
2301 if not self.quiet:
2302 self.check_conditions(verbose=self.verbose)
2304 def describe_parameters(self):
2305 """
2306 Generate a string describing the primitive model parameters that will
2307 be used to calculating limiting values and factors.
2309 Parameters
2310 ----------
2311 None
2313 Returns
2314 -------
2315 param_desc : str
2316 Description of primitive parameters.
2317 """
2318 # Get parameter description from the perfect foresight model
2319 param_desc = super().describe_parameters()
2321 # Make a new entry for weierstrass-p (the weird formatting here is to
2322 # make it easier to adapt into the style of the superclass if we add more
2323 # parameter reports later)
2324 this_entry = [
2325 "WorstPrb",
2326 "probability of worst income shock realization",
2327 "℘",
2328 False,
2329 ]
2330 try:
2331 val = getattr(self, this_entry[0])
2332 except:
2333 val = self.bilt[this_entry[0]]
2334 this_line = (
2335 this_entry[2]
2336 + f"={val:.5f} : "
2337 + this_entry[1]
2338 + " ("
2339 + this_entry[0]
2340 + ")\n"
2341 )
2343 # Add in the new entry and return it
2344 param_desc += this_line
2345 return param_desc
2347 def calc_limiting_values(self):
2348 """
2349 Compute various scalar values that are relevant to characterizing the
2350 solution to an infinite horizon problem. This method should only be called
2351 when T_cycle=1 and cycles=0, otherwise the values generated are meaningless.
2352 This method adds the following values to this instance in the dictionary
2353 attribute called bilt.
2355 APFac : Absolute Patience Factor
2356 GPFacRaw : Growth Patience Factor
2357 GPFacMod : Risk-Modified Growth Patience Factor
2358 GPFacLiv : Mortality-Adjusted Growth Patience Factor
2359 GPFacLivMod : Modigliani Mortality-Adjusted Growth Patience Factor
2360 GPFacSdl : Szeidl Growth Patience Factor
2361 FHWFac : Finite Human Wealth Factor
2362 RPFac : Return Patience Factor
2363 WRPFac : Weak Return Patience Factor
2364 PFVAFac : Perfect Foresight Value of Autarky Factor
2365 VAFac : Value of Autarky Factor
2366 cNrmPDV : Present Discounted Value of Autarky Consumption
2367 MPCmin : Limiting minimum MPC as market resources go to infinity
2368 MPCmax : Limiting maximum MPC as market resources approach minimum level
2369 hNrm : Human wealth divided by permanent income.
2370 ELogPermShk : Expected log permanent income shock
2371 WorstPrb : Probability of worst income shock realization
2372 Delta_mNrm_ZeroFunc : Linear locus where expected change in market resource ratio is zero
2373 BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income
2375 Returns
2376 -------
2377 None
2378 """
2379 super().calc_limiting_values()
2380 aux_dict = self.bilt
2382 # Calculate the risk-modified growth impatience factor
2383 PermShkDstn = self.PermShkDstn[0]
2384 inv_func = lambda x: x ** (-1.0)
2385 Ex_PermShkInv = expected(inv_func, PermShkDstn)[0]
2386 GroCompPermShk = Ex_PermShkInv ** (-1.0)
2387 aux_dict["GPFacMod"] = aux_dict["APFac"] / (self.PermGroFac[0] * GroCompPermShk)
2389 # Calculate the mortality-adjusted growth impatience factor (and version
2390 # with Modigiliani bequests)
2391 aux_dict["GPFacLiv"] = aux_dict["GPFacRaw"] * self.LivPrb[0]
2392 aux_dict["GPFacLivMod"] = aux_dict["GPFacLiv"] * self.LivPrb[0]
2394 # Calculate the risk-modified value of autarky factor
2395 if self.CRRA == 1.0:
2396 UtilCompPermShk = np.exp(expected(np.log, PermShkDstn)[0])
2397 else:
2398 CRRAfunc = lambda x: x ** (1.0 - self.CRRA)
2399 UtilCompPermShk = expected(CRRAfunc, PermShkDstn)[0] ** (
2400 1 / (1.0 - self.CRRA)
2401 )
2402 aux_dict["VAFac"] = self.DiscFac * (self.PermGroFac[0] * UtilCompPermShk) ** (
2403 1.0 - self.CRRA
2404 )
2406 # Calculate the expected log permanent income shock, which will be used
2407 # for the Szeidl variation of the Growth Impatience condition
2408 aux_dict["ELogPermShk"] = expected(np.log, PermShkDstn)[0]
2410 # Calculate the Harmenberg permanent income neutral expected log permanent
2411 # shock and the Harmenberg Growth Patience Factor
2412 Hrm_func = lambda x: x * np.log(x)
2413 PermShk_Hrm = np.exp(expected(Hrm_func, PermShkDstn)[0])
2414 aux_dict["GPFacHrm"] = aux_dict["GPFacRaw"] / PermShk_Hrm
2416 # Calculate the probability of the worst income shock realization
2417 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2418 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2419 ShkPrbsNext = self.IncShkDstn[0].pmv
2420 Ex_IncNext = np.dot(ShkPrbsNext, PermShkValsNext * TranShkValsNext)
2421 PermShkMinNext = np.min(PermShkValsNext)
2422 TranShkMinNext = np.min(TranShkValsNext)
2423 WorstIncNext = PermShkMinNext * TranShkMinNext
2424 WorstIncPrb = np.sum(
2425 ShkPrbsNext[(PermShkValsNext * TranShkValsNext) == WorstIncNext]
2426 )
2427 aux_dict["WorstPrb"] = WorstIncPrb
2429 # Calculate the weak return patience factor
2430 aux_dict["WRPFac"] = WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2432 # Calculate human wealth and the infinite horizon natural borrowing constraint
2433 if aux_dict["FHWFac"] < 1.0:
2434 hNrm = Ex_IncNext / (1.0 - aux_dict["FHWFac"])
2435 else:
2436 hNrm = np.inf
2437 temp = PermShkMinNext * aux_dict["FHWFac"]
2438 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2440 # Find the upper bound of the MPC as market resources approach the minimum
2441 BoroCnstArt = -np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2442 if BoroCnstNat < BoroCnstArt:
2443 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2444 else:
2445 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2446 MPCmax = np.maximum(MPCmax, 0.0)
2448 # Store maximum MPC and human wealth
2449 aux_dict["hNrm"] = hNrm
2450 aux_dict["MPCmax"] = MPCmax
2452 # Generate the "Delta m = 0" function, which is used to find target market resources
2453 # This overwrites the function generated by the perfect foresight version
2454 Ex_Rnrm = self.Rfree[0] / self.PermGroFac[0] * Ex_PermShkInv
2455 aux_dict["Delta_mNrm_ZeroFunc"] = (
2456 lambda m: (1.0 - 1.0 / Ex_Rnrm) * m + 1.0 / Ex_Rnrm
2457 )
2459 self.bilt = aux_dict
2461 self.bilt = aux_dict
2463 def check_GICMod(self, verbose=None):
2464 """
2465 Evaluate and report on the Risk-Modified Growth Impatience Condition.
2466 """
2467 name = "GICMod"
2468 GPFacMod = self.bilt["GPFacMod"]
2469 result = GPFacMod < 1.0
2471 messages = {
2472 True: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor satisfies the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2473 False: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor violates the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2474 }
2475 verbose = self.verbose if verbose is None else verbose
2476 self.log_condition_result(name, result, messages[result], verbose)
2478 def check_GICSdl(self, verbose=None):
2479 """
2480 Evaluate and report on the Szeidl variation of the Growth Impatience Condition.
2481 """
2482 name = "GICSdl"
2483 ELogPermShk = self.bilt["ELogPermShk"]
2484 result = np.log(self.bilt["GPFacRaw"]) < ELogPermShk
2486 messages = {
2487 True: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock satisfies the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2488 False: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock violates the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2489 }
2490 verbose = self.verbose if verbose is None else verbose
2491 self.log_condition_result(name, result, messages[result], verbose)
2493 def check_GICHrm(self, verbose=None):
2494 """
2495 Evaluate and report on the Harmenberg variation of the Growth Impatience Condition.
2496 """
2497 name = "GICHrm"
2498 GPFacHrm = self.bilt["GPFacHrm"]
2499 result = GPFacHrm < 1.0
2501 messages = {
2502 True: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor satisfies the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2503 False: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor violates the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2504 }
2505 verbose = self.verbose if verbose is None else verbose
2506 self.log_condition_result(name, result, messages[result], verbose)
2508 def check_GICLiv(self, verbose=None):
2509 """
2510 Evaluate and report on the Mortality-Adjusted Growth Impatience Condition.
2511 """
2512 name = "GICLiv"
2513 GPFacLiv = self.bilt["GPFacLiv"]
2514 result = GPFacLiv < 1.0
2516 messages = {
2517 True: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor satisfies the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2518 False: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor violates the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2519 }
2520 verbose = self.verbose if verbose is None else verbose
2521 self.log_condition_result(name, result, messages[result], verbose)
2523 def check_FVAC(self, verbose=None):
2524 """
2525 Evaluate and report on the Finite Value of Autarky condition in the presence of income risk.
2526 """
2527 name = "FVAC"
2528 VAFac = self.bilt["VAFac"]
2529 result = VAFac < 1.0
2531 messages = {
2532 True: f"VAFac={VAFac:.5f} : The Risk-Modified Finite Value of Autarky Factor satisfies the Risk-Modified Finite Value of Autarky Condition β(G‖Ψ‖_(1-ρ))^(1-ρ) < 1.",
2533 False: f"VAFac={VAFac:.5f} : The Risk-Modified Finite Value of Autarky Factor violates the Risk-Modified Finite Value of Autarky Condition β(G‖Ψ‖_(1-ρ))^(1-ρ) < 1.",
2534 }
2535 verbose = self.verbose if verbose is None else verbose
2536 self.log_condition_result(name, result, messages[result], verbose)
2538 def check_WRIC(self, verbose=None):
2539 """
2540 Evaluate and report on the Weak Return Impatience Condition.
2541 """
2542 name = "WRIC"
2543 WRPFac = self.bilt["WRPFac"]
2544 result = WRPFac < 1.0
2546 messages = {
2547 True: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor satisfies the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2548 False: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor violates the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2549 }
2550 verbose = self.verbose if verbose is None else verbose
2551 self.log_condition_result(name, result, messages[result], verbose)
2553 def check_conditions(self, verbose=None):
2554 """
2555 This method checks whether the instance's type satisfies various conditions.
2556 When combinations of these conditions are satisfied, the solution to the
2557 problem exhibits different characteristics. (For an exposition of the
2558 conditions, see https://econ-ark.github.io/BufferStockTheory/)
2560 Parameters
2561 ----------
2562 verbose : boolean
2563 Specifies different levels of verbosity of feedback. When False, it only reports whether the
2564 instance's type fails to satisfy a particular condition. When True, it reports all results, i.e.
2565 the factor values for all conditions.
2567 Returns
2568 -------
2569 None
2570 """
2571 self.conditions = {}
2572 self.bilt["conditions_report"] = ""
2573 self.degenerate = False
2574 verbose = self.verbose if verbose is None else verbose
2576 # This method only checks for the conditions for infinite horizon models
2577 # with a 1 period cycle. If these conditions are not met, we exit early.
2578 if self.cycles != 0 or self.T_cycle > 1:
2579 trivial_message = "No conditions report was produced because this functionality is only supported for infinite horizon models with a cycle length of 1."
2580 self.log_condition_result(None, None, trivial_message, verbose)
2581 if not self.quiet:
2582 _log.info(self.bilt["conditions_report"])
2583 return
2585 # Calculate some useful quantities that will be used in the condition checks
2586 self.calc_limiting_values()
2587 param_desc = self.describe_parameters()
2588 self.log_condition_result(None, None, param_desc, verbose)
2590 # Check individual conditions and add their results to the report
2591 self.check_AIC(verbose)
2592 self.check_RIC(verbose)
2593 self.check_WRIC(verbose)
2594 self.check_GICRaw(verbose)
2595 self.check_GICMod(verbose)
2596 self.check_GICLiv(verbose)
2597 self.check_GICSdl(verbose)
2598 self.check_GICHrm(verbose)
2599 super().check_FVAC(verbose)
2600 self.check_FVAC(verbose)
2601 self.check_FHWC(verbose)
2603 # Exit now if verbose output was not requested.
2604 if not verbose:
2605 if not self.quiet:
2606 _log.info(self.bilt["conditions_report"])
2607 return
2609 # Report on the degeneracy of the consumption function solution
2610 if self.conditions["WRIC"] and self.conditions["FVAC"]:
2611 degen_message = "\nBecause both the WRIC and FVAC are satisfied, the recursive solution to the infinite horizon problem represents a contraction mapping on the consumption function. Thus a non-degenerate solution exists."
2612 degenerate = False
2613 elif not self.conditions["WRIC"]:
2614 degen_message = "\nBecause the WRIC is violated, the consumer is so pathologically patient that they will never consume at all. Thus the solution will be degenerate at c(m) = 0 for all m.\n"
2615 degenerate = True
2616 elif not self.conditions["FVAC"]:
2617 degen_message = "\nBecause the FVAC is violated, the recursive solution to the infinite horizon problem might not be a contraction mapping, so the produced solution might not be valid. Proceed with caution."
2618 degenerate = False
2619 self.log_condition_result(None, None, degen_message, verbose)
2620 self.degenerate = degenerate
2622 # Stop here if the solution is degenerate
2623 if degenerate:
2624 if not self.quiet:
2625 _log.info(self.bilt["conditions_report"])
2626 return
2628 # Report on the limiting behavior of the consumption function as m goes to infinity
2629 if self.conditions["RIC"]:
2630 if self.conditions["FHWC"]:
2631 RIC_message = "\nBecause both the RIC and FHWC condition are satisfied, the consumption function will approach the linear perfect foresight solution as m becomes arbitrarily large."
2632 else:
2633 RIC_message = "\nBecause the RIC is satisfied but the FHWC is violated, the GIC is satisfied."
2634 else:
2635 RIC_message = "\nBecause the RIC is violated, the FHWC condition is also violated. The consumer is pathologically impatient but has infinite expected future earnings. Thus the consumption function will not approach any linear limit as m becomes arbitrarily large, and the MPC will asymptote to zero."
2636 self.log_condition_result(None, None, RIC_message, verbose)
2638 # Report on whether a pseudo-steady-state exists at the individual level
2639 if self.conditions["GICRaw"]:
2640 GIC_message = "\nBecause the GICRaw is satisfied, there exists a pseudo-steady-state wealth ratio at which the level of wealth is expected to grow at the same rate as permanent income."
2641 else:
2642 GIC_message = "\nBecause the GICRaw is violated, there might not exist a pseudo-steady-state wealth ratio at which the level of wealth is expected to grow at the same rate as permanent income."
2643 self.log_condition_result(None, None, GIC_message, verbose)
2645 # Report on whether a target wealth ratio exists at the individual level
2646 if self.conditions["GICMod"]:
2647 GICMod_message = "\nBecause the GICMod is satisfied, expected growth of the ratio of market resources to permanent income is less than one as market resources become arbitrarily large. Hence the consumer has a target ratio of market resources to permanent income."
2648 else:
2649 GICMod_message = "\nBecause the GICMod is violated, expected growth of the ratio of market resources to permanent income exceeds one as market resources go to infinity. Hence the consumer might not have a target ratio of market resources to permanent income."
2650 self.log_condition_result(None, None, GICMod_message, verbose)
2652 # Report on whether a target level of wealth exists at the aggregate level
2653 if self.conditions["GICLiv"]:
2654 GICLiv_message = "\nBecause the GICLiv is satisfied, a target ratio of aggregate market resources to aggregate permanent income exists."
2655 else:
2656 GICLiv_message = "\nBecause the GICLiv is violated, a target ratio of aggregate market resources to aggregate permanent income might not exist."
2657 self.log_condition_result(None, None, GICLiv_message, verbose)
2659 # Report on whether invariant distributions exist
2660 if self.conditions["GICSdl"]:
2661 GICSdl_message = "\nBecause the GICSdl is satisfied, there exist invariant distributions of permanent income-normalized variables."
2662 else:
2663 GICSdl_message = "\nBecause the GICSdl is violated, there do not exist invariant distributions of permanent income-normalized variables."
2664 self.log_condition_result(None, None, GICSdl_message, verbose)
2666 # Report on whether blah blah
2667 if self.conditions["GICHrm"]:
2668 GICHrm_message = "\nBecause the GICHrm is satisfied, there exists a target ratio of the individual market resources to permanent income, under the permanent-income-neutral measure."
2669 else:
2670 GICHrm_message = "\nBecause the GICHrm is violated, there does not exist a target ratio of the individual market resources to permanent income, under the permanent-income-neutral measure.."
2671 self.log_condition_result(None, None, GICHrm_message, verbose)
2673 if not self.quiet:
2674 _log.info(self.bilt["conditions_report"])
2677###############################################################################
2679# Specify default parameters used in "kinked R" model
2681KinkedRconsumerType_IncShkDstn_default = IndShockConsumerType_IncShkDstn_default.copy()
2682KinkedRconsumerType_aXtraGrid_default = IndShockConsumerType_aXtraGrid_default.copy()
2683KinkedRconsumerType_kNrmInitDstn_default = (
2684 IndShockConsumerType_kNrmInitDstn_default.copy()
2685)
2686KinkedRconsumerType_pLvlInitDstn_default = (
2687 IndShockConsumerType_pLvlInitDstn_default.copy()
2688)
2690KinkedRconsumerType_solving_default = IndShockConsumerType_solving_default.copy()
2691KinkedRconsumerType_solving_default.update(
2692 {
2693 "Rboro": 1.20, # Interest factor on assets when borrowing, a < 0
2694 "Rsave": 1.02, # Interest factor on assets when saving, a > 0
2695 "BoroCnstArt": None, # Kinked R only matters if borrowing is allowed
2696 }
2697)
2698del KinkedRconsumerType_solving_default["Rfree"]
2700KinkedRconsumerType_simulation_default = IndShockConsumerType_simulation_default.copy()
2702KinkedRconsumerType_defaults = {}
2703KinkedRconsumerType_defaults.update(
2704 KinkedRconsumerType_IncShkDstn_default
2705) # Fill with some parameters
2706KinkedRconsumerType_defaults.update(KinkedRconsumerType_pLvlInitDstn_default)
2707KinkedRconsumerType_defaults.update(KinkedRconsumerType_kNrmInitDstn_default)
2708KinkedRconsumerType_defaults.update(KinkedRconsumerType_aXtraGrid_default)
2709KinkedRconsumerType_defaults.update(KinkedRconsumerType_solving_default)
2710KinkedRconsumerType_defaults.update(KinkedRconsumerType_simulation_default)
2711init_kinked_R = KinkedRconsumerType_defaults
2714class KinkedRconsumerType(IndShockConsumerType):
2715 r"""
2716 A consumer type based on IndShockConsumerType, with different
2717 interest rates for saving (:math:`\mathsf{R}_{save}`) and borrowing
2718 (:math:`\mathsf{R}_{boro}`).
2720 Solver for this class is currently only compatible with linear spline interpolation.
2722 .. math::
2723 \newcommand{\CRRA}{\rho}
2724 \newcommand{\DiePrb}{\mathsf{D}}
2725 \newcommand{\PermGroFac}{\Gamma}
2726 \newcommand{\Rfree}{\mathsf{R}}
2727 \newcommand{\DiscFac}{\beta}
2728 \begin{align*}
2729 v_t(m_t) &= \max_{c_t} u(c_t) + \DiscFac (1-\DiePrb_{t+1}) \mathbb{E}_{t} \left[(\PermGroFac_{t+1}\psi_{t+1})^{1-\CRRA} v_{t+1}(m_{t+1}) \right], \\
2730 & \text{s.t.} \\
2731 a_t &= m_t - c_t, \\
2732 a_t &\geq \underline{a}, \\
2733 m_{t+1} &= \Rfree_t/(\PermGroFac_{t+1} \psi_{t+1}) a_t + \theta_{t+1}, \\
2734 \Rfree_t &= \begin{cases}
2735 \Rfree_{boro} & \text{if } a_t < 0\\
2736 \Rfree_{save} & \text{if } a_t \geq 0,
2737 \end{cases}\\
2738 \Rfree_{boro} &> \Rfree_{save}, \\
2739 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
2740 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1.\\
2741 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA} \\
2742 \end{align*}
2745 Constructors
2746 ------------
2747 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2748 The agent's income shock distributions.
2750 It's default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2751 aXtraGrid: Constructor
2752 The agent's asset grid.
2754 It's default constructor is :func:`HARK.utilities.make_assets_grid`
2756 Solving Parameters
2757 ------------------
2758 cycles: int
2759 0 specifies an infinite horizon model, 1 specifies a finite model.
2760 T_cycle: int
2761 Number of periods in the cycle for this agent type.
2762 CRRA: float, :math:`\rho`
2763 Coefficient of Relative Risk Aversion.
2764 Rboro: float, :math:`\mathsf{R}_{boro}`
2765 Risk Free interest rate when assets are negative.
2766 Rsave: float, :math:`\mathsf{R}_{save}`
2767 Risk Free interest rate when assets are positive.
2768 DiscFac: float, :math:`\beta`
2769 Intertemporal discount factor.
2770 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2771 Survival probability after each period.
2772 PermGroFac: list[float], time varying, :math:`\Gamma`
2773 Permanent income growth factor.
2774 BoroCnstArt: float, :math:`\underline{a}`
2775 The minimum Asset/Perminant Income ratio, None to ignore.
2776 vFuncBool: bool
2777 Whether to calculate the value function during solution.
2778 CubicBool: bool
2779 Whether to use cubic spline interpoliation.
2781 Simulation Parameters
2782 ---------------------
2783 AgentCount: int
2784 Number of agents of this kind that are created during simulations.
2785 T_age: int
2786 Age after which to automatically kill agents, None to ignore.
2787 T_sim: int, required for simulation
2788 Number of periods to simulate.
2789 track_vars: list[strings]
2790 List of variables that should be tracked when running the simulation.
2791 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2793 PermShk is the agent's permanent income shock
2795 TranShk is the agent's transitory income shock
2797 aLvl is the nominal asset level
2799 aNrm is the normalized assets
2801 bNrm is the normalized resources without this period's labor income
2803 cNrm is the normalized consumption
2805 mNrm is the normalized market resources
2807 pLvl is the permanent income level
2809 who_dies is the array of which agents died
2810 aNrmInitMean: float
2811 Mean of Log initial Normalized Assets.
2812 aNrmInitStd: float
2813 Std of Log initial Normalized Assets.
2814 pLvlInitMean: float
2815 Mean of Log initial permanent income.
2816 pLvlInitStd: float
2817 Std of Log initial permanent income.
2818 PermGroFacAgg: float
2819 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2820 PerfMITShk: boolean
2821 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2822 NewbornTransShk: boolean
2823 Whether Newborns have transitory shock.
2825 Attributes
2826 ----------
2827 solution: list[Consumer solution object]
2828 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2829 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2831 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2832 history: Dict[Array]
2833 Created by running the :func:`.simulate()` method.
2834 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2835 Visit :class:`HARK.core.AgentType.simulate` for more information.
2836 """
2838 IncShkDstn_defaults = KinkedRconsumerType_IncShkDstn_default
2839 aXtraGrid_defaults = KinkedRconsumerType_aXtraGrid_default
2840 solving_defaults = KinkedRconsumerType_solving_default
2841 simulation_defaults = KinkedRconsumerType_simulation_default
2842 default_ = {
2843 "params": KinkedRconsumerType_defaults,
2844 "solver": solve_one_period_ConsKinkedR,
2845 "model": "ConsKinkedR.yaml",
2846 }
2848 time_inv_ = copy(IndShockConsumerType.time_inv_)
2849 time_inv_ += ["Rboro", "Rsave"]
2851 def calc_bounding_values(self):
2852 """
2853 Calculate human wealth plus minimum and maximum MPC in an infinite
2854 horizon model with only one period repeated indefinitely. Store results
2855 as attributes of self. Human wealth is the present discounted value of
2856 expected future income after receiving income this period, ignoring mort-
2857 ality. The maximum MPC is the limit of the MPC as m --> mNrmMin. The
2858 minimum MPC is the limit of the MPC as m --> infty. This version deals
2859 with the different interest rates on borrowing vs saving.
2861 Parameters
2862 ----------
2863 None
2865 Returns
2866 -------
2867 None
2868 """
2869 # Unpack the income distribution and get average and worst outcomes
2870 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2871 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2872 ShkPrbsNext = self.IncShkDstn[0].pmv
2873 IncNext = PermShkValsNext * TranShkValsNext
2874 Ex_IncNext = np.dot(ShkPrbsNext, IncNext)
2875 PermShkMinNext = np.min(PermShkValsNext)
2876 TranShkMinNext = np.min(TranShkValsNext)
2877 WorstIncNext = PermShkMinNext * TranShkMinNext
2878 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
2879 # TODO: Check the math above. I think it fails for non-independent shocks
2881 BoroCnstArt = np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2883 # Calculate human wealth and the infinite horizon natural borrowing constraint
2884 hNrm = (Ex_IncNext * self.PermGroFac[0] / self.Rsave) / (
2885 1.0 - self.PermGroFac[0] / self.Rsave
2886 )
2887 temp = self.PermGroFac[0] * PermShkMinNext / self.Rboro
2888 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2890 PatFacTop = (self.DiscFac * self.LivPrb[0] * self.Rsave) ** (
2891 1.0 / self.CRRA
2892 ) / self.Rsave
2893 PatFacBot = (self.DiscFac * self.LivPrb[0] * self.Rboro) ** (
2894 1.0 / self.CRRA
2895 ) / self.Rboro
2896 if BoroCnstNat < BoroCnstArt:
2897 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2898 else:
2899 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * PatFacBot
2900 MPCmin = 1.0 - PatFacTop
2902 # Store the results as attributes of self
2903 self.hNrm = hNrm
2904 self.MPCmin = MPCmin
2905 self.MPCmax = MPCmax
2907 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True): # pragma: nocover
2908 """
2909 Creates a "normalized Euler error" function for this instance, mapping
2910 from market resources to "consumption error per dollar of consumption."
2911 Stores result in attribute eulerErrorFunc as an interpolated function.
2912 Has option to use approximate income distribution stored in self.IncShkDstn
2913 or to use a (temporary) very dense approximation.
2915 SHOULD BE INHERITED FROM ConsIndShockModel
2917 Parameters
2918 ----------
2919 mMax : float
2920 Maximum normalized market resources for the Euler error function.
2921 approx_inc_dstn : Boolean
2922 Indicator for whether to use the approximate discrete income distri-
2923 bution stored in self.IncShkDstn[0], or to use a very accurate
2924 discrete approximation instead. When True, uses approximation in
2925 IncShkDstn; when False, makes and uses a very dense approximation.
2927 Returns
2928 -------
2929 None
2930 """
2931 raise NotImplementedError()
2933 def get_Rfree(self):
2934 """
2935 Returns an array of size self.AgentCount with self.Rboro or self.Rsave in each entry, based
2936 on whether self.aNrmNow >< 0.
2938 Parameters
2939 ----------
2940 None
2942 Returns
2943 -------
2944 RfreeNow : np.array
2945 Array of size self.AgentCount with risk free interest rate for each agent.
2946 """
2947 RfreeNow = self.Rboro * np.ones(self.AgentCount)
2948 RfreeNow[self.state_prev["aNrm"] > 0] = self.Rsave
2949 return RfreeNow
2951 def check_conditions(self, verbose):
2952 """
2953 This empty method overwrites the version inherited from its parent class,
2954 IndShockConsumerType. The condition checks are not appropriate when Rfree
2955 has multiple values.
2957 Parameters
2958 ----------
2959 None
2961 Returns
2962 -------
2963 None
2964 """
2965 pass
2968###############################################################################
2970# Make a dictionary to specify a lifecycle consumer with a finite horizon
2972# Main calibration characteristics
2973birth_age = 25
2974death_age = 90
2975adjust_infl_to = 1992
2976# Use income estimates from Cagetti (2003) for High-school graduates
2977education = "HS"
2978income_calib = Cagetti_income[education]
2980# Income specification
2981income_params = parse_income_spec(
2982 age_min=birth_age,
2983 age_max=death_age,
2984 adjust_infl_to=adjust_infl_to,
2985 **income_calib,
2986 SabelhausSong=True,
2987)
2989# Initial distribution of wealth and permanent income
2990dist_params = income_wealth_dists_from_scf(
2991 base_year=adjust_infl_to, age=birth_age, education=education, wave=1995
2992)
2994# We need survival probabilities only up to death_age-1, because survival
2995# probability at death_age is 1.
2996liv_prb = parse_ssa_life_table(
2997 female=False, cross_sec=True, year=2004, age_min=birth_age, age_max=death_age
2998)
3000# Parameters related to the number of periods implied by the calibration
3001time_params = parse_time_params(age_birth=birth_age, age_death=death_age)
3003# Update all the new parameters
3004init_lifecycle = copy(init_idiosyncratic_shocks)
3005del init_lifecycle["constructors"]
3006init_lifecycle.update(time_params)
3007init_lifecycle.update(dist_params)
3008# Note the income specification overrides the pLvlInitMean from the SCF.
3009init_lifecycle.update(income_params)
3010init_lifecycle.update({"LivPrb": liv_prb})
3011init_lifecycle["Rfree"] = init_lifecycle["T_cycle"] * init_lifecycle["Rfree"]
3013# Make a dictionary to specify an infinite consumer with a four period cycle
3014init_cyclical = copy(init_idiosyncratic_shocks)
3015init_cyclical["PermGroFac"] = [1.1, 1.082251, 2.8, 0.3]
3016init_cyclical["PermShkStd"] = [0.1, 0.1, 0.1, 0.1]
3017init_cyclical["TranShkStd"] = [0.1, 0.1, 0.1, 0.1]
3018init_cyclical["LivPrb"] = 4 * [0.98]
3019init_cyclical["Rfree"] = 4 * [1.03]
3020init_cyclical["T_cycle"] = 4