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