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