Coverage for HARK / ConsumptionSaving / ConsIndShockModel.py: 99%
882 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +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 interpolation.
871 Returns
872 -------
873 solution_now : ConsumerSolution
874 Solution to this period's consumption-saving problem with income risk.
876 """
877 # Verifiy that there is actually a kink in the interest factor
878 assert Rboro >= Rsave, (
879 "Interest factor on debt less than interest factor on savings!"
880 )
881 # If the kink is in the wrong direction, code should break here. If there's
882 # no kink at all, then just use the ConsIndShockModel solver.
883 if Rboro == Rsave:
884 solution_now = solve_one_period_ConsIndShock(
885 solution_next,
886 IncShkDstn,
887 LivPrb,
888 DiscFac,
889 CRRA,
890 Rboro,
891 PermGroFac,
892 BoroCnstArt,
893 aXtraGrid,
894 vFuncBool,
895 CubicBool,
896 )
897 return solution_now
899 # Define the current period utility function and effective discount factor
900 uFunc = UtilityFuncCRRA(CRRA)
901 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
903 # Calculate the probability that we get the worst possible income draw
904 WorstIncPrb = calc_worst_inc_prob(IncShkDstn, use_infimum=False)
905 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
906 Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn)
907 hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext)
909 # Unpack next period's (marginal) value function
910 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
911 vPfuncNext = solution_next.vPfunc
912 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
914 # Calculate the minimum allowable value of money resources in this period
915 BoroCnstNat = calc_boro_const_nat(
916 solution_next.mNrmMin,
917 IncShkDstn,
918 Rboro,
919 PermGroFac,
920 use_infimum=False,
921 )
922 # Set the minimum allowable (normalized) market resources based on the natural
923 # and artificial borrowing constraints
924 mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat)
926 # Update the bounding MPCs and PDV of human wealth:
927 PatFacSave = calc_patience_factor(Rsave, DiscFacEff, CRRA)
928 PatFacBoro = calc_patience_factor(Rboro, DiscFacEff, CRRA)
929 MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave)
930 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
931 # or artificial borrowing constraint actually binds
932 MPCmaxUnc = calc_mpc_max(
933 solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro, BoroCnstNat, BoroCnstArt
934 )
935 MPCmaxNow = 1.0 if BoroCnstNat < mNrmMinNow else MPCmaxUnc
937 cFuncLimitIntercept = MPCminNow * hNrmNow
938 cFuncLimitSlope = MPCminNow
940 # Define the borrowing-constrained consumption function
941 cFuncNowCnst = LinearInterp(
942 np.array([mNrmMinNow, mNrmMinNow + 1.0]),
943 np.array([0.0, 1.0]),
944 )
946 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
947 aNrmNow = np.sort(
948 np.hstack((np.asarray(aXtraGrid) + mNrmMinNow, np.array([0.0, 1e-15]))),
949 )
951 # Make a 1D array of the interest factor at each asset gridpoint
952 Rfree = Rsave * np.ones_like(aNrmNow)
953 Rfree[aNrmNow <= 0] = Rboro
954 i_kink = np.argwhere(aNrmNow == 0.0)[0][0]
956 # Calculate end-of-period marginal value of assets at each gridpoint
957 vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA)
958 EndOfPrdvP = vPfacEff * expected(
959 calc_vp_next,
960 IncShkDstn,
961 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPfuncNext),
962 )
964 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
965 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
966 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
968 # Limiting consumption is zero as m approaches mNrmMin
969 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
970 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
972 # Construct the consumption function as a cubic or linear spline interpolation
973 if CubicBool:
974 # Calculate end-of-period marginal marginal value of assets at each gridpoint
975 vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0)
976 EndOfPrdvPP = vPPfacEff * expected(
977 calc_vpp_next,
978 IncShkDstn,
979 args=(aNrmNow, Rfree, CRRA, PermGroFac, vPPfuncNext),
980 )
981 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
982 MPC = dcda / (dcda + 1.0)
983 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxUnc)
985 # Construct the unconstrained consumption function as a cubic interpolation
986 cFuncNowUnc = CubicInterp(
987 m_for_interpolation,
988 c_for_interpolation,
989 MPC_for_interpolation,
990 cFuncLimitIntercept,
991 cFuncLimitSlope,
992 )
993 # Adjust the coefficients on the kinked portion of the cFunc
994 cFuncNowUnc.coeffs[i_kink + 2] = [
995 c_for_interpolation[i_kink + 1],
996 m_for_interpolation[i_kink + 2] - m_for_interpolation[i_kink + 1],
997 0.0,
998 0.0,
999 ]
1000 else:
1001 # Construct the unconstrained consumption function as a linear interpolation
1002 cFuncNowUnc = LinearInterp(
1003 m_for_interpolation,
1004 c_for_interpolation,
1005 cFuncLimitIntercept,
1006 cFuncLimitSlope,
1007 )
1009 # Combine the constrained and unconstrained functions into the true consumption function.
1010 # LowerEnvelope should only be used when BoroCnstArt is True
1011 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1013 # Make the marginal value function and the marginal marginal value function
1014 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1016 # Define this period's marginal marginal value function
1017 if CubicBool:
1018 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1019 else:
1020 vPPfuncNow = NullFunc() # Dummy object
1022 # Construct this period's value function if requested
1023 if vFuncBool:
1024 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1025 EndOfPrdv = DiscFacEff * expected(
1026 calc_v_next,
1027 IncShkDstn,
1028 args=(aNrmNow, Rfree, CRRA, PermGroFac, vFuncNext),
1029 )
1030 EndOfPrdvNvrs = uFunc.inv(
1031 EndOfPrdv,
1032 ) # value transformed through inverse utility
1033 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1034 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1035 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1036 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1038 # Construct the end-of-period value function
1039 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1040 EndOfPrdvNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1041 EndOfPrdvFunc = ValueFuncCRRA(EndOfPrdvNvrsFunc, CRRA)
1043 # Compute expected value and marginal value on a grid of market resources
1044 mNrm_temp = mNrmMinNow + aXtraGrid
1045 cNrm_temp = cFuncNow(mNrm_temp)
1046 aNrm_temp = mNrm_temp - cNrm_temp
1047 v_temp = uFunc(cNrm_temp) + EndOfPrdvFunc(aNrm_temp)
1048 vP_temp = uFunc.der(cNrm_temp)
1050 # Construct the beginning-of-period value function
1051 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1052 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1053 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1054 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1055 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
1056 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1057 vNvrsFuncNow = CubicInterp(
1058 mNrm_temp,
1059 vNvrs_temp,
1060 vNvrsP_temp,
1061 MPCminNvrs * hNrmNow,
1062 MPCminNvrs,
1063 )
1064 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1065 else:
1066 vFuncNow = NullFunc() # Dummy object
1068 # Create and return this period's solution
1069 solution_now = ConsumerSolution(
1070 cFunc=cFuncNow,
1071 vFunc=vFuncNow,
1072 vPfunc=vPfuncNow,
1073 vPPfunc=vPPfuncNow,
1074 mNrmMin=mNrmMinNow,
1075 hNrm=hNrmNow,
1076 MPCmin=MPCminNow,
1077 MPCmax=MPCmaxNow,
1078 )
1079 return solution_now
1082def make_basic_CRRA_solution_terminal(CRRA):
1083 """
1084 Construct the terminal period solution for a consumption-saving model with
1085 CRRA utility and only one state variable.
1087 Parameters
1088 ----------
1089 CRRA : float
1090 Coefficient of relative risk aversion. This is the only relevant parameter.
1092 Returns
1093 -------
1094 solution_terminal : ConsumerSolution
1095 Terminal period solution for someone with the given CRRA.
1096 """
1097 cFunc_terminal = LinearInterp([0.0, 1.0], [0.0, 1.0]) # c=m at t=T
1098 vFunc_terminal = ValueFuncCRRA(cFunc_terminal, CRRA)
1099 vPfunc_terminal = MargValueFuncCRRA(cFunc_terminal, CRRA)
1100 vPPfunc_terminal = MargMargValueFuncCRRA(cFunc_terminal, CRRA)
1101 solution_terminal = ConsumerSolution(
1102 cFunc=cFunc_terminal,
1103 vFunc=vFunc_terminal,
1104 vPfunc=vPfunc_terminal,
1105 vPPfunc=vPPfunc_terminal,
1106 mNrmMin=0.0,
1107 hNrm=0.0,
1108 MPCmin=1.0,
1109 MPCmax=1.0,
1110 )
1111 return solution_terminal
1114# ============================================================================
1115# == Classes for representing types of consumer agents (and things they do) ==
1116# ============================================================================
1118# Make a dictionary of constructors (very simply for perfect foresight model)
1119PerfForesightConsumerType_constructors_default = {
1120 "solution_terminal": make_basic_CRRA_solution_terminal,
1121 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1122 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1123}
1125# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1126PerfForesightConsumerType_kNrmInitDstn_default = {
1127 "kLogInitMean": -12.0, # Mean of log initial capital
1128 "kLogInitStd": 0.0, # Stdev of log initial capital
1129 "kNrmInitCount": 15, # Number of points in initial capital discretization
1130}
1132# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1133PerfForesightConsumerType_pLvlInitDstn_default = {
1134 "pLogInitMean": 0.0, # Mean of log permanent income
1135 "pLogInitStd": 0.0, # Stdev of log permanent income
1136 "pLvlInitCount": 15, # Number of points in initial capital discretization
1137}
1139# Make a dictionary to specify a perfect foresight consumer type
1140PerfForesightConsumerType_solving_defaults = {
1141 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1142 "cycles": 1, # Finite, non-cyclic model
1143 "T_cycle": 1, # Number of periods in the cycle for this agent type
1144 "pseudo_terminal": False, # Terminal period really does exist
1145 "constructors": PerfForesightConsumerType_constructors_default, # See dictionary above
1146 # PARAMETERS REQUIRED TO SOLVE THE MODEL
1147 "CRRA": 2.0, # Coefficient of relative risk aversion
1148 "Rfree": [1.03], # Interest factor on retained assets
1149 "DiscFac": 0.96, # Intertemporal discount factor
1150 "LivPrb": [0.98], # Survival probability after each period
1151 "PermGroFac": [1.01], # Permanent income growth factor
1152 "BoroCnstArt": None, # Artificial borrowing constraint
1153 "MaxKinks": 400, # Maximum number of grid points to allow in cFunc
1154}
1155PerfForesightConsumerType_simulation_defaults = {
1156 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1157 "AgentCount": 10000, # Number of agents of this type
1158 "T_age": None, # Age after which simulated agents are automatically killed
1159 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1160 # (The portion of PermGroFac attributable to aggregate productivity growth)
1161 # ADDITIONAL OPTIONAL PARAMETERS
1162 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1163 # (Forces Newborns to follow solution path of the agent they replaced if True)
1164}
1165PerfForesightConsumerType_defaults = {}
1166PerfForesightConsumerType_defaults.update(PerfForesightConsumerType_solving_defaults)
1167PerfForesightConsumerType_defaults.update(
1168 PerfForesightConsumerType_kNrmInitDstn_default
1169)
1170PerfForesightConsumerType_defaults.update(
1171 PerfForesightConsumerType_pLvlInitDstn_default
1172)
1173PerfForesightConsumerType_defaults.update(PerfForesightConsumerType_simulation_defaults)
1174init_perfect_foresight = PerfForesightConsumerType_defaults
1177class PerfForesightConsumerType(AgentType):
1178 r"""
1179 A perfect foresight consumer type who has no uncertainty other than mortality.
1180 Their problem is defined by a coefficient of relative risk aversion (:math:`\rho`), intertemporal
1181 discount factor (:math:`\beta`), interest factor (:math:`\mathsf{R}`), an optional artificial borrowing constraint (:math:`\underline{a}`)
1182 and time sequences of the permanent income growth rate (:math:`\Gamma`) and survival probability (:math:`1-\mathsf{D}`).
1183 Their assets and income are normalized by permanent income.
1185 .. math::
1186 \newcommand{\CRRA}{\rho}
1187 \newcommand{\DiePrb}{\mathsf{D}}
1188 \newcommand{\PermGroFac}{\Gamma}
1189 \newcommand{\Rfree}{\mathsf{R}}
1190 \newcommand{\DiscFac}{\beta}
1191 \begin{align*}
1192 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}), \\
1193 & \text{s.t.} \\
1194 a_t &= m_t - c_t, \\
1195 a_t &\geq \underline{a}, \\
1196 m_{t+1} &= \Rfree_{t+1} a_t/\PermGroFac_{t+1} + 1, \\
1197 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}
1198 \end{align*}
1201 Solving Parameters
1202 ------------------
1203 cycles: int
1204 0 specifies an infinite horizon model, 1 specifies a finite model.
1205 T_cycle: int
1206 Number of periods in the cycle for this agent type.
1207 CRRA: float, :math:`\rho`
1208 Coefficient of Relative Risk Aversion.
1209 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1210 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1211 DiscFac: float, :math:`\beta`
1212 Intertemporal discount factor.
1213 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1214 Survival probability after each period.
1215 PermGroFac: list[float], time varying, :math:`\Gamma`
1216 Permanent income growth factor.
1217 BoroCnstArt: float, :math:`\underline{a}`
1218 The minimum Asset/Perminant Income ratio, None to ignore.
1219 MaxKinks: int
1220 Maximum number of gridpoints to allow in cFunc.
1222 Simulation Parameters
1223 ---------------------
1224 AgentCount: int
1225 Number of agents of this kind that are created during simulations.
1226 T_age: int
1227 Age after which to automatically kill agents, None to ignore.
1228 T_sim: int, required for simulation
1229 Number of periods to simulate.
1230 track_vars: list[strings]
1231 List of variables that should be tracked when running the simulation.
1232 For this agent, the options are 'kNrm', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
1234 kNrm is beginning-of-period capital holdings (last period's assets)
1236 aLvl is the nominal asset level
1238 aNrm is the normalized assets
1240 bNrm is the normalized resources without this period's labor income
1242 cNrm is the normalized consumption
1244 mNrm is the normalized market resources
1246 pLvl is the permanent income level
1248 who_dies is the array of which agents died
1249 aNrmInitMean: float
1250 Mean of Log initial Normalized Assets.
1251 aNrmInitStd: float
1252 Std of Log initial Normalized Assets.
1253 pLvlInitMean: float
1254 Mean of Log initial permanent income.
1255 pLvlInitStd: float
1256 Std of Log initial permanent income.
1257 PermGroFacAgg: float
1258 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1259 PerfMITShk: boolean
1260 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1262 Attributes
1263 ----------
1264 solution: list[Consumer solution object]
1265 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1266 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1268 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1269 history: Dict[Array]
1270 Created by running the :func:`.simulate()` method.
1271 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1272 Visit :class:`HARK.core.AgentType.simulate` for more information.
1273 """
1275 solving_defaults = PerfForesightConsumerType_solving_defaults
1276 simulation_defaults = PerfForesightConsumerType_simulation_defaults
1278 default_ = {
1279 "params": PerfForesightConsumerType_defaults,
1280 "solver": solve_one_period_ConsPF,
1281 "model": "ConsPerfForesight.yaml",
1282 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
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_Rport(self):
1422 """
1423 Returns an array of size self.AgentCount with Rfree in every entry,
1424 representing the risk-free portfolio return
1426 Parameters
1427 ----------
1428 None
1430 Returns
1431 -------
1432 RfreeNow : np.array
1433 Array of size self.AgentCount with risk free interest rate for each agent.
1434 """
1435 Rfree_array = np.array(self.Rfree)
1436 return Rfree_array[self.t_cycle - 1]
1438 def transition(self):
1439 pLvlPrev = self.state_prev["pLvl"]
1440 kNrm = self.state_prev["aNrm"]
1441 RportNow = self.get_Rport()
1443 # Calculate new states: normalized market resources and permanent income level
1444 # Updated permanent income level
1445 pLvlNow = pLvlPrev * self.shocks["PermShk"]
1446 # "Effective" interest factor on normalized assets
1447 ReffNow = RportNow / self.shocks["PermShk"]
1448 bNrmNow = ReffNow * kNrm # Bank balances before labor income
1449 # Market resources after income
1450 mNrmNow = bNrmNow + self.shocks["TranShk"]
1452 return kNrm, pLvlNow, bNrmNow, mNrmNow, None
1454 def get_controls(self):
1455 """
1456 Calculates consumption for each consumer of this type using the consumption functions.
1458 Parameters
1459 ----------
1460 None
1462 Returns
1463 -------
1464 None
1465 """
1466 cNrmNow = np.full(self.AgentCount, np.nan)
1467 MPCnow = np.full(self.AgentCount, np.nan)
1468 for t in np.unique(self.t_cycle):
1469 idx = self.t_cycle == t
1470 if np.any(idx):
1471 cNrmNow[idx], MPCnow[idx] = self.solution[t].cFunc.eval_with_derivative(
1472 self.state_now["mNrm"][idx]
1473 )
1474 self.controls["cNrm"] = cNrmNow
1476 # MPCnow is not really a control
1477 self.MPCnow = MPCnow
1479 def get_poststates(self):
1480 """
1481 Calculates end-of-period assets for each consumer of this type.
1483 Parameters
1484 ----------
1485 None
1487 Returns
1488 -------
1489 None
1490 """
1491 self.state_now["aNrm"] = self.state_now["mNrm"] - self.controls["cNrm"]
1492 self.state_now["aLvl"] = self.state_now["aNrm"] * self.state_now["pLvl"]
1493 # Update aggregate permanent productivity level
1494 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1496 def log_condition_result(self, name, result, message, verbose):
1497 """
1498 Records the result of one condition check in the attribute condition_report
1499 of the bilt dictionary, and in the message log.
1501 Parameters
1502 ----------
1503 name : string or None
1504 Name for the condition; if None, no test result is added to conditions.
1505 result : bool
1506 An indicator for whether the condition was passed.
1507 message : str
1508 The messages to record about the condition check.
1509 verbose : bool
1510 Indicator for whether verbose messages should be included in the report.
1511 """
1512 if name is not None:
1513 self.conditions[name] = result
1514 set_verbosity_level((4 - verbose) * 10)
1515 _log.info(message)
1516 self.bilt["conditions_report"] += message + "\n"
1518 def check_AIC(self, verbose=None):
1519 """
1520 Evaluate and report on the Absolute Impatience Condition.
1521 """
1522 name = "AIC"
1523 APFac = self.bilt["APFac"]
1524 result = APFac < 1.0
1526 messages = {
1527 True: f"APFac={APFac:.5f} : The Absolute Patience Factor satisfies the Absolute Impatience Condition (AIC) Þ < 1.",
1528 False: f"APFac={APFac:.5f} : The Absolute Patience Factor violates the Absolute Impatience Condition (AIC) Þ < 1.",
1529 }
1530 verbose = self.verbose if verbose is None else verbose
1531 self.log_condition_result(name, result, messages[result], verbose)
1533 def check_GICRaw(self, verbose=None):
1534 """
1535 Evaluate and report on the Growth Impatience Condition for the Perfect Foresight model.
1536 """
1537 name = "GICRaw"
1538 GPFacRaw = self.bilt["GPFacRaw"]
1539 result = GPFacRaw < 1.0
1541 messages = {
1542 True: f"GPFacRaw={GPFacRaw:.5f} : The Growth Patience Factor satisfies the Growth Impatience Condition (GICRaw) Þ/G < 1.",
1543 False: f"GPFacRaw={GPFacRaw:.5f} : The Growth Patience Factor violates the Growth Impatience Condition (GICRaw) Þ/G < 1.",
1544 }
1545 verbose = self.verbose if verbose is None else verbose
1546 self.log_condition_result(name, result, messages[result], verbose)
1548 def check_RIC(self, verbose=None):
1549 """
1550 Evaluate and report on the Return Impatience Condition.
1551 """
1552 name = "RIC"
1553 RPFac = self.bilt["RPFac"]
1554 result = RPFac < 1.0
1556 messages = {
1557 True: f"RPFac={RPFac:.5f} : The Return Patience Factor satisfies the Return Impatience Condition (RIC) Þ/R < 1.",
1558 False: f"RPFac={RPFac:.5f} : The Return Patience Factor violates the Return Impatience Condition (RIC) Þ/R < 1.",
1559 }
1560 verbose = self.verbose if verbose is None else verbose
1561 self.log_condition_result(name, result, messages[result], verbose)
1563 def check_FHWC(self, verbose=None):
1564 """
1565 Evaluate and report on the Finite Human Wealth Condition.
1566 """
1567 name = "FHWC"
1568 FHWFac = self.bilt["FHWFac"]
1569 result = FHWFac < 1.0
1571 messages = {
1572 True: f"FHWFac={FHWFac:.5f} : The Finite Human Wealth Factor satisfies the Finite Human Wealth Condition (FHWC) G/R < 1.",
1573 False: f"FHWFac={FHWFac:.5f} : The Finite Human Wealth Factor violates the Finite Human Wealth Condition (FHWC) G/R < 1.",
1574 }
1575 verbose = self.verbose if verbose is None else verbose
1576 self.log_condition_result(name, result, messages[result], verbose)
1578 def check_FVAC(self, verbose=None):
1579 """
1580 Evaluate and report on the Finite Value of Autarky Condition under perfect foresight.
1581 """
1582 name = "PFFVAC"
1583 PFVAFac = self.bilt["PFVAFac"]
1584 result = PFVAFac < 1.0
1586 messages = {
1587 True: f"PFVAFac={PFVAFac:.5f} : The Finite Value of Autarky Factor satisfies the Finite Value of Autarky Condition βG^(1-ρ) < 1.",
1588 False: f"PFVAFac={PFVAFac:.5f} : The Finite Value of Autarky Factor violates the Finite Value of Autarky Condition βG^(1-ρ) < 1.",
1589 }
1590 verbose = self.verbose if verbose is None else verbose
1591 self.log_condition_result(name, result, messages[result], verbose)
1593 def describe_parameters(self):
1594 """
1595 Make a string describing this instance's parameter values, including their
1596 representation in code and symbolically.
1598 Returns
1599 -------
1600 param_desc : str
1601 Description of parameters as a unicode string.
1602 """
1603 params_to_describe = [
1604 # [name, description, symbol, time varying]
1605 ["DiscFac", "intertemporal discount factor", "β", False],
1606 ["Rfree", "risk free interest factor", "R", True],
1607 ["PermGroFac", "permanent income growth factor", "G", True],
1608 ["CRRA", "coefficient of relative risk aversion", "ρ", False],
1609 ["LivPrb", "survival probability", "ℒ", True],
1610 ["APFac", "absolute patience factor", "Þ=(βℒR)^(1/ρ)", False],
1611 ]
1613 param_desc = ""
1614 for j in range(len(params_to_describe)):
1615 this_entry = params_to_describe[j]
1616 if this_entry[3]:
1617 val = getattr(self, this_entry[0])[0]
1618 else:
1619 try:
1620 val = getattr(self, this_entry[0])
1621 except:
1622 val = self.bilt[this_entry[0]]
1623 this_line = (
1624 this_entry[2]
1625 + f"={val:.5f} : "
1626 + this_entry[1]
1627 + " ("
1628 + this_entry[0]
1629 + ")\n"
1630 )
1631 param_desc += this_line
1633 return param_desc
1635 def calc_limiting_values(self):
1636 """
1637 Compute various scalar values that are relevant to characterizing the
1638 solution to an infinite horizon problem. This method should only be called
1639 when T_cycle=1 and cycles=0, otherwise the values generated are meaningless.
1640 This method adds the following values to the instance in the dictionary
1641 attribute called bilt.
1643 APFac : Absolute Patience Factor
1644 GPFacRaw : Growth Patience Factor
1645 FHWFac : Finite Human Wealth Factor
1646 RPFac : Return Patience Factor
1647 PFVAFac : Perfect Foresight Value of Autarky Factor
1648 cNrmPDV : Present Discounted Value of Autarky Consumption
1649 MPCmin : Limiting minimum MPC as market resources go to infinity
1650 MPCmax : Limiting maximum MPC as market resources approach minimum level.
1651 hNrm : Human wealth divided by permanent income.
1652 Delta_mNrm_ZeroFunc : Linear consumption function where expected change in market resource ratio is zero
1653 BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income
1655 Returns
1656 -------
1657 None
1658 """
1659 aux_dict = self.bilt
1660 aux_dict["APFac"] = (self.Rfree[0] * self.DiscFac * self.LivPrb[0]) ** (
1661 1 / self.CRRA
1662 )
1663 aux_dict["GPFacRaw"] = aux_dict["APFac"] / self.PermGroFac[0]
1664 aux_dict["FHWFac"] = self.PermGroFac[0] / self.Rfree[0]
1665 aux_dict["RPFac"] = aux_dict["APFac"] / self.Rfree[0]
1666 aux_dict["PFVAFac"] = (self.DiscFac * self.LivPrb[0]) * self.PermGroFac[0] ** (
1667 1.0 - self.CRRA
1668 )
1669 aux_dict["cNrmPDV"] = 1.0 / (1.0 - aux_dict["RPFac"])
1670 aux_dict["MPCmin"] = np.maximum(1.0 - aux_dict["RPFac"], 0.0)
1671 constrained = (
1672 hasattr(self, "BoroCnstArt")
1673 and (self.BoroCnstArt is not None)
1674 and (self.BoroCnstArt > -np.inf)
1675 )
1677 if constrained:
1678 aux_dict["MPCmax"] = 1.0
1679 else:
1680 aux_dict["MPCmax"] = aux_dict["MPCmin"]
1681 if aux_dict["FHWFac"] < 1.0:
1682 aux_dict["hNrm"] = 1.0 / (1.0 - aux_dict["FHWFac"])
1683 else:
1684 aux_dict["hNrm"] = np.inf
1686 # Generate the "Delta m = 0" function, which is used to find target market resources
1687 Ex_Rnrm = self.Rfree[0] / self.PermGroFac[0]
1688 aux_dict["Delta_mNrm_ZeroFunc"] = (
1689 lambda m: (1.0 - 1.0 / Ex_Rnrm) * m + 1.0 / Ex_Rnrm
1690 )
1692 # Generate the "E[M_tp1 / M_t] = G" function, which is used to find balanced growth market resources
1693 PF_Rnrm = self.Rfree[0] / self.PermGroFac[0]
1694 aux_dict["BalGroFunc"] = lambda m: (1.0 - 1.0 / PF_Rnrm) * m + 1.0 / PF_Rnrm
1696 self.bilt = aux_dict
1698 def check_conditions(self, verbose=None):
1699 """
1700 This method checks whether the instance's type satisfies the
1701 Absolute Impatience Condition (AIC), the Return Impatience Condition (RIC),
1702 the Finite Human Wealth Condition (FHWC), the perfect foresight model's
1703 Growth Impatience Condition (GICRaw) and Perfect Foresight Finite Value
1704 of Autarky Condition (FVACPF). Depending on the configuration of parameter
1705 values, somecombination of these conditions must be satisfied in order
1706 for the problem to have a nondegenerate solution. To check which conditions
1707 are required, in the verbose mode a reference to the relevant theoretical
1708 literature is made.
1710 Parameters
1711 ----------
1712 verbose : boolean
1713 Specifies different levels of verbosity of feedback. When False, it
1714 only reports whether the instance's type fails to satisfy a particular
1715 condition. When True, it reports all results, i.e. the factor values
1716 for all conditions.
1718 Returns
1719 -------
1720 None
1721 """
1722 self.conditions = {}
1723 self.bilt["conditions_report"] = ""
1724 self.degenerate = False
1725 verbose = self.verbose if verbose is None else verbose
1727 # This method only checks for the conditions for infinite horizon models
1728 # with a 1 period cycle. If these conditions are not met, we exit early.
1729 if self.cycles != 0 or self.T_cycle > 1:
1730 trivial_message = "No conditions report was produced because this functionality is only supported for infinite horizon models with a cycle length of 1."
1731 self.log_condition_result(None, None, trivial_message, verbose)
1732 if not self.quiet:
1733 _log.info(self.bilt["conditions_report"])
1734 return
1736 # Calculate some useful quantities that will be used in the condition checks
1737 self.calc_limiting_values()
1738 param_desc = self.describe_parameters()
1739 self.log_condition_result(None, None, param_desc, verbose)
1741 # Check individual conditions and add their results to the report
1742 self.check_AIC(verbose)
1743 self.check_RIC(verbose)
1744 self.check_GICRaw(verbose)
1745 self.check_FVAC(verbose)
1746 self.check_FHWC(verbose)
1747 constrained = (
1748 hasattr(self, "BoroCnstArt")
1749 and (self.BoroCnstArt is not None)
1750 and (self.BoroCnstArt > -np.inf)
1751 )
1753 # Exit now if verbose output was not requested.
1754 if not verbose:
1755 if not self.quiet:
1756 _log.info(self.bilt["conditions_report"])
1757 return
1759 # Report on the degeneracy of the consumption function solution
1760 if not constrained:
1761 if self.conditions["FHWC"]:
1762 RIC_message = "\nBecause the FHWC is satisfied, the solution is not c(m)=Infinity."
1763 if self.conditions["RIC"]:
1764 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."
1765 degenerate = False
1766 else:
1767 RIC_message += " However, because the RIC is violated, the solution is degenerate at c(m) = 0 for all m."
1768 degenerate = True
1769 else:
1770 RIC_message = "\nBecause the FHWC condition is violated and the consumer is not constrained, the solution is degenerate at c(m)=Infinity."
1771 degenerate = True
1772 else:
1773 if self.conditions["RIC"]:
1774 RIC_message = "\nBecause the RIC is satisfied and the consumer is constrained, the solution is not c(m)=0 for all m."
1775 if self.conditions["GICRaw"]:
1776 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."
1777 degenerate = False
1778 else:
1779 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."
1780 degenerate = False
1781 else:
1782 if self.conditions["GICRaw"]:
1783 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."
1784 degenerate = False
1785 else:
1786 RIC_message = "\nBecause the RIC is violated but the FHWC is satisfied, the solution is degenerate at c(m)=0 for all m."
1787 degenerate = True
1788 self.log_condition_result(None, None, RIC_message, verbose)
1790 if (
1791 degenerate
1792 ): # All of the other checks are meaningless if the solution is degenerate
1793 if not self.quiet:
1794 _log.info(self.bilt["conditions_report"])
1795 return
1797 # Report on the consequences of the Absolute Impatience Condition
1798 if self.conditions["AIC"]:
1799 AIC_message = "\nBecause the AIC is satisfied, the absolute amount of consumption is expected to fall over time."
1800 else:
1801 AIC_message = "\nBecause the AIC is violated, the absolute amount of consumption is expected to grow over time."
1802 self.log_condition_result(None, None, AIC_message, verbose)
1804 # Report on the consequences of the Growth Impatience Condition
1805 if self.conditions["GICRaw"]:
1806 GIC_message = "\nBecause the GICRaw is satisfed, the ratio of individual wealth to permanent income is expected to fall indefinitely."
1807 elif self.conditions["FHWC"]:
1808 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."
1809 else:
1810 pass # pragma: nocover
1811 # This can never be reached! If GICRaw and FHWC both fail, then the RIC also fails, and we would have exited by this point.
1812 self.log_condition_result(None, None, GIC_message, verbose)
1814 if not self.quiet:
1815 _log.info(self.bilt["conditions_report"])
1817 def calc_stable_points(self, force=False):
1818 """
1819 If the problem is one that satisfies the conditions required for target ratios of different
1820 variables to permanent income to exist, and has been solved to within the self-defined
1821 tolerance, this method calculates the target values of market resources.
1823 Parameters
1824 ----------
1825 force : bool
1826 Indicator for whether the method should be forced to be run even if
1827 the agent seems to be the wrong type. Default is False.
1829 Returns
1830 -------
1831 None
1832 """
1833 # Child classes should not run this method
1834 is_perf_foresight = type(self) is PerfForesightConsumerType
1835 is_ind_shock = type(self) is IndShockConsumerType
1836 if not (is_perf_foresight or is_ind_shock or force):
1837 return
1839 infinite_horizon = self.cycles == 0
1840 single_period = self.T_cycle == 1
1841 if not infinite_horizon:
1842 raise ValueError(
1843 "The calc_stable_points method works only for infinite horizon models."
1844 )
1845 if not single_period:
1846 raise ValueError(
1847 "The calc_stable_points method works only with a single infinitely repeated period."
1848 )
1849 if not hasattr(self, "conditions"):
1850 raise ValueError(
1851 "The check_conditions method must be run before the calc_stable_points method."
1852 )
1853 if not hasattr(self, "solution"):
1854 raise ValueError(
1855 "The solve method must be run before the calc_stable_points method."
1856 )
1858 # Extract balanced growth and delta m_t+1 = 0 functions
1859 BalGroFunc = self.bilt["BalGroFunc"]
1860 Delta_mNrm_ZeroFunc = self.bilt["Delta_mNrm_ZeroFunc"]
1862 # If the GICRaw holds, then there is a balanced growth market resources ratio
1863 if self.conditions["GICRaw"]:
1864 cFunc = self.solution[0].cFunc
1865 func_to_zero = lambda m: BalGroFunc(m) - cFunc(m)
1866 m0 = 1.0
1867 try:
1868 mNrmStE = newton(func_to_zero, m0)
1869 except:
1870 mNrmStE = np.nan
1872 # A target level of assets *might* exist even if the GICMod fails, so check no matter what
1873 func_to_zero = lambda m: Delta_mNrm_ZeroFunc(m) - cFunc(m)
1874 m0 = 1.0 if np.isnan(mNrmStE) else mNrmStE
1875 try:
1876 mNrmTrg = newton(func_to_zero, m0, maxiter=200)
1877 except:
1878 mNrmTrg = np.nan
1879 else:
1880 mNrmStE = np.nan
1881 mNrmTrg = np.nan
1883 self.solution[0].mNrmStE = mNrmStE
1884 self.solution[0].mNrmTrg = mNrmTrg
1885 self.bilt["mNrmStE"] = mNrmStE
1886 self.bilt["mNrmTrg"] = mNrmTrg
1889###############################################################################
1891# Make a dictionary of constructors for the idiosyncratic income shocks model
1892IndShockConsumerType_constructors_default = {
1893 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1894 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1895 "IncShkDstn": construct_lognormal_income_process_unemployment,
1896 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1897 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1898 "aXtraGrid": make_assets_grid,
1899 "solution_terminal": make_basic_CRRA_solution_terminal,
1900}
1902# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1903IndShockConsumerType_kNrmInitDstn_default = {
1904 "kLogInitMean": -12.0, # Mean of log initial capital
1905 "kLogInitStd": 0.0, # Stdev of log initial capital
1906 "kNrmInitCount": 15, # Number of points in initial capital discretization
1907}
1909# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1910IndShockConsumerType_pLvlInitDstn_default = {
1911 "pLogInitMean": 0.0, # Mean of log permanent income
1912 "pLogInitStd": 0.0, # Stdev of log permanent income
1913 "pLvlInitCount": 15, # Number of points in initial capital discretization
1914}
1916# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
1917IndShockConsumerType_IncShkDstn_default = {
1918 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1919 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1920 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1921 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1922 "UnempPrb": 0.05, # Probability of unemployment while working
1923 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1924 "T_retire": 0, # Period of retirement (0 --> no retirement)
1925 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1926 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1927}
1929# Default parameters to make aXtraGrid using make_assets_grid
1930IndShockConsumerType_aXtraGrid_default = {
1931 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1932 "aXtraMax": 20, # Maximum end-of-period "assets above minimum" value
1933 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
1934 "aXtraCount": 48, # Number of points in the grid of "assets above minimum"
1935 "aXtraExtra": None, # Additional other values to add in grid (optional)
1936}
1938# Make a dictionary to specify an idiosyncratic income shocks consumer type
1939IndShockConsumerType_solving_default = {
1940 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1941 "cycles": 1, # Finite, non-cyclic model
1942 "T_cycle": 1, # Number of periods in the cycle for this agent type
1943 "pseudo_terminal": False, # Terminal period really does exist
1944 "constructors": IndShockConsumerType_constructors_default, # See dictionary above
1945 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
1946 "CRRA": 2.0, # Coefficient of relative risk aversion
1947 "Rfree": [1.03], # Interest factor on retained assets
1948 "DiscFac": 0.96, # Intertemporal discount factor
1949 "LivPrb": [0.98], # Survival probability after each period
1950 "PermGroFac": [1.01], # Permanent income growth factor
1951 "BoroCnstArt": 0.0, # Artificial borrowing constraint
1952 "vFuncBool": False, # Whether to calculate the value function during solution
1953 "CubicBool": False, # Whether to use cubic spline interpolation when True
1954 # (Uses linear spline interpolation for cFunc when False)
1955}
1956IndShockConsumerType_simulation_default = {
1957 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1958 "AgentCount": 10000, # Number of agents of this type
1959 "T_age": None, # Age after which simulated agents are automatically killed
1960 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1961 # (The portion of PermGroFac attributable to aggregate productivity growth)
1962 "NewbornTransShk": False, # Whether Newborns have transitory shock
1963 # ADDITIONAL OPTIONAL PARAMETERS
1964 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1965 # (Forces Newborns to follow solution path of the agent they replaced if True)
1966 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1967}
1969IndShockConsumerType_defaults = {}
1970IndShockConsumerType_defaults.update(IndShockConsumerType_IncShkDstn_default)
1971IndShockConsumerType_defaults.update(IndShockConsumerType_kNrmInitDstn_default)
1972IndShockConsumerType_defaults.update(IndShockConsumerType_pLvlInitDstn_default)
1973IndShockConsumerType_defaults.update(IndShockConsumerType_aXtraGrid_default)
1974IndShockConsumerType_defaults.update(IndShockConsumerType_solving_default)
1975IndShockConsumerType_defaults.update(IndShockConsumerType_simulation_default)
1976init_idiosyncratic_shocks = IndShockConsumerType_defaults # Here so that other models which use the old convention don't break
1979class IndShockConsumerType(PerfForesightConsumerType):
1980 r"""
1981 A consumer type with idiosyncratic shocks to permanent and transitory income.
1982 Their problem is defined by a sequence of income distributions, survival probabilities
1983 (:math:`\mathsf{S}`), and permanent income growth rates (:math:`\Gamma`), as well
1984 as time invariant values for risk aversion (:math:`\rho`), discount factor (:math:`\beta`),
1985 the interest rate (:math:`\mathsf{R}`), the grid of end-of-period assets, and an artificial
1986 borrowing constraint (:math:`\underline{a}`).
1988 .. math::
1989 \newcommand{\CRRA}{\rho}
1990 \newcommand{\LivPrb}{\mathsf{S}}
1991 \newcommand{\PermGroFac}{\Gamma}
1992 \newcommand{\Rfree}{\mathsf{R}}
1993 \newcommand{\DiscFac}{\beta}
1994 \begin{align*}
1995 v_t(m_t) &= \max_{c_t}u(c_t) + \DiscFac \LivPrb_t \mathbb{E}_{t} \left[ (\PermGroFac_{t+1} \psi_{t+1})^{1-\CRRA} v_{t+1}(m_{t+1}) \right], \\
1996 & \text{s.t.} \\
1997 a_t &= m_t - c_t, \\
1998 a_t &\geq \underline{a}, \\
1999 m_{t+1} &= a_t \Rfree_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
2000 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
2001 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1, \\
2002 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}.
2003 \end{align*}
2006 Constructors
2007 ------------
2008 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2009 The agent's income shock distributions.
2011 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2012 aXtraGrid: Constructor
2013 The agent's asset grid.
2015 Its default constructor is :func:`HARK.utilities.make_assets_grid`
2017 Solving Parameters
2018 ------------------
2019 cycles: int
2020 0 specifies an infinite horizon model, 1 specifies a finite model.
2021 T_cycle: int
2022 Number of periods in the cycle for this agent type.
2023 CRRA: float, :math:`\rho`
2024 Coefficient of Relative Risk Aversion.
2025 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
2026 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
2027 DiscFac: float, :math:`\beta`
2028 Intertemporal discount factor.
2029 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2030 Survival probability after each period.
2031 PermGroFac: list[float], time varying, :math:`\Gamma`
2032 Permanent income growth factor.
2033 BoroCnstArt: float, :math:`\underline{a}`
2034 The minimum Asset/Perminant Income ratio, None to ignore.
2035 vFuncBool: bool
2036 Whether to calculate the value function during solution.
2037 CubicBool: bool
2038 Whether to use cubic spline interpoliation.
2040 Simulation Parameters
2041 ---------------------
2042 AgentCount: int
2043 Number of agents of this kind that are created during simulations.
2044 T_age: int
2045 Age after which to automatically kill agents, None to ignore.
2046 T_sim: int, required for simulation
2047 Number of periods to simulate.
2048 track_vars: list[strings]
2049 List of variables that should be tracked when running the simulation.
2050 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2052 PermShk is the agent's permanent income shock
2054 TranShk is the agent's transitory income shock
2056 aLvl is the nominal asset level
2058 aNrm is the normalized assets
2060 bNrm is the normalized resources without this period's labor income
2062 cNrm is the normalized consumption
2064 mNrm is the normalized market resources
2066 pLvl is the permanent income level
2068 who_dies is the array of which agents died
2069 aNrmInitMean: float
2070 Mean of Log initial Normalized Assets.
2071 aNrmInitStd: float
2072 Std of Log initial Normalized Assets.
2073 pLvlInitMean: float
2074 Mean of Log initial permanent income.
2075 pLvlInitStd: float
2076 Std of Log initial permanent income.
2077 PermGroFacAgg: float
2078 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2079 PerfMITShk: boolean
2080 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2081 NewbornTransShk: boolean
2082 Whether Newborns have transitory shock.
2084 Attributes
2085 ----------
2086 solution: list[Consumer solution object]
2087 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2088 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2090 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2091 history: Dict[Array]
2092 Created by running the :func:`.simulate()` method.
2093 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2094 Visit :class:`HARK.core.AgentType.simulate` for more information.
2095 """
2097 IncShkDstn_defaults = IndShockConsumerType_IncShkDstn_default
2098 aXtraGrid_defaults = IndShockConsumerType_aXtraGrid_default
2099 solving_defaults = IndShockConsumerType_solving_default
2100 simulation_defaults = IndShockConsumerType_simulation_default
2101 default_ = {
2102 "params": IndShockConsumerType_defaults,
2103 "solver": solve_one_period_ConsIndShock,
2104 "model": "ConsIndShock.yaml",
2105 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
2106 }
2108 time_inv_ = PerfForesightConsumerType.time_inv_ + [
2109 "vFuncBool",
2110 "CubicBool",
2111 "aXtraGrid",
2112 ]
2113 time_vary_ = PerfForesightConsumerType.time_vary_ + [
2114 "IncShkDstn",
2115 "PermShkDstn",
2116 "TranShkDstn",
2117 ]
2118 # This is in the PerfForesight model but not ConsIndShock
2119 time_inv_.remove("MaxKinks")
2120 shock_vars_ = ["PermShk", "TranShk"]
2121 distributions = [
2122 "IncShkDstn",
2123 "PermShkDstn",
2124 "TranShkDstn",
2125 "kNrmInitDstn",
2126 "pLvlInitDstn",
2127 ]
2129 def update_income_process(self):
2130 self.update("IncShkDstn", "PermShkDstn", "TranShkDstn")
2132 def get_shocks(self):
2133 """
2134 Gets permanent and transitory income shocks for this period. Samples from IncShkDstn for
2135 each period in the cycle.
2137 Parameters
2138 ----------
2139 NewbornTransShk : boolean, optional
2140 Whether Newborns have transitory shock. The default is False.
2142 Returns
2143 -------
2144 None
2145 """
2146 # Whether Newborns have transitory shock. The default is False.
2147 NewbornTransShk = self.NewbornTransShk
2149 PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
2150 TranShkNow = np.zeros(self.AgentCount)
2151 newborn = self.t_age == 0
2152 for t in np.unique(self.t_cycle):
2153 idx = self.t_cycle == t
2155 # temporary, see #1022
2156 if self.cycles == 1:
2157 t = t - 1
2159 N = np.sum(idx)
2160 if N > 0:
2161 # set current income distribution
2162 IncShkDstnNow = self.IncShkDstn[t]
2163 # and permanent growth factor
2164 PermGroFacNow = self.PermGroFac[t]
2165 # Get random draws of income shocks from the discrete distribution
2166 IncShks = IncShkDstnNow.draw(N)
2168 PermShkNow[idx] = (
2169 IncShks[0, :] * PermGroFacNow
2170 ) # permanent "shock" includes expected growth
2171 TranShkNow[idx] = IncShks[1, :]
2173 # That procedure used the *last* period in the sequence for newborns, but that's not right
2174 # Redraw shocks for newborns, using the *first* period in the sequence. Approximation.
2175 N = np.sum(newborn)
2176 if N > 0:
2177 idx = newborn
2178 # set current income distribution
2179 IncShkDstnNow = self.IncShkDstn[0]
2180 PermGroFacNow = self.PermGroFac[0] # and permanent growth factor
2182 # Get random draws of income shocks from the discrete distribution
2183 EventDraws = IncShkDstnNow.draw_events(N)
2184 PermShkNow[idx] = (
2185 IncShkDstnNow.atoms[0][EventDraws] * PermGroFacNow
2186 ) # permanent "shock" includes expected growth
2187 TranShkNow[idx] = IncShkDstnNow.atoms[1][EventDraws]
2189 # Whether Newborns have transitory shock. The default is False.
2190 if not NewbornTransShk:
2191 TranShkNow[newborn] = 1.0
2193 # Store the shocks in self
2194 self.shocks["PermShk"] = PermShkNow
2195 self.shocks["TranShk"] = TranShkNow
2197 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
2198 """
2199 Creates a "normalized Euler error" function for this instance, mapping
2200 from market resources to "consumption error per dollar of consumption."
2201 Stores result in attribute eulerErrorFunc as an interpolated function.
2202 Has option to use approximate income distribution stored in self.IncShkDstn
2203 or to use a (temporary) very dense approximation.
2205 Only works on (one period) infinite horizon models at this time, will
2206 be generalized later.
2208 Parameters
2209 ----------
2210 mMax : float
2211 Maximum normalized market resources for the Euler error function.
2212 approx_inc_dstn : Boolean
2213 Indicator for whether to use the approximate discrete income distri-
2214 bution stored in self.IncShkDstn[0], or to use a very accurate
2215 discrete approximation instead. When True, uses approximation in
2216 IncShkDstn; when False, makes and uses a very dense approximation.
2218 Returns
2219 -------
2220 None
2222 Notes
2223 -----
2224 This method is not used by any other code in the library. Rather, it is here
2225 for expository and benchmarking purposes.
2226 """
2227 # Get the income distribution (or make a very dense one)
2228 if approx_inc_dstn:
2229 IncShkDstn = self.IncShkDstn[0]
2230 else:
2231 TranShkDstn = MeanOneLogNormal(sigma=self.TranShkStd[0]).discretize(
2232 N=200,
2233 method="equiprobable",
2234 tail_N=50,
2235 tail_order=1.3,
2236 tail_bound=[0.05, 0.95],
2237 )
2238 TranShkDstn = add_discrete_outcome_constant_mean(
2239 TranShkDstn, p=self.UnempPrb, x=self.IncUnemp
2240 )
2241 PermShkDstn = MeanOneLogNormal(sigma=self.PermShkStd[0]).discretize(
2242 N=200,
2243 method="equiprobable",
2244 tail_N=50,
2245 tail_order=1.3,
2246 tail_bound=[0.05, 0.95],
2247 )
2248 IncShkDstn = combine_indep_dstns(PermShkDstn, TranShkDstn)
2250 # Make a grid of market resources
2251 mNowMin = self.solution[0].mNrmMin + 10 ** (
2252 -15
2253 ) # add tiny bit to get around 0/0 problem
2254 mNowMax = mMax
2255 mNowGrid = np.linspace(mNowMin, mNowMax, 1000)
2257 # Get the consumption function this period and the marginal value function
2258 # for next period. Note that this part assumes a one period cycle.
2259 cFuncNow = self.solution[0].cFunc
2260 vPfuncNext = self.solution[0].vPfunc
2262 # Calculate consumption this period at each gridpoint (and assets)
2263 cNowGrid = cFuncNow(mNowGrid)
2264 aNowGrid = mNowGrid - cNowGrid
2266 # Tile the grids for fast computation
2267 ShkCount = IncShkDstn.pmv.size
2268 aCount = aNowGrid.size
2269 aNowGrid_tiled = np.tile(aNowGrid, (ShkCount, 1))
2270 PermShkVals_tiled = (np.tile(IncShkDstn.atoms[0], (aCount, 1))).transpose()
2271 TranShkVals_tiled = (np.tile(IncShkDstn.atoms[1], (aCount, 1))).transpose()
2272 ShkPrbs_tiled = (np.tile(IncShkDstn.pmv, (aCount, 1))).transpose()
2274 # Calculate marginal value next period for each gridpoint and each shock
2275 mNextArray = (
2276 self.Rfree[0] / (self.PermGroFac[0] * PermShkVals_tiled) * aNowGrid_tiled
2277 + TranShkVals_tiled
2278 )
2279 vPnextArray = vPfuncNext(mNextArray)
2281 # Calculate expected marginal value and implied optimal consumption
2282 ExvPnextGrid = (
2283 self.DiscFac
2284 * self.Rfree[0]
2285 * self.LivPrb[0]
2286 * self.PermGroFac[0] ** (-self.CRRA)
2287 * np.sum(
2288 PermShkVals_tiled ** (-self.CRRA) * vPnextArray * ShkPrbs_tiled, axis=0
2289 )
2290 )
2291 cOptGrid = ExvPnextGrid ** (
2292 -1.0 / self.CRRA
2293 ) # This is the 'Endogenous Gridpoints' step
2295 # Calculate Euler error and store an interpolated function
2296 EulerErrorNrmGrid = (cNowGrid - cOptGrid) / cOptGrid
2297 eulerErrorFunc = LinearInterp(mNowGrid, EulerErrorNrmGrid)
2298 self.eulerErrorFunc = eulerErrorFunc
2300 def pre_solve(self):
2301 self.check_restrictions()
2302 self.construct("solution_terminal")
2303 if not self.quiet:
2304 self.check_conditions(verbose=self.verbose)
2306 def describe_parameters(self):
2307 """
2308 Generate a string describing the primitive model parameters that will
2309 be used to calculating limiting values and factors.
2311 Parameters
2312 ----------
2313 None
2315 Returns
2316 -------
2317 param_desc : str
2318 Description of primitive parameters.
2319 """
2320 # Get parameter description from the perfect foresight model
2321 param_desc = super().describe_parameters()
2323 # Make a new entry for weierstrass-p (the weird formatting here is to
2324 # make it easier to adapt into the style of the superclass if we add more
2325 # parameter reports later)
2326 this_entry = [
2327 "WorstPrb",
2328 "probability of worst income shock realization",
2329 "℘",
2330 False,
2331 ]
2332 try:
2333 val = getattr(self, this_entry[0])
2334 except:
2335 val = self.bilt[this_entry[0]]
2336 this_line = (
2337 this_entry[2]
2338 + f"={val:.5f} : "
2339 + this_entry[1]
2340 + " ("
2341 + this_entry[0]
2342 + ")\n"
2343 )
2345 # Add in the new entry and return it
2346 param_desc += this_line
2347 return param_desc
2349 def calc_limiting_values(self):
2350 """
2351 Compute various scalar values that are relevant to characterizing the
2352 solution to an infinite horizon problem. This method should only be called
2353 when T_cycle=1 and cycles=0, otherwise the values generated are meaningless.
2354 This method adds the following values to this instance in the dictionary
2355 attribute called bilt.
2357 APFac : Absolute Patience Factor
2358 GPFacRaw : Growth Patience Factor
2359 GPFacMod : Risk-Modified Growth Patience Factor
2360 GPFacLiv : Mortality-Adjusted Growth Patience Factor
2361 GPFacLivMod : Modigliani Mortality-Adjusted Growth Patience Factor
2362 GPFacSdl : Szeidl Growth Patience Factor
2363 FHWFac : Finite Human Wealth Factor
2364 RPFac : Return Patience Factor
2365 WRPFac : Weak Return Patience Factor
2366 PFVAFac : Perfect Foresight Value of Autarky Factor
2367 VAFac : Value of Autarky Factor
2368 cNrmPDV : Present Discounted Value of Autarky Consumption
2369 MPCmin : Limiting minimum MPC as market resources go to infinity
2370 MPCmax : Limiting maximum MPC as market resources approach minimum level
2371 hNrm : Human wealth divided by permanent income.
2372 ELogPermShk : Expected log permanent income shock
2373 WorstPrb : Probability of worst income shock realization
2374 Delta_mNrm_ZeroFunc : Linear locus where expected change in market resource ratio is zero
2375 BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income
2377 Returns
2378 -------
2379 None
2380 """
2381 super().calc_limiting_values()
2382 aux_dict = self.bilt
2384 # Calculate the risk-modified growth impatience factor
2385 PermShkDstn = self.PermShkDstn[0]
2386 inv_func = lambda x: x ** (-1.0)
2387 Ex_PermShkInv = expected(inv_func, PermShkDstn)[0]
2388 GroCompPermShk = Ex_PermShkInv ** (-1.0)
2389 aux_dict["GPFacMod"] = aux_dict["APFac"] / (self.PermGroFac[0] * GroCompPermShk)
2391 # Calculate the mortality-adjusted growth impatience factor (and version
2392 # with Modigiliani bequests)
2393 aux_dict["GPFacLiv"] = aux_dict["GPFacRaw"] * self.LivPrb[0]
2394 aux_dict["GPFacLivMod"] = aux_dict["GPFacLiv"] * self.LivPrb[0]
2396 # Calculate the risk-modified value of autarky factor
2397 if self.CRRA == 1.0:
2398 UtilCompPermShk = np.exp(expected(np.log, PermShkDstn)[0])
2399 else:
2400 CRRAfunc = lambda x: x ** (1.0 - self.CRRA)
2401 UtilCompPermShk = expected(CRRAfunc, PermShkDstn)[0] ** (
2402 1 / (1.0 - self.CRRA)
2403 )
2404 aux_dict["VAFac"] = self.DiscFac * (self.PermGroFac[0] * UtilCompPermShk) ** (
2405 1.0 - self.CRRA
2406 )
2408 # Calculate the expected log permanent income shock, which will be used
2409 # for the Szeidl variation of the Growth Impatience condition
2410 aux_dict["ELogPermShk"] = expected(np.log, PermShkDstn)[0]
2412 # Calculate the Harmenberg permanent income neutral expected log permanent
2413 # shock and the Harmenberg Growth Patience Factor
2414 Hrm_func = lambda x: x * np.log(x)
2415 PermShk_Hrm = np.exp(expected(Hrm_func, PermShkDstn)[0])
2416 aux_dict["GPFacHrm"] = aux_dict["GPFacRaw"] / PermShk_Hrm
2418 # Calculate the probability of the worst income shock realization
2419 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2420 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2421 ShkPrbsNext = self.IncShkDstn[0].pmv
2422 Ex_IncNext = np.dot(ShkPrbsNext, PermShkValsNext * TranShkValsNext)
2423 PermShkMinNext = np.min(PermShkValsNext)
2424 TranShkMinNext = np.min(TranShkValsNext)
2425 WorstIncNext = PermShkMinNext * TranShkMinNext
2426 WorstIncPrb = np.sum(
2427 ShkPrbsNext[(PermShkValsNext * TranShkValsNext) == WorstIncNext]
2428 )
2429 aux_dict["WorstPrb"] = WorstIncPrb
2431 # Calculate the weak return patience factor
2432 aux_dict["WRPFac"] = WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2434 # Calculate human wealth and the infinite horizon natural borrowing constraint
2435 if aux_dict["FHWFac"] < 1.0:
2436 hNrm = Ex_IncNext / (1.0 - aux_dict["FHWFac"])
2437 else:
2438 hNrm = np.inf
2439 temp = PermShkMinNext * aux_dict["FHWFac"]
2440 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2442 # Find the upper bound of the MPC as market resources approach the minimum
2443 BoroCnstArt = -np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2444 if BoroCnstNat < BoroCnstArt:
2445 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2446 else:
2447 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2448 MPCmax = np.maximum(MPCmax, 0.0)
2450 # Store maximum MPC and human wealth
2451 aux_dict["hNrm"] = hNrm
2452 aux_dict["MPCmax"] = MPCmax
2454 # Generate the "Delta m = 0" function, which is used to find target market resources
2455 # This overwrites the function generated by the perfect foresight version
2456 Ex_Rnrm = self.Rfree[0] / self.PermGroFac[0] * Ex_PermShkInv
2457 aux_dict["Delta_mNrm_ZeroFunc"] = (
2458 lambda m: (1.0 - 1.0 / Ex_Rnrm) * m + 1.0 / Ex_Rnrm
2459 )
2461 self.bilt = aux_dict
2463 self.bilt = aux_dict
2465 def check_GICMod(self, verbose=None):
2466 """
2467 Evaluate and report on the Risk-Modified Growth Impatience Condition.
2468 """
2469 name = "GICMod"
2470 GPFacMod = self.bilt["GPFacMod"]
2471 result = GPFacMod < 1.0
2473 messages = {
2474 True: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor satisfies the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2475 False: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor violates the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2476 }
2477 verbose = self.verbose if verbose is None else verbose
2478 self.log_condition_result(name, result, messages[result], verbose)
2480 def check_GICSdl(self, verbose=None):
2481 """
2482 Evaluate and report on the Szeidl variation of the Growth Impatience Condition.
2483 """
2484 name = "GICSdl"
2485 ELogPermShk = self.bilt["ELogPermShk"]
2486 result = np.log(self.bilt["GPFacRaw"]) < ELogPermShk
2488 messages = {
2489 True: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock satisfies the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2490 False: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock violates the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2491 }
2492 verbose = self.verbose if verbose is None else verbose
2493 self.log_condition_result(name, result, messages[result], verbose)
2495 def check_GICHrm(self, verbose=None):
2496 """
2497 Evaluate and report on the Harmenberg variation of the Growth Impatience Condition.
2498 """
2499 name = "GICHrm"
2500 GPFacHrm = self.bilt["GPFacHrm"]
2501 result = GPFacHrm < 1.0
2503 messages = {
2504 True: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor satisfies the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2505 False: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor violates the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2506 }
2507 verbose = self.verbose if verbose is None else verbose
2508 self.log_condition_result(name, result, messages[result], verbose)
2510 def check_GICLiv(self, verbose=None):
2511 """
2512 Evaluate and report on the Mortality-Adjusted Growth Impatience Condition.
2513 """
2514 name = "GICLiv"
2515 GPFacLiv = self.bilt["GPFacLiv"]
2516 result = GPFacLiv < 1.0
2518 messages = {
2519 True: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor satisfies the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2520 False: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor violates the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2521 }
2522 verbose = self.verbose if verbose is None else verbose
2523 self.log_condition_result(name, result, messages[result], verbose)
2525 def check_FVAC(self, verbose=None):
2526 """
2527 Evaluate and report on the Finite Value of Autarky condition in the presence of income risk.
2528 """
2529 name = "FVAC"
2530 VAFac = self.bilt["VAFac"]
2531 result = VAFac < 1.0
2533 messages = {
2534 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.",
2535 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.",
2536 }
2537 verbose = self.verbose if verbose is None else verbose
2538 self.log_condition_result(name, result, messages[result], verbose)
2540 def check_WRIC(self, verbose=None):
2541 """
2542 Evaluate and report on the Weak Return Impatience Condition.
2543 """
2544 name = "WRIC"
2545 WRPFac = self.bilt["WRPFac"]
2546 result = WRPFac < 1.0
2548 messages = {
2549 True: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor satisfies the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2550 False: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor violates the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2551 }
2552 verbose = self.verbose if verbose is None else verbose
2553 self.log_condition_result(name, result, messages[result], verbose)
2555 def check_conditions(self, verbose=None):
2556 """
2557 This method checks whether the instance's type satisfies various conditions.
2558 When combinations of these conditions are satisfied, the solution to the
2559 problem exhibits different characteristics. (For an exposition of the
2560 conditions, see https://econ-ark.github.io/BufferStockTheory/)
2562 Parameters
2563 ----------
2564 verbose : boolean
2565 Specifies different levels of verbosity of feedback. When False, it only reports whether the
2566 instance's type fails to satisfy a particular condition. When True, it reports all results, i.e.
2567 the factor values for all conditions.
2569 Returns
2570 -------
2571 None
2572 """
2573 self.conditions = {}
2574 self.bilt["conditions_report"] = ""
2575 self.degenerate = False
2576 verbose = self.verbose if verbose is None else verbose
2578 # This method only checks for the conditions for infinite horizon models
2579 # with a 1 period cycle. If these conditions are not met, we exit early.
2580 if self.cycles != 0 or self.T_cycle > 1:
2581 trivial_message = "No conditions report was produced because this functionality is only supported for infinite horizon models with a cycle length of 1."
2582 self.log_condition_result(None, None, trivial_message, verbose)
2583 if not self.quiet:
2584 _log.info(self.bilt["conditions_report"])
2585 return
2587 # Calculate some useful quantities that will be used in the condition checks
2588 self.calc_limiting_values()
2589 param_desc = self.describe_parameters()
2590 self.log_condition_result(None, None, param_desc, verbose)
2592 # Check individual conditions and add their results to the report
2593 self.check_AIC(verbose)
2594 self.check_RIC(verbose)
2595 self.check_WRIC(verbose)
2596 self.check_GICRaw(verbose)
2597 self.check_GICMod(verbose)
2598 self.check_GICLiv(verbose)
2599 self.check_GICSdl(verbose)
2600 self.check_GICHrm(verbose)
2601 super().check_FVAC(verbose)
2602 self.check_FVAC(verbose)
2603 self.check_FHWC(verbose)
2605 # Exit now if verbose output was not requested.
2606 if not verbose:
2607 if not self.quiet:
2608 _log.info(self.bilt["conditions_report"])
2609 return
2611 # Report on the degeneracy of the consumption function solution
2612 if self.conditions["WRIC"] and self.conditions["FVAC"]:
2613 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."
2614 degenerate = False
2615 elif not self.conditions["WRIC"]:
2616 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"
2617 degenerate = True
2618 elif not self.conditions["FVAC"]:
2619 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."
2620 degenerate = False
2621 self.log_condition_result(None, None, degen_message, verbose)
2622 self.degenerate = degenerate
2624 # Stop here if the solution is degenerate
2625 if degenerate:
2626 if not self.quiet:
2627 _log.info(self.bilt["conditions_report"])
2628 return
2630 # Report on the limiting behavior of the consumption function as m goes to infinity
2631 if self.conditions["RIC"]:
2632 if self.conditions["FHWC"]:
2633 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."
2634 else:
2635 RIC_message = "\nBecause the RIC is satisfied but the FHWC is violated, the GIC is satisfied."
2636 else:
2637 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."
2638 self.log_condition_result(None, None, RIC_message, verbose)
2640 # Report on whether a pseudo-steady-state exists at the individual level
2641 if self.conditions["GICRaw"]:
2642 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."
2643 else:
2644 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."
2645 self.log_condition_result(None, None, GIC_message, verbose)
2647 # Report on whether a target wealth ratio exists at the individual level
2648 if self.conditions["GICMod"]:
2649 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."
2650 else:
2651 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."
2652 self.log_condition_result(None, None, GICMod_message, verbose)
2654 # Report on whether a target level of wealth exists at the aggregate level
2655 if self.conditions["GICLiv"]:
2656 GICLiv_message = "\nBecause the GICLiv is satisfied, a target ratio of aggregate market resources to aggregate permanent income exists."
2657 else:
2658 GICLiv_message = "\nBecause the GICLiv is violated, a target ratio of aggregate market resources to aggregate permanent income might not exist."
2659 self.log_condition_result(None, None, GICLiv_message, verbose)
2661 # Report on whether invariant distributions exist
2662 if self.conditions["GICSdl"]:
2663 GICSdl_message = "\nBecause the GICSdl is satisfied, there exist invariant distributions of permanent income-normalized variables."
2664 else:
2665 GICSdl_message = "\nBecause the GICSdl is violated, there do not exist invariant distributions of permanent income-normalized variables."
2666 self.log_condition_result(None, None, GICSdl_message, verbose)
2668 # Report on whether blah blah
2669 if self.conditions["GICHrm"]:
2670 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."
2671 else:
2672 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.."
2673 self.log_condition_result(None, None, GICHrm_message, verbose)
2675 if not self.quiet:
2676 _log.info(self.bilt["conditions_report"])
2679###############################################################################
2681# Specify default parameters used in "kinked R" model
2683KinkedRconsumerType_IncShkDstn_default = IndShockConsumerType_IncShkDstn_default.copy()
2684KinkedRconsumerType_aXtraGrid_default = IndShockConsumerType_aXtraGrid_default.copy()
2685KinkedRconsumerType_kNrmInitDstn_default = (
2686 IndShockConsumerType_kNrmInitDstn_default.copy()
2687)
2688KinkedRconsumerType_pLvlInitDstn_default = (
2689 IndShockConsumerType_pLvlInitDstn_default.copy()
2690)
2692KinkedRconsumerType_solving_default = IndShockConsumerType_solving_default.copy()
2693KinkedRconsumerType_solving_default.update(
2694 {
2695 "Rboro": 1.20, # Interest factor on assets when borrowing, a < 0
2696 "Rsave": 1.02, # Interest factor on assets when saving, a > 0
2697 "BoroCnstArt": None, # Kinked R only matters if borrowing is allowed
2698 }
2699)
2700del KinkedRconsumerType_solving_default["Rfree"]
2702KinkedRconsumerType_simulation_default = IndShockConsumerType_simulation_default.copy()
2704KinkedRconsumerType_defaults = {}
2705KinkedRconsumerType_defaults.update(
2706 KinkedRconsumerType_IncShkDstn_default
2707) # Fill with some parameters
2708KinkedRconsumerType_defaults.update(KinkedRconsumerType_pLvlInitDstn_default)
2709KinkedRconsumerType_defaults.update(KinkedRconsumerType_kNrmInitDstn_default)
2710KinkedRconsumerType_defaults.update(KinkedRconsumerType_aXtraGrid_default)
2711KinkedRconsumerType_defaults.update(KinkedRconsumerType_solving_default)
2712KinkedRconsumerType_defaults.update(KinkedRconsumerType_simulation_default)
2713init_kinked_R = KinkedRconsumerType_defaults
2716class KinkedRconsumerType(IndShockConsumerType):
2717 r"""
2718 A consumer type based on IndShockConsumerType, with different
2719 interest rates for saving (:math:`\mathsf{R}_{save}`) and borrowing
2720 (:math:`\mathsf{R}_{boro}`).
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*}
2744 Constructors
2745 ------------
2746 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2747 The agent's income shock distributions.
2748 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2749 aXtraGrid: Constructor
2750 The agent's asset grid.
2751 Its default constructor is :func:`HARK.utilities.make_assets_grid`
2753 Solving Parameters
2754 ------------------
2755 cycles: int
2756 0 specifies an infinite horizon model, 1 specifies a finite model.
2757 T_cycle: int
2758 Number of periods in the cycle for this agent type.
2759 CRRA: float, :math:`\rho`
2760 Coefficient of Relative Risk Aversion.
2761 Rboro: float, :math:`\mathsf{R}_{boro}`
2762 Risk Free interest rate when assets are negative.
2763 Rsave: float, :math:`\mathsf{R}_{save}`
2764 Risk Free interest rate when assets are positive.
2765 DiscFac: float, :math:`\beta`
2766 Intertemporal discount factor.
2767 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2768 Survival probability after each period.
2769 PermGroFac: list[float], time varying, :math:`\Gamma`
2770 Permanent income growth factor.
2771 BoroCnstArt: float, :math:`\underline{a}`
2772 The minimum Asset/Perminant Income ratio, None to ignore.
2773 vFuncBool: bool
2774 Whether to calculate the value function during solution.
2775 CubicBool: bool
2776 Whether to use cubic spline interpoliation.
2778 Simulation Parameters
2779 ---------------------
2780 AgentCount: int
2781 Number of agents of this kind that are created during simulations.
2782 T_age: int
2783 Age after which to automatically kill agents, None to ignore.
2784 T_sim: int, required for simulation
2785 Number of periods to simulate.
2786 track_vars: list[strings]
2787 List of variables that should be tracked when running the simulation.
2788 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2790 PermShk is the agent's permanent income shock
2792 TranShk is the agent's transitory income shock
2794 aLvl is the nominal asset level
2796 aNrm is the normalized assets
2798 bNrm is the normalized resources without this period's labor income
2800 cNrm is the normalized consumption
2802 mNrm is the normalized market resources
2804 pLvl is the permanent income level
2806 who_dies is the array of which agents died
2807 aNrmInitMean: float
2808 Mean of Log initial Normalized Assets.
2809 aNrmInitStd: float
2810 Std of Log initial Normalized Assets.
2811 pLvlInitMean: float
2812 Mean of Log initial permanent income.
2813 pLvlInitStd: float
2814 Std of Log initial permanent income.
2815 PermGroFacAgg: float
2816 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2817 PerfMITShk: boolean
2818 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2819 NewbornTransShk: boolean
2820 Whether Newborns have transitory shock.
2822 Attributes
2823 ----------
2824 solution: list[Consumer solution object]
2825 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2826 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2828 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2829 history: Dict[Array]
2830 Created by running the :func:`.simulate()` method.
2831 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2832 Visit :class:`HARK.core.AgentType.simulate` for more information.
2833 """
2835 IncShkDstn_defaults = KinkedRconsumerType_IncShkDstn_default
2836 aXtraGrid_defaults = KinkedRconsumerType_aXtraGrid_default
2837 solving_defaults = KinkedRconsumerType_solving_default
2838 simulation_defaults = KinkedRconsumerType_simulation_default
2839 default_ = {
2840 "params": KinkedRconsumerType_defaults,
2841 "solver": solve_one_period_ConsKinkedR,
2842 "model": "ConsKinkedR.yaml",
2843 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
2844 }
2846 time_inv_ = copy(IndShockConsumerType.time_inv_)
2847 time_inv_ += ["Rboro", "Rsave"]
2849 def calc_bounding_values(self):
2850 """
2851 Calculate human wealth plus minimum and maximum MPC in an infinite
2852 horizon model with only one period repeated indefinitely. Store results
2853 as attributes of self. Human wealth is the present discounted value of
2854 expected future income after receiving income this period, ignoring mort-
2855 ality. The maximum MPC is the limit of the MPC as m --> mNrmMin. The
2856 minimum MPC is the limit of the MPC as m --> infty. This version deals
2857 with the different interest rates on borrowing vs saving.
2859 Parameters
2860 ----------
2861 None
2863 Returns
2864 -------
2865 None
2866 """
2867 # Unpack the income distribution and get average and worst outcomes
2868 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2869 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2870 ShkPrbsNext = self.IncShkDstn[0].pmv
2871 IncNext = PermShkValsNext * TranShkValsNext
2872 Ex_IncNext = np.dot(ShkPrbsNext, IncNext)
2873 PermShkMinNext = np.min(PermShkValsNext)
2874 TranShkMinNext = np.min(TranShkValsNext)
2875 WorstIncNext = PermShkMinNext * TranShkMinNext
2876 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
2877 # TODO: Check the math above. I think it fails for non-independent shocks
2879 BoroCnstArt = np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2881 # Calculate human wealth and the infinite horizon natural borrowing constraint
2882 hNrm = (Ex_IncNext * self.PermGroFac[0] / self.Rsave) / (
2883 1.0 - self.PermGroFac[0] / self.Rsave
2884 )
2885 temp = self.PermGroFac[0] * PermShkMinNext / self.Rboro
2886 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2888 PatFacTop = (self.DiscFac * self.LivPrb[0] * self.Rsave) ** (
2889 1.0 / self.CRRA
2890 ) / self.Rsave
2891 PatFacBot = (self.DiscFac * self.LivPrb[0] * self.Rboro) ** (
2892 1.0 / self.CRRA
2893 ) / self.Rboro
2894 if BoroCnstNat < BoroCnstArt:
2895 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2896 else:
2897 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * PatFacBot
2898 MPCmin = 1.0 - PatFacTop
2900 # Store the results as attributes of self
2901 self.hNrm = hNrm
2902 self.MPCmin = MPCmin
2903 self.MPCmax = MPCmax
2905 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True): # pragma: nocover
2906 """
2907 Creates a "normalized Euler error" function for this instance, mapping
2908 from market resources to "consumption error per dollar of consumption."
2909 Stores result in attribute eulerErrorFunc as an interpolated function.
2910 Has option to use approximate income distribution stored in self.IncShkDstn
2911 or to use a (temporary) very dense approximation.
2913 SHOULD BE INHERITED FROM ConsIndShockModel
2915 Parameters
2916 ----------
2917 mMax : float
2918 Maximum normalized market resources for the Euler error function.
2919 approx_inc_dstn : Boolean
2920 Indicator for whether to use the approximate discrete income distri-
2921 bution stored in self.IncShkDstn[0], or to use a very accurate
2922 discrete approximation instead. When True, uses approximation in
2923 IncShkDstn; when False, makes and uses a very dense approximation.
2925 Returns
2926 -------
2927 None
2928 """
2929 raise NotImplementedError()
2931 def get_Rport(self):
2932 """
2933 Returns an array of size self.AgentCount with self.Rboro or self.Rsave in
2934 each entry, based on whether self.aNrmNow >< 0. This represents the risk-
2935 free portfolio return in this model.
2937 Parameters
2938 ----------
2939 None
2941 Returns
2942 -------
2943 RfreeNow : np.array
2944 Array of size self.AgentCount with risk free interest rate for each agent.
2945 """
2946 RfreeNow = self.Rboro * np.ones(self.AgentCount)
2947 RfreeNow[self.state_prev["aNrm"] > 0] = self.Rsave
2948 return RfreeNow
2950 def check_conditions(self, verbose):
2951 """
2952 This empty method overwrites the version inherited from its parent class,
2953 IndShockConsumerType. The condition checks are not appropriate when Rfree
2954 has multiple values.
2956 Parameters
2957 ----------
2958 None
2960 Returns
2961 -------
2962 None
2963 """
2964 pass
2967###############################################################################
2969# Make a dictionary to specify a lifecycle consumer with a finite horizon
2971# Main calibration characteristics
2972birth_age = 25
2973death_age = 90
2974adjust_infl_to = 1992
2975# Use income estimates from Cagetti (2003) for High-school graduates
2976education = "HS"
2977income_calib = Cagetti_income[education]
2979# Income specification
2980income_params = parse_income_spec(
2981 age_min=birth_age,
2982 age_max=death_age,
2983 adjust_infl_to=adjust_infl_to,
2984 **income_calib,
2985 SabelhausSong=True,
2986)
2988# Initial distribution of wealth and permanent income
2989dist_params = income_wealth_dists_from_scf(
2990 base_year=adjust_infl_to, age=birth_age, education=education, wave=1995
2991)
2993# We need survival probabilities only up to death_age-1, because survival
2994# probability at death_age is 1.
2995liv_prb = parse_ssa_life_table(
2996 female=False, cross_sec=True, year=2004, age_min=birth_age, age_max=death_age
2997)
2999# Parameters related to the number of periods implied by the calibration
3000time_params = parse_time_params(age_birth=birth_age, age_death=death_age)
3002# Update all the new parameters
3003init_lifecycle = copy(init_idiosyncratic_shocks)
3004del init_lifecycle["constructors"]
3005init_lifecycle.update(time_params)
3006init_lifecycle.update(dist_params)
3007# Note the income specification overrides the pLvlInitMean from the SCF.
3008init_lifecycle.update(income_params)
3009init_lifecycle.update({"LivPrb": liv_prb})
3010init_lifecycle["Rfree"] = init_lifecycle["T_cycle"] * init_lifecycle["Rfree"]
3012# Make a dictionary to specify an infinite consumer with a four period cycle
3013init_cyclical = copy(init_idiosyncratic_shocks)
3014init_cyclical["PermGroFac"] = [1.1, 1.082251, 2.8, 0.3]
3015init_cyclical["PermShkStd"] = [0.1, 0.1, 0.1, 0.1]
3016init_cyclical["TranShkStd"] = [0.1, 0.1, 0.1, 0.1]
3017init_cyclical["LivPrb"] = 4 * [0.98]
3018init_cyclical["Rfree"] = 4 * [1.03]
3019init_cyclical["T_cycle"] = 4