From fb74f6868aa074080569595f7c0df0432a76c065 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Fri, 8 Mar 2024 17:46:52 -0500 Subject: [PATCH 01/14] add functional solvers --- HARK/ConsumptionSaving/ConsIndShockModel.py | 280 +++++++++++--------- 1 file changed, 160 insertions(+), 120 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index bc4edd291..66c72ef27 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -16,16 +16,6 @@ from copy import copy, deepcopy import numpy as np -from scipy import sparse as sp -from scipy.optimize import newton - -from HARK import ( - AgentType, - NullFunc, - _log, - make_one_period_oo_solver, - set_verbosity_level, -) from HARK.Calibration.Income.IncomeTools import ( Cagetti_income, parse_income_spec, @@ -72,6 +62,16 @@ jump_to_grid_2D, make_grid_exp_mult, ) +from scipy import sparse as sp +from scipy.optimize import newton + +from HARK import ( + AgentType, + NullFunc, + _log, + make_one_period_oo_solver, + set_verbosity_level, +) __all__ = [ "ConsumerSolution", @@ -209,11 +209,29 @@ def append_solution(self, new_solution): # ===================================================================== +def calc_human_wealth(h_nrm_next, perm_gro_fac, rfree, ex_inc_next): + return (perm_gro_fac / rfree) * (h_nrm_next + ex_inc_next) + + +def calc_pat_fac(rfree, disc_fac_eff, crra): + return ((rfree * disc_fac_eff) ** (1.0 / crra)) / rfree + + +def calc_mpc_min(mpc_min_next, pat_fac): + return 1.0 / (1.0 + pat_fac / mpc_min_next) + + def solve_one_period_ConsPF( - solution_next, DiscFac, LivPrb, CRRA, Rfree, PermGroFac, BoroCnstArt, MaxKinks + solution_next, + DiscFac, + LivPrb, + CRRA, + Rfree, + PermGroFac, + BoroCnstArt, + MaxKinks, ): - """ - Solves one period of a basic perfect foresight consumption-saving model with + """Solves one period of a basic perfect foresight consumption-saving model with a single risk free asset and permanent income growth. Parameters @@ -244,30 +262,29 @@ def solve_one_period_ConsPF( solution_now : ConsumerSolution Solution to the current period of a perfect foresight consumption-saving problem. + """ # Define the utility function and effective discount factor uFunc = UtilityFuncCRRA(CRRA) DiscFacEff = DiscFac * LivPrb # Effective = pure x LivPrb # Prevent comparing None and float if there is no borrowing constraint - if BoroCnstArt is None: - BoroCnstArt = -np.inf # Can borrow as much as we want + # Can borrow as much as we want + BoroCnstArt = -np.inf if BoroCnstArt is None else BoroCnstArt # Calculate human wealth this period - hNrmNow = (PermGroFac / Rfree) * (solution_next.hNrm + 1.0) + hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, 1.0) # Calculate the lower bound of the marginal propensity to consume - PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree - MPCmin = 1.0 / (1.0 + PatFac / solution_next.MPCmin) + PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) + MPCmin = calc_mpc_min(solution_next.MPCmin, PatFac) # Extract the discrete kink points in next period's consumption function; # don't take the last one, as it only defines the extrapolation and is not a kink. mNrmNext = solution_next.cFunc.x_list[:-1] cNrmNext = solution_next.cFunc.y_list[:-1] - vNext = PermGroFac ** (1.0 - CRRA) * uFunc( - solution_next.vFunc.vFuncNvrs.y_list[:-1] - ) - EndOfPrdv = DiscFacEff * vNext + vFuncNvrsNext = solution_next.vFunc.vFuncNvrs.y_list[:-1] + EndOfPrdv = DiscFacEff * PermGroFac ** (1.0 - CRRA) * uFunc(vFuncNvrsNext) # Calculate the end-of-period asset values that would reach those kink points # next period, then invert the first order condition to get consumption. Then @@ -366,6 +383,52 @@ def solve_one_period_ConsPF( return solution_now +def calc_worst_inc_prob(inc_shk_dstn): + probs = inc_shk_dstn.pmv + perm, tran = inc_shk_dstn.atoms + income = perm * tran + worst_inc = np.min(income) + return np.sum(probs[income == worst_inc]) + + +def calc_mpc_max(mpc_max_next, worst_inc_prob, crra, pat_fac): + temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac + return 1.0 / (1.0 + temp_fac / mpc_max_next) + + +def calc_boro_const_nat( + m_nrm_min_next, + inc_shk_dstn, + rfree, + perm_gro_fac, +): + perm, tran = inc_shk_dstn.atoms + temp_fac = (perm_gro_fac * np.min(perm)) / rfree + return (m_nrm_min_next - np.min(tran)) * temp_fac + + +def calc_m_nrm_next(shock, a, rfree, perm_gro_fac): + return rfree / (perm_gro_fac * shock["PermShk"]) * a + shock["TranShk"] + + +def calc_v_next(shock, a, rfree, crra, perm_gro_fac, vfunc_next): + return ( + shock["PermShk"] ** (1.0 - crra) * perm_gro_fac ** (1.0 - crra) + ) * vfunc_next(calc_m_nrm_next(shock, a, rfree, perm_gro_fac)) + + +def calc_vp_next(shock, a, rfree, crra, perm_gro_fac, vp_func_next): + return shock["PermShk"] ** (-crra) * vp_func_next( + calc_m_nrm_next(shock, a, rfree, perm_gro_fac), + ) + + +def calc_vpp_next(shock, a, rfree, crra, perm_gro_fac, vppfunc_next): + return shock["PermShk"] ** (-crra - 1.0) * vppfunc_next( + calc_m_nrm_next(shock, a, rfree, perm_gro_fac), + ) + + def solve_one_period_ConsIndShock( solution_next, IncShkDstn, @@ -379,8 +442,7 @@ def solve_one_period_ConsIndShock( vFuncBool, CubicBool, ): - """ - Solves one period of a consumption-saving model with idiosyncratic shocks to + """Solves one period of a consumption-saving model with idiosyncratic shocks to permanent and transitory income, with one risk free asset and CRRA utility. Parameters @@ -420,23 +482,14 @@ def solve_one_period_ConsIndShock( ------- solution_now : ConsumerSolution Solution to this period's consumption-saving problem with income risk. + """ # Define the current period utility function and effective discount factor uFunc = UtilityFuncCRRA(CRRA) DiscFacEff = DiscFac * LivPrb # "effective" discount factor - # Unpack next period's income shock distribution - ShkPrbsNext = IncShkDstn.pmv - PermShkValsNext = IncShkDstn.atoms[0] - TranShkValsNext = IncShkDstn.atoms[1] - PermShkMinNext = np.min(PermShkValsNext) - TranShkMinNext = np.min(TranShkValsNext) - # Calculate the probability that we get the worst possible income draw - IncNext = PermShkValsNext * TranShkValsNext - WorstIncNext = PermShkMinNext * TranShkMinNext - WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext]) - # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing + WorstIncPrb = calc_worst_inc_prob(IncShkDstn) # Unpack next period's (marginal) value function vFuncNext = solution_next.vFunc # This is None when vFuncBool is False @@ -444,21 +497,19 @@ def solve_one_period_ConsIndShock( vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False # Update the bounding MPCs and PDV of human wealth: - PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree - try: - MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) - except: - MPCminNow = 0.0 - Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext) - hNrmNow = PermGroFac / Rfree * (Ex_IncNext + solution_next.hNrm) - temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac - MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax) - cFuncLimitIntercept = MPCminNow * hNrmNow - cFuncLimitSlope = MPCminNow + PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) + MPCmin = calc_mpc_min(solution_next.MPCmin, PatFac) + Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) + hNrm = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, Ex_IncNext) + MPCmax = calc_mpc_max(solution_next.MPCmax, WorstIncPrb, CRRA, PatFac) + + cFuncLimitIntercept = MPCmin * hNrm + cFuncLimitSlope = MPCmin # Calculate the minimum allowable value of money resources in this period - PermGroFacEffMin = (PermGroFac * PermShkMinNext) / Rfree - BoroCnstNat = (solution_next.mNrmMin - TranShkMinNext) * PermGroFacEffMin + BoroCnstNat = calc_boro_const_nat( + solution_next.mNrmMin, IncShkDstn, Rfree, PermGroFac + ) # Set the minimum allowable (normalized) market resources based on the natural # and artificial borrowing constraints @@ -472,34 +523,24 @@ def solve_one_period_ConsIndShock( if BoroCnstNat < mNrmMinNow: MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1 else: - MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above + MPCmaxEff = MPCmax # Otherwise, it's the MPC calculated above # Define the borrowing-constrained consumption function cFuncNowCnst = LinearInterp( - np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0]) + np.array([mNrmMinNow, mNrmMinNow + 1.0]), + np.array([0.0, 1.0]), ) # Construct the assets grid by adjusting aXtra by the natural borrowing constraint aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat - # Define local functions for taking future expectations - def calc_mNrmNext(S, a, R): - return R / (PermGroFac * S["PermShk"]) * a + S["TranShk"] - - def calc_vNext(S, a, R): - return (S["PermShk"] ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)) * vFuncNext( - calc_mNrmNext(S, a, R) - ) - - def calc_vPnext(S, a, R): - return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a, R)) - - def calc_vPPnext(S, a, R): - return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, a, R)) - # Calculate end-of-period marginal value of assets at each gridpoint vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA) - EndOfPrdvP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(aNrmNow, Rfree)) + EndOfPrdvP = vPfacEff * expected( + calc_vp_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vPfuncNext), + ) # Invert the first order condition to find optimal cNrm from each aNrm gridpoint cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0)) @@ -514,11 +555,13 @@ def calc_vPPnext(S, a, R): # Calculate end-of-period marginal marginal value of assets at each gridpoint vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0) EndOfPrdvPP = vPPfacEff * expected( - calc_vPPnext, IncShkDstn, args=(aNrmNow, Rfree) + calc_vpp_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vPPfuncNext), ) dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2) MPC = dcda / (dcda + 1.0) - MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow) + MPC_for_interpolation = np.insert(MPC, 0, MPCmax) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( @@ -553,9 +596,13 @@ def calc_vPPnext(S, a, R): # Construct this period's value function if requested if vFuncBool: # Calculate end-of-period value, its derivative, and their pseudo-inverse - EndOfPrdv = DiscFacEff * expected(calc_vNext, IncShkDstn, args=(aNrmNow, Rfree)) + EndOfPrdv = DiscFacEff * expected( + calc_v_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vFuncNext), + ) EndOfPrdvNvrs = uFunc.inv( - EndOfPrdv + EndOfPrdv, ) # value transformed through inverse utility EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1)) EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0) @@ -580,9 +627,13 @@ def calc_vPPnext(S, a, R): mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow) vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0) vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA))) - MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA)) + MPCminNvrs = MPCmin ** (-CRRA / (1.0 - CRRA)) vNvrsFuncNow = CubicInterp( - mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs + mNrm_temp, + vNvrs_temp, + vNvrsP_temp, + MPCminNvrs * hNrm, + MPCminNvrs, ) vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA) else: @@ -595,8 +646,8 @@ def calc_vPPnext(S, a, R): vPfunc=vPfuncNow, vPPfunc=vPPfuncNow, mNrmMin=mNrmMinNow, - hNrm=hNrmNow, - MPCmin=MPCminNow, + hNrm=hNrm, + MPCmin=MPCmin, MPCmax=MPCmaxEff, ) return solution_now @@ -616,8 +667,7 @@ def solve_one_period_ConsKinkedR( vFuncBool, CubicBool, ): - """ - Solves one period of a consumption-saving model with idiosyncratic shocks to + """Solves one period of a consumption-saving model with idiosyncratic shocks to permanent and transitory income, with a risk free asset and CRRA utility. In this variation, the interest rate on borrowing Rboro exceeds the interest rate on saving Rsave. @@ -663,6 +713,7 @@ def solve_one_period_ConsKinkedR( ------- solution_now : ConsumerSolution Solution to this period's consumption-saving problem with income risk. + """ # Verifiy that there is actually a kink in the interest factor assert ( @@ -690,17 +741,8 @@ def solve_one_period_ConsKinkedR( uFunc = UtilityFuncCRRA(CRRA) DiscFacEff = DiscFac * LivPrb # "effective" discount factor - # Unpack next period's income shock distribution - ShkPrbsNext = IncShkDstn.pmv - PermShkValsNext = IncShkDstn.atoms[0] - TranShkValsNext = IncShkDstn.atoms[1] - PermShkMinNext = np.min(PermShkValsNext) - TranShkMinNext = np.min(TranShkValsNext) - # Calculate the probability that we get the worst possible income draw - IncNext = PermShkValsNext * TranShkValsNext - WorstIncNext = PermShkMinNext * TranShkMinNext - WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext]) + WorstIncPrb = calc_worst_inc_prob(IncShkDstn) # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing # Unpack next period's (marginal) value function @@ -709,22 +751,20 @@ def solve_one_period_ConsKinkedR( vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False # Update the bounding MPCs and PDV of human wealth: - PatFac = ((Rsave * DiscFacEff) ** (1.0 / CRRA)) / Rsave - PatFacAlt = ((Rboro * DiscFacEff) ** (1.0 / CRRA)) / Rboro - try: - MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) - except: - MPCminNow = 0.0 - Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext) - hNrmNow = (PermGroFac / Rsave) * (Ex_IncNext + solution_next.hNrm) - temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFacAlt - MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax) + PatFacSave = calc_pat_fac(Rsave, DiscFacEff, CRRA) + PatFacBoro = calc_pat_fac(Rboro, DiscFacEff, CRRA) + MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave) + Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) + hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext) + MPCmaxNow = calc_mpc_max(solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro) + cFuncLimitIntercept = MPCminNow * hNrmNow cFuncLimitSlope = MPCminNow # Calculate the minimum allowable value of money resources in this period - PermGroFacEffMin = (PermGroFac * PermShkMinNext) / Rboro - BoroCnstNat = (solution_next.mNrmMin - TranShkMinNext) * PermGroFacEffMin + BoroCnstNat = calc_boro_const_nat( + solution_next.mNrmMin, IncShkDstn, Rboro, PermGroFac + ) # Set the minimum allowable (normalized) market resources based on the natural # and artificial borrowing constraints @@ -742,12 +782,13 @@ def solve_one_period_ConsKinkedR( # Define the borrowing-constrained consumption function cFuncNowCnst = LinearInterp( - np.array([mNrmMinNow, mNrmMinNow + 1.0]), np.array([0.0, 1.0]) + np.array([mNrmMinNow, mNrmMinNow + 1.0]), + np.array([0.0, 1.0]), ) # Construct the assets grid by adjusting aXtra by the natural borrowing constraint aNrmNow = np.sort( - np.hstack((np.asarray(aXtraGrid) + mNrmMinNow, np.array([0.0, 0.0]))) + np.hstack((np.asarray(aXtraGrid) + mNrmMinNow, np.array([0.0, 0.0]))), ) # Make a 1D array of the interest factor at each asset gridpoint @@ -756,24 +797,13 @@ def solve_one_period_ConsKinkedR( i_kink = np.argwhere(aNrmNow == 0.0)[0][0] Rfree[i_kink] = Rboro - # Define local functions for taking future expectations - def calc_mNrmNext(S, a, R): - return R / (PermGroFac * S["PermShk"]) * a + S["TranShk"] - - def calc_vNext(S, a, R): - return (S["PermShk"] ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)) * vFuncNext( - calc_mNrmNext(S, a, R) - ) - - def calc_vPnext(S, a, R): - return S["PermShk"] ** (-CRRA) * vPfuncNext(calc_mNrmNext(S, a, R)) - - def calc_vPPnext(S, a, R): - return S["PermShk"] ** (-CRRA - 1.0) * vPPfuncNext(calc_mNrmNext(S, a, R)) - # Calculate end-of-period marginal value of assets at each gridpoint vPfacEff = DiscFacEff * Rfree * PermGroFac ** (-CRRA) - EndOfPrdvP = vPfacEff * expected(calc_vPnext, IncShkDstn, args=(aNrmNow, Rfree)) + EndOfPrdvP = vPfacEff * expected( + calc_vp_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vPfuncNext), + ) # Invert the first order condition to find optimal cNrm from each aNrm gridpoint cNrmNow = uFunc.derinv(EndOfPrdvP, order=(1, 0)) @@ -788,7 +818,9 @@ def calc_vPPnext(S, a, R): # Calculate end-of-period marginal marginal value of assets at each gridpoint vPPfacEff = DiscFacEff * Rfree * Rfree * PermGroFac ** (-CRRA - 1.0) EndOfPrdvPP = vPPfacEff * expected( - calc_vPPnext, IncShkDstn, args=(aNrmNow, Rfree) + calc_vpp_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vPPfuncNext), ) dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2) MPC = dcda / (dcda + 1.0) @@ -834,9 +866,13 @@ def calc_vPPnext(S, a, R): # Construct this period's value function if requested if vFuncBool: # Calculate end-of-period value, its derivative, and their pseudo-inverse - EndOfPrdv = DiscFacEff * expected(calc_vNext, IncShkDstn, args=(aNrmNow, Rfree)) + EndOfPrdv = DiscFacEff * expected( + calc_v_next, + IncShkDstn, + args=(aNrmNow, Rfree, CRRA, PermGroFac, vFuncNext), + ) EndOfPrdvNvrs = uFunc.inv( - EndOfPrdv + EndOfPrdv, ) # value transformed through inverse utility EndOfPrdvNvrsP = EndOfPrdvP * uFunc.derinv(EndOfPrdv, order=(0, 1)) EndOfPrdvNvrs = np.insert(EndOfPrdvNvrs, 0, 0.0) @@ -863,7 +899,11 @@ def calc_vPPnext(S, a, R): vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA))) MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA)) vNvrsFuncNow = CubicInterp( - mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs + mNrm_temp, + vNvrs_temp, + vNvrsP_temp, + MPCminNvrs * hNrmNow, + MPCminNvrs, ) vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA) else: From 6f9187d0bfc98f5eb75ba24b1502ae8f9422ca0a Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Thu, 14 Mar 2024 13:19:26 -0400 Subject: [PATCH 02/14] take subfunctions out --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 328 +++++++++++-------- 1 file changed, 197 insertions(+), 131 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index c3fb4f5a7..4b9414577 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -329,6 +329,159 @@ def __init__(self, verbose=False, quiet=False, **kwds): self.solve_one_period = make_one_period_oo_solver(ConsSequentialPortfolioSolver) +def calc_m_nrm_next(shocks, b_nrm, perm_gro_fac): + """Calculate future realizations of market resources mNrm from the income + shock distribution S and normalized bank balances b. + """ + return b_nrm / (shocks["PermShk"] * perm_gro_fac) + shocks["TranShk"] + + +def calc_dvdm_next( + shocks, + b_nrm, + share, + crra, + perm_gro_fac, + adjust_prob, + vp_func_adj, + dvdm_func_fxd, +): + """Evaluate realizations of marginal value of market resources next period, + based on the income distribution S, values of bank balances bNrm, and + values of the risky share z. + """ + mNrm_next = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + dvdmAdj_next = vp_func_adj(mNrm_next) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + Share_next_expanded = share + np.zeros_like(mNrm_next) + dvdmFxd_next = dvdm_func_fxd(mNrm_next, Share_next_expanded) + # Combine by adjustment probability + dvdm_next = adjust_prob * dvdmAdj_next + (1.0 - adjust_prob) * dvdmFxd_next + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvdm_next = dvdmAdj_next + + dvdm_next = (shocks["PermShk"] * perm_gro_fac) ** (-crra) * dvdm_next + return dvdm_next + + +def calc_dvds_next( + shocks, + b_nrm, + share, + crra, + perm_gro_fac, + adjust_prob, + dvds_func_fxd, +): + """Evaluate realizations of marginal value of risky share next period, based + on the income distribution S, values of bank balances bNrm, and values of + the risky share z. + """ + mNrm_next = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + + # No marginal value of Share if it's a free choice! + dvdsAdj_next = np.zeros_like(mNrm_next) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + Share_next_expanded = share + np.zeros_like(mNrm_next) + dvdsFxd_next = dvds_func_fxd(mNrm_next, Share_next_expanded) + # Combine by adjustment probability + dvds_next = adjust_prob * dvdsAdj_next + (1.0 - adjust_prob) * dvdsFxd_next + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvds_next = dvdsAdj_next + + dvds_next = (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * dvds_next + return dvds_next + + +# Define functions for calculating end-of-period marginal value +def calc_end_of_prd_dvda(shocks, a_nrm, share, rfree, dvdb_func): + """Compute end-of-period marginal value of assets at values a, conditional + on risky asset return S and risky share z. + """ + # Calculate future realizations of bank balances bNrm + Rxs = shocks - rfree # Excess returns + Rport = rfree + share * Rxs # Portfolio return + bNrm_next = Rport * a_nrm + + # Ensure shape concordance + z_rep = share + np.zeros_like(bNrm_next) + + # Calculate and return dvda + EndOfPrd_dvda = Rport * dvdb_func(bNrm_next, z_rep) + return EndOfPrd_dvda + + +def calc_end_of_prd_dvds( + shocks, + a_nrm, + share, + rfree, + dvdb_func, + dvds_func, +): + """Compute end-of-period marginal value of risky share at values a, conditional + on risky asset return S and risky share z. + """ + # Calculate future realizations of bank balances bNrm + Rxs = shocks - rfree # Excess returns + Rport = rfree + share * Rxs # Portfolio return + bNrm_next = Rport * a_nrm + + # Make the shares match the dimension of b, so that it can be vectorized + z_rep = share + np.zeros_like(bNrm_next) + + # Calculate and return dvds + EndOfPrd_dvds = Rxs * a_nrm * dvdb_func(bNrm_next, z_rep) + dvds_func( + bNrm_next, z_rep + ) + return EndOfPrd_dvds + + +def calc_v_intermed( + shocks, + b_nrm, + share, + crra, + perm_gro_fac, + adjust_prob, + v_func_adj, + v_func_fxd, +): + """Calculate "intermediate" value from next period's bank balances, the + income shocks S, and the risky asset share. + """ + mNrm_next = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + + vAdj_next = v_func_adj(mNrm_next) + if adjust_prob < 1.0: + vFxd_next = v_func_fxd(mNrm_next, share) + # Combine by adjustment probability + v_next = adjust_prob * vAdj_next + (1.0 - adjust_prob) * vFxd_next + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + v_next = vAdj_next + + v_intermed = (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * v_next + return v_intermed + + +def calc_end_of_prd_v(shocks, a_nrm, share, rfree, v_func): + # Calculate future realizations of bank balances bNrm + Rxs = shocks - rfree + Rport = rfree + share * Rxs + bNrm_next = Rport * a_nrm + + # Make an extended share_next of the same dimension as b_nrm so + # that the function can be vectorized + z_rep = share + np.zeros_like(bNrm_next) + + EndOfPrd_v = v_func(bNrm_next, z_rep) + return EndOfPrd_v + + def solve_one_period_ConsPortfolio( solution_next, ShockDstn, @@ -348,8 +501,7 @@ def solve_one_period_ConsPortfolio( DiscreteShareBool, IndepDstnBool, ): - """ - Solve one period of a consumption-saving problem with portfolio allocation + """Solve one period of a consumption-saving problem with portfolio allocation between a riskless and risky asset. This function handles various sub-cases or variations on the problem, including the possibility that the agent does not necessarily get to update their portfolio share in every period, or that @@ -413,6 +565,7 @@ def solve_one_period_ConsPortfolio( ------- solution_now : PortfolioSolution Solution to this period's problem. + """ # Make sure the individual is liquidity constrained. Allowing a consumer to # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix. @@ -423,7 +576,7 @@ def solve_one_period_ConsPortfolio( # the value function is also constructed (else this task would be impossible). if DiscreteShareBool and (not vFuncBool): raise ValueError( - "PortfolioConsumerType requires vFuncBool to be True when DiscreteShareBool is True!" + "PortfolioConsumerType requires vFuncBool to be True when DiscreteShareBool is True!", ) # Define the current period utility function and effective discount factor @@ -468,170 +621,83 @@ def solve_one_period_ConsPortfolio( # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn bNrmNext, ShareNext = np.meshgrid(bNrmGrid, ShareGrid, indexing="ij") - # Define functions that are used internally to evaluate future realizations - def calc_mNrm_next(S, b): - """ - Calculate future realizations of market resources mNrm from the income - shock distribution S and normalized bank balances b. - """ - return b / (S["PermShk"] * PermGroFac) + S["TranShk"] - - def calc_dvdm_next(S, b, z): - """ - Evaluate realizations of marginal value of market resources next period, - based on the income distribution S, values of bank balances bNrm, and - values of the risky share z. - """ - mNrm_next = calc_mNrm_next(S, b) - dvdmAdj_next = vPfuncAdj_next(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - dvdmFxd_next = dvdmFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - dvdm_next = AdjustPrb * dvdmAdj_next + (1.0 - AdjustPrb) * dvdmFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - dvdm_next = dvdmAdj_next - - dvdm_next = (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next - return dvdm_next - - def calc_dvds_next(S, b, z): - """ - Evaluate realizations of marginal value of risky share next period, based - on the income distribution S, values of bank balances bNrm, and values of - the risky share z. - """ - mNrm_next = calc_mNrm_next(S, b) - - # No marginal value of Share if it's a free choice! - dvdsAdj_next = np.zeros_like(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - dvdsFxd_next = dvdsFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - dvds_next = AdjustPrb * dvdsAdj_next + (1.0 - AdjustPrb) * dvdsFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - dvds_next = dvdsAdj_next - - dvds_next = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * dvds_next - return dvds_next - # Calculate end-of-period marginal value of assets and shares at each point # in aNrm and ShareGrid. Does so by taking expectation of next period marginal # values across income and risky return shocks. # Calculate intermediate marginal value of bank balances by taking expectations over income shocks - dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext, ShareNext)) + dvdb_intermed = expected( + calc_dvdm_next, + IncShkDstn, + args=( + bNrmNext, + ShareNext, + CRRA, + PermGroFac, + AdjustPrb, + vPfuncAdj_next, + dvdmFuncFxd_next, + ), + ) dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0)) dvdbNvrsFunc_intermed = BilinearInterp(dvdbNvrs_intermed, bNrmGrid, ShareGrid) dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA) # Calculate intermediate marginal value of risky portfolio share by taking expectations over income shocks - dvds_intermed = expected(calc_dvds_next, IncShkDstn, args=(bNrmNext, ShareNext)) + dvds_intermed = expected( + calc_dvds_next, + IncShkDstn, + args=(bNrmNext, ShareNext, CRRA, PermGroFac, AdjustPrb, dvdsFuncFxd_next), + ) dvdsFunc_intermed = BilinearInterp(dvds_intermed, bNrmGrid, ShareGrid) # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij") - # Define functions for calculating end-of-period marginal value - def calc_EndOfPrd_dvda(S, a, z): - """ - Compute end-of-period marginal value of assets at values a, conditional - on risky asset return S and risky share z. - """ - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree # Excess returns - Rport = Rfree + z * Rxs # Portfolio return - bNrm_next = Rport * a - - # Ensure shape concordance - z_rep = z + np.zeros_like(bNrm_next) - - # Calculate and return dvda - EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_dvda - - def EndOfPrddvds_dist(S, a, z): - """ - Compute end-of-period marginal value of risky share at values a, conditional - on risky asset return S and risky share z. - """ - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree # Excess returns - Rport = Rfree + z * Rxs # Portfolio return - bNrm_next = Rport * a - - # Make the shares match the dimension of b, so that it can be vectorized - z_rep = z + np.zeros_like(bNrm_next) - - # Calculate and return dvds - EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed( - bNrm_next, z_rep - ) + dvdsFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_dvds - # Evaluate realizations of value and marginal value after asset returns are realized # Calculate end-of-period marginal value of assets by taking expectations EndOfPrd_dvda = DiscFacEff * expected( - calc_EndOfPrd_dvda, RiskyDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_dvda, + RiskyDstn, + args=(aNrmNow, ShareNext, Rfree, dvdbFunc_intermed), ) EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) # Calculate end-of-period marginal value of risky portfolio share by taking expectations EndOfPrd_dvds = DiscFacEff * expected( - EndOfPrddvds_dist, RiskyDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_dvds, + RiskyDstn, + args=(aNrmNow, ShareNext, Rfree, dvdbFunc_intermed, dvdsFunc_intermed), ) # Make the end-of-period value function if the value function is requested if vFuncBool: - - def calc_v_intermed(S, b, z): - """ - Calculate "intermediate" value from next period's bank balances, the - income shocks S, and the risky asset share. - """ - mNrm_next = calc_mNrm_next(S, b) - - vAdj_next = vFuncAdj_next(mNrm_next) - if AdjustPrb < 1.0: - vFxd_next = vFuncFxd_next(mNrm_next, z) - # Combine by adjustment probability - v_next = AdjustPrb * vAdj_next + (1.0 - AdjustPrb) * vFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - v_next = vAdj_next - - v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next - return v_intermed - # Calculate intermediate value by taking expectations over income shocks - v_intermed = expected(calc_v_intermed, IncShkDstn, args=(bNrmNext, ShareNext)) + v_intermed = expected( + calc_v_intermed, + IncShkDstn, + args=( + bNrmNext, + ShareNext, + CRRA, + PermGroFac, + AdjustPrb, + vFuncAdj_next, + vFuncFxd_next, + ), + ) # Construct the "intermediate value function" for this period vNvrs_intermed = uFunc.inv(v_intermed) vNvrsFunc_intermed = BilinearInterp(vNvrs_intermed, bNrmGrid, ShareGrid) vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA) - def calc_EndOfPrd_v(S, a, z): - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree - Rport = Rfree + z * Rxs - bNrm_next = Rport * a - - # Make an extended share_next of the same dimension as b_nrm so - # that the function can be vectorized - z_rep = z + np.zeros_like(bNrm_next) - - EndOfPrd_v = vFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_v - # Calculate end-of-period value by taking expectations EndOfPrd_v = DiscFacEff * expected( - calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_v, + RiskyDstn, + args=(aNrmNow, ShareNext, Rfree, vFunc_intermed), ) EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v) @@ -797,7 +863,7 @@ def calc_EndOfPrd_v(S, a, z): np.insert(mNrm_temp[:, 0], 0, 0.0), # x_list np.insert(vNvrs_temp[:, j], 0, 0.0), # f_list np.insert(vNvrsP_temp[:, j], 0, vNvrsP_temp[j, 0]), # dfdx_list - ) + ), ) vNvrsFuncFxd = LinearInterpOnInterp1D(vNvrsFuncFxd_by_Share, ShareGrid) vFuncFxd_now = ValueFuncCRRA(vNvrsFuncFxd, CRRA) From c0af9f4a14fb3475db921fcbdaabea333a717d65 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Thu, 14 Mar 2024 13:48:34 -0400 Subject: [PATCH 03/14] merge expectation functions to make more efficient --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 89 +++++++++++++++----- 1 file changed, 69 insertions(+), 20 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index 4b9414577..a444e925c 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -397,6 +397,45 @@ def calc_dvds_next( return dvds_next +def calc_dv_next( + shocks, + b_nrm, + share, + crra, + perm_gro_fac, + adjust_prob, + vp_func_adj, + dvdm_func_fxd, + dvds_func_fxd, +): + mNrm_next = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + dvdmAdj_next = vp_func_adj(mNrm_next) + + # No marginal value of Share if it's a free choice! + dvdsAdj_next = np.zeros_like(mNrm_next) + + # If there's no chance that portfolio share is fixed, don't bother evaluating + if adjust_prob >= 1.0: + return (shocks["PermShk"] * perm_gro_fac) ** ( + -crra + ) * dvdmAdj_next, dvdsAdj_next + + # Expand to the same dimensions as mNrm + Share_next_expanded = share + np.zeros_like(mNrm_next) + + dvdmFxd_next = dvdm_func_fxd(mNrm_next, Share_next_expanded) + dvdsFxd_next = dvds_func_fxd(mNrm_next, Share_next_expanded) + + # Combine by adjustment probability + dvdm_next = adjust_prob * dvdmAdj_next + (1.0 - adjust_prob) * dvdmFxd_next + dvds_next = adjust_prob * dvdsAdj_next + (1.0 - adjust_prob) * dvdsFxd_next + + dvdm_next = (shocks["PermShk"] * perm_gro_fac) ** (-crra) * dvdm_next + dvds_next = (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * dvds_next + + return dvdm_next, dvds_next + + # Define functions for calculating end-of-period marginal value def calc_end_of_prd_dvda(shocks, a_nrm, share, rfree, dvdb_func): """Compute end-of-period marginal value of assets at values a, conditional @@ -441,6 +480,24 @@ def calc_end_of_prd_dvds( return EndOfPrd_dvds +def calc_end_of_prd_dv(shocks, a_nrm, share, rfree, dvdb_func, dvds_func): + # Calculate future realizations of bank balances bNrm + Rxs = shocks - rfree # Excess returns + Rport = rfree + share * Rxs # Portfolio return + bNrm_next = Rport * a_nrm + + # Ensure shape concordance + z_rep = np.full_like(bNrm_next, share) + + # Calculate and return dvda + dvdb = dvdb_func(bNrm_next, z_rep) + EndOfPrd_dvda = Rport * dvdb + # Calculate and return dvds + EndOfPrd_dvds = Rxs * a_nrm * dvdb + dvds_func(bNrm_next, z_rep) + + return EndOfPrd_dvda, EndOfPrd_dvds + + def calc_v_intermed( shocks, b_nrm, @@ -625,9 +682,10 @@ def solve_one_period_ConsPortfolio( # in aNrm and ShareGrid. Does so by taking expectation of next period marginal # values across income and risky return shocks. - # Calculate intermediate marginal value of bank balances by taking expectations over income shocks - dvdb_intermed = expected( - calc_dvdm_next, + # Calculate intermediate marginal value of bank balances and risky portfolio share + # by taking expectations over income shocks + dvdb_intermed, dvds_intermed = expected( + calc_dv_next, IncShkDstn, args=( bNrmNext, @@ -637,18 +695,14 @@ def solve_one_period_ConsPortfolio( AdjustPrb, vPfuncAdj_next, dvdmFuncFxd_next, + dvdsFuncFxd_next, ), ) + dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0)) dvdbNvrsFunc_intermed = BilinearInterp(dvdbNvrs_intermed, bNrmGrid, ShareGrid) dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA) - # Calculate intermediate marginal value of risky portfolio share by taking expectations over income shocks - dvds_intermed = expected( - calc_dvds_next, - IncShkDstn, - args=(bNrmNext, ShareNext, CRRA, PermGroFac, AdjustPrb, dvdsFuncFxd_next), - ) dvdsFunc_intermed = BilinearInterp(dvds_intermed, bNrmGrid, ShareGrid) # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn @@ -656,21 +710,16 @@ def solve_one_period_ConsPortfolio( # Evaluate realizations of value and marginal value after asset returns are realized - # Calculate end-of-period marginal value of assets by taking expectations - EndOfPrd_dvda = DiscFacEff * expected( - calc_end_of_prd_dvda, - RiskyDstn, - args=(aNrmNow, ShareNext, Rfree, dvdbFunc_intermed), - ) - EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) - - # Calculate end-of-period marginal value of risky portfolio share by taking expectations - EndOfPrd_dvds = DiscFacEff * expected( - calc_end_of_prd_dvds, + # Calculate end-of-period marginal value of assets and risky portfolio share + # by taking expectations + EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected( + calc_end_of_prd_dv, RiskyDstn, args=(aNrmNow, ShareNext, Rfree, dvdbFunc_intermed, dvdsFunc_intermed), ) + EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) + # Make the end-of-period value function if the value function is requested if vFuncBool: # Calculate intermediate value by taking expectations over income shocks From 01745f8c675586dc1984c8371e603e21e81ae4ad Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 20 Mar 2024 11:25:03 -0400 Subject: [PATCH 04/14] all params in bequest are scalars --- HARK/ConsumptionSaving/ConsBequestModel.py | 39 ++++------------------ 1 file changed, 6 insertions(+), 33 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsBequestModel.py b/HARK/ConsumptionSaving/ConsBequestModel.py index 3553ef3ce..9d615dfa8 100644 --- a/HARK/ConsumptionSaving/ConsBequestModel.py +++ b/HARK/ConsumptionSaving/ConsBequestModel.py @@ -43,14 +43,7 @@ class BequestWarmGlowConsumerType(IndShockConsumerType): - time_inv_ = IndShockConsumerType.time_inv_ + [ - "BeqCRRA", - "BeqShift", - ] - - time_vary_ = IndShockConsumerType.time_vary_ + [ - "BeqFac", - ] + time_inv_ = IndShockConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"] def __init__(self, **kwds): params = init_wealth_in_utility.copy() @@ -68,14 +61,8 @@ def update_parameters(self): if not isinstance(self.BeqCRRA, (int, float)): raise ValueError("Bequest CRRA parameter must be a single value.") - if isinstance(self.BeqFac, (int, float)): - self.BeqFac = [self.BeqFac] * self.T_cycle - elif len(self.BeqFac) == 1: - self.BeqFac *= self.T_cycle - elif len(self.BeqFac) != self.T_cycle: - raise ValueError( - "Bequest relative value parameter must be a single value or a list of length T_cycle", - ) + if not isinstance(self.BeqFac, (int, float)): + raise ValueError("Bequest relative value parameter must be a single value.") if not isinstance(self.BeqShift, (int, float)): raise ValueError("Bequest Stone-Geary parameter must be a single value.") @@ -117,15 +104,7 @@ def update_solution_terminal(self): class BequestWarmGlowPortfolioType(PortfolioConsumerType): - time_inv_ = IndShockConsumerType.time_inv_ + [ - "BeqCRRA", - "BeqShift", - "DiscreteShareBool", - ] - - time_vary_ = IndShockConsumerType.time_vary_ + [ - "BeqFac", - ] + time_inv_ = IndShockConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"] def __init__(self, **kwds): params = init_portfolio_bequest.copy() @@ -145,14 +124,8 @@ def update_parameters(self): if not isinstance(self.BeqCRRA, (int, float)): raise ValueError("Bequest CRRA parameter must be a single value.") - if isinstance(self.BeqFac, (int, float)): - self.BeqFac = [self.BeqFac] * self.T_cycle - elif len(self.BeqFac) == 1: - self.BeqFac *= self.T_cycle - elif len(self.BeqFac) != self.T_cycle: - raise ValueError( - "Bequest relative value parameter must be a single value or a list of length T_cycle", - ) + if not isinstance(self.BeqFac, (int, float)): + raise ValueError("Bequest relative value parameter must be a single value.") if not isinstance(self.BeqShift, (int, float)): raise ValueError("Bequest Stone-Geary parameter must be a single value.") From 8f7b3b32348860ec7e234e6c43c109224b117aba Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 27 Mar 2024 09:11:13 -0400 Subject: [PATCH 05/14] fix failing test --- HARK/ConsumptionSaving/ConsBequestModel.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/HARK/ConsumptionSaving/ConsBequestModel.py b/HARK/ConsumptionSaving/ConsBequestModel.py index 4077cffb3..33b080777 100644 --- a/HARK/ConsumptionSaving/ConsBequestModel.py +++ b/HARK/ConsumptionSaving/ConsBequestModel.py @@ -102,7 +102,12 @@ def update_solution_terminal(self): class BequestWarmGlowPortfolioType(PortfolioConsumerType): - time_inv_ = IndShockConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"] + time_inv_ = IndShockConsumerType.time_inv_ + [ + "BeqCRRA", + "BeqShift", + "BeqFac", + "DiscreteShareBool", + ] def __init__(self, **kwds): params = init_portfolio_bequest.copy() From 3bab06261357da5e3dc9943396941a14ef1d6a5f Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 27 Mar 2024 14:25:26 -0400 Subject: [PATCH 06/14] modify calc_mpc_max --- HARK/ConsumptionSaving/ConsIndShockModel.py | 141 ++++++++++---------- 1 file changed, 67 insertions(+), 74 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 06b88c123..ef9fe98f3 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -270,7 +270,7 @@ def solve_one_period_ConsPF( # Calculate the lower bound of the marginal propensity to consume PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) - MPCmin = calc_mpc_min(solution_next.MPCmin, PatFac) + MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac) # Extract the discrete kink points in next period's consumption function; # don't take the last one, as it only defines the extrapolation and is not a kink. @@ -289,12 +289,12 @@ def solve_one_period_ConsPF( # Calculate (pseudo-inverse) value at each consumption kink point vNow = uFunc(cNrmNow) + EndOfPrdv vNvrsNow = uFunc.inverse(vNow) - vNvrsSlopeMin = MPCmin ** (-CRRA / (1.0 - CRRA)) + vNvrsSlopeMin = MPCminNow ** (-CRRA / (1.0 - CRRA)) # Add an additional point to the list of gridpoints for the extrapolation, # using the new value of the lower bound of the MPC. mNrmNow = np.append(mNrmNow, mNrmNow[-1] + 1.0) - cNrmNow = np.append(cNrmNow, cNrmNow[-1] + MPCmin) + cNrmNow = np.append(cNrmNow, cNrmNow[-1] + MPCminNow) vNvrsNow = np.append(vNvrsNow, vNvrsNow[-1] + vNvrsSlopeMin) # If the artificial borrowing constraint binds, combine the constrained and @@ -330,11 +330,11 @@ def solve_one_period_ConsPF( # If it *is* the very last index, then there are only three points # that characterize the consumption function: the artificial borrowing # constraint, the constraint kink, and the extrapolation point. - mXtra = (cNrmNow[-1] - cNrmCnst[-1]) / (1.0 - MPCmin) + mXtra = (cNrmNow[-1] - cNrmCnst[-1]) / (1.0 - MPCminNow) mCrit = mNrmNow[-1] + mXtra cCrit = mCrit - BoroCnstArt mNrmNow = np.array([BoroCnstArt, mCrit, mCrit + 1.0]) - cNrmNow = np.array([0.0, cCrit, cCrit + MPCmin]) + cNrmNow = np.array([0.0, cCrit, cCrit + MPCminNow]) # Adjust vNvrs grid for this three node structure mNextCrit = BoroCnstArt * Rfree + 1.0 @@ -347,31 +347,31 @@ def solve_one_period_ConsPF( # kink point, being sure to adjust the extrapolation. if mNrmNow.size > MaxKinks: mNrmNow = np.concatenate((mNrmNow[:-2], [mNrmNow[-3] + 1.0])) - cNrmNow = np.concatenate((cNrmNow[:-2], [cNrmNow[-3] + MPCmin])) + cNrmNow = np.concatenate((cNrmNow[:-2], [cNrmNow[-3] + MPCminNow])) vNvrsNow = np.concatenate((vNvrsNow[:-2], [vNvrsNow[-3] + vNvrsSlopeMin])) # Construct the consumption function as a linear interpolation. - cFunc = LinearInterp(mNrmNow, cNrmNow) + cFuncNow = LinearInterp(mNrmNow, cNrmNow) # Calculate the upper bound of the MPC as the slope of the bottom segment. - MPCmax = (cNrmNow[1] - cNrmNow[0]) / (mNrmNow[1] - mNrmNow[0]) + MPCmaxNow = (cNrmNow[1] - cNrmNow[0]) / (mNrmNow[1] - mNrmNow[0]) mNrmMinNow = mNrmNow[0] # Construct the (marginal) value function for this period # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations vFuncNvrs = LinearInterp(mNrmNow, vNvrsNow) - vFunc = ValueFuncCRRA(vFuncNvrs, CRRA) - vPfunc = MargValueFuncCRRA(cFunc, CRRA) + vFuncNow = ValueFuncCRRA(vFuncNvrs, CRRA) + vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA) # Construct and return the solution solution_now = ConsumerSolution( - cFunc=cFunc, - vFunc=vFunc, - vPfunc=vPfunc, + cFunc=cFuncNow, + vFunc=vFuncNow, + vPfunc=vPfuncNow, mNrmMin=mNrmMinNow, hNrm=hNrmNow, - MPCmin=MPCmin, - MPCmax=MPCmax, + MPCmin=MPCminNow, + MPCmax=MPCmaxNow, ) return solution_now @@ -384,22 +384,29 @@ def calc_worst_inc_prob(inc_shk_dstn): return np.sum(probs[income == worst_inc]) -def calc_mpc_max(mpc_max_next, worst_inc_prob, crra, pat_fac): - temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac - return 1.0 / (1.0 + temp_fac / mpc_max_next) - - -def calc_boro_const_nat( - m_nrm_min_next, - inc_shk_dstn, - rfree, - perm_gro_fac, -): +def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): perm, tran = inc_shk_dstn.atoms temp_fac = (perm_gro_fac * np.min(perm)) / rfree return (m_nrm_min_next - np.min(tran)) * temp_fac +def calc_m_nrm_min(boro_const_art, boro_const_nat): + return ( + boro_const_nat + if boro_const_art is None + else max(boro_const_nat, boro_const_art) + ) + + +def calc_mpc_max( + mpc_max_next, worst_inc_prob, crra, pat_fac, boro_const_nat, boro_const_art +): + m_nrm_min = calc_m_nrm_min(boro_const_art, boro_const_nat) + temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac + mpc_max = 1.0 / (1.0 + temp_fac / mpc_max_next) + return 1.0 if boro_const_nat < m_nrm_min else mpc_max + + def calc_m_nrm_next(shock, a, rfree, perm_gro_fac): return rfree / (perm_gro_fac * shock["PermShk"]) * a + shock["TranShk"] @@ -482,40 +489,33 @@ def solve_one_period_ConsIndShock( # Calculate the probability that we get the worst possible income draw WorstIncPrb = calc_worst_inc_prob(IncShkDstn) + Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) + hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, Ex_IncNext) # Unpack next period's (marginal) value function vFuncNext = solution_next.vFunc # This is None when vFuncBool is False vPfuncNext = solution_next.vPfunc vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False - # Update the bounding MPCs and PDV of human wealth: - PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) - MPCmin = calc_mpc_min(solution_next.MPCmin, PatFac) - Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) - hNrm = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, Ex_IncNext) - MPCmax = calc_mpc_max(solution_next.MPCmax, WorstIncPrb, CRRA, PatFac) - - cFuncLimitIntercept = MPCmin * hNrm - cFuncLimitSlope = MPCmin - # Calculate the minimum allowable value of money resources in this period BoroCnstNat = calc_boro_const_nat( solution_next.mNrmMin, IncShkDstn, Rfree, PermGroFac ) - # Set the minimum allowable (normalized) market resources based on the natural # and artificial borrowing constraints - if BoroCnstArt is None: - mNrmMinNow = BoroCnstNat - else: - mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt]) + mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat) + # Update the bounding MPCs and PDV of human wealth: + PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) + MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds - if BoroCnstNat < mNrmMinNow: - MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1 - else: - MPCmaxEff = MPCmax # Otherwise, it's the MPC calculated above + MPCmaxNow = calc_mpc_max( + solution_next.MPCmax, WorstIncPrb, CRRA, PatFac, BoroCnstNat, BoroCnstArt + ) + + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow # Define the borrowing-constrained consumption function cFuncNowCnst = LinearInterp( @@ -553,7 +553,7 @@ def solve_one_period_ConsIndShock( ) dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2) MPC = dcda / (dcda + 1.0) - MPC_for_interpolation = np.insert(MPC, 0, MPCmax) + MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( @@ -618,13 +618,13 @@ def solve_one_period_ConsIndShock( vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1)) mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow) vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0) - vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA))) - MPCminNvrs = MPCmin ** (-CRRA / (1.0 - CRRA)) + vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA))) + MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA)) vNvrsFuncNow = CubicInterp( mNrm_temp, vNvrs_temp, vNvrsP_temp, - MPCminNvrs * hNrm, + MPCminNvrs * hNrmNow, MPCminNvrs, ) vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA) @@ -638,9 +638,9 @@ def solve_one_period_ConsIndShock( vPfunc=vPfuncNow, vPPfunc=vPPfuncNow, mNrmMin=mNrmMinNow, - hNrm=hNrm, - MPCmin=MPCmin, - MPCmax=MPCmaxEff, + hNrm=hNrmNow, + MPCmin=MPCminNow, + MPCmax=MPCmaxNow, ) return solution_now @@ -736,41 +736,34 @@ def solve_one_period_ConsKinkedR( # Calculate the probability that we get the worst possible income draw WorstIncPrb = calc_worst_inc_prob(IncShkDstn) # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing + Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) + hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext) # Unpack next period's (marginal) value function vFuncNext = solution_next.vFunc # This is None when vFuncBool is False vPfuncNext = solution_next.vPfunc vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False - # Update the bounding MPCs and PDV of human wealth: - PatFacSave = calc_pat_fac(Rsave, DiscFacEff, CRRA) - PatFacBoro = calc_pat_fac(Rboro, DiscFacEff, CRRA) - MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave) - Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) - hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext) - MPCmaxNow = calc_mpc_max(solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro) - - cFuncLimitIntercept = MPCminNow * hNrmNow - cFuncLimitSlope = MPCminNow - # Calculate the minimum allowable value of money resources in this period BoroCnstNat = calc_boro_const_nat( solution_next.mNrmMin, IncShkDstn, Rboro, PermGroFac ) - # Set the minimum allowable (normalized) market resources based on the natural # and artificial borrowing constraints - if BoroCnstArt is None: - mNrmMinNow = BoroCnstNat - else: - mNrmMinNow = np.max([BoroCnstNat, BoroCnstArt]) + mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat) + # Update the bounding MPCs and PDV of human wealth: + PatFacSave = calc_pat_fac(Rsave, DiscFacEff, CRRA) + PatFacBoro = calc_pat_fac(Rboro, DiscFacEff, CRRA) + MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds - if BoroCnstNat < mNrmMinNow: - MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1 - else: - MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above + MPCmaxNow = calc_mpc_max( + solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro, BoroCnstNat, BoroCnstArt + ) + + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow # Define the borrowing-constrained consumption function cFuncNowCnst = LinearInterp( @@ -888,7 +881,7 @@ def solve_one_period_ConsKinkedR( vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1)) mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow) vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0) - vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA))) + vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA))) MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA)) vNvrsFuncNow = CubicInterp( mNrm_temp, @@ -910,7 +903,7 @@ def solve_one_period_ConsKinkedR( mNrmMin=mNrmMinNow, hNrm=hNrmNow, MPCmin=MPCminNow, - MPCmax=MPCmaxEff, + MPCmax=MPCmaxNow, ) return solution_now From 4de92547dd27a3cc401a6a5844b979cbc280d704 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 27 Mar 2024 16:26:19 -0400 Subject: [PATCH 07/14] pull out helper functions --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 534 +++++++++++-------- 1 file changed, 319 insertions(+), 215 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index d4ecb1b0b..80e015a7b 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -307,6 +307,260 @@ def get_controls(self): self.controls["Share"] = ShareNow +def calc_m_nrm_next(shocks, b_nrm, perm_gro_fac): + """ + Calculate future realizations of market resources mNrm from the income + shock distribution "shocks" and normalized bank balances b. + """ + return b_nrm / (shocks["PermShk"] * perm_gro_fac) + shocks["TranShk"] + + +def calc_dvdm_next( + shocks, b_nrm, share, adjust_prob, perm_gro_fac, crra, vp_func_adj, dvdm_func_fxd +): + """ + Evaluate realizations of marginal value of market resources next period, + based on the income distribution "shocks", values of bank balances bNrm, and + values of the risky share z. + """ + m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + dvdm_adj = vp_func_adj(m_nrm) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + share_exp = np.full_like(m_nrm, share) + dvdm_fxd = dvdm_func_fxd(m_nrm, share_exp) + # Combine by adjustment probability + dvdm_next = adjust_prob * dvdm_adj + (1.0 - adjust_prob) * dvdm_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvdm_next = dvdm_adj + + dvdm_next = (shocks["PermShk"] * perm_gro_fac) ** (-crra) * dvdm_next + return dvdm_next + + +def calc_dvds_next( + shocks, b_nrm, share, adjust_prob, perm_gro_fac, crra, dvds_func_fxd +): + """ + Evaluate realizations of marginal value of risky share next period, based + on the income distribution "shocks", values of bank balances bNrm, and values of + the risky share z. + """ + m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + + # No marginal value of shockshare if it's a free choice! + dvds_adj = np.zeros_like(m_nrm) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + share_exp = np.full_like(m_nrm, share) + dvds_fxd = dvds_func_fxd(m_nrm, share_exp) + # Combine by adjustment probability + dvds_next = adjust_prob * dvds_adj + (1.0 - adjust_prob) * dvds_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvds_next = dvds_adj + + dvds_next = (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * dvds_next + return dvds_next + + +def calc_dvdx_next( + shocks, + b_nrm, + share, + adjust_prob, + perm_gro_fac, + crra, + vp_func_adj, + dvdm_func_fxd, + dvds_func_fxd, +): + m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + dvdm_adj = vp_func_adj(m_nrm) + # No marginal value of shockshare if it's a free choice! + dvds_adj = np.zeros_like(m_nrm) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + share_exp = np.full_like(m_nrm, share) + dvdm_fxd = dvdm_func_fxd(m_nrm, share_exp) + dvds_fxd = dvds_func_fxd(m_nrm, share_exp) + # Combine by adjustment probability + dvdm = adjust_prob * dvdm_adj + (1.0 - adjust_prob) * dvdm_fxd + dvds = adjust_prob * dvds_adj + (1.0 - adjust_prob) * dvds_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvdm = dvdm_adj + dvds = dvds_adj + + perm_shk_fac = shocks["PermShk"] * perm_gro_fac + dvdm = perm_shk_fac ** (-crra) * dvdm + dvds = perm_shk_fac ** (1.0 - crra) * dvds + + return dvdm, dvds + + +def calc_end_of_prd_dvda(shocks, a_nrm, share, rfree, dvdb_func): + """ + Compute end-of-period marginal value of assets at values a, conditional + on risky asset return shocks and risky share z. + """ + # Calculate future realizations of bank balances bNrm + ex_ret = shocks - rfree # Excess returns + r_port = rfree + share * ex_ret # Portfolio return + b_nrm = r_port * a_nrm + + # Ensure shape concordance + share_exp = np.full_like(b_nrm, share) + + # Calculate and return dvda + return r_port * dvdb_func(b_nrm, share_exp) + + +def calc_end_of_prd_dvds(shocks, a_nrm, share, rfree, dvdb_func, dvds_func): + """ + Compute end-of-period marginal value of risky share at values a, conditional + on risky asset return shocks and risky share z. + """ + # Calculate future realizations of bank balances bNrm + ex_ret = shocks - rfree # Excess returns + r_port = rfree + share * ex_ret # Portfolio return + b_nrm = r_port * a_nrm + + # Make the shares match the dimension of b, so that it can be vectorized + share_exp = np.full_like(b_nrm, share) + + # Calculate and return dvds + + return ex_ret * a_nrm * dvdb_func(b_nrm, share_exp) + dvds_func(b_nrm, share_exp) + + +def calc_end_of_prd_dvdx(shocks, a_nrm, share, rfree, dvdb_func, dvds_func): + # Calculate future realizations of bank balances bNrm + ex_ret = shocks - rfree # Excess returns + r_port = rfree + share * ex_ret # Portfolio return + b_nrm = r_port * a_nrm + # Ensure shape concordance + share_exp = np.full_like(b_nrm, share) + + # Calculate and return dvda, dvds + dvda = r_port * dvdb_func(b_nrm, share_exp) + dvds = ex_ret * a_nrm * dvdb_func(b_nrm, share_exp) + dvds_func(b_nrm, share_exp) + return dvda, dvds + + +def calc_v_intermed( + shocks, b_nrm, share, adjust_prob, perm_gro_fac, crra, v_func_adj, v_func_fxd +): + """ + Calculate "intermediate" value from next period's bank balances, the + income shocks shocks, and the risky asset share. + """ + m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) + + v_adj = v_func_adj(m_nrm) + if adjust_prob < 1.0: + v_fxd = v_func_fxd(m_nrm, share) + # Combine by adjustment probability + v_next = adjust_prob * v_adj + (1.0 - adjust_prob) * v_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + v_next = v_adj + + v_intermed = (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * v_next + return v_intermed + + +def calc_end_of_prd_v(shocks, a_nrm, share, rfree, v_func): + # Calculate future realizations of bank balances bNrm + ex_ret = shocks - rfree + r_port = rfree + share * ex_ret + b_rnm = r_port * a_nrm + + # Make an extended share_next of the same dimension as b_nrm so + # that the function can be vectorized + share_exp = np.full_like(b_rnm, share) + + return v_func(b_rnm, share_exp) + + +def calc_m_nrm_next_joint(shocks, a_nrm, share, rfree, perm_gro_fac): + """ + Calculate future realizations of market resources mNrm from the shock + distribution shocks, normalized end-of-period assets a, and risky share z. + """ + # Calculate future realizations of bank balances bNrm + ex_ret = shocks["Risky"] - rfree + r_port = rfree + share * ex_ret + b_nrm = r_port * a_nrm + return b_nrm / (shocks["PermShk"] * perm_gro_fac) + shocks["TranShk"] + + +def calc_end_of_prd_dvdx_joint( + shocks, + a_nrm, + share, + rfree, + adjust_prob, + perm_gro_fac, + crra, + vp_func_adj, + dvdm_func_fxd, + dvds_func_fxd, +): + """ + Evaluate end-of-period marginal value of assets and risky share based + on the shock distribution S, values of bend of period assets a, and + risky share z. + """ + m_nrm = calc_m_nrm_next_joint(shocks, a_nrm, share, rfree, perm_gro_fac) + ex_ret = shocks["Risky"] - rfree + r_port = rfree + share * ex_ret + dvdm_adj = vp_func_adj(m_nrm) + # No marginal value of Share if it's a free choice! + dvds_adj = np.zeros_like(m_nrm) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + share_exp = np.full_like(m_nrm, share) + dvdm_fxd = dvdm_func_fxd(m_nrm, share_exp) + dvds_fxd = dvds_func_fxd(m_nrm, share_exp) + # Combine by adjustment probability + dvdm_next = adjust_prob * dvdm_adj + (1.0 - adjust_prob) * dvdm_fxd + dvds_next = adjust_prob * dvds_adj + (1.0 - adjust_prob) * dvds_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + dvdm_next = dvdm_adj + dvds_next = dvds_adj + + perm_shk_fac = shocks["PermShk"] * perm_gro_fac + temp_fac = perm_shk_fac ** (-crra) * dvdm_next + eop_dvda = r_port * temp_fac + eop_dvds = ex_ret * a_nrm * temp_fac + perm_shk_fac ** (1 - crra) * dvds_next + + return eop_dvda, eop_dvds + + +def calc_end_of_prd_v_joint( + shocks, a_nrm, share, rfree, adjust_prob, perm_gro_fac, crra, v_func_adj, v_func_fxd +): + """ + Evaluate end-of-period value, based on the shock distribution S, values + of bank balances bNrm, and values of the risky share z. + """ + m_nrm = calc_m_nrm_next_joint(shocks, a_nrm, share, rfree, perm_gro_fac) + v_adj = v_func_adj(m_nrm) + + if adjust_prob < 1.0: + # Expand to the same dimensions as mNrm + share_exp = np.full_like(m_nrm, share) + v_fxd = v_func_fxd(m_nrm, share_exp) + # Combine by adjustment probability + v_next = adjust_prob * v_adj + (1.0 - adjust_prob) * v_fxd + else: # Don't bother evaluating if there's no chance that portfolio share is fixed + v_next = v_adj + + return (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * v_next + + def solve_one_period_ConsPortfolio( solution_next, ShockDstn, @@ -453,148 +707,68 @@ def solve_one_period_ConsPortfolio( bNrmNext, ShareNext = np.meshgrid(bNrmGrid, ShareGrid, indexing="ij") # Define functions that are used internally to evaluate future realizations - def calc_mNrm_next(S, b): - """ - Calculate future realizations of market resources mNrm from the income - shock distribution S and normalized bank balances b. - """ - return b / (S["PermShk"] * PermGroFac) + S["TranShk"] - - def calc_dvdm_next(S, b, z): - """ - Evaluate realizations of marginal value of market resources next period, - based on the income distribution S, values of bank balances bNrm, and - values of the risky share z. - """ - mNrm_next = calc_mNrm_next(S, b) - dvdmAdj_next = vPfuncAdj_next(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - dvdmFxd_next = dvdmFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - dvdm_next = AdjustPrb * dvdmAdj_next + (1.0 - AdjustPrb) * dvdmFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - dvdm_next = dvdmAdj_next - - dvdm_next = (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next - return dvdm_next - - def calc_dvds_next(S, b, z): - """ - Evaluate realizations of marginal value of risky share next period, based - on the income distribution S, values of bank balances bNrm, and values of - the risky share z. - """ - mNrm_next = calc_mNrm_next(S, b) - - # No marginal value of Share if it's a free choice! - dvdsAdj_next = np.zeros_like(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - dvdsFxd_next = dvdsFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - dvds_next = AdjustPrb * dvdsAdj_next + (1.0 - AdjustPrb) * dvdsFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - dvds_next = dvdsAdj_next - - dvds_next = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * dvds_next - return dvds_next # Calculate end-of-period marginal value of assets and shares at each point # in aNrm and ShareGrid. Does so by taking expectation of next period marginal # values across income and risky return shocks. - # Calculate intermediate marginal value of bank balances by taking expectations over income shocks - dvdb_intermed = expected(calc_dvdm_next, IncShkDstn, args=(bNrmNext, ShareNext)) + # Calculate intermediate marginal value of bank balances and risky portfolio share + # by taking expectations over income shocks + + dvdb_intermed, dvds_intermed = expected( + calc_dvdx_next, + IncShkDstn, + args=( + bNrmNext, + ShareNext, + AdjustPrb, + PermGroFac, + CRRA, + vPfuncAdj_next, + dvdmFuncFxd_next, + dvdsFuncFxd_next, + ), + ) + dvdbNvrs_intermed = uFunc.derinv(dvdb_intermed, order=(1, 0)) dvdbNvrsFunc_intermed = BilinearInterp(dvdbNvrs_intermed, bNrmGrid, ShareGrid) dvdbFunc_intermed = MargValueFuncCRRA(dvdbNvrsFunc_intermed, CRRA) - # Calculate intermediate marginal value of risky portfolio share by taking expectations over income shocks - dvds_intermed = expected(calc_dvds_next, IncShkDstn, args=(bNrmNext, ShareNext)) dvdsFunc_intermed = BilinearInterp(dvds_intermed, bNrmGrid, ShareGrid) # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij") # Define functions for calculating end-of-period marginal value - def calc_EndOfPrd_dvda(S, a, z): - """ - Compute end-of-period marginal value of assets at values a, conditional - on risky asset return S and risky share z. - """ - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree # Excess returns - Rport = Rfree + z * Rxs # Portfolio return - bNrm_next = Rport * a - - # Ensure shape concordance - z_rep = z + np.zeros_like(bNrm_next) - - # Calculate and return dvda - EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_dvda - - def calc_EndOfPrddvds(S, a, z): - """ - Compute end-of-period marginal value of risky share at values a, conditional - on risky asset return S and risky share z. - """ - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree # Excess returns - Rport = Rfree + z * Rxs # Portfolio return - bNrm_next = Rport * a - - # Make the shares match the dimension of b, so that it can be vectorized - z_rep = z + np.zeros_like(bNrm_next) - - # Calculate and return dvds - EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed( - bNrm_next, z_rep - ) + dvdsFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_dvds # Evaluate realizations of value and marginal value after asset returns are realized - # Calculate end-of-period marginal value of assets by taking expectations - EndOfPrd_dvda = DiscFacEff * expected( - calc_EndOfPrd_dvda, RiskyDstn, args=(aNrmNow, ShareNext) - ) - EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) + # Calculate end-of-period marginal value of assets and risky portfolio share + # by taking expectations - # Calculate end-of-period marginal value of risky portfolio share by taking expectations - EndOfPrd_dvds = DiscFacEff * expected( - calc_EndOfPrddvds, RiskyDstn, args=(aNrmNow, ShareNext) + EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected( + calc_end_of_prd_dvdx, + RiskyDstn, + args=(aNrmNow, ShareNext, Rfree, dvdbFunc_intermed, dvdsFunc_intermed), ) + EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) + # Make the end-of-period value function if the value function is requested if vFuncBool: - - def calc_v_intermed(S, b, z): - """ - Calculate "intermediate" value from next period's bank balances, the - income shocks S, and the risky asset share. - """ - mNrm_next = calc_mNrm_next(S, b) - - vAdj_next = vFuncAdj_next(mNrm_next) - if AdjustPrb < 1.0: - vFxd_next = vFuncFxd_next(mNrm_next, z) - # Combine by adjustment probability - v_next = AdjustPrb * vAdj_next + (1.0 - AdjustPrb) * vFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - v_next = vAdj_next - - v_intermed = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next - return v_intermed - # Calculate intermediate value by taking expectations over income shocks v_intermed = expected( - calc_v_intermed, IncShkDstn, args=(bNrmNext, ShareNext) + calc_v_intermed, + IncShkDstn, + args=( + bNrmNext, + ShareNext, + AdjustPrb, + PermGroFac, + CRRA, + vFuncAdj_next, + vFuncFxd_next, + ), ) # Construct the "intermediate value function" for this period @@ -602,22 +776,11 @@ def calc_v_intermed(S, b, z): vNvrsFunc_intermed = BilinearInterp(vNvrs_intermed, bNrmGrid, ShareGrid) vFunc_intermed = ValueFuncCRRA(vNvrsFunc_intermed, CRRA) - def calc_EndOfPrd_v(S, a, z): - # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree - Rport = Rfree + z * Rxs - bNrm_next = Rport * a - - # Make an extended share_next of the same dimension as b_nrm so - # that the function can be vectorized - z_rep = z + np.zeros_like(bNrm_next) - - EndOfPrd_v = vFunc_intermed(bNrm_next, z_rep) - return EndOfPrd_v - # Calculate end-of-period value by taking expectations EndOfPrd_v = DiscFacEff * expected( - calc_EndOfPrd_v, RiskyDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_v, + RiskyDstn, + args=(aNrmNow, ShareNext, Rfree, vFunc_intermed), ) EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v) @@ -634,76 +797,24 @@ def calc_EndOfPrd_v(S, a, z): aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij") # Define functions that are used internally to evaluate future realizations - def calc_mNrm_next(S, a, z): - """ - Calculate future realizations of market resources mNrm from the shock - distribution S, normalized end-of-period assets a, and risky share z. - """ - # Calculate future realizations of bank balances bNrm - Rxs = S["Risky"] - Rfree - Rport = Rfree + z * Rxs - bNrm_next = Rport * a - mNrm_next = bNrm_next / (S["PermShk"] * PermGroFac) + S["TranShk"] - return mNrm_next - - def calc_EndOfPrd_dvdx(S, a, z): - """ - Evaluate end-of-period marginal value of assets and risky share based - on the shock distribution S, values of bend of period assets a, and - risky share z. - """ - mNrm_next = calc_mNrm_next(S, a, z) - Rxs = S["Risky"] - Rfree - Rport = Rfree + z * Rxs - dvdmAdj_next = vPfuncAdj_next(mNrm_next) - # No marginal value of Share if it's a free choice! - dvdsAdj_next = np.zeros_like(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - dvdmFxd_next = dvdmFuncFxd_next(mNrm_next, Share_next_expanded) - dvdsFxd_next = dvdsFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - dvdm_next = AdjustPrb * dvdmAdj_next + (1.0 - AdjustPrb) * dvdmFxd_next - dvds_next = AdjustPrb * dvdsAdj_next + (1.0 - AdjustPrb) * dvdsFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - dvdm_next = dvdmAdj_next - dvds_next = dvdsAdj_next - - EndOfPrd_dvda = Rport * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next - EndOfPrd_dvds = ( - Rxs * a * (S["PermShk"] * PermGroFac) ** (-CRRA) * dvdm_next - + (S["PermShk"] * PermGroFac) ** (1 - CRRA) * dvds_next - ) - - return EndOfPrd_dvda, EndOfPrd_dvds - - def calc_EndOfPrd_v(S, a, z): - """ - Evaluate end-of-period value, based on the shock distribution S, values - of bank balances bNrm, and values of the risky share z. - """ - mNrm_next = calc_mNrm_next(S, a, z) - vAdj_next = vFuncAdj_next(mNrm_next) - - if AdjustPrb < 1.0: - # Expand to the same dimensions as mNrm - Share_next_expanded = z + np.zeros_like(mNrm_next) - vFxd_next = vFuncFxd_next(mNrm_next, Share_next_expanded) - # Combine by adjustment probability - v_next = AdjustPrb * vAdj_next + (1.0 - AdjustPrb) * vFxd_next - else: # Don't bother evaluating if there's no chance that portfolio share is fixed - v_next = vAdj_next - - EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next - return EndOfPrd_v # Evaluate realizations of value and marginal value after asset returns are realized # Calculate end-of-period marginal value of assets and risky share by taking expectations EndOfPrd_dvda, EndOfPrd_dvds = DiscFacEff * expected( - calc_EndOfPrd_dvdx, ShockDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_dvdx_joint, + ShockDstn, + args=( + aNrmNow, + ShareNext, + Rfree, + AdjustPrb, + PermGroFac, + CRRA, + vPfuncAdj_next, + dvdmFuncFxd_next, + dvdsFuncFxd_next, + ), ) EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) @@ -711,7 +822,18 @@ def calc_EndOfPrd_v(S, a, z): if vFuncBool: # Calculate end-of-period value, its derivative, and their pseudo-inverse EndOfPrd_v = DiscFacEff * expected( - calc_EndOfPrd_v, ShockDstn, args=(aNrmNow, ShareNext) + calc_end_of_prd_v_joint, + ShockDstn, + args=( + aNrmNow, + ShareNext, + Rfree, + AdjustPrb, + PermGroFac, + CRRA, + vFuncAdj_next, + vFuncFxd_next, + ), ) EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v) @@ -841,17 +963,6 @@ def calc_EndOfPrd_v(S, a, z): ShareAdj_now = np.insert(ShareAdj_now, 0, Share_lower_bound) ShareFuncAdj_now = LinearInterp(mNrmAdj_now, ShareAdj_now, ShareLimit, 0.0) - # This is a point at which (a,c,share) have consistent length. Take the - # snapshot for storing the grid and values in the solution. - save_points = { - "a": deepcopy(aNrmGrid), - "eop_dvda_adj": uFunc.der(cNrmAdj_now), - "share_adj": deepcopy(ShareAdj_now), - "share_grid": deepcopy(ShareGrid), - "eop_dvda_fxd": uFunc.der(EndOfPrd_dvda), - "eop_dvds_fxd": EndOfPrd_dvds, - } - # Add the value function if requested if vFuncBool: # Create the value functions for this period, defined over market resources @@ -909,13 +1020,6 @@ def calc_EndOfPrd_v(S, a, z): dvdsFuncFxd=dvdsFuncFxd_now, vFuncFxd=vFuncFxd_now, AdjPrb=AdjustPrb, - # WHAT IS THIS STUFF FOR?? - aGrid=save_points["a"], - Share_adj=save_points["share_adj"], - EndOfPrddvda_adj=save_points["eop_dvda_adj"], - ShareGrid=save_points["share_grid"], - EndOfPrddvda_fxd=save_points["eop_dvda_fxd"], - EndOfPrddvds_fxd=save_points["eop_dvds_fxd"], ) return solution_now From 9764b2ae95618085013d7c373dffc7026c28d9b0 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 27 Mar 2024 19:12:10 -0400 Subject: [PATCH 08/14] fix import legacy solver --- HARK/ConsumptionSaving/ConsWealthPortfolioModel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py b/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py index 27a8c7488..90022901f 100644 --- a/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py @@ -5,10 +5,10 @@ from scipy.optimize import fixed_point, minimize_scalar, root from HARK.ConsumptionSaving.ConsPortfolioModel import ( - ConsPortfolioSolver, PortfolioConsumerType, init_portfolio, ) +from HARK.ConsumptionSaving.LegacyOOsolvers import ConsPortfolioSolver from HARK.core import make_one_period_oo_solver from HARK.distribution import DiscreteDistribution, calc_expectation from HARK.interpolation import ( From 758b3047ab515aadcfbb63cdc3f053f5f510a73a Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Thu, 28 Mar 2024 12:22:09 -0400 Subject: [PATCH 09/14] quick fix for wealthPortfolioModel --- HARK/ConsumptionSaving/ConsWealthPortfolioModel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py b/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py index 90022901f..bb4bfbb16 100644 --- a/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py @@ -97,8 +97,8 @@ def post_solve(self): solution.cFuncAdj = solution.cFunc solution.cFuncFxd = lambda m, s: solution.cFunc(m) share = solution.shareFunc - solution.ShareFuncAdj = lambda m: np.clip(share(m), 0.0, 1.0) - solution.ShareFuncFxd = lambda m, s: np.clip(share(m), 0.0, 1.0) + solution.ShareFuncAdj = share + solution.ShareFuncFxd = lambda m, s: share(m) @dataclass From 43612aa9ba4fafa649ed3c5535dc028b1f2d4765 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Fri, 5 Apr 2024 15:21:25 -0400 Subject: [PATCH 10/14] move functions out of portfolio solver --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 38 +++++++++++++------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index f63b9d4f6..bcd8f35f0 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -309,6 +309,25 @@ def get_controls(self): self.controls["Share"] = ShareNow +def calc_radj(shock, share_limit, rfree, crra): + rport = share_limit * shock + (1.0 - share_limit) * rfree + return rport ** (1.0 - crra) + + +def calc_h_nrm(shocks, perm_gro_fac, share_limit, rfree, crra, h_nrm_next): + perm_shk_fac = perm_gro_fac * shocks["PermShk"] + rport = share_limit * shocks["Risky"] + (1.0 - share_limit) * rfree + hNrm = (perm_shk_fac / rport**crra) * (shocks["TranShk"] + h_nrm_next) + return hNrm + + +def calc_h_nrm_joint(shocks, perm_gro_fac, share_limit, rfree, h_nrm_next): + perm_shk_fac = perm_gro_fac * shocks["PermShk"] + rport = share_limit * shocks["Risky"] + (1.0 - share_limit) * rfree + hNrm = (perm_shk_fac / rport) * (shocks["TranShk"] + h_nrm_next) + return hNrm + + def calc_m_nrm_next(shocks, b_nrm, perm_gro_fac): """ Calculate future realizations of market resources mNrm from the income @@ -688,26 +707,19 @@ def solve_one_period_ConsPortfolio( # Perform an alternate calculation of the absolute patience factor when # returns are risky. This uses the Merton-Samuelson limiting risky share, # which is what's relevant as mNrm goes to infinity. - def calc_Radj(R): - Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree - return Rport ** (1.0 - CRRA) - R_adj = expected(calc_Radj, RiskyDstn)[0] + R_adj = expected(calc_radj, RiskyDstn, args=(ShareLimit, Rfree, CRRA))[0] PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA) MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) # Also perform an alternate calculation for human wealth under risky returns - def calc_hNrm(S): - Risky = S["Risky"] - PermShk = S["PermShk"] - TranShk = S["TranShk"] - G = PermGroFac * PermShk - Rport = ShareLimit * Risky + (1.0 - ShareLimit) * Rfree - hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm) - return hNrm # This correctly accounts for risky returns and risk aversion - hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj + hNrmNow = expected( + calc_h_nrm_joint, + ShockDstn, + args=(PermGroFac, ShareLimit, Rfree, solution_next.hNrm), + ) # This basic equation works if there's no correlation among shocks # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) From 755295cb4158f7c53aac77343d609e8a74f7c78f Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Mon, 15 Apr 2024 20:14:54 -0400 Subject: [PATCH 11/14] update MPCmaxUnc --- HARK/ConsumptionSaving/ConsIndShockModel.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index d75549e36..8cac3037d 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -401,10 +401,8 @@ def calc_m_nrm_min(boro_const_art, boro_const_nat): def calc_mpc_max( mpc_max_next, worst_inc_prob, crra, pat_fac, boro_const_nat, boro_const_art ): - m_nrm_min = calc_m_nrm_min(boro_const_art, boro_const_nat) temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac - mpc_max = 1.0 / (1.0 + temp_fac / mpc_max_next) - return 1.0 if boro_const_nat < m_nrm_min else mpc_max + return 1.0 / (1.0 + temp_fac / mpc_max_next) def calc_m_nrm_next(shock, a, rfree, perm_gro_fac): @@ -510,9 +508,10 @@ def solve_one_period_ConsIndShock( MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds - MPCmaxNow = calc_mpc_max( + MPCmaxUnc = calc_mpc_max( solution_next.MPCmax, WorstIncPrb, CRRA, PatFac, BoroCnstNat, BoroCnstArt ) + MPCmaxNow = 1.0 if BoroCnstNat < mNrmMinNow else MPCmaxUnc cFuncLimitIntercept = MPCminNow * hNrmNow cFuncLimitSlope = MPCminNow @@ -553,7 +552,7 @@ def solve_one_period_ConsIndShock( ) dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2) MPC = dcda / (dcda + 1.0) - MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow) + MPC_for_interpolation = np.insert(MPC, 0, MPCmaxUnc) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( @@ -758,9 +757,10 @@ def solve_one_period_ConsKinkedR( MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds - MPCmaxNow = calc_mpc_max( + MPCmaxUnc = calc_mpc_max( solution_next.MPCmax, WorstIncPrb, CRRA, PatFacBoro, BoroCnstNat, BoroCnstArt ) + MPCmaxNow = 1.0 if BoroCnstNat < mNrmMinNow else MPCmaxUnc cFuncLimitIntercept = MPCminNow * hNrmNow cFuncLimitSlope = MPCminNow @@ -809,7 +809,7 @@ def solve_one_period_ConsKinkedR( ) dcda = EndOfPrdvPP / uFunc.der(np.array(cNrmNow), order=2) MPC = dcda / (dcda + 1.0) - MPC_for_interpolation = np.insert(MPC, 0, MPCmaxNow) + MPC_for_interpolation = np.insert(MPC, 0, MPCmaxUnc) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( From 5e7b6d7be8b8ba3715385e616c154057cb791e1f Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Mon, 15 Apr 2024 20:18:33 -0400 Subject: [PATCH 12/14] remove joint h_nrm and calculate using Matt's method --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index 84940bbb5..eb9259fba 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -316,13 +316,6 @@ def calc_h_nrm(shocks, perm_gro_fac, share_limit, rfree, crra, h_nrm_next): return hNrm -def calc_h_nrm_joint(shocks, perm_gro_fac, share_limit, rfree, h_nrm_next): - perm_shk_fac = perm_gro_fac * shocks["PermShk"] - rport = share_limit * shocks["Risky"] + (1.0 - share_limit) * rfree - hNrm = (perm_shk_fac / rport) * (shocks["TranShk"] + h_nrm_next) - return hNrm - - def calc_m_nrm_next(shocks, b_nrm, perm_gro_fac): """ Calculate future realizations of market resources mNrm from the income @@ -710,10 +703,13 @@ def solve_one_period_ConsPortfolio( # Also perform an alternate calculation for human wealth under risky returns # This correctly accounts for risky returns and risk aversion - hNrmNow = expected( - calc_h_nrm_joint, - ShockDstn, - args=(PermGroFac, ShareLimit, Rfree, solution_next.hNrm), + hNrmNow = ( + expected( + calc_h_nrm, + ShockDstn, + args=(PermGroFac, ShareLimit, Rfree, CRRA, solution_next.hNrm), + ) + / R_adj ) # This basic equation works if there's no correlation among shocks From efa120ac7854d7bd72a4c6bb039600961dbd0377 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Wed, 1 May 2024 10:27:22 -0400 Subject: [PATCH 13/14] add docstrings to ConsIndShockModel --- HARK/ConsumptionSaving/ConsIndShockModel.py | 104 +++++++++++++++++++- 1 file changed, 99 insertions(+), 5 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 8cac3037d..39005d593 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -203,14 +203,37 @@ def append_solution(self, new_solution): def calc_human_wealth(h_nrm_next, perm_gro_fac, rfree, ex_inc_next): + """Calculate human wealth this period given human wealth next period. + + Args: + h_nrm_next (float): Normalized human wealth next period. + perm_gro_fac (float): Permanent income growth factor. + rfree (float): Risk free interest factor. + ex_inc_next (float): Expected income next period. + """ return (perm_gro_fac / rfree) * (h_nrm_next + ex_inc_next) -def calc_pat_fac(rfree, disc_fac_eff, crra): +def calc_patience_factor(rfree, disc_fac_eff, crra): + """Calculate the patience factor for the agent. + + Args: + rfree (float): Risk free interest factor. + disc_fac_eff (float): Effective discount factor. + crra (float): Coefficient of relative risk aversion. + + """ return ((rfree * disc_fac_eff) ** (1.0 / crra)) / rfree def calc_mpc_min(mpc_min_next, pat_fac): + """Calculate the lower bound of the marginal propensity to consume. + + Args: + mpc_min_next (float): Lower bound of the marginal propensity to + consume next period. + pat_fac (float): Patience factor. + """ return 1.0 / (1.0 + pat_fac / mpc_min_next) @@ -269,7 +292,7 @@ def solve_one_period_ConsPF( hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rfree, 1.0) # Calculate the lower bound of the marginal propensity to consume - PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) + PatFac = calc_patience_factor(Rfree, DiscFacEff, CRRA) MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac) # Extract the discrete kink points in next period's consumption function; @@ -377,6 +400,11 @@ def solve_one_period_ConsPF( def calc_worst_inc_prob(inc_shk_dstn): + """Calculate the probability of the worst income shock. + + Args: + inc_shk_dstn (DiscreteDistribution): Distribution of shocks to income. + """ probs = inc_shk_dstn.pmv perm, tran = inc_shk_dstn.atoms income = perm * tran @@ -385,12 +413,26 @@ def calc_worst_inc_prob(inc_shk_dstn): def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): + """Calculate the natural borrowing constraint. + + Args: + m_nrm_min_next (float): Minimum normalized market resources next period. + inc_shk_dstn (DiscreteDstn): Distribution of shocks to income. + rfree (float): Risk free interest factor. + perm_gro_fac (float): Permanent income growth factor. + """ perm, tran = inc_shk_dstn.atoms temp_fac = (perm_gro_fac * np.min(perm)) / rfree return (m_nrm_min_next - np.min(tran)) * temp_fac def calc_m_nrm_min(boro_const_art, boro_const_nat): + """Calculate the minimum normalized market resources this period. + + Args: + boro_const_art (float): Artificial borrowing constraint. + boro_const_nat (float): Natural borrowing constraint. + """ return ( boro_const_nat if boro_const_art is None @@ -401,27 +443,79 @@ def calc_m_nrm_min(boro_const_art, boro_const_nat): def calc_mpc_max( mpc_max_next, worst_inc_prob, crra, pat_fac, boro_const_nat, boro_const_art ): + """Calculate the upper bound of the marginal propensity to consume. + + Args: + mpc_max_next (float): Upper bound of the marginal propensity to + consume next period. + worst_inc_prob (float): Probability of the worst income shock. + crra (float): Coefficient of relative risk aversion. + pat_fac (float): Patience factor. + boro_const_nat (float): Natural borrowing constraint. + boro_const_art (float): Artificial borrowing constraint. + """ temp_fac = (worst_inc_prob ** (1.0 / crra)) * pat_fac return 1.0 / (1.0 + temp_fac / mpc_max_next) def calc_m_nrm_next(shock, a, rfree, perm_gro_fac): + """Calculate normalized market resources next period. + + Args: + shock (float): Realization of shocks to income. + a (np.ndarray): Exogenous grid of end-of-period assets. + rfree (float): Risk free interest factor. + perm_gro_fac (float): Permanent income growth factor. + """ return rfree / (perm_gro_fac * shock["PermShk"]) * a + shock["TranShk"] def calc_v_next(shock, a, rfree, crra, perm_gro_fac, vfunc_next): + """Calculate continuation value function with respect to + end-of-period assets. + + Args: + shock (float): Realization of shocks to income. + a (np.ndarray): Exogenous grid of end-of-period assets. + rfree (float): Risk free interest factor. + crra (float): Coefficient of relative risk aversion. + perm_gro_fac (float): Permanent income growth factor. + vfunc_next (Callable): Value function next period. + """ return ( shock["PermShk"] ** (1.0 - crra) * perm_gro_fac ** (1.0 - crra) ) * vfunc_next(calc_m_nrm_next(shock, a, rfree, perm_gro_fac)) def calc_vp_next(shock, a, rfree, crra, perm_gro_fac, vp_func_next): + """Calculate the continuation marginal value function with respect to + end-of-period assets. + + Args: + shock (float): Realization of shocks to income. + a (np.ndarray): Exogenous grid of end-of-period assets. + rfree (float): Risk free interest factor. + crra (float): Coefficient of relative risk aversion. + perm_gro_fac (float): Permanent income growth factor. + vp_func_next (Callable): Marginal value function next period. + """ return shock["PermShk"] ** (-crra) * vp_func_next( calc_m_nrm_next(shock, a, rfree, perm_gro_fac), ) def calc_vpp_next(shock, a, rfree, crra, perm_gro_fac, vppfunc_next): + """Calculate the continuation marginal marginal value function + with respect to end-of-period assets. + + Args: + shock (float): Realization of shocks to income. + a (np.ndarray): Exogenous grid of end-of-period assets. + rfree (float): Risk free interest factor. + crra (float): Coefficient of relative risk aversion. + perm_gro_fac (float): Permanent income growth factor. + vppfunc_next (Callable): Marginal marginal value function next period. + """ return shock["PermShk"] ** (-crra - 1.0) * vppfunc_next( calc_m_nrm_next(shock, a, rfree, perm_gro_fac), ) @@ -504,7 +598,7 @@ def solve_one_period_ConsIndShock( mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat) # Update the bounding MPCs and PDV of human wealth: - PatFac = calc_pat_fac(Rfree, DiscFacEff, CRRA) + PatFac = calc_patience_factor(Rfree, DiscFacEff, CRRA) MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFac) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds @@ -752,8 +846,8 @@ def solve_one_period_ConsKinkedR( mNrmMinNow = calc_m_nrm_min(BoroCnstArt, BoroCnstNat) # Update the bounding MPCs and PDV of human wealth: - PatFacSave = calc_pat_fac(Rsave, DiscFacEff, CRRA) - PatFacBoro = calc_pat_fac(Rboro, DiscFacEff, CRRA) + PatFacSave = calc_patience_factor(Rsave, DiscFacEff, CRRA) + PatFacBoro = calc_patience_factor(Rboro, DiscFacEff, CRRA) MPCminNow = calc_mpc_min(solution_next.MPCmin, PatFacSave) # Set the upper limit of the MPC (at mNrmMinNow) based on whether the natural # or artificial borrowing constraint actually binds From fbbaeeff5128b9bcfff4d4646519427e4b229834 Mon Sep 17 00:00:00 2001 From: alanlujan91 Date: Thu, 2 May 2024 15:37:20 -0400 Subject: [PATCH 14/14] update ConsPortfolioModel docs --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 32 ++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index eb9259fba..26f00d4b9 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -305,11 +305,29 @@ def get_controls(self): def calc_radj(shock, share_limit, rfree, crra): + """Expected rate of return adjusted by CRRA + + Args: + shock (DiscreteDistribution): Distribution of risky asset returns + share_limit (float): limiting lower bound of risky portfolio share + rfree (float): Risk free interest rate + crra (float): Coefficient of relative risk aversion + """ rport = share_limit * shock + (1.0 - share_limit) * rfree return rport ** (1.0 - crra) -def calc_h_nrm(shocks, perm_gro_fac, share_limit, rfree, crra, h_nrm_next): +def calc_human_wealth(shocks, perm_gro_fac, share_limit, rfree, crra, h_nrm_next): + """Calculate human wealth this period given human wealth next period. + + Args: + shocks (DiscreteDistribution): Joint distribution of shocks to income and returns. + perm_gro_fac (float): Permanent income growth factor + share_limit (float): limiting lower bound of risky portfolio share + rfree (float): Risk free interest rate + crra (float): Coefficient of relative risk aversion + h_nrm_next (float): Human wealth next period + """ perm_shk_fac = perm_gro_fac * shocks["PermShk"] rport = share_limit * shocks["Risky"] + (1.0 - share_limit) * rfree hNrm = (perm_shk_fac / rport**crra) * (shocks["TranShk"] + h_nrm_next) @@ -385,6 +403,11 @@ def calc_dvdx_next( dvdm_func_fxd, dvds_func_fxd, ): + """ + Evaluate realizations of marginal values next period, based + on the income distribution "shocks", values of bank balances bNrm, and values of + the risky share z. + """ m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac) dvdm_adj = vp_func_adj(m_nrm) # No marginal value of shockshare if it's a free choice! @@ -445,6 +468,10 @@ def calc_end_of_prd_dvds(shocks, a_nrm, share, rfree, dvdb_func, dvds_func): def calc_end_of_prd_dvdx(shocks, a_nrm, share, rfree, dvdb_func, dvds_func): + """ + Compute end-of-period marginal values at values a, conditional + on risky asset return shocks and risky share z. + """ # Calculate future realizations of bank balances bNrm ex_ret = shocks - rfree # Excess returns r_port = rfree + share * ex_ret # Portfolio return @@ -480,6 +507,7 @@ def calc_v_intermed( def calc_end_of_prd_v(shocks, a_nrm, share, rfree, v_func): + """Compute end-of-period values.""" # Calculate future realizations of bank balances bNrm ex_ret = shocks - rfree r_port = rfree + share * ex_ret @@ -705,7 +733,7 @@ def solve_one_period_ConsPortfolio( # This correctly accounts for risky returns and risk aversion hNrmNow = ( expected( - calc_h_nrm, + calc_human_wealth, ShockDstn, args=(PermGroFac, ShareLimit, Rfree, CRRA, solution_next.hNrm), )