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