Coverage for HARK / ConsumptionSaving / ConsRiskyAssetModel.py: 96%
492 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +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(PortfolioBool):
67 """
68 Trivial constructor function that chooses between two solvers.
69 """
70 if PortfolioBool:
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": 20, # Maximum end-of-period "assets above minimum" value
159 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
160 "aXtraCount": 48, # 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've (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 when True
195 # (Uses linear spline interpolation for cFunc when False)
196 "RiskyShareFixed": 1.0, # Fixed share of risky asset when PortfolioBool is False
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 "PortfolioBool": False, # Whether this instance can choose portfolio shares
200 "PortfolioBisect": False, # What does this do?
201 "pseudo_terminal": False,
202}
203IndShockRiskyAssetConsumerType_simulation_default = {
204 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
205 "AgentCount": 10000, # Number of agents of this type
206 "T_age": None, # Age after which simulated agents are automatically killed
207 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
208 # (The portion of PermGroFac attributable to aggregate productivity growth)
209 "NewbornTransShk": False, # Whether Newborns have transitory shock
210 # ADDITIONAL OPTIONAL PARAMETERS
211 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
212 # (Forces Newborns to follow solution path of the agent they replaced if True)
213 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
214 "sim_common_Rrisky": True, # Whether risky returns have a shared/common value across agents
215}
216IndShockRiskyAssetConsumerType_default = {}
217IndShockRiskyAssetConsumerType_default.update(
218 IndShockRiskyAssetConsumerType_IncShkDstn_default
219)
220IndShockRiskyAssetConsumerType_default.update(
221 IndShockRiskyAssetConsumerType_RiskyDstn_default
222)
223IndShockRiskyAssetConsumerType_default.update(
224 IndShockRiskyAssetConsumerType_aXtraGrid_default
225)
226IndShockRiskyAssetConsumerType_default.update(
227 IndShockRiskyAssetConsumerType_ShareGrid_default
228)
229IndShockRiskyAssetConsumerType_default.update(
230 IndShockRiskyAssetConsumerType_solving_default
231)
232IndShockRiskyAssetConsumerType_default.update(
233 IndShockRiskyAssetConsumerType_simulation_default
234)
235IndShockRiskyAssetConsumerType_default.update(
236 IndShockRiskyAssetConsumerType_kNrmInitDstn_default
237)
238IndShockRiskyAssetConsumerType_default.update(
239 IndShockRiskyAssetConsumerType_pLvlInitDstn_default
240)
241init_risky_asset = IndShockRiskyAssetConsumerType_default
244class IndShockRiskyAssetConsumerType(IndShockConsumerType):
245 r"""
246 A consumer type based on IndShockConsumerType, that has access to a risky asset for their savings. The
247 risky asset has lognormal returns that are possibly correlated with his
248 income shocks.
250 If PortfolioBool is False, then the risky asset share is always one.
251 Otherwise the agent can optimize 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.
313 PortfolioBool: Boolean
314 Determines whether agent will use portfolio optimization or they only have access to risky assets. If false, the risky share is always one.
316 Simulation Parameters
317 ---------------------
318 sim_common_Rrisky: Boolean
319 Whether risky returns have a shared/common value across agents. If True, Risky return's can't be time varying.
320 AgentCount: int
321 Number of agents of this kind that are created during simulations.
322 T_age: int
323 Age after which to automatically kill agents, None to ignore.
324 T_sim: int, required for simulation
325 Number of periods to simulate.
326 track_vars: list[strings]
327 List of variables that should be tracked when running the simulation.
328 For this agent, the options are 'Adjust', 'PermShk', 'Risky', 'TranShk', 'aLvl', 'aNrm', 'bNrm', 'cNrm', 'mNrm', 'pLvl', and 'who_dies'.
330 Adjust is the array of which agents can adjust
332 PermShk is the agent's permanent income shock
334 Risky is the agent's risky asset shock
336 TranShk is the agent's transitory income shock
338 aLvl is the nominal asset level
340 aNrm is the normalized assets
342 bNrm is the normalized resources without this period's labor income
344 cNrm is the normalized consumption
346 mNrm is the normalized market resources
348 pLvl is the permanent income level
350 who_dies is the array of which agents died
351 aNrmInitMean: float
352 Mean of Log initial Normalized Assets.
353 aNrmInitStd: float
354 Std of Log initial Normalized Assets.
355 pLvlInitMean: float
356 Mean of Log initial permanent income.
357 pLvlInitStd: float
358 Std of Log initial permanent income.
359 PermGroFacAgg: float
360 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
361 PerfMITShk: boolean
362 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
363 NewbornTransShk: boolean
364 Whether Newborns have transitory shock.
366 Attributes
367 ----------
368 solution: list[Consumer solution object]
369 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
370 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle. If PortfolioBool is True, the solution also contains ShareFunc.
372 If PortfolioBool is True, the solution also contains:
373 ShareFunc - The asset share function for this period, defined over normalized market resources :math:`S=ShareFunc(mNrm)`.
375 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
376 history: Dict[Array]
377 Created by running the :func:`.simulate()` method.
378 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
379 Visit :class:`HARK.core.AgentType.simulate` for more information.
380 """
382 IncShkDstn_default = IndShockRiskyAssetConsumerType_IncShkDstn_default
383 RiskyDstn_default = IndShockRiskyAssetConsumerType_RiskyDstn_default
384 aXtraGrid_default = IndShockRiskyAssetConsumerType_aXtraGrid_default
385 ShareGrid_default = IndShockRiskyAssetConsumerType_ShareGrid_default
386 solving_default = IndShockRiskyAssetConsumerType_solving_default
387 simulation_default = IndShockRiskyAssetConsumerType_simulation_default # So sphinx documents defaults
388 default_ = {
389 "params": IndShockRiskyAssetConsumerType_default,
390 "solver": NullFunc(),
391 "model": "ConsRiskyAsset.yaml",
392 "track_vars": ["aNrm", "cNrm", "Share", "mNrm", "pLvl"],
393 }
395 time_inv_ = IndShockConsumerType.time_inv_ + [
396 "PortfolioBisect",
397 "ShareGrid",
398 "PortfolioBool",
399 "IndepDstnBool",
400 "RiskyShareFixed",
401 ]
402 time_vary_ = IndShockConsumerType.time_vary_ + ["ShockDstn", "ShareLimit"]
403 shock_vars_ = IndShockConsumerType.shock_vars_ + ["Adjust", "Risky"]
404 distributions = [
405 "IncShkDstn",
406 "PermShkDstn",
407 "TranShkDstn",
408 "RiskyDstn",
409 "ShockDstn",
410 "kNrmInitDstn",
411 "pLvlInitDstn",
412 "RiskyDstn",
413 ]
415 def pre_solve(self):
416 self.construct("solution_terminal")
417 self.update_timing()
418 self.solution_terminal.ShareFunc = ConstantFunction(1.0)
420 def update_timing(self):
421 """
422 This method simply ensures that a few attributes that could be in either
423 time_inv or time_vary are appropriately labeled.
424 """
425 if type(self.AdjustDstn) is IndexDistribution:
426 self.add_to_time_vary("AdjustPrb")
427 self.del_from_time_inv("AdjustPrb")
428 else:
429 self.add_to_time_inv("AdjustPrb")
430 self.del_from_time_vary("AdjustPrb")
431 if hasattr(self.RiskyDstn, "__getitem__"):
432 self.add_to_time_vary("RiskyDstn")
433 else:
434 self.add_to_time_inv("RiskyDstn")
435 if type(self.ShareLimit) is list:
436 self.add_to_time_vary("ShareLimit")
437 self.del_from_time_inv("ShareLimit")
438 else:
439 self.add_to_time_inv("ShareLimit")
440 self.del_from_time_vary("ShareLimit")
442 def get_Rport(self):
443 """
444 Calculates realized return factor for each agent, using the attributes Rfree,
445 RiskyNow, and ShareNow.
447 Parameters
448 ----------
449 None
451 Returns
452 -------
453 Rport : np.array
454 Array of size AgentCount with each simulated agent's realized portfolio
455 return factor. Will be used by get_states() to calculate mNrmNow, where it
456 will be mislabeled as "Rfree".
457 """
459 RfreeNow = super().get_Rport()
460 RiskyNow = self.shocks["Risky"]
461 if self.PortfolioBool:
462 ShareNow = self.controls["Share"]
463 else:
464 ShareNow = np.ones_like(RiskyNow) # Only asset is risky asset
466 Rport = ShareNow * RiskyNow + (1.0 - ShareNow) * RfreeNow
467 self.Rport = Rport
468 return Rport
470 def get_Risky(self):
471 """
472 Draws a new risky return factor.
474 Parameters
475 ----------
476 None
478 Returns
479 -------
480 None
481 """
483 # How we draw the shocks depends on whether their distribution is time-varying
484 if "RiskyDstn" in self.time_vary:
485 if self.sim_common_Rrisky:
486 raise AttributeError(
487 "If sim_common_Rrisky is True, RiskyDstn cannot be time-varying!"
488 )
490 else:
491 # Make use of the IndexDistribution.draw() method
492 self.shocks["Risky"] = self.RiskyDstn.draw(
493 np.maximum(self.t_cycle - 1, 0)
494 if self.cycles == 1
495 else self.t_cycle
496 )
498 else:
499 # Draw either a common economy-wide return, or one for each agent
500 if self.sim_common_Rrisky:
501 self.shocks["Risky"] = self.RiskyDstn.draw(1)
502 else:
503 self.shocks["Risky"] = self.RiskyDstn.draw(self.AgentCount)
505 def get_Adjust(self):
506 """
507 Sets the attribute Adjust as a boolean array of size AgentCount, indicating
508 whether each agent is able to adjust their risky portfolio share this period.
509 Uses the attribute AdjustPrb to draw from a Bernoulli distribution.
511 Parameters
512 ----------
513 None
515 Returns
516 -------
517 None
518 """
519 if "AdjustPrb" in self.time_vary:
520 self.shocks["Adjust"] = self.AdjustDstn.draw(
521 np.maximum(self.t_cycle - 1, 0) if self.cycles == 1 else self.t_cycle
522 )
523 else:
524 self.shocks["Adjust"] = self.AdjustDstn.draw(self.AgentCount)
526 def initialize_sim(self):
527 """
528 Initialize the state of simulation attributes. Simply calls the same
529 method for IndShockConsumerType, then initializes the new states/shocks
530 Adjust and Share.
532 Parameters
533 ----------
534 None
536 Returns
537 -------
538 None
539 """
540 self.shocks["Adjust"] = np.zeros(self.AgentCount, dtype=bool)
541 # Initialize Share to default value; will be updated in get_controls()
542 self.controls["Share"] = np.ones(self.AgentCount)
543 IndShockConsumerType.initialize_sim(self)
545 def get_shocks(self):
546 """
547 Draw idiosyncratic income shocks, just as for IndShockConsumerType, then draw
548 a single common value for the risky asset return. Also draws whether each
549 agent is able to adjust their portfolio this period.
551 Parameters
552 ----------
553 None
555 Returns
556 -------
557 None
558 """
559 IndShockConsumerType.get_shocks(self)
560 self.get_Risky()
561 self.get_Adjust()
563 def get_controls(self):
564 """
565 Calculates consumption for each consumer of this type using the consumption functions;
566 also calculates risky asset share when PortfolioBool=True
568 Parameters
569 ----------
570 None
572 Returns
573 -------
574 None
575 """
576 cNrmNow = np.full(self.AgentCount, np.nan)
577 MPCnow = np.full(self.AgentCount, np.nan)
578 ShareNow = np.full(self.AgentCount, np.nan)
579 for t in np.unique(self.t_cycle):
580 idx = self.t_cycle == t
581 if np.any(idx):
582 mNrm = self.state_now["mNrm"][idx]
583 cNrmNow[idx], MPCnow[idx] = self.solution[t].cFunc.eval_with_derivative(
584 mNrm
585 )
586 if self.PortfolioBool:
587 ShareNow[idx] = self.solution[t].ShareFunc(mNrm)
588 else:
589 ShareNow[idx] = self.RiskyShareFixed
590 self.controls["cNrm"] = cNrmNow
591 self.controls["Share"] = ShareNow
592 self.MPCnow = MPCnow
594 def check_conditions(self, verbose=None):
595 raise NotImplementedError() # pragma: nocover
597 def calc_limiting_values(self):
598 raise NotImplementedError() # pragma: nocover
601# This is to preserve compatibility with old name
602RiskyAssetConsumerType = IndShockRiskyAssetConsumerType
605###############################################################################
608def solve_one_period_ConsPortChoice(
609 solution_next,
610 ShockDstn,
611 IncShkDstn,
612 RiskyDstn,
613 LivPrb,
614 DiscFac,
615 CRRA,
616 Rfree,
617 PermGroFac,
618 BoroCnstArt,
619 aXtraGrid,
620 ShareGrid,
621 ShareLimit,
622 vFuncBool,
623 IndepDstnBool,
624):
625 """
626 Solve one period of a consumption-saving problem with portfolio allocation
627 between a riskless and risky asset. This function handles only the most
628 fundamental portfolio choice problem: frictionless reallocation of the
629 portfolio each period as a continuous choice.
631 Parameters
632 ----------
633 solution_next : PortfolioSolution
634 Solution to next period's problem.
635 ShockDstn : Distribution
636 Joint distribution of permanent income shocks, transitory income shocks,
637 and risky returns. This is only used if the input IndepDstnBool is False,
638 indicating that income and return distributions can't be assumed to be
639 independent.
640 IncShkDstn : Distribution
641 Discrete distribution of permanent income shocks and transitory income
642 shocks. This is only used if the input IndepDstnBool is True, indicating
643 that income and return distributions are independent.
644 RiskyDstn : Distribution
645 Distribution of risky asset returns. This is only used if the input
646 IndepDstnBool is True, indicating that income and return distributions
647 are independent.
648 LivPrb : float
649 Survival probability; likelihood of being alive at the beginning of
650 the succeeding period.
651 DiscFac : float
652 Intertemporal discount factor for future utility.
653 CRRA : float
654 Coefficient of relative risk aversion.
655 Rfree : float
656 Risk free interest factor on end-of-period assets.
657 PermGroFac : float
658 Expected permanent income growth factor at the end of this period.
659 BoroCnstArt: float or None
660 Borrowing constraint for the minimum allowable assets to end the
661 period with. In this model, it is *required* to be zero.
662 aXtraGrid: np.array
663 Array of "extra" end-of-period asset values-- assets above the
664 absolute minimum acceptable level.
665 ShareGrid : np.array
666 Array of risky portfolio shares on which to define the interpolation
667 of the consumption function when Share is fixed. Also used when the
668 risky share choice is specified as discrete rather than continuous.
669 ShareLimit : float
670 Limiting lower bound of risky portfolio share as mNrm approaches infinity.
671 vFuncBool: boolean
672 An indicator for whether the value function should be computed and
673 included in the reported solution.
674 IndepDstnBool : bool
675 Indicator for whether the income and risky return distributions are in-
676 dependent of each other, which can speed up the expectations step.
678 Returns
679 -------
680 solution_now : PortfolioSolution
681 Solution to this period's problem.
683 :meta private:
684 """
685 # Make sure the individual is liquidity constrained. Allowing a consumer to
686 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
687 if BoroCnstArt != 0.0:
688 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
690 # Define the current period utility function and effective discount factor
691 uFunc = UtilityFuncCRRA(CRRA)
692 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
694 # Unpack next period's solution for easier access
695 vPfunc_next = solution_next.vPfunc
696 vFunc_next = solution_next.vFunc
698 # Set a flag for whether the natural borrowing constraint is zero, which
699 # depends on whether the smallest transitory income shock is zero
700 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
702 # Prepare to calculate end-of-period marginal values by creating an array
703 # of market resources that the agent could have next period, considering
704 # the grid of end-of-period assets and the distribution of shocks he might
705 # experience next period.
707 # Unpack the risky return shock distribution
708 Risky_next = RiskyDstn.atoms
709 RiskyMax = np.max(Risky_next)
710 RiskyMin = np.min(Risky_next)
712 # Perform an alternate calculation of the absolute patience factor when
713 # returns are risky. This uses the Merton-Samuelson limiting risky share,
714 # which is what's relevant as mNrm goes to infinity.
715 def calc_Radj(R):
716 Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree
717 return Rport ** (1.0 - CRRA)
719 R_adj = expected(calc_Radj, RiskyDstn)[0]
720 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
721 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
723 # Also perform an alternate calculation for human wealth under risky returns
724 def calc_hNrm(S):
725 Risky = S["Risky"]
726 PermShk = S["PermShk"]
727 TranShk = S["TranShk"]
728 G = PermGroFac * PermShk
729 Rport = ShareLimit * Risky + (1.0 - ShareLimit) * Rfree
730 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
731 return hNrm
733 # This correctly accounts for risky returns and risk aversion
734 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
736 # This basic equation works if there's no correlation among shocks
737 # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm)
739 # Define the terms for the limiting linear consumption function as m gets very big
740 cFuncLimitIntercept = MPCminNow * hNrmNow
741 cFuncLimitSlope = MPCminNow
743 # bNrm represents R*a, balances after asset return shocks but before income.
744 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
745 if BoroCnstNat_iszero:
746 aNrmGrid = aXtraGrid
747 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
748 else:
749 # Add an asset point at exactly zero
750 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
751 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
753 # Get grid and shock sizes, for easier indexing
754 aNrmCount = aNrmGrid.size
755 ShareCount = ShareGrid.size
757 # If the income shock distribution is independent from the risky return distribution,
758 # then taking end-of-period expectations can proceed in a two part process: First,
759 # construct an "intermediate" value function by integrating out next period's income
760 # shocks, *then* compute end-of-period expectations by integrating out return shocks.
761 # This method is lengthy to code, but can be significantly faster.
762 if IndepDstnBool:
763 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
764 bNrmNext = bNrmGrid
766 # Define functions that are used internally to evaluate future realizations
767 def calc_mNrm_next(S, b):
768 """
769 Calculate future realizations of market resources mNrm from the income
770 shock distribution S and normalized bank balances b.
771 """
772 return b / (S["PermShk"] * PermGroFac) + S["TranShk"]
774 def calc_dvdm_next(S, b):
775 """
776 Evaluate realizations of marginal value of market resources next period,
777 based on the income distribution S and values of bank balances bNrm
778 """
779 mNrm_next = calc_mNrm_next(S, b)
780 G = S["PermShk"] * PermGroFac
781 dvdm_next = G ** (-CRRA) * vPfunc_next(mNrm_next)
782 return dvdm_next
784 # Calculate end-of-period marginal value of assets and shares at each point
785 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
786 # values across income and risky return shocks.
788 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
789 dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext))
790 dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0))
792 dvdbNvrsFunc_intermed = LinearInterp(bNrmGrid, dvdbNvrs_intermed)
793 dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA)
795 # The intermediate marginal value of risky portfolio share is zero in this
796 # model because risky share is flexible: we can just change it next period,
797 # so there is no marginal value of Share once the return is realized.
798 dvdsFunc_intermed = ConstantFunction(0.0) # all zeros
800 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
801 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
803 # Define functions for calculating end-of-period marginal value
804 def calc_EndOfPrd_dvda(R, a, z):
805 """
806 Compute end-of-period marginal value of assets at values a, conditional
807 on risky asset return R and risky share z.
808 """
809 # Calculate future realizations of bank balances bNrm
810 Rxs = R - Rfree # Excess returns
811 Rport = Rfree + z * Rxs # Portfolio return
812 bNrm_next = Rport * a
814 # Calculate and return dvda
815 EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next)
816 return EndOfPrd_dvda
818 def calc_EndOfPrd_dvds(R, a, z):
819 """
820 Compute end-of-period marginal value of risky share at values a, conditional
821 on risky asset return S and risky share z.
822 """
823 # Calculate future realizations of bank balances bNrm
824 Rxs = R - Rfree # Excess returns
825 Rport = Rfree + z * Rxs # Portfolio return
826 bNrm_next = Rport * a
828 # Calculate and return dvds (second term is all zeros)
829 EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(
830 bNrm_next
831 )
832 return EndOfPrd_dvds
834 TempDstn = RiskyDstn # relabeling for below
836 # Evaluate realizations of value and marginal value after asset returns are realized
838 # Calculate end-of-period marginal value of risky portfolio share by taking expectations
839 EndOfPrd_dvds = DiscFacEff * expected(
840 calc_EndOfPrd_dvds, RiskyDstn, args=(aNrmNow, ShareNext)
841 )
843 # Make the end-of-period value function if the value function is requested
844 if vFuncBool:
846 def calc_v_intermed(S, b):
847 """
848 Calculate "intermediate" value from next period's bank balances, the
849 income shocks S, and the risky asset share.
850 """
851 mNrm_next = calc_mNrm_next(S, b)
852 v_next = vFunc_next(mNrm_next)
853 v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
854 return v_intermed
856 # Calculate intermediate value by taking expectations over income shocks
857 v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext))
859 # Construct the "intermediate value function" for this period
860 vNvrs_intermed = uFunc.inv(v_intermed)
861 vNvrsFunc_intermed = LinearInterp(bNrmGrid, vNvrs_intermed)
862 vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA)
864 def calc_EndOfPrd_v(S, a, z):
865 # Calculate future realizations of bank balances bNrm
866 Rxs = S - Rfree
867 Rport = Rfree + z * Rxs
868 bNrm_next = Rport * a
870 EndOfPrd_v = vFunc_intermed(bNrm_next)
871 return EndOfPrd_v
873 # Calculate end-of-period value by taking expectations
874 EndOfPrd_v = DiscFacEff * expected(
875 calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext)
876 )
877 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
879 # Now make an end-of-period value function over aNrm and Share
880 EndOfPrd_vNvrsFunc = BilinearInterp(EndOfPrd_vNvrs, aNrmGrid, ShareGrid)
881 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
882 # This will be used later to make the value function for this period
884 # If the income shock distribution and risky return distribution are *NOT*
885 # independent, then computation of end-of-period expectations are simpler in
886 # code, but might take longer to execute
887 else:
888 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
889 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
891 # Define functions that are used internally to evaluate future realizations
892 def calc_mNrm_next(S, a, z):
893 """
894 Calculate future realizations of market resources mNrm from the shock
895 distribution S, normalized end-of-period assets a, and risky share z.
896 """
897 # Calculate future realizations of bank balances bNrm
898 Rxs = S["Risky"] - Rfree
899 Rport = Rfree + z * Rxs
900 bNrm_next = Rport * a
901 mNrm_next = bNrm_next / (S["PermShk"] * PermGroFac) + S["TranShk"]
902 return mNrm_next
904 def calc_EndOfPrd_dvdx(S, a, z):
905 """
906 Evaluate end-of-period marginal value of assets and risky share based
907 on the shock distribution S, normalized end-of-period assets a, and
908 risky share z.
909 """
910 mNrm_next = calc_mNrm_next(S, a, z)
911 Rxs = S["Risky"] - Rfree
912 Rport = Rfree + z * Rxs
913 dvdm_next = vPfunc_next(mNrm_next)
914 # No marginal value of Share if it's a free choice!
915 dvds_next = np.zeros_like(mNrm_next)
917 EndOfPrd_dvda = Rport * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
918 EndOfPrd_dvds = (
919 Rxs * a * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
920 + (S["PermShk"] * PermGroFac) ** (1 - CRRA) * dvds_next
921 )
923 return EndOfPrd_dvda, EndOfPrd_dvds
925 def calc_EndOfPrd_v(S, a, z):
926 """
927 Evaluate end-of-period value, based on the shock distribution S, values
928 of bank balances bNrm, and values of the risky share z.
929 """
930 mNrm_next = calc_mNrm_next(S, a, z)
931 v_next = vFunc_next(mNrm_next)
932 EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
933 return EndOfPrd_v
935 calc_EndOfPrd_dvda = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[0]
936 calc_EndOfPrd_dvds = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[1]
937 TempDstn = ShockDstn
939 # Evaluate realizations of value and marginal value after asset returns are realized
941 # Calculate end-of-period marginal value of assets and risky share by taking expectations
942 EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected(
943 calc_EndOfPrd_dvdx, ShockDstn, args=(aNrmNow, ShareNext)
944 )
945 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
947 # Construct the end-of-period value function if requested
948 if vFuncBool:
949 # Calculate end-of-period value, its derivative, and their pseudo-inverse
950 EndOfPrd_v = DiscFacEff * expected(
951 calc_EndOfPrd_v, ShockDstn, args=(aNrmNow, ShareNext)
952 )
953 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
955 # value transformed through inverse utility
956 EndOfPrd_vNvrsP = EndOfPrd_dvda * uFunc.derinv(EndOfPrd_v, order=(0, 1))
958 # Construct the end-of-period value function
959 EndOfPrd_vNvrsFunc_by_Share = []
960 for j in range(ShareCount):
961 EndOfPrd_vNvrsFunc_by_Share.append(
962 CubicInterp(
963 aNrmNow[:, j], EndOfPrd_vNvrs[:, j], EndOfPrd_vNvrsP[:, j]
964 )
965 )
966 EndOfPrd_vNvrsFunc = LinearInterpOnInterp1D(
967 EndOfPrd_vNvrsFunc_by_Share, ShareGrid
968 )
969 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
971 # Now find the optimal (continuous) risky share on [0,1] by solving the first
972 # order condition EndOfPrd_dvds == 0.
973 FOC_s = EndOfPrd_dvds # Relabel for convenient typing
975 # If agent wants to put more than 100% into risky asset, he is constrained.
976 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
977 constrained_top = FOC_s[:, -1] > 0.0
978 constrained_bot = FOC_s[:, 0] < 0.0
979 constrained = np.logical_or(constrained_top, constrained_bot)
980 a_idx = np.arange(aNrmCount)
982 # For each value of aNrm, find the value of Share such that FOC_s == 0
983 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0)
984 share_idx = np.argmax(crossing, axis=1)
986 for k in range(3):
987 # This represents the index of the segment of the share grid where dvds flips
988 # from positive to negative, indicating that there's a zero *on* the segment.
989 # The exception is for aNrm values that are flagged as constrained, for which
990 # there will be no crossing point and we can just use the boundary value.
992 # Now that we have a *range* for the location of the optimal share, we can
993 # do a refined search for the optimal share at each aNrm value where there
994 # is an interior solution (not constrained). We now make a refined ShareGrid
995 # that has *different* values on it for each aNrm value.
996 bot_s = ShareNext[a_idx, share_idx]
997 top_s = ShareNext[a_idx, share_idx + 1]
998 for j in range(aNrmCount):
999 if constrained[j]:
1000 continue
1001 ShareNext[j, :] = np.linspace(bot_s[j], top_s[j], ShareCount)
1003 # Now evaluate end-of-period marginal value of risky share on the refined grid
1004 EndOfPrd_dvds = DiscFacEff * expected(
1005 calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)
1006 )
1007 these = np.logical_not(constrained)
1008 FOC_s[these, :] = EndOfPrd_dvds[these, :] # Relabel for convenient typing
1010 # Look for "crossing points" again
1011 crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0)
1012 share_idx[these] = np.argmax(crossing, axis=1)
1014 # Calculate end-of-period marginal value of assets on the refined grid
1015 EndOfPrd_dvda = DiscFacEff * expected(
1016 calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)
1017 )
1018 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
1020 # Calculate the fractional distance between those share gridpoints where the
1021 # zero should be found, assuming a linear function; call it alpha
1022 bot_s = ShareNext[a_idx, share_idx]
1023 top_s = ShareNext[a_idx, share_idx + 1]
1024 bot_f = FOC_s[a_idx, share_idx]
1025 top_f = FOC_s[a_idx, share_idx + 1]
1026 bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx]
1027 top_c = EndOfPrd_dvdaNvrs[a_idx, share_idx + 1]
1028 alpha = 1.0 - top_f / (top_f - bot_f)
1030 # Calculate the continuous optimal risky share and optimal consumption
1031 Share_now = (1.0 - alpha) * bot_s + alpha * top_s
1032 cNrm_now = (1.0 - alpha) * bot_c + alpha * top_c
1034 # If agent wants to put more than 100% into risky asset, he is constrained.
1035 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
1036 constrained_top = FOC_s[:, -1] > 0.0
1037 constrained_bot = FOC_s[:, 0] < 0.0
1039 # Apply the constraints to both risky share and consumption (but lower
1040 # constraint should never be relevant)
1041 Share_now[constrained_top] = 1.0
1042 Share_now[constrained_bot] = 0.0
1043 cNrm_now[constrained_top] = EndOfPrd_dvdaNvrs[constrained_top, -1]
1044 cNrm_now[constrained_bot] = EndOfPrd_dvdaNvrs[constrained_bot, 0]
1046 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
1047 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
1048 # can't depend on the risky share if it doesn't meaningfully exist. Apply
1049 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
1050 if not BoroCnstNat_iszero:
1051 Share_now[0] = 1.0
1052 cNrm_now[0] = EndOfPrd_dvdaNvrs[0, -1]
1054 # Construct functions characterizing the solution for this period
1056 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
1057 # then construct the consumption function when the agent can adjust his share
1058 mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0)
1059 cNrm_now = np.insert(cNrm_now, 0, 0.0)
1060 cFunc_now = LinearInterp(mNrm_now, cNrm_now, cFuncLimitIntercept, cFuncLimitSlope)
1062 # Construct the marginal value (of mNrm) function
1063 vPfunc_now = MargValueFuncCRRA(cFunc_now, CRRA)
1065 # If the share choice is continuous, just make an ordinary interpolating function
1066 if BoroCnstNat_iszero:
1067 Share_lower_bound = ShareLimit
1068 else:
1069 Share_lower_bound = 1.0
1070 Share_now = np.insert(Share_now, 0, Share_lower_bound)
1071 ShareFunc_now = LinearInterp(mNrm_now, Share_now, ShareLimit, 0.0)
1073 # Add the value function if requested
1074 if vFuncBool:
1075 # Create the value functions for this period, defined over market resources
1076 # mNrm when agent can adjust his portfolio, and over market resources and
1077 # fixed share when agent can not adjust his portfolio.
1079 # Construct the value function
1080 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
1081 cNrm_temp = cFunc_now(mNrm_temp)
1082 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # Fix tiny violations
1083 Share_temp = ShareFunc_now(mNrm_temp)
1084 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
1085 vNvrs_temp = uFunc.inv(v_temp)
1086 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
1087 vNvrsFunc = CubicInterp(
1088 np.insert(mNrm_temp, 0, 0.0), # x_list
1089 np.insert(vNvrs_temp, 0, 0.0), # f_list
1090 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
1091 )
1092 # Re-curve the pseudo-inverse value function
1093 vFunc_now = ValueFuncCRRA(vNvrsFunc, CRRA)
1095 else: # If vFuncBool is False, fill in dummy values
1096 vFunc_now = NullFunc()
1098 # Package and return the solution
1099 solution_now = ConsumerSolution(
1100 cFunc=cFunc_now,
1101 vPfunc=vPfunc_now,
1102 vFunc=vFunc_now,
1103 hNrm=hNrmNow,
1104 MPCmin=MPCminNow,
1105 )
1106 solution_now.ShareFunc = ShareFunc_now
1107 return solution_now
1110###############################################################################
1113def solve_one_period_ConsIndShockRiskyAsset(
1114 solution_next,
1115 IncShkDstn,
1116 RiskyDstn,
1117 ShockDstn,
1118 LivPrb,
1119 DiscFac,
1120 Rfree,
1121 CRRA,
1122 PermGroFac,
1123 BoroCnstArt,
1124 aXtraGrid,
1125 RiskyShareFixed,
1126 vFuncBool,
1127 CubicBool,
1128 IndepDstnBool,
1129):
1130 """
1131 Solves one period of a consumption-saving model with idiosyncratic shocks to
1132 permanent and transitory income, with one risky asset and CRRA utility.
1134 Parameters
1135 ----------
1136 solution_next : ConsumerSolution
1137 The solution to next period's one period problem.
1138 IncShkDstn : Distribution
1139 Discrete distribution of permanent income shocks and transitory income
1140 shocks. This is only used if the input IndepDstnBool is True, indicating
1141 that income and return distributions are independent.
1142 RiskyDstn : Distribution
1143 Distribution of risky asset returns. This is only used if the input
1144 IndepDstnBool is True, indicating that income and return distributions
1145 are independent.
1146 ShockDstn : Distribution
1147 Joint distribution of permanent income shocks, transitory income shocks,
1148 and risky returns. This is only used if the input IndepDstnBool is False,
1149 indicating that income and return distributions can't be assumed to be
1150 independent.
1151 LivPrb : float
1152 Survival probability; likelihood of being alive at the beginning of
1153 the succeeding period.
1154 DiscFac : float
1155 Intertemporal discount factor for future utility.
1156 Rfree : float
1157 Risk free interest factor on end-of-period assets.
1158 CRRA : float
1159 Coefficient of relative risk aversion.
1160 PermGroFac : float
1161 Expected permanent income growth factor at the end of this period.
1162 BoroCnstArt: float or None
1163 Borrowing constraint for the minimum allowable assets to end the
1164 period with. If it is less than the natural borrowing constraint,
1165 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
1166 rowing constraint.
1167 aXtraGrid: np.array
1168 Array of "extra" end-of-period asset values-- assets above the
1169 absolute minimum acceptable level.
1170 RiskyShareFixed : float
1171 Fixed fraction of end-of-period assets that are allocated to the risky asset.
1172 vFuncBool: boolean
1173 An indicator for whether the value function should be computed and
1174 included in the reported solution.
1175 CubicBool: boolean
1176 An indicator for whether the solver should use cubic or linear interpolation.
1177 IndepDstnBool : bool
1178 Indicator for whether the income and risky return distributions are in-
1179 dependent of each other, which can speed up the expectations step.
1181 Returns
1182 -------
1183 solution_now : ConsumerSolution
1184 Solution to this period's consumption-saving problem with income risk.
1186 :meta private:
1187 """
1188 # Do a quick validity check; don't want to allow borrowing with risky returns
1189 if BoroCnstArt != 0.0:
1190 raise ValueError("RiskyAssetConsumerType must have BoroCnstArt=0.0!")
1192 # Define the current period utility function and effective discount factor
1193 uFunc = UtilityFuncCRRA(CRRA)
1194 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
1196 # Unpack next period's income shock distribution
1197 ShkPrbsNext = ShockDstn.pmv
1198 PermShkValsNext = ShockDstn.atoms[0]
1199 TranShkValsNext = ShockDstn.atoms[1]
1200 RiskyValsNext = ShockDstn.atoms[2]
1201 PermShkMinNext = np.min(PermShkValsNext)
1202 TranShkMinNext = np.min(TranShkValsNext)
1203 RiskyMinNext = np.min(RiskyValsNext)
1204 RiskyMaxNext = np.max(RiskyValsNext)
1206 # Unpack next period's (marginal) value function
1207 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
1208 vPfuncNext = solution_next.vPfunc
1209 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
1211 # Perform an alternate calculation of the absolute patience factor when returns are risky
1212 def calc_Radj(R):
1213 Rport = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree
1214 return Rport ** (1.0 - CRRA)
1216 R_adj = expected(calc_Radj, RiskyDstn)
1217 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
1218 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
1219 MPCminNow = MPCminNow[0]
1221 # Also perform an alternate calculation for human wealth under risky returns
1222 def calc_hNrm(S):
1223 Risky = S["Risky"]
1224 PermShk = S["PermShk"]
1225 TranShk = S["TranShk"]
1226 G = PermGroFac * PermShk
1227 Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree
1228 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
1229 return hNrm
1231 # This correctly accounts for risky returns and risk aversion
1232 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
1233 hNrmNow = hNrmNow[0]
1235 # The above attempts to pin down the limiting consumption function for this
1236 # model, however it is not clear why it creates bugs, so for now we allow
1237 # for a linear extrapolation beyond the last asset point
1238 cFuncLimitIntercept = MPCminNow * hNrmNow
1239 cFuncLimitSlope = MPCminNow
1241 # Calculate the minimum allowable value of market resources in this period
1242 BoroCnstNat_cand = (
1243 (solution_next.mNrmMin - TranShkValsNext)
1244 * (PermGroFac * PermShkValsNext)
1245 / RiskyValsNext
1246 )
1247 BoroCnstNat = np.max(BoroCnstNat_cand) # Must be at least this
1249 # Set a flag for whether the natural borrowing constraint is zero, which
1250 # depends on whether the smallest transitory income shock is zero
1251 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
1253 # Set the minimum allowable (normalized) market resources based on the natural
1254 # and artificial borrowing constraints
1255 if BoroCnstArt is None:
1256 mNrmMinNow = BoroCnstNat
1257 else:
1258 mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt])
1260 # The MPCmax code is a bit unusual here, and possibly "harmlessly wrong".
1261 # The "worst event" should depend on the risky return factor as well as
1262 # income shocks. However, the natural borrowing constraint is only ever
1263 # relevant in this model when it's zero, so the MPC at mNrm is only relevant
1264 # in the case where risky returns don't matter at all (because a=0).
1266 # Calculate the probability that we get the worst possible income draw
1267 IncNext = PermShkValsNext * TranShkValsNext
1268 WorstIncNext = PermShkMinNext * TranShkMinNext
1269 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
1270 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
1272 # Update the upper bounding MPC as market resources approach the lower bound
1273 temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
1274 MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
1276 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
1277 # or artificial borrowing constraint actually binds
1278 if BoroCnstNat < mNrmMinNow:
1279 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
1280 else:
1281 MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above
1283 # Define the borrowing-constrained consumption function
1284 cFuncNowCnst = LinearInterp(
1285 np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0])
1286 )
1288 # Big methodological split here: whether the income and return distributions are independent.
1289 # Calculation of end-of-period marginal (marginal) value uses different approaches
1290 if IndepDstnBool:
1291 # bNrm represents R*a, balances after asset return shocks but before income.
1292 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
1293 if BoroCnstNat_iszero:
1294 bNrmNow = np.insert(
1295 RiskyMaxNext * aXtraGrid, 0, RiskyMinNext * aXtraGrid[0]
1296 )
1297 aNrmNow = aXtraGrid.copy()
1298 else:
1299 # Add a bank balances point at exactly zero
1300 bNrmNow = RiskyMaxNext * np.insert(aXtraGrid, 0, 0.0)
1301 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1303 # Define local functions for taking future expectations when the interest
1304 # factor *is* independent from the income shock distribution. These go
1305 # from "bank balances" bNrm = R * aNrm to t+1 realizations.
1306 def calc_mNrmNext(S, b):
1307 return b / (PermGroFac * S["PermShk"]) + S["TranShk"]
1309 def calc_vNext(S, b):
1310 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, b))
1312 def calc_vPnext(S, b):
1313 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, b))
1315 def calc_vPPnext(S, b):
1316 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, b))
1318 # Calculate marginal value of bank balances at each gridpoint
1319 vPfacEff = PermGroFac ** (-CRRA)
1320 Intermed_vP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(bNrmNow))
1321 Intermed_vPnvrs = uFunc.derinv(Intermed_vP, order=(1, 0))
1323 if BoroCnstNat_iszero:
1324 Intermed_vPnvrs = np.insert(Intermed_vPnvrs, 0, 0.0)
1325 bNrm_temp = np.insert(bNrmNow, 0, 0.0)
1326 else:
1327 bNrm_temp = bNrmNow.copy()
1329 # If using cubic spline interpolation, also calculate "intermediate"
1330 # marginal marginal value of bank balances
1331 if CubicBool:
1332 vPPfacEff = PermGroFac ** (-CRRA - 1.0)
1333 Intermed_vPP = vPPfacEff * expected(
1334 calc_vPPnext, IncShkDstn, args=(bNrmNow)
1335 )
1336 Intermed_vPnvrsP = Intermed_vPP * uFunc.derinv(Intermed_vP, order=(1, 1))
1337 if BoroCnstNat_iszero:
1338 Intermed_vPnvrsP = np.insert(Intermed_vPnvrsP, 0, Intermed_vPnvrsP[0])
1340 # Make a cubic spline intermediate pseudo-inverse marginal value function
1341 Intermed_vPnvrsFunc = CubicInterp(
1342 bNrm_temp,
1343 Intermed_vPnvrs,
1344 Intermed_vPnvrsP,
1345 lower_extrap=True,
1346 )
1347 Intermed_vPPfunc = MargMargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1348 else:
1349 # Make a linear interpolation intermediate pseudo-inverse marginal value function
1350 Intermed_vPnvrsFunc = LinearInterp(
1351 bNrm_temp, Intermed_vPnvrs, lower_extrap=True
1352 )
1354 # "Recurve" the intermediate pseudo-inverse marginal value function
1355 Intermed_vPfunc = MargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1357 # If the value function is requested, calculate "intermediate" value
1358 if vFuncBool:
1359 vFacEff = PermGroFac ** (1.0 - CRRA)
1360 Intermed_v = vFacEff * expected(calc_vNext, IncShkDstn, args=(bNrmNow))
1361 Intermed_vNvrs = uFunc.inv(Intermed_v)
1362 # value transformed through inverse utility
1363 Intermed_vNvrsP = Intermed_vP * uFunc.derinv(Intermed_v, order=(0, 1))
1364 if BoroCnstNat_iszero:
1365 Intermed_vNvrs = np.insert(Intermed_vNvrs, 0, 0.0)
1366 Intermed_vNvrsP = np.insert(Intermed_vNvrsP, 0, Intermed_vNvrsP[0])
1367 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1369 # Make a cubic spline intermediate pseudo-inverse value function
1370 Intermed_vNvrsFunc = CubicInterp(bNrm_temp, Intermed_vNvrs, Intermed_vNvrsP)
1372 # "Recurve" the intermediate pseudo-inverse value function
1373 Intermed_vFunc = ValueFuncCRRA(Intermed_vNvrsFunc, CRRA)
1375 # We have "intermediate" (marginal) value functions defined over bNrm,
1376 # so now we want to take expectations over Risky realizations at each aNrm.
1378 # Begin by re-defining transition functions for taking expectations, which are all very simple!
1379 Z = RiskyShareFixed # for shorter notation
1381 def calc_bNrmNext(R, a):
1382 Rport = Z * R + (1 - Z) * Rfree
1383 return Rport * a
1385 def calc_vNext(R, a):
1386 return Intermed_vFunc(calc_bNrmNext(R, a))
1388 def calc_vPnext(R, a):
1389 Rport = Z * R + (1 - Z) * Rfree
1390 return Rport * Intermed_vPfunc(calc_bNrmNext(R, a))
1392 def calc_vPPnext(R, a):
1393 Rport = Z * R + (1 - Z) * Rfree
1394 return Rport * Rport * Intermed_vPPfunc(calc_bNrmNext(R, a))
1396 # Calculate end-of-period marginal value of assets at each gridpoint
1397 EndOfPrdvP = DiscFacEff * expected(calc_vPnext, RiskyDstn, args=(aNrmNow))
1399 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1400 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1401 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1403 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1404 if CubicBool:
1405 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1406 EndOfPrdvPP = DiscFacEff * expected(calc_vPPnext, RiskyDstn, args=(aNrmNow))
1407 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1408 MPC = dcda / (dcda + 1.0)
1409 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1411 # Limiting consumption is zero as m approaches mNrmMin
1412 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1413 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1415 # Construct the end-of-period value function if requested
1416 if vFuncBool:
1417 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1418 EndOfPrdv = DiscFacEff * expected(calc_vNext, RiskyDstn, args=(aNrmNow))
1419 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1420 # value transformed through inverse utility
1421 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1423 # Construct the end-of-period value function
1424 if BoroCnstNat_iszero:
1425 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1426 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1427 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1428 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1429 else:
1430 aNrm_temp = aNrmNow.copy()
1432 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1433 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1435 # NON-INDEPENDENT METHOD BEGINS HERE
1436 else:
1437 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
1438 # aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
1439 if BoroCnstNat_iszero:
1440 aNrmNow = aXtraGrid
1441 else:
1442 # Add an asset point at exactly zero
1443 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1445 # Define local functions for taking future expectations when the interest
1446 # factor is *not* independent from the income shock distribution
1447 Z = RiskyShareFixed # for shorter notation
1449 def calc_mNrmNext(S, a):
1450 Risky = S["Risky"]
1451 Rport = Z * Risky + (1 - Z) * Rfree
1452 return Rport / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
1454 def calc_vNext(S, a):
1455 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, a))
1457 def calc_vPnext(S, a):
1458 Risky = S["Risky"]
1459 Rport = Z * Risky + (1 - Z) * Rfree
1460 return Rport * S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a))
1462 def calc_vPPnext(S, a):
1463 Risky = S["Risky"]
1464 Rport = Z * Risky + (1 - Z) * Rfree
1465 return (
1466 (Rport**2)
1467 * S["PermShk"] ** (-CRRA - 1.0)
1468 * vPPfuncNext(calc_mNrmNext(S, a))
1469 )
1471 # Calculate end-of-period marginal value of assets at each gridpoint
1472 vPfacEff = DiscFacEff * PermGroFac ** (-CRRA)
1473 EndOfPrdvP = vPfacEff * expected(calc_vPnext, ShockDstn, args=(aNrmNow))
1475 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1476 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1477 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1479 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1480 if CubicBool:
1481 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1482 vPPfacEff = DiscFacEff * PermGroFac ** (-CRRA - 1.0)
1483 EndOfPrdvPP = vPPfacEff * expected(calc_vPPnext, ShockDstn, args=(aNrmNow))
1484 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1485 MPC = dcda / (dcda + 1.0)
1486 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1488 # Limiting consumption is zero as m approaches mNrmMin
1489 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1490 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1492 # Construct the end-of-period value function if requested
1493 if vFuncBool:
1494 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1495 vFacEff = DiscFacEff * PermGroFac ** (1.0 - CRRA)
1496 EndOfPrdv = vFacEff * expected(calc_vNext, ShockDstn, args=(aNrmNow))
1497 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1498 # value transformed through inverse utility
1499 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1501 # Construct the end-of-period value function
1502 if BoroCnstNat_iszero:
1503 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1504 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1505 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1506 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1507 else:
1508 aNrm_temp = aNrmNow.copy()
1510 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1511 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1513 # Construct the consumption function; this uses the same method whether the
1514 # income distribution is independent from the return distribution
1515 if CubicBool:
1516 # Construct the unconstrained consumption function as a cubic interpolation
1517 cFuncNowUnc = CubicInterp(
1518 m_for_interpolation,
1519 c_for_interpolation,
1520 MPC_for_interpolation,
1521 cFuncLimitIntercept,
1522 cFuncLimitSlope,
1523 )
1524 else:
1525 # Construct the unconstrained consumption function as a linear interpolation
1526 cFuncNowUnc = LinearInterp(
1527 m_for_interpolation,
1528 c_for_interpolation,
1529 cFuncLimitIntercept,
1530 cFuncLimitSlope,
1531 )
1533 # Combine the constrained and unconstrained functions into the true consumption function.
1534 # LowerEnvelope should only be used when BoroCnstArt is True
1535 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1537 # Make the marginal value function and the marginal marginal value function
1538 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1540 # Define this period's marginal marginal value function
1541 if CubicBool:
1542 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1543 else:
1544 vPPfuncNow = NullFunc() # Dummy object
1546 # Construct this period's value function if requested. This version is set
1547 # up for the non-independent distributions, need to write a faster version.
1548 if vFuncBool:
1549 # Compute expected value and marginal value on a grid of market resources
1550 mNrm_temp = mNrmMinNow + aXtraGrid
1551 cNrm_temp = cFuncNow(mNrm_temp)
1552 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # fix tiny errors
1553 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
1554 vP_temp = uFunc.der(cNrm_temp)
1556 # Construct the beginning-of-period value function
1557 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1558 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1559 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1560 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1561 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
1562 # MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1563 vNvrsFuncNow = CubicInterp(mNrm_temp, vNvrs_temp, vNvrsP_temp)
1564 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1565 else:
1566 vFuncNow = NullFunc() # Dummy object
1568 # Create and return this period's solution
1569 solution_now = ConsumerSolution(
1570 cFunc=cFuncNow,
1571 vFunc=vFuncNow,
1572 vPfunc=vPfuncNow,
1573 vPPfunc=vPPfuncNow,
1574 mNrmMin=mNrmMinNow,
1575 hNrm=hNrmNow,
1576 MPCmin=MPCminNow,
1577 MPCmax=MPCmaxEff,
1578 )
1579 solution_now.ShareFunc = ConstantFunction(RiskyShareFixed)
1580 return solution_now