From 821f0fd721a0090ecd8a52db9ac56e760d829e7f Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Thu, 28 Mar 2024 15:40:44 -0400 Subject: [PATCH 1/8] Found correct human wealth formula with risky returns Through a combination of math and experimentation, I've found the correct computation of "human wealth" in the presence of risky returns. I've implemented it for RiskyAssetModel and FixedShareModel. On the next commit I'll get it working for the PortfolioChoice models, and then turn to ConsMarkovModel. --- HARK/ConsumptionSaving/ConsRiskyAssetModel.py | 38 +++++++++++-------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py index 965569533..9356dd8ae 100644 --- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py @@ -528,24 +528,26 @@ def solve_one_period_ConsIndShockRiskyAsset( def calc_Radj(R): return R ** (1.0 - CRRA) - PatFac = (DiscFacEff * expected(calc_Radj, RiskyDstn)) ** (1.0 / CRRA) + Radj = expected(calc_Radj, RiskyDstn) + PatFac = (DiscFacEff * Radj) ** (1.0 / CRRA) MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) + MPCminNow = MPCminNow[0] # Returns as one element array, extract # Also perform an alternate calculation for human wealth under risky returns def calc_hNrm(S): Risky = S["Risky"] PermShk = S["PermShk"] TranShk = S["TranShk"] - hNrm = (PermGroFac / Risky) * (PermShk * TranShk + solution_next.hNrm) + hNrm = (PermGroFac / Risky**CRRA) * (PermShk * TranShk + solution_next.hNrm) return hNrm - hNrmNow = expected(calc_hNrm, ShockDstn) + # This correctly incorporates risk aversion and risky returns + hNrmNow = expected(calc_hNrm, ShockDstn) / Radj + hNrmNow = hNrmNow[0] - # The above attempts to pin down the limiting consumption function for this - # model, however it is not clear why it creates bugs, so for now we allow - # for a linear extrapolation beyond the last asset point - cFuncLimitIntercept = None - cFuncLimitSlope = None + # Use adjusted MPCmin and hNrm to specify limiting linear behavior of cFunc + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow # Returns as one element array, extract # Calculate the minimum allowable value of market resources in this period BoroCnstNat_cand = ( @@ -812,6 +814,7 @@ def calc_vPPnext(S, a): # income distribution is independent from the return distribution if CubicBool: # Construct the unconstrained consumption function as a cubic interpolation + cFuncNowUnc = CubicInterp( m_for_interpolation, c_for_interpolation, @@ -1410,27 +1413,32 @@ def solve_one_period_FixedShareRiskyAsset( # Perform an alternate calculation of the absolute patience factor when returns are risky def calc_Radj(R): - R_temp = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree - return R_temp ** (1.0 - CRRA) + Rport = RiskyShareFixed * R + (1.0 - RiskyShareFixed) * Rfree + return Rport ** (1.0 - CRRA) - PatFac = (DiscFacEff * expected(calc_Radj, RiskyDstn)) ** (1.0 / CRRA) + R_adj = expected(calc_Radj, RiskyDstn) + PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA) MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) + MPCminNow = MPCminNow[0] # Also perform an alternate calculation for human wealth under risky returns def calc_hNrm(S): Risky = S["Risky"] PermShk = S["PermShk"] TranShk = S["TranShk"] - hNrm = (PermGroFac / Risky) * (PermShk * TranShk + solution_next.hNrm) + Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree + hNrm = (PermGroFac / Rport**CRRA) * (PermShk * TranShk + solution_next.hNrm) return hNrm - hNrmNow = expected(calc_hNrm, ShockDstn) + # This correctly accounts for risky returns and risk aversion + hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj + hNrmNow = hNrmNow[0] # The above attempts to pin down the limiting consumption function for this # model, however it is not clear why it creates bugs, so for now we allow # for a linear extrapolation beyond the last asset point - cFuncLimitIntercept = None - cFuncLimitSlope = None + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow # Calculate the minimum allowable value of market resources in this period BoroCnstNat_cand = ( From f36700e97204a96cd48cb2aa11460ce4786efa29 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Fri, 29 Mar 2024 16:24:41 -0400 Subject: [PATCH 2/8] Added cFunc extrapolation to basic portfolio model Seems to mostly work, *except* that MPCmin is now off by 0.00001 or so. This only becomes evident if you have an enormous top gridpoint, like m=10000. I have a guess why. --- HARK/ConsumptionSaving/ConsRiskyAssetModel.py | 42 +++++++++++++++++-- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py index 9356dd8ae..d106550ad 100644 --- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py @@ -538,7 +538,8 @@ def calc_hNrm(S): Risky = S["Risky"] PermShk = S["PermShk"] TranShk = S["TranShk"] - hNrm = (PermGroFac / Risky**CRRA) * (PermShk * TranShk + solution_next.hNrm) + G = PermGroFac * PermShk + hNrm = (G / Risky**CRRA) * (TranShk + solution_next.hNrm) return hNrm # This correctly incorporates risk aversion and risky returns @@ -985,6 +986,38 @@ def solve_one_period_ConsPortChoice( RiskyMax = np.max(Risky_next) RiskyMin = np.min(Risky_next) + # 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) + PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA) + MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) + MPCminNow = MPCminNow[0] + + # 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 = hNrmNow[0] + + # The above attempts to pin down the limiting consumption function for this + # model, however it is not clear why it creates bugs, so for now we allow + # for a linear extrapolation beyond the last asset point + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow + # bNrm represents R*a, balances after asset return shocks but before income. # This just uses the highest risky return as a rough shifter for the aXtraGrid. if BoroCnstNat_iszero: @@ -1264,7 +1297,7 @@ def calc_EndOfPrd_v(S, a, z): # then construct the consumption function when the agent can adjust his share mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0) cNrm_now = np.insert(cNrm_now, 0, 0.0) - cFunc_now = LinearInterp(mNrm_now, cNrm_now) + cFunc_now = LinearInterp(mNrm_now, cNrm_now, cFuncLimitIntercept, cFuncLimitSlope) # Construct the marginal value (of mNrm) function vPfunc_now = MargValueFuncCRRA(cFunc_now, CRRA) @@ -1307,6 +1340,8 @@ def calc_EndOfPrd_v(S, a, z): cFunc=cFunc_now, vPfunc=vPfunc_now, vFunc=vFunc_now, + hNrm=hNrmNow, + MPCmin=MPCminNow, ) solution_now.ShareFunc = ShareFunc_now return solution_now @@ -1426,8 +1461,9 @@ def calc_hNrm(S): Risky = S["Risky"] PermShk = S["PermShk"] TranShk = S["TranShk"] + G = PermGroFac * PermShk Rport = RiskyShareFixed * Risky + (1.0 - RiskyShareFixed) * Rfree - hNrm = (PermGroFac / Rport**CRRA) * (PermShk * TranShk + solution_next.hNrm) + hNrm = (G / Rport**CRRA) * (TranShk + solution_next.hNrm) return hNrm # This correctly accounts for risky returns and risk aversion From fbb84d846f484280941614df917bc66e494d8ab9 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 2 Apr 2024 16:37:00 -0400 Subject: [PATCH 3/8] Add "refinement loop" to basic portfolio solver The proper human wealth and lower bound of MPC calculations are now in the "advanced" portfolio solver in ConsPortfolioModel, but not actually used in the consumption function. Calculation of the limiting share has been tweaked to be more accurate. The "basic" solver in ConsRiskyAssetModel now goes through a loop to refine its search for optimal portfolio share at each aNrm gridpoint by "zooming in". We will probably want to add a parameter to specify how many times to do this loop, but it's currently set to 3. --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 37 ++++++++ HARK/ConsumptionSaving/ConsRiskyAssetModel.py | 86 +++++++++++++------ 2 files changed, 98 insertions(+), 25 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index d4ecb1b0b..e2006f5ec 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -227,6 +227,8 @@ def update_solution_terminal(self): dvdmFuncFxd=dvdmFuncFxd_terminal, dvdsFuncFxd=dvdsFuncFxd_terminal, ) + self.solution_terminal.hNrm = 0.0 + self.solution_terminal.MPCmin = 1.0 def initialize_sim(self): """ @@ -428,6 +430,39 @@ def solve_one_period_ConsPortfolio( Risky_next = RiskyDstn.atoms RiskyMax = np.max(Risky_next) RiskyMin = np.min(Risky_next) + + # 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] + 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 + + # This basic equation works if there's no correlation among shocks + # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) + + # The above attempts to pin down the limiting consumption function for this + # model, however it is not clear why it creates bugs, so for now we allow + # for a linear extrapolation beyond the last asset point + cFuncLimitIntercept = MPCminNow * hNrmNow + cFuncLimitSlope = MPCminNow # bNrm represents R*a, balances after asset return shocks but before income. # This just uses the highest risky return as a rough shifter for the aXtraGrid. @@ -917,6 +952,8 @@ def calc_EndOfPrd_v(S, a, z): EndOfPrddvda_fxd=save_points["eop_dvda_fxd"], EndOfPrddvds_fxd=save_points["eop_dvds_fxd"], ) + solution_now.hNrm = hNrmNow + solution_now.MPCmin = MPCminNow return solution_now diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py index d106550ad..8e9bde4e4 100644 --- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py @@ -269,7 +269,7 @@ def temp_f(s): RiskyDstn.pmv, ) - SharePF = minimize_scalar(temp_f, bounds=(0.0, 1.0), method="bounded").x + SharePF = minimize_scalar(temp_f, bracket=(0.0, 1.0), method="golden", tol=1e-10).x[0] self.ShareLimit = SharePF self.add_to_time_inv("ShareLimit") @@ -993,10 +993,9 @@ def calc_Radj(R): Rport = ShareLimit * R + (1.0 - ShareLimit) * Rfree return Rport ** (1.0 - CRRA) - R_adj = expected(calc_Radj, RiskyDstn) + R_adj = expected(calc_Radj, RiskyDstn)[0] PatFac = (DiscFacEff * R_adj) ** (1.0 / CRRA) MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin) - MPCminNow = MPCminNow[0] # Also perform an alternate calculation for human wealth under risky returns def calc_hNrm(S): @@ -1010,11 +1009,11 @@ def calc_hNrm(S): # This correctly accounts for risky returns and risk aversion hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj - hNrmNow = hNrmNow[0] + + # This basic equation works if there's no correlation among shocks + # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) - # The above attempts to pin down the limiting consumption function for this - # model, however it is not clear why it creates bugs, so for now we allow - # for a linear extrapolation beyond the last asset point + # Define the terms for the limiting linear consumption function as m gets very big cFuncLimitIntercept = MPCminNow * hNrmNow cFuncLimitSlope = MPCminNow @@ -1055,8 +1054,8 @@ def calc_dvdm_next(S, b): based on the income distribution S and values of bank balances bNrm """ mNrm_next = calc_mNrm_next(S, b) - Gamma = S["PermShk"] * PermGroFac - dvdm_next = Gamma ** (-CRRA) * vPfunc_next(mNrm_next) + G = S["PermShk"] * PermGroFac + dvdm_next = G ** (-CRRA) * vPfunc_next(mNrm_next) return dvdm_next # Calculate end-of-period marginal value of assets and shares at each point @@ -1079,13 +1078,13 @@ def calc_dvdm_next(S, b): aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij") # Define functions for calculating end-of-period marginal value - def calc_EndOfPrd_dvda(S, a, z): + def calc_EndOfPrd_dvda(R, a, z): """ Compute end-of-period marginal value of assets at values a, conditional - on risky asset return S and risky share z. + on risky asset return R and risky share z. """ # Calculate future realizations of bank balances bNrm - Rxs = S - Rfree # Excess returns + Rxs = R - Rfree # Excess returns Rport = Rfree + z * Rxs # Portfolio return bNrm_next = Rport * a @@ -1093,23 +1092,22 @@ def calc_EndOfPrd_dvda(S, a, z): EndOfPrd_dvda = Rport * dvdbFunc_intermed(bNrm_next) return EndOfPrd_dvda - def calc_EndOfPrddvds(S, a, z): + def calc_EndOfPrd_dvds(R, 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 + Rxs = R - Rfree # Excess returns Rport = Rfree + z * Rxs # Portfolio return bNrm_next = Rport * a # Calculate and return dvds (second term is all zeros) - EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed( - bNrm_next - ) + EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(bNrm_next) return EndOfPrd_dvds # Evaluate realizations of value and marginal value after asset returns are realized + TempDstn = RiskyDstn # Calculate end-of-period marginal value of assets by taking expectations EndOfPrd_dvda = DiscFacEff * expected( @@ -1119,7 +1117,7 @@ def calc_EndOfPrddvds(S, a, z): # Calculate end-of-period marginal value of risky portfolio share by taking expectations EndOfPrd_dvds = DiscFacEff * expected( - calc_EndOfPrddvds, RiskyDstn, args=(aNrmNow, ShareNext) + calc_EndOfPrd_dvds, RiskyDstn, args=(aNrmNow, ShareNext) ) # Make the end-of-period value function if the value function is requested @@ -1186,7 +1184,7 @@ def calc_mNrm_next(S, a, z): 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 + on the shock distribution S, normalized end-of-period assets a, and risky share z. """ mNrm_next = calc_mNrm_next(S, a, z) @@ -1213,6 +1211,10 @@ def calc_EndOfPrd_v(S, a, z): v_next = vFunc_next(mNrm_next) EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next return EndOfPrd_v + + calc_EndOfPrd_dvda = lambda S, a, z : calc_EndOfPrd_dvdx(S, a, z)[0] + calc_EndOfPrd_dvds = lambda S, a, z : calc_EndOfPrd_dvdx(S, a, z)[1] + TempDstn = ShockDstn # Evaluate realizations of value and marginal value after asset returns are realized @@ -1249,18 +1251,52 @@ def calc_EndOfPrd_v(S, a, z): # Now find the optimal (continuous) risky share on [0,1] by solving the first # order condition EndOfPrd_dvds == 0. FOC_s = EndOfPrd_dvds # Relabel for convenient typing + + # If agent wants to put more than 100% into risky asset, he is constrained. + # Likewise if he wants to put less than 0% into risky asset, he is constrained. + constrained_top = FOC_s[:, -1] > 0.0 + constrained_bot = FOC_s[:, 0] < 0.0 + constrained = np.logical_or(constrained_top, constrained_bot) + a_idx = np.arange(aNrmCount) # For each value of aNrm, find the value of Share such that FOC_s == 0 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0) share_idx = np.argmax(crossing, axis=1) - # This represents the index of the segment of the share grid where dvds flips - # from positive to negative, indicating that there's a zero *on* the segment + + for k in range(3): + # This represents the index of the segment of the share grid where dvds flips + # from positive to negative, indicating that there's a zero *on* the segment. + # The exception is for aNrm values that are flagged as constrained, for which + # there will be no crossing point and we can just use the boundary value. + + # Now that we have a *range* for the location of the optimal share, we can + # do a refined search for the optimal share at each aNrm value where there + # is an interior solution (not constrained). We now make a refined ShareGrid + # that has *different* values on it for each aNrm value. + bot_s = ShareNext[a_idx,share_idx] + top_s = ShareNext[a_idx,share_idx + 1] + for j in range(aNrmCount): + if constrained[j]: + continue + ShareNext[j,:] = np.linspace(bot_s[j],top_s[j],ShareCount) + + # Now evaluate end-of-period marginal value of risky share on the refined grid + EndOfPrd_dvds = DiscFacEff * expected(calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)) + these = np.logical_not(constrained) + FOC_s[these,:] = EndOfPrd_dvds[these,:] # Relabel for convenient typing + + # Look for "crossing points" again + crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0) + share_idx[these] = np.argmax(crossing, axis=1) + + # Recalculate end-of-period marginal value of assets on the refined grid + EndOfPrd_dvda = DiscFacEff * expected(calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)) + EndOfPrd_dvdaNvrs[these,:] = uFunc.derinv(EndOfPrd_dvda[these,:]) # Calculate the fractional distance between those share gridpoints where the # zero should be found, assuming a linear function; call it alpha - a_idx = np.arange(aNrmCount) - bot_s = ShareGrid[share_idx] - top_s = ShareGrid[share_idx + 1] + bot_s = ShareNext[a_idx, share_idx] + top_s = ShareNext[a_idx, share_idx + 1] bot_f = FOC_s[a_idx, share_idx] top_f = FOC_s[a_idx, share_idx + 1] bot_c = EndOfPrd_dvdaNvrs[a_idx, share_idx] @@ -1276,7 +1312,7 @@ def calc_EndOfPrd_v(S, a, z): constrained_top = FOC_s[:, -1] > 0.0 constrained_bot = FOC_s[:, 0] < 0.0 - # Apply those constraints to both risky share and consumption (but lower + # Apply the constraints to both risky share and consumption (but lower # constraint should never be relevant) Share_now[constrained_top] = 1.0 Share_now[constrained_bot] = 0.0 From a8c7b244cc52b86c06506003786a5f88a937c7de Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 2 Apr 2024 16:51:14 -0400 Subject: [PATCH 4/8] Reformatting run --- HARK/ConsumptionSaving/ConsRiskyAssetModel.py | 57 ++++++++++--------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py index 8e9bde4e4..bd53fee12 100644 --- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py @@ -269,7 +269,9 @@ def temp_f(s): RiskyDstn.pmv, ) - SharePF = minimize_scalar(temp_f, bracket=(0.0, 1.0), method="golden", tol=1e-10).x[0] + SharePF = minimize_scalar( + temp_f, bracket=(0.0, 1.0), method="golden", tol=1e-10 + ).x[0] self.ShareLimit = SharePF self.add_to_time_inv("ShareLimit") @@ -1009,7 +1011,7 @@ def calc_hNrm(S): # This correctly accounts for risky returns and risk aversion hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj - + # This basic equation works if there's no correlation among shocks # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) @@ -1103,17 +1105,14 @@ def calc_EndOfPrd_dvds(R, a, z): bNrm_next = Rport * a # Calculate and return dvds (second term is all zeros) - EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed(bNrm_next) + EndOfPrd_dvds = Rxs * a * dvdbFunc_intermed(bNrm_next) + dvdsFunc_intermed( + bNrm_next + ) return EndOfPrd_dvds - # Evaluate realizations of value and marginal value after asset returns are realized - TempDstn = RiskyDstn + TempDstn = RiskyDstn # relabeling for below - # 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) + # Evaluate realizations of value and marginal value after asset returns are realized # Calculate end-of-period marginal value of risky portfolio share by taking expectations EndOfPrd_dvds = DiscFacEff * expected( @@ -1211,9 +1210,9 @@ def calc_EndOfPrd_v(S, a, z): v_next = vFunc_next(mNrm_next) EndOfPrd_v = (S["PermShk"] * PermGroFac) ** (1.0 - CRRA) * v_next return EndOfPrd_v - - calc_EndOfPrd_dvda = lambda S, a, z : calc_EndOfPrd_dvdx(S, a, z)[0] - calc_EndOfPrd_dvds = lambda S, a, z : calc_EndOfPrd_dvdx(S, a, z)[1] + + calc_EndOfPrd_dvda = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[0] + calc_EndOfPrd_dvds = lambda S, a, z: calc_EndOfPrd_dvdx(S, a, z)[1] TempDstn = ShockDstn # Evaluate realizations of value and marginal value after asset returns are realized @@ -1251,7 +1250,7 @@ def calc_EndOfPrd_v(S, a, z): # Now find the optimal (continuous) risky share on [0,1] by solving the first # order condition EndOfPrd_dvds == 0. FOC_s = EndOfPrd_dvds # Relabel for convenient typing - + # If agent wants to put more than 100% into risky asset, he is constrained. # Likewise if he wants to put less than 0% into risky asset, he is constrained. constrained_top = FOC_s[:, -1] > 0.0 @@ -1262,36 +1261,40 @@ def calc_EndOfPrd_v(S, a, z): # For each value of aNrm, find the value of Share such that FOC_s == 0 crossing = np.logical_and(FOC_s[:, 1:] <= 0.0, FOC_s[:, :-1] >= 0.0) share_idx = np.argmax(crossing, axis=1) - + for k in range(3): # This represents the index of the segment of the share grid where dvds flips # from positive to negative, indicating that there's a zero *on* the segment. # The exception is for aNrm values that are flagged as constrained, for which # there will be no crossing point and we can just use the boundary value. - + # Now that we have a *range* for the location of the optimal share, we can # do a refined search for the optimal share at each aNrm value where there # is an interior solution (not constrained). We now make a refined ShareGrid # that has *different* values on it for each aNrm value. - bot_s = ShareNext[a_idx,share_idx] - top_s = ShareNext[a_idx,share_idx + 1] + bot_s = ShareNext[a_idx, share_idx] + top_s = ShareNext[a_idx, share_idx + 1] for j in range(aNrmCount): if constrained[j]: continue - ShareNext[j,:] = np.linspace(bot_s[j],top_s[j],ShareCount) - + ShareNext[j, :] = np.linspace(bot_s[j], top_s[j], ShareCount) + # Now evaluate end-of-period marginal value of risky share on the refined grid - EndOfPrd_dvds = DiscFacEff * expected(calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext)) + EndOfPrd_dvds = DiscFacEff * expected( + calc_EndOfPrd_dvds, TempDstn, args=(aNrmNow, ShareNext) + ) these = np.logical_not(constrained) - FOC_s[these,:] = EndOfPrd_dvds[these,:] # Relabel for convenient typing - + FOC_s[these, :] = EndOfPrd_dvds[these, :] # Relabel for convenient typing + # Look for "crossing points" again crossing = np.logical_and(FOC_s[these, 1:] <= 0.0, FOC_s[these, :-1] >= 0.0) share_idx[these] = np.argmax(crossing, axis=1) - - # Recalculate end-of-period marginal value of assets on the refined grid - EndOfPrd_dvda = DiscFacEff * expected(calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext)) - EndOfPrd_dvdaNvrs[these,:] = uFunc.derinv(EndOfPrd_dvda[these,:]) + + # Calculate end-of-period marginal value of assets on the refined grid + EndOfPrd_dvda = DiscFacEff * expected( + calc_EndOfPrd_dvda, TempDstn, args=(aNrmNow, ShareNext) + ) + EndOfPrd_dvdaNvrs = uFunc.derinv(EndOfPrd_dvda) # Calculate the fractional distance between those share gridpoints where the # zero should be found, assuming a linear function; call it alpha From 10070379d244b42241bd82a34c50aa3094cdcfbb Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 2 Apr 2024 17:09:52 -0400 Subject: [PATCH 5/8] Fix indexing error I broke everything --- HARK/ConsumptionSaving/ConsRiskyAssetModel.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py index bd53fee12..aafd824f4 100644 --- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py +++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py @@ -271,7 +271,9 @@ def temp_f(s): SharePF = minimize_scalar( temp_f, bracket=(0.0, 1.0), method="golden", tol=1e-10 - ).x[0] + ).x + if type(SharePF) is np.array: + SharePF = SharePF[0] self.ShareLimit = SharePF self.add_to_time_inv("ShareLimit") From 3287467a2545ecf9a7ee033c8199c61206893545 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 2 Apr 2024 17:14:02 -0400 Subject: [PATCH 6/8] Missed two tabs --- HARK/ConsumptionSaving/ConsPortfolioModel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index e2006f5ec..fe8e04d35 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -430,7 +430,7 @@ def solve_one_period_ConsPortfolio( Risky_next = RiskyDstn.atoms RiskyMax = np.max(Risky_next) RiskyMin = np.min(Risky_next) - + # 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. @@ -454,7 +454,7 @@ def calc_hNrm(S): # This correctly accounts for risky returns and risk aversion hNrmNow = expected(calc_hNrm, ShockDstn) / R_adj - + # This basic equation works if there's no correlation among shocks # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) From 352e1f1b5d58d9293dd05d607e460b60766f6de4 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 9 Apr 2024 13:29:51 -0400 Subject: [PATCH 7/8] Correct human wealth calculation in ConsMarkovModel Incorporate proper discounting for human wealth into ConsMarkovModel so that extrapolation works properly. --- HARK/ConsumptionSaving/ConsMarkovModel.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py b/HARK/ConsumptionSaving/ConsMarkovModel.py index 9a5196113..586dfba27 100644 --- a/HARK/ConsumptionSaving/ConsMarkovModel.py +++ b/HARK/ConsumptionSaving/ConsMarkovModel.py @@ -381,9 +381,14 @@ def calc_vPPnext(S, a, R): MPCmaxEff = MPCmaxNow MPCmaxEff[BoroCnstNat_list < mNrmMin_list] = 1.0 - # Calculate the current Markov-state-conditional PDV of human wealth + # Calculate the current Markov-state-conditional PDV of human wealth, correctly + # accounting for risky returns and risk aversion hNrmPlusIncNext = Ex_IncNextAll + solution_next.hNrm - hNrmNow = np.dot(MrkvArray, (PermGroFac_list / Rfree_list) * hNrmPlusIncNext) + R_adj = np.dot(MrkvArray, Rfree_list ** (1.0 - CRRA)) + hNrmNow = ( + np.dot(MrkvArray, (PermGroFac_list / Rfree_list**CRRA) * hNrmPlusIncNext) + / R_adj + ) # Calculate the lower bound on MPC as m gets arbitrarily large temp = ( From 1a27df54053c99cf5053b53f99e8e22ab01aa2ed Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 9 Apr 2024 13:56:41 -0400 Subject: [PATCH 8/8] Edit CHANGELOG, fix one comment --- Documentation/CHANGELOG.md | 3 ++- HARK/ConsumptionSaving/ConsPortfolioModel.py | 4 +--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index ac280fc55..552d4b12a 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -22,8 +22,9 @@ Release Date: TBA - Add option to pass pre-built grid to `LinearFast`. [1388](https://github.com/econ-ark/HARK/pull/1388) - Moves calculation of stable points out of ConsIndShock solver, into method called by post_solve [#1349](https://github.com/econ-ark/HARK/pull/1349) - Adds cubic spline interpolation and value function construction to "warm glow bequest" models. -- Fixes cubic spline interpolation for ConsMedShockModel +- Fixes cubic spline interpolation for ConsMedShockModel. - Moves computation of "stable points" from inside of ConsIndShock solver to a post-solution method. [1349](https://github.com/econ-ark/HARK/pull/1349) +- Corrects calculation of "human wealth" under risky returns, providing correct limiting linear consumption function. [1403](https://github.com/econ-ark/HARK/pull/1403) ### 0.14.1 diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index fe8e04d35..a88be9f17 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -458,9 +458,7 @@ def calc_hNrm(S): # This basic equation works if there's no correlation among shocks # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) - # The above attempts to pin down the limiting consumption function for this - # model, however it is not clear why it creates bugs, so for now we allow - # for a linear extrapolation beyond the last asset point + # Set the terms of the limiting linear consumption function as mNrm goes to infinity cFuncLimitIntercept = MPCminNow * hNrmNow cFuncLimitSlope = MPCminNow