diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index c322ea8be..c7d142f25 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -22,6 +22,7 @@ Release Data: TBD frictions to moving funds between them. To circumvent these frictions, he has access to an income-deduction scheme to accumulate risky assets. PR: [#832](https://github.com/econ-ark/HARK/pull/832). See [this forthcoming REMARK](https://github.com/Mv77/RiskyContrib) for the model's details. * 'cycles' agent property moved from constructor argument to parameter [#1031](https://github.com/econ-ark/HARK/pull/1031) +* Uses iterated expectations to speed-up the solution of `RiskyContrib` when income and returns are independent [#1058](https://github.com/econ-ark/HARK/pull/1058). #### Minor Changes diff --git a/HARK/ConsumptionSaving/ConsRiskyContribModel.py b/HARK/ConsumptionSaving/ConsRiskyContribModel.py index 603bd0719..93edcb1e1 100644 --- a/HARK/ConsumptionSaving/ConsRiskyContribModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyContribModel.py @@ -1,7 +1,25 @@ """ This file contains classes and functions for representing, solving, and simulating -agents who must allocate their resources among consumption, saving in a risk-free -asset (with a low return), and saving in a risky asset (with higher average return). +a consumer type with idiosyncratic shocks to permanent and transitory income, +who can save in both a risk-free and a risky asset but faces frictions to +moving funds between them. The agent can only consume out of his risk-free +asset. + +The model is described in detail in the REMARK: +https://econ-ark.org/materials/riskycontrib + +@software{mateo_velasquez_giraldo_2021_4977915, + author = {Mateo Velásquez-Giraldo}, + title = {{Mv77/RiskyContrib: A Two-Asset Savings Model with + an Income-Contribution Scheme}}, + month = jun, + year = 2021, + publisher = {Zenodo}, + version = {v1.0.1}, + doi = {10.5281/zenodo.4977915}, + url = {https://doi.org/10.5281/zenodo.4977915} +} + """ import numpy as np from copy import deepcopy @@ -55,7 +73,7 @@ class RiskyContribConsumerType(RiskyAssetConsumerType): """ time_inv_ = deepcopy(RiskyAssetConsumerType.time_inv_) - time_inv_ = time_inv_ + ["DiscreteShareBool"] + time_inv_ = time_inv_ + ["DiscreteShareBool", "joint_dist_solver"] # The new state variables (over those in ConsIndShock) are: # - nMrm: start-of-period risky resources. @@ -73,7 +91,7 @@ class RiskyContribConsumerType(RiskyAssetConsumerType): ] shock_vars_ = RiskyAssetConsumerType.shock_vars_ - def __init__(self, verbose=False, quiet=False, **kwds): + def __init__(self, verbose=False, quiet=False, joint_dist_solver=False, **kwds): params = init_risky_contrib.copy() params.update(kwds) @@ -103,7 +121,12 @@ def __init__(self, verbose=False, quiet=False, **kwds): "Sha": self.get_controls_Sha, "Cns": self.get_controls_Cns, } - + + # The model can be solved more quickly if income and risky returns are + # independent. However, people might want to use the general solver + # even when they are independent for debugging and testing. + self.joint_dist_solver=joint_dist_solver + # Set the solver for the portfolio model, and update various constructed attributes self.solve_one_period = solveRiskyContrib self.update() @@ -476,8 +499,9 @@ def get_controls_Reb(self): ) ) - # Store controls as attributes of self - self.controls["dfrac"] = dfrac + # Limit dfrac to [-1,1] to prevent negative balances. Values outside + # the range can come from extrapolation. + self.controls["dfrac"] = np.minimum(np.maximum(dfrac,-1),1.0) def get_states_Sha(self): """ @@ -996,6 +1020,9 @@ def n_nrm_next(shocks, nNrm, Share, PermGroFac): def solve_RiskyContrib_Cns( solution_next, ShockDstn, + IncShkDstn, + RiskyDstn, + IndepDstnBool, LivPrb, DiscFac, CRRA, @@ -1009,6 +1036,7 @@ def solve_RiskyContrib_Cns( vFuncBool, AdjustPrb, DiscreteShareBool, + joint_dist_solver, **unused_params ): """ @@ -1021,6 +1049,14 @@ def solve_RiskyContrib_Cns( ShockDstn : DiscreteDistribution Joint distribution of next period's (0) permanent income shock, (1) transitory income shock, and (2) risky asset return factor. + IncShkDstn : DiscreteDistribution + Joint distribution of next period's (0) permanent income shock and (1) + transitory income shock. + RiskyDstn : DiscreteDistribution + Distribution of next period's risky asset return factor. + IndepDstnBool : bool + Indicates whether the income and risky return distributions are + independent. LivPrb : float Probability of surviving until next period. DiscFac : float @@ -1050,6 +1086,9 @@ def solve_RiskyContrib_Cns( DiscreteShareBool : bool Boolean that determines whether only a discrete set of contribution shares (ShareGrid) is allowed. + joint_dist_solver: bool + Should the general solver be used even if income and returns are + independent? Returns ------- @@ -1115,60 +1154,194 @@ def solve_RiskyContrib_Cns( if vFuncBool: v_next = lambda m, n, s: vFunc_Reb_Adj_next(m, n) - # Now construct a function that evaluates and discounts them given a - # vector of return and income shocks and an end-of-period state - def end_of_period_derivs(shocks, a, nTil, s): - """ - Computes the end-of-period derivatives (and optionally the value) of the - continuation value function, conditional on shocks. This is so that the - expectations can be calculated by integrating over shocks. - Parameters - ---------- - shocks : np.array - Length-3 array with the stochastic shocks that get realized between the - end of the current period and the start of next period. Their order is - (0) permanent income shock, (1) transitory income shock, (2) risky - asset return. - a : float - end-of-period risk-free assets. - nTil : float - end-of-period risky assets. - s : float - end-of-period income deduction share. - """ - temp_fac_A = utilityP(shocks[0] * PermGroFac, CRRA) - temp_fac_B = (shocks[0] * PermGroFac) ** (1.0 - CRRA) - - # Find next-period asset balances - m_next = m_nrm_next(shocks, a, s, Rfree, PermGroFac) - n_next = n_nrm_next(shocks, nTil, s, PermGroFac) - - # Interpolate next-period-value derivatives - dvdm_tp1 = dvdm_next(m_next, n_next, s) - dvdn_tp1 = dvdn_next(m_next, n_next, s) - if shocks[1] == 0: - dvds_tp1 = dvds_next(m_next, n_next, s) - else: - dvds_tp1 = shocks[1] * (dvdn_tp1 - dvdm_tp1) + dvds_next(m_next, n_next, s) + if IndepDstnBool and not joint_dist_solver: + + # If income and returns are independent we can use the law of iterated + # expectations to speed up the computation of end-of-period derivatives + + # Define "post-return variables" + # b_aux = aNrm * R + # g_aux = nNrmTilde * Rtilde + # and create a function that interpolates end-of-period marginal values + # as functions of those and the contribution share + + def post_return_derivs(inc_shocks, b_aux, g_aux, s): + + perm_shk = inc_shocks[0] + tran_shk = inc_shocks[1] + + temp_fac_A = utilityP(perm_shk * PermGroFac, CRRA) + temp_fac_B = (perm_shk * PermGroFac) ** (1.0 - CRRA) + + # Find next-period asset balances + m_next = b_aux / (perm_shk * PermGroFac) + (1.0 - s) * tran_shk + n_next = g_aux / (perm_shk * PermGroFac) + s * tran_shk + + # Interpolate next-period-value derivatives + dvdm_tp1 = dvdm_next(m_next, n_next, s) + dvdn_tp1 = dvdn_next(m_next, n_next, s) + if tran_shk == 0: + dvds_tp1 = dvds_next(m_next, n_next, s) + else: + dvds_tp1 = tran_shk * (dvdn_tp1 - dvdm_tp1) + dvds_next( + m_next, n_next, s + ) + + # Discount next-period-value derivatives to current period + + # Liquid resources + pr_dvda = temp_fac_A * dvdm_tp1 + # Iliquid resources + pr_dvdn = temp_fac_A * dvdn_tp1 + # Contribution share + pr_dvds = temp_fac_B * dvds_tp1 + + # End of period value function, if needed + if vFuncBool: + pr_v = temp_fac_B * v_next(m_next, n_next, s) + return np.stack([pr_dvda, pr_dvdn, pr_dvds, pr_v]) + else: + return np.stack([pr_dvda, pr_dvdn, pr_dvds]) - # Discount next-period-value derivatives to current period + # Define grids + b_aux_grid = np.concatenate([np.array([0.0]), Rfree * aXtraGrid]) + g_aux_grid = np.concatenate([np.array([0.0]), max(RiskyDstn.X) * nNrmGrid]) - # Liquid resources - end_of_prd_dvda = DiscFac * Rfree * LivPrb * temp_fac_A * dvdm_tp1 - # Iliquid resources - end_of_prd_dvdn = DiscFac * shocks[2] * LivPrb * temp_fac_A * dvdn_tp1 - # Contribution share - end_of_prd_dvds = DiscFac * LivPrb * temp_fac_B * dvds_tp1 + # Create tiled arrays with conforming dimensions. + b_aux_tiled, g_aux_tiled, Share_tiled = np.meshgrid( + b_aux_grid, g_aux_grid, ShareGrid, indexing="ij" + ) + + # Find end of period derivatives and value as expectations of (discounted) + # next period's derivatives and value. + pr_derivs = calc_expectation( + IncShkDstn, post_return_derivs, b_aux_tiled, g_aux_tiled, Share_tiled + )[:, :, :, :, 0] + + # Unpack results and create interpolators + pr_dvdb_func = MargValueFuncCRRA( + TrilinearInterp(uPinv(pr_derivs[0]), b_aux_grid, g_aux_grid, ShareGrid), CRRA + ) + pr_dvdg_func = MargValueFuncCRRA( + TrilinearInterp(uPinv(pr_derivs[1]), b_aux_grid, g_aux_grid, ShareGrid), CRRA + ) + pr_dvds_func = TrilinearInterp(pr_derivs[2], b_aux_grid, g_aux_grid, ShareGrid) - # End of period value function, i11f needed if vFuncBool: - end_of_prd_v = DiscFac * LivPrb * temp_fac_B * v_next(m_next, n_next, s) - return np.stack( - [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v] + + pr_vFunc = ValueFuncCRRA( + TrilinearInterp(uInv(pr_derivs[3]), b_aux_grid, g_aux_grid, ShareGrid), CRRA ) - else: - return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds]) + + # Now construct a function that produces end-of-period derivatives + # given the risky return draw + def end_of_period_derivs(risky_ret, a, nTil, s): + """ + Computes the end-of-period derivatives (and optionally the value) of the + continuation value function, conditional on risky returns. This is so that the + expectations can be calculated by integrating over risky returns. + + Parameters + ---------- + risky_ret : float + Risky return factor + a : float + end-of-period risk-free assets. + nTil : float + end-of-period risky assets. + s : float + end-of-period income deduction share. + """ + + # Find next-period asset balances + b_aux = a * Rfree + g_aux = nTil * risky_ret + + # Interpolate post-return derivatives + pr_dvdb = pr_dvdb_func(b_aux, g_aux, s) + pr_dvdg = pr_dvdg_func(b_aux, g_aux, s) + pr_dvds = pr_dvds_func(b_aux, g_aux, s) + + # Discount + + # Liquid resources + end_of_prd_dvda = DiscFac * Rfree * LivPrb * pr_dvdb + # Iliquid resources + end_of_prd_dvdn = DiscFac * risky_ret * LivPrb * pr_dvdg + # Contribution share + end_of_prd_dvds = DiscFac * LivPrb * pr_dvds + + # End of period value function, i11f needed + if vFuncBool: + end_of_prd_v = DiscFac * LivPrb * pr_vFunc(b_aux, g_aux, s) + return np.stack( + [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v] + ) + else: + return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds]) + + else: + + # If income and returns are not independent, we just integrate over + # them jointly. + + # Construct a function that evaluates and discounts them given a + # vector of return and income shocks and an end-of-period state + def end_of_period_derivs(shocks, a, nTil, s): + """ + Computes the end-of-period derivatives (and optionally the value) of the + continuation value function, conditional on shocks. This is so that the + expectations can be calculated by integrating over shocks. + + Parameters + ---------- + shocks : np.array + Length-3 array with the stochastic shocks that get realized between the + end of the current period and the start of next period. Their order is + (0) permanent income shock, (1) transitory income shock, (2) risky + asset return. + a : float + end-of-period risk-free assets. + nTil : float + end-of-period risky assets. + s : float + end-of-period income deduction share. + """ + temp_fac_A = utilityP(shocks[0] * PermGroFac, CRRA) + temp_fac_B = (shocks[0] * PermGroFac) ** (1.0 - CRRA) + + # Find next-period asset balances + m_next = m_nrm_next(shocks, a, s, Rfree, PermGroFac) + n_next = n_nrm_next(shocks, nTil, s, PermGroFac) + + # Interpolate next-period-value derivatives + dvdm_tp1 = dvdm_next(m_next, n_next, s) + dvdn_tp1 = dvdn_next(m_next, n_next, s) + if shocks[1] == 0: + dvds_tp1 = dvds_next(m_next, n_next, s) + else: + dvds_tp1 = shocks[1] * (dvdn_tp1 - dvdm_tp1) + dvds_next( + m_next, n_next, s + ) + + # Discount next-period-value derivatives to current period + + # Liquid resources + end_of_prd_dvda = DiscFac * Rfree * LivPrb * temp_fac_A * dvdm_tp1 + # Iliquid resources + end_of_prd_dvdn = DiscFac * shocks[2] * LivPrb * temp_fac_A * dvdn_tp1 + # Contribution share + end_of_prd_dvds = DiscFac * LivPrb * temp_fac_B * dvds_tp1 + + # End of period value function, i11f needed + if vFuncBool: + end_of_prd_v = DiscFac * LivPrb * temp_fac_B * v_next(m_next, n_next, s) + return np.stack( + [end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds, end_of_prd_v] + ) + else: + return np.stack([end_of_prd_dvda, end_of_prd_dvdn, end_of_prd_dvds]) # Now find the expected values on a (a, nTil, s) grid @@ -1185,7 +1358,11 @@ def end_of_period_derivs(shocks, a, nTil, s): # Find end of period derivatives and value as expectations of (discounted) # next period's derivatives and value. eop_derivs = calc_expectation( - ShockDstn, end_of_period_derivs, aNrm_tiled, nNrm_tiled, Share_tiled + RiskyDstn if IndepDstnBool and not joint_dist_solver else ShockDstn, + end_of_period_derivs, + aNrm_tiled, + nNrm_tiled, + Share_tiled, )[:, :, :, :, 0] # Unpack results @@ -1684,6 +1861,9 @@ def solve_RiskyContrib_Reb( def solveRiskyContrib( solution_next, ShockDstn, + IncShkDstn, + RiskyDstn, + IndepDstnBool, LivPrb, DiscFac, CRRA, @@ -1699,6 +1879,7 @@ def solveRiskyContrib( vFuncBool, AdjustPrb, DiscreteShareBool, + joint_dist_solver, ): """ Solve a full period (with its three stages) of the agent's problem @@ -1710,6 +1891,14 @@ def solveRiskyContrib( ShockDstn : DiscreteDistribution Joint distribution of next period's (0) permanent income shock, (1) transitory income shock, and (2) risky asset return factor. + IncShkDstn : DiscreteDistribution + Joint distribution of next period's (0) permanent income shock and (1) + transitory income shock. + RiskyDstn : DiscreteDistribution + Distribution of next period's risky asset return factor. + IndepDstnBool : bool + Indicates whether the income and risky return distributions are + independent. LivPrb : float Probability of surviving until next period. DiscFac : float @@ -1743,6 +1932,9 @@ def solveRiskyContrib( DiscreteShareBool : bool Boolean that determines whether only a discrete set of contribution shares (ShareGrid) is allowed. + joint_dist_solver: bool + Should the general solver be used even if income and returns are + independent? Returns ------- @@ -1753,6 +1945,9 @@ def solveRiskyContrib( # Pack parameters to be passed to stage-specific solvers kws = { "ShockDstn": ShockDstn, + "IncShkDstn": IncShkDstn, + "RiskyDstn": RiskyDstn, + "IndepDstnBool": IndepDstnBool, "LivPrb": LivPrb, "DiscFac": DiscFac, "CRRA": CRRA, @@ -1768,6 +1963,7 @@ def solveRiskyContrib( "vFuncBool": vFuncBool, "AdjustPrb": AdjustPrb, "DiscreteShareBool": DiscreteShareBool, + "joint_dist_solver": joint_dist_solver, } # Stages of the problem in chronological order diff --git a/HARK/ConsumptionSaving/tests/test_ConsRiskyContribModel.py b/HARK/ConsumptionSaving/tests/test_ConsRiskyContribModel.py index 9b247bfdc..cf6e5e179 100644 --- a/HARK/ConsumptionSaving/tests/test_ConsRiskyContribModel.py +++ b/HARK/ConsumptionSaving/tests/test_ConsRiskyContribModel.py @@ -45,10 +45,25 @@ def test_finite_cont_share(self): cont_params = copy(self.par_finite) cont_params["DiscreteShareBool"] = False cont_params["vFuncBool"] = False - + + fin_cont_agent = RiskyContribConsumerType(**cont_params) + + # Independent solver + fin_cont_agent.solve() + self.assertAlmostEqual( + fin_cont_agent.solution[0].stage_sols["Reb"].dfracFunc_Adj(3., 4.), -0.87671241 + ) + self.assertAlmostEqual( + fin_cont_agent.solution[0].stage_sols["Sha"].ShareFunc_Adj(5., 0.1), 0.14641409 + ) + self.assertAlmostEqual( + fin_cont_agent.solution[0].stage_sols["Cns"].cFunc(3., 4., 0.1), 2.4560881 + ) + + # General correlated solver + fin_cont_agent.joint_dist_solver = True fin_cont_agent.solve() - self.assertAlmostEqual( fin_cont_agent.solution[0].stage_sols["Reb"].dfracFunc_Adj(3, 4), -0.87848691 ) @@ -66,6 +81,22 @@ def test_finite_disc_share(self): disc_params["vFuncBool"] = True fin_disc_agent = RiskyContribConsumerType(**disc_params) + + # Independent solver + fin_disc_agent.solve() + + self.assertAlmostEqual( + fin_disc_agent.solution[0].stage_sols["Reb"].dfracFunc_Adj(3., 4.), -0.8767603 + ) + self.assertAlmostEqual( + fin_disc_agent.solution[0].stage_sols["Sha"].ShareFunc_Adj(5., 0.1), 0.1 + ) + self.assertAlmostEqual( + fin_disc_agent.solution[0].stage_sols["Cns"].cFunc(3., 4., 0.1), 2.45608803 + ) + + # General correlated solver + fin_disc_agent.joint_dist_solver = True fin_disc_agent.solve() self.assertAlmostEqual(