Coverage for HARK / ConsumptionSaving / ConsBequestModel.py: 96%
360 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-07 05:16 +0000
1"""
2Classes to solve consumption-saving models with a bequest motive and
3idiosyncratic shocks to income and wealth. All models here assume
4separable CRRA utility of consumption and Stone-Geary utility of
5savings with geometric discounting of the continuation value and
6shocks to income that have transitory and/or permanent components.
8It currently solves two types of models:
9 1) A standard lifecycle model with a terminal and/or accidental bequest motive.
10 2) A portfolio choice model with a terminal and/or accidental bequest motive.
11"""
13import numpy as np
15from HARK import NullFunc
16from HARK.Calibration.Income.IncomeProcesses import (
17 construct_lognormal_income_process_unemployment,
18 get_PermShkDstn_from_IncShkDstn,
19 get_TranShkDstn_from_IncShkDstn,
20)
21from HARK.ConsumptionSaving.ConsIndShockModel import (
22 ConsumerSolution,
23 IndShockConsumerType,
24 make_basic_CRRA_solution_terminal,
25 make_lognormal_kNrm_init_dstn,
26 make_lognormal_pLvl_init_dstn,
27)
28from HARK.Calibration.Assets.AssetProcesses import (
29 make_lognormal_RiskyDstn,
30 combine_IncShkDstn_and_RiskyDstn,
31 calc_ShareLimit_for_CRRA,
32)
33from HARK.ConsumptionSaving.ConsPortfolioModel import (
34 PortfolioConsumerType,
35 PortfolioSolution,
36 make_portfolio_solution_terminal,
37 make_AdjustDstn,
38)
39from HARK.ConsumptionSaving.ConsRiskyAssetModel import make_simple_ShareGrid
40from HARK.distributions import expected
41from HARK.interpolation import (
42 BilinearInterp,
43 ConstantFunction,
44 CubicInterp,
45 IdentityFunction,
46 LinearInterp,
47 LinearInterpOnInterp1D,
48 LowerEnvelope,
49 MargMargValueFuncCRRA,
50 MargValueFuncCRRA,
51 ValueFuncCRRA,
52)
53from HARK.rewards import UtilityFuncCRRA, UtilityFuncStoneGeary
54from HARK.utilities import make_assets_grid
57def make_bequest_solution_terminal(
58 CRRA, BeqCRRATerm, BeqFacTerm, BeqShiftTerm, aXtraGrid
59):
60 """
61 Make the terminal period solution when there is a warm glow bequest motive with
62 Stone-Geary form utility. If there is no warm glow bequest motive (BeqFacTerm = 0),
63 then the terminal period solution is identical to ConsIndShock.
65 Parameters
66 ----------
67 CRRA : float
68 Coefficient on relative risk aversion over consumption.
69 BeqCRRATerm : float
70 Coefficient on relative risk aversion in the terminal warm glow bequest motive.
71 BeqFacTerm : float
72 Scaling factor for the terminal warm glow bequest motive.
73 BeqShiftTerm : float
74 Stone-Geary shifter term for the terminal warm glow bequest motive.
75 aXtraGrid : np.array
76 Set of assets-above-minimum to be used in the solution.
78 Returns
79 -------
80 solution_terminal : ConsumerSolution
81 Terminal period solution when there is a warm glow bequest.
82 """
83 if BeqFacTerm == 0.0: # No terminal bequest
84 solution_terminal = make_basic_CRRA_solution_terminal(CRRA)
85 return solution_terminal
87 utility = UtilityFuncCRRA(CRRA)
88 warm_glow = UtilityFuncStoneGeary(
89 BeqCRRATerm,
90 factor=BeqFacTerm,
91 shifter=BeqShiftTerm,
92 )
94 aNrmGrid = np.append(0.0, aXtraGrid) if BeqShiftTerm != 0.0 else aXtraGrid
95 cNrmGrid = utility.derinv(warm_glow.der(aNrmGrid))
96 vGrid = utility(cNrmGrid) + warm_glow(aNrmGrid)
97 cNrmGridW0 = np.append(0.0, cNrmGrid)
98 mNrmGridW0 = np.append(0.0, aNrmGrid + cNrmGrid)
99 vNvrsGridW0 = np.append(0.0, utility.inv(vGrid))
101 cFunc_term = LinearInterp(mNrmGridW0, cNrmGridW0)
102 vNvrsFunc_term = LinearInterp(mNrmGridW0, vNvrsGridW0)
103 vFunc_term = ValueFuncCRRA(vNvrsFunc_term, CRRA)
104 vPfunc_term = MargValueFuncCRRA(cFunc_term, CRRA)
105 vPPfunc_term = MargMargValueFuncCRRA(cFunc_term, CRRA)
107 solution_terminal = ConsumerSolution(
108 cFunc=cFunc_term,
109 vFunc=vFunc_term,
110 vPfunc=vPfunc_term,
111 vPPfunc=vPPfunc_term,
112 mNrmMin=0.0,
113 hNrm=0.0,
114 MPCmax=1.0,
115 MPCmin=1.0,
116 )
117 return solution_terminal
120def make_warmglow_portfolio_solution_terminal(
121 CRRA, BeqCRRATerm, BeqFacTerm, BeqShiftTerm, aXtraGrid
122):
123 """
124 Make the terminal period solution when there is a warm glow bequest motive with
125 Stone-Geary form utility and portfolio choice. If there is no warm glow bequest
126 motive (BeqFacTerm = 0), then the terminal period solution is identical to ConsPortfolio.
128 Parameters
129 ----------
130 CRRA : float
131 Coefficient on relative risk aversion over consumption.
132 BeqCRRATerm : float
133 Coefficient on relative risk aversion in the terminal warm glow bequest motive.
134 BeqFacTerm : float
135 Scaling factor for the terminal warm glow bequest motive.
136 BeqShiftTerm : float
137 Stone-Geary shifter term for the terminal warm glow bequest motive.
138 aXtraGrid : np.array
139 Set of assets-above-minimum to be used in the solution.
141 Returns
142 -------
143 solution_terminal : ConsumerSolution
144 Terminal period solution when there is a warm glow bequest and portfolio choice.
145 """
146 if BeqFacTerm == 0.0: # No terminal bequest
147 solution_terminal = make_portfolio_solution_terminal(CRRA)
148 return solution_terminal
150 # Solve the terminal period problem when there is no portfolio choice
151 solution_terminal_no_port = make_bequest_solution_terminal(
152 CRRA, BeqCRRATerm, BeqFacTerm, BeqShiftTerm, aXtraGrid
153 )
155 # Take consumption function from the no portfolio choice solution
156 cFuncAdj_terminal = solution_terminal_no_port.cFunc
157 cFuncFxd_terminal = lambda m, s: solution_terminal_no_port(m)
159 # Risky share is irrelevant-- no end-of-period assets; set to zero
160 ShareFuncAdj_terminal = ConstantFunction(0.0)
161 ShareFuncFxd_terminal = IdentityFunction(i_dim=1, n_dims=2)
163 # Value function is simply utility from consuming market resources
164 vFuncAdj_terminal = solution_terminal_no_port.vFunc
165 vFuncFxd_terminal = lambda m, s: solution_terminal_no_port.cFunc(m)
167 # Marginal value of market resources is marg utility at the consumption function
168 vPfuncAdj_terminal = solution_terminal_no_port.vPfunc
169 dvdmFuncFxd_terminal = lambda m, s: solution_terminal_no_port.vPfunc(m)
170 # No future, no marg value of Share
171 dvdsFuncFxd_terminal = ConstantFunction(0.0)
173 # Construct the terminal period solution
174 solution_terminal = PortfolioSolution(
175 cFuncAdj=cFuncAdj_terminal,
176 ShareFuncAdj=ShareFuncAdj_terminal,
177 vFuncAdj=vFuncAdj_terminal,
178 vPfuncAdj=vPfuncAdj_terminal,
179 cFuncFxd=cFuncFxd_terminal,
180 ShareFuncFxd=ShareFuncFxd_terminal,
181 vFuncFxd=vFuncFxd_terminal,
182 dvdmFuncFxd=dvdmFuncFxd_terminal,
183 dvdsFuncFxd=dvdsFuncFxd_terminal,
184 )
185 return solution_terminal
188###############################################################################
191def solve_one_period_ConsWarmBequest(
192 solution_next,
193 IncShkDstn,
194 LivPrb,
195 DiscFac,
196 CRRA,
197 Rfree,
198 PermGroFac,
199 BoroCnstArt,
200 aXtraGrid,
201 BeqCRRA,
202 BeqFac,
203 BeqShift,
204 CubicBool,
205 vFuncBool,
206):
207 """
208 Solves one period of a consumption-saving model with idiosyncratic shocks to
209 permanent and transitory income, with one risk free asset and CRRA utility.
210 The consumer also has a "warm glow" bequest motive in which they gain additional
211 utility based on their terminal wealth upon death.
213 Parameters
214 ----------
215 solution_next : ConsumerSolution
216 The solution to next period's one period problem.
217 IncShkDstn : distribution.Distribution
218 A discrete approximation to the income process between the period being
219 solved and the one immediately following (in solution_next).
220 LivPrb : float
221 Survival probability; likelihood of being alive at the beginning of
222 the succeeding period.
223 DiscFac : float
224 Intertemporal discount factor for future utility.
225 CRRA : float
226 Coefficient of relative risk aversion.
227 Rfree : float
228 Risk free interest factor on end-of-period assets.
229 PermGroFac : float
230 Expected permanent income growth factor at the end of this period.
231 BoroCnstArt: float or None
232 Borrowing constraint for the minimum allowable assets to end the
233 period with. If it is less than the natural borrowing constraint,
234 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
235 rowing constraint.
236 aXtraGrid : np.array
237 Array of "extra" end-of-period asset values-- assets above the
238 absolute minimum acceptable level.
239 BeqCRRA : float
240 Coefficient of relative risk aversion for warm glow bequest motive.
241 BeqFac : float
242 Multiplicative intensity factor for the warm glow bequest motive.
243 BeqShift : float
244 Stone-Geary shifter in the warm glow bequest motive.
245 CubicBool : bool
246 An indicator for whether the solver should use cubic or linear interpolation.
247 vFuncBool: boolean
248 An indicator for whether the value function should be computed and
249 included in the reported solution.
251 Returns
252 -------
253 solution_now : ConsumerSolution
254 Solution to this period's consumption-saving problem with income risk.
255 """
256 # Define the current period utility function and effective discount factor
257 uFunc = UtilityFuncCRRA(CRRA)
258 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
259 BeqFacEff = (1.0 - LivPrb) * BeqFac # "effective" bequest factor
260 warm_glow = UtilityFuncStoneGeary(BeqCRRA, BeqFacEff, BeqShift)
262 # Unpack next period's income shock distribution
263 ShkPrbsNext = IncShkDstn.pmv
264 PermShkValsNext = IncShkDstn.atoms[0]
265 TranShkValsNext = IncShkDstn.atoms[1]
266 PermShkMinNext = np.min(PermShkValsNext)
267 TranShkMinNext = np.min(TranShkValsNext)
269 # Calculate the probability that we get the worst possible income draw
270 IncNext = PermShkValsNext * TranShkValsNext
271 WorstIncNext = PermShkMinNext * TranShkMinNext
272 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
273 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
275 # Unpack next period's (marginal) value function
276 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
277 vPfuncNext = solution_next.vPfunc
278 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
280 # Update the bounding MPCs and PDV of human wealth:
281 PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
282 try:
283 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
284 except:
285 MPCminNow = 0.0
286 Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
287 hNrmNow = PermGroFac / Rfree * (Ex_IncNext + solution_next.hNrm)
288 temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
289 MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
290 cFuncLimitIntercept = MPCminNow * hNrmNow
291 cFuncLimitSlope = MPCminNow
293 # Calculate the minimum allowable value of money resources in this period
294 PermGroFacEffMin = (PermGroFac * PermShkMinNext) / Rfree
295 BoroCnstNat = (solution_next.mNrmMin - TranShkMinNext) * PermGroFacEffMin
296 BoroCnstNat = np.max([BoroCnstNat, -BeqShift])
298 # Set the minimum allowable (normalized) market resources based on the natural
299 # and artificial borrowing constraints
300 if BoroCnstArt is None:
301 mNrmMinNow = BoroCnstNat
302 else:
303 mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt])
305 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
306 # or artificial borrowing constraint actually binds
307 if BoroCnstNat < mNrmMinNow:
308 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
309 else:
310 MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above
312 # Define the borrowing-constrained consumption function
313 cFuncNowCnst = LinearInterp(
314 np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0])
315 )
317 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
318 aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
320 # Define local functions for taking future expectations
321 def calc_mNrmNext(S, a, R):
322 return R / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
324 def calc_vNext(S, a, R):
325 return (S["PermShk"] ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)) * vFuncNext(
326 calc_mNrmNext(S, a, R)
327 )
329 def calc_vPnext(S, a, R):
330 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a, R))
332 def calc_vPPnext(S, a, R):
333 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, a, R))
335 # Calculate end-of-period marginal value of assets at each gridpoint
336 vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA)
337 EndOfPrdvP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(aNrmNow, Rfree))
338 EndOfPrdvP += warm_glow.der(aNrmNow)
340 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
341 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
342 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
344 # Limiting consumption is zero as m approaches mNrmMin
345 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
346 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
348 # Construct the consumption function as a cubic or linear spline interpolation
349 if CubicBool:
350 # Calculate end-of-period marginal marginal value of assets at each gridpoint
351 vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0)
352 EndOfPrdvPP = vPPfacEff * expected(
353 calc_vPPnext, IncShkDstn, args=(aNrmNow, Rfree)
354 )
355 EndOfPrdvPP += warm_glow.der(aNrmNow, order=2)
356 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
357 MPC = dcda / (dcda + 1.0)
358 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
360 # Construct the unconstrained consumption function as a cubic interpolation
361 cFuncNowUnc = CubicInterp(
362 m_for_interpolation,
363 c_for_interpolation,
364 MPC_for_interpolation,
365 cFuncLimitIntercept,
366 cFuncLimitSlope,
367 )
368 else:
369 # Construct the unconstrained consumption function as a linear interpolation
370 cFuncNowUnc = LinearInterp(
371 m_for_interpolation,
372 c_for_interpolation,
373 cFuncLimitIntercept,
374 cFuncLimitSlope,
375 )
377 # Combine the constrained and unconstrained functions into the true consumption function.
378 # LowerEnvelope should only be used when BoroCnstArt is True
379 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
381 # Make the marginal value function
382 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
384 # Define this period's marginal marginal value function
385 if CubicBool:
386 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
387 else:
388 vPPfuncNow = NullFunc() # Dummy object
390 # Construct this period's value function if requested
391 if vFuncBool:
392 # Calculate end-of-period value, its derivative, and their pseudo-inverse
393 EndOfPrdv = DiscFacEff * expected(calc_vNext, IncShkDstn, args=(aNrmNow, Rfree))
394 EndOfPrdv += warm_glow(aNrmNow)
395 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
396 # value transformed through inverse utility
397 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
398 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
399 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
400 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
402 # Construct the end-of-period value function
403 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
404 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
405 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
407 # Compute expected value and marginal value on a grid of market resources
408 mNrm_temp = mNrmMinNow + aXtraGrid
409 cNrm_temp = cFuncNow(mNrm_temp)
410 aNrm_temp = mNrm_temp - cNrm_temp
411 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
412 vP_temp = uFunc.der(cNrm_temp)
414 # Construct the beginning-of-period value function
415 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
416 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
417 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
418 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
419 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
420 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
421 vNvrsFuncNow = CubicInterp(
422 mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
423 )
424 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
425 else:
426 vFuncNow = NullFunc() # Dummy object
428 # Create and return this period's solution
429 solution_now = ConsumerSolution(
430 cFunc=cFuncNow,
431 vFunc=vFuncNow,
432 vPfunc=vPfuncNow,
433 vPPfunc=vPPfuncNow,
434 mNrmMin=mNrmMinNow,
435 hNrm=hNrmNow,
436 MPCmin=MPCminNow,
437 MPCmax=MPCmaxEff,
438 )
439 return solution_now
442###############################################################################
445def solve_one_period_ConsPortfolioWarmGlow(
446 solution_next,
447 IncShkDstn,
448 RiskyDstn,
449 LivPrb,
450 DiscFac,
451 CRRA,
452 Rfree,
453 PermGroFac,
454 BoroCnstArt,
455 aXtraGrid,
456 ShareGrid,
457 AdjustPrb,
458 ShareLimit,
459 vFuncBool,
460 DiscreteShareBool,
461 BeqCRRA,
462 BeqFac,
463 BeqShift,
464):
465 """
466 Solve one period of a consumption-saving problem with portfolio allocation
467 between a riskless and risky asset. This function handles various sub-cases
468 or variations on the problem, including the possibility that the agent does
469 not necessarily get to update their portfolio share in every period, or that
470 they must choose a discrete rather than continuous risky share.
472 Parameters
473 ----------
474 solution_next : PortfolioSolution
475 Solution to next period's problem.
476 ShockDstn : Distribution
477 Joint distribution of permanent income shocks, transitory income shocks,
478 and risky returns. This is only used if the input IndepDstnBool is False,
479 indicating that income and return distributions can't be assumed to be
480 independent.
481 IncShkDstn : Distribution
482 Discrete distribution of permanent income shocks and transitory income
483 shocks. This is only used if the input IndepDstnBool is True, indicating
484 that income and return distributions are independent.
485 RiskyDstn : Distribution
486 Distribution of risky asset returns. This is only used if the input
487 IndepDstnBool is True, indicating that income and return distributions
488 are independent.
489 LivPrb : float
490 Survival probability; likelihood of being alive at the beginning of
491 the succeeding period.
492 DiscFac : float
493 Intertemporal discount factor for future utility.
494 CRRA : float
495 Coefficient of relative risk aversion.
496 Rfree : float
497 Risk free interest factor on end-of-period assets.
498 PermGroFac : float
499 Expected permanent income growth factor at the end of this period.
500 BoroCnstArt: float or None
501 Borrowing constraint for the minimum allowable assets to end the
502 period with. In this model, it is *required* to be zero.
503 aXtraGrid: np.array
504 Array of "extra" end-of-period asset values-- assets above the
505 absolute minimum acceptable level.
506 ShareGrid : np.array
507 Array of risky portfolio shares on which to define the interpolation
508 of the consumption function when Share is fixed. Also used when the
509 risky share choice is specified as discrete rather than continuous.
510 AdjustPrb : float
511 Probability that the agent will be able to update his portfolio share.
512 ShareLimit : float
513 Limiting lower bound of risky portfolio share as mNrm approaches infinity.
514 vFuncBool: boolean
515 An indicator for whether the value function should be computed and
516 included in the reported solution.
517 DiscreteShareBool : bool
518 Indicator for whether risky portfolio share should be optimized on the
519 continuous [0,1] interval using the FOC (False), or instead only selected
520 from the discrete set of values in ShareGrid (True). If True, then
521 vFuncBool must also be True.
522 IndepDstnBool : bool
523 Indicator for whether the income and risky return distributions are in-
524 dependent of each other, which can speed up the expectations step.
525 BeqCRRA : float
526 Coefficient of relative risk aversion for warm glow bequest motive.
527 BeqFac : float
528 Multiplicative intensity factor for the warm glow bequest motive.
529 BeqShift : float
530 Stone-Geary shifter in the warm glow bequest motive.
532 Returns
533 -------
534 solution_now : PortfolioSolution
535 Solution to this period's problem.
536 """
537 # Make sure the individual is liquidity constrained. Allowing a consumer to
538 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
539 if BoroCnstArt != 0.0:
540 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
542 # Make sure that if risky portfolio share is optimized only discretely, then
543 # the value function is also constructed (else this task would be impossible).
544 if DiscreteShareBool and (not vFuncBool):
545 raise ValueError(
546 "PortfolioConsumerType requires vFuncBool to be True when DiscreteShareBool is True!"
547 )
549 # Define the current period utility function and effective discount factor
550 uFunc = UtilityFuncCRRA(CRRA)
551 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
552 BeqFacEff = (1.0 - LivPrb) * BeqFac # "effective" bequest factor
553 warm_glow = UtilityFuncStoneGeary(BeqCRRA, BeqFacEff, BeqShift)
555 # Unpack next period's solution for easier access
556 vPfuncAdj_next = solution_next.vPfuncAdj
557 dvdmFuncFxd_next = solution_next.dvdmFuncFxd
558 dvdsFuncFxd_next = solution_next.dvdsFuncFxd
559 vFuncAdj_next = solution_next.vFuncAdj
560 vFuncFxd_next = solution_next.vFuncFxd
562 # Set a flag for whether the natural borrowing constraint is zero, which
563 # depends on whether the smallest transitory income shock is zero
564 BoroCnstNat_iszero = (np.min(IncShkDstn.atoms[1]) == 0.0) or (
565 BeqFac != 0.0 and BeqShift == 0.0
566 )
568 # Prepare to calculate end-of-period marginal values by creating an array
569 # of market resources that the agent could have next period, considering
570 # the grid of end-of-period assets and the distribution of shocks he might
571 # experience next period.
573 # Unpack the risky return shock distribution
574 Risky_next = RiskyDstn.atoms
575 RiskyMax = np.max(Risky_next)
576 RiskyMin = np.min(Risky_next)
578 # bNrm represents R*a, balances after asset return shocks but before income.
579 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
580 if BoroCnstNat_iszero:
581 aNrmGrid = aXtraGrid
582 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
583 else:
584 # Add an asset point at exactly zero
585 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
586 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
588 # Get grid and shock sizes, for easier indexing
589 aNrmCount = aNrmGrid.size
590 ShareCount = ShareGrid.size
592 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
593 bNrmNext, ShareNext = np.meshgrid(bNrmGrid, ShareGrid, indexing="ij")
595 # Define functions that are used internally to evaluate future realizations
596 def calc_mNrm_next(S, b):
597 """
598 Calculate future realizations of market resources mNrm from the income
599 shock distribution S and normalized bank balances b.
600 """
601 return b / (S["PermShk"] * PermGroFac) + S["TranShk"]
603 def calc_dvdm_next(S, b, z):
604 """
605 Evaluate realizations of marginal value of market resources next period,
606 based on the income distribution S, values of bank balances bNrm, and
607 values of the risky share z.
608 """
609 mNrm_next = calc_mNrm_next(S, b)
610 dvdmAdj_next = vPfuncAdj_next(mNrm_next)
612 if AdjustPrb < 1.0:
613 # Expand to the same dimensions as mNrm
614 Share_next_expanded = z + np.zeros_like(mNrm_next)
615 dvdmFxd_next = dvdmFuncFxd_next(mNrm_next, Share_next_expanded)
616 # Combine by adjustment probability
617 dvdm_next = AdjustPrb * dvdmAdj_next + (1.0 - AdjustPrb) * dvdmFxd_next
618 else: # Don't bother evaluating if there's no chance that portfolio share is fixed
619 dvdm_next = dvdmAdj_next
621 dvdm_next = (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
622 return dvdm_next
624 def calc_dvds_next(S, b, z):
625 """
626 Evaluate realizations of marginal value of risky share next period, based
627 on the income distribution S, values of bank balances bNrm, and values of
628 the risky share z.
629 """
630 mNrm_next = calc_mNrm_next(S, b)
632 # No marginal value of Share if it's a free choice!
633 dvdsAdj_next = np.zeros_like(mNrm_next)
635 if AdjustPrb < 1.0:
636 # Expand to the same dimensions as mNrm
637 Share_next_expanded = z + np.zeros_like(mNrm_next)
638 dvdsFxd_next = dvdsFuncFxd_next(mNrm_next, Share_next_expanded)
639 # Combine by adjustment probability
640 dvds_next = AdjustPrb * dvdsAdj_next + (1.0 - AdjustPrb) * dvdsFxd_next
641 else: # Don't bother evaluating if there's no chance that portfolio share is fixed
642 dvds_next = dvdsAdj_next
644 dvds_next = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * dvds_next
645 return dvds_next
647 # Calculate end-of-period marginal value of assets and shares at each point
648 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
649 # values across income and risky return shocks.
651 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
652 dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext, ShareNext))
653 dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0))
654 dvdbNvrsFunc_intermed = BilinearInterp(dvdbNvrs_intermed, bNrmGrid, ShareGrid)
655 dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA)
657 # Calculate intermediate marginal value of risky portfolio share by taking expectations over income shocks
658 dvds_intermed = expected(calc_dvds_next, IncShkDstn, args=(bNrmNext, ShareNext))
659 dvdsFunc_intermed = BilinearInterp(dvds_intermed, bNrmGrid, ShareGrid)
661 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
662 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
664 # Define functions for calculating end-of-period marginal value
665 def calc_EndOfPrd_dvda(S, a, z):
666 """
667 Compute end-of-period marginal value of assets at values a, conditional
668 on risky asset return S and risky share z.
669 """
670 # Calculate future realizations of bank balances bNrm
671 Rxs = S - Rfree # Excess returns
672 Rport = Rfree + z * Rxs # Portfolio return
673 bNrm_next = Rport * a
675 # Ensure shape concordance
676 z_rep = z + np.zeros_like(bNrm_next)
678 # Calculate and return dvda
679 EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next, z_rep)
680 return EndOfPrd_dvda
682 def EndOfPrddvds_dist(S, a, z):
683 """
684 Compute end-of-period marginal value of risky share at values a, conditional
685 on risky asset return S and risky share z.
686 """
687 # Calculate future realizations of bank balances bNrm
688 Rxs = S - Rfree # Excess returns
689 Rport = Rfree + z * Rxs # Portfolio return
690 bNrm_next = Rport * a
692 # Make the shares match the dimension of b, so that it can be vectorized
693 z_rep = z + np.zeros_like(bNrm_next)
695 # Calculate and return dvds
696 EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(
697 bNrm_next, z_rep
698 ) + dvdsFunc_intermed(bNrm_next, z_rep)
699 return EndOfPrd_dvds
701 # Evaluate realizations of value and marginal value after asset returns are realized
703 # Calculate end-of-period marginal value of assets by taking expectations
704 EndOfPrd_dvda = DiscFacEff * expected(
705 calc_EndOfPrd_dvda, RiskyDstn, args=(aNrmNow, ShareNext)
706 )
707 warm_glow_der = warm_glow.der(aNrmNow)
708 EndOfPrd_dvda += np.where(np.isnan(warm_glow_der), 0.0, warm_glow_der)
709 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
711 # Calculate end-of-period marginal value of risky portfolio share by taking expectations
712 EndOfPrd_dvds = DiscFacEff * expected(
713 EndOfPrddvds_dist, RiskyDstn, args=(aNrmNow, ShareNext)
714 )
716 # Make the end-of-period value function if the value function is requested
717 if vFuncBool:
719 def calc_v_intermed(S, b, z):
720 """
721 Calculate "intermediate" value from next period's bank balances, the
722 income shocks S, and the risky asset share.
723 """
724 mNrm_next = calc_mNrm_next(S, b)
726 vAdj_next = vFuncAdj_next(mNrm_next)
727 if AdjustPrb < 1.0:
728 vFxd_next = vFuncFxd_next(mNrm_next, z)
729 # Combine by adjustment probability
730 v_next = AdjustPrb * vAdj_next + (1.0 - AdjustPrb) * vFxd_next
731 else: # Don't bother evaluating if there's no chance that portfolio share is fixed
732 v_next = vAdj_next
734 v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
735 return v_intermed
737 # Calculate intermediate value by taking expectations over income shocks
738 v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext, ShareNext))
740 # Construct the "intermediate value function" for this period
741 vNvrs_intermed = uFunc.inv(v_intermed)
742 vNvrsFunc_intermed = BilinearInterp(vNvrs_intermed, bNrmGrid, ShareGrid)
743 vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA)
745 def calc_EndOfPrd_v(S, a, z):
746 # Calculate future realizations of bank balances bNrm
747 Rxs = S - Rfree
748 Rport = Rfree + z * Rxs
749 bNrm_next = Rport * a
751 # Make an extended share_next of the same dimension as b_nrm so
752 # that the function can be vectorized
753 z_rep = z + np.zeros_like(bNrm_next)
755 EndOfPrd_v = vFunc_intermed(bNrm_next, z_rep)
756 return EndOfPrd_v
758 # Calculate end-of-period value by taking expectations
759 EndOfPrd_v = DiscFacEff * expected(
760 calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext)
761 )
762 EndOfPrd_v += warm_glow(aNrmNow)
763 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
765 # Now make an end-of-period value function over aNrm and Share
766 EndOfPrd_vNvrsFunc = BilinearInterp(EndOfPrd_vNvrs, aNrmGrid, ShareGrid)
767 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
768 # This will be used later to make the value function for this period
770 # Find the optimal risky asset share either by choosing the best value among
771 # the discrete grid choices, or by satisfying the FOC with equality (continuous)
772 if DiscreteShareBool:
773 # If we're restricted to discrete choices, then portfolio share is
774 # the one with highest value for each aNrm gridpoint
775 opt_idx = np.argmax(EndOfPrd_v, axis=1)
776 ShareAdj_now = ShareGrid[opt_idx]
778 # Take cNrm at that index as well... and that's it!
779 cNrmAdj_now = EndOfPrd_dvdaNvrs[np.arange(aNrmCount), opt_idx]
781 else:
782 # Now find the optimal (continuous) risky share on [0,1] by solving the first
783 # order condition EndOfPrd_dvds == 0.
784 FOC_s = EndOfPrd_dvds # Relabel for convenient typing
786 # For each value of aNrm, find the value of Share such that FOC_s == 0
787 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0)
788 share_idx = np.argmax(crossing, axis=1)
789 # This represents the index of the segment of the share grid where dvds flips
790 # from positive to negative, indicating that there's a zero *on* the segment
792 # Calculate the fractional distance between those share gridpoints where the
793 # zero should be found, assuming a linear function; call it alpha
794 a_idx = np.arange(aNrmCount)
795 bot_s = ShareGrid[share_idx]
796 top_s = ShareGrid[share_idx + 1]
797 bot_f = FOC_s[a_idx, share_idx]
798 top_f = FOC_s[a_idx, share_idx + 1]
799 bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx]
800 top_c = EndOfPrd_dvdaNvrs[a_idx, share_idx + 1]
801 alpha = 1.0 - top_f / (top_f - bot_f)
803 # Calculate the continuous optimal risky share and optimal consumption
804 ShareAdj_now = (1.0 - alpha) * bot_s + alpha * top_s
805 cNrmAdj_now = (1.0 - alpha) * bot_c + alpha * top_c
807 # If agent wants to put more than 100% into risky asset, he is constrained.
808 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
809 constrained_top = FOC_s[:, -1] > 0.0
810 constrained_bot = FOC_s[:, 0] < 0.0
812 # Apply those constraints to both risky share and consumption (but lower
813 # constraint should never be relevant)
814 ShareAdj_now[constrained_top] = 1.0
815 ShareAdj_now[constrained_bot] = 0.0
816 cNrmAdj_now[constrained_top] = EndOfPrd_dvdaNvrs[constrained_top, -1]
817 cNrmAdj_now[constrained_bot] = EndOfPrd_dvdaNvrs[constrained_bot, 0]
819 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
820 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
821 # can't depend on the risky share if it doesn't meaningfully exist. Apply
822 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
823 if not BoroCnstNat_iszero:
824 ShareAdj_now[0] = 1.0
825 cNrmAdj_now[0] = EndOfPrd_dvdaNvrs[0, -1]
827 # Construct functions characterizing the solution for this period
829 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
830 # then construct the consumption function when the agent can adjust his share
831 mNrmAdj_now = np.insert(aNrmGrid + cNrmAdj_now, 0, 0.0)
832 cNrmAdj_now = np.insert(cNrmAdj_now, 0, 0.0)
833 cFuncAdj_now = LinearInterp(mNrmAdj_now, cNrmAdj_now)
835 # Construct the marginal value (of mNrm) function when the agent can adjust
836 vPfuncAdj_now = MargValueFuncCRRA(cFuncAdj_now, CRRA)
838 # Construct the consumption function when the agent *can't* adjust the risky
839 # share, as well as the marginal value of Share function
840 cFuncFxd_by_Share = []
841 dvdsFuncFxd_by_Share = []
842 for j in range(ShareCount):
843 cNrmFxd_temp = np.insert(EndOfPrd_dvdaNvrs[:, j], 0, 0.0)
844 mNrmFxd_temp = np.insert(aNrmGrid + cNrmFxd_temp[1:], 0, 0.0)
845 dvdsFxd_temp = np.insert(EndOfPrd_dvds[:, j], 0, EndOfPrd_dvds[0, j])
846 cFuncFxd_by_Share.append(LinearInterp(mNrmFxd_temp, cNrmFxd_temp))
847 dvdsFuncFxd_by_Share.append(LinearInterp(mNrmFxd_temp, dvdsFxd_temp))
848 cFuncFxd_now = LinearInterpOnInterp1D(cFuncFxd_by_Share, ShareGrid)
849 dvdsFuncFxd_now = LinearInterpOnInterp1D(dvdsFuncFxd_by_Share, ShareGrid)
851 # The share function when the agent can't adjust his portfolio is trivial
852 ShareFuncFxd_now = IdentityFunction(i_dim=1, n_dims=2)
854 # Construct the marginal value of mNrm function when the agent can't adjust his share
855 dvdmFuncFxd_now = MargValueFuncCRRA(cFuncFxd_now, CRRA)
857 # Construct the optimal risky share function when adjusting is possible.
858 # The interpolation method depends on whether the choice is discrete or continuous.
859 if DiscreteShareBool:
860 # If the share choice is discrete, the "interpolated" share function acts
861 # like a step function, with jumps at the midpoints of mNrm gridpoints.
862 # Because an actual step function would break our (assumed continuous) linear
863 # interpolator, there's a *tiny* region with extremely high slope.
864 mNrmAdj_mid = (mNrmAdj_now[2:] + mNrmAdj_now[1:-1]) / 2
865 mNrmAdj_plus = mNrmAdj_mid * (1.0 + 1e-12)
866 mNrmAdj_comb = (np.transpose(np.vstack((mNrmAdj_mid, mNrmAdj_plus)))).flatten()
867 mNrmAdj_comb = np.append(np.insert(mNrmAdj_comb, 0, 0.0), mNrmAdj_now[-1])
868 Share_comb = (np.transpose(np.vstack((ShareAdj_now, ShareAdj_now)))).flatten()
869 ShareFuncAdj_now = LinearInterp(mNrmAdj_comb, Share_comb)
871 else:
872 # If the share choice is continuous, just make an ordinary interpolating function
873 if BoroCnstNat_iszero:
874 Share_lower_bound = ShareLimit
875 else:
876 Share_lower_bound = 1.0
877 ShareAdj_now = np.insert(ShareAdj_now, 0, Share_lower_bound)
878 ShareFuncAdj_now = LinearInterp(mNrmAdj_now, ShareAdj_now, ShareLimit, 0.0)
880 # Add the value function if requested
881 if vFuncBool:
882 # Create the value functions for this period, defined over market resources
883 # mNrm when agent can adjust his portfolio, and over market resources and
884 # fixed share when agent can not adjust his portfolio.
886 # Construct the value function when the agent can adjust his portfolio
887 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
888 cNrm_temp = cFuncAdj_now(mNrm_temp)
889 aNrm_temp = mNrm_temp - cNrm_temp
890 Share_temp = ShareFuncAdj_now(mNrm_temp)
891 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
892 vNvrs_temp = uFunc.inv(v_temp)
893 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
894 vNvrsFuncAdj = CubicInterp(
895 np.insert(mNrm_temp, 0, 0.0), # x_list
896 np.insert(vNvrs_temp, 0, 0.0), # f_list
897 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
898 )
899 # Re-curve the pseudo-inverse value function
900 vFuncAdj_now = ValueFuncCRRA(vNvrsFuncAdj, CRRA)
902 # Construct the value function when the agent *can't* adjust his portfolio
903 mNrm_temp, Share_temp = np.meshgrid(aXtraGrid, ShareGrid)
904 cNrm_temp = cFuncFxd_now(mNrm_temp, Share_temp)
905 aNrm_temp = mNrm_temp - cNrm_temp
906 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
907 vNvrs_temp = uFunc.inv(v_temp)
908 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
909 vNvrsFuncFxd_by_Share = []
910 for j in range(ShareCount):
911 vNvrsFuncFxd_by_Share.append(
912 CubicInterp(
913 np.insert(mNrm_temp[:, 0], 0, 0.0), # x_list
914 np.insert(vNvrs_temp[:, j], 0, 0.0), # f_list
915 np.insert(vNvrsP_temp[:, j], 0, vNvrsP_temp[j, 0]), # dfdx_list
916 )
917 )
918 vNvrsFuncFxd = LinearInterpOnInterp1D(vNvrsFuncFxd_by_Share, ShareGrid)
919 vFuncFxd_now = ValueFuncCRRA(vNvrsFuncFxd, CRRA)
921 else: # If vFuncBool is False, fill in dummy values
922 vFuncAdj_now = NullFunc()
923 vFuncFxd_now = NullFunc()
925 # Package and return the solution
926 solution_now = PortfolioSolution(
927 cFuncAdj=cFuncAdj_now,
928 ShareFuncAdj=ShareFuncAdj_now,
929 vPfuncAdj=vPfuncAdj_now,
930 vFuncAdj=vFuncAdj_now,
931 cFuncFxd=cFuncFxd_now,
932 ShareFuncFxd=ShareFuncFxd_now,
933 dvdmFuncFxd=dvdmFuncFxd_now,
934 dvdsFuncFxd=dvdsFuncFxd_now,
935 vFuncFxd=vFuncFxd_now,
936 AdjPrb=AdjustPrb,
937 )
938 return solution_now
941###############################################################################
943# Make a dictionary of constructors for the warm glow bequest model
944warmglow_constructor_dict = {
945 "IncShkDstn": construct_lognormal_income_process_unemployment,
946 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
947 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
948 "aXtraGrid": make_assets_grid,
949 "solution_terminal": make_bequest_solution_terminal,
950 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
951 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
952}
954# Make a dictionary with parameters for the default constructor for kNrmInitDstn
955default_kNrmInitDstn_params = {
956 "kLogInitMean": -12.0, # Mean of log initial capital
957 "kLogInitStd": 0.0, # Stdev of log initial capital
958 "kNrmInitCount": 15, # Number of points in initial capital discretization
959}
961# Make a dictionary with parameters for the default constructor for pLvlInitDstn
962default_pLvlInitDstn_params = {
963 "pLogInitMean": 0.0, # Mean of log permanent income
964 "pLogInitStd": 0.0, # Stdev of log permanent income
965 "pLvlInitCount": 15, # Number of points in initial capital discretization
966}
968# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
969default_IncShkDstn_params = {
970 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
971 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
972 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
973 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
974 "UnempPrb": 0.05, # Probability of unemployment while working
975 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
976 "T_retire": 0, # Period of retirement (0 --> no retirement)
977 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
978 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
979}
981# Default parameters to make aXtraGrid using make_assets_grid
982default_aXtraGrid_params = {
983 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
984 "aXtraMax": 20, # Maximum end-of-period "assets above minimum" value
985 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
986 "aXtraCount": 48, # Number of points in the grid of "assets above minimum"
987 "aXtraExtra": None, # Additional other values to add in grid (optional)
988}
990# Make a dictionary to specify awarm glow bequest consumer type
991init_warm_glow = {
992 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
993 "cycles": 1, # Finite, non-cyclic model
994 "T_cycle": 1, # Number of periods in the cycle for this agent type
995 "constructors": warmglow_constructor_dict, # See dictionary above
996 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
997 "CRRA": 2.0, # Coefficient of relative risk aversion on consumption
998 "Rfree": [1.03], # Interest factor on retained assets
999 "DiscFac": 0.96, # Intertemporal discount factor
1000 "LivPrb": [0.98], # Survival probability after each period
1001 "PermGroFac": [1.01], # Permanent income growth factor
1002 "BoroCnstArt": 0.0, # Artificial borrowing constraint
1003 "BeqCRRA": 2.0, # Coefficient of relative risk aversion for bequest motive
1004 "BeqFac": 40.0, # Scaling factor for bequest motive
1005 "BeqShift": 0.0, # Stone-Geary shifter term for bequest motive
1006 "BeqCRRATerm": 2.0, # Coefficient of relative risk aversion for bequest motive, terminal period only
1007 "BeqFacTerm": 40.0, # Scaling factor for bequest motive, terminal period only
1008 "BeqShiftTerm": 0.0, # Stone-Geary shifter term for bequest motive, terminal period only
1009 "vFuncBool": False, # Whether to calculate the value function during solution
1010 "CubicBool": False, # Whether to use cubic spline interpolation when True
1011 # (Uses linear spline interpolation for cFunc when False)
1012 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1013 "AgentCount": 10000, # Number of agents of this type
1014 "T_age": None, # Age after which simulated agents are automatically killed
1015 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1016 # (The portion of PermGroFac attributable to aggregate productivity growth)
1017 "NewbornTransShk": False, # Whether Newborns have transitory shock
1018 # ADDITIONAL OPTIONAL PARAMETERS
1019 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1020 # (Forces Newborns to follow solution path of the agent they replaced if True)
1021 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1022}
1023init_warm_glow.update(default_IncShkDstn_params)
1024init_warm_glow.update(default_aXtraGrid_params)
1025init_warm_glow.update(default_kNrmInitDstn_params)
1026init_warm_glow.update(default_pLvlInitDstn_params)
1028# Make a dictionary with bequest motives turned off
1029init_accidental_bequest = init_warm_glow.copy()
1030init_accidental_bequest["BeqFac"] = 0.0
1031init_accidental_bequest["BeqShift"] = 0.0
1032init_accidental_bequest["BeqFacTerm"] = 0.0
1033init_accidental_bequest["BeqShiftTerm"] = 0.0
1035# Make a dictionary that has *only* a terminal period bequest
1036init_warm_glow_terminal_only = init_warm_glow.copy()
1037init_warm_glow_terminal_only["BeqFac"] = 0.0
1038init_warm_glow_terminal_only["BeqShift"] = 0.0
1041class BequestWarmGlowConsumerType(IndShockConsumerType):
1042 r"""
1043 A consumer type with based on IndShockConsumerType, with an additional bequest motive.
1044 They gain utility for any wealth they leave when they die, according to a Stone-Geary utility.
1046 .. math::
1047 \newcommand{\CRRA}{\rho}
1048 \newcommand{\DiePrb}{\mathsf{D}}
1049 \newcommand{\PermGroFac}{\Gamma}
1050 \newcommand{\Rfree}{\mathsf{R}}
1051 \newcommand{\DiscFac}{\beta}
1052 \begin{align*}
1053 v_t(m_t) &= \max_{c_t}u(c_t) + \DiePrb_{t+1} u_{Beq}(a_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], \\
1054 & \text{s.t.} \\
1055 a_t &= m_t - c_t, \\
1056 a_t &\geq \underline{a}, \\
1057 m_{t+1} &= a_t \Rfree_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
1058 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
1059 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1, \\
1060 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA} \\
1061 u_{Beq} (a) &= \textbf{BeqFac} \frac{(a+\textbf{BeqShift})^{1-\CRRA_{Beq}}}{1-\CRRA_{Beq}} \\
1062 \end{align*}
1065 Constructors
1066 ------------
1067 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
1068 The agent's income shock distributions.
1070 It's default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
1071 aXtraGrid: Constructor
1072 The agent's asset grid.
1074 It's default constructor is :func:`HARK.utilities.make_assets_grid`
1076 Solving Parameters
1077 ------------------
1078 cycles: int
1079 0 specifies an infinite horizon model, 1 specifies a finite model.
1080 T_cycle: int
1081 Number of periods in the cycle for this agent type.
1082 CRRA: float, :math:`\rho`
1083 Coefficient of Relative Risk Aversion.
1084 BeqCRRA: float, :math:`\rho_{Beq}`
1085 Coefficient of Relative Risk Aversion for the bequest motive.
1086 If this value isn't the same as CRRA, then the model can only be represented as a Bellman equation.
1087 This may cause unintented behavior.
1088 BeqCRRATerm: float, :math:`\rho_{Beq}`
1089 The Coefficient of Relative Risk Aversion for the bequest motive, but only in the terminal period.
1090 In most cases this should be the same as beqCRRA.
1091 BeqShift: float, :math:`\textbf{BeqShift}`
1092 The Shift term from the bequest motive's utility function.
1093 If this value isn't 0, then the model can only be represented as a Bellman equation.
1094 This may cause unintented behavior.
1095 BeqShiftTerm: float, :math:`\textbf{BeqShift}`
1096 The shift term from the bequest motive's utility function, in the terminal period.
1097 In most cases this should be the same as beqShift
1098 BeqFac: float, :math:`\textbf{BeqFac}`
1099 The weight for the bequest's utility function.
1100 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1101 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1102 DiscFac: float, :math:`\beta`
1103 Intertemporal discount factor.
1104 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1105 Survival probability after each period.
1106 PermGroFac: list[float], time varying, :math:`\Gamma`
1107 Permanent income growth factor.
1108 BoroCnstArt: float, :math:`\underline{a}`
1109 The minimum Asset/Perminant Income ratio, None to ignore.
1110 vFuncBool: bool
1111 Whether to calculate the value function during solution.
1112 CubicBool: bool
1113 Whether to use cubic spline interpoliation.
1115 Simulation Parameters
1116 ---------------------
1117 AgentCount: int
1118 Number of agents of this kind that are created during simulations.
1119 T_age: int
1120 Age after which to automatically kill agents, None to ignore.
1121 T_sim: int, required for simulation
1122 Number of periods to simulate.
1123 track_vars: list[strings]
1124 List of variables that should be tracked when running the simulation.
1125 For this agent, the options are 'PermShk', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
1127 PermShk is the agent's permanent income shock
1129 TranShk is the agent's transitory income shock
1131 aLvl is the nominal asset level
1133 aNrm is the normalized assets
1135 bNrm is the normalized resources without this period's labor income
1137 cNrm is the normalized consumption
1139 mNrm is the normalized market resources
1141 pLvl is the permanent income level
1143 who_dies is the array of which agents died
1144 aNrmInitMean: float
1145 Mean of Log initial Normalized Assets.
1146 aNrmInitStd: float
1147 Std of Log initial Normalized Assets.
1148 pLvlInitMean: float
1149 Mean of Log initial permanent income.
1150 pLvlInitStd: float
1151 Std of Log initial permanent income.
1152 PermGroFacAgg: float
1153 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1154 PerfMITShk: boolean
1155 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1156 NewbornTransShk: boolean
1157 Whether Newborns have transitory shock.
1159 Attributes
1160 ----------
1161 solution: list[Consumer solution object]
1162 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1163 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1165 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1166 history: Dict[Array]
1167 Created by running the :func:`.simulate()` method.
1168 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1169 Visit :class:`HARK.core.AgentType.simulate` for more information.
1170 """
1172 time_inv_ = IndShockConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"]
1173 default_ = {
1174 "params": init_accidental_bequest,
1175 "solver": solve_one_period_ConsWarmBequest,
1176 "model": "ConsIndShock.yaml",
1177 }
1179 def pre_solve(self):
1180 self.construct("solution_terminal")
1182 def check_conditions(self, verbose=None):
1183 raise NotImplementedError()
1185 def calc_limiting_values(self):
1186 raise NotImplementedError()
1189###############################################################################
1192# Make a dictionary of constructors for the portfolio choice consumer type
1193portfolio_bequest_constructor_dict = {
1194 "IncShkDstn": construct_lognormal_income_process_unemployment,
1195 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1196 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1197 "aXtraGrid": make_assets_grid,
1198 "RiskyDstn": make_lognormal_RiskyDstn,
1199 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
1200 "ShareLimit": calc_ShareLimit_for_CRRA,
1201 "ShareGrid": make_simple_ShareGrid,
1202 "AdjustDstn": make_AdjustDstn,
1203 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1204 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1205 "solution_terminal": make_warmglow_portfolio_solution_terminal,
1206}
1208# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1209default_kNrmInitDstn_params = {
1210 "kLogInitMean": -12.0, # Mean of log initial capital
1211 "kLogInitStd": 0.0, # Stdev of log initial capital
1212 "kNrmInitCount": 15, # Number of points in initial capital discretization
1213}
1215# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1216default_pLvlInitDstn_params = {
1217 "pLogInitMean": 0.0, # Mean of log permanent income
1218 "pLogInitStd": 0.0, # Stdev of log permanent income
1219 "pLvlInitCount": 15, # Number of points in initial capital discretization
1220}
1223# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
1224default_IncShkDstn_params = {
1225 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1226 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1227 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1228 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1229 "UnempPrb": 0.05, # Probability of unemployment while working
1230 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1231 "T_retire": 0, # Period of retirement (0 --> no retirement)
1232 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1233 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1234}
1236# Default parameters to make aXtraGrid using make_assets_grid
1237default_aXtraGrid_params = {
1238 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1239 "aXtraMax": 100, # Maximum end-of-period "assets above minimum" value
1240 "aXtraNestFac": 1, # Exponential nesting factor for aXtraGrid
1241 "aXtraCount": 200, # Number of points in the grid of "assets above minimum"
1242 "aXtraExtra": None, # Additional other values to add in grid (optional)
1243}
1245# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn (and uniform ShareGrid)
1246default_RiskyDstn_and_ShareGrid_params = {
1247 "RiskyAvg": 1.08, # Mean return factor of risky asset
1248 "RiskyStd": 0.18362634887, # Stdev of log returns on risky asset
1249 "RiskyCount": 5, # Number of integration nodes to use in approximation of risky returns
1250 "ShareCount": 25, # Number of discrete points in the risky share approximation
1251}
1253# Make a dictionary to specify a risky asset consumer type
1254init_portfolio_bequest = {
1255 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1256 "cycles": 1, # Finite, non-cyclic model
1257 "T_cycle": 1, # Number of periods in the cycle for this agent type
1258 "constructors": portfolio_bequest_constructor_dict, # See dictionary above
1259 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
1260 "CRRA": 5.0, # Coefficient of relative risk aversion
1261 "Rfree": [1.03], # Return factor on risk free asset
1262 "DiscFac": 0.90, # Intertemporal discount factor
1263 "LivPrb": [0.98], # Survival probability after each period
1264 "PermGroFac": [1.01], # Permanent income growth factor
1265 "BoroCnstArt": 0.0, # Artificial borrowing constraint
1266 "BeqCRRA": 2.0, # Coefficient of relative risk aversion for bequest motive
1267 "BeqFac": 40.0, # Scaling factor for bequest motive
1268 "BeqShift": 0.0, # Stone-Geary shifter term for bequest motive
1269 "BeqCRRATerm": 2.0, # Coefficient of relative risk aversion for bequest motive, terminal period only
1270 "BeqFacTerm": 40.0, # Scaling factor for bequest motive, terminal period only
1271 "BeqShiftTerm": 0.0, # Stone-Geary shifter term for bequest motive, terminal period only
1272 "DiscreteShareBool": False, # Whether risky asset share is restricted to discrete values
1273 "vFuncBool": False, # Whether to calculate the value function during solution
1274 "CubicBool": False, # Whether to use cubic spline interpolation when True
1275 # (Uses linear spline interpolation for cFunc when False)
1276 "IndepDstnBool": True, # Indicator for whether return & income shocks are independent
1277 "PortfolioBool": True, # Whether this agent has portfolio choice
1278 "PortfolioBisect": False, # What does this do?
1279 "AdjustPrb": 1.0, # Probability that the agent can update their risky portfolio share each period
1280 "RiskyShareFixed": None, # This just needs to exist because of inheritance, does nothing
1281 "sim_common_Rrisky": True, # Whether risky returns have a shared/common value across agents
1282 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1283 "AgentCount": 10000, # Number of agents of this type
1284 "T_age": None, # Age after which simulated agents are automatically killed
1285 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1286 # (The portion of PermGroFac attributable to aggregate productivity growth)
1287 "NewbornTransShk": False, # Whether Newborns have transitory shock
1288 # ADDITIONAL OPTIONAL PARAMETERS
1289 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1290 # (Forces Newborns to follow solution path of the agent they replaced if True)
1291 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1292}
1293init_portfolio_bequest.update(default_kNrmInitDstn_params)
1294init_portfolio_bequest.update(default_pLvlInitDstn_params)
1295init_portfolio_bequest.update(default_IncShkDstn_params)
1296init_portfolio_bequest.update(default_aXtraGrid_params)
1297init_portfolio_bequest.update(default_RiskyDstn_and_ShareGrid_params)
1300class BequestWarmGlowPortfolioType(PortfolioConsumerType):
1301 r"""
1302 A consumer type with based on PortfolioConsumerType, with an additional bequest motive.
1303 They gain utility for any wealth they leave when they die, according to a Stone-Geary utility.
1305 .. math::
1306 \newcommand{\CRRA}{\rho}
1307 \newcommand{\DiePrb}{\mathsf{D}}
1308 \newcommand{\PermGroFac}{\Gamma}
1309 \newcommand{\Rfree}{\mathsf{R}}
1310 \newcommand{\DiscFac}{\beta}
1311 \begin{align*}
1312 v_t(m_t,S_t) &= \max_{c_t,S^{*}_t} u(c_t) + \DiePrb_{t+1} u_{Beq}(a_t)+ \DiscFac (1-\DiePrb_{t+1}) \mathbb{E}_{t} \left[(\PermGroFac_{t+1}\psi_{t+1})^{1-\CRRA} v_{t+1}(m_{t+1},S_{t+1}) \right], \\
1313 & \text{s.t.} \\
1314 a_t &= m_t - c_t, \\
1315 a_t &\geq \underline{a}, \\
1316 m_{t+1} &= \mathsf{R}_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) a_t + \theta_{t+1}, \\
1317 \mathsf{R}_{t+1} &=S_t\phi_{t+1}\mathbf{R}_{t+1}+ (1-S_t)\mathsf{R}_{t+1}, \\
1318 S_{t+1} &= \begin{cases}
1319 S^{*}_t & \text{if } p_t < \wp\\
1320 S_t & \text{if } p_t \geq \wp,
1321 \end{cases}\\
1322 (\psi_{t+1},\theta_{t+1},\phi_{t+1},p_t) &\sim F_{t+1}, \\
1323 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1. \\
1324 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA} \\
1325 u_{Beq} (a) &= \textbf{BeqFac} \frac{(a+\textbf{BeqShift})^{1-\CRRA_{Beq}}}{1-\CRRA_{Beq}} \\
1326 \end{align*}
1329 Constructors
1330 ------------
1331 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
1332 The agent's income shock distributions.
1334 It's default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
1335 aXtraGrid: Constructor
1336 The agent's asset grid.
1338 It's default constructor is :func:`HARK.utilities.make_assets_grid`
1339 ShareGrid: Constructor
1340 The agent's risky asset share grid
1342 It's default constructor is :func:`HARK.ConsumptionSaving.ConsRiskyAssetModel.make_simple_ShareGrid`
1343 RiskyDstn: Constructor, :math:`\phi`
1344 The agent's asset shock distribution for risky assets.
1346 It's default constructor is :func:`HARK.Calibration.Assets.AssetProcesses.make_lognormal_RiskyDstn`
1348 Solving Parameters
1349 ------------------
1350 cycles: int
1351 0 specifies an infinite horizon model, 1 specifies a finite model.
1352 T_cycle: int
1353 Number of periods in the cycle for this agent type.
1354 CRRA: float, :math:`\rho`
1355 Coefficient of Relative Risk Aversion.
1356 BeqCRRA: float, :math:`\rho_{Beq}`
1357 Coefficient of Relative Risk Aversion for the bequest motive.
1358 If this value isn't the same as CRRA, then the model can only be represented as a Bellman equation.
1359 This may cause unintented behavior.
1360 BeqCRRATerm: float, :math:`\rho_{Beq}`
1361 The Coefficient of Relative Risk Aversion for the bequest motive, but only in the terminal period.
1362 In most cases this should be the same as beqCRRA.
1363 BeqShift: float, :math:`\textbf{BeqShift}`
1364 The Shift term from the bequest motive's utility function.
1365 If this value isn't 0, then the model can only be represented as a Bellman equation.
1366 This may cause unintented behavior.
1367 BeqShiftTerm: float, :math:`\textbf{BeqShift}`
1368 The shift term from the bequest motive's utility function, in the terminal period.
1369 In most cases this should be the same as beqShift
1370 BeqFac: float, :math:`\textbf{BeqFac}`
1371 The weight for the bequest's utility function.
1372 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1373 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1374 DiscFac: float, :math:`\beta`
1375 Intertemporal discount factor.
1376 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1377 Survival probability after each period.
1378 PermGroFac: list[float], time varying, :math:`\Gamma`
1379 Permanent income growth factor.
1380 BoroCnstArt: float, default=0.0, :math:`\underline{a}`
1381 The minimum Asset/Perminant Income ratio. for this agent, BoroCnstArt must be 0.
1382 vFuncBool: bool
1383 Whether to calculate the value function during solution.
1384 CubicBool: bool
1385 Whether to use cubic spline interpoliation.
1386 AdjustPrb: float or list[float], time varying
1387 Must be between 0 and 1. Probability that the agent can update their risky portfolio share each period. Pass a list of floats to make AdjustPrb time varying.
1389 Simulation Parameters
1390 ---------------------
1391 sim_common_Rrisky: Boolean
1392 Whether risky returns have a shared/common value across agents. If True, Risky return's can't be time varying.
1393 AgentCount: int
1394 Number of agents of this kind that are created during simulations.
1395 T_age: int
1396 Age after which to automatically kill agents, None to ignore.
1397 T_sim: int, required for simulation
1398 Number of periods to simulate.
1399 track_vars: list[strings]
1400 List of variables that should be tracked when running the simulation.
1401 For this agent, the options are 'Adjust', 'PermShk', 'Risky', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
1403 Adjust is the array of which agents can adjust
1405 PermShk is the agent's permanent income shock
1407 Risky is the agent's risky asset shock
1409 TranShk is the agent's transitory income shock
1411 aLvl is the nominal asset level
1413 aNrm is the normalized assets
1415 bNrm is the normalized resources without this period's labor income
1417 cNrm is the normalized consumption
1419 mNrm is the normalized market resources
1421 pLvl is the permanent income level
1423 who_dies is the array of which agents died
1424 aNrmInitMean: float
1425 Mean of Log initial Normalized Assets.
1426 aNrmInitStd: float
1427 Std of Log initial Normalized Assets.
1428 pLvlInitMean: float
1429 Mean of Log initial permanent income.
1430 pLvlInitStd: float
1431 Std of Log initial permanent income.
1432 PermGroFacAgg: float
1433 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1434 PerfMITShk: boolean
1435 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1436 NewbornTransShk: boolean
1437 Whether Newborns have transitory shock.
1439 Attributes
1440 ----------
1441 solution: list[Consumer solution object]
1442 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1443 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1445 Visit :class:`HARK.ConsumptionSaving.ConsPortfolioModel.PortfolioSolution` for more information about the solution.
1447 history: Dict[Array]
1448 Created by running the :func:`.simulate()` method.
1449 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1450 Visit :class:`HARK.core.AgentType.simulate` for more information.
1451 """
1453 time_inv_ = PortfolioConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"]
1454 default_ = {
1455 "params": init_portfolio_bequest,
1456 "solver": solve_one_period_ConsPortfolioWarmGlow,
1457 "model": "ConsRiskyAsset.yaml",
1458 }