Coverage for HARK / ConsumptionSaving / ConsIndShockModel.py: 99%
892 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-08 05:31 +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 when True
1957 # (Uses linear spline interpolation for cFunc when False)
1958}
1959IndShockConsumerType_simulation_default = {
1960 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1961 "AgentCount": 10000, # Number of agents of this type
1962 "T_age": None, # Age after which simulated agents are automatically killed
1963 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1964 # (The portion of PermGroFac attributable to aggregate productivity growth)
1965 "NewbornTransShk": False, # Whether Newborns have transitory shock
1966 # ADDITIONAL OPTIONAL PARAMETERS
1967 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1968 # (Forces Newborns to follow solution path of the agent they replaced if True)
1969 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1970}
1972IndShockConsumerType_defaults = {}
1973IndShockConsumerType_defaults.update(IndShockConsumerType_IncShkDstn_default)
1974IndShockConsumerType_defaults.update(IndShockConsumerType_kNrmInitDstn_default)
1975IndShockConsumerType_defaults.update(IndShockConsumerType_pLvlInitDstn_default)
1976IndShockConsumerType_defaults.update(IndShockConsumerType_aXtraGrid_default)
1977IndShockConsumerType_defaults.update(IndShockConsumerType_solving_default)
1978IndShockConsumerType_defaults.update(IndShockConsumerType_simulation_default)
1979init_idiosyncratic_shocks = IndShockConsumerType_defaults # Here so that other models which use the old convention don't break
1982class IndShockConsumerType(PerfForesightConsumerType):
1983 r"""
1984 A consumer type with idiosyncratic shocks to permanent and transitory income.
1985 Their problem is defined by a sequence of income distributions, survival probabilities
1986 (:math:`\mathsf{S}`), and permanent income growth rates (:math:`\Gamma`), as well
1987 as time invariant values for risk aversion (:math:`\rho`), discount factor (:math:`\beta`),
1988 the interest rate (:math:`\mathsf{R}`), the grid of end-of-period assets, and an artificial
1989 borrowing constraint (:math:`\underline{a}`).
1991 .. math::
1992 \newcommand{\CRRA}{\rho}
1993 \newcommand{\LivPrb}{\mathsf{S}}
1994 \newcommand{\PermGroFac}{\Gamma}
1995 \newcommand{\Rfree}{\mathsf{R}}
1996 \newcommand{\DiscFac}{\beta}
1997 \begin{align*}
1998 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], \\
1999 & \text{s.t.} \\
2000 a_t &= m_t - c_t, \\
2001 a_t &\geq \underline{a}, \\
2002 m_{t+1} &= a_t \Rfree_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
2003 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
2004 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1, \\
2005 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}.
2006 \end{align*}
2009 Constructors
2010 ------------
2011 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2012 The agent's income shock distributions.
2014 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2015 aXtraGrid: Constructor
2016 The agent's asset grid.
2018 Its default constructor is :func:`HARK.utilities.make_assets_grid`
2020 Solving Parameters
2021 ------------------
2022 cycles: int
2023 0 specifies an infinite horizon model, 1 specifies a finite model.
2024 T_cycle: int
2025 Number of periods in the cycle for this agent type.
2026 CRRA: float, :math:`\rho`
2027 Coefficient of Relative Risk Aversion.
2028 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
2029 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
2030 DiscFac: float, :math:`\beta`
2031 Intertemporal discount factor.
2032 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2033 Survival probability after each period.
2034 PermGroFac: list[float], time varying, :math:`\Gamma`
2035 Permanent income growth factor.
2036 BoroCnstArt: float, :math:`\underline{a}`
2037 The minimum Asset/Perminant Income ratio, None to ignore.
2038 vFuncBool: bool
2039 Whether to calculate the value function during solution.
2040 CubicBool: bool
2041 Whether to use cubic spline interpoliation.
2043 Simulation Parameters
2044 ---------------------
2045 AgentCount: int
2046 Number of agents of this kind that are created during simulations.
2047 T_age: int
2048 Age after which to automatically kill agents, None to ignore.
2049 T_sim: int, required for simulation
2050 Number of periods to simulate.
2051 track_vars: list[strings]
2052 List of variables that should be tracked when running the simulation.
2053 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2055 PermShk is the agent's permanent income shock
2057 TranShk is the agent's transitory income shock
2059 aLvl is the nominal asset level
2061 aNrm is the normalized assets
2063 bNrm is the normalized resources without this period's labor income
2065 cNrm is the normalized consumption
2067 mNrm is the normalized market resources
2069 pLvl is the permanent income level
2071 who_dies is the array of which agents died
2072 kLogInitMean: float
2073 Mean of Log initial Normalized Assets.
2074 kLogInitStd: float
2075 Std of Log initial Normalized Assets.
2076 pLogInitMean: float
2077 Mean of Log initial permanent income.
2078 pLogInitStd: float
2079 Std of Log initial permanent income.
2080 PermGroFacAgg: float
2081 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2082 PerfMITShk: boolean
2083 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2084 NewbornTransShk: boolean
2085 Whether Newborns have transitory shock.
2087 Attributes
2088 ----------
2089 solution: list[Consumer solution object]
2090 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2091 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2093 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2094 history: Dict[Array]
2095 Created by running the :func:`.simulate()` method.
2096 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2097 Visit :class:`HARK.core.AgentType.simulate` for more information.
2098 """
2100 IncShkDstn_defaults = IndShockConsumerType_IncShkDstn_default
2101 aXtraGrid_defaults = IndShockConsumerType_aXtraGrid_default
2102 solving_defaults = IndShockConsumerType_solving_default
2103 simulation_defaults = IndShockConsumerType_simulation_default
2104 default_ = {
2105 "params": IndShockConsumerType_defaults,
2106 "solver": solve_one_period_ConsIndShock,
2107 "model": "ConsIndShock.yaml",
2108 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
2109 }
2111 time_inv_ = PerfForesightConsumerType.time_inv_ + [
2112 "vFuncBool",
2113 "CubicBool",
2114 "aXtraGrid",
2115 ]
2116 time_vary_ = PerfForesightConsumerType.time_vary_ + [
2117 "IncShkDstn",
2118 "PermShkDstn",
2119 "TranShkDstn",
2120 ]
2121 # This is in the PerfForesight model but not ConsIndShock
2122 time_inv_.remove("MaxKinks")
2123 shock_vars_ = ["PermShk", "TranShk"]
2124 distributions = [
2125 "IncShkDstn",
2126 "PermShkDstn",
2127 "TranShkDstn",
2128 "kNrmInitDstn",
2129 "pLvlInitDstn",
2130 ]
2132 def update_income_process(self):
2133 self.update("IncShkDstn", "PermShkDstn", "TranShkDstn")
2135 def get_shocks(self):
2136 """
2137 Gets permanent and transitory income shocks for this period. Samples from IncShkDstn for
2138 each period in the cycle.
2140 Parameters
2141 ----------
2142 NewbornTransShk : boolean, optional
2143 Whether Newborns have transitory shock. The default is False.
2145 Returns
2146 -------
2147 None
2148 """
2149 # Whether Newborns have transitory shock. The default is False.
2150 NewbornTransShk = self.NewbornTransShk
2152 PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
2153 TranShkNow = np.zeros(self.AgentCount)
2154 newborn = self.t_age == 0
2155 for t in np.unique(self.t_cycle):
2156 idx = self.t_cycle == t
2158 # temporary, see #1022
2159 if self.cycles == 1:
2160 t = t - 1
2162 N = np.sum(idx)
2163 if N > 0:
2164 # set current income distribution
2165 IncShkDstnNow = self.IncShkDstn[t]
2166 # and permanent growth factor
2167 PermGroFacNow = self.PermGroFac[t]
2168 # Get random draws of income shocks from the discrete distribution
2169 IncShks = IncShkDstnNow.draw(N)
2171 PermShkNow[idx] = (
2172 IncShks[0, :] * PermGroFacNow
2173 ) # permanent "shock" includes expected growth
2174 TranShkNow[idx] = IncShks[1, :]
2176 # That procedure used the *last* period in the sequence for newborns, but that's not right
2177 # Redraw shocks for newborns, using the *first* period in the sequence. Approximation.
2178 N = np.sum(newborn)
2179 if N > 0:
2180 idx = newborn
2181 # set current income distribution
2182 IncShkDstnNow = self.IncShkDstn[0]
2183 PermGroFacNow = self.PermGroFac[0] # and permanent growth factor
2185 # Get random draws of income shocks from the discrete distribution
2186 EventDraws = IncShkDstnNow.draw_events(N)
2187 PermShkNow[idx] = (
2188 IncShkDstnNow.atoms[0][EventDraws] * PermGroFacNow
2189 ) # permanent "shock" includes expected growth
2190 TranShkNow[idx] = IncShkDstnNow.atoms[1][EventDraws]
2192 # Whether Newborns have transitory shock. The default is False.
2193 if not NewbornTransShk:
2194 TranShkNow[newborn] = 1.0
2196 # Store the shocks in self
2197 self.shocks["PermShk"] = PermShkNow
2198 self.shocks["TranShk"] = TranShkNow
2200 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
2201 """
2202 Creates a "normalized Euler error" function for this instance, mapping
2203 from market resources to "consumption error per dollar of consumption."
2204 Stores result in attribute eulerErrorFunc as an interpolated function.
2205 Has option to use approximate income distribution stored in self.IncShkDstn
2206 or to use a (temporary) very dense approximation.
2208 Only works on (one period) infinite horizon models at this time, will
2209 be generalized later.
2211 Parameters
2212 ----------
2213 mMax : float
2214 Maximum normalized market resources for the Euler error function.
2215 approx_inc_dstn : Boolean
2216 Indicator for whether to use the approximate discrete income distri-
2217 bution stored in self.IncShkDstn[0], or to use a very accurate
2218 discrete approximation instead. When True, uses approximation in
2219 IncShkDstn; when False, makes and uses a very dense approximation.
2221 Returns
2222 -------
2223 None
2225 Notes
2226 -----
2227 This method is not used by any other code in the library. Rather, it is here
2228 for expository and benchmarking purposes.
2229 """
2230 # Get the income distribution (or make a very dense one)
2231 if approx_inc_dstn:
2232 IncShkDstn = self.IncShkDstn[0]
2233 else:
2234 TranShkDstn = MeanOneLogNormal(sigma=self.TranShkStd[0]).discretize(
2235 N=200,
2236 method="equiprobable",
2237 tail_N=50,
2238 tail_order=1.3,
2239 tail_bound=[0.05, 0.95],
2240 )
2241 TranShkDstn = add_discrete_outcome_constant_mean(
2242 TranShkDstn, p=self.UnempPrb, x=self.IncUnemp
2243 )
2244 PermShkDstn = MeanOneLogNormal(sigma=self.PermShkStd[0]).discretize(
2245 N=200,
2246 method="equiprobable",
2247 tail_N=50,
2248 tail_order=1.3,
2249 tail_bound=[0.05, 0.95],
2250 )
2251 IncShkDstn = combine_indep_dstns(PermShkDstn, TranShkDstn)
2253 # Make a grid of market resources
2254 mNowMin = self.solution[0].mNrmMin + 10 ** (
2255 -15
2256 ) # add tiny bit to get around 0/0 problem
2257 mNowMax = mMax
2258 mNowGrid = np.linspace(mNowMin, mNowMax, 1000)
2260 # Get the consumption function this period and the marginal value function
2261 # for next period. Note that this part assumes a one period cycle.
2262 cFuncNow = self.solution[0].cFunc
2263 vPfuncNext = self.solution[0].vPfunc
2265 # Calculate consumption this period at each gridpoint (and assets)
2266 cNowGrid = cFuncNow(mNowGrid)
2267 aNowGrid = mNowGrid - cNowGrid
2269 # Tile the grids for fast computation
2270 ShkCount = IncShkDstn.pmv.size
2271 aCount = aNowGrid.size
2272 aNowGrid_tiled = np.tile(aNowGrid, (ShkCount, 1))
2273 PermShkVals_tiled = (np.tile(IncShkDstn.atoms[0], (aCount, 1))).transpose()
2274 TranShkVals_tiled = (np.tile(IncShkDstn.atoms[1], (aCount, 1))).transpose()
2275 ShkPrbs_tiled = (np.tile(IncShkDstn.pmv, (aCount, 1))).transpose()
2277 # Calculate marginal value next period for each gridpoint and each shock
2278 mNextArray = (
2279 self.Rfree[0] / (self.PermGroFac[0] * PermShkVals_tiled) * aNowGrid_tiled
2280 + TranShkVals_tiled
2281 )
2282 vPnextArray = vPfuncNext(mNextArray)
2284 # Calculate expected marginal value and implied optimal consumption
2285 ExvPnextGrid = (
2286 self.DiscFac
2287 * self.Rfree[0]
2288 * self.LivPrb[0]
2289 * self.PermGroFac[0] ** (-self.CRRA)
2290 * np.sum(
2291 PermShkVals_tiled ** (-self.CRRA) * vPnextArray * ShkPrbs_tiled, axis=0
2292 )
2293 )
2294 cOptGrid = ExvPnextGrid ** (
2295 -1.0 / self.CRRA
2296 ) # This is the 'Endogenous Gridpoints' step
2298 # Calculate Euler error and store an interpolated function
2299 EulerErrorNrmGrid = (cNowGrid - cOptGrid) / cOptGrid
2300 eulerErrorFunc = LinearInterp(mNowGrid, EulerErrorNrmGrid)
2301 self.eulerErrorFunc = eulerErrorFunc
2303 def pre_solve(self):
2304 self.check_restrictions()
2305 self.construct("solution_terminal")
2306 if not self.quiet:
2307 self.check_conditions(verbose=self.verbose)
2309 def describe_parameters(self):
2310 """
2311 Generate a string describing the primitive model parameters that will
2312 be used to calculating limiting values and factors.
2314 Parameters
2315 ----------
2316 None
2318 Returns
2319 -------
2320 param_desc : str
2321 Description of primitive parameters.
2322 """
2323 # Get parameter description from the perfect foresight model
2324 param_desc = super().describe_parameters()
2326 # Make a new entry for weierstrass-p (the weird formatting here is to
2327 # make it easier to adapt into the style of the superclass if we add more
2328 # parameter reports later)
2329 this_entry = [
2330 "WorstPrb",
2331 "probability of worst income shock realization",
2332 "℘",
2333 False,
2334 ]
2335 try:
2336 val = getattr(self, this_entry[0])
2337 except:
2338 val = self.bilt[this_entry[0]]
2339 this_line = (
2340 this_entry[2]
2341 + f"={val:.5f} : "
2342 + this_entry[1]
2343 + " ("
2344 + this_entry[0]
2345 + ")\n"
2346 )
2348 # Add in the new entry and return it
2349 param_desc += this_line
2350 return param_desc
2352 def calc_limiting_values(self):
2353 """
2354 Compute various scalar values that are relevant to characterizing the
2355 solution to an infinite horizon problem. This method should only be called
2356 when T_cycle=1 and cycles=0, otherwise the values generated are meaningless.
2357 This method adds the following values to this instance in the dictionary
2358 attribute called bilt.
2360 APFac : Absolute Patience Factor
2361 GPFacRaw : Growth Patience Factor
2362 GPFacMod : Risk-Modified Growth Patience Factor
2363 GPFacLiv : Mortality-Adjusted Growth Patience Factor
2364 GPFacLivMod : Modigliani Mortality-Adjusted Growth Patience Factor
2365 GPFacSdl : Szeidl Growth Patience Factor
2366 FHWFac : Finite Human Wealth Factor
2367 RPFac : Return Patience Factor
2368 WRPFac : Weak Return Patience Factor
2369 PFVAFac : Perfect Foresight Value of Autarky Factor
2370 VAFac : Value of Autarky Factor
2371 cNrmPDV : Present Discounted Value of Autarky Consumption
2372 MPCmin : Limiting minimum MPC as market resources go to infinity
2373 MPCmax : Limiting maximum MPC as market resources approach minimum level
2374 hNrm : Human wealth divided by permanent income.
2375 ELogPermShk : Expected log permanent income shock
2376 WorstPrb : Probability of worst income shock realization
2377 Delta_mNrm_ZeroFunc : Linear locus where expected change in market resource ratio is zero
2378 BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income
2380 Returns
2381 -------
2382 None
2383 """
2384 super().calc_limiting_values()
2385 aux_dict = self.bilt
2387 # Calculate the risk-modified growth impatience factor
2388 PermShkDstn = self.PermShkDstn[0]
2389 inv_func = lambda x: x ** (-1.0)
2390 Ex_PermShkInv = expected(inv_func, PermShkDstn)[0]
2391 GroCompPermShk = Ex_PermShkInv ** (-1.0)
2392 aux_dict["GPFacMod"] = aux_dict["APFac"] / (self.PermGroFac[0] * GroCompPermShk)
2394 # Calculate the mortality-adjusted growth impatience factor (and version
2395 # with Modigiliani bequests)
2396 aux_dict["GPFacLiv"] = aux_dict["GPFacRaw"] * self.LivPrb[0]
2397 aux_dict["GPFacLivMod"] = aux_dict["GPFacLiv"] * self.LivPrb[0]
2399 # Calculate the risk-modified value of autarky factor
2400 if self.CRRA == 1.0:
2401 UtilCompPermShk = np.exp(expected(np.log, PermShkDstn)[0])
2402 else:
2403 CRRAfunc = lambda x: x ** (1.0 - self.CRRA)
2404 UtilCompPermShk = expected(CRRAfunc, PermShkDstn)[0] ** (
2405 1 / (1.0 - self.CRRA)
2406 )
2407 aux_dict["VAFac"] = self.DiscFac * (self.PermGroFac[0] * UtilCompPermShk) ** (
2408 1.0 - self.CRRA
2409 )
2411 # Calculate the expected log permanent income shock, which will be used
2412 # for the Szeidl variation of the Growth Impatience condition
2413 aux_dict["ELogPermShk"] = expected(np.log, PermShkDstn)[0]
2415 # Calculate the Harmenberg permanent income neutral expected log permanent
2416 # shock and the Harmenberg Growth Patience Factor
2417 Hrm_func = lambda x: x * np.log(x)
2418 PermShk_Hrm = np.exp(expected(Hrm_func, PermShkDstn)[0])
2419 aux_dict["GPFacHrm"] = aux_dict["GPFacRaw"] / PermShk_Hrm
2421 # Calculate the probability of the worst income shock realization
2422 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2423 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2424 ShkPrbsNext = self.IncShkDstn[0].pmv
2425 Ex_IncNext = np.dot(ShkPrbsNext, PermShkValsNext * TranShkValsNext)
2426 PermShkMinNext = np.min(PermShkValsNext)
2427 TranShkMinNext = np.min(TranShkValsNext)
2428 WorstIncNext = PermShkMinNext * TranShkMinNext
2429 WorstIncPrb = np.sum(
2430 ShkPrbsNext[(PermShkValsNext * TranShkValsNext) == WorstIncNext]
2431 )
2432 aux_dict["WorstPrb"] = WorstIncPrb
2434 # Calculate the weak return patience factor
2435 aux_dict["WRPFac"] = WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2437 # Calculate human wealth and the infinite horizon natural borrowing constraint
2438 if aux_dict["FHWFac"] < 1.0:
2439 hNrm = Ex_IncNext / (1.0 - aux_dict["FHWFac"])
2440 else:
2441 hNrm = np.inf
2442 temp = PermShkMinNext * aux_dict["FHWFac"]
2443 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2445 # Find the upper bound of the MPC as market resources approach the minimum
2446 BoroCnstArt = -np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2447 if BoroCnstNat < BoroCnstArt:
2448 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2449 else:
2450 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * aux_dict["RPFac"]
2451 MPCmax = np.maximum(MPCmax, 0.0)
2453 # Store maximum MPC and human wealth
2454 aux_dict["hNrm"] = hNrm
2455 aux_dict["MPCmax"] = MPCmax
2457 # Generate the "Delta m = 0" function, which is used to find target market resources
2458 # This overwrites the function generated by the perfect foresight version
2459 Ex_Rnrm = self.Rfree[0] / self.PermGroFac[0] * Ex_PermShkInv
2460 aux_dict["Delta_mNrm_ZeroFunc"] = (
2461 lambda m: (1.0 - 1.0 / Ex_Rnrm) * m + 1.0 / Ex_Rnrm
2462 )
2464 self.bilt = aux_dict
2466 self.bilt = aux_dict
2468 def check_GICMod(self, verbose=None):
2469 """
2470 Evaluate and report on the Risk-Modified Growth Impatience Condition.
2471 """
2472 name = "GICMod"
2473 GPFacMod = self.bilt["GPFacMod"]
2474 result = GPFacMod < 1.0
2476 messages = {
2477 True: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor satisfies the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2478 False: f"GPFacMod={GPFacMod:.5f} : The Risk-Modified Growth Patience Factor violates the Risk-Modified Growth Impatience Condition (GICMod) Þ/(G‖Ψ‖_(-1)) < 1.",
2479 }
2480 verbose = self.verbose if verbose is None else verbose
2481 self.log_condition_result(name, result, messages[result], verbose)
2483 def check_GICSdl(self, verbose=None):
2484 """
2485 Evaluate and report on the Szeidl variation of the Growth Impatience Condition.
2486 """
2487 name = "GICSdl"
2488 ELogPermShk = self.bilt["ELogPermShk"]
2489 result = np.log(self.bilt["GPFacRaw"]) < ELogPermShk
2491 messages = {
2492 True: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock satisfies the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2493 False: f"E[log Ψ]={ELogPermShk:.5f} : The expected log permanent income shock violates the Szeidl Growth Impatience Condition (GICSdl) log(Þ/G) < E[log Ψ].",
2494 }
2495 verbose = self.verbose if verbose is None else verbose
2496 self.log_condition_result(name, result, messages[result], verbose)
2498 def check_GICHrm(self, verbose=None):
2499 """
2500 Evaluate and report on the Harmenberg variation of the Growth Impatience Condition.
2501 """
2502 name = "GICHrm"
2503 GPFacHrm = self.bilt["GPFacHrm"]
2504 result = GPFacHrm < 1.0
2506 messages = {
2507 True: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor satisfies the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2508 False: f"GPFacHrm={GPFacHrm:.5f} : The Harmenberg Expected Growth Patience Factor violates the Harmenberg Growth Normalized Impatience Condition (GICHrm) Þ/G < exp(E[Ψlog Ψ]).",
2509 }
2510 verbose = self.verbose if verbose is None else verbose
2511 self.log_condition_result(name, result, messages[result], verbose)
2513 def check_GICLiv(self, verbose=None):
2514 """
2515 Evaluate and report on the Mortality-Adjusted Growth Impatience Condition.
2516 """
2517 name = "GICLiv"
2518 GPFacLiv = self.bilt["GPFacLiv"]
2519 result = GPFacLiv < 1.0
2521 messages = {
2522 True: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor satisfies the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2523 False: f"GPFacLiv={GPFacLiv:.5f} : The Mortality-Adjusted Growth Patience Factor violates the Mortality-Adjusted Growth Impatience Condition (GICLiv) ℒÞ/G < 1.",
2524 }
2525 verbose = self.verbose if verbose is None else verbose
2526 self.log_condition_result(name, result, messages[result], verbose)
2528 def check_FVAC(self, verbose=None):
2529 """
2530 Evaluate and report on the Finite Value of Autarky condition in the presence of income risk.
2531 """
2532 name = "FVAC"
2533 VAFac = self.bilt["VAFac"]
2534 result = VAFac < 1.0
2536 messages = {
2537 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.",
2538 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.",
2539 }
2540 verbose = self.verbose if verbose is None else verbose
2541 self.log_condition_result(name, result, messages[result], verbose)
2543 def check_WRIC(self, verbose=None):
2544 """
2545 Evaluate and report on the Weak Return Impatience Condition.
2546 """
2547 name = "WRIC"
2548 WRPFac = self.bilt["WRPFac"]
2549 result = WRPFac < 1.0
2551 messages = {
2552 True: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor satisfies the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2553 False: f"WRPFac={WRPFac:.5f} : The Weak Return Patience Factor violates the Weak Return Impatience Condition (WRIC) ℘ Þ/R < 1.",
2554 }
2555 verbose = self.verbose if verbose is None else verbose
2556 self.log_condition_result(name, result, messages[result], verbose)
2558 def check_conditions(self, verbose=None):
2559 """
2560 This method checks whether the instance's type satisfies various conditions.
2561 When combinations of these conditions are satisfied, the solution to the
2562 problem exhibits different characteristics. (For an exposition of the
2563 conditions, see https://econ-ark.github.io/BufferStockTheory/)
2565 Parameters
2566 ----------
2567 verbose : boolean
2568 Specifies different levels of verbosity of feedback. When False, it only reports whether the
2569 instance's type fails to satisfy a particular condition. When True, it reports all results, i.e.
2570 the factor values for all conditions.
2572 Returns
2573 -------
2574 None
2575 """
2576 self.conditions = {}
2577 self.bilt["conditions_report"] = ""
2578 self.degenerate = False
2579 verbose = self.verbose if verbose is None else verbose
2581 # This method only checks for the conditions for infinite horizon models
2582 # with a 1 period cycle. If these conditions are not met, we exit early.
2583 if self.cycles != 0 or self.T_cycle > 1:
2584 trivial_message = "No conditions report was produced because this functionality is only supported for infinite horizon models with a cycle length of 1."
2585 self.log_condition_result(None, None, trivial_message, verbose)
2586 if not self.quiet:
2587 _log.info(self.bilt["conditions_report"])
2588 return
2590 # Calculate some useful quantities that will be used in the condition checks
2591 self.calc_limiting_values()
2592 param_desc = self.describe_parameters()
2593 self.log_condition_result(None, None, param_desc, verbose)
2595 # Check individual conditions and add their results to the report
2596 self.check_AIC(verbose)
2597 self.check_RIC(verbose)
2598 self.check_WRIC(verbose)
2599 self.check_GICRaw(verbose)
2600 self.check_GICMod(verbose)
2601 self.check_GICLiv(verbose)
2602 self.check_GICSdl(verbose)
2603 self.check_GICHrm(verbose)
2604 super().check_FVAC(verbose)
2605 self.check_FVAC(verbose)
2606 self.check_FHWC(verbose)
2608 # Exit now if verbose output was not requested.
2609 if not verbose:
2610 if not self.quiet:
2611 _log.info(self.bilt["conditions_report"])
2612 return
2614 # Report on the degeneracy of the consumption function solution
2615 if self.conditions["WRIC"] and self.conditions["FVAC"]:
2616 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."
2617 degenerate = False
2618 elif not self.conditions["WRIC"]:
2619 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"
2620 degenerate = True
2621 elif not self.conditions["FVAC"]:
2622 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."
2623 degenerate = False
2624 self.log_condition_result(None, None, degen_message, verbose)
2625 self.degenerate = degenerate
2627 # Stop here if the solution is degenerate
2628 if degenerate:
2629 if not self.quiet:
2630 _log.info(self.bilt["conditions_report"])
2631 return
2633 # Report on the limiting behavior of the consumption function as m goes to infinity
2634 if self.conditions["RIC"]:
2635 if self.conditions["FHWC"]:
2636 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."
2637 else:
2638 RIC_message = "\nBecause the RIC is satisfied but the FHWC is violated, the GIC is satisfied."
2639 else:
2640 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."
2641 self.log_condition_result(None, None, RIC_message, verbose)
2643 # Report on whether a pseudo-steady-state exists at the individual level
2644 if self.conditions["GICRaw"]:
2645 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."
2646 else:
2647 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."
2648 self.log_condition_result(None, None, GIC_message, verbose)
2650 # Report on whether a target wealth ratio exists at the individual level
2651 if self.conditions["GICMod"]:
2652 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."
2653 else:
2654 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."
2655 self.log_condition_result(None, None, GICMod_message, verbose)
2657 # Report on whether a target level of wealth exists at the aggregate level
2658 if self.conditions["GICLiv"]:
2659 GICLiv_message = "\nBecause the GICLiv is satisfied, a target ratio of aggregate market resources to aggregate permanent income exists."
2660 else:
2661 GICLiv_message = "\nBecause the GICLiv is violated, a target ratio of aggregate market resources to aggregate permanent income might not exist."
2662 self.log_condition_result(None, None, GICLiv_message, verbose)
2664 # Report on whether invariant distributions exist
2665 if self.conditions["GICSdl"]:
2666 GICSdl_message = "\nBecause the GICSdl is satisfied, there exist invariant distributions of permanent income-normalized variables."
2667 else:
2668 GICSdl_message = "\nBecause the GICSdl is violated, there do not exist invariant distributions of permanent income-normalized variables."
2669 self.log_condition_result(None, None, GICSdl_message, verbose)
2671 # Report on whether blah blah
2672 if self.conditions["GICHrm"]:
2673 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."
2674 else:
2675 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.."
2676 self.log_condition_result(None, None, GICHrm_message, verbose)
2678 if not self.quiet:
2679 _log.info(self.bilt["conditions_report"])
2682###############################################################################
2684# Specify default parameters used in "kinked R" model
2686KinkedRconsumerType_IncShkDstn_default = IndShockConsumerType_IncShkDstn_default.copy()
2687KinkedRconsumerType_aXtraGrid_default = IndShockConsumerType_aXtraGrid_default.copy()
2688KinkedRconsumerType_kNrmInitDstn_default = (
2689 IndShockConsumerType_kNrmInitDstn_default.copy()
2690)
2691KinkedRconsumerType_pLvlInitDstn_default = (
2692 IndShockConsumerType_pLvlInitDstn_default.copy()
2693)
2695KinkedRconsumerType_solving_default = IndShockConsumerType_solving_default.copy()
2696KinkedRconsumerType_solving_default.update(
2697 {
2698 "Rboro": 1.20, # Interest factor on assets when borrowing, a < 0
2699 "Rsave": 1.02, # Interest factor on assets when saving, a > 0
2700 "BoroCnstArt": None, # Kinked R only matters if borrowing is allowed
2701 }
2702)
2703del KinkedRconsumerType_solving_default["Rfree"]
2705KinkedRconsumerType_simulation_default = IndShockConsumerType_simulation_default.copy()
2707KinkedRconsumerType_defaults = {}
2708KinkedRconsumerType_defaults.update(
2709 KinkedRconsumerType_IncShkDstn_default
2710) # Fill with some parameters
2711KinkedRconsumerType_defaults.update(KinkedRconsumerType_pLvlInitDstn_default)
2712KinkedRconsumerType_defaults.update(KinkedRconsumerType_kNrmInitDstn_default)
2713KinkedRconsumerType_defaults.update(KinkedRconsumerType_aXtraGrid_default)
2714KinkedRconsumerType_defaults.update(KinkedRconsumerType_solving_default)
2715KinkedRconsumerType_defaults.update(KinkedRconsumerType_simulation_default)
2716init_kinked_R = KinkedRconsumerType_defaults
2719class KinkedRconsumerType(IndShockConsumerType):
2720 r"""
2721 A consumer type based on IndShockConsumerType, with different
2722 interest rates for saving (:math:`\mathsf{R}_{save}`) and borrowing
2723 (:math:`\mathsf{R}_{boro}`).
2725 .. math::
2726 \newcommand{\CRRA}{\rho}
2727 \newcommand{\DiePrb}{\mathsf{D}}
2728 \newcommand{\PermGroFac}{\Gamma}
2729 \newcommand{\Rfree}{\mathsf{R}}
2730 \newcommand{\DiscFac}{\beta}
2731 \begin{align*}
2732 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], \\
2733 & \text{s.t.} \\
2734 a_t &= m_t - c_t, \\
2735 a_t &\geq \underline{a}, \\
2736 m_{t+1} &= \Rfree_t/(\PermGroFac_{t+1} \psi_{t+1}) a_t + \theta_{t+1}, \\
2737 \Rfree_t &= \begin{cases}
2738 \Rfree_{boro} & \text{if } a_t < 0\\
2739 \Rfree_{save} & \text{if } a_t \geq 0,
2740 \end{cases}\\
2741 \Rfree_{boro} &> \Rfree_{save}, \\
2742 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
2743 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1.\\
2744 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA} \\
2745 \end{align*}
2747 Constructors
2748 ------------
2749 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
2750 The agent's income shock distributions.
2751 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
2752 aXtraGrid: Constructor
2753 The agent's asset grid.
2754 Its default constructor is :func:`HARK.utilities.make_assets_grid`
2756 Solving Parameters
2757 ------------------
2758 cycles: int
2759 0 specifies an infinite horizon model, 1 specifies a finite model.
2760 T_cycle: int
2761 Number of periods in the cycle for this agent type.
2762 CRRA: float, :math:`\rho`
2763 Coefficient of Relative Risk Aversion.
2764 Rboro: float, :math:`\mathsf{R}_{boro}`
2765 Risk Free interest rate when assets are negative.
2766 Rsave: float, :math:`\mathsf{R}_{save}`
2767 Risk Free interest rate when assets are positive.
2768 DiscFac: float, :math:`\beta`
2769 Intertemporal discount factor.
2770 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
2771 Survival probability after each period.
2772 PermGroFac: list[float], time varying, :math:`\Gamma`
2773 Permanent income growth factor.
2774 BoroCnstArt: float, :math:`\underline{a}`
2775 The minimum Asset/Perminant Income ratio, None to ignore.
2776 vFuncBool: bool
2777 Whether to calculate the value function during solution.
2778 CubicBool: bool
2779 Whether to use cubic spline interpoliation.
2781 Simulation Parameters
2782 ---------------------
2783 AgentCount: int
2784 Number of agents of this kind that are created during simulations.
2785 T_age: int
2786 Age after which to automatically kill agents, None to ignore.
2787 T_sim: int, required for simulation
2788 Number of periods to simulate.
2789 track_vars: list[strings]
2790 List of variables that should be tracked when running the simulation.
2791 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
2793 PermShk is the agent's permanent income shock
2795 TranShk is the agent's transitory income shock
2797 aLvl is the nominal asset level
2799 aNrm is the normalized assets
2801 bNrm is the normalized resources without this period's labor income
2803 cNrm is the normalized consumption
2805 mNrm is the normalized market resources
2807 pLvl is the permanent income level
2809 who_dies is the array of which agents died
2810 kLogInitMean: float
2811 Mean of Log initial Normalized Assets.
2812 kLogInitStd: float
2813 Std of Log initial Normalized Assets.
2814 pLogInitMean: float
2815 Mean of Log initial permanent income.
2816 pLogInitStd: float
2817 Std of Log initial permanent income.
2818 PermGroFacAgg: float
2819 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
2820 PerfMITShk: boolean
2821 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
2822 NewbornTransShk: boolean
2823 Whether Newborns have transitory shock.
2825 Attributes
2826 ----------
2827 solution: list[Consumer solution object]
2828 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
2829 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
2831 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
2832 history: Dict[Array]
2833 Created by running the :func:`.simulate()` method.
2834 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
2835 Visit :class:`HARK.core.AgentType.simulate` for more information.
2836 """
2838 IncShkDstn_defaults = KinkedRconsumerType_IncShkDstn_default
2839 aXtraGrid_defaults = KinkedRconsumerType_aXtraGrid_default
2840 solving_defaults = KinkedRconsumerType_solving_default
2841 simulation_defaults = KinkedRconsumerType_simulation_default
2842 default_ = {
2843 "params": KinkedRconsumerType_defaults,
2844 "solver": solve_one_period_ConsKinkedR,
2845 "model": "ConsKinkedR.yaml",
2846 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
2847 }
2849 time_inv_ = copy(IndShockConsumerType.time_inv_)
2850 time_inv_ += ["Rboro", "Rsave"]
2852 def calc_bounding_values(self):
2853 """
2854 Calculate human wealth plus minimum and maximum MPC in an infinite
2855 horizon model with only one period repeated indefinitely. Store results
2856 as attributes of self. Human wealth is the present discounted value of
2857 expected future income after receiving income this period, ignoring mort-
2858 ality. The maximum MPC is the limit of the MPC as m --> mNrmMin. The
2859 minimum MPC is the limit of the MPC as m --> infty. This version deals
2860 with the different interest rates on borrowing vs saving.
2862 Parameters
2863 ----------
2864 None
2866 Returns
2867 -------
2868 None
2869 """
2870 # Unpack the income distribution and get average and worst outcomes
2871 PermShkValsNext = self.IncShkDstn[0].atoms[0]
2872 TranShkValsNext = self.IncShkDstn[0].atoms[1]
2873 ShkPrbsNext = self.IncShkDstn[0].pmv
2874 IncNext = PermShkValsNext * TranShkValsNext
2875 Ex_IncNext = np.dot(ShkPrbsNext, IncNext)
2876 PermShkMinNext = np.min(PermShkValsNext)
2877 TranShkMinNext = np.min(TranShkValsNext)
2878 WorstIncNext = PermShkMinNext * TranShkMinNext
2879 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
2880 # TODO: Check the math above. I think it fails for non-independent shocks
2882 BoroCnstArt = np.inf if self.BoroCnstArt is None else self.BoroCnstArt
2884 # Calculate human wealth and the infinite horizon natural borrowing constraint
2885 hNrm = (Ex_IncNext * self.PermGroFac[0] / self.Rsave) / (
2886 1.0 - self.PermGroFac[0] / self.Rsave
2887 )
2888 temp = self.PermGroFac[0] * PermShkMinNext / self.Rboro
2889 BoroCnstNat = -TranShkMinNext * temp / (1.0 - temp)
2891 PatFacTop = (self.DiscFac * self.LivPrb[0] * self.Rsave) ** (
2892 1.0 / self.CRRA
2893 ) / self.Rsave
2894 PatFacBot = (self.DiscFac * self.LivPrb[0] * self.Rboro) ** (
2895 1.0 / self.CRRA
2896 ) / self.Rboro
2897 if BoroCnstNat < BoroCnstArt:
2898 MPCmax = 1.0 # if natural borrowing constraint is overridden by artificial one, MPCmax is 1
2899 else:
2900 MPCmax = 1.0 - WorstIncPrb ** (1.0 / self.CRRA) * PatFacBot
2901 MPCmin = 1.0 - PatFacTop
2903 # Store the results as attributes of self
2904 self.hNrm = hNrm
2905 self.MPCmin = MPCmin
2906 self.MPCmax = MPCmax
2908 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True): # pragma: nocover
2909 """
2910 Creates a "normalized Euler error" function for this instance, mapping
2911 from market resources to "consumption error per dollar of consumption."
2912 Stores result in attribute eulerErrorFunc as an interpolated function.
2913 Has option to use approximate income distribution stored in self.IncShkDstn
2914 or to use a (temporary) very dense approximation.
2916 SHOULD BE INHERITED FROM ConsIndShockModel
2918 Parameters
2919 ----------
2920 mMax : float
2921 Maximum normalized market resources for the Euler error function.
2922 approx_inc_dstn : Boolean
2923 Indicator for whether to use the approximate discrete income distri-
2924 bution stored in self.IncShkDstn[0], or to use a very accurate
2925 discrete approximation instead. When True, uses approximation in
2926 IncShkDstn; when False, makes and uses a very dense approximation.
2928 Returns
2929 -------
2930 None
2931 """
2932 raise NotImplementedError()
2934 def get_Rport(self):
2935 """
2936 Returns an array of size self.AgentCount with self.Rboro or self.Rsave in
2937 each entry, based on whether self.aNrmNow >< 0. This represents the risk-
2938 free portfolio return in this model.
2940 Parameters
2941 ----------
2942 None
2944 Returns
2945 -------
2946 RfreeNow : np.array
2947 Array of size self.AgentCount with risk free interest rate for each agent.
2948 """
2949 RfreeNow = self.Rboro * np.ones(self.AgentCount)
2950 RfreeNow[self.state_prev["aNrm"] > 0] = self.Rsave
2951 return RfreeNow
2953 def check_conditions(self, verbose):
2954 """
2955 This empty method overwrites the version inherited from its parent class,
2956 IndShockConsumerType. The condition checks are not appropriate when Rfree
2957 has multiple values.
2959 Parameters
2960 ----------
2961 None
2963 Returns
2964 -------
2965 None
2966 """
2967 pass
2970###############################################################################
2972# Make a dictionary to specify a lifecycle consumer with a finite horizon
2974# Main calibration characteristics
2975birth_age = 25
2976death_age = 90
2977adjust_infl_to = 1992
2978# Use income estimates from Cagetti (2003) for High-school graduates
2979education = "HS"
2980income_calib = Cagetti_income[education]
2982# Income specification
2983income_params = parse_income_spec(
2984 age_min=birth_age,
2985 age_max=death_age,
2986 adjust_infl_to=adjust_infl_to,
2987 **income_calib,
2988 SabelhausSong=True,
2989)
2991# Initial distribution of wealth and permanent income
2992dist_params = income_wealth_dists_from_scf(
2993 base_year=adjust_infl_to, age=birth_age, education=education, wave=1995
2994)
2996# We need survival probabilities only up to death_age-1, because survival
2997# probability at death_age is 1.
2998liv_prb = parse_ssa_life_table(
2999 female=False, cross_sec=True, year=2004, age_min=birth_age, age_max=death_age
3000)
3002# Parameters related to the number of periods implied by the calibration
3003time_params = parse_time_params(age_birth=birth_age, age_death=death_age)
3005# Update all the new parameters
3006init_lifecycle = copy(init_idiosyncratic_shocks)
3007del init_lifecycle["constructors"]
3008init_lifecycle.update(time_params)
3009init_lifecycle.update(dist_params)
3010# Note the income specification overrides the pLvlInitMean from the SCF.
3011init_lifecycle.update(income_params)
3012init_lifecycle.update({"LivPrb": liv_prb})
3013init_lifecycle["Rfree"] = init_lifecycle["T_cycle"] * init_lifecycle["Rfree"]
3015# Make a dictionary to specify an infinite consumer with a four period cycle
3016init_cyclical = copy(init_idiosyncratic_shocks)
3017init_cyclical["PermGroFac"] = [1.1, 1.082251, 2.8, 0.3]
3018init_cyclical["PermShkStd"] = [0.1, 0.1, 0.1, 0.1]
3019init_cyclical["TranShkStd"] = [0.1, 0.1, 0.1, 0.1]
3020init_cyclical["LivPrb"] = 4 * [0.98]
3021init_cyclical["Rfree"] = 4 * [1.03]
3022init_cyclical["T_cycle"] = 4
3024# Make dictionaries based on Carroll QJE (1997) lifecycle specifications
3025buffer_stock_lifecycle_base = {
3026 "CRRA": 2.0,
3027 "DiscFac": 0.96,
3028 "PermGroFacAgg": 1.02,
3029 "kLogInitMean": -1000.0,
3030 "kLogInitStd": 0.0,
3031 "pLogInitStd": 0.0,
3032 "Rfree": 49 * [1.00],
3033 "PermShkStd": 40 * [0.1] + 9 * [0.0],
3034 "TranShkStd": 40 * [0.1] + 9 * [0.0],
3035 "UnempPrb": 0.005,
3036 "IncUnemp": 0.0,
3037 "UnempPrbRet": 0.0005,
3038 "IncUnempRet": 0.0,
3039 "LivPrb": 49 * [1.0],
3040 "T_cycle": 49,
3041 "T_retire": 40,
3042 "T_age": 50,
3043 "AgentCount": 10000,
3044 "cycles": 1,
3045}
3047unskilled_update = {
3048 "pLogInitMean": np.log(1 / 1.03),
3049 "PermGroFac": 14 * [1.03] + 25 * [1.0] + [0.7] + 9 * [1.0],
3050}
3052operative_update = {
3053 "pLogInitMean": np.log(1 / 1.025),
3054 "PermGroFac": 24 * [1.025] + 15 * [1.01] + [0.7] + 9 * [1.0],
3055}
3057manager_update = {
3058 "pLogInitMean": np.log(1 / 1.03),
3059 "PermGroFac": 29 * [1.03] + 10 * [0.99] + [0.7] + 9 * [1.0],
3060}
3062# Carroll QJE (1997) life-cycle calibration for unskilled workers
3063buffer_stock_lifecycle_unskilled = copy(buffer_stock_lifecycle_base)
3064buffer_stock_lifecycle_unskilled.update(unskilled_update)
3066# Carroll QJE (1997) life-cycle calibration for operative workers
3067buffer_stock_lifecycle_operative = copy(buffer_stock_lifecycle_base)
3068buffer_stock_lifecycle_operative.update(operative_update)
3070# Carroll QJE (1997) life-cycle calibration for managerial workers
3071buffer_stock_lifecycle_manager = copy(buffer_stock_lifecycle_base)
3072buffer_stock_lifecycle_manager.update(manager_update)