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