Coverage for HARK / ConsumptionSaving / ConsRiskyAssetModel.py: 95%
494 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-07 05:16 +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 }
394 time_inv_ = IndShockConsumerType.time_inv_ + [
395 "PortfolioBisect",
396 "ShareGrid",
397 "PortfolioBool",
398 "IndepDstnBool",
399 "RiskyShareFixed",
400 ]
401 time_vary_ = IndShockConsumerType.time_vary_ + ["ShockDstn", "ShareLimit"]
402 shock_vars_ = IndShockConsumerType.shock_vars_ + ["Adjust", "Risky"]
403 distributions = [
404 "IncShkDstn",
405 "PermShkDstn",
406 "TranShkDstn",
407 "RiskyDstn",
408 "ShockDstn",
409 "kNrmInitDstn",
410 "pLvlInitDstn",
411 "RiskyDstn",
412 ]
414 def pre_solve(self):
415 self.construct("solution_terminal")
416 self.update_timing()
417 self.solution_terminal.ShareFunc = ConstantFunction(1.0)
419 def update_timing(self):
420 """
421 This method simply ensures that a few attributes that could be in either
422 time_inv or time_vary are appropriately labeled.
423 """
424 if type(self.AdjustDstn) is IndexDistribution:
425 self.add_to_time_vary("AdjustPrb")
426 self.del_from_time_inv("AdjustPrb")
427 else:
428 self.add_to_time_inv("AdjustPrb")
429 self.del_from_time_vary("AdjustPrb")
430 if hasattr(self.RiskyDstn, "__getitem__"):
431 self.add_to_time_vary("RiskyDstn")
432 else:
433 self.add_to_time_inv("RiskyDstn")
434 if type(self.ShareLimit) is list:
435 self.add_to_time_vary("ShareLimit")
436 self.del_from_time_inv("ShareLimit")
437 else:
438 self.add_to_time_inv("ShareLimit")
439 self.del_from_time_vary("ShareLimit")
441 def get_Rfree(self):
442 """
443 Calculates realized return factor for each agent, using the attributes Rfree,
444 RiskyNow, and ShareNow. This method is a bit of a misnomer, as the return
445 factor is not riskless, but would more accurately be labeled as Rport. However,
446 this method makes the portfolio model compatible with its parent class.
448 Parameters
449 ----------
450 None
452 Returns
453 -------
454 Rport : np.array
455 Array of size AgentCount with each simulated agent's realized portfolio
456 return factor. Will be used by get_states() to calculate mNrmNow, where it
457 will be mislabeled as "Rfree".
458 """
460 RfreeNow = super().get_Rfree()
461 RiskyNow = self.shocks["Risky"]
462 if self.PortfolioBool:
463 ShareNow = self.controls["Share"]
464 else:
465 ShareNow = np.ones_like(RiskyNow) # Only asset is risky asset
467 Rport = ShareNow * RiskyNow + (1.0 - ShareNow) * RfreeNow
468 self.Rport = Rport
469 return Rport
471 def get_Risky(self):
472 """
473 Draws a new risky return factor.
475 Parameters
476 ----------
477 None
479 Returns
480 -------
481 None
482 """
484 # How we draw the shocks depends on whether their distribution is time-varying
485 if "RiskyDstn" in self.time_vary:
486 if self.sim_common_Rrisky:
487 raise AttributeError(
488 "If sim_common_Rrisky is True, RiskyDstn cannot be time-varying!"
489 )
491 else:
492 # Make use of the IndexDistribution.draw() method
493 self.shocks["Risky"] = self.RiskyDstn.draw(
494 np.maximum(self.t_cycle - 1, 0)
495 if self.cycles == 1
496 else self.t_cycle
497 )
499 else:
500 # Draw either a common economy-wide return, or one for each agent
501 if self.sim_common_Rrisky:
502 self.shocks["Risky"] = self.RiskyDstn.draw(1)
503 else:
504 self.shocks["Risky"] = self.RiskyDstn.draw(self.AgentCount)
506 def get_Adjust(self):
507 """
508 Sets the attribute Adjust as a boolean array of size AgentCount, indicating
509 whether each agent is able to adjust their risky portfolio share this period.
510 Uses the attribute AdjustPrb to draw from a Bernoulli distribution.
512 Parameters
513 ----------
514 None
516 Returns
517 -------
518 None
519 """
520 if "AdjustPrb" in self.time_vary:
521 self.shocks["Adjust"] = self.AdjustDstn.draw(
522 np.maximum(self.t_cycle - 1, 0) if self.cycles == 1 else self.t_cycle
523 )
524 else:
525 self.shocks["Adjust"] = self.AdjustDstn.draw(self.AgentCount)
527 def initialize_sim(self):
528 """
529 Initialize the state of simulation attributes. Simply calls the same
530 method for IndShockConsumerType, then initializes the new states/shocks
531 Adjust and Share.
533 Parameters
534 ----------
535 None
537 Returns
538 -------
539 None
540 """
541 self.shocks["Adjust"] = np.zeros(self.AgentCount, dtype=bool)
542 # Initialize Share to default value; will be updated in get_controls()
543 self.controls["Share"] = np.ones(self.AgentCount)
544 IndShockConsumerType.initialize_sim(self)
546 def get_shocks(self):
547 """
548 Draw idiosyncratic income shocks, just as for IndShockConsumerType, then draw
549 a single common value for the risky asset return. Also draws whether each
550 agent is able to adjust their portfolio this period.
552 Parameters
553 ----------
554 None
556 Returns
557 -------
558 None
559 """
560 IndShockConsumerType.get_shocks(self)
561 self.get_Risky()
562 self.get_Adjust()
564 def get_controls(self):
565 """
566 Calculates consumption for each consumer of this type using the consumption functions;
567 also calculates risky asset share when PortfolioBool=True
569 Parameters
570 ----------
571 None
573 Returns
574 -------
575 None
576 """
577 cNrmNow = np.full(self.AgentCount, np.nan)
578 MPCnow = np.full(self.AgentCount, np.nan)
579 ShareNow = np.full(self.AgentCount, np.nan)
580 for t in np.unique(self.t_cycle):
581 idx = self.t_cycle == t
582 if np.any(idx):
583 mNrm = self.state_now["mNrm"][idx]
584 cNrmNow[idx], MPCnow[idx] = self.solution[t].cFunc.eval_with_derivative(
585 mNrm
586 )
587 if self.PortfolioBool:
588 ShareNow[idx] = self.solution[t].ShareFunc(mNrm)
589 else:
590 ShareNow[idx] = self.RiskyShareFixed
591 self.controls["cNrm"] = cNrmNow
592 self.controls["Share"] = ShareNow
593 self.MPCnow = MPCnow
595 def check_conditions(self, verbose=None):
596 raise NotImplementedError()
598 def calc_limiting_values(self):
599 raise NotImplementedError()
602# This is to preserve compatibility with old name
603RiskyAssetConsumerType = IndShockRiskyAssetConsumerType
606###############################################################################
609def solve_one_period_ConsPortChoice(
610 solution_next,
611 ShockDstn,
612 IncShkDstn,
613 RiskyDstn,
614 LivPrb,
615 DiscFac,
616 CRRA,
617 Rfree,
618 PermGroFac,
619 BoroCnstArt,
620 aXtraGrid,
621 ShareGrid,
622 ShareLimit,
623 vFuncBool,
624 IndepDstnBool,
625):
626 """
627 Solve one period of a consumption-saving problem with portfolio allocation
628 between a riskless and risky asset. This function handles only the most
629 fundamental portfolio choice problem: frictionless reallocation of the
630 portfolio each period as a continuous choice.
632 Parameters
633 ----------
634 solution_next : PortfolioSolution
635 Solution to next period's problem.
636 ShockDstn : Distribution
637 Joint distribution of permanent income shocks, transitory income shocks,
638 and risky returns. This is only used if the input IndepDstnBool is False,
639 indicating that income and return distributions can't be assumed to be
640 independent.
641 IncShkDstn : Distribution
642 Discrete distribution of permanent income shocks and transitory income
643 shocks. This is only used if the input IndepDstnBool is True, indicating
644 that income and return distributions are independent.
645 RiskyDstn : Distribution
646 Distribution of risky asset returns. This is only used if the input
647 IndepDstnBool is True, indicating that income and return distributions
648 are independent.
649 LivPrb : float
650 Survival probability; likelihood of being alive at the beginning of
651 the succeeding period.
652 DiscFac : float
653 Intertemporal discount factor for future utility.
654 CRRA : float
655 Coefficient of relative risk aversion.
656 Rfree : float
657 Risk free interest factor on end-of-period assets.
658 PermGroFac : float
659 Expected permanent income growth factor at the end of this period.
660 BoroCnstArt: float or None
661 Borrowing constraint for the minimum allowable assets to end the
662 period with. In this model, it is *required* to be zero.
663 aXtraGrid: np.array
664 Array of "extra" end-of-period asset values-- assets above the
665 absolute minimum acceptable level.
666 ShareGrid : np.array
667 Array of risky portfolio shares on which to define the interpolation
668 of the consumption function when Share is fixed. Also used when the
669 risky share choice is specified as discrete rather than continuous.
670 ShareLimit : float
671 Limiting lower bound of risky portfolio share as mNrm approaches infinity.
672 vFuncBool: boolean
673 An indicator for whether the value function should be computed and
674 included in the reported solution.
675 IndepDstnBool : bool
676 Indicator for whether the income and risky return distributions are in-
677 dependent of each other, which can speed up the expectations step.
679 Returns
680 -------
681 solution_now : PortfolioSolution
682 Solution to this period's problem.
684 :meta private:
685 """
686 # Make sure the individual is liquidity constrained. Allowing a consumer to
687 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
688 if BoroCnstArt != 0.0:
689 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
691 # Define the current period utility function and effective discount factor
692 uFunc = UtilityFuncCRRA(CRRA)
693 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
695 # Unpack next period's solution for easier access
696 vPfunc_next = solution_next.vPfunc
697 vFunc_next = solution_next.vFunc
699 # Set a flag for whether the natural borrowing constraint is zero, which
700 # depends on whether the smallest transitory income shock is zero
701 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
703 # Prepare to calculate end-of-period marginal values by creating an array
704 # of market resources that the agent could have next period, considering
705 # the grid of end-of-period assets and the distribution of shocks he might
706 # experience next period.
708 # Unpack the risky return shock distribution
709 Risky_next = RiskyDstn.atoms
710 RiskyMax = np.max(Risky_next)
711 RiskyMin = np.min(Risky_next)
713 # Perform an alternate calculation of the absolute patience factor when
714 # returns are risky. This uses the Merton-Samuelson limiting risky share,
715 # which is what's relevant as mNrm goes to infinity.
716 def calc_Radj(R):
717 Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree
718 return Rport ** (1.0 - CRRA)
720 R_adj = expected(calc_Radj, RiskyDstn)[0]
721 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
722 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
724 # Also perform an alternate calculation for human wealth under risky returns
725 def calc_hNrm(S):
726 Risky = S["Risky"]
727 PermShk = S["PermShk"]
728 TranShk = S["TranShk"]
729 G = PermGroFac * PermShk
730 Rport = ShareLimit * Risky + (1.0 - ShareLimit) * Rfree
731 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
732 return hNrm
734 # This correctly accounts for risky returns and risk aversion
735 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
737 # This basic equation works if there's no correlation among shocks
738 # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm)
740 # Define the terms for the limiting linear consumption function as m gets very big
741 cFuncLimitIntercept = MPCminNow * hNrmNow
742 cFuncLimitSlope = MPCminNow
744 # bNrm represents R*a, balances after asset return shocks but before income.
745 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
746 if BoroCnstNat_iszero:
747 aNrmGrid = aXtraGrid
748 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
749 else:
750 # Add an asset point at exactly zero
751 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
752 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
754 # Get grid and shock sizes, for easier indexing
755 aNrmCount = aNrmGrid.size
756 ShareCount = ShareGrid.size
758 # If the income shock distribution is independent from the risky return distribution,
759 # then taking end-of-period expectations can proceed in a two part process: First,
760 # construct an "intermediate" value function by integrating out next period's income
761 # shocks, *then* compute end-of-period expectations by integrating out return shocks.
762 # This method is lengthy to code, but can be significantly faster.
763 if IndepDstnBool:
764 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
765 bNrmNext = bNrmGrid
767 # Define functions that are used internally to evaluate future realizations
768 def calc_mNrm_next(S, b):
769 """
770 Calculate future realizations of market resources mNrm from the income
771 shock distribution S and normalized bank balances b.
772 """
773 return b / (S["PermShk"] * PermGroFac) + S["TranShk"]
775 def calc_dvdm_next(S, b):
776 """
777 Evaluate realizations of marginal value of market resources next period,
778 based on the income distribution S and values of bank balances bNrm
779 """
780 mNrm_next = calc_mNrm_next(S, b)
781 G = S["PermShk"] * PermGroFac
782 dvdm_next = G ** (-CRRA) * vPfunc_next(mNrm_next)
783 return dvdm_next
785 # Calculate end-of-period marginal value of assets and shares at each point
786 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
787 # values across income and risky return shocks.
789 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
790 dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext))
791 dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0))
793 dvdbNvrsFunc_intermed = LinearInterp(bNrmGrid, dvdbNvrs_intermed)
794 dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA)
796 # The intermediate marginal value of risky portfolio share is zero in this
797 # model because risky share is flexible: we can just change it next period,
798 # so there is no marginal value of Share once the return is realized.
799 dvdsFunc_intermed = ConstantFunction(0.0) # all zeros
801 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
802 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
804 # Define functions for calculating end-of-period marginal value
805 def calc_EndOfPrd_dvda(R, a, z):
806 """
807 Compute end-of-period marginal value of assets at values a, conditional
808 on risky asset return R and risky share z.
809 """
810 # Calculate future realizations of bank balances bNrm
811 Rxs = R - Rfree # Excess returns
812 Rport = Rfree + z * Rxs # Portfolio return
813 bNrm_next = Rport * a
815 # Calculate and return dvda
816 EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next)
817 return EndOfPrd_dvda
819 def calc_EndOfPrd_dvds(R, a, z):
820 """
821 Compute end-of-period marginal value of risky share at values a, conditional
822 on risky asset return S and risky share z.
823 """
824 # Calculate future realizations of bank balances bNrm
825 Rxs = R - Rfree # Excess returns
826 Rport = Rfree + z * Rxs # Portfolio return
827 bNrm_next = Rport * a
829 # Calculate and return dvds (second term is all zeros)
830 EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(
831 bNrm_next
832 )
833 return EndOfPrd_dvds
835 TempDstn = RiskyDstn # relabeling for below
837 # Evaluate realizations of value and marginal value after asset returns are realized
839 # Calculate end-of-period marginal value of risky portfolio share by taking expectations
840 EndOfPrd_dvds = DiscFacEff * expected(
841 calc_EndOfPrd_dvds, RiskyDstn, args=(aNrmNow, ShareNext)
842 )
844 # Make the end-of-period value function if the value function is requested
845 if vFuncBool:
847 def calc_v_intermed(S, b):
848 """
849 Calculate "intermediate" value from next period's bank balances, the
850 income shocks S, and the risky asset share.
851 """
852 mNrm_next = calc_mNrm_next(S, b)
853 v_next = vFunc_next(mNrm_next)
854 v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
855 return v_intermed
857 # Calculate intermediate value by taking expectations over income shocks
858 v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext))
860 # Construct the "intermediate value function" for this period
861 vNvrs_intermed = uFunc.inv(v_intermed)
862 vNvrsFunc_intermed = LinearInterp(bNrmGrid, vNvrs_intermed)
863 vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA)
865 def calc_EndOfPrd_v(S, a, z):
866 # Calculate future realizations of bank balances bNrm
867 Rxs = S - Rfree
868 Rport = Rfree + z * Rxs
869 bNrm_next = Rport * a
871 EndOfPrd_v = vFunc_intermed(bNrm_next)
872 return EndOfPrd_v
874 # Calculate end-of-period value by taking expectations
875 EndOfPrd_v = DiscFacEff * expected(
876 calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext)
877 )
878 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
880 # Now make an end-of-period value function over aNrm and Share
881 EndOfPrd_vNvrsFunc = BilinearInterp(EndOfPrd_vNvrs, aNrmGrid, ShareGrid)
882 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
883 # This will be used later to make the value function for this period
885 # If the income shock distribution and risky return distribution are *NOT*
886 # independent, then computation of end-of-period expectations are simpler in
887 # code, but might take longer to execute
888 else:
889 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
890 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
892 # Define functions that are used internally to evaluate future realizations
893 def calc_mNrm_next(S, a, z):
894 """
895 Calculate future realizations of market resources mNrm from the shock
896 distribution S, normalized end-of-period assets a, and risky share z.
897 """
898 # Calculate future realizations of bank balances bNrm
899 Rxs = S["Risky"] - Rfree
900 Rport = Rfree + z * Rxs
901 bNrm_next = Rport * a
902 mNrm_next = bNrm_next / (S["PermShk"] * PermGroFac) + S["TranShk"]
903 return mNrm_next
905 def calc_EndOfPrd_dvdx(S, a, z):
906 """
907 Evaluate end-of-period marginal value of assets and risky share based
908 on the shock distribution S, normalized end-of-period assets a, and
909 risky share z.
910 """
911 mNrm_next = calc_mNrm_next(S, a, z)
912 Rxs = S["Risky"] - Rfree
913 Rport = Rfree + z * Rxs
914 dvdm_next = vPfunc_next(mNrm_next)
915 # No marginal value of Share if it's a free choice!
916 dvds_next = np.zeros_like(mNrm_next)
918 EndOfPrd_dvda = Rport * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
919 EndOfPrd_dvds = (
920 Rxs * a * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
921 + (S["PermShk"] * PermGroFac) ** (1 - CRRA) * dvds_next
922 )
924 return EndOfPrd_dvda, EndOfPrd_dvds
926 def calc_EndOfPrd_v(S, a, z):
927 """
928 Evaluate end-of-period value, based on the shock distribution S, values
929 of bank balances bNrm, and values of the risky share z.
930 """
931 mNrm_next = calc_mNrm_next(S, a, z)
932 v_next = vFunc_next(mNrm_next)
933 EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
934 return EndOfPrd_v
936 calc_EndOfPrd_dvda = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[0]
937 calc_EndOfPrd_dvds = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[1]
938 TempDstn = ShockDstn
940 # Evaluate realizations of value and marginal value after asset returns are realized
942 # Calculate end-of-period marginal value of assets and risky share by taking expectations
943 EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected(
944 calc_EndOfPrd_dvdx, ShockDstn, args=(aNrmNow, ShareNext)
945 )
946 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
948 # Construct the end-of-period value function if requested
949 if vFuncBool:
950 # Calculate end-of-period value, its derivative, and their pseudo-inverse
951 EndOfPrd_v = DiscFacEff * expected(
952 calc_EndOfPrd_v, ShockDstn, args=(aNrmNow, ShareNext)
953 )
954 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
956 # value transformed through inverse utility
957 EndOfPrd_vNvrsP = EndOfPrd_dvda * uFunc.derinv(EndOfPrd_v, order=(0, 1))
959 # Construct the end-of-period value function
960 EndOfPrd_vNvrsFunc_by_Share = []
961 for j in range(ShareCount):
962 EndOfPrd_vNvrsFunc_by_Share.append(
963 CubicInterp(
964 aNrmNow[:, j], EndOfPrd_vNvrs[:, j], EndOfPrd_vNvrsP[:, j]
965 )
966 )
967 EndOfPrd_vNvrsFunc = LinearInterpOnInterp1D(
968 EndOfPrd_vNvrsFunc_by_Share, ShareGrid
969 )
970 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
972 # Now find the optimal (continuous) risky share on [0,1] by solving the first
973 # order condition EndOfPrd_dvds == 0.
974 FOC_s = EndOfPrd_dvds # Relabel for convenient typing
976 # If agent wants to put more than 100% into risky asset, he is constrained.
977 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
978 constrained_top = FOC_s[:, -1] > 0.0
979 constrained_bot = FOC_s[:, 0] < 0.0
980 constrained = np.logical_or(constrained_top, constrained_bot)
981 a_idx = np.arange(aNrmCount)
983 # For each value of aNrm, find the value of Share such that FOC_s == 0
984 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0)
985 share_idx = np.argmax(crossing, axis=1)
987 for k in range(3):
988 # This represents the index of the segment of the share grid where dvds flips
989 # from positive to negative, indicating that there's a zero *on* the segment.
990 # The exception is for aNrm values that are flagged as constrained, for which
991 # there will be no crossing point and we can just use the boundary value.
993 # Now that we have a *range* for the location of the optimal share, we can
994 # do a refined search for the optimal share at each aNrm value where there
995 # is an interior solution (not constrained). We now make a refined ShareGrid
996 # that has *different* values on it for each aNrm value.
997 bot_s = ShareNext[a_idx, share_idx]
998 top_s = ShareNext[a_idx, share_idx + 1]
999 for j in range(aNrmCount):
1000 if constrained[j]:
1001 continue
1002 ShareNext[j, :] = np.linspace(bot_s[j], top_s[j], ShareCount)
1004 # Now evaluate end-of-period marginal value of risky share on the refined grid
1005 EndOfPrd_dvds = DiscFacEff * expected(
1006 calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)
1007 )
1008 these = np.logical_not(constrained)
1009 FOC_s[these, :] = EndOfPrd_dvds[these, :] # Relabel for convenient typing
1011 # Look for "crossing points" again
1012 crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0)
1013 share_idx[these] = np.argmax(crossing, axis=1)
1015 # Calculate end-of-period marginal value of assets on the refined grid
1016 EndOfPrd_dvda = DiscFacEff * expected(
1017 calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)
1018 )
1019 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
1021 # Calculate the fractional distance between those share gridpoints where the
1022 # zero should be found, assuming a linear function; call it alpha
1023 bot_s = ShareNext[a_idx, share_idx]
1024 top_s = ShareNext[a_idx, share_idx + 1]
1025 bot_f = FOC_s[a_idx, share_idx]
1026 top_f = FOC_s[a_idx, share_idx + 1]
1027 bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx]
1028 top_c = EndOfPrd_dvdaNvrs[a_idx, share_idx + 1]
1029 alpha = 1.0 - top_f / (top_f - bot_f)
1031 # Calculate the continuous optimal risky share and optimal consumption
1032 Share_now = (1.0 - alpha) * bot_s + alpha * top_s
1033 cNrm_now = (1.0 - alpha) * bot_c + alpha * top_c
1035 # If agent wants to put more than 100% into risky asset, he is constrained.
1036 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
1037 constrained_top = FOC_s[:, -1] > 0.0
1038 constrained_bot = FOC_s[:, 0] < 0.0
1040 # Apply the constraints to both risky share and consumption (but lower
1041 # constraint should never be relevant)
1042 Share_now[constrained_top] = 1.0
1043 Share_now[constrained_bot] = 0.0
1044 cNrm_now[constrained_top] = EndOfPrd_dvdaNvrs[constrained_top, -1]
1045 cNrm_now[constrained_bot] = EndOfPrd_dvdaNvrs[constrained_bot, 0]
1047 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
1048 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
1049 # can't depend on the risky share if it doesn't meaningfully exist. Apply
1050 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
1051 if not BoroCnstNat_iszero:
1052 Share_now[0] = 1.0
1053 cNrm_now[0] = EndOfPrd_dvdaNvrs[0, -1]
1055 # Construct functions characterizing the solution for this period
1057 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
1058 # then construct the consumption function when the agent can adjust his share
1059 mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0)
1060 cNrm_now = np.insert(cNrm_now, 0, 0.0)
1061 cFunc_now = LinearInterp(mNrm_now, cNrm_now, cFuncLimitIntercept, cFuncLimitSlope)
1063 # Construct the marginal value (of mNrm) function
1064 vPfunc_now = MargValueFuncCRRA(cFunc_now, CRRA)
1066 # If the share choice is continuous, just make an ordinary interpolating function
1067 if BoroCnstNat_iszero:
1068 Share_lower_bound = ShareLimit
1069 else:
1070 Share_lower_bound = 1.0
1071 Share_now = np.insert(Share_now, 0, Share_lower_bound)
1072 ShareFunc_now = LinearInterp(mNrm_now, Share_now, ShareLimit, 0.0)
1074 # Add the value function if requested
1075 if vFuncBool:
1076 # Create the value functions for this period, defined over market resources
1077 # mNrm when agent can adjust his portfolio, and over market resources and
1078 # fixed share when agent can not adjust his portfolio.
1080 # Construct the value function
1081 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
1082 cNrm_temp = cFunc_now(mNrm_temp)
1083 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # Fix tiny violations
1084 Share_temp = ShareFunc_now(mNrm_temp)
1085 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
1086 vNvrs_temp = uFunc.inv(v_temp)
1087 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
1088 vNvrsFunc = CubicInterp(
1089 np.insert(mNrm_temp, 0, 0.0), # x_list
1090 np.insert(vNvrs_temp, 0, 0.0), # f_list
1091 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
1092 )
1093 # Re-curve the pseudo-inverse value function
1094 vFunc_now = ValueFuncCRRA(vNvrsFunc, CRRA)
1096 else: # If vFuncBool is False, fill in dummy values
1097 vFunc_now = NullFunc()
1099 # Package and return the solution
1100 solution_now = ConsumerSolution(
1101 cFunc=cFunc_now,
1102 vPfunc=vPfunc_now,
1103 vFunc=vFunc_now,
1104 hNrm=hNrmNow,
1105 MPCmin=MPCminNow,
1106 )
1107 solution_now.ShareFunc = ShareFunc_now
1108 return solution_now
1111###############################################################################
1114def solve_one_period_ConsIndShockRiskyAsset(
1115 solution_next,
1116 IncShkDstn,
1117 RiskyDstn,
1118 ShockDstn,
1119 LivPrb,
1120 DiscFac,
1121 Rfree,
1122 CRRA,
1123 PermGroFac,
1124 BoroCnstArt,
1125 aXtraGrid,
1126 RiskyShareFixed,
1127 vFuncBool,
1128 CubicBool,
1129 IndepDstnBool,
1130):
1131 """
1132 Solves one period of a consumption-saving model with idiosyncratic shocks to
1133 permanent and transitory income, with one risky asset and CRRA utility.
1135 Parameters
1136 ----------
1137 solution_next : ConsumerSolution
1138 The solution to next period's one period problem.
1139 IncShkDstn : Distribution
1140 Discrete distribution of permanent income shocks and transitory income
1141 shocks. This is only used if the input IndepDstnBool is True, indicating
1142 that income and return distributions are independent.
1143 RiskyDstn : Distribution
1144 Distribution of risky asset returns. This is only used if the input
1145 IndepDstnBool is True, indicating that income and return distributions
1146 are independent.
1147 ShockDstn : Distribution
1148 Joint distribution of permanent income shocks, transitory income shocks,
1149 and risky returns. This is only used if the input IndepDstnBool is False,
1150 indicating that income and return distributions can't be assumed to be
1151 independent.
1152 LivPrb : float
1153 Survival probability; likelihood of being alive at the beginning of
1154 the succeeding period.
1155 DiscFac : float
1156 Intertemporal discount factor for future utility.
1157 Rfree : float
1158 Risk free interest factor on end-of-period assets.
1159 CRRA : float
1160 Coefficient of relative risk aversion.
1161 PermGroFac : float
1162 Expected permanent income growth factor at the end of this period.
1163 BoroCnstArt: float or None
1164 Borrowing constraint for the minimum allowable assets to end the
1165 period with. If it is less than the natural borrowing constraint,
1166 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
1167 rowing constraint.
1168 aXtraGrid: np.array
1169 Array of "extra" end-of-period asset values-- assets above the
1170 absolute minimum acceptable level.
1171 RiskyShareFixed : float
1172 Fixed fraction of end-of-period assets that are allocated to the risky asset.
1173 vFuncBool: boolean
1174 An indicator for whether the value function should be computed and
1175 included in the reported solution.
1176 CubicBool: boolean
1177 An indicator for whether the solver should use cubic or linear interpolation.
1178 IndepDstnBool : bool
1179 Indicator for whether the income and risky return distributions are in-
1180 dependent of each other, which can speed up the expectations step.
1182 Returns
1183 -------
1184 solution_now : ConsumerSolution
1185 Solution to this period's consumption-saving problem with income risk.
1187 :meta private:
1188 """
1189 # Do a quick validity check; don't want to allow borrowing with risky returns
1190 if BoroCnstArt != 0.0:
1191 raise ValueError("RiskyAssetConsumerType must have BoroCnstArt=0.0!")
1193 # Define the current period utility function and effective discount factor
1194 uFunc = UtilityFuncCRRA(CRRA)
1195 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
1197 # Unpack next period's income shock distribution
1198 ShkPrbsNext = ShockDstn.pmv
1199 PermShkValsNext = ShockDstn.atoms[0]
1200 TranShkValsNext = ShockDstn.atoms[1]
1201 RiskyValsNext = ShockDstn.atoms[2]
1202 PermShkMinNext = np.min(PermShkValsNext)
1203 TranShkMinNext = np.min(TranShkValsNext)
1204 RiskyMinNext = np.min(RiskyValsNext)
1205 RiskyMaxNext = np.max(RiskyValsNext)
1207 # Unpack next period's (marginal) value function
1208 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
1209 vPfuncNext = solution_next.vPfunc
1210 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
1212 # Perform an alternate calculation of the absolute patience factor when returns are risky
1213 def calc_Radj(R):
1214 Rport = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree
1215 return Rport ** (1.0 - CRRA)
1217 R_adj = expected(calc_Radj, RiskyDstn)
1218 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
1219 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
1220 MPCminNow = MPCminNow[0]
1222 # Also perform an alternate calculation for human wealth under risky returns
1223 def calc_hNrm(S):
1224 Risky = S["Risky"]
1225 PermShk = S["PermShk"]
1226 TranShk = S["TranShk"]
1227 G = PermGroFac * PermShk
1228 Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree
1229 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
1230 return hNrm
1232 # This correctly accounts for risky returns and risk aversion
1233 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
1234 hNrmNow = hNrmNow[0]
1236 # The above attempts to pin down the limiting consumption function for this
1237 # model, however it is not clear why it creates bugs, so for now we allow
1238 # for a linear extrapolation beyond the last asset point
1239 cFuncLimitIntercept = MPCminNow * hNrmNow
1240 cFuncLimitSlope = MPCminNow
1242 # Calculate the minimum allowable value of market resources in this period
1243 BoroCnstNat_cand = (
1244 (solution_next.mNrmMin - TranShkValsNext)
1245 * (PermGroFac * PermShkValsNext)
1246 / RiskyValsNext
1247 )
1248 BoroCnstNat = np.max(BoroCnstNat_cand) # Must be at least this
1250 # Set a flag for whether the natural borrowing constraint is zero, which
1251 # depends on whether the smallest transitory income shock is zero
1252 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
1254 # Set the minimum allowable (normalized) market resources based on the natural
1255 # and artificial borrowing constraints
1256 if BoroCnstArt is None:
1257 mNrmMinNow = BoroCnstNat
1258 else:
1259 mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt])
1261 # The MPCmax code is a bit unusual here, and possibly "harmlessly wrong".
1262 # The "worst event" should depend on the risky return factor as well as
1263 # income shocks. However, the natural borrowing constraint is only ever
1264 # relevant in this model when it's zero, so the MPC at mNrm is only relevant
1265 # in the case where risky returns don't matter at all (because a=0).
1267 # Calculate the probability that we get the worst possible income draw
1268 IncNext = PermShkValsNext * TranShkValsNext
1269 WorstIncNext = PermShkMinNext * TranShkMinNext
1270 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
1271 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
1273 # Update the upper bounding MPC as market resources approach the lower bound
1274 temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
1275 MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
1277 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
1278 # or artificial borrowing constraint actually binds
1279 if BoroCnstNat < mNrmMinNow:
1280 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
1281 else:
1282 MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above
1284 # Define the borrowing-constrained consumption function
1285 cFuncNowCnst = LinearInterp(
1286 np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0])
1287 )
1289 # Big methodological split here: whether the income and return distributions are independent.
1290 # Calculation of end-of-period marginal (marginal) value uses different approaches
1291 if IndepDstnBool:
1292 # bNrm represents R*a, balances after asset return shocks but before income.
1293 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
1294 if BoroCnstNat_iszero:
1295 bNrmNow = np.insert(
1296 RiskyMaxNext * aXtraGrid, 0, RiskyMinNext * aXtraGrid[0]
1297 )
1298 aNrmNow = aXtraGrid.copy()
1299 else:
1300 # Add a bank balances point at exactly zero
1301 bNrmNow = RiskyMaxNext * np.insert(aXtraGrid, 0, 0.0)
1302 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1304 # Define local functions for taking future expectations when the interest
1305 # factor *is* independent from the income shock distribution. These go
1306 # from "bank balances" bNrm = R * aNrm to t+1 realizations.
1307 def calc_mNrmNext(S, b):
1308 return b / (PermGroFac * S["PermShk"]) + S["TranShk"]
1310 def calc_vNext(S, b):
1311 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, b))
1313 def calc_vPnext(S, b):
1314 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, b))
1316 def calc_vPPnext(S, b):
1317 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, b))
1319 # Calculate marginal value of bank balances at each gridpoint
1320 vPfacEff = PermGroFac ** (-CRRA)
1321 Intermed_vP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(bNrmNow))
1322 Intermed_vPnvrs = uFunc.derinv(Intermed_vP, order=(1, 0))
1324 if BoroCnstNat_iszero:
1325 Intermed_vPnvrs = np.insert(Intermed_vPnvrs, 0, 0.0)
1326 bNrm_temp = np.insert(bNrmNow, 0, 0.0)
1327 else:
1328 bNrm_temp = bNrmNow.copy()
1330 # If using cubic spline interpolation, also calculate "intermediate"
1331 # marginal marginal value of bank balances
1332 if CubicBool:
1333 vPPfacEff = PermGroFac ** (-CRRA - 1.0)
1334 Intermed_vPP = vPPfacEff * expected(
1335 calc_vPPnext, IncShkDstn, args=(bNrmNow)
1336 )
1337 Intermed_vPnvrsP = Intermed_vPP * uFunc.derinv(Intermed_vP, order=(1, 1))
1338 if BoroCnstNat_iszero:
1339 Intermed_vPnvrsP = np.insert(Intermed_vPnvrsP, 0, Intermed_vPnvrsP[0])
1341 # Make a cubic spline intermediate pseudo-inverse marginal value function
1342 Intermed_vPnvrsFunc = CubicInterp(
1343 bNrm_temp,
1344 Intermed_vPnvrs,
1345 Intermed_vPnvrsP,
1346 lower_extrap=True,
1347 )
1348 Intermed_vPPfunc = MargMargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1349 else:
1350 # Make a linear interpolation intermediate pseudo-inverse marginal value function
1351 Intermed_vPnvrsFunc = LinearInterp(
1352 bNrm_temp, Intermed_vPnvrs, lower_extrap=True
1353 )
1355 # "Recurve" the intermediate pseudo-inverse marginal value function
1356 Intermed_vPfunc = MargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1358 # If the value function is requested, calculate "intermediate" value
1359 if vFuncBool:
1360 vFacEff = PermGroFac ** (1.0 - CRRA)
1361 Intermed_v = vFacEff * expected(calc_vNext, IncShkDstn, args=(bNrmNow))
1362 Intermed_vNvrs = uFunc.inv(Intermed_v)
1363 # value transformed through inverse utility
1364 Intermed_vNvrsP = Intermed_vP * uFunc.derinv(Intermed_v, order=(0, 1))
1365 if BoroCnstNat_iszero:
1366 Intermed_vNvrs = np.insert(Intermed_vNvrs, 0, 0.0)
1367 Intermed_vNvrsP = np.insert(Intermed_vNvrsP, 0, Intermed_vNvrsP[0])
1368 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1370 # Make a cubic spline intermediate pseudo-inverse value function
1371 Intermed_vNvrsFunc = CubicInterp(bNrm_temp, Intermed_vNvrs, Intermed_vNvrsP)
1373 # "Recurve" the intermediate pseudo-inverse value function
1374 Intermed_vFunc = ValueFuncCRRA(Intermed_vNvrsFunc, CRRA)
1376 # We have "intermediate" (marginal) value functions defined over bNrm,
1377 # so now we want to take expectations over Risky realizations at each aNrm.
1379 # Begin by re-defining transition functions for taking expectations, which are all very simple!
1380 Z = RiskyShareFixed # for shorter notation
1382 def calc_bNrmNext(R, a):
1383 Rport = Z * R + (1 - Z) * Rfree
1384 return Rport * a
1386 def calc_vNext(R, a):
1387 return Intermed_vFunc(calc_bNrmNext(R, a))
1389 def calc_vPnext(R, a):
1390 Rport = Z * R + (1 - Z) * Rfree
1391 return Rport * Intermed_vPfunc(calc_bNrmNext(R, a))
1393 def calc_vPPnext(R, a):
1394 Rport = Z * R + (1 - Z) * Rfree
1395 return Rport * Rport * Intermed_vPPfunc(calc_bNrmNext(R, a))
1397 # Calculate end-of-period marginal value of assets at each gridpoint
1398 EndOfPrdvP = DiscFacEff * expected(calc_vPnext, RiskyDstn, args=(aNrmNow))
1400 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1401 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1402 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1404 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1405 if CubicBool:
1406 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1407 EndOfPrdvPP = DiscFacEff * expected(calc_vPPnext, RiskyDstn, args=(aNrmNow))
1408 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1409 MPC = dcda / (dcda + 1.0)
1410 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1412 # Limiting consumption is zero as m approaches mNrmMin
1413 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1414 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1416 # Construct the end-of-period value function if requested
1417 if vFuncBool:
1418 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1419 EndOfPrdv = DiscFacEff * expected(calc_vNext, RiskyDstn, args=(aNrmNow))
1420 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1421 # value transformed through inverse utility
1422 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1424 # Construct the end-of-period value function
1425 if BoroCnstNat_iszero:
1426 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1427 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1428 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1429 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1430 else:
1431 aNrm_temp = aNrmNow.copy()
1433 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1434 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1436 # NON-INDEPENDENT METHOD BEGINS HERE
1437 else:
1438 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
1439 # aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
1440 if BoroCnstNat_iszero:
1441 aNrmNow = aXtraGrid
1442 else:
1443 # Add an asset point at exactly zero
1444 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1446 # Define local functions for taking future expectations when the interest
1447 # factor is *not* independent from the income shock distribution
1448 Z = RiskyShareFixed # for shorter notation
1450 def calc_mNrmNext(S, a):
1451 Risky = S["Risky"]
1452 Rport = Z * Risky + (1 - Z) * Rfree
1453 return Rport / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
1455 def calc_vNext(S, a):
1456 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, a))
1458 def calc_vPnext(S, a):
1459 Risky = S["Risky"]
1460 Rport = Z * Risky + (1 - Z) * Rfree
1461 return Rport * S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a))
1463 def calc_vPPnext(S, a):
1464 Risky = S["Risky"]
1465 Rport = Z * Risky + (1 - Z) * Rfree
1466 return (
1467 (Rport**2)
1468 * S["PermShk"] ** (-CRRA - 1.0)
1469 * vPPfuncNext(calc_mNrmNext(S, a))
1470 )
1472 # Calculate end-of-period marginal value of assets at each gridpoint
1473 vPfacEff = DiscFacEff * PermGroFac ** (-CRRA)
1474 EndOfPrdvP = vPfacEff * expected(calc_vPnext, ShockDstn, args=(aNrmNow))
1476 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1477 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1478 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1480 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1481 if CubicBool:
1482 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1483 vPPfacEff = DiscFacEff * PermGroFac ** (-CRRA - 1.0)
1484 EndOfPrdvPP = vPPfacEff * expected(calc_vPPnext, ShockDstn, args=(aNrmNow))
1485 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1486 MPC = dcda / (dcda + 1.0)
1487 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1489 # Limiting consumption is zero as m approaches mNrmMin
1490 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1491 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1493 # Construct the end-of-period value function if requested
1494 if vFuncBool:
1495 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1496 vFacEff = DiscFacEff * PermGroFac ** (1.0 - CRRA)
1497 EndOfPrdv = vFacEff * expected(calc_vNext, ShockDstn, args=(aNrmNow))
1498 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1499 # value transformed through inverse utility
1500 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1502 # Construct the end-of-period value function
1503 if BoroCnstNat_iszero:
1504 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1505 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1506 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1507 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1508 else:
1509 aNrm_temp = aNrmNow.copy()
1511 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1512 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1514 # Construct the consumption function; this uses the same method whether the
1515 # income distribution is independent from the return distribution
1516 if CubicBool:
1517 # Construct the unconstrained consumption function as a cubic interpolation
1518 cFuncNowUnc = CubicInterp(
1519 m_for_interpolation,
1520 c_for_interpolation,
1521 MPC_for_interpolation,
1522 cFuncLimitIntercept,
1523 cFuncLimitSlope,
1524 )
1525 else:
1526 # Construct the unconstrained consumption function as a linear interpolation
1527 cFuncNowUnc = LinearInterp(
1528 m_for_interpolation,
1529 c_for_interpolation,
1530 cFuncLimitIntercept,
1531 cFuncLimitSlope,
1532 )
1534 # Combine the constrained and unconstrained functions into the true consumption function.
1535 # LowerEnvelope should only be used when BoroCnstArt is True
1536 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1538 # Make the marginal value function and the marginal marginal value function
1539 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1541 # Define this period's marginal marginal value function
1542 if CubicBool:
1543 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1544 else:
1545 vPPfuncNow = NullFunc() # Dummy object
1547 # Construct this period's value function if requested. This version is set
1548 # up for the non-independent distributions, need to write a faster version.
1549 if vFuncBool:
1550 # Compute expected value and marginal value on a grid of market resources
1551 mNrm_temp = mNrmMinNow + aXtraGrid
1552 cNrm_temp = cFuncNow(mNrm_temp)
1553 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # fix tiny errors
1554 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
1555 vP_temp = uFunc.der(cNrm_temp)
1557 # Construct the beginning-of-period value function
1558 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1559 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1560 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1561 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1562 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
1563 # MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1564 vNvrsFuncNow = CubicInterp(mNrm_temp, vNvrs_temp, vNvrsP_temp)
1565 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1566 else:
1567 vFuncNow = NullFunc() # Dummy object
1569 # Create and return this period's solution
1570 solution_now = ConsumerSolution(
1571 cFunc=cFuncNow,
1572 vFunc=vFuncNow,
1573 vPfunc=vPfuncNow,
1574 vPPfunc=vPPfuncNow,
1575 mNrmMin=mNrmMinNow,
1576 hNrm=hNrmNow,
1577 MPCmin=MPCminNow,
1578 MPCmax=MPCmaxEff,
1579 )
1580 solution_now.ShareFunc = ConstantFunction(RiskyShareFixed)
1581 return solution_now