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