Coverage for HARK/ConsumptionSaving/ConsRiskyAssetModel.py: 95%
490 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +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
596# This is to preserve compatibility with old name
597RiskyAssetConsumerType = IndShockRiskyAssetConsumerType
600###############################################################################
603def solve_one_period_ConsPortChoice(
604 solution_next,
605 ShockDstn,
606 IncShkDstn,
607 RiskyDstn,
608 LivPrb,
609 DiscFac,
610 CRRA,
611 Rfree,
612 PermGroFac,
613 BoroCnstArt,
614 aXtraGrid,
615 ShareGrid,
616 ShareLimit,
617 vFuncBool,
618 IndepDstnBool,
619):
620 """
621 Solve one period of a consumption-saving problem with portfolio allocation
622 between a riskless and risky asset. This function handles only the most
623 fundamental portfolio choice problem: frictionless reallocation of the
624 portfolio each period as a continuous choice.
626 Parameters
627 ----------
628 solution_next : PortfolioSolution
629 Solution to next period's problem.
630 ShockDstn : Distribution
631 Joint distribution of permanent income shocks, transitory income shocks,
632 and risky returns. This is only used if the input IndepDstnBool is False,
633 indicating that income and return distributions can't be assumed to be
634 independent.
635 IncShkDstn : Distribution
636 Discrete distribution of permanent income shocks and transitory income
637 shocks. This is only used if the input IndepDstnBool is True, indicating
638 that income and return distributions are independent.
639 RiskyDstn : Distribution
640 Distribution of risky asset returns. This is only used if the input
641 IndepDstnBool is True, indicating that income and return distributions
642 are independent.
643 LivPrb : float
644 Survival probability; likelihood of being alive at the beginning of
645 the succeeding period.
646 DiscFac : float
647 Intertemporal discount factor for future utility.
648 CRRA : float
649 Coefficient of relative risk aversion.
650 Rfree : float
651 Risk free interest factor on end-of-period assets.
652 PermGroFac : float
653 Expected permanent income growth factor at the end of this period.
654 BoroCnstArt: float or None
655 Borrowing constraint for the minimum allowable assets to end the
656 period with. In this model, it is *required* to be zero.
657 aXtraGrid: np.array
658 Array of "extra" end-of-period asset values-- assets above the
659 absolute minimum acceptable level.
660 ShareGrid : np.array
661 Array of risky portfolio shares on which to define the interpolation
662 of the consumption function when Share is fixed. Also used when the
663 risky share choice is specified as discrete rather than continuous.
664 ShareLimit : float
665 Limiting lower bound of risky portfolio share as mNrm approaches infinity.
666 vFuncBool: boolean
667 An indicator for whether the value function should be computed and
668 included in the reported solution.
669 IndepDstnBool : bool
670 Indicator for whether the income and risky return distributions are in-
671 dependent of each other, which can speed up the expectations step.
673 Returns
674 -------
675 solution_now : PortfolioSolution
676 Solution to this period's problem.
678 :meta private:
679 """
680 # Make sure the individual is liquidity constrained. Allowing a consumer to
681 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
682 if BoroCnstArt != 0.0:
683 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
685 # Define the current period utility function and effective discount factor
686 uFunc = UtilityFuncCRRA(CRRA)
687 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
689 # Unpack next period's solution for easier access
690 vPfunc_next = solution_next.vPfunc
691 vFunc_next = solution_next.vFunc
693 # Set a flag for whether the natural borrowing constraint is zero, which
694 # depends on whether the smallest transitory income shock is zero
695 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
697 # Prepare to calculate end-of-period marginal values by creating an array
698 # of market resources that the agent could have next period, considering
699 # the grid of end-of-period assets and the distribution of shocks he might
700 # experience next period.
702 # Unpack the risky return shock distribution
703 Risky_next = RiskyDstn.atoms
704 RiskyMax = np.max(Risky_next)
705 RiskyMin = np.min(Risky_next)
707 # Perform an alternate calculation of the absolute patience factor when
708 # returns are risky. This uses the Merton-Samuelson limiting risky share,
709 # which is what's relevant as mNrm goes to infinity.
710 def calc_Radj(R):
711 Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree
712 return Rport ** (1.0 - CRRA)
714 R_adj = expected(calc_Radj, RiskyDstn)[0]
715 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
716 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
718 # Also perform an alternate calculation for human wealth under risky returns
719 def calc_hNrm(S):
720 Risky = S["Risky"]
721 PermShk = S["PermShk"]
722 TranShk = S["TranShk"]
723 G = PermGroFac * PermShk
724 Rport = ShareLimit * Risky + (1.0 - ShareLimit) * Rfree
725 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
726 return hNrm
728 # This correctly accounts for risky returns and risk aversion
729 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
731 # This basic equation works if there's no correlation among shocks
732 # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm)
734 # Define the terms for the limiting linear consumption function as m gets very big
735 cFuncLimitIntercept = MPCminNow * hNrmNow
736 cFuncLimitSlope = MPCminNow
738 # bNrm represents R*a, balances after asset return shocks but before income.
739 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
740 if BoroCnstNat_iszero:
741 aNrmGrid = aXtraGrid
742 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
743 else:
744 # Add an asset point at exactly zero
745 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
746 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
748 # Get grid and shock sizes, for easier indexing
749 aNrmCount = aNrmGrid.size
750 ShareCount = ShareGrid.size
752 # If the income shock distribution is independent from the risky return distribution,
753 # then taking end-of-period expectations can proceed in a two part process: First,
754 # construct an "intermediate" value function by integrating out next period's income
755 # shocks, *then* compute end-of-period expectations by integrating out return shocks.
756 # This method is lengthy to code, but can be significantly faster.
757 if IndepDstnBool:
758 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
759 bNrmNext = bNrmGrid
761 # Define functions that are used internally to evaluate future realizations
762 def calc_mNrm_next(S, b):
763 """
764 Calculate future realizations of market resources mNrm from the income
765 shock distribution S and normalized bank balances b.
766 """
767 return b / (S["PermShk"] * PermGroFac) + S["TranShk"]
769 def calc_dvdm_next(S, b):
770 """
771 Evaluate realizations of marginal value of market resources next period,
772 based on the income distribution S and values of bank balances bNrm
773 """
774 mNrm_next = calc_mNrm_next(S, b)
775 G = S["PermShk"] * PermGroFac
776 dvdm_next = G ** (-CRRA) * vPfunc_next(mNrm_next)
777 return dvdm_next
779 # Calculate end-of-period marginal value of assets and shares at each point
780 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
781 # values across income and risky return shocks.
783 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
784 dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext))
785 dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0))
787 dvdbNvrsFunc_intermed = LinearInterp(bNrmGrid, dvdbNvrs_intermed)
788 dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA)
790 # The intermediate marginal value of risky portfolio share is zero in this
791 # model because risky share is flexible: we can just change it next period,
792 # so there is no marginal value of Share once the return is realized.
793 dvdsFunc_intermed = ConstantFunction(0.0) # all zeros
795 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
796 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
798 # Define functions for calculating end-of-period marginal value
799 def calc_EndOfPrd_dvda(R, a, z):
800 """
801 Compute end-of-period marginal value of assets at values a, conditional
802 on risky asset return R and risky share z.
803 """
804 # Calculate future realizations of bank balances bNrm
805 Rxs = R - Rfree # Excess returns
806 Rport = Rfree + z * Rxs # Portfolio return
807 bNrm_next = Rport * a
809 # Calculate and return dvda
810 EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next)
811 return EndOfPrd_dvda
813 def calc_EndOfPrd_dvds(R, a, z):
814 """
815 Compute end-of-period marginal value of risky share at values a, conditional
816 on risky asset return S and risky share z.
817 """
818 # Calculate future realizations of bank balances bNrm
819 Rxs = R - Rfree # Excess returns
820 Rport = Rfree + z * Rxs # Portfolio return
821 bNrm_next = Rport * a
823 # Calculate and return dvds (second term is all zeros)
824 EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(
825 bNrm_next
826 )
827 return EndOfPrd_dvds
829 TempDstn = RiskyDstn # relabeling for below
831 # Evaluate realizations of value and marginal value after asset returns are realized
833 # Calculate end-of-period marginal value of risky portfolio share by taking expectations
834 EndOfPrd_dvds = DiscFacEff * expected(
835 calc_EndOfPrd_dvds, RiskyDstn, args=(aNrmNow, ShareNext)
836 )
838 # Make the end-of-period value function if the value function is requested
839 if vFuncBool:
841 def calc_v_intermed(S, b):
842 """
843 Calculate "intermediate" value from next period's bank balances, the
844 income shocks S, and the risky asset share.
845 """
846 mNrm_next = calc_mNrm_next(S, b)
847 v_next = vFunc_next(mNrm_next)
848 v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
849 return v_intermed
851 # Calculate intermediate value by taking expectations over income shocks
852 v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext))
854 # Construct the "intermediate value function" for this period
855 vNvrs_intermed = uFunc.inv(v_intermed)
856 vNvrsFunc_intermed = LinearInterp(bNrmGrid, vNvrs_intermed)
857 vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA)
859 def calc_EndOfPrd_v(S, a, z):
860 # Calculate future realizations of bank balances bNrm
861 Rxs = S - Rfree
862 Rport = Rfree + z * Rxs
863 bNrm_next = Rport * a
865 EndOfPrd_v = vFunc_intermed(bNrm_next)
866 return EndOfPrd_v
868 # Calculate end-of-period value by taking expectations
869 EndOfPrd_v = DiscFacEff * expected(
870 calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext)
871 )
872 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
874 # Now make an end-of-period value function over aNrm and Share
875 EndOfPrd_vNvrsFunc = BilinearInterp(EndOfPrd_vNvrs, aNrmGrid, ShareGrid)
876 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
877 # This will be used later to make the value function for this period
879 # If the income shock distribution and risky return distribution are *NOT*
880 # independent, then computation of end-of-period expectations are simpler in
881 # code, but might take longer to execute
882 else:
883 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
884 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
886 # Define functions that are used internally to evaluate future realizations
887 def calc_mNrm_next(S, a, z):
888 """
889 Calculate future realizations of market resources mNrm from the shock
890 distribution S, normalized end-of-period assets a, and risky share z.
891 """
892 # Calculate future realizations of bank balances bNrm
893 Rxs = S["Risky"] - Rfree
894 Rport = Rfree + z * Rxs
895 bNrm_next = Rport * a
896 mNrm_next = bNrm_next / (S["PermShk"] * PermGroFac) + S["TranShk"]
897 return mNrm_next
899 def calc_EndOfPrd_dvdx(S, a, z):
900 """
901 Evaluate end-of-period marginal value of assets and risky share based
902 on the shock distribution S, normalized end-of-period assets a, and
903 risky share z.
904 """
905 mNrm_next = calc_mNrm_next(S, a, z)
906 Rxs = S["Risky"] - Rfree
907 Rport = Rfree + z * Rxs
908 dvdm_next = vPfunc_next(mNrm_next)
909 # No marginal value of Share if it's a free choice!
910 dvds_next = np.zeros_like(mNrm_next)
912 EndOfPrd_dvda = Rport * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
913 EndOfPrd_dvds = (
914 Rxs * a * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next
915 + (S["PermShk"] * PermGroFac) ** (1 - CRRA) * dvds_next
916 )
918 return EndOfPrd_dvda, EndOfPrd_dvds
920 def calc_EndOfPrd_v(S, a, z):
921 """
922 Evaluate end-of-period value, based on the shock distribution S, values
923 of bank balances bNrm, and values of the risky share z.
924 """
925 mNrm_next = calc_mNrm_next(S, a, z)
926 v_next = vFunc_next(mNrm_next)
927 EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next
928 return EndOfPrd_v
930 calc_EndOfPrd_dvda = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[0]
931 calc_EndOfPrd_dvds = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[1]
932 TempDstn = ShockDstn
934 # Evaluate realizations of value and marginal value after asset returns are realized
936 # Calculate end-of-period marginal value of assets and risky share by taking expectations
937 EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected(
938 calc_EndOfPrd_dvdx, ShockDstn, args=(aNrmNow, ShareNext)
939 )
940 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
942 # Construct the end-of-period value function if requested
943 if vFuncBool:
944 # Calculate end-of-period value, its derivative, and their pseudo-inverse
945 EndOfPrd_v = DiscFacEff * expected(
946 calc_EndOfPrd_v, ShockDstn, args=(aNrmNow, ShareNext)
947 )
948 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
950 # value transformed through inverse utility
951 EndOfPrd_vNvrsP = EndOfPrd_dvda * uFunc.derinv(EndOfPrd_v, order=(0, 1))
953 # Construct the end-of-period value function
954 EndOfPrd_vNvrsFunc_by_Share = []
955 for j in range(ShareCount):
956 EndOfPrd_vNvrsFunc_by_Share.append(
957 CubicInterp(
958 aNrmNow[:, j], EndOfPrd_vNvrs[:, j], EndOfPrd_vNvrsP[:, j]
959 )
960 )
961 EndOfPrd_vNvrsFunc = LinearInterpOnInterp1D(
962 EndOfPrd_vNvrsFunc_by_Share, ShareGrid
963 )
964 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
966 # Now find the optimal (continuous) risky share on [0,1] by solving the first
967 # order condition EndOfPrd_dvds == 0.
968 FOC_s = EndOfPrd_dvds # Relabel for convenient typing
970 # If agent wants to put more than 100% into risky asset, he is constrained.
971 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
972 constrained_top = FOC_s[:, -1] > 0.0
973 constrained_bot = FOC_s[:, 0] < 0.0
974 constrained = np.logical_or(constrained_top, constrained_bot)
975 a_idx = np.arange(aNrmCount)
977 # For each value of aNrm, find the value of Share such that FOC_s == 0
978 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0)
979 share_idx = np.argmax(crossing, axis=1)
981 for k in range(3):
982 # This represents the index of the segment of the share grid where dvds flips
983 # from positive to negative, indicating that there's a zero *on* the segment.
984 # The exception is for aNrm values that are flagged as constrained, for which
985 # there will be no crossing point and we can just use the boundary value.
987 # Now that we have a *range* for the location of the optimal share, we can
988 # do a refined search for the optimal share at each aNrm value where there
989 # is an interior solution (not constrained). We now make a refined ShareGrid
990 # that has *different* values on it for each aNrm value.
991 bot_s = ShareNext[a_idx, share_idx]
992 top_s = ShareNext[a_idx, share_idx + 1]
993 for j in range(aNrmCount):
994 if constrained[j]:
995 continue
996 ShareNext[j, :] = np.linspace(bot_s[j], top_s[j], ShareCount)
998 # Now evaluate end-of-period marginal value of risky share on the refined grid
999 EndOfPrd_dvds = DiscFacEff * expected(
1000 calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)
1001 )
1002 these = np.logical_not(constrained)
1003 FOC_s[these, :] = EndOfPrd_dvds[these, :] # Relabel for convenient typing
1005 # Look for "crossing points" again
1006 crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0)
1007 share_idx[these] = np.argmax(crossing, axis=1)
1009 # Calculate end-of-period marginal value of assets on the refined grid
1010 EndOfPrd_dvda = DiscFacEff * expected(
1011 calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)
1012 )
1013 EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda)
1015 # Calculate the fractional distance between those share gridpoints where the
1016 # zero should be found, assuming a linear function; call it alpha
1017 bot_s = ShareNext[a_idx, share_idx]
1018 top_s = ShareNext[a_idx, share_idx + 1]
1019 bot_f = FOC_s[a_idx, share_idx]
1020 top_f = FOC_s[a_idx, share_idx + 1]
1021 bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx]
1022 top_c = EndOfPrd_dvdaNvrs[a_idx, share_idx + 1]
1023 alpha = 1.0 - top_f / (top_f - bot_f)
1025 # Calculate the continuous optimal risky share and optimal consumption
1026 Share_now = (1.0 - alpha) * bot_s + alpha * top_s
1027 cNrm_now = (1.0 - alpha) * bot_c + alpha * top_c
1029 # If agent wants to put more than 100% into risky asset, he is constrained.
1030 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
1031 constrained_top = FOC_s[:, -1] > 0.0
1032 constrained_bot = FOC_s[:, 0] < 0.0
1034 # Apply the constraints to both risky share and consumption (but lower
1035 # constraint should never be relevant)
1036 Share_now[constrained_top] = 1.0
1037 Share_now[constrained_bot] = 0.0
1038 cNrm_now[constrained_top] = EndOfPrd_dvdaNvrs[constrained_top, -1]
1039 cNrm_now[constrained_bot] = EndOfPrd_dvdaNvrs[constrained_bot, 0]
1041 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
1042 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
1043 # can't depend on the risky share if it doesn't meaningfully exist. Apply
1044 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
1045 if not BoroCnstNat_iszero:
1046 Share_now[0] = 1.0
1047 cNrm_now[0] = EndOfPrd_dvdaNvrs[0, -1]
1049 # Construct functions characterizing the solution for this period
1051 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
1052 # then construct the consumption function when the agent can adjust his share
1053 mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0)
1054 cNrm_now = np.insert(cNrm_now, 0, 0.0)
1055 cFunc_now = LinearInterp(mNrm_now, cNrm_now, cFuncLimitIntercept, cFuncLimitSlope)
1057 # Construct the marginal value (of mNrm) function
1058 vPfunc_now = MargValueFuncCRRA(cFunc_now, CRRA)
1060 # If the share choice is continuous, just make an ordinary interpolating function
1061 if BoroCnstNat_iszero:
1062 Share_lower_bound = ShareLimit
1063 else:
1064 Share_lower_bound = 1.0
1065 Share_now = np.insert(Share_now, 0, Share_lower_bound)
1066 ShareFunc_now = LinearInterp(mNrm_now, Share_now, ShareLimit, 0.0)
1068 # Add the value function if requested
1069 if vFuncBool:
1070 # Create the value functions for this period, defined over market resources
1071 # mNrm when agent can adjust his portfolio, and over market resources and
1072 # fixed share when agent can not adjust his portfolio.
1074 # Construct the value function
1075 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
1076 cNrm_temp = cFunc_now(mNrm_temp)
1077 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # Fix tiny violations
1078 Share_temp = ShareFunc_now(mNrm_temp)
1079 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp, Share_temp)
1080 vNvrs_temp = uFunc.inv(v_temp)
1081 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
1082 vNvrsFunc = CubicInterp(
1083 np.insert(mNrm_temp, 0, 0.0), # x_list
1084 np.insert(vNvrs_temp, 0, 0.0), # f_list
1085 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
1086 )
1087 # Re-curve the pseudo-inverse value function
1088 vFunc_now = ValueFuncCRRA(vNvrsFunc, CRRA)
1090 else: # If vFuncBool is False, fill in dummy values
1091 vFunc_now = NullFunc()
1093 # Package and return the solution
1094 solution_now = ConsumerSolution(
1095 cFunc=cFunc_now,
1096 vPfunc=vPfunc_now,
1097 vFunc=vFunc_now,
1098 hNrm=hNrmNow,
1099 MPCmin=MPCminNow,
1100 )
1101 solution_now.ShareFunc = ShareFunc_now
1102 return solution_now
1105###############################################################################
1108def solve_one_period_ConsIndShockRiskyAsset(
1109 solution_next,
1110 IncShkDstn,
1111 RiskyDstn,
1112 ShockDstn,
1113 LivPrb,
1114 DiscFac,
1115 Rfree,
1116 CRRA,
1117 PermGroFac,
1118 BoroCnstArt,
1119 aXtraGrid,
1120 RiskyShareFixed,
1121 vFuncBool,
1122 CubicBool,
1123 IndepDstnBool,
1124):
1125 """
1126 Solves one period of a consumption-saving model with idiosyncratic shocks to
1127 permanent and transitory income, with one risky asset and CRRA utility.
1129 Parameters
1130 ----------
1131 solution_next : ConsumerSolution
1132 The solution to next period's one period problem.
1133 IncShkDstn : Distribution
1134 Discrete distribution of permanent income shocks and transitory income
1135 shocks. This is only used if the input IndepDstnBool is True, indicating
1136 that income and return distributions are independent.
1137 RiskyDstn : Distribution
1138 Distribution of risky asset returns. This is only used if the input
1139 IndepDstnBool is True, indicating that income and return distributions
1140 are independent.
1141 ShockDstn : Distribution
1142 Joint distribution of permanent income shocks, transitory income shocks,
1143 and risky returns. This is only used if the input IndepDstnBool is False,
1144 indicating that income and return distributions can't be assumed to be
1145 independent.
1146 LivPrb : float
1147 Survival probability; likelihood of being alive at the beginning of
1148 the succeeding period.
1149 DiscFac : float
1150 Intertemporal discount factor for future utility.
1151 Rfree : float
1152 Risk free interest factor on end-of-period assets.
1153 CRRA : float
1154 Coefficient of relative risk aversion.
1155 PermGroFac : float
1156 Expected permanent income growth factor at the end of this period.
1157 BoroCnstArt: float or None
1158 Borrowing constraint for the minimum allowable assets to end the
1159 period with. If it is less than the natural borrowing constraint,
1160 then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
1161 rowing constraint.
1162 aXtraGrid: np.array
1163 Array of "extra" end-of-period asset values-- assets above the
1164 absolute minimum acceptable level.
1165 RiskyShareFixed : float
1166 Fixed fraction of end-of-period assets that are allocated to the risky asset.
1167 vFuncBool: boolean
1168 An indicator for whether the value function should be computed and
1169 included in the reported solution.
1170 CubicBool: boolean
1171 An indicator for whether the solver should use cubic or linear interpolation.
1172 IndepDstnBool : bool
1173 Indicator for whether the income and risky return distributions are in-
1174 dependent of each other, which can speed up the expectations step.
1176 Returns
1177 -------
1178 solution_now : ConsumerSolution
1179 Solution to this period's consumption-saving problem with income risk.
1181 :meta private:
1182 """
1183 # Do a quick validity check; don't want to allow borrowing with risky returns
1184 if BoroCnstArt != 0.0:
1185 raise ValueError("RiskyAssetConsumerType must have BoroCnstArt=0.0!")
1187 # Define the current period utility function and effective discount factor
1188 uFunc = UtilityFuncCRRA(CRRA)
1189 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
1191 # Unpack next period's income shock distribution
1192 ShkPrbsNext = ShockDstn.pmv
1193 PermShkValsNext = ShockDstn.atoms[0]
1194 TranShkValsNext = ShockDstn.atoms[1]
1195 RiskyValsNext = ShockDstn.atoms[2]
1196 PermShkMinNext = np.min(PermShkValsNext)
1197 TranShkMinNext = np.min(TranShkValsNext)
1198 RiskyMinNext = np.min(RiskyValsNext)
1199 RiskyMaxNext = np.max(RiskyValsNext)
1201 # Unpack next period's (marginal) value function
1202 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
1203 vPfuncNext = solution_next.vPfunc
1204 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
1206 # Perform an alternate calculation of the absolute patience factor when returns are risky
1207 def calc_Radj(R):
1208 Rport = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree
1209 return Rport ** (1.0 - CRRA)
1211 R_adj = expected(calc_Radj, RiskyDstn)
1212 PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA)
1213 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
1214 MPCminNow = MPCminNow[0]
1216 # Also perform an alternate calculation for human wealth under risky returns
1217 def calc_hNrm(S):
1218 Risky = S["Risky"]
1219 PermShk = S["PermShk"]
1220 TranShk = S["TranShk"]
1221 G = PermGroFac * PermShk
1222 Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree
1223 hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm)
1224 return hNrm
1226 # This correctly accounts for risky returns and risk aversion
1227 hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj
1228 hNrmNow = hNrmNow[0]
1230 # The above attempts to pin down the limiting consumption function for this
1231 # model, however it is not clear why it creates bugs, so for now we allow
1232 # for a linear extrapolation beyond the last asset point
1233 cFuncLimitIntercept = MPCminNow * hNrmNow
1234 cFuncLimitSlope = MPCminNow
1236 # Calculate the minimum allowable value of market resources in this period
1237 BoroCnstNat_cand = (
1238 (solution_next.mNrmMin - TranShkValsNext)
1239 * (PermGroFac * PermShkValsNext)
1240 / RiskyValsNext
1241 )
1242 BoroCnstNat = np.max(BoroCnstNat_cand) # Must be at least this
1244 # Set a flag for whether the natural borrowing constraint is zero, which
1245 # depends on whether the smallest transitory income shock is zero
1246 BoroCnstNat_iszero = np.min(IncShkDstn.atoms[1]) == 0.0
1248 # Set the minimum allowable (normalized) market resources based on the natural
1249 # and artificial borrowing constraints
1250 if BoroCnstArt is None:
1251 mNrmMinNow = BoroCnstNat
1252 else:
1253 mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt])
1255 # The MPCmax code is a bit unusual here, and possibly "harmlessly wrong".
1256 # The "worst event" should depend on the risky return factor as well as
1257 # income shocks. However, the natural borrowing constraint is only ever
1258 # relevant in this model when it's zero, so the MPC at mNrm is only relevant
1259 # in the case where risky returns don't matter at all (because a=0).
1261 # Calculate the probability that we get the worst possible income draw
1262 IncNext = PermShkValsNext * TranShkValsNext
1263 WorstIncNext = PermShkMinNext * TranShkMinNext
1264 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
1265 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
1267 # Update the upper bounding MPC as market resources approach the lower bound
1268 temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
1269 MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
1271 # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural
1272 # or artificial borrowing constraint actually binds
1273 if BoroCnstNat < mNrmMinNow:
1274 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
1275 else:
1276 MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above
1278 # Define the borrowing-constrained consumption function
1279 cFuncNowCnst = LinearInterp(
1280 np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0])
1281 )
1283 # Big methodological split here: whether the income and return distributions are independent.
1284 # Calculation of end-of-period marginal (marginal) value uses different approaches
1285 if IndepDstnBool:
1286 # bNrm represents R*a, balances after asset return shocks but before income.
1287 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
1288 if BoroCnstNat_iszero:
1289 bNrmNow = np.insert(
1290 RiskyMaxNext * aXtraGrid, 0, RiskyMinNext * aXtraGrid[0]
1291 )
1292 aNrmNow = aXtraGrid.copy()
1293 else:
1294 # Add a bank balances point at exactly zero
1295 bNrmNow = RiskyMaxNext * np.insert(aXtraGrid, 0, 0.0)
1296 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1298 # Define local functions for taking future expectations when the interest
1299 # factor *is* independent from the income shock distribution. These go
1300 # from "bank balances" bNrm = R * aNrm to t+1 realizations.
1301 def calc_mNrmNext(S, b):
1302 return b / (PermGroFac * S["PermShk"]) + S["TranShk"]
1304 def calc_vNext(S, b):
1305 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, b))
1307 def calc_vPnext(S, b):
1308 return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, b))
1310 def calc_vPPnext(S, b):
1311 return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, b))
1313 # Calculate marginal value of bank balances at each gridpoint
1314 vPfacEff = PermGroFac ** (-CRRA)
1315 Intermed_vP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(bNrmNow))
1316 Intermed_vPnvrs = uFunc.derinv(Intermed_vP, order=(1, 0))
1318 if BoroCnstNat_iszero:
1319 Intermed_vPnvrs = np.insert(Intermed_vPnvrs, 0, 0.0)
1320 bNrm_temp = np.insert(bNrmNow, 0, 0.0)
1321 else:
1322 bNrm_temp = bNrmNow.copy()
1324 # If using cubic spline interpolation, also calculate "intermediate"
1325 # marginal marginal value of bank balances
1326 if CubicBool:
1327 vPPfacEff = PermGroFac ** (-CRRA - 1.0)
1328 Intermed_vPP = vPPfacEff * expected(
1329 calc_vPPnext, IncShkDstn, args=(bNrmNow)
1330 )
1331 Intermed_vPnvrsP = Intermed_vPP * uFunc.derinv(Intermed_vP, order=(1, 1))
1332 if BoroCnstNat_iszero:
1333 Intermed_vPnvrsP = np.insert(Intermed_vPnvrsP, 0, Intermed_vPnvrsP[0])
1335 # Make a cubic spline intermediate pseudo-inverse marginal value function
1336 Intermed_vPnvrsFunc = CubicInterp(
1337 bNrm_temp,
1338 Intermed_vPnvrs,
1339 Intermed_vPnvrsP,
1340 lower_extrap=True,
1341 )
1342 Intermed_vPPfunc = MargMargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1343 else:
1344 # Make a linear interpolation intermediate pseudo-inverse marginal value function
1345 Intermed_vPnvrsFunc = LinearInterp(
1346 bNrm_temp, Intermed_vPnvrs, lower_extrap=True
1347 )
1349 # "Recurve" the intermediate pseudo-inverse marginal value function
1350 Intermed_vPfunc = MargValueFuncCRRA(Intermed_vPnvrsFunc, CRRA)
1352 # If the value function is requested, calculate "intermediate" value
1353 if vFuncBool:
1354 vFacEff = PermGroFac ** (1.0 - CRRA)
1355 Intermed_v = vFacEff * expected(calc_vNext, IncShkDstn, args=(bNrmNow))
1356 Intermed_vNvrs = uFunc.inv(Intermed_v)
1357 # value transformed through inverse utility
1358 Intermed_vNvrsP = Intermed_vP * uFunc.derinv(Intermed_v, order=(0, 1))
1359 if BoroCnstNat_iszero:
1360 Intermed_vNvrs = np.insert(Intermed_vNvrs, 0, 0.0)
1361 Intermed_vNvrsP = np.insert(Intermed_vNvrsP, 0, Intermed_vNvrsP[0])
1362 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1364 # Make a cubic spline intermediate pseudo-inverse value function
1365 Intermed_vNvrsFunc = CubicInterp(bNrm_temp, Intermed_vNvrs, Intermed_vNvrsP)
1367 # "Recurve" the intermediate pseudo-inverse value function
1368 Intermed_vFunc = ValueFuncCRRA(Intermed_vNvrsFunc, CRRA)
1370 # We have "intermediate" (marginal) value functions defined over bNrm,
1371 # so now we want to take expectations over Risky realizations at each aNrm.
1373 # Begin by re-defining transition functions for taking expectations, which are all very simple!
1374 Z = RiskyShareFixed # for shorter notation
1376 def calc_bNrmNext(R, a):
1377 Rport = Z * R + (1 - Z) * Rfree
1378 return Rport * a
1380 def calc_vNext(R, a):
1381 return Intermed_vFunc(calc_bNrmNext(R, a))
1383 def calc_vPnext(R, a):
1384 Rport = Z * R + (1 - Z) * Rfree
1385 return Rport * Intermed_vPfunc(calc_bNrmNext(R, a))
1387 def calc_vPPnext(R, a):
1388 Rport = Z * R + (1 - Z) * Rfree
1389 return Rport * Rport * Intermed_vPPfunc(calc_bNrmNext(R, a))
1391 # Calculate end-of-period marginal value of assets at each gridpoint
1392 EndOfPrdvP = DiscFacEff * expected(calc_vPnext, RiskyDstn, args=(aNrmNow))
1394 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1395 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1396 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1398 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1399 if CubicBool:
1400 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1401 EndOfPrdvPP = DiscFacEff * expected(calc_vPPnext, RiskyDstn, args=(aNrmNow))
1402 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1403 MPC = dcda / (dcda + 1.0)
1404 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1406 # Limiting consumption is zero as m approaches mNrmMin
1407 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1408 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1410 # Construct the end-of-period value function if requested
1411 if vFuncBool:
1412 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1413 EndOfPrdv = DiscFacEff * expected(calc_vNext, RiskyDstn, args=(aNrmNow))
1414 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1415 # value transformed through inverse utility
1416 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1418 # Construct the end-of-period value function
1419 if BoroCnstNat_iszero:
1420 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1421 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1422 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1423 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1424 else:
1425 aNrm_temp = aNrmNow.copy()
1427 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1428 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1430 # NON-INDEPENDENT METHOD BEGINS HERE
1431 else:
1432 # Construct the assets grid by adjusting aXtra by the natural borrowing constraint
1433 # aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
1434 if BoroCnstNat_iszero:
1435 aNrmNow = aXtraGrid
1436 else:
1437 # Add an asset point at exactly zero
1438 aNrmNow = np.insert(aXtraGrid, 0, 0.0)
1440 # Define local functions for taking future expectations when the interest
1441 # factor is *not* independent from the income shock distribution
1442 Z = RiskyShareFixed # for shorter notation
1444 def calc_mNrmNext(S, a):
1445 Risky = S["Risky"]
1446 Rport = Z * Risky + (1 - Z) * Rfree
1447 return Rport / (PermGroFac * S["PermShk"]) * a + S["TranShk"]
1449 def calc_vNext(S, a):
1450 return S["PermShk"] ** (1.0 - CRRA) * vFuncNext(calc_mNrmNext(S, a))
1452 def calc_vPnext(S, a):
1453 Risky = S["Risky"]
1454 Rport = Z * Risky + (1 - Z) * Rfree
1455 return Rport * S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a))
1457 def calc_vPPnext(S, a):
1458 Risky = S["Risky"]
1459 Rport = Z * Risky + (1 - Z) * Rfree
1460 return (
1461 (Rport**2)
1462 * S["PermShk"] ** (-CRRA - 1.0)
1463 * vPPfuncNext(calc_mNrmNext(S, a))
1464 )
1466 # Calculate end-of-period marginal value of assets at each gridpoint
1467 vPfacEff = DiscFacEff * PermGroFac ** (-CRRA)
1468 EndOfPrdvP = vPfacEff * expected(calc_vPnext, ShockDstn, args=(aNrmNow))
1470 # Invert the first order condition to find optimal cNrm from each aNrm gridpoint
1471 cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0))
1472 mNrmNow = cNrmNow + aNrmNow # Endogenous mNrm gridpoints
1474 # Calculate the MPC at each gridpoint if using cubic spline interpolation
1475 if CubicBool:
1476 # Calculate end-of-period marginal marginal value of assets at each gridpoint
1477 vPPfacEff = DiscFacEff * PermGroFac ** (-CRRA - 1.0)
1478 EndOfPrdvPP = vPPfacEff * expected(calc_vPPnext, ShockDstn, args=(aNrmNow))
1479 dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2)
1480 MPC = dcda / (dcda + 1.0)
1481 MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow)
1483 # Limiting consumption is zero as m approaches mNrmMin
1484 c_for_interpolation = np.insert(cNrmNow, 0, 0.0)
1485 m_for_interpolation = np.insert(mNrmNow, 0, BoroCnstNat)
1487 # Construct the end-of-period value function if requested
1488 if vFuncBool:
1489 # Calculate end-of-period value, its derivative, and their pseudo-inverse
1490 vFacEff = DiscFacEff * PermGroFac ** (1.0 - CRRA)
1491 EndOfPrdv = vFacEff * expected(calc_vNext, ShockDstn, args=(aNrmNow))
1492 EndOfPrdvNvrs = uFunc.inv(EndOfPrdv)
1493 # value transformed through inverse utility
1494 EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1))
1496 # Construct the end-of-period value function
1497 if BoroCnstNat_iszero:
1498 EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0)
1499 EndOfPrdvNvrsP = np.insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
1500 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1501 aNrm_temp = np.insert(aNrmNow, 0, BoroCnstNat)
1502 else:
1503 aNrm_temp = aNrmNow.copy()
1505 EndOfPrd_vNvrsFunc = CubicInterp(aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP)
1506 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1508 # Construct the consumption function; this uses the same method whether the
1509 # income distribution is independent from the return distribution
1510 if CubicBool:
1511 # Construct the unconstrained consumption function as a cubic interpolation
1512 cFuncNowUnc = CubicInterp(
1513 m_for_interpolation,
1514 c_for_interpolation,
1515 MPC_for_interpolation,
1516 cFuncLimitIntercept,
1517 cFuncLimitSlope,
1518 )
1519 else:
1520 # Construct the unconstrained consumption function as a linear interpolation
1521 cFuncNowUnc = LinearInterp(
1522 m_for_interpolation,
1523 c_for_interpolation,
1524 cFuncLimitIntercept,
1525 cFuncLimitSlope,
1526 )
1528 # Combine the constrained and unconstrained functions into the true consumption function.
1529 # LowerEnvelope should only be used when BoroCnstArt is True
1530 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst, nan_bool=False)
1532 # Make the marginal value function and the marginal marginal value function
1533 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
1535 # Define this period's marginal marginal value function
1536 if CubicBool:
1537 vPPfuncNow = MargMargValueFuncCRRA(cFuncNow, CRRA)
1538 else:
1539 vPPfuncNow = NullFunc() # Dummy object
1541 # Construct this period's value function if requested. This version is set
1542 # up for the non-independent distributions, need to write a faster version.
1543 if vFuncBool:
1544 # Compute expected value and marginal value on a grid of market resources
1545 mNrm_temp = mNrmMinNow + aXtraGrid
1546 cNrm_temp = cFuncNow(mNrm_temp)
1547 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # fix tiny errors
1548 v_temp = uFunc(cNrm_temp) + EndOfPrd_vFunc(aNrm_temp)
1549 vP_temp = uFunc.der(cNrm_temp)
1551 # Construct the beginning-of-period value function
1552 vNvrs_temp = uFunc.inv(v_temp) # value transformed through inv utility
1553 vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
1554 mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
1555 vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
1556 vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
1557 # MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
1558 vNvrsFuncNow = CubicInterp(mNrm_temp, vNvrs_temp, vNvrsP_temp)
1559 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
1560 else:
1561 vFuncNow = NullFunc() # Dummy object
1563 # Create and return this period's solution
1564 solution_now = ConsumerSolution(
1565 cFunc=cFuncNow,
1566 vFunc=vFuncNow,
1567 vPfunc=vPfuncNow,
1568 vPPfunc=vPPfuncNow,
1569 mNrmMin=mNrmMinNow,
1570 hNrm=hNrmNow,
1571 MPCmin=MPCminNow,
1572 MPCmax=MPCmaxEff,
1573 )
1574 solution_now.ShareFunc = ConstantFunction(RiskyShareFixed)
1575 return solution_now