Coverage for HARK / ConsumptionSaving / ConsRiskyAssetModel.py: 98%
488 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 06:00 +0000
1"""
2This file contains a class that adds a risky asset with a log-normal return
3factor to IndShockConsumerType. It is meant as a container of methods for dealing
4with risky assets that will be useful to models what will inherit from it.
5"""
7import numpy as np
9from HARK import NullFunc
10from HARK.ConsumptionSaving.ConsIndShockModel import (
11 ConsumerSolution,
12 IndShockConsumerType,
13 make_basic_CRRA_solution_terminal,
14 make_lognormal_kNrm_init_dstn,
15 make_lognormal_pLvl_init_dstn,
16)
17from HARK.Calibration.Assets.AssetProcesses import (
18 make_lognormal_RiskyDstn,
19 combine_IncShkDstn_and_RiskyDstn,
20 calc_ShareLimit_for_CRRA,
21)
22from HARK.Calibration.Income.IncomeProcesses import (
23 construct_lognormal_income_process_unemployment,
24 get_PermShkDstn_from_IncShkDstn,
25 get_TranShkDstn_from_IncShkDstn,
26)
27from HARK.distributions import (
28 Bernoulli,
29 expected,
30 IndexDistribution,
31)
32from HARK.interpolation import (
33 BilinearInterp,
34 ConstantFunction,
35 LinearInterp,
36 LinearInterpOnInterp1D,
37 LowerEnvelope,
38 CubicInterp,
39 MargMargValueFuncCRRA,
40 MargValueFuncCRRA,
41 ValueFuncCRRA,
42)
43from HARK.rewards import UtilityFuncCRRA
44from HARK.utilities import make_assets_grid
46###############################################################################
49def make_simple_ShareGrid(ShareCount):
50 """
51 Make a uniformly spaced grid on the unit interval, representing risky asset shares.
53 Parameters
54 ----------
55 ShareCount : int
56 Number of points in the grid.
58 Returns
59 -------
60 ShareGrid : np.array
61 """
62 ShareGrid = np.linspace(0.0, 1.0, ShareCount)
63 return ShareGrid
66def select_risky_solver(RiskyShareFixed):
67 """
68 Trivial constructor function that chooses between two solvers.
69 """
70 if RiskyShareFixed is None:
71 solve_one_period = solve_one_period_ConsPortChoice
72 else:
73 solve_one_period = solve_one_period_ConsIndShockRiskyAsset
74 return solve_one_period
77def make_AdjustDstn(AdjustPrb, T_cycle, RNG):
78 """
79 Make the distribution of "allowed to adjust" outcomes (a Bernoulli dstn) that
80 could depend on age.
82 Parameters
83 ----------
84 AdjustPrb : float or [float]
85 Probability of being allowed to adjust portfolio allocation, by period of cycle.
86 T_cycle : int
87 Number of periods in the cycle.
88 RNG : RandomState
89 Instance's own random number generator.
91 Returns
92 -------
93 AdjustDstn : BernoulliDistribution or IndexDistribution
94 Distribution object for whether agents can update their portfolios.
95 """
96 if type(AdjustPrb) is list and (len(AdjustPrb) == T_cycle):
97 AdjustDstn = IndexDistribution(
98 Bernoulli, {"p": AdjustPrb}, seed=RNG.integers(0, 2**31 - 1)
99 )
100 elif type(AdjustPrb) is list:
101 raise AttributeError(
102 "If AdjustPrb is time-varying, it must have length of T_cycle!"
103 )
104 else:
105 AdjustDstn = Bernoulli(p=AdjustPrb, seed=RNG.integers(0, 2**31 - 1))
106 return AdjustDstn
109###############################################################################
111# Make a dictionary of constructors for the risky asset model
112IndShockRiskyAssetConsumerType_constructor_default = {
113 "IncShkDstn": construct_lognormal_income_process_unemployment,
114 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
115 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
116 "aXtraGrid": make_assets_grid,
117 "RiskyDstn": make_lognormal_RiskyDstn,
118 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
119 "ShareLimit": calc_ShareLimit_for_CRRA,
120 "ShareGrid": make_simple_ShareGrid,
121 "AdjustDstn": make_AdjustDstn,
122 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
123 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
124 "solution_terminal": make_basic_CRRA_solution_terminal,
125 "solve_one_period": select_risky_solver,
126}
128# Make a dictionary with parameters for the default constructor for kNrmInitDstn
129IndShockRiskyAssetConsumerType_kNrmInitDstn_default = {
130 "kLogInitMean": -12.0, # Mean of log initial capital
131 "kLogInitStd": 0.0, # Stdev of log initial capital
132 "kNrmInitCount": 15, # Number of points in initial capital discretization
133}
135# Make a dictionary with parameters for the default constructor for pLvlInitDstn
136IndShockRiskyAssetConsumerType_pLvlInitDstn_default = {
137 "pLogInitMean": 0.0, # Mean of log permanent income
138 "pLogInitStd": 0.0, # Stdev of log permanent income
139 "pLvlInitCount": 15, # Number of points in initial capital discretization
140}
142# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
143IndShockRiskyAssetConsumerType_IncShkDstn_default = {
144 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
145 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
146 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
147 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
148 "UnempPrb": 0.05, # Probability of unemployment while working
149 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
150 "T_retire": 0, # Period of retirement (0 --> no retirement)
151 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
152 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
153}
155# Default parameters to make aXtraGrid using make_assets_grid
156IndShockRiskyAssetConsumerType_aXtraGrid_default = {
157 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
158 "aXtraMax": 100.0, # Maximum end-of-period "assets above minimum" value
159 "aXtraNestFac": 2, # Exponential nesting factor for aXtraGrid
160 "aXtraCount": 200, # Number of points in the grid of "assets above minimum"
161 "aXtraExtra": None, # Additional other values to add in grid (optional)
162}
164# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn
165IndShockRiskyAssetConsumerType_RiskyDstn_default = {
166 "RiskyAvg": 1.0803701891, # Mean return factor of risky asset
167 "RiskyStd": 0.162927447983, # Stdev of log returns on risky asset
168 "RiskyCount": 5, # Number of integration nodes to use in approximation of risky returns
169}
170# Risky return factor moments are based on SP500 real returns from Shiller's
171# "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
172# Access it through the internet archive
173# We have (or will) rounded them to the nearest .01
175# Default parameters to make RiskyDstn with make_simple_ShareGrid
176IndShockRiskyAssetConsumerType_ShareGrid_default = {
177 "ShareCount": 25, # Number of discrete points in the risky share approximation
178}
180# Make a dictionary to specify a risky asset consumer type
181IndShockRiskyAssetConsumerType_solving_default = {
182 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
183 "cycles": 1, # Finite, non-cyclic model
184 "T_cycle": 1, # Number of periods in the cycle for this agent type
185 "constructors": IndShockRiskyAssetConsumerType_constructor_default, # See dictionary above
186 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
187 "CRRA": 2.0, # Coefficient of relative risk aversion
188 "Rfree": [1.03], # Return factor on risk free asset (not used by this type)
189 "DiscFac": 0.96, # Intertemporal discount factor
190 "LivPrb": [0.98], # Survival probability after each period
191 "PermGroFac": [1.01], # Permanent income growth factor
192 "BoroCnstArt": 0.0, # Artificial borrowing constraint
193 "vFuncBool": False, # Whether to calculate the value function during solution
194 "CubicBool": False, # Whether to use cubic spline interpolation
195 "RiskyShareFixed": 1.0, # Fixed share of risky asset; set to None for portfolio choice
196 "ShareAugFac": 0, # Number of times to "zoom in" for an "augmented" search for optimal risky share
197 "AdjustPrb": 1.0, # Probability that the agent can update their risky portfolio share each period
198 "IndepDstnBool": True, # Whether return and income shocks are independent
199 "PortfolioBisect": False, # What does this do?
200 "pseudo_terminal": False,
201}
202IndShockRiskyAssetConsumerType_simulation_default = {
203 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
204 "AgentCount": 10000, # Number of agents of this type
205 "T_age": None, # Age after which simulated agents are automatically killed
206 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
207 # (The portion of PermGroFac attributable to aggregate productivity growth)
208 "NewbornTransShk": False, # Whether Newborns have transitory shock
209 # ADDITIONAL OPTIONAL PARAMETERS
210 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
211 # (Forces Newborns to follow solution path of the agent they replaced if True)
212 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
213 "sim_common_Rrisky": True, # Whether risky returns have a shared/common value across agents
214}
215IndShockRiskyAssetConsumerType_default = {}
216IndShockRiskyAssetConsumerType_default.update(
217 IndShockRiskyAssetConsumerType_IncShkDstn_default
218)
219IndShockRiskyAssetConsumerType_default.update(
220 IndShockRiskyAssetConsumerType_RiskyDstn_default
221)
222IndShockRiskyAssetConsumerType_default.update(
223 IndShockRiskyAssetConsumerType_aXtraGrid_default
224)
225IndShockRiskyAssetConsumerType_default.update(
226 IndShockRiskyAssetConsumerType_ShareGrid_default
227)
228IndShockRiskyAssetConsumerType_default.update(
229 IndShockRiskyAssetConsumerType_solving_default
230)
231IndShockRiskyAssetConsumerType_default.update(
232 IndShockRiskyAssetConsumerType_simulation_default
233)
234IndShockRiskyAssetConsumerType_default.update(
235 IndShockRiskyAssetConsumerType_kNrmInitDstn_default
236)
237IndShockRiskyAssetConsumerType_default.update(
238 IndShockRiskyAssetConsumerType_pLvlInitDstn_default
239)
240init_risky_asset = IndShockRiskyAssetConsumerType_default
243class IndShockRiskyAssetConsumerType(IndShockConsumerType):
244 r"""
245 A consumer type based on IndShockConsumerType, that has access to a risky asset
246 for their savings, as well as a risk-free asset. The risky asset has lognormal
247 returns that are possibly correlated with his income shocks.
249 If RiskyShareFixed is set to a number (on the unit interval), then the agent's
250 portfolio share is fixed at that value. If it is instead set to None, then the
251 agent can flexibly choose their risky asset share.
253 .. math::
254 \newcommand{\CRRA}{\rho}
255 \newcommand{\DiePrb}{\mathsf{D}}
256 \newcommand{\PermGroFac}{\Gamma}
257 \newcommand{\Rfree}{\mathsf{R}}
258 \newcommand{\DiscFac}{\beta}
259 \begin{align*}
260 v_t(m_t) &= \max_{c_t,S_t} u(c_t) + \DiscFac (1-\DiePrb_{t+1}) \mathbb{E}_{t} \left[(\PermGroFac_{t+1}\psi_{t+1})^{1-\CRRA} v_{t+1}(m_{t+1}) \right], \\
261 & \text{s.t.} \\
262 a_t &= m_t - c_t, \\
263 a_t &\geq \underline{a}, \\
264 m_{t+1} &= \mathsf{R}_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) a_t + \theta_{t+1}, \\
265 \mathsf{R}_{t+1} &=S_t\phi_{t+1}\mathbf{R}_{t+1}+ (1-S_t)\mathsf{R}_{t+1}, \\
266 (\psi_{t+1},\theta_{t+1},\phi_{t+1}) &\sim F_{t+1}, \\
267 \mathbb{E}[\psi]=\mathbb{E}[\theta] &= 1. \\
268 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA} \\
269 \end{align*}
272 Constructors
273 ------------
274 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
275 The agent's income shock distributions.
277 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
278 aXtraGrid: Constructor
279 The agent's asset grid.
281 Its default constructor is :func:`HARK.utilities.make_assets_grid`
282 ShareGrid: Constructor
283 The agent's risky asset share grid
285 Its default constructor is :func:`HARK.ConsumptionSaving.ConsRiskyAssetModel.make_simple_ShareGrid`
286 RiskyDstn: Constructor, :math:`\phi`
287 The agent's asset shock distribution for risky assets.
289 Its default constructor is :func:`HARK.Calibration.Assets.AssetProcesses.make_lognormal_RiskyDstn`
291 Solving Parameters
292 ------------------
293 cycles: int
294 0 specifies an infinite horizon model, 1 specifies a finite model.
295 T_cycle: int
296 Number of periods in the cycle for this agent type.
297 CRRA: float, :math:`\rho`
298 Coefficient of Relative Risk Aversion.
299 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
300 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
301 DiscFac: float, :math:`\beta`
302 Intertemporal discount factor.
303 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
304 Survival probability after each period.
305 PermGroFac: list[float], time varying, :math:`\Gamma`
306 Permanent income growth factor.
307 BoroCnstArt: float, default=0.0, :math:`\underline{a}`
308 The minimum Asset/Perminant Income ratio. for this agent, BoroCnstArt must be 0.
309 vFuncBool: bool
310 Whether to calculate the value function during solution.
311 CubicBool: bool
312 Whether to use cubic spline interpoliation.
314 Simulation Parameters
315 ---------------------
316 sim_common_Rrisky: Boolean
317 Whether risky returns have a shared/common value across agents. If True, Risky return's can't be time varying.
318 AgentCount: int
319 Number of agents of this kind that are created during simulations.
320 T_age: int
321 Age after which to automatically kill agents, None to ignore.
322 T_sim: int, required for simulation
323 Number of periods to simulate.
324 track_vars: list[strings]
325 List of variables that should be tracked when running the simulation.
326 For this agent, the options are 'Adjust', 'PermShk', 'Risky', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
328 Adjust is the array of which agents can adjust
330 PermShk is the agent's permanent income shock
332 Risky is the agent's risky asset shock
334 TranShk is the agent's transitory income shock
336 aLvl is the nominal asset level
338 aNrm is the normalized assets
340 bNrm is the normalized resources without this period's labor income
342 cNrm is the normalized consumption
344 mNrm is the normalized market resources
346 pLvl is the permanent income level
348 who_dies is the array of which agents died
349 aNrmInitMean: float
350 Mean of Log initial Normalized Assets.
351 aNrmInitStd: float
352 Std of Log initial Normalized Assets.
353 pLvlInitMean: float
354 Mean of Log initial permanent income.
355 pLvlInitStd: float
356 Std of Log initial permanent income.
357 PermGroFacAgg: float
358 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
359 PerfMITShk: boolean
360 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
361 NewbornTransShk: boolean
362 Whether Newborns have transitory shock.
364 Attributes
365 ----------
366 solution: list[Consumer solution object]
367 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
368 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
370 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
371 history: Dict[Array]
372 Created by running the :func:`.simulate()` method.
373 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
374 Visit :class:`HARK.core.AgentType.simulate` for more information.
375 """
377 IncShkDstn_default = IndShockRiskyAssetConsumerType_IncShkDstn_default
378 RiskyDstn_default = IndShockRiskyAssetConsumerType_RiskyDstn_default
379 aXtraGrid_default = IndShockRiskyAssetConsumerType_aXtraGrid_default
380 ShareGrid_default = IndShockRiskyAssetConsumerType_ShareGrid_default
381 solving_default = IndShockRiskyAssetConsumerType_solving_default
382 simulation_default = IndShockRiskyAssetConsumerType_simulation_default # So sphinx documents defaults
383 default_ = {
384 "params": IndShockRiskyAssetConsumerType_default,
385 "solver": NullFunc(),
386 "model": "ConsRiskyAsset.yaml",
387 "track_vars": ["aNrm", "cNrm", "Share", "mNrm", "pLvl"],
388 }
390 time_inv_ = IndShockConsumerType.time_inv_ + [
391 "PortfolioBisect",
392 "ShareGrid",
393 "IndepDstnBool",
394 "RiskyShareFixed",
395 "ShareAugFac",
396 ]
397 time_vary_ = IndShockConsumerType.time_vary_ + ["ShockDstn", "ShareLimit"]
398 shock_vars_ = IndShockConsumerType.shock_vars_ + ["Adjust", "Risky"]
399 distributions = [
400 "IncShkDstn",
401 "PermShkDstn",
402 "TranShkDstn",
403 "RiskyDstn",
404 "ShockDstn",
405 "kNrmInitDstn",
406 "pLvlInitDstn",
407 "RiskyDstn",
408 ]
410 def pre_solve(self):
411 self.construct("solution_terminal")
412 self.update_timing()
413 self.solution_terminal.ShareFunc = ConstantFunction(1.0)
415 def update_timing(self):
416 """
417 This method simply ensures that a few attributes that could be in either
418 time_inv or time_vary are appropriately labeled.
419 """
420 if type(self.AdjustDstn) is IndexDistribution:
421 self.add_to_time_vary("AdjustPrb")
422 self.del_from_time_inv("AdjustPrb")
423 else:
424 self.add_to_time_inv("AdjustPrb")
425 self.del_from_time_vary("AdjustPrb")
426 if hasattr(self.RiskyDstn, "__getitem__"):
427 self.add_to_time_vary("RiskyDstn")
428 else:
429 self.add_to_time_inv("RiskyDstn")
430 if type(self.ShareLimit) is list:
431 self.add_to_time_vary("ShareLimit")
432 self.del_from_time_inv("ShareLimit")
433 else:
434 self.add_to_time_inv("ShareLimit")
435 self.del_from_time_vary("ShareLimit")
437 def get_Rport(self):
438 """
439 Calculates realized return factor for each agent, using the attributes Rfree,
440 RiskyNow, and ShareNow.
442 Parameters
443 ----------
444 None
446 Returns
447 -------
448 Rport : np.array
449 Array of size AgentCount with each simulated agent's realized portfolio
450 return factor. Will be used by get_states() to calculate mNrmNow.
451 """
452 RfreeNow = super().get_Rport()
453 RiskyNow = self.shocks["Risky"]
454 ShareNow = self.controls["Share"]
455 Rport = ShareNow * RiskyNow + (1.0 - ShareNow) * RfreeNow
456 self.Rport = Rport
457 return Rport
459 def get_Risky(self):
460 """
461 Draws a new risky return factor.
463 Parameters
464 ----------
465 None
467 Returns
468 -------
469 None
470 """
472 # How we draw the shocks depends on whether their distribution is time-varying
473 if "RiskyDstn" in self.time_vary:
474 if self.sim_common_Rrisky:
475 raise AttributeError(
476 "If sim_common_Rrisky is True, RiskyDstn cannot be time-varying!"
477 )
479 else:
480 # Make use of the IndexDistribution.draw() method
481 self.shocks["Risky"] = self.RiskyDstn.draw(
482 np.maximum(self.t_cycle - 1, 0)
483 if self.cycles == 1
484 else self.t_cycle
485 )
487 else:
488 # Draw either a common economy-wide return, or one for each agent
489 if self.sim_common_Rrisky:
490 self.shocks["Risky"] = self.RiskyDstn.draw(1)
491 else:
492 self.shocks["Risky"] = self.RiskyDstn.draw(self.AgentCount)
494 def get_Adjust(self):
495 """
496 Sets the attribute Adjust as a boolean array of size AgentCount, indicating
497 whether each agent is able to adjust their risky portfolio share this period.
498 Uses the attribute AdjustPrb to draw from a Bernoulli distribution.
500 Parameters
501 ----------
502 None
504 Returns
505 -------
506 None
507 """
508 if "AdjustPrb" in self.time_vary:
509 self.shocks["Adjust"] = self.AdjustDstn.draw(
510 np.maximum(self.t_cycle - 1, 0) if self.cycles == 1 else self.t_cycle
511 )
512 else:
513 self.shocks["Adjust"] = self.AdjustDstn.draw(self.AgentCount)
515 def initialize_sim(self):
516 """
517 Initialize the state of simulation attributes. Simply calls the same
518 method for IndShockConsumerType, then initializes the new states/shocks
519 Adjust and Share.
521 Parameters
522 ----------
523 None
525 Returns
526 -------
527 None
528 """
529 self.shocks["Adjust"] = np.zeros(self.AgentCount, dtype=bool)
530 # Initialize Share to default value; will be updated in get_controls()
531 self.controls["Share"] = np.ones(self.AgentCount)
532 IndShockConsumerType.initialize_sim(self)
534 def get_shocks(self):
535 """
536 Draw idiosyncratic income shocks, just as for IndShockConsumerType, then draw
537 a single common value for the risky asset return. Also draws whether each
538 agent is able to adjust their portfolio this period.
540 Parameters
541 ----------
542 None
544 Returns
545 -------
546 None
547 """
548 IndShockConsumerType.get_shocks(self)
549 self.get_Risky()
550 self.get_Adjust()
552 def get_controls(self):
553 """
554 Calculates consumption for each consumer of this type using the consumption functions;
555 also calculates risky asset share when there is portfolio share.
557 Parameters
558 ----------
559 None
561 Returns
562 -------
563 None
564 """
565 cNrmNow = np.full(self.AgentCount, np.nan)
566 MPCnow = np.full(self.AgentCount, np.nan)
567 ShareNow = np.full(self.AgentCount, np.nan)
568 for t in np.unique(self.t_cycle):
569 idx = self.t_cycle == t
570 if np.any(idx):
571 mNrm = self.state_now["mNrm"][idx]
572 cNrmNow[idx], MPCnow[idx] = self.solution[t].cFunc.eval_with_derivative(
573 mNrm
574 )
575 if self.RiskyShareFixed is None:
576 ShareNow[idx] = self.solution[t].ShareFunc(mNrm)
577 else:
578 ShareNow[idx] = self.RiskyShareFixed
579 self.controls["cNrm"] = cNrmNow
580 self.controls["Share"] = ShareNow
581 self.MPCnow = MPCnow
583 def check_conditions(self, verbose=None):
584 raise NotImplementedError() # pragma: nocover
586 def calc_limiting_values(self):
587 raise NotImplementedError() # pragma: nocover
590# This is to preserve compatibility with old name
591RiskyAssetConsumerType = IndShockRiskyAssetConsumerType
594###############################################################################
597def solve_one_period_ConsPortChoice(
598 solution_next,
599 ShockDstn,
600 IncShkDstn,
601 RiskyDstn,
602 LivPrb,
603 DiscFac,
604 CRRA,
605 Rfree,
606 PermGroFac,
607 BoroCnstArt,
608 aXtraGrid,
609 ShareGrid,
610 ShareLimit,
611 ShareAugFac,
612 vFuncBool,
613 IndepDstnBool,
614):
615 """
616 Solve one period of a consumption-saving problem with portfolio allocation
617 between a riskless and risky asset. This function handles only the most
618 fundamental portfolio choice problem: frictionless reallocation of the
619 portfolio each period as a continuous choice.
621 Parameters
622 ----------
623 solution_next : PortfolioSolution
624 Solution to next period's problem.
625 ShockDstn : Distribution
626 Joint distribution of permanent income shocks, transitory income shocks,
627 and risky returns. This is only used if the input IndepDstnBool is False,
628 indicating that income and return distributions can't be assumed to be
629 independent.
630 IncShkDstn : Distribution
631 Discrete distribution of permanent income shocks and transitory income
632 shocks. This is only used if the input IndepDstnBool is True, indicating
633 that income and return distributions are independent.
634 RiskyDstn : Distribution
635 Distribution of risky asset returns. This is only used if the input
636 IndepDstnBool is True, indicating that income and return distributions
637 are independent.
638 LivPrb : float
639 Survival probability; likelihood of being alive at the beginning of
640 the succeeding period.
641 DiscFac : float
642 Intertemporal discount factor for future utility.
643 CRRA : float
644 Coefficient of relative risk aversion.
645 Rfree : float
646 Risk free interest factor on end-of-period assets.
647 PermGroFac : float
648 Expected permanent income growth factor at the end of this period.
649 BoroCnstArt: float or None
650 Borrowing constraint for the minimum allowable assets to end the
651 period with. In this model, it is *required* to be zero.
652 aXtraGrid: np.array
653 Array of "extra" end-of-period asset values-- assets above the
654 absolute minimum acceptable level.
655 ShareGrid : np.array
656 Array of risky portfolio shares on which to define the interpolation
657 of the consumption function when Share is fixed. Also used when the
658 risky share choice is specified as discrete rather than continuous.
659 ShareLimit : float
660 Limiting lower bound of risky portfolio share as mNrm approaches infinity.
661 ShareAugFac : int
662 Number of times to perform an "augmented" search for optimal risky asset
663 shares by "zooming in" on on FOC-crossing region. Setting this above zero
664 will make the solution slightly more accurate and can aid stability for
665 "unusual" or extreme parameter sets.
666 vFuncBool: boolean
667 An indicator for whether the value function should be computed and
668 included in the reported solution.
669 IndepDstnBool : bool
670 Indicator for whether the income and risky return distributions are in-
671 dependent of each other, which can speed up the expectations step.
673 Returns
674 -------
675 solution_now : PortfolioSolution
676 Solution to this period's problem.
678 :meta private:
679 """
680 # Make sure the individual is liquidity constrained. Allowing a consumer to
681 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
682 if BoroCnstArt != 0.0:
683 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
685 # Define the current period utility function and effective discount factor
686 uFunc = UtilityFuncCRRA(CRRA)
687 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
689 # Unpack next period's solution for easier access
690 vPfunc_next = solution_next.vPfunc
691 vFunc_next = solution_next.vFunc
693 # Set a flag for whether the natural borrowing constraint is zero, which
694 # depends on whether the smallest transitory income shock is zero
695 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
697 # Prepare to calculate end-of-period marginal values by creating an array
698 # of market resources that the agent could have next period, considering
699 # the grid of end-of-period assets and the distribution of shocks he might
700 # experience next period.
702 # Unpack the risky return shock distribution
703 Risky_next = RiskyDstn.atoms
704 RiskyMax = np.max(Risky_next)
705 RiskyMin = np.min(Risky_next)
707 # Perform an alternate calculation of the absolute patience factor when
708 # returns are risky. This uses the Merton-Samuelson limiting risky share,
709 # which is what's relevant as mNrm goes to infinity.
710 def calc_Radj(R):
711 Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree
712 return Rport ** (1.0 - CRRA)
714 R_adj = expected(calc_Radj, RiskyDstn)[0]
715 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
716 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
718 # Also perform an alternate calculation for human wealth under risky returns
719 def calc_hNrm(S):
720 Risky = S["Risky"]
721 PermShk = S["PermShk"]
722 TranShk = S["TranShk"]
723 G = PermGroFac * PermShk
724 Rport = ShareLimit * Risky + (1.0 - ShareLimit) * Rfree
725 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
726 return hNrm
728 # This correctly accounts for risky returns and risk aversion
729 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
731 # This basic equation works if there's no correlation among shocks
732 # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm)
734 # Define the terms for the limiting linear consumption function as m gets very big
735 cFuncLimitIntercept = MPCminNow * hNrmNow
736 cFuncLimitSlope = MPCminNow
738 # bNrm represents R*a, balances after asset return shocks but before income.
739 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
740 if BoroCnstNat_iszero:
741 aNrmGrid = aXtraGrid
742 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
743 else:
744 # Add an asset point at exactly zero
745 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
746 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
748 # Get grid and shock sizes, for easier indexing
749 aNrmCount = aNrmGrid.size
750 ShareCount = ShareGrid.size
752 # If the income shock distribution is independent from the risky return distribution,
753 # then taking end-of-period expectations can proceed in a two part process: First,
754 # construct an "intermediate" value function by integrating out next period's income
755 # shocks, *then* compute end-of-period expectations by integrating out return shocks.
756 # This method is lengthy to code, but can be significantly faster.
757 if IndepDstnBool:
758 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
759 bNrmNext = bNrmGrid
761 # Define functions that are used internally to evaluate future realizations
762 def calc_mNrm_next(S, b):
763 """
764 Calculate future realizations of market resources mNrm from the income
765 shock distribution S and normalized bank balances b.
766 """
767 return b / (S["PermShk"] * PermGroFac) + S["TranShk"]
769 def calc_dvdm_next(S, b):
770 """
771 Evaluate realizations of marginal value of market resources next period,
772 based on the income distribution S and values of bank balances bNrm
773 """
774 mNrm_next = calc_mNrm_next(S, b)
775 G = S["PermShk"] * PermGroFac
776 dvdm_next = G ** (-CRRA) * vPfunc_next(mNrm_next)
777 return dvdm_next
779 # Calculate end-of-period marginal value of assets and shares at each point
780 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
781 # values across income and risky return shocks.
783 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
784 dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext))
785 dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0))
787 dvdbNvrsFunc_intermed = LinearInterp(bNrmGrid, dvdbNvrs_intermed)
788 dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA)
790 # The intermediate marginal value of risky portfolio share is zero in this
791 # model because risky share is flexible: we can just change it next period,
792 # so there is no marginal value of Share once the return is realized.
793 dvdsFunc_intermed = ConstantFunction(0.0) # all zeros
795 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
796 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
798 # Define functions for calculating end-of-period marginal value
799 def calc_EndOfPrd_dvda(R, a, z):
800 """
801 Compute end-of-period marginal value of assets at values a, conditional
802 on risky asset return R and risky share z.
803 """
804 # Calculate future realizations of bank balances bNrm
805 Rxs = R - Rfree # Excess returns
806 Rport = Rfree + z * Rxs # Portfolio return
807 bNrm_next = Rport * a
809 # Calculate and return dvda
810 EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next)
811 return EndOfPrd_dvda
813 def calc_EndOfPrd_dvds(R, a, z):
814 """
815 Compute end-of-period marginal value of risky share at values a, conditional
816 on risky asset return S and risky share z.
817 """
818 # Calculate future realizations of bank balances bNrm
819 Rxs = R - Rfree # Excess returns
820 Rport = Rfree + z * Rxs # Portfolio return
821 bNrm_next = Rport * a
823 # Calculate and return dvds (second term is all zeros)
824 EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(
825 bNrm_next
826 )
827 return EndOfPrd_dvds
829 TempDstn = RiskyDstn # relabeling for below
831 # Evaluate realizations of value and marginal value after asset returns are realized
833 # Calculate end-of-period marginal value of risky portfolio share by taking expectations
834 EndOfPrd_dvds = DiscFacEff * expected(
835 calc_EndOfPrd_dvds, RiskyDstn, args=(aNrmNow, ShareNext)
836 )
838 # Make the end-of-period value function if the value function is requested
839 if vFuncBool:
841 def calc_v_intermed(S, b):
842 """
843 Calculate "intermediate" value from next period's bank balances, the
844 income shocks S, and the risky asset share.
845 """
846 mNrm_next = calc_mNrm_next(S, b)
847 v_next = vFunc_next(mNrm_next)
848 v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
849 return v_intermed
851 # Calculate intermediate value by taking expectations over income shocks
852 v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext))
854 # Construct the "intermediate value function" for this period
855 vNvrs_intermed = uFunc.inv(v_intermed)
856 vNvrsFunc_intermed = LinearInterp(bNrmGrid, vNvrs_intermed)
857 vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA)
859 def calc_EndOfPrd_v(S, a, z):
860 # Calculate future realizations of bank balances bNrm
861 Rxs = S - Rfree
862 Rport = Rfree + z * Rxs
863 bNrm_next = Rport * a
865 EndOfPrd_v = vFunc_intermed(bNrm_next)
866 return EndOfPrd_v
868 # Calculate end-of-period value by taking expectations
869 EndOfPrd_v = DiscFacEff * expected(
870 calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext)
871 )
872 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
874 # Now make an end-of-period value function over aNrm and Share
875 EndOfPrd_vNvrsFunc = BilinearInterp(EndOfPrd_vNvrs, aNrmGrid, ShareGrid)
876 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
877 # This will be used later to make the value function for this period
879 # If the income shock distribution and risky return distribution are *NOT*
880 # independent, then computation of end-of-period expectations are simpler in
881 # code, but might take longer to execute
882 else:
883 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
884 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
886 # Define functions that are used internally to evaluate future realizations
887 def calc_mNrm_next(S, a, z):
888 """
889 Calculate future realizations of market resources mNrm from the shock
890 distribution S, normalized end-of-period assets a, and risky share z.
891 """
892 # Calculate future realizations of bank balances bNrm
893 Rxs = S["Risky"] - Rfree
894 Rport = Rfree + z * Rxs
895 bNrm_next = Rport * a
896 mNrm_next = bNrm_next / (S["PermShk"] * PermGroFac) + S["TranShk"]
897 return mNrm_next
899 def calc_EndOfPrd_dvdx(S, a, z):
900 """
901 Evaluate end-of-period marginal value of assets and risky share based
902 on the shock distribution S, normalized end-of-period assets a, and
903 risky share z.
904 """
905 mNrm_next = calc_mNrm_next(S, a, z)
906 Rxs = S["Risky"] - Rfree
907 Rport = Rfree + z * Rxs
908 dvdm_next = vPfunc_next(mNrm_next)
909 # No marginal value of Share if it's a free choice!
910 dvds_next = np.zeros_like(mNrm_next)
912 EndOfPrd_dvda = Rport * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
913 EndOfPrd_dvds = (
914 Rxs * a * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
915 + (S["PermShk"] * PermGroFac) ** (1 - CRRA) * dvds_next
916 )
918 return EndOfPrd_dvda, EndOfPrd_dvds
920 def calc_EndOfPrd_v(S, a, z):
921 """
922 Evaluate end-of-period value, based on the shock distribution S, values
923 of bank balances bNrm, and values of the risky share z.
924 """
925 mNrm_next = calc_mNrm_next(S, a, z)
926 v_next = vFunc_next(mNrm_next)
927 EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
928 return EndOfPrd_v
930 calc_EndOfPrd_dvda = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[0]
931 calc_EndOfPrd_dvds = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[1]
932 TempDstn = ShockDstn
934 # Evaluate realizations of value and marginal value after asset returns are realized
936 # Calculate end-of-period marginal value of assets and risky share by taking expectations
937 EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected(
938 calc_EndOfPrd_dvdx, ShockDstn, args=(aNrmNow, ShareNext)
939 )
940 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
942 # Construct the end-of-period value function if requested
943 if vFuncBool:
944 # Calculate end-of-period value, its derivative, and their pseudo-inverse
945 EndOfPrd_v = DiscFacEff * expected(
946 calc_EndOfPrd_v, ShockDstn, args=(aNrmNow, ShareNext)
947 )
948 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
950 # value transformed through inverse utility
951 EndOfPrd_vNvrsP = EndOfPrd_dvda * uFunc.derinv(EndOfPrd_v, order=(0, 1))
953 # Construct the end-of-period value function
954 EndOfPrd_vNvrsFunc_by_Share = []
955 for j in range(ShareCount):
956 EndOfPrd_vNvrsFunc_by_Share.append(
957 CubicInterp(
958 aNrmNow[:, j], EndOfPrd_vNvrs[:, j], EndOfPrd_vNvrsP[:, j]
959 )
960 )
961 EndOfPrd_vNvrsFunc = LinearInterpOnInterp1D(
962 EndOfPrd_vNvrsFunc_by_Share, ShareGrid
963 )
964 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
966 # Now find the optimal (continuous) risky share on [0,1] by solving the first
967 # order condition EndOfPrd_dvds == 0.
968 FOC_s = EndOfPrd_dvds # Relabel for convenient typing
970 # If agent wants to put more than 100% into risky asset, he is constrained.
971 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
972 constrained_top = FOC_s[:, -1] > 0.0
973 constrained_bot = FOC_s[:, 0] < 0.0
974 constrained = np.logical_or(constrained_top, constrained_bot)
975 a_idx = np.arange(aNrmCount)
977 # For each value of aNrm, find the value of Share such that FOC_s == 0
978 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0)
979 share_idx = np.argmax(crossing, axis=1)
981 # share_idx represents the index of the segment of the share grid where dvds flips
982 # from positive to negative, indicating that there's a zero *on* the segment.
983 # The exception is for aNrm values that are flagged as constrained, for which
984 # there will be no crossing point and we can just use the boundary value.
986 # Now that we have a *range* for the location of the optimal share, we can
987 # do a refined search for the optimal share at each aNrm value where there
988 # is an interior solution (not constrained). We now make a refined ShareGrid
989 # that has *different* values on it for each aNrm value.
991 for k in range(ShareAugFac):
992 bot_s = ShareNext[a_idx, share_idx]
993 top_s = ShareNext[a_idx, share_idx + 1]
994 for j in range(aNrmCount):
995 if constrained[j]:
996 continue
997 ShareNext[j, :] = np.linspace(bot_s[j], top_s[j], ShareCount)
999 # Now evaluate end-of-period marginal value of risky share on the refined grid
1000 EndOfPrd_dvds = DiscFacEff * expected(
1001 calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)
1002 )
1003 these = np.logical_not(constrained)
1004 FOC_s[these, :] = EndOfPrd_dvds[these, :] # Relabel for convenient typing
1006 # Look for "crossing points" again
1007 crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0)
1008 share_idx[these] = np.argmax(crossing, axis=1)
1010 # Calculate end-of-period marginal value of assets on the refined grid
1011 EndOfPrd_dvda = DiscFacEff * expected(
1012 calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)
1013 )
1014 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
1016 # Calculate the fractional distance between those share gridpoints where the
1017 # zero should be found, assuming a linear function; call it alpha
1018 bot_s = ShareNext[a_idx, share_idx]
1019 top_s = ShareNext[a_idx, share_idx + 1]
1020 bot_f = FOC_s[a_idx, share_idx]
1021 top_f = FOC_s[a_idx, share_idx + 1]
1022 bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx]
1023 top_c = EndOfPrd_dvdaNvrs[a_idx, share_idx + 1]
1024 alpha = 1.0 - top_f / (top_f - bot_f)
1026 # Calculate the continuous optimal risky share and optimal consumption
1027 Share_now = (1.0 - alpha) * bot_s + alpha * top_s
1028 cNrm_now = (1.0 - alpha) * bot_c + alpha * top_c
1030 # If agent wants to put more than 100% into risky asset, he is constrained.
1031 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
1032 constrained_top = FOC_s[:, -1] > 0.0
1033 constrained_bot = FOC_s[:, 0] < 0.0
1035 # Apply the constraints to both risky share and consumption (but lower
1036 # constraint should never be relevant)
1037 Share_now[constrained_top] = 1.0
1038 Share_now[constrained_bot] = 0.0
1039 cNrm_now[constrained_top] = EndOfPrd_dvdaNvrs[constrained_top, -1]
1040 cNrm_now[constrained_bot] = EndOfPrd_dvdaNvrs[constrained_bot, 0]
1042 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
1043 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
1044 # can't depend on the risky share if it doesn't meaningfully exist. Apply
1045 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
1046 if not BoroCnstNat_iszero:
1047 Share_now[0] = 1.0
1048 cNrm_now[0] = EndOfPrd_dvdaNvrs[0, -1]
1050 # Construct functions characterizing the solution for this period
1052 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
1053 # then construct the consumption function when the agent can adjust his share
1054 mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0)
1055 cNrm_now = np.insert(cNrm_now, 0, 0.0)
1056 cFunc_now = LinearInterp(mNrm_now, cNrm_now, cFuncLimitIntercept, cFuncLimitSlope)
1058 # Construct the marginal value (of mNrm) function
1059 vPfunc_now = MargValueFuncCRRA(cFunc_now, CRRA)
1061 # If the share choice is continuous, just make an ordinary interpolating function
1062 if BoroCnstNat_iszero:
1063 Share_lower_bound = ShareLimit
1064 else:
1065 Share_lower_bound = 1.0
1066 Share_now = np.insert(Share_now, 0, Share_lower_bound)
1067 ShareFunc_now = LinearInterp(mNrm_now, Share_now, ShareLimit, 0.0)
1069 # Add the value function if requested
1070 if vFuncBool:
1071 # Create the value functions for this period, defined over market resources
1072 # mNrm when agent can adjust his portfolio, and over market resources and
1073 # fixed share when agent can not adjust his portfolio.
1075 # Construct the value function
1076 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
1077 cNrm_temp = cFunc_now(mNrm_temp)
1078 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # Fix tiny violations
1079 Share_temp = ShareFunc_now(mNrm_temp)
1080 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
1081 vNvrs_temp = uFunc.inv(v_temp)
1082 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
1083 vNvrsFunc = CubicInterp(
1084 np.insert(mNrm_temp, 0, 0.0), # x_list
1085 np.insert(vNvrs_temp, 0, 0.0), # f_list
1086 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
1087 )
1088 # Re-curve the pseudo-inverse value function
1089 vFunc_now = ValueFuncCRRA(vNvrsFunc, CRRA)
1091 else: # If vFuncBool is False, fill in dummy values
1092 vFunc_now = NullFunc()
1094 # Package and return the solution
1095 solution_now = ConsumerSolution(
1096 cFunc=cFunc_now,
1097 vPfunc=vPfunc_now,
1098 vFunc=vFunc_now,
1099 hNrm=hNrmNow,
1100 MPCmin=MPCminNow,
1101 )
1102 solution_now.ShareFunc = ShareFunc_now
1103 return solution_now
1106###############################################################################
1109def solve_one_period_ConsIndShockRiskyAsset(
1110 solution_next,
1111 IncShkDstn,
1112 RiskyDstn,
1113 ShockDstn,
1114 LivPrb,
1115 DiscFac,
1116 Rfree,
1117 CRRA,
1118 PermGroFac,
1119 BoroCnstArt,
1120 aXtraGrid,
1121 RiskyShareFixed,
1122 vFuncBool,
1123 CubicBool,
1124 IndepDstnBool,
1125):
1126 """
1127 Solves one period of a consumption-saving model with idiosyncratic shocks to
1128 permanent and transitory income, with a risky and riskless asset and CRRA utility.
1130 Parameters
1131 ----------
1132 solution_next : ConsumerSolution
1133 The solution to next period's one period problem.
1134 IncShkDstn : Distribution
1135 Discrete distribution of permanent income shocks and transitory income
1136 shocks. This is only used if the input IndepDstnBool is True, indicating
1137 that income and return distributions are independent.
1138 RiskyDstn : Distribution
1139 Distribution of risky asset returns. This is only used if the input
1140 IndepDstnBool is True, indicating that income and return distributions
1141 are independent.
1142 ShockDstn : Distribution
1143 Joint distribution of permanent income shocks, transitory income shocks,
1144 and risky returns. This is only used if the input IndepDstnBool is False,
1145 indicating that income and return distributions can't be assumed to be
1146 independent.
1147 LivPrb : float
1148 Survival probability; likelihood of being alive at the beginning of
1149 the succeeding period.
1150 DiscFac : float
1151 Intertemporal discount factor for future utility.
1152 Rfree : float
1153 Risk free interest factor on end-of-period assets.
1154 CRRA : float
1155 Coefficient of relative risk aversion.
1156 PermGroFac : float
1157 Expected permanent income growth factor at the end of this period.
1158 BoroCnstArt: float or None
1159 Borrowing constraint for the minimum allowable assets to end the
1160 period with. If it is less than the natural borrowing constraint,
1161 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
1162 rowing constraint.
1163 aXtraGrid: np.array
1164 Array of "extra" end-of-period asset values-- assets above the
1165 absolute minimum acceptable level.
1166 RiskyShareFixed : float
1167 Fixed fraction of end-of-period assets that are allocated to the risky asset.
1168 vFuncBool: boolean
1169 An indicator for whether the value function should be computed and
1170 included in the reported solution.
1171 CubicBool: boolean
1172 An indicator for whether the solver should use cubic or linear interpolation.
1173 IndepDstnBool : bool
1174 Indicator for whether the income and risky return distributions are in-
1175 dependent of each other, which can speed up the expectations step.
1177 Returns
1178 -------
1179 solution_now : ConsumerSolution
1180 Solution to this period's consumption-saving problem with income risk.
1182 :meta private:
1183 """
1184 # Do a quick validity check; don't want to allow borrowing with risky returns
1185 if BoroCnstArt != 0.0:
1186 raise ValueError("RiskyAssetConsumerType must have BoroCnstArt=0.0!")
1188 # Define the current period utility function and effective discount factor
1189 uFunc = UtilityFuncCRRA(CRRA)
1190 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
1192 # Unpack next period's income shock distribution
1193 ShkPrbsNext = ShockDstn.pmv
1194 PermShkValsNext = ShockDstn.atoms[0]
1195 TranShkValsNext = ShockDstn.atoms[1]
1196 RiskyValsNext = ShockDstn.atoms[2]
1197 PermShkMinNext = np.min(PermShkValsNext)
1198 TranShkMinNext = np.min(TranShkValsNext)
1199 RiskyMinNext = np.min(RiskyValsNext)
1200 RiskyMaxNext = np.max(RiskyValsNext)
1202 # Unpack next period's (marginal) value function
1203 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
1204 vPfuncNext = solution_next.vPfunc
1205 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
1207 # Perform an alternate calculation of the absolute patience factor when returns are risky
1208 def calc_Radj(R):
1209 Rport = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree
1210 return Rport ** (1.0 - CRRA)
1212 R_adj = expected(calc_Radj, RiskyDstn)
1213 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
1214 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
1215 MPCminNow = MPCminNow[0]
1217 # Also perform an alternate calculation for human wealth under risky returns
1218 def calc_hNrm(S):
1219 Risky = S["Risky"]
1220 PermShk = S["PermShk"]
1221 TranShk = S["TranShk"]
1222 G = PermGroFac * PermShk
1223 Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree
1224 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
1225 return hNrm
1227 # This correctly accounts for risky returns and risk aversion
1228 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
1229 hNrmNow = hNrmNow[0]
1231 # The above attempts to pin down the limiting consumption function for this
1232 # model, however it is not clear why it creates bugs, so for now we allow
1233 # for a linear extrapolation beyond the last asset point
1234 cFuncLimitIntercept = MPCminNow * hNrmNow
1235 cFuncLimitSlope = MPCminNow
1237 # Calculate the minimum allowable value of market resources in this period
1238 BoroCnstNat_cand = (
1239 (solution_next.mNrmMin - TranShkValsNext)
1240 * (PermGroFac * PermShkValsNext)
1241 / RiskyValsNext
1242 )
1243 BoroCnstNat = np.max(BoroCnstNat_cand) # Must be at least this
1245 # Set a flag for whether the natural borrowing constraint is zero, which
1246 # depends on whether the smallest transitory income shock is zero
1247 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
1249 # Set the minimum allowable (normalized) market resources based on the natural
1250 # and artificial borrowing constraints
1251 mNrmMinNow = (
1252 np.max([BoroCnstNat, BoroCnstArt]) if BoroCnstArt is not None else BoroCnstNat
1253 )
1255 # The MPCmax code is a bit unusual here, and possibly "harmlessly wrong".
1256 # The "worst event" should depend on the risky return factor as well as
1257 # income shocks. However, the natural borrowing constraint is only ever
1258 # relevant in this model when it's zero, so the MPC at mNrm is only relevant
1259 # in the case where risky returns don't matter at all (because a=0).
1261 # Calculate the probability that we get the worst possible income draw
1262 IncNext = PermShkValsNext * TranShkValsNext
1263 WorstIncNext = PermShkMinNext * TranShkMinNext
1264 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
1265 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
1267 # Update the upper bounding MPC as market resources approach the lower bound
1268 temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
1269 MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
1271 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
1272 # or artificial borrowing constraint actually binds
1273 if BoroCnstNat < mNrmMinNow:
1274 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
1275 else:
1276 MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above
1278 # Define the borrowing-constrained consumption function
1279 cFuncNowCnst = LinearInterp(
1280 np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0])
1281 )
1283 # Big methodological split here: whether the income and return distributions are independent.
1284 # Calculation of end-of-period marginal (marginal) value uses different approaches
1285 if IndepDstnBool:
1286 # bNrm represents R*a, balances after asset return shocks but before income.
1287 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
1288 if BoroCnstNat_iszero:
1289 bNrmNow = np.insert(
1290 RiskyMaxNext * aXtraGrid, 0, RiskyMinNext * aXtraGrid[0]
1291 )
1292 aNrmNow = aXtraGrid.copy()
1293 else:
1294 # Add a bank balances point at exactly zero
1295 bNrmNow = RiskyMaxNext * np.insert(aXtraGrid, 0, 0.0)
1296 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1298 # Define local functions for taking future expectations when the interest
1299 # factor *is* independent from the income shock distribution. These go
1300 # from "bank balances" bNrm = R * aNrm to t+1 realizations.
1301 def calc_mNrmNext(S, b):
1302 return b / (PermGroFac * S["PermShk"]) + S["TranShk"]
1304 def calc_vNext(S, b):
1305 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, b))
1307 def calc_vPnext(S, b):
1308 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, b))
1310 def calc_vPPnext(S, b):
1311 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, b))
1313 # Calculate marginal value of bank balances at each gridpoint
1314 vPfacEff = PermGroFac ** (-CRRA)
1315 Intermed_vP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(bNrmNow))
1316 Intermed_vPnvrs = uFunc.derinv(Intermed_vP, order=(1, 0))
1318 if BoroCnstNat_iszero:
1319 Intermed_vPnvrs = np.insert(Intermed_vPnvrs, 0, 0.0)
1320 bNrm_temp = np.insert(bNrmNow, 0, 0.0)
1321 else:
1322 bNrm_temp = bNrmNow.copy()
1324 # If using cubic spline interpolation, also calculate "intermediate"
1325 # marginal marginal value of bank balances
1326 if CubicBool:
1327 vPPfacEff = PermGroFac ** (-CRRA - 1.0)
1328 Intermed_vPP = vPPfacEff * expected(
1329 calc_vPPnext, IncShkDstn, args=(bNrmNow)
1330 )
1331 Intermed_vPnvrsP = Intermed_vPP * uFunc.derinv(Intermed_vP, order=(1, 1))
1332 if BoroCnstNat_iszero:
1333 Intermed_vPnvrsP = np.insert(Intermed_vPnvrsP, 0, Intermed_vPnvrsP[0])
1335 # Make a cubic spline intermediate pseudo-inverse marginal value function
1336 Intermed_vPnvrsFunc = CubicInterp(
1337 bNrm_temp,
1338 Intermed_vPnvrs,
1339 Intermed_vPnvrsP,
1340 lower_extrap=True,
1341 )
1342 Intermed_vPPfunc = MargMargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1343 else:
1344 # Make a linear interpolation intermediate pseudo-inverse marginal value function
1345 Intermed_vPnvrsFunc = LinearInterp(
1346 bNrm_temp, Intermed_vPnvrs, lower_extrap=True
1347 )
1349 # "Recurve" the intermediate pseudo-inverse marginal value function
1350 Intermed_vPfunc = MargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1352 # If the value function is requested, calculate "intermediate" value
1353 if vFuncBool:
1354 vFacEff = PermGroFac ** (1.0 - CRRA)
1355 Intermed_v = vFacEff * expected(calc_vNext, IncShkDstn, args=(bNrmNow))
1356 Intermed_vNvrs = uFunc.inv(Intermed_v)
1357 # value transformed through inverse utility
1358 Intermed_vNvrsP = Intermed_vP * uFunc.derinv(Intermed_v, order=(0, 1))
1359 if BoroCnstNat_iszero:
1360 Intermed_vNvrs = np.insert(Intermed_vNvrs, 0, 0.0)
1361 Intermed_vNvrsP = np.insert(Intermed_vNvrsP, 0, Intermed_vNvrsP[0])
1362 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1364 # Make a cubic spline intermediate pseudo-inverse value function
1365 Intermed_vNvrsFunc = CubicInterp(bNrm_temp, Intermed_vNvrs, Intermed_vNvrsP)
1367 # "Recurve" the intermediate pseudo-inverse value function
1368 Intermed_vFunc = ValueFuncCRRA(Intermed_vNvrsFunc, CRRA)
1370 # We have "intermediate" (marginal) value functions defined over bNrm,
1371 # so now we want to take expectations over Risky realizations at each aNrm.
1373 # Begin by re-defining transition functions for taking expectations, which are all very simple!
1374 Z = RiskyShareFixed # for shorter notation
1376 def calc_bNrmNext(R, a):
1377 Rport = Z * R + (1 - Z) * Rfree
1378 return Rport * a
1380 def calc_vNext(R, a):
1381 return Intermed_vFunc(calc_bNrmNext(R, a))
1383 def calc_vPnext(R, a):
1384 Rport = Z * R + (1 - Z) * Rfree
1385 return Rport * Intermed_vPfunc(calc_bNrmNext(R, a))
1387 def calc_vPPnext(R, a):
1388 Rport = Z * R + (1 - Z) * Rfree
1389 return Rport * Rport * Intermed_vPPfunc(calc_bNrmNext(R, a))
1391 # Calculate end-of-period marginal value of assets at each gridpoint
1392 EndOfPrdvP = DiscFacEff * expected(calc_vPnext, RiskyDstn, args=(aNrmNow))
1394 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1395 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1396 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1398 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1399 if CubicBool:
1400 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1401 EndOfPrdvPP = DiscFacEff * expected(calc_vPPnext, RiskyDstn, args=(aNrmNow))
1402 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1403 MPC = dcda / (dcda + 1.0)
1404 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1406 # Limiting consumption is zero as m approaches mNrmMin
1407 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1408 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1410 # Construct the end-of-period value function if requested
1411 if vFuncBool:
1412 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1413 EndOfPrdv = DiscFacEff * expected(calc_vNext, RiskyDstn, args=(aNrmNow))
1414 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1415 # value transformed through inverse utility
1416 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1418 # Construct the end-of-period value function
1419 if BoroCnstNat_iszero:
1420 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1421 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1422 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1423 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1424 else:
1425 aNrm_temp = aNrmNow.copy()
1427 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1428 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1430 # NON-INDEPENDENT METHOD BEGINS HERE
1431 else:
1432 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
1433 # aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
1434 if BoroCnstNat_iszero:
1435 aNrmNow = aXtraGrid
1436 else:
1437 # Add an asset point at exactly zero
1438 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1440 # Define local functions for taking future expectations when the interest
1441 # factor is *not* independent from the income shock distribution
1442 Z = RiskyShareFixed # for shorter notation
1444 def calc_mNrmNext(S, a):
1445 Risky = S["Risky"]
1446 Rport = Z * Risky + (1 - Z) * Rfree
1447 return Rport / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
1449 def calc_vNext(S, a):
1450 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, a))
1452 def calc_vPnext(S, a):
1453 Risky = S["Risky"]
1454 Rport = Z * Risky + (1 - Z) * Rfree
1455 return Rport * S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a))
1457 def calc_vPPnext(S, a):
1458 Risky = S["Risky"]
1459 Rport = Z * Risky + (1 - Z) * Rfree
1460 return (
1461 (Rport**2)
1462 * S["PermShk"] ** (-CRRA - 1.0)
1463 * vPPfuncNext(calc_mNrmNext(S, a))
1464 )
1466 # Calculate end-of-period marginal value of assets at each gridpoint
1467 vPfacEff = DiscFacEff * PermGroFac ** (-CRRA)
1468 EndOfPrdvP = vPfacEff * expected(calc_vPnext, ShockDstn, args=(aNrmNow))
1470 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1471 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1472 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1474 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1475 if CubicBool:
1476 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1477 vPPfacEff = DiscFacEff * PermGroFac ** (-CRRA - 1.0)
1478 EndOfPrdvPP = vPPfacEff * expected(calc_vPPnext, ShockDstn, args=(aNrmNow))
1479 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1480 MPC = dcda / (dcda + 1.0)
1481 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1483 # Limiting consumption is zero as m approaches mNrmMin
1484 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1485 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1487 # Construct the end-of-period value function if requested
1488 if vFuncBool:
1489 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1490 vFacEff = DiscFacEff * PermGroFac ** (1.0 - CRRA)
1491 EndOfPrdv = vFacEff * expected(calc_vNext, ShockDstn, args=(aNrmNow))
1492 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1493 # value transformed through inverse utility
1494 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1496 # Construct the end-of-period value function
1497 if BoroCnstNat_iszero:
1498 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1499 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1500 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1501 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1502 else:
1503 aNrm_temp = aNrmNow.copy()
1505 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1506 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1508 # Construct the consumption function; this uses the same method whether the
1509 # income distribution is independent from the return distribution
1510 if CubicBool:
1511 # Construct the unconstrained consumption function as a cubic interpolation
1512 cFuncNowUnc = CubicInterp(
1513 m_for_interpolation,
1514 c_for_interpolation,
1515 MPC_for_interpolation,
1516 cFuncLimitIntercept,
1517 cFuncLimitSlope,
1518 )
1519 else:
1520 # Construct the unconstrained consumption function as a linear interpolation
1521 cFuncNowUnc = LinearInterp(
1522 m_for_interpolation,
1523 c_for_interpolation,
1524 cFuncLimitIntercept,
1525 cFuncLimitSlope,
1526 )
1528 # Combine the constrained and unconstrained functions into the true consumption function.
1529 # LowerEnvelope should only be used when BoroCnstArt is True
1530 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1532 # Make the marginal value function and the marginal marginal value function
1533 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1535 # Define this period's marginal marginal value function
1536 if CubicBool:
1537 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1538 else:
1539 vPPfuncNow = NullFunc() # Dummy object
1541 # Construct this period's value function if requested. This version is set
1542 # up for the non-independent distributions, need to write a faster version.
1543 if vFuncBool:
1544 # Compute expected value and marginal value on a grid of market resources
1545 mNrm_temp = mNrmMinNow + aXtraGrid
1546 cNrm_temp = cFuncNow(mNrm_temp)
1547 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # fix tiny errors
1548 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
1549 vP_temp = uFunc.der(cNrm_temp)
1551 # Construct the beginning-of-period value function
1552 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1553 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1554 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1555 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1556 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
1557 # MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1558 vNvrsFuncNow = CubicInterp(mNrm_temp, vNvrs_temp, vNvrsP_temp)
1559 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1560 else:
1561 vFuncNow = NullFunc() # Dummy object
1563 # Create and return this period's solution
1564 solution_now = ConsumerSolution(
1565 cFunc=cFuncNow,
1566 vFunc=vFuncNow,
1567 vPfunc=vPfuncNow,
1568 vPPfunc=vPPfuncNow,
1569 mNrmMin=mNrmMinNow,
1570 hNrm=hNrmNow,
1571 MPCmin=MPCminNow,
1572 MPCmax=MPCmaxEff,
1573 )
1574 solution_now.ShareFunc = ConstantFunction(RiskyShareFixed)
1575 return solution_now