Coverage for HARK / ConsumptionSaving / ConsRiskyContribModel.py: 98%
506 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-07 05:16 +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]
236 else:
237 WithdrawTax = WithdrawTax
239 # Value and marginal value function of the adjusting agent
240 vFunc_Reb_Adj_term = ValueFuncCRRA(lambda m, n: m + n / (1 + WithdrawTax), CRRA)
241 dvdmFunc_Reb_Adj_term = MargValueFuncCRRA(
242 lambda m, n: m + n / (1 + WithdrawTax), CRRA
243 )
244 # A marginal unit of n will be withdrawn and put into m. Then consumed.
245 dvdnFunc_Reb_Adj_term = lambda m, n: dvdmFunc_Reb_Adj_term(m, n) / (1 + WithdrawTax)
247 Reb_stage_sol = RiskyContribRebSolution(
248 # Rebalancing stage
249 vFunc_Adj=vFunc_Reb_Adj_term,
250 dfracFunc_Adj=dfracFunc_Adj_term,
251 dvdmFunc_Adj=dvdmFunc_Reb_Adj_term,
252 dvdnFunc_Adj=dvdnFunc_Reb_Adj_term,
253 # Adjusting stage
254 vFunc_Fxd=vFunc_Cns_term,
255 dfracFunc_Fxd=ConstantFunction(0.0),
256 dvdmFunc_Fxd=dvdmFunc_Cns_term,
257 dvdnFunc_Fxd=dvdnFunc_Cns_term,
258 dvdsFunc_Fxd=dvdsFunc_Cns_term,
259 )
261 # Construct the terminal period solution
262 solution_terminal = RiskyContribSolution(
263 Reb_stage_sol, Sha_stage_sol, Cns_stage_sol
264 )
265 return solution_terminal
268###############################################################################
270# %% Classes for RiskyContrib type solution objects
273# Class for asset adjustment stage solution
274class RiskyContribRebSolution(MetricObject):
275 """
276 A class for representing the solution to the asset-rebalancing stage of
277 the 'RiskyContrib' model.
279 Parameters
280 ----------
281 vFunc_Adj : ValueFunc2D
282 Stage value function over normalized liquid resources and normalized
283 iliquid resources when the agent is able to adjust his portfolio.
284 dfracFunc_Adj : Interp2D
285 Deposit function over normalized liquid resources and normalized
286 iliquid resources when the agent is able to adjust his portfolio.
287 dvdmFunc_Adj : MargValueFunc2D
288 Marginal value over normalized liquid resources when the agent is able
289 to adjust his portfolio.
290 dvdnFunc_Adj : MargValueFunc2D
291 Marginal value over normalized liquid resources when the agent is able
292 to adjust his portfolio.
293 vFunc_Fxd : ValueFunc3D
294 Stage value function over normalized liquid resources, normalized
295 iliquid resources, and income contribution share when the agent is
296 not able to adjust his portfolio.
297 dfracFunc_Fxd : Interp2D
298 Deposit function over normalized liquid resources, normalized iliquid
299 resources, and income contribution share when the agent is not able to
300 adjust his portfolio.
301 Must be ConstantFunction(0.0)
302 dvdmFunc_Fxd : MargValueFunc3D
303 Marginal value over normalized liquid resources when the agent is not
304 able to adjust his portfolio.
305 dvdnFunc_Fxd : MargValueFunc3D
306 Marginal value over normalized iliquid resources when the agent is not
307 able to adjust his portfolio.
308 dvdsFunc_Fxd : Interp3D
309 Marginal value function over income contribution share when the agent
310 is not able to ajust his portfolio.
311 """
313 distance_criteria = ["dvdmFunc_Adj", "dvdnFunc_Adj"]
315 def __init__(
316 self,
317 # Rebalancing stage, adjusting
318 vFunc_Adj=None,
319 dfracFunc_Adj=None,
320 dvdmFunc_Adj=None,
321 dvdnFunc_Adj=None,
322 # Rebalancing stage, fixed
323 vFunc_Fxd=None,
324 dfracFunc_Fxd=None,
325 dvdmFunc_Fxd=None,
326 dvdnFunc_Fxd=None,
327 dvdsFunc_Fxd=None,
328 ):
329 # Rebalancing stage
330 if vFunc_Adj is None:
331 vFunc_Adj = NullFunc()
332 if dfracFunc_Adj is None:
333 dfracFunc_Adj = NullFunc()
334 if dvdmFunc_Adj is None:
335 dvdmFunc_Adj = NullFunc()
336 if dvdnFunc_Adj is None:
337 dvdnFunc_Adj = NullFunc()
339 if vFunc_Fxd is None:
340 vFunc_Fxd = NullFunc()
341 if dfracFunc_Fxd is None:
342 dfracFunc_Fxd = NullFunc()
343 if dvdmFunc_Fxd is None:
344 dvdmFunc_Fxd = NullFunc()
345 if dvdnFunc_Fxd is None:
346 dvdnFunc_Fxd = NullFunc()
347 if dvdsFunc_Fxd is None:
348 dvdsFunc_Fxd = NullFunc()
350 # Components of the adjusting problem
351 self.vFunc_Adj = vFunc_Adj
352 self.dfracFunc_Adj = dfracFunc_Adj
353 self.dvdmFunc_Adj = dvdmFunc_Adj
354 self.dvdnFunc_Adj = dvdnFunc_Adj
356 # Components of the fixed problem
357 self.vFunc_Fxd = vFunc_Fxd
358 self.dfracFunc_Fxd = dfracFunc_Fxd
359 self.dvdmFunc_Fxd = dvdmFunc_Fxd
360 self.dvdnFunc_Fxd = dvdnFunc_Fxd
361 self.dvdsFunc_Fxd = dvdsFunc_Fxd
364# Class for the contribution share stage solution
365class RiskyContribShaSolution(MetricObject):
366 """
367 A class for representing the solution to the contribution-share stage of
368 the 'RiskyContrib' model.
370 Parameters
371 ----------
372 vFunc_Adj : ValueFunc2D
373 Stage value function over normalized liquid resources and normalized
374 iliquid resources when the agent is able to adjust his portfolio.
375 ShareFunc_Adj : Interp2D
376 Income contribution share function over normalized liquid resources
377 and normalized iliquid resources when the agent is able to adjust his
378 portfolio.
379 dvdmFunc_Adj : MargValueFunc2D
380 Marginal value function over normalized liquid resources when the agent
381 is able to adjust his portfolio.
382 dvdnFunc_Adj : MargValueFunc2D
383 Marginal value function over normalized iliquid resources when the
384 agent is able to adjust his portfolio.
385 vFunc_Fxd : ValueFunc3D
386 Stage value function over normalized liquid resources, normalized
387 iliquid resources, and income contribution share when the agent is not
388 able to adjust his portfolio.
389 ShareFunc_Fxd : Interp3D
390 Income contribution share function over normalized liquid resources,
391 iliquid resources, and income contribution share when the agent is not
392 able to adjust his portfolio.
393 Should be an IdentityFunc.
394 dvdmFunc_Fxd : MargValueFunc3D
395 Marginal value function over normalized liquid resources when the agent
396 is not able to adjust his portfolio.
397 dvdnFunc_Fxd : MargValueFunc3D
398 Marginal value function over normalized iliquid resources when the
399 agent is not able to adjust his portfolio.
400 dvdsFunc_Fxd : Interp3D
401 Marginal value function over income contribution share when the agent
402 is not able to adjust his portfolio
403 """
405 distance_criteria = ["dvdmFunc_Adj", "dvdnFunc_Adj"]
407 def __init__(
408 self,
409 # Contribution stage, adjust
410 vFunc_Adj=None,
411 ShareFunc_Adj=None,
412 dvdmFunc_Adj=None,
413 dvdnFunc_Adj=None,
414 # Contribution stage, fixed
415 vFunc_Fxd=None,
416 ShareFunc_Fxd=None,
417 dvdmFunc_Fxd=None,
418 dvdnFunc_Fxd=None,
419 dvdsFunc_Fxd=None,
420 ):
421 # Contribution stage, adjust
422 if vFunc_Adj is None:
423 vFunc_Adj = NullFunc()
424 if ShareFunc_Adj is None:
425 ShareFunc_Adj = NullFunc()
426 if dvdmFunc_Adj is None:
427 dvdmFunc_Adj = NullFunc()
428 if dvdnFunc_Adj is None:
429 dvdnFunc_Adj = NullFunc()
431 # Contribution stage, fixed
432 if vFunc_Fxd is None:
433 vFunc_Fxd = NullFunc()
434 if ShareFunc_Fxd is None:
435 ShareFunc_Fxd = NullFunc()
436 if dvdmFunc_Fxd is None:
437 dvdmFunc_Fxd = NullFunc()
438 if dvdnFunc_Fxd is None:
439 dvdnFunc_Fxd = NullFunc()
440 if dvdsFunc_Fxd is None:
441 dvdsFunc_Fxd = NullFunc()
443 # Set attributes of self
444 self.vFunc_Adj = vFunc_Adj
445 self.ShareFunc_Adj = ShareFunc_Adj
446 self.dvdmFunc_Adj = dvdmFunc_Adj
447 self.dvdnFunc_Adj = dvdnFunc_Adj
449 self.vFunc_Fxd = vFunc_Fxd
450 self.ShareFunc_Fxd = ShareFunc_Fxd
451 self.dvdmFunc_Fxd = dvdmFunc_Fxd
452 self.dvdnFunc_Fxd = dvdnFunc_Fxd
453 self.dvdsFunc_Fxd = dvdsFunc_Fxd
456# Class for the consumption stage solution
457class RiskyContribCnsSolution(MetricObject):
458 """
459 A class for representing the solution to the consumption stage of the
460 'RiskyContrib' model.
462 Parameters
463 ----------
464 vFunc : ValueFunc3D
465 Stage-value function over normalized liquid resources, normalized
466 iliquid resources, and income contribution share.
467 cFunc : Interp3D
468 Consumption function over normalized liquid resources, normalized
469 iliquid resources, and income contribution share.
470 dvdmFunc : MargValueFunc3D
471 Marginal value function over normalized liquid resources.
472 dvdnFunc : MargValueFunc3D
473 Marginal value function over normalized iliquid resources.
474 dvdsFunc : Interp3D
475 Marginal value function over income contribution share.
476 """
478 distance_criteria = ["dvdmFunc", "dvdnFunc"]
480 def __init__(
481 self,
482 # Consumption stage
483 vFunc=None,
484 cFunc=None,
485 dvdmFunc=None,
486 dvdnFunc=None,
487 dvdsFunc=None,
488 ):
489 if vFunc is None:
490 vFunc = NullFunc()
491 if cFunc is None:
492 cFunc = NullFunc()
493 if dvdmFunc is None:
494 dvdmFunc = NullFunc()
495 if dvdnFunc is None:
496 dvdmFunc = NullFunc()
497 if dvdsFunc is None:
498 dvdsFunc = NullFunc()
500 self.vFunc = vFunc
501 self.cFunc = cFunc
502 self.dvdmFunc = dvdmFunc
503 self.dvdnFunc = dvdnFunc
504 self.dvdsFunc = dvdsFunc
507# Class for the solution of a whole period
508class RiskyContribSolution(MetricObject):
509 """
510 A class for representing the solution to a full time-period of the
511 'RiskyContrib' agent type's problem.
513 Parameters
514 ----------
515 Reb : RiskyContribRebSolution
516 Solution to the period's rebalancing stage.
517 Sha : RiskyContribShaSolution
518 Solution to the period's contribution-share stage.
519 Cns : RiskyContribCnsSolution
520 Solution to the period's consumption stage.
521 """
523 # Solutions are to be compared on the basis of their sub-period solutions
524 distance_criteria = ["stage_sols"]
526 def __init__(self, Reb, Sha, Cns):
527 # Dictionary of stage solutions
528 self.stage_sols = {"Reb": Reb, "Sha": Sha, "Cns": Cns}
531# %% Auxiliary functions and transition equations for the RiskyContrib model.
534def rebalance_assets(d, m, n, tau):
535 """
536 A function that produces post-rebalancing assets for given initial assets,
537 rebalancing action, and tax rate.
539 Parameters
540 ----------
541 d : np.array
542 Array with rebalancing decisions. d > 0 represents depositing d*m into
543 the risky asset account. d<0 represents withdrawing ``|d|*n`` (pre-tax)
544 from the risky account into the risky account.
545 m : np.array
546 Initial risk-free assets.
547 n : np.array
548 Initial risky assets.
549 tau : float
550 Tax rate on flows from the risky to the risk-free asset.
552 Returns
553 -------
554 mTil : np.array
555 Post-rebalancing risk-free assets.
556 nTil : np.arrat
557 Post-rebalancing risky assets.
559 """
560 # Initialize
561 mTil = np.zeros_like(m) + np.nan
562 nTil = np.zeros_like(m) + np.nan
564 # Contributions
565 inds = d >= 0
566 mTil[inds] = m[inds] * (1 - d[inds])
567 nTil[inds] = n[inds] + m[inds] * d[inds]
569 # Withdrawals
570 inds = d < 0
571 mTil[inds] = m[inds] - d[inds] * n[inds] * (1 - tau)
572 nTil[inds] = n[inds] * (1 + d[inds])
574 return (mTil, nTil)
577# Transition equations for the consumption stage
578def m_nrm_next(shocks, aNrm, Share, Rfree, PermGroFac):
579 """
580 Given end-of-period balances and contribution share and the
581 start-of-next-period shocks, figure out next period's normalized riskless
582 assets
584 Parameters
585 ----------
586 shocks : np.array
587 Length-3 array with the stochastic shocks that get realized between the
588 end of the current period and the start of next period. Their order is
589 (0) permanent income shock, (1) transitory income shock, (2) risky
590 asset return.
591 aNrm : float
592 End-of-period risk-free asset balances.
593 Share : float
594 End-of-period income deduction share.
595 Rfree : float
596 Risk-free return factor.
597 PermGroFac : float
598 Permanent income growth factor.
600 Returns
601 -------
602 m_nrm_tp1 : float
603 Next-period normalized riskless balance.
605 """
606 # Extract shocks
607 perm_shk = shocks[0]
608 tran_shk = shocks[1]
610 m_nrm_tp1 = Rfree * aNrm / (perm_shk * PermGroFac) + (1.0 - Share) * tran_shk
612 return m_nrm_tp1
615def n_nrm_next(shocks, nNrm, Share, PermGroFac):
616 """
617 Given end-of-period balances and contribution share and the
618 start-of-next-period shocks, figure out next period's normalized risky
619 assets
621 Parameters
622 ----------
623 shocks : np.array
624 Length-3 array with the stochastic shocks that get realized between the
625 end of the current period and the start of next period. Their order is
626 (0) permanent income shock, (1) transitory income shock, (2) risky
627 asset return.
628 nNrm : float
629 End-of-period risky asset balances.
630 Share : float
631 End-of-period income deduction share.
632 PermGroFac : float
633 Permanent income growth factor.
635 Returns
636 -------
637 n_nrm_tp1 : float
638 Next-period normalized risky balance.
640 """
642 # Extract shocks
643 perm_shk = shocks[0]
644 tran_shk = shocks[1]
645 R_risky = shocks[2]
647 n_nrm_tp1 = R_risky * nNrm / (perm_shk * PermGroFac) + Share * tran_shk
649 return n_nrm_tp1
652# %% RiskyContrib solvers
655# Consumption stage solver
656def solve_RiskyContrib_Cns(
657 solution_next,
658 ShockDstn,
659 IncShkDstn,
660 RiskyDstn,
661 IndepDstnBool,
662 LivPrb,
663 DiscFac,
664 CRRA,
665 Rfree,
666 PermGroFac,
667 BoroCnstArt,
668 aXtraGrid,
669 nNrmGrid,
670 mNrmGrid,
671 ShareGrid,
672 vFuncBool,
673 AdjustPrb,
674 DiscreteShareBool,
675 joint_dist_solver,
676 **unused_params,
677):
678 """
679 Solves the consumption stage of the agent's problem
681 Parameters
682 ----------
683 solution_next : RiskyContribRebSolution
684 Solution to the first stage of the next period in the agent's problem.
685 ShockDstn : DiscreteDistribution
686 Joint distribution of next period's (0) permanent income shock, (1)
687 transitory income shock, and (2) risky asset return factor.
688 IncShkDstn : DiscreteDistribution
689 Joint distribution of next period's (0) permanent income shock and (1)
690 transitory income shock.
691 RiskyDstn : DiscreteDistribution
692 Distribution of next period's risky asset return factor.
693 IndepDstnBool : bool
694 Indicates whether the income and risky return distributions are
695 independent.
696 LivPrb : float
697 Probability of surviving until next period.
698 DiscFac : float
699 Time-preference discount factor.
700 CRRA : float
701 Coefficient of relative risk aversion.
702 Rfree : float
703 Risk-free return factor.
704 PermGroFac : float
705 Deterministic permanent income growth factor.
706 BoroCnstArt : float
707 Minimum allowed market resources (must be 0).
708 aXtraGrid : numpy array
709 Exogenous grid for end-of-period risk free resources.
710 nNrmGrid : numpy array
711 Exogenous grid for risky resources.
712 mNrmGrid : numpy array
713 Exogenous grid for risk-free resources.
714 ShareGrid : numpt array
715 Exogenous grid for the income contribution share.
716 vFuncBool : bool
717 Boolean that determines wether the value function's level needs to be
718 computed.
719 AdjustPrb : float
720 Probability thet the agent will be able to adjust his portfolio next
721 period.
722 DiscreteShareBool : bool
723 Boolean that determines whether only a discrete set of contribution
724 shares (ShareGrid) is allowed.
725 joint_dist_solver: bool
726 Should the general solver be used even if income and returns are
727 independent?
729 Returns
730 -------
731 solution : RiskyContribCnsSolution
732 Solution to the agent's consumption stage problem.
734 """
735 # Make sure the individual is liquidity constrained. Allowing a consumer to
736 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
737 if BoroCnstArt != 0.0:
738 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
740 # Make sure that if risky portfolio share is optimized only discretely, then
741 # the value function is also constructed (else this task would be impossible).
742 if DiscreteShareBool and (not vFuncBool):
743 raise ValueError(
744 "PortfolioConsumerType requires vFuncBool to be True when DiscreteShareBool is True!"
745 )
747 # Define temporary functions for utility and its derivative and inverse
748 u = lambda x: utility(x, CRRA)
749 uPinv = lambda x: utilityP_inv(x, CRRA)
750 uInv = lambda x: utility_inv(x, CRRA)
752 # Unpack next period's solution
753 vFunc_Reb_Adj_next = solution_next.vFunc_Adj
754 dvdmFunc_Reb_Adj_next = solution_next.dvdmFunc_Adj
755 dvdnFunc_Reb_Adj_next = solution_next.dvdnFunc_Adj
757 vFunc_Reb_Fxd_next = solution_next.vFunc_Fxd
758 dvdmFunc_Reb_Fxd_next = solution_next.dvdmFunc_Fxd
759 dvdnFunc_Reb_Fxd_next = solution_next.dvdnFunc_Fxd
760 dvdsFunc_Reb_Fxd_next = solution_next.dvdsFunc_Fxd
762 # STEP ONE
763 # Find end-of-period (continuation) value function and its derivatives.
765 # Start by constructing functions for next-period's pre-adjustment-shock
766 # expected value functions
767 if AdjustPrb < 1.0:
768 dvdm_next = lambda m, n, s: AdjustPrb * dvdmFunc_Reb_Adj_next(m, n) + (
769 1.0 - AdjustPrb
770 ) * dvdmFunc_Reb_Fxd_next(m, n, s)
771 dvdn_next = lambda m, n, s: AdjustPrb * dvdnFunc_Reb_Adj_next(m, n) + (
772 1.0 - AdjustPrb
773 ) * dvdnFunc_Reb_Fxd_next(m, n, s)
774 dvds_next = lambda m, n, s: (1.0 - AdjustPrb) * dvdsFunc_Reb_Fxd_next(m, n, s)
776 # Value function if needed
777 if vFuncBool:
778 v_next = lambda m, n, s: AdjustPrb * vFunc_Reb_Adj_next(m, n) + (
779 1.0 - AdjustPrb
780 ) * vFunc_Reb_Fxd_next(m, n, s)
782 else:
783 dvdm_next = lambda m, n, s: dvdmFunc_Reb_Adj_next(m, n)
784 dvdn_next = lambda m, n, s: dvdnFunc_Reb_Adj_next(m, n)
785 dvds_next = ConstantFunction(0.0)
787 if vFuncBool:
788 v_next = lambda m, n, s: vFunc_Reb_Adj_next(m, n)
790 if IndepDstnBool and not joint_dist_solver:
791 # If income and returns are independent we can use the law of iterated
792 # expectations to speed up the computation of end-of-period derivatives
794 # Define "post-return variables"
795 # b_aux = aNrm * R
796 # g_aux = nNrmTilde * Rtilde
797 # and create a function that interpolates end-of-period marginal values
798 # as functions of those and the contribution share
800 def post_return_derivs(inc_shocks, b_aux, g_aux, s):
801 perm_shk = inc_shocks[0]
802 tran_shk = inc_shocks[1]
804 temp_fac_A = utilityP(perm_shk * PermGroFac, CRRA)
805 temp_fac_B = (perm_shk * PermGroFac) ** (1.0 - CRRA)
807 # Find next-period asset balances
808 m_next = b_aux / (perm_shk * PermGroFac) + (1.0 - s) * tran_shk
809 n_next = g_aux / (perm_shk * PermGroFac) + s * tran_shk
811 # Interpolate next-period-value derivatives
812 dvdm_tp1 = dvdm_next(m_next, n_next, s)
813 dvdn_tp1 = dvdn_next(m_next, n_next, s)
814 if tran_shk == 0:
815 dvds_tp1 = dvds_next(m_next, n_next, s)
816 else:
817 dvds_tp1 = tran_shk * (dvdn_tp1 - dvdm_tp1) + dvds_next(
818 m_next, n_next, s
819 )
821 # Discount next-period-value derivatives to current period
823 # Liquid resources
824 pr_dvda = temp_fac_A * dvdm_tp1
825 # Iliquid resources
826 pr_dvdn = temp_fac_A * dvdn_tp1
827 # Contribution share
828 pr_dvds = temp_fac_B * dvds_tp1
830 # End of period value function, if needed
831 if vFuncBool:
832 pr_v = temp_fac_B * v_next(m_next, n_next, s)
833 return np.stack([pr_dvda, pr_dvdn, pr_dvds, pr_v])
834 else:
835 return np.stack([pr_dvda, pr_dvdn, pr_dvds])
837 # Define grids
838 b_aux_grid = np.concatenate([np.array([0.0]), Rfree * aXtraGrid])
839 g_aux_grid = np.concatenate(
840 [np.array([0.0]), max(RiskyDstn.atoms.flatten()) * nNrmGrid]
841 )
843 # Create tiled arrays with conforming dimensions.
844 b_aux_tiled, g_aux_tiled, Share_tiled = np.meshgrid(
845 b_aux_grid, g_aux_grid, ShareGrid, indexing="ij"
846 )
848 # Find end of period derivatives and value as expectations of (discounted)
849 # next period's derivatives and value.
850 pr_derivs = calc_expectation(
851 IncShkDstn, post_return_derivs, b_aux_tiled, g_aux_tiled, Share_tiled
852 )
854 # Unpack results and create interpolators
855 pr_dvdb_func = MargValueFuncCRRA(
856 TrilinearInterp(uPinv(pr_derivs[0]), b_aux_grid, g_aux_grid, ShareGrid),
857 CRRA,
858 )
859 pr_dvdg_func = MargValueFuncCRRA(
860 TrilinearInterp(uPinv(pr_derivs[1]), b_aux_grid, g_aux_grid, ShareGrid),
861 CRRA,
862 )
863 pr_dvds_func = TrilinearInterp(pr_derivs[2], b_aux_grid, g_aux_grid, ShareGrid)
865 if vFuncBool:
866 pr_vFunc = ValueFuncCRRA(
867 TrilinearInterp(uInv(pr_derivs[3]), b_aux_grid, g_aux_grid, ShareGrid),
868 CRRA,
869 )
871 # Now construct a function that produces end-of-period derivatives
872 # given the risky return draw
873 def end_of_period_derivs(risky_ret, a, nTil, s):
874 """
875 Computes the end-of-period derivatives (and optionally the value) of the
876 continuation value function, conditional on risky returns. This is so that the
877 expectations can be calculated by integrating over risky returns.
879 Parameters
880 ----------
881 risky_ret : float
882 Risky return factor
883 a : float
884 end-of-period risk-free assets.
885 nTil : float
886 end-of-period risky assets.
887 s : float
888 end-of-period income deduction share.
889 """
891 # Find next-period asset balances
892 b_aux = a * Rfree
893 g_aux = nTil * risky_ret
895 # Interpolate post-return derivatives
896 pr_dvdb = pr_dvdb_func(b_aux, g_aux, s)
897 pr_dvdg = pr_dvdg_func(b_aux, g_aux, s)
898 pr_dvds = pr_dvds_func(b_aux, g_aux, s)
900 # Discount
902 # Liquid resources
903 end_of_prd_dvda = DiscFac * Rfree * LivPrb * pr_dvdb
904 # Iliquid resources
905 end_of_prd_dvdn = DiscFac * risky_ret * LivPrb * pr_dvdg
906 # Contribution share
907 end_of_prd_dvds = DiscFac * LivPrb * pr_dvds
909 # End of period value function, i11f needed
910 if vFuncBool:
911 end_of_prd_v = DiscFac * LivPrb * pr_vFunc(b_aux, g_aux, s)
912 return np.stack(
913 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
914 )
915 else:
916 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
918 else:
919 # If income and returns are not independent, we just integrate over
920 # them jointly.
922 # Construct a function that evaluates and discounts them given a
923 # vector of return and income shocks and an end-of-period state
924 def end_of_period_derivs(shocks, a, nTil, s):
925 """
926 Computes the end-of-period derivatives (and optionally the value) of the
927 continuation value function, conditional on shocks. This is so that the
928 expectations can be calculated by integrating over shocks.
930 Parameters
931 ----------
932 shocks : np.array
933 Length-3 array with the stochastic shocks that get realized between the
934 end of the current period and the start of next period. Their order is
935 (0) permanent income shock, (1) transitory income shock, (2) risky
936 asset return.
937 a : float
938 end-of-period risk-free assets.
939 nTil : float
940 end-of-period risky assets.
941 s : float
942 end-of-period income deduction share.
943 """
944 temp_fac_A = utilityP(shocks[0] * PermGroFac, CRRA)
945 temp_fac_B = (shocks[0] * PermGroFac) ** (1.0 - CRRA)
947 # Find next-period asset balances
948 m_next = m_nrm_next(shocks, a, s, Rfree, PermGroFac)
949 n_next = n_nrm_next(shocks, nTil, s, PermGroFac)
951 # Interpolate next-period-value derivatives
952 dvdm_tp1 = dvdm_next(m_next, n_next, s)
953 dvdn_tp1 = dvdn_next(m_next, n_next, s)
954 if shocks[1] == 0:
955 dvds_tp1 = dvds_next(m_next, n_next, s)
956 else:
957 dvds_tp1 = shocks[1] * (dvdn_tp1 - dvdm_tp1) + dvds_next(
958 m_next, n_next, s
959 )
961 # Discount next-period-value derivatives to current period
963 # Liquid resources
964 end_of_prd_dvda = DiscFac * Rfree * LivPrb * temp_fac_A * dvdm_tp1
965 # Iliquid resources
966 end_of_prd_dvdn = DiscFac * shocks[2] * LivPrb * temp_fac_A * dvdn_tp1
967 # Contribution share
968 end_of_prd_dvds = DiscFac * LivPrb * temp_fac_B * dvds_tp1
970 # End of period value function, i11f needed
971 if vFuncBool:
972 end_of_prd_v = DiscFac * LivPrb * temp_fac_B * v_next(m_next, n_next, s)
973 return np.stack(
974 [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v]
975 )
976 else:
977 return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds])
979 # Now find the expected values on a (a, nTil, s) grid
981 # The "inversion" machinery can deal with assets of 0 even if there is a
982 # natural borrowing constraint, so include zeros.
983 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
984 aNrmGrid = np.concatenate([np.array([0.0]), aXtraGrid])
986 # Create tiled arrays with conforming dimensions.
987 aNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
988 aNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
989 )
991 # Find end of period derivatives and value as expectations of (discounted)
992 # next period's derivatives and value.
993 eop_derivs = calc_expectation(
994 RiskyDstn if IndepDstnBool and not joint_dist_solver else ShockDstn,
995 end_of_period_derivs,
996 aNrm_tiled,
997 nNrm_tiled,
998 Share_tiled,
999 )
1001 # Unpack results
1002 eop_dvdaNvrs = uPinv(eop_derivs[0])
1003 eop_dvdnNvrs = uPinv(eop_derivs[1])
1004 eop_dvds = eop_derivs[2]
1005 if vFuncBool:
1006 eop_vNvrs = uInv(eop_derivs[3])
1008 # Construct an interpolator for eop_V. It will be used later.
1009 eop_vFunc = ValueFuncCRRA(
1010 TrilinearInterp(eop_vNvrs, aNrmGrid, nNrmGrid, ShareGrid), CRRA
1011 )
1013 # STEP TWO:
1014 # Solve the consumption problem and create interpolators for c, vCns,
1015 # and its derivatives.
1017 # Apply EGM over liquid resources at every (n,s) to find consumption.
1018 c_end = eop_dvdaNvrs
1019 mNrm_end = aNrm_tiled + c_end
1021 # Now construct interpolators for c and the derivatives of vCns.
1022 # The m grid is different for every (n,s). We interpolate the object of
1023 # interest on the regular m grid for every (n,s). At the end we will have
1024 # values of the functions of interest on a regular (m,n,s) grid. We use
1025 # trilinear interpolation on those points.
1027 # Expand the exogenous m grid to contain 0.
1028 mNrmGrid = np.insert(mNrmGrid, 0, 0)
1030 # Dimensions might have changed, so re-create tiled arrays
1031 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1032 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1033 )
1035 # Initialize arrays
1036 c_vals = np.zeros_like(mNrm_tiled)
1037 dvdnNvrs_vals = np.zeros_like(mNrm_tiled)
1038 dvds_vals = np.zeros_like(mNrm_tiled)
1040 nNrm_N = nNrmGrid.size
1041 Share_N = ShareGrid.size
1042 for nInd in range(nNrm_N):
1043 for sInd in range(Share_N):
1044 # Extract the endogenous m grid for particular (n,s).
1045 m_ns = mNrm_end[:, nInd, sInd]
1047 # Check if there is a natural constraint
1048 if m_ns[0] == 0.0:
1049 # There's no need to insert points since we have m==0.0
1051 # c
1052 c_vals[:, nInd, sInd] = LinearInterp(m_ns, c_end[:, nInd, sInd])(
1053 mNrmGrid
1054 )
1056 # dvdnNvrs
1057 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1058 m_ns, eop_dvdnNvrs[:, nInd, sInd]
1059 )(mNrmGrid)
1061 # dvds
1062 dvds_vals[:, nInd, sInd] = LinearInterp(m_ns, eop_dvds[:, nInd, sInd])(
1063 mNrmGrid
1064 )
1066 else:
1067 # We know that:
1068 # -The lowest gridpoints of both a and n are 0.
1069 # -Consumption at m < m0 is m.
1070 # -dvdn_Fxd at (m,n) for m < m0(n) is dvdn_Fxd(m0,n)
1071 # -Same is true for dvds_Fxd
1073 m_ns = np.concatenate([np.array([0]), m_ns])
1075 # c
1076 c_vals[:, nInd, sInd] = LinearInterp(
1077 m_ns, np.concatenate([np.array([0]), c_end[:, nInd, sInd]])
1078 )(mNrmGrid)
1080 # dvdnNvrs
1081 dvdnNvrs_vals[:, nInd, sInd] = LinearInterp(
1082 m_ns,
1083 np.concatenate(
1084 [
1085 np.array([eop_dvdnNvrs[0, nInd, sInd]]),
1086 eop_dvdnNvrs[:, nInd, sInd],
1087 ]
1088 ),
1089 )(mNrmGrid)
1091 # dvds
1092 dvds_vals[:, nInd, sInd] = LinearInterp(
1093 m_ns,
1094 np.concatenate(
1095 [
1096 np.array([eop_dvds[0, nInd, sInd]]),
1097 eop_dvds[:, nInd, sInd],
1098 ]
1099 ),
1100 )(mNrmGrid)
1102 # With the arrays filled, create 3D interpolators
1104 # Consumption interpolator
1105 cFunc = TrilinearInterp(c_vals, mNrmGrid, nNrmGrid, ShareGrid)
1106 # dvdmCns interpolator
1107 dvdmFunc_Cns = MargValueFuncCRRA(cFunc, CRRA)
1108 # dvdnCns interpolator
1109 dvdnNvrsFunc = TrilinearInterp(dvdnNvrs_vals, mNrmGrid, nNrmGrid, ShareGrid)
1110 dvdnFunc_Cns = MargValueFuncCRRA(dvdnNvrsFunc, CRRA)
1111 # dvdsCns interpolator
1112 dvdsFunc_Cns = TrilinearInterp(dvds_vals, mNrmGrid, nNrmGrid, ShareGrid)
1114 # Compute value function if needed
1115 if vFuncBool:
1116 # Consumption in the regular grid
1117 aNrm_reg = mNrm_tiled - c_vals
1118 vCns = u(c_vals) + eop_vFunc(aNrm_reg, nNrm_tiled, Share_tiled)
1119 vNvrsCns = uInv(vCns)
1120 vNvrsFunc_Cns = TrilinearInterp(vNvrsCns, mNrmGrid, nNrmGrid, ShareGrid)
1121 vFunc_Cns = ValueFuncCRRA(vNvrsFunc_Cns, CRRA)
1122 else:
1123 vFunc_Cns = NullFunc()
1125 # Assemble solution
1126 solution = RiskyContribCnsSolution(
1127 vFunc=vFunc_Cns,
1128 cFunc=cFunc,
1129 dvdmFunc=dvdmFunc_Cns,
1130 dvdnFunc=dvdnFunc_Cns,
1131 dvdsFunc=dvdsFunc_Cns,
1132 )
1134 return solution
1137# Solver for the contribution stage
1138def solve_RiskyContrib_Sha(
1139 solution_next,
1140 CRRA,
1141 AdjustPrb,
1142 mNrmGrid,
1143 nNrmGrid,
1144 ShareGrid,
1145 DiscreteShareBool,
1146 vFuncBool,
1147 **unused_params,
1148):
1149 """
1150 Solves the income-contribution-share stag of the agent's problem
1152 Parameters
1153 ----------
1154 solution_next : RiskyContribCnsSolution
1155 Solution to the agent's consumption stage problem that follows.
1156 CRRA : float
1157 Coefficient of relative risk aversion.
1158 AdjustPrb : float
1159 Probability that the agent will be able to rebalance his portfolio
1160 next period.
1161 mNrmGrid : numpy array
1162 Exogenous grid for risk-free resources.
1163 nNrmGrid : numpy array
1164 Exogenous grid for risky resources.
1165 ShareGrid : numpy array
1166 Exogenous grid for the income contribution share.
1167 DiscreteShareBool : bool
1168 Boolean that determines whether only a discrete set of contribution
1169 shares (ShareGrid) is allowed.
1170 vFuncBool : bool
1171 Determines whether the level of the value function is computed.
1173 Yields
1174 ------
1175 solution : RiskyContribShaSolution
1176 Solution to the income-contribution-share stage of the agent's problem.
1178 """
1179 # Unpack solution from the next sub-stage
1180 vFunc_Cns_next = solution_next.vFunc
1181 cFunc_next = solution_next.cFunc
1182 dvdmFunc_Cns_next = solution_next.dvdmFunc
1183 dvdnFunc_Cns_next = solution_next.dvdnFunc
1184 dvdsFunc_Cns_next = solution_next.dvdsFunc
1186 uPinv = lambda x: utilityP_inv(x, CRRA)
1188 # Create tiled grids
1190 # Add 0 to the m and n grids
1191 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1192 nNrm_N = len(nNrmGrid)
1193 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1194 mNrm_N = len(mNrmGrid)
1196 if AdjustPrb == 1.0:
1197 # If the readjustment probability is 1, set the share to 0:
1198 # - If there is a withdrawal tax: better for the agent to observe
1199 # income before rebalancing.
1200 # - If there is no tax: all shares should yield the same value.
1201 mNrm_tiled, nNrm_tiled = np.meshgrid(mNrmGrid, nNrmGrid, indexing="ij")
1203 opt_idx = np.zeros_like(mNrm_tiled, dtype=int)
1204 opt_Share = ShareGrid[opt_idx]
1206 if vFuncBool:
1207 vNvrsSha = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, opt_Share)
1209 else:
1210 # Figure out optimal share by evaluating all alternatives at all
1211 # (m,n) combinations
1212 m_idx_tiled, n_idx_tiled = np.meshgrid(
1213 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1214 )
1216 mNrm_tiled, nNrm_tiled, Share_tiled = np.meshgrid(
1217 mNrmGrid, nNrmGrid, ShareGrid, indexing="ij"
1218 )
1220 if DiscreteShareBool:
1221 # Evaluate value function to optimize over shares.
1222 # Do it in inverse space
1223 vNvrs = vFunc_Cns_next.vFuncNvrs(mNrm_tiled, nNrm_tiled, Share_tiled)
1225 # Find the optimal share at each (m,n).
1226 opt_idx = np.argmax(vNvrs, axis=2)
1228 # Compute objects needed for the value function and its derivatives
1229 vNvrsSha = vNvrs[m_idx_tiled, n_idx_tiled, opt_idx]
1230 opt_Share = ShareGrid[opt_idx]
1232 # Project grids
1233 mNrm_tiled = mNrm_tiled[:, :, 0]
1234 nNrm_tiled = nNrm_tiled[:, :, 0]
1236 else:
1237 # Evaluate the marginal value of the contribution share at
1238 # every (m,n,s) gridpoint
1239 dvds = dvdsFunc_Cns_next(mNrm_tiled, nNrm_tiled, Share_tiled)
1241 # If the derivative is negative at the lowest share, then s[0] is optimal
1242 constrained_bot = dvds[:, :, 0] <= 0.0
1243 # If it is poitive at the highest share, then s[-1] is optimal
1244 constrained_top = dvds[:, :, -1] >= 0.0
1246 # Find indices at which the derivative crosses 0 for the 1st time
1247 # will be 0 if it never does, but "constrained_top/bot" deals with that
1248 crossings = np.logical_and(dvds[:, :, :-1] >= 0.0, dvds[:, :, 1:] <= 0.0)
1249 idx = np.argmax(crossings, axis=2)
1251 # Linearly interpolate the optimal share
1252 idx1 = idx + 1
1253 slopes = (
1254 dvds[m_idx_tiled, n_idx_tiled, idx1]
1255 - dvds[m_idx_tiled, n_idx_tiled, idx]
1256 ) / (ShareGrid[idx1] - ShareGrid[idx])
1257 opt_Share = ShareGrid[idx] - dvds[m_idx_tiled, n_idx_tiled, idx] / slopes
1259 # Replace the ones we knew were constrained
1260 opt_Share[constrained_bot] = ShareGrid[0]
1261 opt_Share[constrained_top] = ShareGrid[-1]
1263 # Project grids
1264 mNrm_tiled = mNrm_tiled[:, :, 0]
1265 nNrm_tiled = nNrm_tiled[:, :, 0]
1267 # Evaluate the inverse value function at the optimal shares
1268 if vFuncBool:
1269 vNvrsSha = vFunc_Cns_next.func(mNrm_tiled, nNrm_tiled, opt_Share)
1271 dvdmNvrsSha = cFunc_next(mNrm_tiled, nNrm_tiled, opt_Share)
1272 dvdnSha = dvdnFunc_Cns_next(mNrm_tiled, nNrm_tiled, opt_Share)
1273 dvdnNvrsSha = uPinv(dvdnSha)
1275 # Interpolators
1277 # Value function if needed
1278 if vFuncBool:
1279 vNvrsFunc_Sha = BilinearInterp(vNvrsSha, mNrmGrid, nNrmGrid)
1280 vFunc_Sha = ValueFuncCRRA(vNvrsFunc_Sha, CRRA)
1281 else:
1282 vFunc_Sha = NullFunc()
1284 # Contribution share function
1285 if DiscreteShareBool:
1286 ShareFunc = DiscreteInterp(
1287 BilinearInterp(opt_idx, mNrmGrid, nNrmGrid), ShareGrid
1288 )
1289 else:
1290 ShareFunc = BilinearInterp(opt_Share, mNrmGrid, nNrmGrid)
1292 # Derivatives
1293 dvdmNvrsFunc_Sha = BilinearInterp(dvdmNvrsSha, mNrmGrid, nNrmGrid)
1294 dvdmFunc_Sha = MargValueFuncCRRA(dvdmNvrsFunc_Sha, CRRA)
1295 dvdnNvrsFunc_Sha = BilinearInterp(dvdnNvrsSha, mNrmGrid, nNrmGrid)
1296 dvdnFunc_Sha = MargValueFuncCRRA(dvdnNvrsFunc_Sha, CRRA)
1298 solution = RiskyContribShaSolution(
1299 vFunc_Adj=vFunc_Sha,
1300 ShareFunc_Adj=ShareFunc,
1301 dvdmFunc_Adj=dvdmFunc_Sha,
1302 dvdnFunc_Adj=dvdnFunc_Sha,
1303 # The fixed agent does nothing at this stage,
1304 # so his value functions are the next problem's
1305 vFunc_Fxd=vFunc_Cns_next,
1306 ShareFunc_Fxd=IdentityFunction(i_dim=2, n_dims=3),
1307 dvdmFunc_Fxd=dvdmFunc_Cns_next,
1308 dvdnFunc_Fxd=dvdnFunc_Cns_next,
1309 dvdsFunc_Fxd=dvdsFunc_Cns_next,
1310 )
1312 return solution
1315# Solver for the asset rebalancing stage
1316def solve_RiskyContrib_Reb(
1317 solution_next,
1318 CRRA,
1319 WithdrawTax,
1320 nNrmGrid,
1321 mNrmGrid,
1322 dfracGrid,
1323 vFuncBool,
1324 **unused_params,
1325):
1326 """
1327 Solves the asset-rebalancing-stage of the agent's problem
1329 Parameters
1330 ----------
1331 solution_next : RiskyContribShaSolution
1332 Solution to the income-contribution-share stage problem that follows.
1333 CRRA : float
1334 Coefficient of relative risk aversion.
1335 WithdrawTax : float
1336 Tax rate on risky asset withdrawals.
1337 nNrmGrid : numpy array
1338 Exogenous grid for risky resources.
1339 mNrmGrid : numpy array
1340 Exogenous grid for risk-free resources.
1341 dfracGrid : numpy array
1342 Grid for rebalancing flows. The final grid will be equivalent to
1343 [-nNrm*dfracGrid, dfracGrid*mNrm].
1344 vFuncBool : bool
1345 Determines whether the level of th value function must be computed.
1347 Returns
1348 -------
1349 solution : RiskyContribShaSolution
1350 Solution to the asset-rebalancing stage of the agent's problem.
1352 """
1353 # Extract next stage's solution
1354 vFunc_Adj_next = solution_next.vFunc_Adj
1355 dvdmFunc_Adj_next = solution_next.dvdmFunc_Adj
1356 dvdnFunc_Adj_next = solution_next.dvdnFunc_Adj
1358 vFunc_Fxd_next = solution_next.vFunc_Fxd
1359 dvdmFunc_Fxd_next = solution_next.dvdmFunc_Fxd
1360 dvdnFunc_Fxd_next = solution_next.dvdnFunc_Fxd
1361 dvdsFunc_Fxd_next = solution_next.dvdsFunc_Fxd
1363 uPinv = lambda x: utilityP_inv(x, CRRA)
1365 # Create tiled grids
1367 # Add 0 to the m and n grids
1368 nNrmGrid = np.concatenate([np.array([0.0]), nNrmGrid])
1369 nNrm_N = len(nNrmGrid)
1370 mNrmGrid = np.concatenate([np.array([0.0]), mNrmGrid])
1371 mNrm_N = len(mNrmGrid)
1372 d_N = len(dfracGrid)
1374 # Duplicate d so that possible values are -dfracGrid,dfracGrid. Duplicate 0 is
1375 # intentional since the tax causes a discontinuity. We need the value
1376 # from the left and right.
1377 dfracGrid = np.concatenate((-1 * np.flip(dfracGrid), dfracGrid))
1379 # It will be useful to pre-evaluate marginals at every (m,n,d) combination
1381 # Create tiled arrays for every d,m,n option
1382 d_N2 = len(dfracGrid)
1383 d_tiled, mNrm_tiled, nNrm_tiled = np.meshgrid(
1384 dfracGrid, mNrmGrid, nNrmGrid, indexing="ij"
1385 )
1387 # Get post-rebalancing assets.
1388 m_tilde, n_tilde = rebalance_assets(d_tiled, mNrm_tiled, nNrm_tiled, WithdrawTax)
1390 # Now the marginals, in inverse space
1391 dvdmNvrs = dvdmFunc_Adj_next.cFunc(m_tilde, n_tilde)
1392 dvdnNvrs = dvdnFunc_Adj_next.cFunc(m_tilde, n_tilde)
1394 # Pre-evaluate the inverse of (1-WithdrawTax)
1395 taxNvrs = uPinv(1 - WithdrawTax)
1396 # Create a tiled array of the tax
1397 taxNvrs_tiled = np.tile(
1398 np.reshape(
1399 np.concatenate([np.repeat(taxNvrs, d_N), np.ones(d_N, dtype=np.double)]),
1400 (d_N2, 1, 1),
1401 ),
1402 (1, mNrm_N, nNrm_N),
1403 )
1405 # The FOC is dvdn = tax*dvdm or dvdnNvrs = taxNvrs*dvdmNvrs
1406 dvdDNvrs = dvdnNvrs - taxNvrs_tiled * dvdmNvrs
1407 # The optimal d will be at the first point where dvdD < 0. The inverse
1408 # transformation flips the sign.
1410 # If the derivative is negative (inverse positive) at the lowest d,
1411 # then d == -1.0 is optimal
1412 constrained_bot = dvdDNvrs[0, :, :] >= 0.0
1413 # If it is positive (inverse negative) at the highest d, then d[-1] = 1.0
1414 # is optimal
1415 constrained_top = (
1416 dvdDNvrs[
1417 -1,
1418 :,
1419 :,
1420 ]
1421 <= 0.0
1422 )
1424 # Find indices at which the derivative crosses 0 for the 1st time
1425 # will be 0 if it never does, but "constrained_top/bot" deals with that
1426 crossings = np.logical_and(dvdDNvrs[:-1, :, :] <= 0.0, dvdDNvrs[1:, :, :] >= 0.0)
1427 idx = np.argmax(crossings, axis=0)
1429 m_idx_tiled, n_idx_tiled = np.meshgrid(
1430 np.arange(mNrm_N), np.arange(nNrm_N), indexing="ij"
1431 )
1433 # Linearly interpolate the optimal withdrawal percentage d
1434 idx1 = idx + 1
1435 slopes = (
1436 dvdDNvrs[idx1, m_idx_tiled, n_idx_tiled]
1437 - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled]
1438 ) / (dfracGrid[idx1] - dfracGrid[idx])
1439 dfrac_opt = dfracGrid[idx] - dvdDNvrs[idx, m_idx_tiled, n_idx_tiled] / slopes
1441 # Replace the ones we knew were constrained
1442 dfrac_opt[constrained_bot] = dfracGrid[0]
1443 dfrac_opt[constrained_top] = dfracGrid[-1]
1445 # Find m_tilde and n_tilde
1446 mtil_opt, ntil_opt = rebalance_assets(
1447 dfrac_opt, mNrm_tiled[0], nNrm_tiled[0], WithdrawTax
1448 )
1450 # Now the derivatives. These are not straight forward because of corner
1451 # solutions with partial derivatives that change the limits. The idea then
1452 # is to evaluate the possible uses of the marginal unit of resources and
1453 # take the maximum.
1455 # An additional unit of m
1456 marg_m = dvdmFunc_Adj_next(mtil_opt, ntil_opt)
1457 # An additional unit of n kept in n
1458 marg_n = dvdnFunc_Adj_next(mtil_opt, ntil_opt)
1459 # An additional unit of n withdrawn to m
1460 marg_n_to_m = marg_m * (1 - WithdrawTax)
1462 # Marginal value is the maximum of the marginals in their possible uses
1463 dvdm_Adj = np.maximum(marg_m, marg_n)
1464 dvdmNvrs_Adj = uPinv(dvdm_Adj)
1465 dvdn_Adj = np.maximum(marg_n, marg_n_to_m)
1466 dvdnNvrs_Adj = uPinv(dvdn_Adj)
1468 # Interpolators
1470 # Value
1471 if vFuncBool:
1472 vNvrs_Adj = vFunc_Adj_next.vFuncNvrs(mtil_opt, ntil_opt)
1473 vNvrsFunc_Adj = BilinearInterp(vNvrs_Adj, mNrmGrid, nNrmGrid)
1474 vFunc_Adj = ValueFuncCRRA(vNvrsFunc_Adj, CRRA)
1475 else:
1476 vFunc_Adj = NullFunc()
1478 # Marginals
1479 dvdmFunc_Adj = MargValueFuncCRRA(
1480 BilinearInterp(dvdmNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1481 )
1482 dvdnFunc_Adj = MargValueFuncCRRA(
1483 BilinearInterp(dvdnNvrs_Adj, mNrmGrid, nNrmGrid), CRRA
1484 )
1486 # Decison
1487 dfracFunc_Adj = BilinearInterp(dfrac_opt, mNrmGrid, nNrmGrid)
1489 solution = RiskyContribRebSolution(
1490 # Rebalancing stage adjusting
1491 vFunc_Adj=vFunc_Adj,
1492 dfracFunc_Adj=dfracFunc_Adj,
1493 dvdmFunc_Adj=dvdmFunc_Adj,
1494 dvdnFunc_Adj=dvdnFunc_Adj,
1495 # Rebalancing stage fixed (nothing happens, so value functions are
1496 # the ones from the next stage)
1497 vFunc_Fxd=vFunc_Fxd_next,
1498 dfracFunc_Fxd=ConstantFunction(0.0),
1499 dvdmFunc_Fxd=dvdmFunc_Fxd_next,
1500 dvdnFunc_Fxd=dvdnFunc_Fxd_next,
1501 dvdsFunc_Fxd=dvdsFunc_Fxd_next,
1502 )
1504 return solution
1507def solveRiskyContrib(
1508 solution_next,
1509 ShockDstn,
1510 IncShkDstn,
1511 RiskyDstn,
1512 IndepDstnBool,
1513 LivPrb,
1514 DiscFac,
1515 CRRA,
1516 Rfree,
1517 PermGroFac,
1518 WithdrawTax,
1519 BoroCnstArt,
1520 aXtraGrid,
1521 nNrmGrid,
1522 mNrmGrid,
1523 ShareGrid,
1524 dfracGrid,
1525 vFuncBool,
1526 AdjustPrb,
1527 DiscreteShareBool,
1528 joint_dist_solver,
1529):
1530 """
1531 Solve a full period (with its three stages) of the agent's problem
1533 Parameters
1534 ----------
1535 solution_next : RiskyContribSolution
1536 Solution to next period's problem.
1537 ShockDstn : DiscreteDistribution
1538 Joint distribution of next period's (0) permanent income shock, (1)
1539 transitory income shock, and (2) risky asset return factor.
1540 IncShkDstn : DiscreteDistribution
1541 Joint distribution of next period's (0) permanent income shock and (1)
1542 transitory income shock.
1543 RiskyDstn : DiscreteDistribution
1544 Distribution of next period's risky asset return factor.
1545 IndepDstnBool : bool
1546 Indicates whether the income and risky return distributions are
1547 independent.
1548 LivPrb : float
1549 Probability of surviving until next period.
1550 DiscFac : float
1551 Time-preference discount factor.
1552 CRRA : float
1553 Coefficient of relative risk aversion.
1554 Rfree : float
1555 Risk-free return factor.
1556 PermGroFac : float
1557 Deterministic permanent income growth factor.
1558 WithdrawTax : float
1559 Tax rate on risky asset withdrawals.
1560 BoroCnstArt : float
1561 Minimum allowed market resources (must be 0).
1562 aXtraGrid : numpy array
1563 Exogenous grid for end-of-period risk free resources.
1564 nNrmGrid : numpy array
1565 Exogenous grid for risky resources.
1566 mNrmGrid : numpy array
1567 Exogenous grid for risk-free resources.
1568 ShareGrid : numpy array
1569 Exogenous grid for the income contribution share.
1570 dfracGrid : numpy array
1571 Grid for rebalancing flows. The final grid will be equivalent to
1572 [-nNrm*dfracGrid, dfracGrid*mNrm].
1573 vFuncBool : bool
1574 Determines whether the level of th value function must be computed.
1575 AdjustPrb : float
1576 Probability that the agent will be able to rebalance his portfolio
1577 next period.
1578 DiscreteShareBool : bool
1579 Boolean that determines whether only a discrete set of contribution
1580 shares (ShareGrid) is allowed.
1581 joint_dist_solver: bool
1582 Should the general solver be used even if income and returns are
1583 independent?
1585 Returns
1586 -------
1587 periodSol : RiskyContribSolution
1588 Solution to the agent's current-period problem.
1590 """
1591 # Pack parameters to be passed to stage-specific solvers
1592 kws = {
1593 "ShockDstn": ShockDstn,
1594 "IncShkDstn": IncShkDstn,
1595 "RiskyDstn": RiskyDstn,
1596 "IndepDstnBool": IndepDstnBool,
1597 "LivPrb": LivPrb,
1598 "DiscFac": DiscFac,
1599 "CRRA": CRRA,
1600 "Rfree": Rfree,
1601 "PermGroFac": PermGroFac,
1602 "WithdrawTax": WithdrawTax,
1603 "BoroCnstArt": BoroCnstArt,
1604 "aXtraGrid": aXtraGrid,
1605 "nNrmGrid": nNrmGrid,
1606 "mNrmGrid": mNrmGrid,
1607 "ShareGrid": ShareGrid,
1608 "dfracGrid": dfracGrid,
1609 "vFuncBool": vFuncBool,
1610 "AdjustPrb": AdjustPrb,
1611 "DiscreteShareBool": DiscreteShareBool,
1612 "joint_dist_solver": joint_dist_solver,
1613 }
1615 # Stages of the problem in chronological order
1616 Stages = ["Reb", "Sha", "Cns"]
1617 n_stages = len(Stages)
1618 # Solvers, indexed by stage names
1619 Solvers = {
1620 "Reb": solve_RiskyContrib_Reb,
1621 "Sha": solve_RiskyContrib_Sha,
1622 "Cns": solve_RiskyContrib_Cns,
1623 }
1625 # Initialize empty solution
1626 stage_sols = {}
1627 # Solve stages backwards
1628 for i in reversed(range(n_stages)):
1629 stage = Stages[i]
1631 # In the last stage, the next solution is the first stage of the next
1632 # period. Otherwise, its the next stage of his period.
1633 if i == n_stages - 1:
1634 sol_next_stage = solution_next.stage_sols[Stages[0]]
1635 else:
1636 sol_next_stage = stage_sols[Stages[i + 1]]
1638 # Solve
1639 stage_sols[stage] = Solvers[stage](sol_next_stage, **kws)
1641 # Assemble stage solutions into period solution
1642 periodSol = RiskyContribSolution(**stage_sols)
1644 return periodSol
1647# %% Base risky-contrib dictionaries
1649risky_contrib_constructor_dict = {
1650 "IncShkDstn": construct_lognormal_income_process_unemployment,
1651 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1652 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1653 "aXtraGrid": make_assets_grid,
1654 "RiskyDstn": make_lognormal_RiskyDstn,
1655 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
1656 "ShareLimit": calc_ShareLimit_for_CRRA,
1657 "AdjustDstn": make_AdjustDstn,
1658 "solution_terminal": make_solution_terminal_risky_contrib,
1659 "ShareGrid": make_bounded_ShareGrid,
1660 "dfracGrid": make_simple_dGrid,
1661 "mNrmGrid": make_mNrm_grid,
1662 "nNrmGrid": make_nNrm_grid,
1663 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1664 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1665}
1667risky_contrib_params = {
1668 "constructors": risky_contrib_constructor_dict,
1669 # Preferences. The points of the model are more evident for more risk
1670 # averse and impatient agents
1671 "CRRA": 5.0,
1672 "DiscFac": 0.90,
1673 "WithdrawTax": [0.1],
1674 # Artificial borrowing constraint must be on
1675 "BoroCnstArt": 0.0,
1676 # Grids go up high wealth/P ratios and are less clustered at the bottom.
1677 "aXtraMax": 250,
1678 "aXtraCount": 50,
1679 "aXtraNestFac": 1,
1680 # Same goes for the new grids of the model
1681 "mNrmMin": 1e-6,
1682 "mNrmMax": 250,
1683 "mNrmCount": 50,
1684 "mNrmNestFac": 1,
1685 "nNrmMin": 1e-6,
1686 "nNrmMax": 250,
1687 "nNrmCount": 50,
1688 "nNrmNestFac": 1,
1689 # Income deduction/contribution share grid
1690 "ShareCount": 10,
1691 "ShareMax": 0.9,
1692 "DiscreteShareBool": False,
1693 # Grid for finding the optimal rebalancing flow
1694 "dCount": 20,
1695 "joint_dist_solver": False,
1696}
1697risky_asset_params = {
1698 # Risky return factor moments. Based on SP500 real returns from Shiller's
1699 # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
1700 "RiskyAvg": 1.080370891,
1701 "RiskyStd": 0.177196585,
1702 "ShareCount": 25, # Number of discrete points in the risky share approximation
1703 # Number of integration nodes to use in approximation of risky returns
1704 "RiskyCount": 5,
1705 # Probability that the agent can adjust their portfolio each period
1706 "AdjustPrb": [1.0],
1707 # When simulating the model, should all agents get the same risky return in
1708 # a given period?
1709 "sim_common_Rrisky": True,
1710}
1712# Infinite horizon version
1713init_risky_contrib = init_risky_asset.copy()
1714init_risky_contrib.update(risky_contrib_params)
1716# Lifecycle version
1717init_risky_contrib_lifecycle = init_lifecycle.copy()
1718init_risky_contrib_lifecycle.update(risky_asset_params)
1719init_risky_contrib_lifecycle.update(risky_contrib_params)
1721###############################################################################
1724class RiskyContribConsumerType(RiskyAssetConsumerType):
1725 """
1726 A consumer type with idiosyncratic shocks to permanent and transitory income,
1727 who can save in both a risk-free and a risky asset but faces frictions to
1728 moving funds between them. The agent can only consume out of his risk-free
1729 asset.
1731 The frictions are:
1733 - A proportional tax on funds moved from the risky to the risk-free
1734 asset.
1735 - A stochastic inability to move funds between his accounts.
1737 To partially avoid the second friction, the agent can commit to have a
1738 fraction of his labor income, which is usually deposited in his risk-free
1739 account, diverted to his risky account. He can change this fraction
1740 only in periods where he is able to move funds between accounts.
1741 """
1743 # The model is solved and simulated spliting each of the agent's
1744 # decisions into its own "stage". The stages in chronological order
1745 # are
1746 # - Reb: asset-rebalancing stage.
1747 # - Sha: definition of the income contribution share.
1748 # - Cns: consumption stage.
1749 stages = ["Reb", "Sha", "Cns"]
1750 # Each stage has its own states and controls, and its methods to find them.
1752 time_inv_ = RiskyAssetConsumerType.time_inv_ + [
1753 "DiscreteShareBool",
1754 "joint_dist_solver",
1755 "ShareGrid",
1756 "nNrmGrid",
1757 "mNrmGrid",
1758 "RiskyDstn",
1759 "dfracGrid",
1760 ]
1761 time_vary_ = RiskyAssetConsumerType.time_vary_ + ["WithdrawTax", "AdjustPrb"]
1763 # The new state variables (over those in ConsIndShock) are:
1764 # - nNrm: start-of-period risky resources.
1765 # - mNrmTilde: post-rebalancing risk-free resources.
1766 # - nNrmTilde: post-rebalancing risky resources.
1767 # - Share: income-deduction share.
1768 # For details, see
1769 # https://github.com/Mv77/RiskyContrib/blob/main/RiskyContrib.pdf
1770 state_vars = RiskyAssetConsumerType.state_vars + [
1771 "gNrm",
1772 "nNrm",
1773 "mNrmTilde",
1774 "nNrmTilde",
1775 "Share",
1776 ]
1777 shock_vars_ = RiskyAssetConsumerType.shock_vars_
1778 default_ = {"params": init_risky_contrib, "solver": solveRiskyContrib}
1780 def __init__(self, **kwds):
1781 super().__init__(**kwds)
1782 # It looks like I can't assign this at the class level, unfortunately
1783 self.get_states = {
1784 "Reb": self.get_states_Reb,
1785 "Sha": self.get_states_Sha,
1786 "Cns": self.get_states_Cns,
1787 }
1788 self.get_controls = {
1789 "Reb": self.get_controls_Reb,
1790 "Sha": self.get_controls_Sha,
1791 "Cns": self.get_controls_Cns,
1792 }
1794 def pre_solve(self):
1795 self.construct("solution_terminal")
1797 def initialize_sim(self):
1798 """
1799 Initialize the state of simulation attributes.
1801 Parameters
1802 ----------
1803 None
1805 Returns
1806 -------
1807 None
1808 """
1809 RiskyAssetConsumerType.initialize_sim(self)
1810 self.state_now["Share"] = np.zeros(self.AgentCount)
1812 def sim_birth(self, which_agents):
1813 """
1814 Create new agents to replace ones who have recently died; takes draws of
1815 initial aNrm and pLvl, as in ConsIndShockModel, then sets Share, Adjust
1816 and post-rebalancing risky asset nNrmTilde to zero as initial values.
1817 Parameters
1818 ----------
1819 which_agents : np.array
1820 Boolean array of size AgentCount indicating which agents should be "born".
1822 Returns
1823 -------
1824 None
1825 """
1827 RiskyAssetConsumerType.sim_birth(self, which_agents)
1828 self.state_now["Share"][which_agents] = 0.0
1829 self.state_now["nNrmTilde"][which_agents] = 0.0
1831 def sim_one_period(self):
1832 """
1833 Simulates one period for this type.
1835 Has to be re-defined instead of using AgentType.sim_one_period() because
1836 of the "stages" structure.
1838 Parameters
1839 ----------
1840 None
1841 Returns
1842 -------
1843 None
1844 """
1846 if not hasattr(self, "solution"):
1847 raise Exception(
1848 "Model instance does not have a solution stored. To simulate, it is necessary"
1849 " to run the `solve()` method of the class first."
1850 )
1852 # Mortality adjusts the agent population
1853 self.get_mortality() # Replace some agents with "newborns"
1855 # Make state_now into state_prev, clearing state_now
1856 for var in self.state_now:
1857 self.state_prev[var] = self.state_now[var]
1859 if isinstance(self.state_now[var], np.ndarray):
1860 self.state_now[var] = np.empty(self.AgentCount)
1861 else:
1862 # Probably an aggregate variable. It may be getting set by the Market.
1863 pass
1865 if self.read_shocks: # If shock histories have been pre-specified, use those
1866 self.read_shocks_from_history()
1867 else: # Otherwise, draw shocks as usual according to subclass-specific method
1868 self.get_shocks()
1870 # Sequentially get states and controls of every stage
1871 for s in self.stages:
1872 self.get_states[s]()
1873 self.get_controls[s]()
1875 self.get_post_states()
1877 # Advance time for all agents
1878 self.t_age = self.t_age + 1 # Age all consumers by one period
1879 self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle
1880 self.t_cycle[self.t_cycle == self.T_cycle] = (
1881 0 # Resetting to zero for those who have reached the end
1882 )
1884 def get_states_Reb(self):
1885 """
1886 Get states for the first "stage": rebalancing.
1887 """
1889 pLvlPrev = self.state_prev["pLvl"]
1890 aNrmPrev = self.state_prev["aNrm"]
1891 SharePrev = self.state_prev["Share"]
1892 nNrmTildePrev = self.state_prev["nNrmTilde"]
1893 Rfree = self.get_Rfree()
1894 Rrisk = self.shocks["Risky"]
1896 # Calculate new states:
1898 # Permanent income
1899 self.state_now["pLvl"] = pLvlPrev * self.shocks["PermShk"]
1900 self.state_now["PlvlAgg"] = self.state_prev["PlvlAgg"] * self.PermShkAggNow
1902 # Assets: mNrm and nNrm
1904 # Compute the effective growth factor of each asset
1905 RfEff = Rfree / self.shocks["PermShk"]
1906 RrEff = Rrisk / self.shocks["PermShk"]
1908 self.state_now["bNrm"] = RfEff * aNrmPrev # Liquid balances before labor income
1909 self.state_now["gNrm"] = (
1910 RrEff * nNrmTildePrev
1911 ) # Iliquid balances before labor income
1913 # Liquid balances after labor income
1914 self.state_now["mNrm"] = self.state_now["bNrm"] + self.shocks["TranShk"] * (
1915 1 - SharePrev
1916 )
1917 # Iliquid balances after labor income
1918 self.state_now["nNrm"] = (
1919 self.state_now["gNrm"] + self.shocks["TranShk"] * SharePrev
1920 )
1922 return None
1924 def get_controls_Reb(self):
1925 """
1926 Get controls for the first stage: rebalancing
1927 """
1928 dfrac = np.zeros(self.AgentCount) + np.nan
1930 # Loop over each period of the cycle, getting controls separately depending on "age"
1931 for t in range(self.T_cycle):
1932 # Find agents in this period-stage
1933 these = t == self.t_cycle
1935 # Get controls for agents who *can* adjust.
1936 those = np.logical_and(these, self.shocks["Adjust"])
1937 dfrac[those] = (
1938 self.solution[t]
1939 .stage_sols["Reb"]
1940 .dfracFunc_Adj(
1941 self.state_now["mNrm"][those], self.state_now["nNrm"][those]
1942 )
1943 )
1945 # Get Controls for agents who *can't* adjust.
1946 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
1947 dfrac[those] = (
1948 self.solution[t]
1949 .stage_sols["Reb"]
1950 .dfracFunc_Fxd(
1951 self.state_now["mNrm"][those],
1952 self.state_now["nNrm"][those],
1953 self.state_prev["Share"][those],
1954 )
1955 )
1957 # Limit dfrac to [-1,1] to prevent negative balances. Values outside
1958 # the range can come from extrapolation.
1959 self.controls["dfrac"] = np.minimum(np.maximum(dfrac, -1), 1.0)
1961 def get_states_Sha(self):
1962 """
1963 Get states for the second "stage": choosing the contribution share.
1964 """
1966 # Post-states are assets after rebalancing
1968 if "WithdrawTax" not in self.time_vary:
1969 mNrmTilde, nNrmTilde = rebalance_assets(
1970 self.controls["dfrac"],
1971 self.state_now["mNrm"],
1972 self.state_now["nNrm"],
1973 self.WithdrawTax,
1974 )
1976 else:
1977 # Initialize
1978 mNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1979 nNrmTilde = np.zeros_like(self.state_now["mNrm"]) + np.nan
1981 # Loop over each period of the cycle, getting controls separately depending on "age"
1982 for t in range(self.T_cycle):
1983 # Find agents in this period-stage
1984 these = t == self.t_cycle
1986 if np.sum(these) > 0:
1987 WithdrawTax = self.WithdrawTax[t]
1989 mNrmTilde[these], nNrmTilde[these] = rebalance_assets(
1990 self.controls["dfrac"][these],
1991 self.state_now["mNrm"][these],
1992 self.state_now["nNrm"][these],
1993 WithdrawTax,
1994 )
1996 self.state_now["mNrmTilde"] = mNrmTilde
1997 self.state_now["nNrmTilde"] = nNrmTilde
1999 def get_controls_Sha(self):
2000 """
2001 Get controls for the second "stage": choosing the contribution share.
2002 """
2004 Share = np.zeros(self.AgentCount) + np.nan
2006 # Loop over each period of the cycle, getting controls separately depending on "age"
2007 for t in range(self.T_cycle):
2008 # Find agents in this period-stage
2009 these = t == self.t_cycle
2011 # Get controls for agents who *can* adjust.
2012 those = np.logical_and(these, self.shocks["Adjust"])
2013 Share[those] = (
2014 self.solution[t]
2015 .stage_sols["Sha"]
2016 .ShareFunc_Adj(
2017 self.state_now["mNrmTilde"][those],
2018 self.state_now["nNrmTilde"][those],
2019 )
2020 )
2022 # Get Controls for agents who *can't* adjust.
2023 those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
2024 Share[those] = (
2025 self.solution[t]
2026 .stage_sols["Sha"]
2027 .ShareFunc_Fxd(
2028 self.state_now["mNrmTilde"][those],
2029 self.state_now["nNrmTilde"][those],
2030 self.state_prev["Share"][those],
2031 )
2032 )
2034 # Store controls as attributes of self
2035 self.controls["Share"] = Share
2037 def get_states_Cns(self):
2038 """
2039 Get states for the third "stage": consumption.
2040 """
2042 # Contribution share becomes a state in the consumption problem
2043 self.state_now["Share"] = self.controls["Share"]
2045 def get_controls_Cns(self):
2046 """
2047 Get controls for the third "stage": consumption.
2048 """
2050 cNrm = np.zeros(self.AgentCount) + np.nan
2052 # Loop over each period of the cycle, getting controls separately depending on "age"
2053 for t in range(self.T_cycle):
2054 # Find agents in this period-stage
2055 these = t == self.t_cycle
2057 # Get consumption
2058 cNrm[these] = (
2059 self.solution[t]
2060 .stage_sols["Cns"]
2061 .cFunc(
2062 self.state_now["mNrmTilde"][these],
2063 self.state_now["nNrmTilde"][these],
2064 self.state_now["Share"][these],
2065 )
2066 )
2068 # Store controls as attributes of self
2069 # Since agents might be willing to end the period with a = 0, make
2070 # sure consumption does not go over m because of some numerical error.
2071 self.controls["cNrm"] = np.minimum(cNrm, self.state_now["mNrmTilde"])
2073 def get_post_states(self):
2074 """
2075 Set variables that are not a state to any problem but need to be
2076 computed in order to interact with shocks and produce next period's
2077 states.
2078 """
2079 self.state_now["aNrm"] = self.state_now["mNrmTilde"] - self.controls["cNrm"]