Coverage for HARK / ConsumptionSaving / ConsRiskyContribModel.py: 98%
508 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-10 06:19 +0000
1"""
2This file contains classes and functions for representing, solving, and simulating
3a consumer type with idiosyncratic shocks to permanent and transitory income,
4who can save in both a risk-free and a risky asset but faces frictions to
5moving funds between them. The agent can only consume out of his risk-free
6asset.
8The model is described in detail in the REMARK:
9https://econ-ark.org/materials/riskycontrib
11.. code:: bibtex
13 @software{mateo_velasquez_giraldo_2021_4977915,
14 author = {Mateo Velásquez-Giraldo},
15 title = {{Mv77/RiskyContrib: A Two-Asset Savings Model with
16 an Income-Contribution Scheme}},
17 month = jun,
18 year = 2021,
19 publisher = {Zenodo},
20 version = {v1.0.1},
21 doi = {10.5281/zenodo.4977915},
22 url = {https://doi.org/10.5281/zenodo.4977915}
23 }
25"""
27import numpy as np
29from HARK import NullFunc # Basic HARK features
30from HARK.ConsumptionSaving.ConsIndShockModel import utility # CRRA utility function
31from HARK.ConsumptionSaving.ConsIndShockModel import (
32 utility_inv, # Inverse CRRA utility function
33)
34from HARK.ConsumptionSaving.ConsIndShockModel import (
35 utilityP, # CRRA marginal utility function
36)
37from HARK.ConsumptionSaving.ConsIndShockModel import (
38 utilityP_inv, # Inverse CRRA marginal utility function
39)
40from HARK.Calibration.Income.IncomeProcesses import (
41 construct_lognormal_income_process_unemployment,
42 get_PermShkDstn_from_IncShkDstn,
43 get_TranShkDstn_from_IncShkDstn,
44)
45from HARK.Calibration.Assets.AssetProcesses import (
46 make_lognormal_RiskyDstn,
47 combine_IncShkDstn_and_RiskyDstn,
48 calc_ShareLimit_for_CRRA,
49)
50from HARK.ConsumptionSaving.ConsIndShockModel import (
51 init_lifecycle,
52 make_lognormal_kNrm_init_dstn,
53 make_lognormal_pLvl_init_dstn,
54)
55from HARK.ConsumptionSaving.ConsRiskyAssetModel import (
56 RiskyAssetConsumerType,
57 init_risky_asset,
58 make_AdjustDstn,
59)
60from HARK.distributions import expected
61from HARK.interpolation import BilinearInterp # 2D interpolator
62from HARK.interpolation import (
63 ConstantFunction, # Interpolator-like class that returns constant value
64)
65from HARK.interpolation import (
66 IdentityFunction, # Interpolator-like class that returns one of its arguments
67)
68from HARK.interpolation import LinearInterp # Piecewise linear interpolation
69from HARK.interpolation import TrilinearInterp # 3D interpolator
70from HARK.interpolation import DiscreteInterp, MargValueFuncCRRA, ValueFuncCRRA
71from HARK.metric import MetricObject
72from HARK.utilities import make_grid_exp_mult, make_assets_grid
74###############################################################################
77def make_bounded_ShareGrid(ShareCount, ShareMax):
78 """
79 Make a uniformly spaced grid on the unit interval, representing shares
80 contributed toward the risky asset.
82 Parameters
83 ----------
84 ShareCount : int
85 Number of points in the grid.
86 ShareMax : float
87 Highest risky fraction allowed.
89 Returns
90 -------
91 ShareGrid : np.array
92 """
93 ShareGrid = np.linspace(0.0, ShareMax, ShareCount)
94 return ShareGrid
97def make_simple_dGrid(dCount):
98 """
99 Make a uniformly spaced grid on the unit interval, representing rebalancing rates.
101 Parameters
102 ----------
103 dCount : int
104 Number of points in the grid.
106 Returns
107 -------
108 dGrid : np.array
109 """
110 dGrid = np.linspace(0.0, 1.0, dCount)
111 return dGrid
114def make_nNrm_grid(nNrmMin, nNrmMax, nNrmCount, nNrmNestFac):
115 """
116 Creates the agent's illiquid assets grid by constructing a multi-exponentially
117 spaced grid of nNrm values.
119 Parameters
120 ----------
121 nNrmMin : float
122 Minimum value in the illiquid assets grid.
123 nNrmMax : float
124 Maximum value in the illiquid assets grid.
125 nNrmCount : float
126 Number of gridpoints in the illiquid assets grid.
127 nNrmNestFac : int
128 Degree of exponential nesting for illiquid assets.
130 Returns
131 -------
132 nNrmGrid : np.array
133 Constructed grid of illiquid assets.
134 """
135 nNrmGrid = make_grid_exp_mult(
136 ming=nNrmMin, maxg=nNrmMax, ng=nNrmCount, timestonest=nNrmNestFac
137 )
138 return nNrmGrid
141def make_mNrm_grid(mNrmMin, mNrmMax, mNrmCount, mNrmNestFac):
142 """
143 Creates the agent's liquid assets grid by constructing a multi-exponentially
144 spaced grid of mNrm values.
146 Parameters
147 ----------
148 mNrmMin : float
149 Minimum value in the liquid assets grid.
150 mNrmMax : float
151 Maximum value in the liquid assets grid.
152 mNrmCount : float
153 Number of gridpoints in the liquid assets grid.
154 mNrmNestFac : int
155 Degree of exponential nesting for liquid assets.
157 Returns
158 -------
159 mNrmGrid : np.array
160 Constructed grid of liquid assets.
161 """
162 mNrmGrid = make_grid_exp_mult(
163 ming=mNrmMin, maxg=mNrmMax, ng=mNrmCount, timestonest=mNrmNestFac
164 )
165 return mNrmGrid
168def make_solution_terminal_risky_contrib(CRRA, WithdrawTax):
169 """
170 Solves the terminal period. The solution is trivial.
171 Cns: agent will consume all of his liquid resources.
172 Sha: irrelevant as there is no "next" period.
173 Reb: agent will shift all of his resources to the risk-free asset.
175 Parameters
176 ----------
177 CRRA : float
178 Coefficient of relative risk aversion.
179 WithdrawTax : float
180 Tax penalty for withdrawing from the risky asset.
182 Returns
183 -------
184 solution_terminal : RiskyContribSolution
185 Terminal period solution object
186 """
188 # Construct the terminal solution backwards.
190 # Start with the consumption stage. All liquid resources are consumed.
191 cFunc_term = IdentityFunction(i_dim=0, n_dims=3)
192 vFunc_Cns_term = ValueFuncCRRA(cFunc_term, CRRA=CRRA)
193 # Marginal values
194 dvdmFunc_Cns_term = MargValueFuncCRRA(cFunc_term, CRRA=CRRA)
195 dvdnFunc_Cns_term = ConstantFunction(0.0)
196 dvdsFunc_Cns_term = ConstantFunction(0.0)
198 Cns_stage_sol = RiskyContribCnsSolution(
199 # Consumption stage
200 vFunc=vFunc_Cns_term,
201 cFunc=cFunc_term,
202 dvdmFunc=dvdmFunc_Cns_term,
203 dvdnFunc=dvdnFunc_Cns_term,
204 dvdsFunc=dvdsFunc_Cns_term,
205 )
207 # Share stage
209 # It's irrelevant because there is no future period. Set share to 0.
210 # Create a dummy 2-d consumption function to get value function and marginal
211 c2d = IdentityFunction(i_dim=0, n_dims=2)
212 Sha_stage_sol = RiskyContribShaSolution(
213 # Adjust
214 vFunc_Adj=ValueFuncCRRA(c2d, CRRA=CRRA),
215 ShareFunc_Adj=ConstantFunction(0.0),
216 dvdmFunc_Adj=MargValueFuncCRRA(c2d, CRRA=CRRA),
217 dvdnFunc_Adj=ConstantFunction(0.0),
218 # Fixed
219 vFunc_Fxd=vFunc_Cns_term,
220 ShareFunc_Fxd=IdentityFunction(i_dim=2, n_dims=3),
221 dvdmFunc_Fxd=dvdmFunc_Cns_term,
222 dvdnFunc_Fxd=dvdnFunc_Cns_term,
223 dvdsFunc_Fxd=dvdsFunc_Cns_term,
224 )
226 # Rebalancing stage
228 # Adjusting agent:
229 # Withdraw everything from the pension fund and consume everything
230 dfracFunc_Adj_term = ConstantFunction(-1.0)
232 # Find the withdrawal penalty. If it is time-varying, assume it takes
233 # the same value as in the last non-terminal period
234 if type(WithdrawTax) is list:
235 WithdrawTax = WithdrawTax[-1]
237 # Value and marginal value function of the adjusting agent
238 vFunc_Reb_Adj_term = ValueFuncCRRA(lambda m, n: m + n / (1 + WithdrawTax), CRRA)
239 dvdmFunc_Reb_Adj_term = MargValueFuncCRRA(
240 lambda m, n: m + n / (1 + WithdrawTax), CRRA
241 )
242 # A marginal unit of n will be withdrawn and put into m. Then consumed.
243 dvdnFunc_Reb_Adj_term = lambda m, n: dvdmFunc_Reb_Adj_term(m, n) / (1 + WithdrawTax)
245 Reb_stage_sol = RiskyContribRebSolution(
246 # Rebalancing stage
247 vFunc_Adj=vFunc_Reb_Adj_term,
248 dfracFunc_Adj=dfracFunc_Adj_term,
249 dvdmFunc_Adj=dvdmFunc_Reb_Adj_term,
250 dvdnFunc_Adj=dvdnFunc_Reb_Adj_term,
251 # Adjusting stage
252 vFunc_Fxd=vFunc_Cns_term,
253 dfracFunc_Fxd=ConstantFunction(0.0),
254 dvdmFunc_Fxd=dvdmFunc_Cns_term,
255 dvdnFunc_Fxd=dvdnFunc_Cns_term,
256 dvdsFunc_Fxd=dvdsFunc_Cns_term,
257 )
259 # Construct the terminal period solution
260 solution_terminal = RiskyContribSolution(
261 Reb_stage_sol, Sha_stage_sol, Cns_stage_sol
262 )
263 return solution_terminal
266###############################################################################
268# %% Classes for RiskyContrib type solution objects
271# Class for asset adjustment stage solution
272class RiskyContribRebSolution(MetricObject):
273 """
274 A class for representing the solution to the asset-rebalancing stage of
275 the 'RiskyContrib' model.
277 Parameters
278 ----------
279 vFunc_Adj : ValueFunc2D
280 Stage value function over normalized liquid resources and normalized
281 iliquid resources when the agent is able to adjust his portfolio.
282 dfracFunc_Adj : Interp2D
283 Deposit function over normalized liquid resources and normalized
284 iliquid resources when the agent is able to adjust his portfolio.
285 dvdmFunc_Adj : MargValueFunc2D
286 Marginal value over normalized liquid resources when the agent is able
287 to adjust his portfolio.
288 dvdnFunc_Adj : MargValueFunc2D
289 Marginal value over normalized liquid resources when the agent is able
290 to adjust his portfolio.
291 vFunc_Fxd : ValueFunc3D
292 Stage value function over normalized liquid resources, normalized
293 iliquid resources, and income contribution share when the agent is
294 not able to adjust his portfolio.
295 dfracFunc_Fxd : Interp2D
296 Deposit function over normalized liquid resources, normalized iliquid
297 resources, and income contribution share when the agent is not able to
298 adjust his portfolio.
299 Must be ConstantFunction(0.0)
300 dvdmFunc_Fxd : MargValueFunc3D
301 Marginal value over normalized liquid resources when the agent is not
302 able to adjust his portfolio.
303 dvdnFunc_Fxd : MargValueFunc3D
304 Marginal value over normalized iliquid resources when the agent is not
305 able to adjust his portfolio.
306 dvdsFunc_Fxd : Interp3D
307 Marginal value function over income contribution share when the agent
308 is not able to ajust his portfolio.
309 """
311 distance_criteria = ["dvdmFunc_Adj", "dvdnFunc_Adj"]
313 def __init__(
314 self,
315 # Rebalancing stage, adjusting
316 vFunc_Adj=None,
317 dfracFunc_Adj=None,
318 dvdmFunc_Adj=None,
319 dvdnFunc_Adj=None,
320 # Rebalancing stage, fixed
321 vFunc_Fxd=None,
322 dfracFunc_Fxd=None,
323 dvdmFunc_Fxd=None,
324 dvdnFunc_Fxd=None,
325 dvdsFunc_Fxd=None,
326 ):
327 # Rebalancing stage
328 if vFunc_Adj is None:
329 vFunc_Adj = NullFunc()
330 if dfracFunc_Adj is None:
331 dfracFunc_Adj = NullFunc()
332 if dvdmFunc_Adj is None:
333 dvdmFunc_Adj = NullFunc()
334 if dvdnFunc_Adj is None:
335 dvdnFunc_Adj = NullFunc()
337 if vFunc_Fxd is None:
338 vFunc_Fxd = NullFunc()
339 if dfracFunc_Fxd is None:
340 dfracFunc_Fxd = NullFunc()
341 if dvdmFunc_Fxd is None:
342 dvdmFunc_Fxd = NullFunc()
343 if dvdnFunc_Fxd is None:
344 dvdnFunc_Fxd = NullFunc()
345 if dvdsFunc_Fxd is None:
346 dvdsFunc_Fxd = NullFunc()
348 # Components of the adjusting problem
349 self.vFunc_Adj = vFunc_Adj
350 self.dfracFunc_Adj = dfracFunc_Adj
351 self.dvdmFunc_Adj = dvdmFunc_Adj
352 self.dvdnFunc_Adj = dvdnFunc_Adj
354 # Components of the fixed problem
355 self.vFunc_Fxd = vFunc_Fxd
356 self.dfracFunc_Fxd = dfracFunc_Fxd
357 self.dvdmFunc_Fxd = dvdmFunc_Fxd
358 self.dvdnFunc_Fxd = dvdnFunc_Fxd
359 self.dvdsFunc_Fxd = dvdsFunc_Fxd
362# Class for the contribution share stage solution
363class RiskyContribShaSolution(MetricObject):
364 """
365 A class for representing the solution to the contribution-share stage of
366 the 'RiskyContrib' model.
368 Parameters
369 ----------
370 vFunc_Adj : ValueFunc2D
371 Stage value function over normalized liquid resources and normalized
372 iliquid resources when the agent is able to adjust his portfolio.
373 ShareFunc_Adj : Interp2D
374 Income contribution share function over normalized liquid resources
375 and normalized iliquid resources when the agent is able to adjust his
376 portfolio.
377 dvdmFunc_Adj : MargValueFunc2D
378 Marginal value function over normalized liquid resources when the agent
379 is able to adjust his portfolio.
380 dvdnFunc_Adj : MargValueFunc2D
381 Marginal value function over normalized iliquid resources when the
382 agent is able to adjust his portfolio.
383 vFunc_Fxd : ValueFunc3D
384 Stage value function over normalized liquid resources, normalized
385 iliquid resources, and income contribution share when the agent is not
386 able to adjust his portfolio.
387 ShareFunc_Fxd : Interp3D
388 Income contribution share function over normalized liquid resources,
389 iliquid resources, and income contribution share when the agent is not
390 able to adjust his portfolio.
391 Should be an IdentityFunc.
392 dvdmFunc_Fxd : MargValueFunc3D
393 Marginal value function over normalized liquid resources when the agent
394 is not able to adjust his portfolio.
395 dvdnFunc_Fxd : MargValueFunc3D
396 Marginal value function over normalized iliquid resources when the
397 agent is not able to adjust his portfolio.
398 dvdsFunc_Fxd : Interp3D
399 Marginal value function over income contribution share when the agent
400 is not able to adjust his portfolio
401 """
403 distance_criteria = ["dvdmFunc_Adj", "dvdnFunc_Adj"]
405 def __init__(
406 self,
407 # Contribution stage, adjust
408 vFunc_Adj=None,
409 ShareFunc_Adj=None,
410 dvdmFunc_Adj=None,
411 dvdnFunc_Adj=None,
412 # Contribution stage, fixed
413 vFunc_Fxd=None,
414 ShareFunc_Fxd=None,
415 dvdmFunc_Fxd=None,
416 dvdnFunc_Fxd=None,
417 dvdsFunc_Fxd=None,
418 ):
419 # Contribution stage, adjust
420 if vFunc_Adj is None:
421 vFunc_Adj = NullFunc()
422 if ShareFunc_Adj is None:
423 ShareFunc_Adj = NullFunc()
424 if dvdmFunc_Adj is None:
425 dvdmFunc_Adj = NullFunc()
426 if dvdnFunc_Adj is None:
427 dvdnFunc_Adj = NullFunc()
429 # Contribution stage, fixed
430 if vFunc_Fxd is None:
431 vFunc_Fxd = NullFunc()
432 if ShareFunc_Fxd is None:
433 ShareFunc_Fxd = NullFunc()
434 if dvdmFunc_Fxd is None:
435 dvdmFunc_Fxd = NullFunc()
436 if dvdnFunc_Fxd is None:
437 dvdnFunc_Fxd = NullFunc()
438 if dvdsFunc_Fxd is None:
439 dvdsFunc_Fxd = NullFunc()
441 # Set attributes of self
442 self.vFunc_Adj = vFunc_Adj
443 self.ShareFunc_Adj = ShareFunc_Adj
444 self.dvdmFunc_Adj = dvdmFunc_Adj
445 self.dvdnFunc_Adj = dvdnFunc_Adj
447 self.vFunc_Fxd = vFunc_Fxd
448 self.ShareFunc_Fxd = ShareFunc_Fxd
449 self.dvdmFunc_Fxd = dvdmFunc_Fxd
450 self.dvdnFunc_Fxd = dvdnFunc_Fxd
451 self.dvdsFunc_Fxd = dvdsFunc_Fxd
454# Class for the consumption stage solution
455class RiskyContribCnsSolution(MetricObject):
456 """
457 A class for representing the solution to the consumption stage of the
458 'RiskyContrib' model.
460 Parameters
461 ----------
462 vFunc : ValueFunc3D
463 Stage-value function over normalized liquid resources, normalized
464 iliquid resources, and income contribution share.
465 cFunc : Interp3D
466 Consumption function over normalized liquid resources, normalized
467 iliquid resources, and income contribution share.
468 dvdmFunc : MargValueFunc3D
469 Marginal value function over normalized liquid resources.
470 dvdnFunc : MargValueFunc3D
471 Marginal value function over normalized iliquid resources.
472 dvdsFunc : Interp3D
473 Marginal value function over income contribution share.
474 """
476 distance_criteria = ["dvdmFunc", "dvdnFunc"]
478 def __init__(
479 self,
480 # Consumption stage
481 vFunc=None,
482 cFunc=None,
483 dvdmFunc=None,
484 dvdnFunc=None,
485 dvdsFunc=None,
486 ):
487 if vFunc is None:
488 vFunc = NullFunc()
489 if cFunc is None:
490 cFunc = NullFunc()
491 if dvdmFunc is None:
492 dvdmFunc = NullFunc()
493 if dvdnFunc is None:
494 dvdmFunc = NullFunc()
495 if dvdsFunc is None:
496 dvdsFunc = NullFunc()
498 self.vFunc = vFunc
499 self.cFunc = cFunc
500 self.dvdmFunc = dvdmFunc
501 self.dvdnFunc = dvdnFunc
502 self.dvdsFunc = dvdsFunc
505# Class for the solution of a whole period
506class RiskyContribSolution(MetricObject):
507 """
508 A class for representing the solution to a full time-period of the
509 'RiskyContrib' agent type's problem.
511 Parameters
512 ----------
513 Reb : RiskyContribRebSolution
514 Solution to the period's rebalancing stage.
515 Sha : RiskyContribShaSolution
516 Solution to the period's contribution-share stage.
517 Cns : RiskyContribCnsSolution
518 Solution to the period's consumption stage.
519 """
521 # Solutions are to be compared on the basis of their sub-period solutions
522 distance_criteria = ["stage_sols"]
524 def __init__(self, Reb, Sha, Cns):
525 # Dictionary of stage solutions
526 self.stage_sols = {"Reb": Reb, "Sha": Sha, "Cns": Cns}
529# %% Auxiliary functions and transition equations for the RiskyContrib model.
532def rebalance_assets(d, m, n, tau):
533 """
534 A function that produces post-rebalancing assets for given initial assets,
535 rebalancing action, and tax rate.
537 Parameters
538 ----------
539 d : np.array
540 Array with rebalancing decisions. d > 0 represents depositing d*m into
541 the risky asset account. d<0 represents withdrawing ``|d|*n`` (pre-tax)
542 from the risky account into the risky account.
543 m : np.array
544 Initial risk-free assets.
545 n : np.array
546 Initial risky assets.
547 tau : float
548 Tax rate on flows from the risky to the risk-free asset.
550 Returns
551 -------
552 mTil : np.array
553 Post-rebalancing risk-free assets.
554 nTil : np.arrat
555 Post-rebalancing risky assets.
557 """
558 # Initialize
559 mTil = np.zeros_like(m) + np.nan
560 nTil = np.zeros_like(m) + np.nan
562 # Contributions
563 inds = d >= 0
564 mTil[inds] = m[inds] * (1 - d[inds])
565 nTil[inds] = n[inds] + m[inds] * d[inds]
567 # Withdrawals
568 inds = d < 0
569 mTil[inds] = m[inds] - d[inds] * n[inds] * (1 - tau)
570 nTil[inds] = n[inds] * (1 + d[inds])
572 return (mTil, nTil)
575# Transition equations for the consumption stage
576def m_nrm_next(shocks, aNrm, Share, Rfree, PermGroFac):
577 """
578 Given end-of-period balances and contribution share and the
579 start-of-next-period shocks, figure out next period's normalized riskless
580 assets
582 Parameters
583 ----------
584 shocks : np.array
585 Length-3 array with the stochastic shocks that get realized between the
586 end of the current period and the start of next period. Their order is
587 (0) permanent income shock, (1) transitory income shock, (2) risky
588 asset return.
589 aNrm : float
590 End-of-period risk-free asset balances.
591 Share : float
592 End-of-period income deduction share.
593 Rfree : float
594 Risk-free return factor.
595 PermGroFac : float
596 Permanent income growth factor.
598 Returns
599 -------
600 m_nrm_tp1 : float
601 Next-period normalized riskless balance.
603 """
604 # Extract shocks
605 perm_shk = shocks["PermShk"]
606 tran_shk = shocks["TranShk"]
608 m_nrm_tp1 = Rfree * aNrm / (perm_shk * PermGroFac) + (1.0 - Share) * tran_shk
610 return m_nrm_tp1
613def n_nrm_next(shocks, nNrm, Share, PermGroFac):
614 """
615 Given end-of-period balances and contribution share and the
616 start-of-next-period shocks, figure out next period's normalized risky
617 assets
619 Parameters
620 ----------
621 shocks : dict
622 Dictionary with the stochastic shocks that get realized between the
623 end of the current period and the start of next period: "PermShk" is the
624 permanent income shock, "TranShk" is the transitory income shock, and
625 "Risky" is the risky asset return.
626 nNrm : float
627 End-of-period risky asset balances.
628 Share : float
629 End-of-period income deduction share.
630 PermGroFac : float
631 Permanent income growth factor.
633 Returns
634 -------
635 n_nrm_tp1 : float
636 Next-period normalized risky balance.
638 """
640 # Extract shocks
641 perm_shk = shocks["PermShk"]
642 tran_shk = shocks["TranShk"]
643 R_risky = shocks["Risky"]
645 n_nrm_tp1 = R_risky * nNrm / (perm_shk * PermGroFac) + Share * tran_shk
647 return n_nrm_tp1
650# %% RiskyContrib solvers
653# Consumption stage solver
654def solve_RiskyContrib_Cns(
655 solution_next,
656 ShockDstn,
657 IncShkDstn,
658 RiskyDstn,
659 IndepDstnBool,
660 LivPrb,
661 DiscFac,
662 CRRA,
663 Rfree,
664 PermGroFac,
665 BoroCnstArt,
666 aXtraGrid,
667 nNrmGrid,
668 mNrmGrid,
669 ShareGrid,
670 vFuncBool,
671 AdjustPrb,
672 DiscreteShareBool,
673 joint_dist_solver,
674 **unused_params,
675):
676 """
677 Solves the consumption stage of the agent's problem
679 Parameters
680 ----------
681 solution_next : RiskyContribRebSolution
682 Solution to the first stage of the next period in the agent's problem.
683 ShockDstn : DiscreteDistribution
684 Joint distribution of next period's (0) permanent income shock, (1)
685 transitory income shock, and (2) risky asset return factor.
686 IncShkDstn : DiscreteDistribution
687 Joint distribution of next period's (0) permanent income shock and (1)
688 transitory income shock.
689 RiskyDstn : DiscreteDistribution
690 Distribution of next period's risky asset return factor.
691 IndepDstnBool : bool
692 Indicates whether the income and risky return distributions are
693 independent.
694 LivPrb : float
695 Probability of surviving until next period.
696 DiscFac : float
697 Time-preference discount factor.
698 CRRA : float
699 Coefficient of relative risk aversion.
700 Rfree : float
701 Risk-free return factor.
702 PermGroFac : float
703 Deterministic permanent income growth factor.
704 BoroCnstArt : float
705 Minimum allowed market resources (must be 0).
706 aXtraGrid : numpy array
707 Exogenous grid for end-of-period risk free resources.
708 nNrmGrid : numpy array
709 Exogenous grid for risky resources.
710 mNrmGrid : numpy array
711 Exogenous grid for risk-free resources.
712 ShareGrid : numpt array
713 Exogenous grid for the income contribution share.
714 vFuncBool : bool
715 Boolean that determines wether the value function's level needs to be
716 computed.
717 AdjustPrb : float
718 Probability thet the agent will be able to adjust his portfolio next
719 period.
720 DiscreteShareBool : bool
721 Boolean that determines whether only a discrete set of contribution
722 shares (ShareGrid) is allowed.
723 joint_dist_solver: bool
724 Should the general solver be used even if income and returns are
725 independent?
727 Returns
728 -------
729 solution : RiskyContribCnsSolution
730 Solution to the agent's consumption stage problem.
732 """
733 # Make sure the individual is liquidity constrained. Allowing a consumer to
734 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
735 if BoroCnstArt != 0.0:
736 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
738 # Make sure that if risky portfolio share is optimized only discretely, then
739 # the value function is also constructed (else this task would be impossible).
740 if DiscreteShareBool and (not vFuncBool):
741 raise ValueError(
742 "PortfolioConsumerType requires vFuncBool to be True when DiscreteShareBool is True!"
743 )
745 # Define temporary functions for utility and its derivative and inverse
746 u = lambda x: utility(x, CRRA)
747 uPinv = lambda x: utilityP_inv(x, CRRA)
748 uInv = lambda x: utility_inv(x, CRRA)
750 # Unpack next period's solution
751 vFunc_Reb_Adj_next = solution_next.vFunc_Adj
752 dvdmFunc_Reb_Adj_next = solution_next.dvdmFunc_Adj
753 dvdnFunc_Reb_Adj_next = solution_next.dvdnFunc_Adj
755 vFunc_Reb_Fxd_next = solution_next.vFunc_Fxd
756 dvdmFunc_Reb_Fxd_next = solution_next.dvdmFunc_Fxd
757 dvdnFunc_Reb_Fxd_next = solution_next.dvdnFunc_Fxd
758 dvdsFunc_Reb_Fxd_next = solution_next.dvdsFunc_Fxd
760 # STEP ONE
761 # Find end-of-period (continuation) value function and its derivatives.
763 # Start by constructing functions for next-period's pre-adjustment-shock
764 # expected value functions
765 if AdjustPrb < 1.0:
766 dvdm_next = lambda m, n, s: AdjustPrb * dvdmFunc_Reb_Adj_next(m, n) + (
767 1.0 - AdjustPrb
768 ) * dvdmFunc_Reb_Fxd_next(m, n, s)
769 dvdn_next = lambda m, n, s: AdjustPrb * dvdnFunc_Reb_Adj_next(m, n) + (
770 1.0 - AdjustPrb
771 ) * dvdnFunc_Reb_Fxd_next(m, n, s)
772 dvds_next = lambda m, n, s: (1.0 - AdjustPrb) * dvdsFunc_Reb_Fxd_next(m, n, s)
774 # Value function if needed
775 if vFuncBool:
776 v_next = lambda m, n, s: AdjustPrb * vFunc_Reb_Adj_next(m, n) + (
777 1.0 - AdjustPrb
778 ) * vFunc_Reb_Fxd_next(m, n, s)
780 else:
781 dvdm_next = lambda m, n, s: dvdmFunc_Reb_Adj_next(m, n)
782 dvdn_next = lambda m, n, s: dvdnFunc_Reb_Adj_next(m, n)
783 dvds_next = ConstantFunction(0.0)
785 if vFuncBool:
786 v_next = lambda m, n, s: vFunc_Reb_Adj_next(m, n)
788 if IndepDstnBool and not joint_dist_solver:
789 # If income and returns are independent we can use the law of iterated
790 # expectations to speed up the computation of end-of-period derivatives
792 # Define "post-return variables"
793 # b_aux = aNrm * R
794 # g_aux = nNrmTilde * Rtilde
795 # and create a function that interpolates end-of-period marginal values
796 # as functions of those and the contribution share
798 def post_return_derivs(inc_shocks, b_aux, g_aux, s):
799 perm_shk = inc_shocks["PermShk"]
800 tran_shk = inc_shocks["TranShk"]
802 temp_fac_A = utilityP(perm_shk * PermGroFac, CRRA)
803 temp_fac_B = (perm_shk * PermGroFac) ** (1.0 - CRRA)
805 # Find next-period asset balances
806 m_next = b_aux / (perm_shk * PermGroFac) + (1.0 - s) * tran_shk
807 n_next = g_aux / (perm_shk * PermGroFac) + s * tran_shk
809 # Interpolate next-period-value derivatives
810 dvdm_tp1 = dvdm_next(m_next, n_next, s)
811 dvdn_tp1 = dvdn_next(m_next, n_next, s)
812 if tran_shk == 0:
813 dvds_tp1 = dvds_next(m_next, n_next, s)
814 else:
815 dvds_tp1 = tran_shk * (dvdn_tp1 - dvdm_tp1) + dvds_next(
816 m_next, n_next, s
817 )
819 # Discount next-period-value derivatives to current period
821 # Liquid resources
822 pr_dvda = temp_fac_A * dvdm_tp1
823 # Iliquid resources
824 pr_dvdn = temp_fac_A * dvdn_tp1
825 # Contribution share
826 pr_dvds = temp_fac_B * dvds_tp1
828 # End of period value function, if needed
829 if vFuncBool:
830 pr_v = temp_fac_B * v_next(m_next, n_next, s)
831 return np.stack([pr_dvda, pr_dvdn, pr_dvds, pr_v])
832 else:
833 return np.stack([pr_dvda, pr_dvdn, pr_dvds])
835 # Define grids
836 b_aux_grid = np.concatenate([np.array([0.0]), Rfree * aXtraGrid])
837 g_aux_grid = np.concatenate(
838 [np.array([0.0]), max(RiskyDstn.atoms.flatten()) * nNrmGrid]
839 )
841 # Create tiled arrays with conforming dimensions.
842 b_aux_tiled, g_aux_tiled, Share_tiled = np.meshgrid(
843 b_aux_grid, g_aux_grid, ShareGrid, indexing="ij"
844 )
846 # Find end of period derivatives and value as expectations of (discounted)
847 # next period's derivatives and value.
848 pr_derivs = expected(
849 post_return_derivs,
850 IncShkDstn,
851 (b_aux_tiled, g_aux_tiled, Share_tiled),
852 vectorized=False,
853 )
855 # Unpack results and create interpolators
856 pr_dvdb_func = MargValueFuncCRRA(
857 TrilinearInterp(uPinv(pr_derivs[0]), b_aux_grid, g_aux_grid, ShareGrid),
858 CRRA,
859 )
860 pr_dvdg_func = MargValueFuncCRRA(
861 TrilinearInterp(uPinv(pr_derivs[1]), b_aux_grid, g_aux_grid, ShareGrid),
862 CRRA,
863 )
864 pr_dvds_func = TrilinearInterp(pr_derivs[2], b_aux_grid, g_aux_grid, ShareGrid)
866 if vFuncBool:
867 pr_vFunc = ValueFuncCRRA(
868 TrilinearInterp(uInv(pr_derivs[3]), b_aux_grid, g_aux_grid, ShareGrid),
869 CRRA,
870 )
872 # Now construct a function that produces end-of-period derivatives
873 # given the risky return draw
874 def end_of_period_derivs(risky_ret, a, nTil, s):
875 """
876 Computes the end-of-period derivatives (and optionally the value) of the
877 continuation value function, conditional on risky returns. This is so that the
878 expectations can be calculated by integrating over risky returns.
880 Parameters
881 ----------
882 risky_ret : float
883 Risky return factor
884 a : float
885 end-of-period risk-free assets.
886 nTil : float
887 end-of-period risky assets.
888 s : float
889 end-of-period income deduction share.
890 """
892 # Find next-period asset balances
893 b_aux = a * Rfree
894 g_aux = nTil * risky_ret
896 # Interpolate post-return derivatives
897 pr_dvdb = pr_dvdb_func(b_aux, g_aux, s)
898 pr_dvdg = pr_dvdg_func(b_aux, g_aux, s)
899 pr_dvds = pr_dvds_func(b_aux, g_aux, s)
901 # Discount
903 # Liquid resources
904 end_of_prd_dvda = DiscFac * Rfree * LivPrb * pr_dvdb
905 # Iliquid resources
906 end_of_prd_dvdn = DiscFac * risky_ret * LivPrb * pr_dvdg
907 # Contribution share
908 end_of_prd_dvds = DiscFac * LivPrb * pr_dvds
910 # End of period value function, i11f needed
911 if vFuncBool:
912 end_of_prd_v = DiscFac * LivPrb * pr_vFunc(b_aux, g_aux, s)
913 return np.stack(
914 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
915 )
916 else:
917 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
919 else:
920 # If income and returns are not independent, we just integrate over
921 # them jointly.
923 # Construct a function that evaluates and discounts them given a
924 # vector of return and income shocks and an end-of-period state
925 def end_of_period_derivs(shocks, a, nTil, s):
926 """
927 Computes the end-of-period derivatives (and optionally the value) of the
928 continuation value function, conditional on shocks. This is so that the
929 expectations can be calculated by integrating over shocks.
931 Parameters
932 ----------
933 shocks : dict
934 Dictionary with the stochastic shocks that get realized between the
935 end of the current period and the start of next period: "PermShk" is the
936 permanent income shock, "TranShk" is the transitory income shock, and
937 "Risky" is the risky asset return.
938 a : float
939 end-of-period risk-free assets.
940 nTil : float
941 end-of-period risky assets.
942 s : float
943 end-of-period income deduction share.
944 """
945 perm_shk = shocks["PermShk"]
946 tran_shk = shocks["TranShk"]
947 risky = shocks["Risky"]
949 temp_fac_A = utilityP(perm_shk * PermGroFac, CRRA)
950 temp_fac_B = (perm_shk * PermGroFac) ** (1.0 - CRRA)
952 # Find next-period asset balances
953 m_next = m_nrm_next(shocks, a, s, Rfree, PermGroFac)
954 n_next = n_nrm_next(shocks, nTil, s, PermGroFac)
956 # Interpolate next-period-value derivatives
957 dvdm_tp1 = dvdm_next(m_next, n_next, s)
958 dvdn_tp1 = dvdn_next(m_next, n_next, s)
959 if tran_shk == 0:
960 dvds_tp1 = dvds_next(m_next, n_next, s)
961 else:
962 dvds_tp1 = tran_shk * (dvdn_tp1 - dvdm_tp1) + dvds_next(
963 m_next, n_next, s
964 )
966 # Discount next-period-value derivatives to current period
968 # Liquid resources
969 end_of_prd_dvda = DiscFac * Rfree * LivPrb * temp_fac_A * dvdm_tp1
970 # Iliquid resources
971 end_of_prd_dvdn = DiscFac * risky * LivPrb * temp_fac_A * dvdn_tp1
972 # Contribution share
973 end_of_prd_dvds = DiscFac * LivPrb * temp_fac_B * dvds_tp1
975 # End of period value function, i11f needed
976 if vFuncBool:
977 end_of_prd_v = DiscFac * LivPrb * temp_fac_B * v_next(m_next, n_next, s)
978 return np.stack(
979 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
980 )
981 else:
982 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
984 # Now find the expected values on a (a, nTil, s) grid
986 # The "inversion" machinery can deal with assets of 0 even if there is a
987 # natural borrowing constraint, so include zeros.
988 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
989 aNrmGrid = np.concatenate([np.array([0.0]), aXtraGrid])
991 # Create tiled arrays with conforming dimensions.
992 aNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
993 aNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
994 )
996 # Find end of period derivatives and value as expectations of (discounted)
997 # next period's derivatives and value.
998 eop_derivs = expected(
999 end_of_period_derivs,
1000 RiskyDstn if IndepDstnBool and not joint_dist_solver else ShockDstn,
1001 (aNrm_tiled, nNrm_tiled, Share_tiled),
1002 vectorized=False,
1003 )
1005 # Unpack results
1006 eop_dvdaNvrs = uPinv(eop_derivs[0])
1007 eop_dvdnNvrs = uPinv(eop_derivs[1])
1008 eop_dvds = eop_derivs[2]
1009 if vFuncBool:
1010 eop_vNvrs = uInv(eop_derivs[3])
1012 # Construct an interpolator for eop_V. It will be used later.
1013 eop_vFunc = ValueFuncCRRA(
1014 TrilinearInterp(eop_vNvrs, aNrmGrid, nNrmGrid, ShareGrid), CRRA
1015 )
1017 # STEP TWO:
1018 # Solve the consumption problem and create interpolators for c, vCns,
1019 # and its derivatives.
1021 # Apply EGM over liquid resources at every (n,s) to find consumption.
1022 c_end = eop_dvdaNvrs
1023 mNrm_end = aNrm_tiled + c_end
1025 # Now construct interpolators for c and the derivatives of vCns.
1026 # The m grid is different for every (n,s). We interpolate the object of
1027 # interest on the regular m grid for every (n,s). At the end we will have
1028 # values of the functions of interest on a regular (m,n,s) grid. We use
1029 # trilinear interpolation on those points.
1031 # Expand the exogenous m grid to contain 0.
1032 mNrmGrid = np.insert(mNrmGrid, 0, 0)
1034 # Dimensions might have changed, so re-create tiled arrays
1035 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1036 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1037 )
1039 # Initialize arrays
1040 c_vals = np.zeros_like(mNrm_tiled)
1041 dvdnNvrs_vals = np.zeros_like(mNrm_tiled)
1042 dvds_vals = np.zeros_like(mNrm_tiled)
1044 nNrm_N = nNrmGrid.size
1045 Share_N = ShareGrid.size
1046 for nInd in range(nNrm_N):
1047 for sInd in range(Share_N):
1048 # Extract the endogenous m grid for particular (n,s).
1049 m_ns = mNrm_end[:, nInd, sInd]
1051 # Check if there is a natural constraint
1052 if m_ns[0] == 0.0:
1053 # There's no need to insert points since we have m==0.0
1055 # c
1056 c_vals[:, nInd, sInd] = LinearInterp(m_ns, c_end[:, nInd, sInd])(
1057 mNrmGrid
1058 )
1060 # dvdnNvrs
1061 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1062 m_ns, eop_dvdnNvrs[:, nInd, sInd]
1063 )(mNrmGrid)
1065 # dvds
1066 dvds_vals[:, nInd, sInd] = LinearInterp(m_ns, eop_dvds[:, nInd, sInd])(
1067 mNrmGrid
1068 )
1070 else:
1071 # We know that:
1072 # -The lowest gridpoints of both a and n are 0.
1073 # -Consumption at m < m0 is m.
1074 # -dvdn_Fxd at (m,n) for m < m0(n) is dvdn_Fxd(m0,n)
1075 # -Same is true for dvds_Fxd
1077 m_ns = np.concatenate([np.array([0]), m_ns])
1079 # c
1080 c_vals[:, nInd, sInd] = LinearInterp(
1081 m_ns, np.concatenate([np.array([0]), c_end[:, nInd, sInd]])
1082 )(mNrmGrid)
1084 # dvdnNvrs
1085 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1086 m_ns,
1087 np.concatenate(
1088 [
1089 np.array([eop_dvdnNvrs[0, nInd, sInd]]),
1090 eop_dvdnNvrs[:, nInd, sInd],
1091 ]
1092 ),
1093 )(mNrmGrid)
1095 # dvds
1096 dvds_vals[:, nInd, sInd] = LinearInterp(
1097 m_ns,
1098 np.concatenate(
1099 [
1100 np.array([eop_dvds[0, nInd, sInd]]),
1101 eop_dvds[:, nInd, sInd],
1102 ]
1103 ),
1104 )(mNrmGrid)
1106 # With the arrays filled, create 3D interpolators
1108 # Consumption interpolator
1109 cFunc = TrilinearInterp(c_vals, mNrmGrid, nNrmGrid, ShareGrid)
1110 # dvdmCns interpolator
1111 dvdmFunc_Cns = MargValueFuncCRRA(cFunc, CRRA)
1112 # dvdnCns interpolator
1113 dvdnNvrsFunc = TrilinearInterp(dvdnNvrs_vals, mNrmGrid, nNrmGrid, ShareGrid)
1114 dvdnFunc_Cns = MargValueFuncCRRA(dvdnNvrsFunc, CRRA)
1115 # dvdsCns interpolator
1116 dvdsFunc_Cns = TrilinearInterp(dvds_vals, mNrmGrid, nNrmGrid, ShareGrid)
1118 # Compute value function if needed
1119 if vFuncBool:
1120 # Consumption in the regular grid
1121 aNrm_reg = mNrm_tiled - c_vals
1122 vCns = u(c_vals) + eop_vFunc(aNrm_reg, nNrm_tiled, Share_tiled)
1123 vNvrsCns = uInv(vCns)
1124 vNvrsFunc_Cns = TrilinearInterp(vNvrsCns, mNrmGrid, nNrmGrid, ShareGrid)
1125 vFunc_Cns = ValueFuncCRRA(vNvrsFunc_Cns, CRRA)
1126 else:
1127 vFunc_Cns = NullFunc()
1129 # Assemble solution
1130 solution = RiskyContribCnsSolution(
1131 vFunc=vFunc_Cns,
1132 cFunc=cFunc,
1133 dvdmFunc=dvdmFunc_Cns,
1134 dvdnFunc=dvdnFunc_Cns,
1135 dvdsFunc=dvdsFunc_Cns,
1136 )
1138 return solution
1141# Solver for the contribution stage
1142def solve_RiskyContrib_Sha(
1143 solution_next,
1144 CRRA,
1145 AdjustPrb,
1146 mNrmGrid,
1147 nNrmGrid,
1148 ShareGrid,
1149 DiscreteShareBool,
1150 vFuncBool,
1151 **unused_params,
1152):
1153 """
1154 Solves the income-contribution-share stag of the agent's problem
1156 Parameters
1157 ----------
1158 solution_next : RiskyContribCnsSolution
1159 Solution to the agent's consumption stage problem that follows.
1160 CRRA : float
1161 Coefficient of relative risk aversion.
1162 AdjustPrb : float
1163 Probability that the agent will be able to rebalance his portfolio
1164 next period.
1165 mNrmGrid : numpy array
1166 Exogenous grid for risk-free resources.
1167 nNrmGrid : numpy array
1168 Exogenous grid for risky resources.
1169 ShareGrid : numpy array
1170 Exogenous grid for the income contribution share.
1171 DiscreteShareBool : bool
1172 Boolean that determines whether only a discrete set of contribution
1173 shares (ShareGrid) is allowed.
1174 vFuncBool : bool
1175 Determines whether the level of the value function is computed.
1177 Yields
1178 ------
1179 solution : RiskyContribShaSolution
1180 Solution to the income-contribution-share stage of the agent's problem.
1182 """
1183 # Unpack solution from the next sub-stage
1184 vFunc_Cns_next = solution_next.vFunc
1185 cFunc_next = solution_next.cFunc
1186 dvdmFunc_Cns_next = solution_next.dvdmFunc
1187 dvdnFunc_Cns_next = solution_next.dvdnFunc
1188 dvdsFunc_Cns_next = solution_next.dvdsFunc
1190 uPinv = lambda x: utilityP_inv(x, CRRA)
1192 # Create tiled grids
1194 # Add 0 to the m and n grids
1195 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1196 nNrm_N = len(nNrmGrid)
1197 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1198 mNrm_N = len(mNrmGrid)
1200 if AdjustPrb == 1.0:
1201 # If the readjustment probability is 1, set the share to 0:
1202 # - If there is a withdrawal tax: better for the agent to observe
1203 # income before rebalancing.
1204 # - If there is no tax: all shares should yield the same value.
1205 mNrm_tiled, nNrm_tiled = np.meshgrid(mNrmGrid, nNrmGrid, indexing="ij")
1207 opt_idx = np.zeros_like(mNrm_tiled, dtype=int)
1208 opt_Share = ShareGrid[opt_idx]
1210 if vFuncBool:
1211 vNvrsSha = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, opt_Share)
1213 else:
1214 # Figure out optimal share by evaluating all alternatives at all
1215 # (m,n) combinations
1216 m_idx_tiled, n_idx_tiled = np.meshgrid(
1217 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1218 )
1220 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1221 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1222 )
1224 if DiscreteShareBool:
1225 # Evaluate value function to optimize over shares.
1226 # Do it in inverse space
1227 vNvrs = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, Share_tiled)
1229 # Find the optimal share at each (m,n).
1230 opt_idx = np.argmax(vNvrs, axis=2)
1232 # Compute objects needed for the value function and its derivatives
1233 vNvrsSha = vNvrs[m_idx_tiled, n_idx_tiled, opt_idx]
1234 opt_Share = ShareGrid[opt_idx]
1236 # Project grids
1237 mNrm_tiled = mNrm_tiled[:, :, 0]
1238 nNrm_tiled = nNrm_tiled[:, :, 0]
1240 else:
1241 # Evaluate the marginal value of the contribution share at
1242 # every (m,n,s) gridpoint
1243 dvds = dvdsFunc_Cns_next(mNrm_tiled, nNrm_tiled, Share_tiled)
1245 # If the derivative is negative at the lowest share, then s[0] is optimal
1246 constrained_bot = dvds[:, :, 0] <= 0.0
1247 # If it is poitive at the highest share, then s[-1] is optimal
1248 constrained_top = dvds[:, :, -1] >= 0.0
1250 # Find indices at which the derivative crosses 0 for the 1st time
1251 # will be 0 if it never does, but "constrained_top/bot" deals with that
1252 crossings = np.logical_and(dvds[:, :, :-1] >= 0.0, dvds[:, :, 1:] <= 0.0)
1253 idx = np.argmax(crossings, axis=2)
1255 # Linearly interpolate the optimal share
1256 idx1 = idx + 1
1257 slopes = (
1258 dvds[m_idx_tiled, n_idx_tiled, idx1]
1259 - dvds[m_idx_tiled, n_idx_tiled, idx]
1260 ) / (ShareGrid[idx1] - ShareGrid[idx])
1261 opt_Share = ShareGrid[idx] - dvds[m_idx_tiled, n_idx_tiled, idx] / slopes
1263 # Replace the ones we knew were constrained
1264 opt_Share[constrained_bot] = ShareGrid[0]
1265 opt_Share[constrained_top] = ShareGrid[-1]
1267 # Project grids
1268 mNrm_tiled = mNrm_tiled[:, :, 0]
1269 nNrm_tiled = nNrm_tiled[:, :, 0]
1271 # Evaluate the inverse value function at the optimal shares
1272 if vFuncBool:
1273 vNvrsSha = vFunc_Cns_next.func(mNrm_tiled, nNrm_tiled, opt_Share)
1275 dvdmNvrsSha = cFunc_next(mNrm_tiled, nNrm_tiled, opt_Share)
1276 dvdnSha = dvdnFunc_Cns_next(mNrm_tiled, nNrm_tiled, opt_Share)
1277 dvdnNvrsSha = uPinv(dvdnSha)
1279 # Interpolators
1281 # Value function if needed
1282 if vFuncBool:
1283 vNvrsFunc_Sha = BilinearInterp(vNvrsSha, mNrmGrid, nNrmGrid)
1284 vFunc_Sha = ValueFuncCRRA(vNvrsFunc_Sha, CRRA)
1285 else:
1286 vFunc_Sha = NullFunc()
1288 # Contribution share function
1289 if DiscreteShareBool:
1290 ShareFunc = DiscreteInterp(
1291 BilinearInterp(opt_idx, mNrmGrid, nNrmGrid), ShareGrid
1292 )
1293 else:
1294 ShareFunc = BilinearInterp(opt_Share, mNrmGrid, nNrmGrid)
1296 # Derivatives
1297 dvdmNvrsFunc_Sha = BilinearInterp(dvdmNvrsSha, mNrmGrid, nNrmGrid)
1298 dvdmFunc_Sha = MargValueFuncCRRA(dvdmNvrsFunc_Sha, CRRA)
1299 dvdnNvrsFunc_Sha = BilinearInterp(dvdnNvrsSha, mNrmGrid, nNrmGrid)
1300 dvdnFunc_Sha = MargValueFuncCRRA(dvdnNvrsFunc_Sha, CRRA)
1302 solution = RiskyContribShaSolution(
1303 vFunc_Adj=vFunc_Sha,
1304 ShareFunc_Adj=ShareFunc,
1305 dvdmFunc_Adj=dvdmFunc_Sha,
1306 dvdnFunc_Adj=dvdnFunc_Sha,
1307 # The fixed agent does nothing at this stage,
1308 # so his value functions are the next problem's
1309 vFunc_Fxd=vFunc_Cns_next,
1310 ShareFunc_Fxd=IdentityFunction(i_dim=2, n_dims=3),
1311 dvdmFunc_Fxd=dvdmFunc_Cns_next,
1312 dvdnFunc_Fxd=dvdnFunc_Cns_next,
1313 dvdsFunc_Fxd=dvdsFunc_Cns_next,
1314 )
1316 return solution
1319# Solver for the asset rebalancing stage
1320def solve_RiskyContrib_Reb(
1321 solution_next,
1322 CRRA,
1323 WithdrawTax,
1324 nNrmGrid,
1325 mNrmGrid,
1326 dfracGrid,
1327 vFuncBool,
1328 **unused_params,
1329):
1330 """
1331 Solves the asset-rebalancing-stage of the agent's problem
1333 Parameters
1334 ----------
1335 solution_next : RiskyContribShaSolution
1336 Solution to the income-contribution-share stage problem that follows.
1337 CRRA : float
1338 Coefficient of relative risk aversion.
1339 WithdrawTax : float
1340 Tax rate on risky asset withdrawals.
1341 nNrmGrid : numpy array
1342 Exogenous grid for risky resources.
1343 mNrmGrid : numpy array
1344 Exogenous grid for risk-free resources.
1345 dfracGrid : numpy array
1346 Grid for rebalancing flows. The final grid will be equivalent to
1347 [-nNrm*dfracGrid, dfracGrid*mNrm].
1348 vFuncBool : bool
1349 Determines whether the level of th value function must be computed.
1351 Returns
1352 -------
1353 solution : RiskyContribShaSolution
1354 Solution to the asset-rebalancing stage of the agent's problem.
1356 """
1357 # Extract next stage's solution
1358 vFunc_Adj_next = solution_next.vFunc_Adj
1359 dvdmFunc_Adj_next = solution_next.dvdmFunc_Adj
1360 dvdnFunc_Adj_next = solution_next.dvdnFunc_Adj
1362 vFunc_Fxd_next = solution_next.vFunc_Fxd
1363 dvdmFunc_Fxd_next = solution_next.dvdmFunc_Fxd
1364 dvdnFunc_Fxd_next = solution_next.dvdnFunc_Fxd
1365 dvdsFunc_Fxd_next = solution_next.dvdsFunc_Fxd
1367 uPinv = lambda x: utilityP_inv(x, CRRA)
1369 # Create tiled grids
1371 # Add 0 to the m and n grids
1372 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1373 nNrm_N = len(nNrmGrid)
1374 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1375 mNrm_N = len(mNrmGrid)
1376 d_N = len(dfracGrid)
1378 # Duplicate d so that possible values are -dfracGrid,dfracGrid. Duplicate 0 is
1379 # intentional since the tax causes a discontinuity. We need the value
1380 # from the left and right.
1381 dfracGrid = np.concatenate((-1 * np.flip(dfracGrid), dfracGrid))
1383 # It will be useful to pre-evaluate marginals at every (m,n,d) combination
1385 # Create tiled arrays for every d,m,n option
1386 d_N2 = len(dfracGrid)
1387 d_tiled, mNrm_tiled, nNrm_tiled = np.meshgrid(
1388 dfracGrid, mNrmGrid, nNrmGrid, indexing="ij"
1389 )
1391 # Get post-rebalancing assets.
1392 m_tilde, n_tilde = rebalance_assets(d_tiled, mNrm_tiled, nNrm_tiled, WithdrawTax)
1394 # Now the marginals, in inverse space
1395 dvdmNvrs = dvdmFunc_Adj_next.cFunc(m_tilde, n_tilde)
1396 dvdnNvrs = dvdnFunc_Adj_next.cFunc(m_tilde, n_tilde)
1398 # Pre-evaluate the inverse of (1-WithdrawTax)
1399 taxNvrs = uPinv(1 - WithdrawTax)
1400 # Create a tiled array of the tax
1401 taxNvrs_tiled = np.tile(
1402 np.reshape(
1403 np.concatenate([np.repeat(taxNvrs, d_N), np.ones(d_N, dtype=np.double)]),
1404 (d_N2, 1, 1),
1405 ),
1406 (1, mNrm_N, nNrm_N),
1407 )
1409 # The FOC is dvdn = tax*dvdm or dvdnNvrs = taxNvrs*dvdmNvrs
1410 dvdDNvrs = dvdnNvrs - taxNvrs_tiled * dvdmNvrs
1411 # The optimal d will be at the first point where dvdD < 0. The inverse
1412 # transformation flips the sign.
1414 # If the derivative is negative (inverse positive) at the lowest d,
1415 # then d == -1.0 is optimal
1416 constrained_bot = dvdDNvrs[0, :, :] >= 0.0
1417 # If it is positive (inverse negative) at the highest d, then d[-1] = 1.0
1418 # is optimal
1419 constrained_top = (
1420 dvdDNvrs[
1421 -1,
1422 :,
1423 :,
1424 ]
1425 <= 0.0
1426 )
1428 # Find indices at which the derivative crosses 0 for the 1st time
1429 # will be 0 if it never does, but "constrained_top/bot" deals with that
1430 crossings = np.logical_and(dvdDNvrs[:-1, :, :] <= 0.0, dvdDNvrs[1:, :, :] >= 0.0)
1431 idx = np.argmax(crossings, axis=0)
1433 m_idx_tiled, n_idx_tiled = np.meshgrid(
1434 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1435 )
1437 # Linearly interpolate the optimal withdrawal percentage d
1438 idx1 = idx + 1
1439 slopes = (
1440 dvdDNvrs[idx1, m_idx_tiled, n_idx_tiled]
1441 - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled]
1442 ) / (dfracGrid[idx1] - dfracGrid[idx])
1443 dfrac_opt = dfracGrid[idx] - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled] / slopes
1445 # Replace the ones we knew were constrained
1446 dfrac_opt[constrained_bot] = dfracGrid[0]
1447 dfrac_opt[constrained_top] = dfracGrid[-1]
1449 # Find m_tilde and n_tilde
1450 mtil_opt, ntil_opt = rebalance_assets(
1451 dfrac_opt, mNrm_tiled[0], nNrm_tiled[0], WithdrawTax
1452 )
1454 # Now the derivatives. These are not straight forward because of corner
1455 # solutions with partial derivatives that change the limits. The idea then
1456 # is to evaluate the possible uses of the marginal unit of resources and
1457 # take the maximum.
1459 # An additional unit of m
1460 marg_m = dvdmFunc_Adj_next(mtil_opt, ntil_opt)
1461 # An additional unit of n kept in n
1462 marg_n = dvdnFunc_Adj_next(mtil_opt, ntil_opt)
1463 # An additional unit of n withdrawn to m
1464 marg_n_to_m = marg_m * (1 - WithdrawTax)
1466 # Marginal value is the maximum of the marginals in their possible uses
1467 dvdm_Adj = np.maximum(marg_m, marg_n)
1468 dvdmNvrs_Adj = uPinv(dvdm_Adj)
1469 dvdn_Adj = np.maximum(marg_n, marg_n_to_m)
1470 dvdnNvrs_Adj = uPinv(dvdn_Adj)
1472 # Interpolators
1474 # Value
1475 if vFuncBool:
1476 vNvrs_Adj = vFunc_Adj_next.vFuncNvrs(mtil_opt, ntil_opt)
1477 vNvrsFunc_Adj = BilinearInterp(vNvrs_Adj, mNrmGrid, nNrmGrid)
1478 vFunc_Adj = ValueFuncCRRA(vNvrsFunc_Adj, CRRA)
1479 else:
1480 vFunc_Adj = NullFunc()
1482 # Marginals
1483 dvdmFunc_Adj = MargValueFuncCRRA(
1484 BilinearInterp(dvdmNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1485 )
1486 dvdnFunc_Adj = MargValueFuncCRRA(
1487 BilinearInterp(dvdnNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1488 )
1490 # Decison
1491 dfracFunc_Adj = BilinearInterp(dfrac_opt, mNrmGrid, nNrmGrid)
1493 solution = RiskyContribRebSolution(
1494 # Rebalancing stage adjusting
1495 vFunc_Adj=vFunc_Adj,
1496 dfracFunc_Adj=dfracFunc_Adj,
1497 dvdmFunc_Adj=dvdmFunc_Adj,
1498 dvdnFunc_Adj=dvdnFunc_Adj,
1499 # Rebalancing stage fixed (nothing happens, so value functions are
1500 # the ones from the next stage)
1501 vFunc_Fxd=vFunc_Fxd_next,
1502 dfracFunc_Fxd=ConstantFunction(0.0),
1503 dvdmFunc_Fxd=dvdmFunc_Fxd_next,
1504 dvdnFunc_Fxd=dvdnFunc_Fxd_next,
1505 dvdsFunc_Fxd=dvdsFunc_Fxd_next,
1506 )
1508 return solution
1511def solveRiskyContrib(
1512 solution_next,
1513 ShockDstn,
1514 IncShkDstn,
1515 RiskyDstn,
1516 IndepDstnBool,
1517 LivPrb,
1518 DiscFac,
1519 CRRA,
1520 Rfree,
1521 PermGroFac,
1522 WithdrawTax,
1523 BoroCnstArt,
1524 aXtraGrid,
1525 nNrmGrid,
1526 mNrmGrid,
1527 ShareGrid,
1528 dfracGrid,
1529 vFuncBool,
1530 AdjustPrb,
1531 DiscreteShareBool,
1532 joint_dist_solver,
1533):
1534 """
1535 Solve a full period (with its three stages) of the agent's problem
1537 Parameters
1538 ----------
1539 solution_next : RiskyContribSolution
1540 Solution to next period's problem.
1541 ShockDstn : DiscreteDistribution
1542 Joint distribution of next period's (0) permanent income shock, (1)
1543 transitory income shock, and (2) risky asset return factor.
1544 IncShkDstn : DiscreteDistribution
1545 Joint distribution of next period's (0) permanent income shock and (1)
1546 transitory income shock.
1547 RiskyDstn : DiscreteDistribution
1548 Distribution of next period's risky asset return factor.
1549 IndepDstnBool : bool
1550 Indicates whether the income and risky return distributions are
1551 independent.
1552 LivPrb : float
1553 Probability of surviving until next period.
1554 DiscFac : float
1555 Time-preference discount factor.
1556 CRRA : float
1557 Coefficient of relative risk aversion.
1558 Rfree : float
1559 Risk-free return factor.
1560 PermGroFac : float
1561 Deterministic permanent income growth factor.
1562 WithdrawTax : float
1563 Tax rate on risky asset withdrawals.
1564 BoroCnstArt : float
1565 Minimum allowed market resources (must be 0).
1566 aXtraGrid : numpy array
1567 Exogenous grid for end-of-period risk free resources.
1568 nNrmGrid : numpy array
1569 Exogenous grid for risky resources.
1570 mNrmGrid : numpy array
1571 Exogenous grid for risk-free resources.
1572 ShareGrid : numpy array
1573 Exogenous grid for the income contribution share.
1574 dfracGrid : numpy array
1575 Grid for rebalancing flows. The final grid will be equivalent to
1576 [-nNrm*dfracGrid, dfracGrid*mNrm].
1577 vFuncBool : bool
1578 Determines whether the level of th value function must be computed.
1579 AdjustPrb : float
1580 Probability that the agent will be able to rebalance his portfolio
1581 next period.
1582 DiscreteShareBool : bool
1583 Boolean that determines whether only a discrete set of contribution
1584 shares (ShareGrid) is allowed.
1585 joint_dist_solver: bool
1586 Should the general solver be used even if income and returns are
1587 independent?
1589 Returns
1590 -------
1591 periodSol : RiskyContribSolution
1592 Solution to the agent's current-period problem.
1594 """
1595 # Pack parameters to be passed to stage-specific solvers
1596 kws = {
1597 "ShockDstn": ShockDstn,
1598 "IncShkDstn": IncShkDstn,
1599 "RiskyDstn": RiskyDstn,
1600 "IndepDstnBool": IndepDstnBool,
1601 "LivPrb": LivPrb,
1602 "DiscFac": DiscFac,
1603 "CRRA": CRRA,
1604 "Rfree": Rfree,
1605 "PermGroFac": PermGroFac,
1606 "WithdrawTax": WithdrawTax,
1607 "BoroCnstArt": BoroCnstArt,
1608 "aXtraGrid": aXtraGrid,
1609 "nNrmGrid": nNrmGrid,
1610 "mNrmGrid": mNrmGrid,
1611 "ShareGrid": ShareGrid,
1612 "dfracGrid": dfracGrid,
1613 "vFuncBool": vFuncBool,
1614 "AdjustPrb": AdjustPrb,
1615 "DiscreteShareBool": DiscreteShareBool,
1616 "joint_dist_solver": joint_dist_solver,
1617 }
1619 # Stages of the problem in chronological order
1620 Stages = ["Reb", "Sha", "Cns"]
1621 n_stages = len(Stages)
1622 # Solvers, indexed by stage names
1623 Solvers = {
1624 "Reb": solve_RiskyContrib_Reb,
1625 "Sha": solve_RiskyContrib_Sha,
1626 "Cns": solve_RiskyContrib_Cns,
1627 }
1629 # Initialize empty solution
1630 stage_sols = {}
1631 # Solve stages backwards
1632 for i in reversed(range(n_stages)):
1633 stage = Stages[i]
1635 # In the last stage, the next solution is the first stage of the next
1636 # period. Otherwise, its the next stage of his period.
1637 if i == n_stages - 1:
1638 sol_next_stage = solution_next.stage_sols[Stages[0]]
1639 else:
1640 sol_next_stage = stage_sols[Stages[i + 1]]
1642 # Solve
1643 stage_sols[stage] = Solvers[stage](sol_next_stage, **kws)
1645 # Assemble stage solutions into period solution
1646 periodSol = RiskyContribSolution(**stage_sols)
1648 return periodSol
1651# %% Base risky-contrib dictionaries
1653risky_contrib_constructor_dict = {
1654 "IncShkDstn": construct_lognormal_income_process_unemployment,
1655 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1656 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1657 "aXtraGrid": make_assets_grid,
1658 "RiskyDstn": make_lognormal_RiskyDstn,
1659 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
1660 "ShareLimit": calc_ShareLimit_for_CRRA,
1661 "AdjustDstn": make_AdjustDstn,
1662 "solution_terminal": make_solution_terminal_risky_contrib,
1663 "ShareGrid": make_bounded_ShareGrid,
1664 "dfracGrid": make_simple_dGrid,
1665 "mNrmGrid": make_mNrm_grid,
1666 "nNrmGrid": make_nNrm_grid,
1667 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1668 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1669}
1671risky_contrib_params = {
1672 "constructors": risky_contrib_constructor_dict,
1673 # Preferences. The points of the model are more evident for more risk
1674 # averse and impatient agents
1675 "CRRA": 5.0,
1676 "DiscFac": 0.90,
1677 "WithdrawTax": [0.1],
1678 # Artificial borrowing constraint must be on
1679 "BoroCnstArt": 0.0,
1680 # Grids go up high wealth/P ratios and are less clustered at the bottom.
1681 "aXtraMax": 250,
1682 "aXtraCount": 50,
1683 "aXtraNestFac": 1,
1684 # Same goes for the new grids of the model
1685 "mNrmMin": 1e-6,
1686 "mNrmMax": 250,
1687 "mNrmCount": 50,
1688 "mNrmNestFac": 1,
1689 "nNrmMin": 1e-6,
1690 "nNrmMax": 250,
1691 "nNrmCount": 50,
1692 "nNrmNestFac": 1,
1693 # Income deduction/contribution share grid
1694 "ShareCount": 10,
1695 "ShareMax": 0.9,
1696 "DiscreteShareBool": False,
1697 # Grid for finding the optimal rebalancing flow
1698 "dCount": 20,
1699 "joint_dist_solver": False,
1700}
1701risky_asset_params = {
1702 # Risky return factor moments. Based on SP500 real returns from Shiller's
1703 # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
1704 "RiskyAvg": 1.080370891,
1705 "RiskyStd": 0.177196585,
1706 "ShareCount": 25, # Number of discrete points in the risky share approximation
1707 # Number of integration nodes to use in approximation of risky returns
1708 "RiskyCount": 5,
1709 # Probability that the agent can adjust their portfolio each period
1710 "AdjustPrb": [1.0],
1711 # When simulating the model, should all agents get the same risky return in
1712 # a given period?
1713 "sim_common_Rrisky": True,
1714}
1716# Infinite horizon version
1717init_risky_contrib = init_risky_asset.copy()
1718init_risky_contrib.update(risky_contrib_params)
1720# Lifecycle version
1721init_risky_contrib_lifecycle = init_lifecycle.copy()
1722init_risky_contrib_lifecycle.update(risky_asset_params)
1723init_risky_contrib_lifecycle.update(risky_contrib_params)
1725###############################################################################
1728class RiskyContribConsumerType(RiskyAssetConsumerType):
1729 """
1730 A consumer type with idiosyncratic shocks to permanent and transitory income,
1731 who can save in both a risk-free and a risky asset but faces frictions to
1732 moving funds between them. The agent can only consume out of his risk-free
1733 asset.
1735 The frictions are:
1737 - A proportional tax on funds moved from the risky to the risk-free
1738 asset.
1739 - A stochastic inability to move funds between his accounts.
1741 To partially avoid the second friction, the agent can commit to have a
1742 fraction of his labor income, which is usually deposited in his risk-free
1743 account, diverted to his risky account. He can change this fraction
1744 only in periods where he is able to move funds between accounts.
1745 """
1747 # The model is solved and simulated spliting each of the agent's
1748 # decisions into its own "stage". The stages in chronological order
1749 # are
1750 # - Reb: asset-rebalancing stage.
1751 # - Sha: definition of the income contribution share.
1752 # - Cns: consumption stage.
1753 stages = ["Reb", "Sha", "Cns"]
1754 # Each stage has its own states and controls, and its methods to find them.
1756 time_inv_ = RiskyAssetConsumerType.time_inv_ + [
1757 "DiscreteShareBool",
1758 "joint_dist_solver",
1759 "ShareGrid",
1760 "nNrmGrid",
1761 "mNrmGrid",
1762 "RiskyDstn",
1763 "dfracGrid",
1764 ]
1765 time_vary_ = RiskyAssetConsumerType.time_vary_ + ["WithdrawTax", "AdjustPrb"]
1767 # The new state variables (over those in ConsIndShock) are:
1768 # - nNrm: start-of-period risky resources.
1769 # - mNrmTilde: post-rebalancing risk-free resources.
1770 # - nNrmTilde: post-rebalancing risky resources.
1771 # - Share: income-deduction share.
1772 # For details, see
1773 # https://github.com/Mv77/RiskyContrib/blob/main/RiskyContrib.pdf
1774 state_vars = RiskyAssetConsumerType.state_vars + [
1775 "gNrm",
1776 "nNrm",
1777 "mNrmTilde",
1778 "nNrmTilde",
1779 "Share",
1780 ]
1781 shock_vars_ = RiskyAssetConsumerType.shock_vars_
1782 default_ = {
1783 "params": init_risky_contrib,
1784 "solver": solveRiskyContrib,
1785 "track_vars": ["aNrm", "cNrm", "mNrm", "nNrm", "dfrac", "Share", "pLvl"],
1786 }
1788 def __init__(self, **kwds):
1789 super().__init__(**kwds)
1790 # It looks like I can't assign this at the class level, unfortunately
1791 self.get_states = {
1792 "Reb": self.get_states_Reb,
1793 "Sha": self.get_states_Sha,
1794 "Cns": self.get_states_Cns,
1795 }
1796 self.get_controls = {
1797 "Reb": self.get_controls_Reb,
1798 "Sha": self.get_controls_Sha,
1799 "Cns": self.get_controls_Cns,
1800 }
1802 def pre_solve(self):
1803 self.construct("solution_terminal")
1805 def initialize_sim(self):
1806 """
1807 Initialize the state of simulation attributes.
1809 Parameters
1810 ----------
1811 None
1813 Returns
1814 -------
1815 None
1816 """
1817 RiskyAssetConsumerType.initialize_sim(self)
1818 self.state_now["Share"] = np.zeros(self.AgentCount)
1820 def sim_birth(self, which_agents):
1821 """
1822 Create new agents to replace ones who have recently died; takes draws of
1823 initial aNrm and pLvl, as in ConsIndShockModel, then sets Share, Adjust
1824 and post-rebalancing risky asset nNrmTilde to zero as initial values.
1825 Parameters
1826 ----------
1827 which_agents : np.array
1828 Boolean array of size AgentCount indicating which agents should be "born".
1830 Returns
1831 -------
1832 None
1833 """
1835 RiskyAssetConsumerType.sim_birth(self, which_agents)
1836 self.state_now["Share"][which_agents] = 0.0
1837 self.state_now["nNrmTilde"][which_agents] = 0.0
1839 def sim_one_period(self):
1840 """
1841 Simulates one period for this type.
1843 Has to be re-defined instead of using AgentType.sim_one_period() because
1844 of the "stages" structure.
1846 Parameters
1847 ----------
1848 None
1849 Returns
1850 -------
1851 None
1852 """
1854 if not hasattr(self, "solution"):
1855 raise Exception(
1856 "Model instance does not have a solution stored. To simulate, it is necessary"
1857 " to run the `solve()` method of the class first."
1858 )
1860 # Mortality adjusts the agent population
1861 self.get_mortality() # Replace some agents with "newborns"
1863 # Make state_now into state_prev, clearing state_now
1864 for var in self.state_now:
1865 self.state_prev[var] = self.state_now[var]
1867 if isinstance(self.state_now[var], np.ndarray):
1868 self.state_now[var] = np.empty(self.AgentCount)
1869 else:
1870 # Probably an aggregate variable. It may be getting set by the Market.
1871 pass
1873 if self.read_shocks: # If shock histories have been pre-specified, use those
1874 self.read_shocks_from_history()
1875 else: # Otherwise, draw shocks as usual according to subclass-specific method
1876 self.get_shocks()
1878 # Sequentially get states and controls of every stage
1879 for s in self.stages:
1880 self.get_states[s]()
1881 self.get_controls[s]()
1883 self.get_post_states()
1885 # Advance time for all agents
1886 self.t_age = self.t_age + 1 # Age all consumers by one period
1887 self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle
1888 self.t_cycle[self.t_cycle == self.T_cycle] = (
1889 0 # Resetting to zero for those who have reached the end
1890 )
1892 def get_states_Reb(self):
1893 """
1894 Get states for the first "stage": rebalancing.
1895 """
1897 pLvlPrev = self.state_prev["pLvl"]
1898 aNrmPrev = self.state_prev["aNrm"]
1899 SharePrev = self.state_prev["Share"]
1900 nNrmTildePrev = self.state_prev["nNrmTilde"]
1901 Rfree = np.array(self.Rfree)[self.t_cycle - 1]
1902 Rrisk = self.shocks["Risky"]
1904 # Calculate new states:
1906 # Permanent income
1907 self.state_now["pLvl"] = pLvlPrev * self.shocks["PermShk"]
1908 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1910 # Assets: mNrm and nNrm
1912 # Compute the effective growth factor of each asset
1913 RfEff = Rfree / self.shocks["PermShk"]
1914 RrEff = Rrisk / self.shocks["PermShk"]
1916 self.state_now["bNrm"] = RfEff * aNrmPrev # Liquid balances before labor income
1917 self.state_now["gNrm"] = (
1918 RrEff * nNrmTildePrev
1919 ) # Iliquid balances before labor income
1921 # Liquid balances after labor income
1922 self.state_now["mNrm"] = self.state_now["bNrm"] + self.shocks["TranShk"] * (
1923 1 - SharePrev
1924 )
1925 # Iliquid balances after labor income
1926 self.state_now["nNrm"] = (
1927 self.state_now["gNrm"] + self.shocks["TranShk"] * SharePrev
1928 )
1930 return None
1932 def get_controls_Reb(self):
1933 """
1934 Get controls for the first stage: rebalancing
1935 """
1936 dfrac = np.zeros(self.AgentCount) + np.nan
1938 # Loop over each period of the cycle, getting controls separately depending on "age"
1939 for t in range(self.T_cycle):
1940 # Find agents in this period-stage
1941 these = t == self.t_cycle
1943 # Get controls for agents who *can* adjust.
1944 those = np.logical_and(these, self.shocks["Adjust"])
1945 dfrac[those] = (
1946 self.solution[t]
1947 .stage_sols["Reb"]
1948 .dfracFunc_Adj(
1949 self.state_now["mNrm"][those], self.state_now["nNrm"][those]
1950 )
1951 )
1953 # Get Controls for agents who *can't* adjust.
1954 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
1955 dfrac[those] = (
1956 self.solution[t]
1957 .stage_sols["Reb"]
1958 .dfracFunc_Fxd(
1959 self.state_now["mNrm"][those],
1960 self.state_now["nNrm"][those],
1961 self.state_prev["Share"][those],
1962 )
1963 )
1965 # Limit dfrac to [-1,1] to prevent negative balances. Values outside
1966 # the range can come from extrapolation.
1967 self.controls["dfrac"] = np.minimum(np.maximum(dfrac, -1), 1.0)
1969 def get_states_Sha(self):
1970 """
1971 Get states for the second "stage": choosing the contribution share.
1972 """
1974 # Post-states are assets after rebalancing
1976 if "WithdrawTax" not in self.time_vary:
1977 mNrmTilde, nNrmTilde = rebalance_assets(
1978 self.controls["dfrac"],
1979 self.state_now["mNrm"],
1980 self.state_now["nNrm"],
1981 self.WithdrawTax,
1982 )
1984 else:
1985 # Initialize
1986 mNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1987 nNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1989 # Loop over each period of the cycle, getting controls separately depending on "age"
1990 for t in range(self.T_cycle):
1991 # Find agents in this period-stage
1992 these = t == self.t_cycle
1994 if np.sum(these) > 0:
1995 WithdrawTax = self.WithdrawTax[t]
1997 mNrmTilde[these], nNrmTilde[these] = rebalance_assets(
1998 self.controls["dfrac"][these],
1999 self.state_now["mNrm"][these],
2000 self.state_now["nNrm"][these],
2001 WithdrawTax,
2002 )
2004 self.state_now["mNrmTilde"] = mNrmTilde
2005 self.state_now["nNrmTilde"] = nNrmTilde
2007 def get_controls_Sha(self):
2008 """
2009 Get controls for the second "stage": choosing the contribution share.
2010 """
2012 Share = np.zeros(self.AgentCount) + np.nan
2014 # Loop over each period of the cycle, getting controls separately depending on "age"
2015 for t in range(self.T_cycle):
2016 # Find agents in this period-stage
2017 these = t == self.t_cycle
2019 # Get controls for agents who *can* adjust.
2020 those = np.logical_and(these, self.shocks["Adjust"])
2021 Share[those] = (
2022 self.solution[t]
2023 .stage_sols["Sha"]
2024 .ShareFunc_Adj(
2025 self.state_now["mNrmTilde"][those],
2026 self.state_now["nNrmTilde"][those],
2027 )
2028 )
2030 # Get Controls for agents who *can't* adjust.
2031 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
2032 Share[those] = (
2033 self.solution[t]
2034 .stage_sols["Sha"]
2035 .ShareFunc_Fxd(
2036 self.state_now["mNrmTilde"][those],
2037 self.state_now["nNrmTilde"][those],
2038 self.state_prev["Share"][those],
2039 )
2040 )
2042 # Store controls as attributes of self
2043 self.controls["Share"] = Share
2045 def get_states_Cns(self):
2046 """
2047 Get states for the third "stage": consumption.
2048 """
2050 # Contribution share becomes a state in the consumption problem
2051 self.state_now["Share"] = self.controls["Share"]
2053 def get_controls_Cns(self):
2054 """
2055 Get controls for the third "stage": consumption.
2056 """
2058 cNrm = np.zeros(self.AgentCount) + np.nan
2060 # Loop over each period of the cycle, getting controls separately depending on "age"
2061 for t in range(self.T_cycle):
2062 # Find agents in this period-stage
2063 these = t == self.t_cycle
2065 # Get consumption
2066 cNrm[these] = (
2067 self.solution[t]
2068 .stage_sols["Cns"]
2069 .cFunc(
2070 self.state_now["mNrmTilde"][these],
2071 self.state_now["nNrmTilde"][these],
2072 self.state_now["Share"][these],
2073 )
2074 )
2076 # Store controls as attributes of self
2077 # Since agents might be willing to end the period with a = 0, make
2078 # sure consumption does not go over m because of some numerical error.
2079 self.controls["cNrm"] = np.minimum(cNrm, self.state_now["mNrmTilde"])
2081 def get_post_states(self):
2082 """
2083 Set variables that are not a state to any problem but need to be
2084 computed in order to interact with shocks and produce next period's
2085 states.
2086 """
2087 self.state_now["aNrm"] = self.state_now["mNrmTilde"] - self.controls["cNrm"]