Coverage for HARK / ConsumptionSaving / ConsRiskyContribModel.py: 98%
508 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-08 05:31 +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 calc_expectation
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 = calc_expectation(
849 IncShkDstn, post_return_derivs, b_aux_tiled, g_aux_tiled, Share_tiled
850 )
852 # Unpack results and create interpolators
853 pr_dvdb_func = MargValueFuncCRRA(
854 TrilinearInterp(uPinv(pr_derivs[0]), b_aux_grid, g_aux_grid, ShareGrid),
855 CRRA,
856 )
857 pr_dvdg_func = MargValueFuncCRRA(
858 TrilinearInterp(uPinv(pr_derivs[1]), b_aux_grid, g_aux_grid, ShareGrid),
859 CRRA,
860 )
861 pr_dvds_func = TrilinearInterp(pr_derivs[2], b_aux_grid, g_aux_grid, ShareGrid)
863 if vFuncBool:
864 pr_vFunc = ValueFuncCRRA(
865 TrilinearInterp(uInv(pr_derivs[3]), b_aux_grid, g_aux_grid, ShareGrid),
866 CRRA,
867 )
869 # Now construct a function that produces end-of-period derivatives
870 # given the risky return draw
871 def end_of_period_derivs(risky_ret, a, nTil, s):
872 """
873 Computes the end-of-period derivatives (and optionally the value) of the
874 continuation value function, conditional on risky returns. This is so that the
875 expectations can be calculated by integrating over risky returns.
877 Parameters
878 ----------
879 risky_ret : float
880 Risky return factor
881 a : float
882 end-of-period risk-free assets.
883 nTil : float
884 end-of-period risky assets.
885 s : float
886 end-of-period income deduction share.
887 """
889 # Find next-period asset balances
890 b_aux = a * Rfree
891 g_aux = nTil * risky_ret
893 # Interpolate post-return derivatives
894 pr_dvdb = pr_dvdb_func(b_aux, g_aux, s)
895 pr_dvdg = pr_dvdg_func(b_aux, g_aux, s)
896 pr_dvds = pr_dvds_func(b_aux, g_aux, s)
898 # Discount
900 # Liquid resources
901 end_of_prd_dvda = DiscFac * Rfree * LivPrb * pr_dvdb
902 # Iliquid resources
903 end_of_prd_dvdn = DiscFac * risky_ret * LivPrb * pr_dvdg
904 # Contribution share
905 end_of_prd_dvds = DiscFac * LivPrb * pr_dvds
907 # End of period value function, i11f needed
908 if vFuncBool:
909 end_of_prd_v = DiscFac * LivPrb * pr_vFunc(b_aux, g_aux, s)
910 return np.stack(
911 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
912 )
913 else:
914 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
916 else:
917 # If income and returns are not independent, we just integrate over
918 # them jointly.
920 # Construct a function that evaluates and discounts them given a
921 # vector of return and income shocks and an end-of-period state
922 def end_of_period_derivs(shocks, a, nTil, s):
923 """
924 Computes the end-of-period derivatives (and optionally the value) of the
925 continuation value function, conditional on shocks. This is so that the
926 expectations can be calculated by integrating over shocks.
928 Parameters
929 ----------
930 shocks : dict
931 Dictionary with the stochastic shocks that get realized between the
932 end of the current period and the start of next period: "PermShk" is the
933 permanent income shock, "TranShk" is the transitory income shock, and
934 "Risky" is the risky asset return.
935 a : float
936 end-of-period risk-free assets.
937 nTil : float
938 end-of-period risky assets.
939 s : float
940 end-of-period income deduction share.
941 """
942 perm_shk = shocks["PermShk"]
943 tran_shk = shocks["TranShk"]
944 risky = shocks["Risky"]
946 temp_fac_A = utilityP(perm_shk * PermGroFac, CRRA)
947 temp_fac_B = (perm_shk * PermGroFac) ** (1.0 - CRRA)
949 # Find next-period asset balances
950 m_next = m_nrm_next(shocks, a, s, Rfree, PermGroFac)
951 n_next = n_nrm_next(shocks, nTil, s, PermGroFac)
953 # Interpolate next-period-value derivatives
954 dvdm_tp1 = dvdm_next(m_next, n_next, s)
955 dvdn_tp1 = dvdn_next(m_next, n_next, s)
956 if tran_shk == 0:
957 dvds_tp1 = dvds_next(m_next, n_next, s)
958 else:
959 dvds_tp1 = tran_shk * (dvdn_tp1 - dvdm_tp1) + dvds_next(
960 m_next, n_next, s
961 )
963 # Discount next-period-value derivatives to current period
965 # Liquid resources
966 end_of_prd_dvda = DiscFac * Rfree * LivPrb * temp_fac_A * dvdm_tp1
967 # Iliquid resources
968 end_of_prd_dvdn = DiscFac * risky * LivPrb * temp_fac_A * dvdn_tp1
969 # Contribution share
970 end_of_prd_dvds = DiscFac * LivPrb * temp_fac_B * dvds_tp1
972 # End of period value function, i11f needed
973 if vFuncBool:
974 end_of_prd_v = DiscFac * LivPrb * temp_fac_B * v_next(m_next, n_next, s)
975 return np.stack(
976 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
977 )
978 else:
979 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
981 # Now find the expected values on a (a, nTil, s) grid
983 # The "inversion" machinery can deal with assets of 0 even if there is a
984 # natural borrowing constraint, so include zeros.
985 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
986 aNrmGrid = np.concatenate([np.array([0.0]), aXtraGrid])
988 # Create tiled arrays with conforming dimensions.
989 aNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
990 aNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
991 )
993 # Find end of period derivatives and value as expectations of (discounted)
994 # next period's derivatives and value.
995 eop_derivs = calc_expectation(
996 RiskyDstn if IndepDstnBool and not joint_dist_solver else ShockDstn,
997 end_of_period_derivs,
998 aNrm_tiled,
999 nNrm_tiled,
1000 Share_tiled,
1001 )
1003 # Unpack results
1004 eop_dvdaNvrs = uPinv(eop_derivs[0])
1005 eop_dvdnNvrs = uPinv(eop_derivs[1])
1006 eop_dvds = eop_derivs[2]
1007 if vFuncBool:
1008 eop_vNvrs = uInv(eop_derivs[3])
1010 # Construct an interpolator for eop_V. It will be used later.
1011 eop_vFunc = ValueFuncCRRA(
1012 TrilinearInterp(eop_vNvrs, aNrmGrid, nNrmGrid, ShareGrid), CRRA
1013 )
1015 # STEP TWO:
1016 # Solve the consumption problem and create interpolators for c, vCns,
1017 # and its derivatives.
1019 # Apply EGM over liquid resources at every (n,s) to find consumption.
1020 c_end = eop_dvdaNvrs
1021 mNrm_end = aNrm_tiled + c_end
1023 # Now construct interpolators for c and the derivatives of vCns.
1024 # The m grid is different for every (n,s). We interpolate the object of
1025 # interest on the regular m grid for every (n,s). At the end we will have
1026 # values of the functions of interest on a regular (m,n,s) grid. We use
1027 # trilinear interpolation on those points.
1029 # Expand the exogenous m grid to contain 0.
1030 mNrmGrid = np.insert(mNrmGrid, 0, 0)
1032 # Dimensions might have changed, so re-create tiled arrays
1033 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1034 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1035 )
1037 # Initialize arrays
1038 c_vals = np.zeros_like(mNrm_tiled)
1039 dvdnNvrs_vals = np.zeros_like(mNrm_tiled)
1040 dvds_vals = np.zeros_like(mNrm_tiled)
1042 nNrm_N = nNrmGrid.size
1043 Share_N = ShareGrid.size
1044 for nInd in range(nNrm_N):
1045 for sInd in range(Share_N):
1046 # Extract the endogenous m grid for particular (n,s).
1047 m_ns = mNrm_end[:, nInd, sInd]
1049 # Check if there is a natural constraint
1050 if m_ns[0] == 0.0:
1051 # There's no need to insert points since we have m==0.0
1053 # c
1054 c_vals[:, nInd, sInd] = LinearInterp(m_ns, c_end[:, nInd, sInd])(
1055 mNrmGrid
1056 )
1058 # dvdnNvrs
1059 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1060 m_ns, eop_dvdnNvrs[:, nInd, sInd]
1061 )(mNrmGrid)
1063 # dvds
1064 dvds_vals[:, nInd, sInd] = LinearInterp(m_ns, eop_dvds[:, nInd, sInd])(
1065 mNrmGrid
1066 )
1068 else:
1069 # We know that:
1070 # -The lowest gridpoints of both a and n are 0.
1071 # -Consumption at m < m0 is m.
1072 # -dvdn_Fxd at (m,n) for m < m0(n) is dvdn_Fxd(m0,n)
1073 # -Same is true for dvds_Fxd
1075 m_ns = np.concatenate([np.array([0]), m_ns])
1077 # c
1078 c_vals[:, nInd, sInd] = LinearInterp(
1079 m_ns, np.concatenate([np.array([0]), c_end[:, nInd, sInd]])
1080 )(mNrmGrid)
1082 # dvdnNvrs
1083 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1084 m_ns,
1085 np.concatenate(
1086 [
1087 np.array([eop_dvdnNvrs[0, nInd, sInd]]),
1088 eop_dvdnNvrs[:, nInd, sInd],
1089 ]
1090 ),
1091 )(mNrmGrid)
1093 # dvds
1094 dvds_vals[:, nInd, sInd] = LinearInterp(
1095 m_ns,
1096 np.concatenate(
1097 [
1098 np.array([eop_dvds[0, nInd, sInd]]),
1099 eop_dvds[:, nInd, sInd],
1100 ]
1101 ),
1102 )(mNrmGrid)
1104 # With the arrays filled, create 3D interpolators
1106 # Consumption interpolator
1107 cFunc = TrilinearInterp(c_vals, mNrmGrid, nNrmGrid, ShareGrid)
1108 # dvdmCns interpolator
1109 dvdmFunc_Cns = MargValueFuncCRRA(cFunc, CRRA)
1110 # dvdnCns interpolator
1111 dvdnNvrsFunc = TrilinearInterp(dvdnNvrs_vals, mNrmGrid, nNrmGrid, ShareGrid)
1112 dvdnFunc_Cns = MargValueFuncCRRA(dvdnNvrsFunc, CRRA)
1113 # dvdsCns interpolator
1114 dvdsFunc_Cns = TrilinearInterp(dvds_vals, mNrmGrid, nNrmGrid, ShareGrid)
1116 # Compute value function if needed
1117 if vFuncBool:
1118 # Consumption in the regular grid
1119 aNrm_reg = mNrm_tiled - c_vals
1120 vCns = u(c_vals) + eop_vFunc(aNrm_reg, nNrm_tiled, Share_tiled)
1121 vNvrsCns = uInv(vCns)
1122 vNvrsFunc_Cns = TrilinearInterp(vNvrsCns, mNrmGrid, nNrmGrid, ShareGrid)
1123 vFunc_Cns = ValueFuncCRRA(vNvrsFunc_Cns, CRRA)
1124 else:
1125 vFunc_Cns = NullFunc()
1127 # Assemble solution
1128 solution = RiskyContribCnsSolution(
1129 vFunc=vFunc_Cns,
1130 cFunc=cFunc,
1131 dvdmFunc=dvdmFunc_Cns,
1132 dvdnFunc=dvdnFunc_Cns,
1133 dvdsFunc=dvdsFunc_Cns,
1134 )
1136 return solution
1139# Solver for the contribution stage
1140def solve_RiskyContrib_Sha(
1141 solution_next,
1142 CRRA,
1143 AdjustPrb,
1144 mNrmGrid,
1145 nNrmGrid,
1146 ShareGrid,
1147 DiscreteShareBool,
1148 vFuncBool,
1149 **unused_params,
1150):
1151 """
1152 Solves the income-contribution-share stag of the agent's problem
1154 Parameters
1155 ----------
1156 solution_next : RiskyContribCnsSolution
1157 Solution to the agent's consumption stage problem that follows.
1158 CRRA : float
1159 Coefficient of relative risk aversion.
1160 AdjustPrb : float
1161 Probability that the agent will be able to rebalance his portfolio
1162 next period.
1163 mNrmGrid : numpy array
1164 Exogenous grid for risk-free resources.
1165 nNrmGrid : numpy array
1166 Exogenous grid for risky resources.
1167 ShareGrid : numpy array
1168 Exogenous grid for the income contribution share.
1169 DiscreteShareBool : bool
1170 Boolean that determines whether only a discrete set of contribution
1171 shares (ShareGrid) is allowed.
1172 vFuncBool : bool
1173 Determines whether the level of the value function is computed.
1175 Yields
1176 ------
1177 solution : RiskyContribShaSolution
1178 Solution to the income-contribution-share stage of the agent's problem.
1180 """
1181 # Unpack solution from the next sub-stage
1182 vFunc_Cns_next = solution_next.vFunc
1183 cFunc_next = solution_next.cFunc
1184 dvdmFunc_Cns_next = solution_next.dvdmFunc
1185 dvdnFunc_Cns_next = solution_next.dvdnFunc
1186 dvdsFunc_Cns_next = solution_next.dvdsFunc
1188 uPinv = lambda x: utilityP_inv(x, CRRA)
1190 # Create tiled grids
1192 # Add 0 to the m and n grids
1193 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1194 nNrm_N = len(nNrmGrid)
1195 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1196 mNrm_N = len(mNrmGrid)
1198 if AdjustPrb == 1.0:
1199 # If the readjustment probability is 1, set the share to 0:
1200 # - If there is a withdrawal tax: better for the agent to observe
1201 # income before rebalancing.
1202 # - If there is no tax: all shares should yield the same value.
1203 mNrm_tiled, nNrm_tiled = np.meshgrid(mNrmGrid, nNrmGrid, indexing="ij")
1205 opt_idx = np.zeros_like(mNrm_tiled, dtype=int)
1206 opt_Share = ShareGrid[opt_idx]
1208 if vFuncBool:
1209 vNvrsSha = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, opt_Share)
1211 else:
1212 # Figure out optimal share by evaluating all alternatives at all
1213 # (m,n) combinations
1214 m_idx_tiled, n_idx_tiled = np.meshgrid(
1215 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1216 )
1218 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1219 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1220 )
1222 if DiscreteShareBool:
1223 # Evaluate value function to optimize over shares.
1224 # Do it in inverse space
1225 vNvrs = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, Share_tiled)
1227 # Find the optimal share at each (m,n).
1228 opt_idx = np.argmax(vNvrs, axis=2)
1230 # Compute objects needed for the value function and its derivatives
1231 vNvrsSha = vNvrs[m_idx_tiled, n_idx_tiled, opt_idx]
1232 opt_Share = ShareGrid[opt_idx]
1234 # Project grids
1235 mNrm_tiled = mNrm_tiled[:, :, 0]
1236 nNrm_tiled = nNrm_tiled[:, :, 0]
1238 else:
1239 # Evaluate the marginal value of the contribution share at
1240 # every (m,n,s) gridpoint
1241 dvds = dvdsFunc_Cns_next(mNrm_tiled, nNrm_tiled, Share_tiled)
1243 # If the derivative is negative at the lowest share, then s[0] is optimal
1244 constrained_bot = dvds[:, :, 0] <= 0.0
1245 # If it is poitive at the highest share, then s[-1] is optimal
1246 constrained_top = dvds[:, :, -1] >= 0.0
1248 # Find indices at which the derivative crosses 0 for the 1st time
1249 # will be 0 if it never does, but "constrained_top/bot" deals with that
1250 crossings = np.logical_and(dvds[:, :, :-1] >= 0.0, dvds[:, :, 1:] <= 0.0)
1251 idx = np.argmax(crossings, axis=2)
1253 # Linearly interpolate the optimal share
1254 idx1 = idx + 1
1255 slopes = (
1256 dvds[m_idx_tiled, n_idx_tiled, idx1]
1257 - dvds[m_idx_tiled, n_idx_tiled, idx]
1258 ) / (ShareGrid[idx1] - ShareGrid[idx])
1259 opt_Share = ShareGrid[idx] - dvds[m_idx_tiled, n_idx_tiled, idx] / slopes
1261 # Replace the ones we knew were constrained
1262 opt_Share[constrained_bot] = ShareGrid[0]
1263 opt_Share[constrained_top] = ShareGrid[-1]
1265 # Project grids
1266 mNrm_tiled = mNrm_tiled[:, :, 0]
1267 nNrm_tiled = nNrm_tiled[:, :, 0]
1269 # Evaluate the inverse value function at the optimal shares
1270 if vFuncBool:
1271 vNvrsSha = vFunc_Cns_next.func(mNrm_tiled, nNrm_tiled, opt_Share)
1273 dvdmNvrsSha = cFunc_next(mNrm_tiled, nNrm_tiled, opt_Share)
1274 dvdnSha = dvdnFunc_Cns_next(mNrm_tiled, nNrm_tiled, opt_Share)
1275 dvdnNvrsSha = uPinv(dvdnSha)
1277 # Interpolators
1279 # Value function if needed
1280 if vFuncBool:
1281 vNvrsFunc_Sha = BilinearInterp(vNvrsSha, mNrmGrid, nNrmGrid)
1282 vFunc_Sha = ValueFuncCRRA(vNvrsFunc_Sha, CRRA)
1283 else:
1284 vFunc_Sha = NullFunc()
1286 # Contribution share function
1287 if DiscreteShareBool:
1288 ShareFunc = DiscreteInterp(
1289 BilinearInterp(opt_idx, mNrmGrid, nNrmGrid), ShareGrid
1290 )
1291 else:
1292 ShareFunc = BilinearInterp(opt_Share, mNrmGrid, nNrmGrid)
1294 # Derivatives
1295 dvdmNvrsFunc_Sha = BilinearInterp(dvdmNvrsSha, mNrmGrid, nNrmGrid)
1296 dvdmFunc_Sha = MargValueFuncCRRA(dvdmNvrsFunc_Sha, CRRA)
1297 dvdnNvrsFunc_Sha = BilinearInterp(dvdnNvrsSha, mNrmGrid, nNrmGrid)
1298 dvdnFunc_Sha = MargValueFuncCRRA(dvdnNvrsFunc_Sha, CRRA)
1300 solution = RiskyContribShaSolution(
1301 vFunc_Adj=vFunc_Sha,
1302 ShareFunc_Adj=ShareFunc,
1303 dvdmFunc_Adj=dvdmFunc_Sha,
1304 dvdnFunc_Adj=dvdnFunc_Sha,
1305 # The fixed agent does nothing at this stage,
1306 # so his value functions are the next problem's
1307 vFunc_Fxd=vFunc_Cns_next,
1308 ShareFunc_Fxd=IdentityFunction(i_dim=2, n_dims=3),
1309 dvdmFunc_Fxd=dvdmFunc_Cns_next,
1310 dvdnFunc_Fxd=dvdnFunc_Cns_next,
1311 dvdsFunc_Fxd=dvdsFunc_Cns_next,
1312 )
1314 return solution
1317# Solver for the asset rebalancing stage
1318def solve_RiskyContrib_Reb(
1319 solution_next,
1320 CRRA,
1321 WithdrawTax,
1322 nNrmGrid,
1323 mNrmGrid,
1324 dfracGrid,
1325 vFuncBool,
1326 **unused_params,
1327):
1328 """
1329 Solves the asset-rebalancing-stage of the agent's problem
1331 Parameters
1332 ----------
1333 solution_next : RiskyContribShaSolution
1334 Solution to the income-contribution-share stage problem that follows.
1335 CRRA : float
1336 Coefficient of relative risk aversion.
1337 WithdrawTax : float
1338 Tax rate on risky asset withdrawals.
1339 nNrmGrid : numpy array
1340 Exogenous grid for risky resources.
1341 mNrmGrid : numpy array
1342 Exogenous grid for risk-free resources.
1343 dfracGrid : numpy array
1344 Grid for rebalancing flows. The final grid will be equivalent to
1345 [-nNrm*dfracGrid, dfracGrid*mNrm].
1346 vFuncBool : bool
1347 Determines whether the level of th value function must be computed.
1349 Returns
1350 -------
1351 solution : RiskyContribShaSolution
1352 Solution to the asset-rebalancing stage of the agent's problem.
1354 """
1355 # Extract next stage's solution
1356 vFunc_Adj_next = solution_next.vFunc_Adj
1357 dvdmFunc_Adj_next = solution_next.dvdmFunc_Adj
1358 dvdnFunc_Adj_next = solution_next.dvdnFunc_Adj
1360 vFunc_Fxd_next = solution_next.vFunc_Fxd
1361 dvdmFunc_Fxd_next = solution_next.dvdmFunc_Fxd
1362 dvdnFunc_Fxd_next = solution_next.dvdnFunc_Fxd
1363 dvdsFunc_Fxd_next = solution_next.dvdsFunc_Fxd
1365 uPinv = lambda x: utilityP_inv(x, CRRA)
1367 # Create tiled grids
1369 # Add 0 to the m and n grids
1370 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1371 nNrm_N = len(nNrmGrid)
1372 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1373 mNrm_N = len(mNrmGrid)
1374 d_N = len(dfracGrid)
1376 # Duplicate d so that possible values are -dfracGrid,dfracGrid. Duplicate 0 is
1377 # intentional since the tax causes a discontinuity. We need the value
1378 # from the left and right.
1379 dfracGrid = np.concatenate((-1 * np.flip(dfracGrid), dfracGrid))
1381 # It will be useful to pre-evaluate marginals at every (m,n,d) combination
1383 # Create tiled arrays for every d,m,n option
1384 d_N2 = len(dfracGrid)
1385 d_tiled, mNrm_tiled, nNrm_tiled = np.meshgrid(
1386 dfracGrid, mNrmGrid, nNrmGrid, indexing="ij"
1387 )
1389 # Get post-rebalancing assets.
1390 m_tilde, n_tilde = rebalance_assets(d_tiled, mNrm_tiled, nNrm_tiled, WithdrawTax)
1392 # Now the marginals, in inverse space
1393 dvdmNvrs = dvdmFunc_Adj_next.cFunc(m_tilde, n_tilde)
1394 dvdnNvrs = dvdnFunc_Adj_next.cFunc(m_tilde, n_tilde)
1396 # Pre-evaluate the inverse of (1-WithdrawTax)
1397 taxNvrs = uPinv(1 - WithdrawTax)
1398 # Create a tiled array of the tax
1399 taxNvrs_tiled = np.tile(
1400 np.reshape(
1401 np.concatenate([np.repeat(taxNvrs, d_N), np.ones(d_N, dtype=np.double)]),
1402 (d_N2, 1, 1),
1403 ),
1404 (1, mNrm_N, nNrm_N),
1405 )
1407 # The FOC is dvdn = tax*dvdm or dvdnNvrs = taxNvrs*dvdmNvrs
1408 dvdDNvrs = dvdnNvrs - taxNvrs_tiled * dvdmNvrs
1409 # The optimal d will be at the first point where dvdD < 0. The inverse
1410 # transformation flips the sign.
1412 # If the derivative is negative (inverse positive) at the lowest d,
1413 # then d == -1.0 is optimal
1414 constrained_bot = dvdDNvrs[0, :, :] >= 0.0
1415 # If it is positive (inverse negative) at the highest d, then d[-1] = 1.0
1416 # is optimal
1417 constrained_top = (
1418 dvdDNvrs[
1419 -1,
1420 :,
1421 :,
1422 ]
1423 <= 0.0
1424 )
1426 # Find indices at which the derivative crosses 0 for the 1st time
1427 # will be 0 if it never does, but "constrained_top/bot" deals with that
1428 crossings = np.logical_and(dvdDNvrs[:-1, :, :] <= 0.0, dvdDNvrs[1:, :, :] >= 0.0)
1429 idx = np.argmax(crossings, axis=0)
1431 m_idx_tiled, n_idx_tiled = np.meshgrid(
1432 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1433 )
1435 # Linearly interpolate the optimal withdrawal percentage d
1436 idx1 = idx + 1
1437 slopes = (
1438 dvdDNvrs[idx1, m_idx_tiled, n_idx_tiled]
1439 - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled]
1440 ) / (dfracGrid[idx1] - dfracGrid[idx])
1441 dfrac_opt = dfracGrid[idx] - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled] / slopes
1443 # Replace the ones we knew were constrained
1444 dfrac_opt[constrained_bot] = dfracGrid[0]
1445 dfrac_opt[constrained_top] = dfracGrid[-1]
1447 # Find m_tilde and n_tilde
1448 mtil_opt, ntil_opt = rebalance_assets(
1449 dfrac_opt, mNrm_tiled[0], nNrm_tiled[0], WithdrawTax
1450 )
1452 # Now the derivatives. These are not straight forward because of corner
1453 # solutions with partial derivatives that change the limits. The idea then
1454 # is to evaluate the possible uses of the marginal unit of resources and
1455 # take the maximum.
1457 # An additional unit of m
1458 marg_m = dvdmFunc_Adj_next(mtil_opt, ntil_opt)
1459 # An additional unit of n kept in n
1460 marg_n = dvdnFunc_Adj_next(mtil_opt, ntil_opt)
1461 # An additional unit of n withdrawn to m
1462 marg_n_to_m = marg_m * (1 - WithdrawTax)
1464 # Marginal value is the maximum of the marginals in their possible uses
1465 dvdm_Adj = np.maximum(marg_m, marg_n)
1466 dvdmNvrs_Adj = uPinv(dvdm_Adj)
1467 dvdn_Adj = np.maximum(marg_n, marg_n_to_m)
1468 dvdnNvrs_Adj = uPinv(dvdn_Adj)
1470 # Interpolators
1472 # Value
1473 if vFuncBool:
1474 vNvrs_Adj = vFunc_Adj_next.vFuncNvrs(mtil_opt, ntil_opt)
1475 vNvrsFunc_Adj = BilinearInterp(vNvrs_Adj, mNrmGrid, nNrmGrid)
1476 vFunc_Adj = ValueFuncCRRA(vNvrsFunc_Adj, CRRA)
1477 else:
1478 vFunc_Adj = NullFunc()
1480 # Marginals
1481 dvdmFunc_Adj = MargValueFuncCRRA(
1482 BilinearInterp(dvdmNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1483 )
1484 dvdnFunc_Adj = MargValueFuncCRRA(
1485 BilinearInterp(dvdnNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1486 )
1488 # Decison
1489 dfracFunc_Adj = BilinearInterp(dfrac_opt, mNrmGrid, nNrmGrid)
1491 solution = RiskyContribRebSolution(
1492 # Rebalancing stage adjusting
1493 vFunc_Adj=vFunc_Adj,
1494 dfracFunc_Adj=dfracFunc_Adj,
1495 dvdmFunc_Adj=dvdmFunc_Adj,
1496 dvdnFunc_Adj=dvdnFunc_Adj,
1497 # Rebalancing stage fixed (nothing happens, so value functions are
1498 # the ones from the next stage)
1499 vFunc_Fxd=vFunc_Fxd_next,
1500 dfracFunc_Fxd=ConstantFunction(0.0),
1501 dvdmFunc_Fxd=dvdmFunc_Fxd_next,
1502 dvdnFunc_Fxd=dvdnFunc_Fxd_next,
1503 dvdsFunc_Fxd=dvdsFunc_Fxd_next,
1504 )
1506 return solution
1509def solveRiskyContrib(
1510 solution_next,
1511 ShockDstn,
1512 IncShkDstn,
1513 RiskyDstn,
1514 IndepDstnBool,
1515 LivPrb,
1516 DiscFac,
1517 CRRA,
1518 Rfree,
1519 PermGroFac,
1520 WithdrawTax,
1521 BoroCnstArt,
1522 aXtraGrid,
1523 nNrmGrid,
1524 mNrmGrid,
1525 ShareGrid,
1526 dfracGrid,
1527 vFuncBool,
1528 AdjustPrb,
1529 DiscreteShareBool,
1530 joint_dist_solver,
1531):
1532 """
1533 Solve a full period (with its three stages) of the agent's problem
1535 Parameters
1536 ----------
1537 solution_next : RiskyContribSolution
1538 Solution to next period's problem.
1539 ShockDstn : DiscreteDistribution
1540 Joint distribution of next period's (0) permanent income shock, (1)
1541 transitory income shock, and (2) risky asset return factor.
1542 IncShkDstn : DiscreteDistribution
1543 Joint distribution of next period's (0) permanent income shock and (1)
1544 transitory income shock.
1545 RiskyDstn : DiscreteDistribution
1546 Distribution of next period's risky asset return factor.
1547 IndepDstnBool : bool
1548 Indicates whether the income and risky return distributions are
1549 independent.
1550 LivPrb : float
1551 Probability of surviving until next period.
1552 DiscFac : float
1553 Time-preference discount factor.
1554 CRRA : float
1555 Coefficient of relative risk aversion.
1556 Rfree : float
1557 Risk-free return factor.
1558 PermGroFac : float
1559 Deterministic permanent income growth factor.
1560 WithdrawTax : float
1561 Tax rate on risky asset withdrawals.
1562 BoroCnstArt : float
1563 Minimum allowed market resources (must be 0).
1564 aXtraGrid : numpy array
1565 Exogenous grid for end-of-period risk free resources.
1566 nNrmGrid : numpy array
1567 Exogenous grid for risky resources.
1568 mNrmGrid : numpy array
1569 Exogenous grid for risk-free resources.
1570 ShareGrid : numpy array
1571 Exogenous grid for the income contribution share.
1572 dfracGrid : numpy array
1573 Grid for rebalancing flows. The final grid will be equivalent to
1574 [-nNrm*dfracGrid, dfracGrid*mNrm].
1575 vFuncBool : bool
1576 Determines whether the level of th value function must be computed.
1577 AdjustPrb : float
1578 Probability that the agent will be able to rebalance his portfolio
1579 next period.
1580 DiscreteShareBool : bool
1581 Boolean that determines whether only a discrete set of contribution
1582 shares (ShareGrid) is allowed.
1583 joint_dist_solver: bool
1584 Should the general solver be used even if income and returns are
1585 independent?
1587 Returns
1588 -------
1589 periodSol : RiskyContribSolution
1590 Solution to the agent's current-period problem.
1592 """
1593 # Pack parameters to be passed to stage-specific solvers
1594 kws = {
1595 "ShockDstn": ShockDstn,
1596 "IncShkDstn": IncShkDstn,
1597 "RiskyDstn": RiskyDstn,
1598 "IndepDstnBool": IndepDstnBool,
1599 "LivPrb": LivPrb,
1600 "DiscFac": DiscFac,
1601 "CRRA": CRRA,
1602 "Rfree": Rfree,
1603 "PermGroFac": PermGroFac,
1604 "WithdrawTax": WithdrawTax,
1605 "BoroCnstArt": BoroCnstArt,
1606 "aXtraGrid": aXtraGrid,
1607 "nNrmGrid": nNrmGrid,
1608 "mNrmGrid": mNrmGrid,
1609 "ShareGrid": ShareGrid,
1610 "dfracGrid": dfracGrid,
1611 "vFuncBool": vFuncBool,
1612 "AdjustPrb": AdjustPrb,
1613 "DiscreteShareBool": DiscreteShareBool,
1614 "joint_dist_solver": joint_dist_solver,
1615 }
1617 # Stages of the problem in chronological order
1618 Stages = ["Reb", "Sha", "Cns"]
1619 n_stages = len(Stages)
1620 # Solvers, indexed by stage names
1621 Solvers = {
1622 "Reb": solve_RiskyContrib_Reb,
1623 "Sha": solve_RiskyContrib_Sha,
1624 "Cns": solve_RiskyContrib_Cns,
1625 }
1627 # Initialize empty solution
1628 stage_sols = {}
1629 # Solve stages backwards
1630 for i in reversed(range(n_stages)):
1631 stage = Stages[i]
1633 # In the last stage, the next solution is the first stage of the next
1634 # period. Otherwise, its the next stage of his period.
1635 if i == n_stages - 1:
1636 sol_next_stage = solution_next.stage_sols[Stages[0]]
1637 else:
1638 sol_next_stage = stage_sols[Stages[i + 1]]
1640 # Solve
1641 stage_sols[stage] = Solvers[stage](sol_next_stage, **kws)
1643 # Assemble stage solutions into period solution
1644 periodSol = RiskyContribSolution(**stage_sols)
1646 return periodSol
1649# %% Base risky-contrib dictionaries
1651risky_contrib_constructor_dict = {
1652 "IncShkDstn": construct_lognormal_income_process_unemployment,
1653 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1654 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1655 "aXtraGrid": make_assets_grid,
1656 "RiskyDstn": make_lognormal_RiskyDstn,
1657 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
1658 "ShareLimit": calc_ShareLimit_for_CRRA,
1659 "AdjustDstn": make_AdjustDstn,
1660 "solution_terminal": make_solution_terminal_risky_contrib,
1661 "ShareGrid": make_bounded_ShareGrid,
1662 "dfracGrid": make_simple_dGrid,
1663 "mNrmGrid": make_mNrm_grid,
1664 "nNrmGrid": make_nNrm_grid,
1665 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1666 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1667}
1669risky_contrib_params = {
1670 "constructors": risky_contrib_constructor_dict,
1671 # Preferences. The points of the model are more evident for more risk
1672 # averse and impatient agents
1673 "CRRA": 5.0,
1674 "DiscFac": 0.90,
1675 "WithdrawTax": [0.1],
1676 # Artificial borrowing constraint must be on
1677 "BoroCnstArt": 0.0,
1678 # Grids go up high wealth/P ratios and are less clustered at the bottom.
1679 "aXtraMax": 250,
1680 "aXtraCount": 50,
1681 "aXtraNestFac": 1,
1682 # Same goes for the new grids of the model
1683 "mNrmMin": 1e-6,
1684 "mNrmMax": 250,
1685 "mNrmCount": 50,
1686 "mNrmNestFac": 1,
1687 "nNrmMin": 1e-6,
1688 "nNrmMax": 250,
1689 "nNrmCount": 50,
1690 "nNrmNestFac": 1,
1691 # Income deduction/contribution share grid
1692 "ShareCount": 10,
1693 "ShareMax": 0.9,
1694 "DiscreteShareBool": False,
1695 # Grid for finding the optimal rebalancing flow
1696 "dCount": 20,
1697 "joint_dist_solver": False,
1698}
1699risky_asset_params = {
1700 # Risky return factor moments. Based on SP500 real returns from Shiller's
1701 # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
1702 "RiskyAvg": 1.080370891,
1703 "RiskyStd": 0.177196585,
1704 "ShareCount": 25, # Number of discrete points in the risky share approximation
1705 # Number of integration nodes to use in approximation of risky returns
1706 "RiskyCount": 5,
1707 # Probability that the agent can adjust their portfolio each period
1708 "AdjustPrb": [1.0],
1709 # When simulating the model, should all agents get the same risky return in
1710 # a given period?
1711 "sim_common_Rrisky": True,
1712}
1714# Infinite horizon version
1715init_risky_contrib = init_risky_asset.copy()
1716init_risky_contrib.update(risky_contrib_params)
1718# Lifecycle version
1719init_risky_contrib_lifecycle = init_lifecycle.copy()
1720init_risky_contrib_lifecycle.update(risky_asset_params)
1721init_risky_contrib_lifecycle.update(risky_contrib_params)
1723###############################################################################
1726class RiskyContribConsumerType(RiskyAssetConsumerType):
1727 """
1728 A consumer type with idiosyncratic shocks to permanent and transitory income,
1729 who can save in both a risk-free and a risky asset but faces frictions to
1730 moving funds between them. The agent can only consume out of his risk-free
1731 asset.
1733 The frictions are:
1735 - A proportional tax on funds moved from the risky to the risk-free
1736 asset.
1737 - A stochastic inability to move funds between his accounts.
1739 To partially avoid the second friction, the agent can commit to have a
1740 fraction of his labor income, which is usually deposited in his risk-free
1741 account, diverted to his risky account. He can change this fraction
1742 only in periods where he is able to move funds between accounts.
1743 """
1745 # The model is solved and simulated spliting each of the agent's
1746 # decisions into its own "stage". The stages in chronological order
1747 # are
1748 # - Reb: asset-rebalancing stage.
1749 # - Sha: definition of the income contribution share.
1750 # - Cns: consumption stage.
1751 stages = ["Reb", "Sha", "Cns"]
1752 # Each stage has its own states and controls, and its methods to find them.
1754 time_inv_ = RiskyAssetConsumerType.time_inv_ + [
1755 "DiscreteShareBool",
1756 "joint_dist_solver",
1757 "ShareGrid",
1758 "nNrmGrid",
1759 "mNrmGrid",
1760 "RiskyDstn",
1761 "dfracGrid",
1762 ]
1763 time_vary_ = RiskyAssetConsumerType.time_vary_ + ["WithdrawTax", "AdjustPrb"]
1765 # The new state variables (over those in ConsIndShock) are:
1766 # - nNrm: start-of-period risky resources.
1767 # - mNrmTilde: post-rebalancing risk-free resources.
1768 # - nNrmTilde: post-rebalancing risky resources.
1769 # - Share: income-deduction share.
1770 # For details, see
1771 # https://github.com/Mv77/RiskyContrib/blob/main/RiskyContrib.pdf
1772 state_vars = RiskyAssetConsumerType.state_vars + [
1773 "gNrm",
1774 "nNrm",
1775 "mNrmTilde",
1776 "nNrmTilde",
1777 "Share",
1778 ]
1779 shock_vars_ = RiskyAssetConsumerType.shock_vars_
1780 default_ = {
1781 "params": init_risky_contrib,
1782 "solver": solveRiskyContrib,
1783 "track_vars": ["aNrm", "cNrm", "mNrm", "nNrm", "dfrac", "Share", "pLvl"],
1784 }
1786 def __init__(self, **kwds):
1787 super().__init__(**kwds)
1788 # It looks like I can't assign this at the class level, unfortunately
1789 self.get_states = {
1790 "Reb": self.get_states_Reb,
1791 "Sha": self.get_states_Sha,
1792 "Cns": self.get_states_Cns,
1793 }
1794 self.get_controls = {
1795 "Reb": self.get_controls_Reb,
1796 "Sha": self.get_controls_Sha,
1797 "Cns": self.get_controls_Cns,
1798 }
1800 def pre_solve(self):
1801 self.construct("solution_terminal")
1803 def initialize_sim(self):
1804 """
1805 Initialize the state of simulation attributes.
1807 Parameters
1808 ----------
1809 None
1811 Returns
1812 -------
1813 None
1814 """
1815 RiskyAssetConsumerType.initialize_sim(self)
1816 self.state_now["Share"] = np.zeros(self.AgentCount)
1818 def sim_birth(self, which_agents):
1819 """
1820 Create new agents to replace ones who have recently died; takes draws of
1821 initial aNrm and pLvl, as in ConsIndShockModel, then sets Share, Adjust
1822 and post-rebalancing risky asset nNrmTilde to zero as initial values.
1823 Parameters
1824 ----------
1825 which_agents : np.array
1826 Boolean array of size AgentCount indicating which agents should be "born".
1828 Returns
1829 -------
1830 None
1831 """
1833 RiskyAssetConsumerType.sim_birth(self, which_agents)
1834 self.state_now["Share"][which_agents] = 0.0
1835 self.state_now["nNrmTilde"][which_agents] = 0.0
1837 def sim_one_period(self):
1838 """
1839 Simulates one period for this type.
1841 Has to be re-defined instead of using AgentType.sim_one_period() because
1842 of the "stages" structure.
1844 Parameters
1845 ----------
1846 None
1847 Returns
1848 -------
1849 None
1850 """
1852 if not hasattr(self, "solution"):
1853 raise Exception(
1854 "Model instance does not have a solution stored. To simulate, it is necessary"
1855 " to run the `solve()` method of the class first."
1856 )
1858 # Mortality adjusts the agent population
1859 self.get_mortality() # Replace some agents with "newborns"
1861 # Make state_now into state_prev, clearing state_now
1862 for var in self.state_now:
1863 self.state_prev[var] = self.state_now[var]
1865 if isinstance(self.state_now[var], np.ndarray):
1866 self.state_now[var] = np.empty(self.AgentCount)
1867 else:
1868 # Probably an aggregate variable. It may be getting set by the Market.
1869 pass
1871 if self.read_shocks: # If shock histories have been pre-specified, use those
1872 self.read_shocks_from_history()
1873 else: # Otherwise, draw shocks as usual according to subclass-specific method
1874 self.get_shocks()
1876 # Sequentially get states and controls of every stage
1877 for s in self.stages:
1878 self.get_states[s]()
1879 self.get_controls[s]()
1881 self.get_post_states()
1883 # Advance time for all agents
1884 self.t_age = self.t_age + 1 # Age all consumers by one period
1885 self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle
1886 self.t_cycle[self.t_cycle == self.T_cycle] = (
1887 0 # Resetting to zero for those who have reached the end
1888 )
1890 def get_states_Reb(self):
1891 """
1892 Get states for the first "stage": rebalancing.
1893 """
1895 pLvlPrev = self.state_prev["pLvl"]
1896 aNrmPrev = self.state_prev["aNrm"]
1897 SharePrev = self.state_prev["Share"]
1898 nNrmTildePrev = self.state_prev["nNrmTilde"]
1899 Rfree = np.array(self.Rfree)[self.t_cycle - 1]
1900 Rrisk = self.shocks["Risky"]
1902 # Calculate new states:
1904 # Permanent income
1905 self.state_now["pLvl"] = pLvlPrev * self.shocks["PermShk"]
1906 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1908 # Assets: mNrm and nNrm
1910 # Compute the effective growth factor of each asset
1911 RfEff = Rfree / self.shocks["PermShk"]
1912 RrEff = Rrisk / self.shocks["PermShk"]
1914 self.state_now["bNrm"] = RfEff * aNrmPrev # Liquid balances before labor income
1915 self.state_now["gNrm"] = (
1916 RrEff * nNrmTildePrev
1917 ) # Iliquid balances before labor income
1919 # Liquid balances after labor income
1920 self.state_now["mNrm"] = self.state_now["bNrm"] + self.shocks["TranShk"] * (
1921 1 - SharePrev
1922 )
1923 # Iliquid balances after labor income
1924 self.state_now["nNrm"] = (
1925 self.state_now["gNrm"] + self.shocks["TranShk"] * SharePrev
1926 )
1928 return None
1930 def get_controls_Reb(self):
1931 """
1932 Get controls for the first stage: rebalancing
1933 """
1934 dfrac = np.zeros(self.AgentCount) + np.nan
1936 # Loop over each period of the cycle, getting controls separately depending on "age"
1937 for t in range(self.T_cycle):
1938 # Find agents in this period-stage
1939 these = t == self.t_cycle
1941 # Get controls for agents who *can* adjust.
1942 those = np.logical_and(these, self.shocks["Adjust"])
1943 dfrac[those] = (
1944 self.solution[t]
1945 .stage_sols["Reb"]
1946 .dfracFunc_Adj(
1947 self.state_now["mNrm"][those], self.state_now["nNrm"][those]
1948 )
1949 )
1951 # Get Controls for agents who *can't* adjust.
1952 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
1953 dfrac[those] = (
1954 self.solution[t]
1955 .stage_sols["Reb"]
1956 .dfracFunc_Fxd(
1957 self.state_now["mNrm"][those],
1958 self.state_now["nNrm"][those],
1959 self.state_prev["Share"][those],
1960 )
1961 )
1963 # Limit dfrac to [-1,1] to prevent negative balances. Values outside
1964 # the range can come from extrapolation.
1965 self.controls["dfrac"] = np.minimum(np.maximum(dfrac, -1), 1.0)
1967 def get_states_Sha(self):
1968 """
1969 Get states for the second "stage": choosing the contribution share.
1970 """
1972 # Post-states are assets after rebalancing
1974 if "WithdrawTax" not in self.time_vary:
1975 mNrmTilde, nNrmTilde = rebalance_assets(
1976 self.controls["dfrac"],
1977 self.state_now["mNrm"],
1978 self.state_now["nNrm"],
1979 self.WithdrawTax,
1980 )
1982 else:
1983 # Initialize
1984 mNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1985 nNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1987 # Loop over each period of the cycle, getting controls separately depending on "age"
1988 for t in range(self.T_cycle):
1989 # Find agents in this period-stage
1990 these = t == self.t_cycle
1992 if np.sum(these) > 0:
1993 WithdrawTax = self.WithdrawTax[t]
1995 mNrmTilde[these], nNrmTilde[these] = rebalance_assets(
1996 self.controls["dfrac"][these],
1997 self.state_now["mNrm"][these],
1998 self.state_now["nNrm"][these],
1999 WithdrawTax,
2000 )
2002 self.state_now["mNrmTilde"] = mNrmTilde
2003 self.state_now["nNrmTilde"] = nNrmTilde
2005 def get_controls_Sha(self):
2006 """
2007 Get controls for the second "stage": choosing the contribution share.
2008 """
2010 Share = np.zeros(self.AgentCount) + np.nan
2012 # Loop over each period of the cycle, getting controls separately depending on "age"
2013 for t in range(self.T_cycle):
2014 # Find agents in this period-stage
2015 these = t == self.t_cycle
2017 # Get controls for agents who *can* adjust.
2018 those = np.logical_and(these, self.shocks["Adjust"])
2019 Share[those] = (
2020 self.solution[t]
2021 .stage_sols["Sha"]
2022 .ShareFunc_Adj(
2023 self.state_now["mNrmTilde"][those],
2024 self.state_now["nNrmTilde"][those],
2025 )
2026 )
2028 # Get Controls for agents who *can't* adjust.
2029 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
2030 Share[those] = (
2031 self.solution[t]
2032 .stage_sols["Sha"]
2033 .ShareFunc_Fxd(
2034 self.state_now["mNrmTilde"][those],
2035 self.state_now["nNrmTilde"][those],
2036 self.state_prev["Share"][those],
2037 )
2038 )
2040 # Store controls as attributes of self
2041 self.controls["Share"] = Share
2043 def get_states_Cns(self):
2044 """
2045 Get states for the third "stage": consumption.
2046 """
2048 # Contribution share becomes a state in the consumption problem
2049 self.state_now["Share"] = self.controls["Share"]
2051 def get_controls_Cns(self):
2052 """
2053 Get controls for the third "stage": consumption.
2054 """
2056 cNrm = np.zeros(self.AgentCount) + np.nan
2058 # Loop over each period of the cycle, getting controls separately depending on "age"
2059 for t in range(self.T_cycle):
2060 # Find agents in this period-stage
2061 these = t == self.t_cycle
2063 # Get consumption
2064 cNrm[these] = (
2065 self.solution[t]
2066 .stage_sols["Cns"]
2067 .cFunc(
2068 self.state_now["mNrmTilde"][these],
2069 self.state_now["nNrmTilde"][these],
2070 self.state_now["Share"][these],
2071 )
2072 )
2074 # Store controls as attributes of self
2075 # Since agents might be willing to end the period with a = 0, make
2076 # sure consumption does not go over m because of some numerical error.
2077 self.controls["cNrm"] = np.minimum(cNrm, self.state_now["mNrmTilde"])
2079 def get_post_states(self):
2080 """
2081 Set variables that are not a state to any problem but need to be
2082 computed in order to interact with shocks and produce next period's
2083 states.
2084 """
2085 self.state_now["aNrm"] = self.state_now["mNrmTilde"] - self.controls["cNrm"]