Coverage for HARK/ConsumptionSaving/ConsRiskyContribModel.py: 98%
506 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +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, tau):
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 tau : float
180 Tax rate of some kind.
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(tau) is list:
235 tau = tau[-1]
236 else:
237 tau = tau
239 # Value and marginal value function of the adjusting agent
240 vFunc_Reb_Adj_term = ValueFuncCRRA(lambda m, n: m + n / (1 + tau), CRRA)
241 dvdmFunc_Reb_Adj_term = MargValueFuncCRRA(lambda m, n: m + n / (1 + tau), CRRA)
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 + tau)
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, CRRA, tau, nNrmGrid, mNrmGrid, dfracGrid, vFuncBool, **unused_params
1316):
1317 """
1318 Solves the asset-rebalancing-stage of the agent's problem
1320 Parameters
1321 ----------
1322 solution_next : RiskyContribShaSolution
1323 Solution to the income-contribution-share stage problem that follows.
1324 CRRA : float
1325 Coefficient of relative risk aversion.
1326 tau : float
1327 Tax rate on risky asset withdrawals.
1328 nNrmGrid : numpy array
1329 Exogenous grid for risky resources.
1330 mNrmGrid : numpy array
1331 Exogenous grid for risk-free resources.
1332 dfracGrid : numpy array
1333 Grid for rebalancing flows. The final grid will be equivalent to
1334 [-nNrm*dfracGrid, dfracGrid*mNrm].
1335 vFuncBool : bool
1336 Determines whether the level of th value function must be computed.
1338 Returns
1339 -------
1340 solution : RiskyContribShaSolution
1341 Solution to the asset-rebalancing stage of the agent's problem.
1343 """
1344 # Extract next stage's solution
1345 vFunc_Adj_next = solution_next.vFunc_Adj
1346 dvdmFunc_Adj_next = solution_next.dvdmFunc_Adj
1347 dvdnFunc_Adj_next = solution_next.dvdnFunc_Adj
1349 vFunc_Fxd_next = solution_next.vFunc_Fxd
1350 dvdmFunc_Fxd_next = solution_next.dvdmFunc_Fxd
1351 dvdnFunc_Fxd_next = solution_next.dvdnFunc_Fxd
1352 dvdsFunc_Fxd_next = solution_next.dvdsFunc_Fxd
1354 uPinv = lambda x: utilityP_inv(x, CRRA)
1356 # Create tiled grids
1358 # Add 0 to the m and n grids
1359 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1360 nNrm_N = len(nNrmGrid)
1361 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1362 mNrm_N = len(mNrmGrid)
1363 d_N = len(dfracGrid)
1365 # Duplicate d so that possible values are -dfracGrid,dfracGrid. Duplicate 0 is
1366 # intentional since the tax causes a discontinuity. We need the value
1367 # from the left and right.
1368 dfracGrid = np.concatenate((-1 * np.flip(dfracGrid), dfracGrid))
1370 # It will be useful to pre-evaluate marginals at every (m,n,d) combination
1372 # Create tiled arrays for every d,m,n option
1373 d_N2 = len(dfracGrid)
1374 d_tiled, mNrm_tiled, nNrm_tiled = np.meshgrid(
1375 dfracGrid, mNrmGrid, nNrmGrid, indexing="ij"
1376 )
1378 # Get post-rebalancing assets.
1379 m_tilde, n_tilde = rebalance_assets(d_tiled, mNrm_tiled, nNrm_tiled, tau)
1381 # Now the marginals, in inverse space
1382 dvdmNvrs = dvdmFunc_Adj_next.cFunc(m_tilde, n_tilde)
1383 dvdnNvrs = dvdnFunc_Adj_next.cFunc(m_tilde, n_tilde)
1385 # Pre-evaluate the inverse of (1-tau)
1386 taxNvrs = uPinv(1 - tau)
1387 # Create a tiled array of the tax
1388 taxNvrs_tiled = np.tile(
1389 np.reshape(
1390 np.concatenate([np.repeat(taxNvrs, d_N), np.ones(d_N, dtype=np.double)]),
1391 (d_N2, 1, 1),
1392 ),
1393 (1, mNrm_N, nNrm_N),
1394 )
1396 # The FOC is dvdn = tax*dvdm or dvdnNvrs = taxNvrs*dvdmNvrs
1397 dvdDNvrs = dvdnNvrs - taxNvrs_tiled * dvdmNvrs
1398 # The optimal d will be at the first point where dvdD < 0. The inverse
1399 # transformation flips the sign.
1401 # If the derivative is negative (inverse positive) at the lowest d,
1402 # then d == -1.0 is optimal
1403 constrained_bot = dvdDNvrs[0, :, :] >= 0.0
1404 # If it is positive (inverse negative) at the highest d, then d[-1] = 1.0
1405 # is optimal
1406 constrained_top = (
1407 dvdDNvrs[
1408 -1,
1409 :,
1410 :,
1411 ]
1412 <= 0.0
1413 )
1415 # Find indices at which the derivative crosses 0 for the 1st time
1416 # will be 0 if it never does, but "constrained_top/bot" deals with that
1417 crossings = np.logical_and(dvdDNvrs[:-1, :, :] <= 0.0, dvdDNvrs[1:, :, :] >= 0.0)
1418 idx = np.argmax(crossings, axis=0)
1420 m_idx_tiled, n_idx_tiled = np.meshgrid(
1421 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1422 )
1424 # Linearly interpolate the optimal withdrawal percentage d
1425 idx1 = idx + 1
1426 slopes = (
1427 dvdDNvrs[idx1, m_idx_tiled, n_idx_tiled]
1428 - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled]
1429 ) / (dfracGrid[idx1] - dfracGrid[idx])
1430 dfrac_opt = dfracGrid[idx] - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled] / slopes
1432 # Replace the ones we knew were constrained
1433 dfrac_opt[constrained_bot] = dfracGrid[0]
1434 dfrac_opt[constrained_top] = dfracGrid[-1]
1436 # Find m_tilde and n_tilde
1437 mtil_opt, ntil_opt = rebalance_assets(dfrac_opt, mNrm_tiled[0], nNrm_tiled[0], tau)
1439 # Now the derivatives. These are not straight forward because of corner
1440 # solutions with partial derivatives that change the limits. The idea then
1441 # is to evaluate the possible uses of the marginal unit of resources and
1442 # take the maximum.
1444 # An additional unit of m
1445 marg_m = dvdmFunc_Adj_next(mtil_opt, ntil_opt)
1446 # An additional unit of n kept in n
1447 marg_n = dvdnFunc_Adj_next(mtil_opt, ntil_opt)
1448 # An additional unit of n withdrawn to m
1449 marg_n_to_m = marg_m * (1 - tau)
1451 # Marginal value is the maximum of the marginals in their possible uses
1452 dvdm_Adj = np.maximum(marg_m, marg_n)
1453 dvdmNvrs_Adj = uPinv(dvdm_Adj)
1454 dvdn_Adj = np.maximum(marg_n, marg_n_to_m)
1455 dvdnNvrs_Adj = uPinv(dvdn_Adj)
1457 # Interpolators
1459 # Value
1460 if vFuncBool:
1461 vNvrs_Adj = vFunc_Adj_next.vFuncNvrs(mtil_opt, ntil_opt)
1462 vNvrsFunc_Adj = BilinearInterp(vNvrs_Adj, mNrmGrid, nNrmGrid)
1463 vFunc_Adj = ValueFuncCRRA(vNvrsFunc_Adj, CRRA)
1464 else:
1465 vFunc_Adj = NullFunc()
1467 # Marginals
1468 dvdmFunc_Adj = MargValueFuncCRRA(
1469 BilinearInterp(dvdmNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1470 )
1471 dvdnFunc_Adj = MargValueFuncCRRA(
1472 BilinearInterp(dvdnNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1473 )
1475 # Decison
1476 dfracFunc_Adj = BilinearInterp(dfrac_opt, mNrmGrid, nNrmGrid)
1478 solution = RiskyContribRebSolution(
1479 # Rebalancing stage adjusting
1480 vFunc_Adj=vFunc_Adj,
1481 dfracFunc_Adj=dfracFunc_Adj,
1482 dvdmFunc_Adj=dvdmFunc_Adj,
1483 dvdnFunc_Adj=dvdnFunc_Adj,
1484 # Rebalancing stage fixed (nothing happens, so value functions are
1485 # the ones from the next stage)
1486 vFunc_Fxd=vFunc_Fxd_next,
1487 dfracFunc_Fxd=ConstantFunction(0.0),
1488 dvdmFunc_Fxd=dvdmFunc_Fxd_next,
1489 dvdnFunc_Fxd=dvdnFunc_Fxd_next,
1490 dvdsFunc_Fxd=dvdsFunc_Fxd_next,
1491 )
1493 return solution
1496def solveRiskyContrib(
1497 solution_next,
1498 ShockDstn,
1499 IncShkDstn,
1500 RiskyDstn,
1501 IndepDstnBool,
1502 LivPrb,
1503 DiscFac,
1504 CRRA,
1505 Rfree,
1506 PermGroFac,
1507 tau,
1508 BoroCnstArt,
1509 aXtraGrid,
1510 nNrmGrid,
1511 mNrmGrid,
1512 ShareGrid,
1513 dfracGrid,
1514 vFuncBool,
1515 AdjustPrb,
1516 DiscreteShareBool,
1517 joint_dist_solver,
1518):
1519 """
1520 Solve a full period (with its three stages) of the agent's problem
1522 Parameters
1523 ----------
1524 solution_next : RiskyContribSolution
1525 Solution to next period's problem.
1526 ShockDstn : DiscreteDistribution
1527 Joint distribution of next period's (0) permanent income shock, (1)
1528 transitory income shock, and (2) risky asset return factor.
1529 IncShkDstn : DiscreteDistribution
1530 Joint distribution of next period's (0) permanent income shock and (1)
1531 transitory income shock.
1532 RiskyDstn : DiscreteDistribution
1533 Distribution of next period's risky asset return factor.
1534 IndepDstnBool : bool
1535 Indicates whether the income and risky return distributions are
1536 independent.
1537 LivPrb : float
1538 Probability of surviving until next period.
1539 DiscFac : float
1540 Time-preference discount factor.
1541 CRRA : float
1542 Coefficient of relative risk aversion.
1543 Rfree : float
1544 Risk-free return factor.
1545 PermGroFac : float
1546 Deterministic permanent income growth factor.
1547 tau : float
1548 Tax rate on risky asset withdrawals.
1549 BoroCnstArt : float
1550 Minimum allowed market resources (must be 0).
1551 aXtraGrid : numpy array
1552 Exogenous grid for end-of-period risk free resources.
1553 nNrmGrid : numpy array
1554 Exogenous grid for risky resources.
1555 mNrmGrid : numpy array
1556 Exogenous grid for risk-free resources.
1557 ShareGrid : numpy array
1558 Exogenous grid for the income contribution share.
1559 dfracGrid : numpy array
1560 Grid for rebalancing flows. The final grid will be equivalent to
1561 [-nNrm*dfracGrid, dfracGrid*mNrm].
1562 vFuncBool : bool
1563 Determines whether the level of th value function must be computed.
1564 AdjustPrb : float
1565 Probability that the agent will be able to rebalance his portfolio
1566 next period.
1567 DiscreteShareBool : bool
1568 Boolean that determines whether only a discrete set of contribution
1569 shares (ShareGrid) is allowed.
1570 joint_dist_solver: bool
1571 Should the general solver be used even if income and returns are
1572 independent?
1574 Returns
1575 -------
1576 periodSol : RiskyContribSolution
1577 Solution to the agent's current-period problem.
1579 """
1580 # Pack parameters to be passed to stage-specific solvers
1581 kws = {
1582 "ShockDstn": ShockDstn,
1583 "IncShkDstn": IncShkDstn,
1584 "RiskyDstn": RiskyDstn,
1585 "IndepDstnBool": IndepDstnBool,
1586 "LivPrb": LivPrb,
1587 "DiscFac": DiscFac,
1588 "CRRA": CRRA,
1589 "Rfree": Rfree,
1590 "PermGroFac": PermGroFac,
1591 "tau": tau,
1592 "BoroCnstArt": BoroCnstArt,
1593 "aXtraGrid": aXtraGrid,
1594 "nNrmGrid": nNrmGrid,
1595 "mNrmGrid": mNrmGrid,
1596 "ShareGrid": ShareGrid,
1597 "dfracGrid": dfracGrid,
1598 "vFuncBool": vFuncBool,
1599 "AdjustPrb": AdjustPrb,
1600 "DiscreteShareBool": DiscreteShareBool,
1601 "joint_dist_solver": joint_dist_solver,
1602 }
1604 # Stages of the problem in chronological order
1605 Stages = ["Reb", "Sha", "Cns"]
1606 n_stages = len(Stages)
1607 # Solvers, indexed by stage names
1608 Solvers = {
1609 "Reb": solve_RiskyContrib_Reb,
1610 "Sha": solve_RiskyContrib_Sha,
1611 "Cns": solve_RiskyContrib_Cns,
1612 }
1614 # Initialize empty solution
1615 stage_sols = {}
1616 # Solve stages backwards
1617 for i in reversed(range(n_stages)):
1618 stage = Stages[i]
1620 # In the last stage, the next solution is the first stage of the next
1621 # period. Otherwise, its the next stage of his period.
1622 if i == n_stages - 1:
1623 sol_next_stage = solution_next.stage_sols[Stages[0]]
1624 else:
1625 sol_next_stage = stage_sols[Stages[i + 1]]
1627 # Solve
1628 stage_sols[stage] = Solvers[stage](sol_next_stage, **kws)
1630 # Assemble stage solutions into period solution
1631 periodSol = RiskyContribSolution(**stage_sols)
1633 return periodSol
1636# %% Base risky-contrib dictionaries
1638risky_contrib_constructor_dict = {
1639 "IncShkDstn": construct_lognormal_income_process_unemployment,
1640 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1641 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1642 "aXtraGrid": make_assets_grid,
1643 "RiskyDstn": make_lognormal_RiskyDstn,
1644 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
1645 "ShareLimit": calc_ShareLimit_for_CRRA,
1646 "AdjustDstn": make_AdjustDstn,
1647 "solution_terminal": make_solution_terminal_risky_contrib,
1648 "ShareGrid": make_bounded_ShareGrid,
1649 "dfracGrid": make_simple_dGrid,
1650 "mNrmGrid": make_mNrm_grid,
1651 "nNrmGrid": make_nNrm_grid,
1652 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1653 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1654}
1656risky_contrib_params = {
1657 "constructors": risky_contrib_constructor_dict,
1658 # Preferences. The points of the model are more evident for more risk
1659 # averse and impatient agents
1660 "CRRA": 5.0,
1661 "DiscFac": 0.90,
1662 # Artificial borrowing constraint must be on
1663 "BoroCnstArt": 0.0,
1664 # Grids go up high wealth/P ratios and are less clustered at the bottom.
1665 "aXtraMax": 250,
1666 "aXtraCount": 50,
1667 "aXtraNestFac": 1,
1668 # Same goes for the new grids of the model
1669 "mNrmMin": 1e-6,
1670 "mNrmMax": 250,
1671 "mNrmCount": 50,
1672 "mNrmNestFac": 1,
1673 "nNrmMin": 1e-6,
1674 "nNrmMax": 250,
1675 "nNrmCount": 50,
1676 "nNrmNestFac": 1,
1677 # Income deduction/contribution share grid
1678 "ShareCount": 10,
1679 "ShareMax": 0.9,
1680 "DiscreteShareBool": False,
1681 # Grid for finding the optimal rebalancing flow
1682 "dCount": 20,
1683 "joint_dist_solver": False,
1684}
1685risky_asset_params = {
1686 # Risky return factor moments. Based on SP500 real returns from Shiller's
1687 # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
1688 "RiskyAvg": 1.080370891,
1689 "RiskyStd": 0.177196585,
1690 "ShareCount": 25, # Number of discrete points in the risky share approximation
1691 # Number of integration nodes to use in approximation of risky returns
1692 "RiskyCount": 5,
1693 # Probability that the agent can adjust their portfolio each period
1694 "AdjustPrb": [1.0],
1695 # When simulating the model, should all agents get the same risky return in
1696 # a given period?
1697 "sim_common_Rrisky": True,
1698}
1700# Infinite horizon version
1701init_risky_contrib = init_risky_asset.copy()
1702init_risky_contrib.update(risky_contrib_params)
1704# Lifecycle version
1705init_risky_contrib_lifecycle = init_lifecycle.copy()
1706init_risky_contrib_lifecycle.update(risky_asset_params)
1707init_risky_contrib_lifecycle.update(risky_contrib_params)
1709###############################################################################
1712class RiskyContribConsumerType(RiskyAssetConsumerType):
1713 """
1714 A consumer type with idiosyncratic shocks to permanent and transitory income,
1715 who can save in both a risk-free and a risky asset but faces frictions to
1716 moving funds between them. The agent can only consume out of his risk-free
1717 asset.
1719 The frictions are:
1721 - A proportional tax on funds moved from the risky to the risk-free
1722 asset.
1723 - A stochastic inability to move funds between his accounts.
1725 To partially avoid the second friction, the agent can commit to have a
1726 fraction of his labor income, which is usually deposited in his risk-free
1727 account, diverted to his risky account. He can change this fraction
1728 only in periods where he is able to move funds between accounts.
1729 """
1731 # The model is solved and simulated spliting each of the agent's
1732 # decisions into its own "stage". The stages in chronological order
1733 # are
1734 # - Reb: asset-rebalancing stage.
1735 # - Sha: definition of the income contribution share.
1736 # - Cns: consumption stage.
1737 stages = ["Reb", "Sha", "Cns"]
1738 # Each stage has its own states and controls, and its methods to find them.
1740 time_inv_ = RiskyAssetConsumerType.time_inv_ + [
1741 "DiscreteShareBool",
1742 "joint_dist_solver",
1743 "ShareGrid",
1744 "nNrmGrid",
1745 "mNrmGrid",
1746 "RiskyDstn",
1747 "dfracGrid",
1748 ]
1749 time_vary_ = RiskyAssetConsumerType.time_vary_ + ["tau", "AdjustPrb"]
1751 # The new state variables (over those in ConsIndShock) are:
1752 # - nNrm: start-of-period risky resources.
1753 # - mNrmTilde: post-rebalancing risk-free resources.
1754 # - nNrmTilde: post-rebalancing risky resources.
1755 # - Share: income-deduction share.
1756 # For details, see
1757 # https://github.com/Mv77/RiskyContrib/blob/main/RiskyContrib.pdf
1758 state_vars = RiskyAssetConsumerType.state_vars + [
1759 "gNrm",
1760 "nNrm",
1761 "mNrmTilde",
1762 "nNrmTilde",
1763 "Share",
1764 ]
1765 shock_vars_ = RiskyAssetConsumerType.shock_vars_
1766 default_ = {"params": init_risky_contrib, "solver": solveRiskyContrib}
1768 def __init__(self, **kwds):
1769 super().__init__(**kwds)
1770 # It looks like I can't assign this at the class level, unfortunately
1771 self.get_states = {
1772 "Reb": self.get_states_Reb,
1773 "Sha": self.get_states_Sha,
1774 "Cns": self.get_states_Cns,
1775 }
1776 self.get_controls = {
1777 "Reb": self.get_controls_Reb,
1778 "Sha": self.get_controls_Sha,
1779 "Cns": self.get_controls_Cns,
1780 }
1782 def pre_solve(self):
1783 self.construct("solution_terminal")
1785 def initialize_sim(self):
1786 """
1787 Initialize the state of simulation attributes.
1789 Parameters
1790 ----------
1791 None
1793 Returns
1794 -------
1795 None
1796 """
1797 RiskyAssetConsumerType.initialize_sim(self)
1798 self.state_now["Share"] = np.zeros(self.AgentCount)
1800 def sim_birth(self, which_agents):
1801 """
1802 Create new agents to replace ones who have recently died; takes draws of
1803 initial aNrm and pLvl, as in ConsIndShockModel, then sets Share, Adjust
1804 and post-rebalancing risky asset nNrmTilde to zero as initial values.
1805 Parameters
1806 ----------
1807 which_agents : np.array
1808 Boolean array of size AgentCount indicating which agents should be "born".
1810 Returns
1811 -------
1812 None
1813 """
1815 RiskyAssetConsumerType.sim_birth(self, which_agents)
1816 self.state_now["Share"][which_agents] = 0.0
1817 self.state_now["nNrmTilde"][which_agents] = 0.0
1819 def sim_one_period(self):
1820 """
1821 Simulates one period for this type.
1823 Has to be re-defined instead of using AgentType.sim_one_period() because
1824 of the "stages" structure.
1826 Parameters
1827 ----------
1828 None
1829 Returns
1830 -------
1831 None
1832 """
1834 if not hasattr(self, "solution"):
1835 raise Exception(
1836 "Model instance does not have a solution stored. To simulate, it is necessary"
1837 " to run the `solve()` method of the class first."
1838 )
1840 # Mortality adjusts the agent population
1841 self.get_mortality() # Replace some agents with "newborns"
1843 # Make state_now into state_prev, clearing state_now
1844 for var in self.state_now:
1845 self.state_prev[var] = self.state_now[var]
1847 if isinstance(self.state_now[var], np.ndarray):
1848 self.state_now[var] = np.empty(self.AgentCount)
1849 else:
1850 # Probably an aggregate variable. It may be getting set by the Market.
1851 pass
1853 if self.read_shocks: # If shock histories have been pre-specified, use those
1854 self.read_shocks_from_history()
1855 else: # Otherwise, draw shocks as usual according to subclass-specific method
1856 self.get_shocks()
1858 # Sequentially get states and controls of every stage
1859 for s in self.stages:
1860 self.get_states[s]()
1861 self.get_controls[s]()
1863 self.get_post_states()
1865 # Advance time for all agents
1866 self.t_age = self.t_age + 1 # Age all consumers by one period
1867 self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle
1868 self.t_cycle[self.t_cycle == self.T_cycle] = (
1869 0 # Resetting to zero for those who have reached the end
1870 )
1872 def get_states_Reb(self):
1873 """
1874 Get states for the first "stage": rebalancing.
1875 """
1877 pLvlPrev = self.state_prev["pLvl"]
1878 aNrmPrev = self.state_prev["aNrm"]
1879 SharePrev = self.state_prev["Share"]
1880 nNrmTildePrev = self.state_prev["nNrmTilde"]
1881 Rfree = self.get_Rfree()
1882 Rrisk = self.shocks["Risky"]
1884 # Calculate new states:
1886 # Permanent income
1887 self.state_now["pLvl"] = pLvlPrev * self.shocks["PermShk"]
1888 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1890 # Assets: mNrm and nNrm
1892 # Compute the effective growth factor of each asset
1893 RfEff = Rfree / self.shocks["PermShk"]
1894 RrEff = Rrisk / self.shocks["PermShk"]
1896 self.state_now["bNrm"] = RfEff * aNrmPrev # Liquid balances before labor income
1897 self.state_now["gNrm"] = (
1898 RrEff * nNrmTildePrev
1899 ) # Iliquid balances before labor income
1901 # Liquid balances after labor income
1902 self.state_now["mNrm"] = self.state_now["bNrm"] + self.shocks["TranShk"] * (
1903 1 - SharePrev
1904 )
1905 # Iliquid balances after labor income
1906 self.state_now["nNrm"] = (
1907 self.state_now["gNrm"] + self.shocks["TranShk"] * SharePrev
1908 )
1910 return None
1912 def get_controls_Reb(self):
1913 """
1914 Get controls for the first stage: rebalancing
1915 """
1916 dfrac = np.zeros(self.AgentCount) + np.nan
1918 # Loop over each period of the cycle, getting controls separately depending on "age"
1919 for t in range(self.T_cycle):
1920 # Find agents in this period-stage
1921 these = t == self.t_cycle
1923 # Get controls for agents who *can* adjust.
1924 those = np.logical_and(these, self.shocks["Adjust"])
1925 dfrac[those] = (
1926 self.solution[t]
1927 .stage_sols["Reb"]
1928 .dfracFunc_Adj(
1929 self.state_now["mNrm"][those], self.state_now["nNrm"][those]
1930 )
1931 )
1933 # Get Controls for agents who *can't* adjust.
1934 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
1935 dfrac[those] = (
1936 self.solution[t]
1937 .stage_sols["Reb"]
1938 .dfracFunc_Fxd(
1939 self.state_now["mNrm"][those],
1940 self.state_now["nNrm"][those],
1941 self.state_prev["Share"][those],
1942 )
1943 )
1945 # Limit dfrac to [-1,1] to prevent negative balances. Values outside
1946 # the range can come from extrapolation.
1947 self.controls["dfrac"] = np.minimum(np.maximum(dfrac, -1), 1.0)
1949 def get_states_Sha(self):
1950 """
1951 Get states for the second "stage": choosing the contribution share.
1952 """
1954 # Post-states are assets after rebalancing
1956 if "tau" not in self.time_vary:
1957 mNrmTilde, nNrmTilde = rebalance_assets(
1958 self.controls["dfrac"],
1959 self.state_now["mNrm"],
1960 self.state_now["nNrm"],
1961 self.tau,
1962 )
1964 else:
1965 # Initialize
1966 mNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1967 nNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1969 # Loop over each period of the cycle, getting controls separately depending on "age"
1970 for t in range(self.T_cycle):
1971 # Find agents in this period-stage
1972 these = t == self.t_cycle
1974 if np.sum(these) > 0:
1975 tau = self.tau[t]
1977 mNrmTilde[these], nNrmTilde[these] = rebalance_assets(
1978 self.controls["dfrac"][these],
1979 self.state_now["mNrm"][these],
1980 self.state_now["nNrm"][these],
1981 tau,
1982 )
1984 self.state_now["mNrmTilde"] = mNrmTilde
1985 self.state_now["nNrmTilde"] = nNrmTilde
1987 def get_controls_Sha(self):
1988 """
1989 Get controls for the second "stage": choosing the contribution share.
1990 """
1992 Share = np.zeros(self.AgentCount) + np.nan
1994 # Loop over each period of the cycle, getting controls separately depending on "age"
1995 for t in range(self.T_cycle):
1996 # Find agents in this period-stage
1997 these = t == self.t_cycle
1999 # Get controls for agents who *can* adjust.
2000 those = np.logical_and(these, self.shocks["Adjust"])
2001 Share[those] = (
2002 self.solution[t]
2003 .stage_sols["Sha"]
2004 .ShareFunc_Adj(
2005 self.state_now["mNrmTilde"][those],
2006 self.state_now["nNrmTilde"][those],
2007 )
2008 )
2010 # Get Controls for agents who *can't* adjust.
2011 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
2012 Share[those] = (
2013 self.solution[t]
2014 .stage_sols["Sha"]
2015 .ShareFunc_Fxd(
2016 self.state_now["mNrmTilde"][those],
2017 self.state_now["nNrmTilde"][those],
2018 self.state_prev["Share"][those],
2019 )
2020 )
2022 # Store controls as attributes of self
2023 self.controls["Share"] = Share
2025 def get_states_Cns(self):
2026 """
2027 Get states for the third "stage": consumption.
2028 """
2030 # Contribution share becomes a state in the consumption problem
2031 self.state_now["Share"] = self.controls["Share"]
2033 def get_controls_Cns(self):
2034 """
2035 Get controls for the third "stage": consumption.
2036 """
2038 cNrm = np.zeros(self.AgentCount) + np.nan
2040 # Loop over each period of the cycle, getting controls separately depending on "age"
2041 for t in range(self.T_cycle):
2042 # Find agents in this period-stage
2043 these = t == self.t_cycle
2045 # Get consumption
2046 cNrm[these] = (
2047 self.solution[t]
2048 .stage_sols["Cns"]
2049 .cFunc(
2050 self.state_now["mNrmTilde"][these],
2051 self.state_now["nNrmTilde"][these],
2052 self.state_now["Share"][these],
2053 )
2054 )
2056 # Store controls as attributes of self
2057 # Since agents might be willing to end the period with a = 0, make
2058 # sure consumption does not go over m because of some numerical error.
2059 self.controls["cNrm"] = np.minimum(cNrm, self.state_now["mNrmTilde"])
2061 def get_post_states(self):
2062 """
2063 Set variables that are not a state to any problem but need to be
2064 computed in order to interact with shocks and produce next period's
2065 states.
2066 """
2067 self.state_now["aNrm"] = self.state_now["mNrmTilde"] - self.controls["cNrm"]