diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index dad06a173..39005d593 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -202,11 +202,52 @@ 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_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) + + 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 @@ -237,30 +278,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_patience_factor(Rfree, DiscFacEff, CRRA) + 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. 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 @@ -272,12 +312,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 @@ -313,11 +353,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 @@ -330,35 +370,157 @@ 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 +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 + worst_inc = np.min(income) + return np.sum(probs[income == worst_inc]) + + +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 + 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 +): + """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), + ) + + def solve_one_period_ConsIndShock( solution_next, IncShkDstn, @@ -372,8 +534,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 @@ -412,86 +573,59 @@ 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) + 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 = ((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 - # 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 - 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_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 - if BoroCnstNat < mNrmMinNow: - MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1 - else: - MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above + 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 # 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)) @@ -506,11 +640,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, MPCmaxUnc) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( @@ -545,9 +681,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) @@ -571,10 +711,14 @@ def calc_vPPnext(S, a, R): 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, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs + mNrm_temp, + vNvrs_temp, + vNvrsP_temp, + MPCminNvrs * hNrmNow, + MPCminNvrs, ) vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA) else: @@ -589,7 +733,7 @@ def calc_vPPnext(S, a, R): mNrmMin=mNrmMinNow, hNrm=hNrmNow, MPCmin=MPCminNow, - MPCmax=MPCmaxEff, + MPCmax=MPCmaxNow, ) return solution_now @@ -608,8 +752,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. @@ -655,6 +798,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 ( @@ -682,64 +826,48 @@ 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 + 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: - 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) - 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 - 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_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 - if BoroCnstNat < mNrmMinNow: - MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1 - else: - MPCmaxEff = MPCmaxNow # Otherwise, it's the MPC calculated above + 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 # 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 @@ -748,24 +876,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)) @@ -780,11 +897,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, MPCmaxUnc) # Construct the unconstrained consumption function as a cubic interpolation cFuncNowUnc = CubicInterp( @@ -826,9 +945,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) @@ -852,10 +975,14 @@ def calc_vPPnext(S, a, R): 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, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs + mNrm_temp, + vNvrs_temp, + vNvrsP_temp, + MPCminNvrs * hNrmNow, + MPCminNvrs, ) vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA) else: @@ -870,7 +997,7 @@ def calc_vPPnext(S, a, R): mNrmMin=mNrmMinNow, hNrm=hNrmNow, MPCmin=MPCminNow, - MPCmax=MPCmaxEff, + MPCmax=MPCmaxNow, ) return solution_now diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index f061f1952..26f00d4b9 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -304,6 +304,300 @@ def get_controls(self): self.controls["Share"] = ShareNow +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_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) + return hNrm + + +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, +): + """ + 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! + 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): + """ + 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 + 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): + """Compute end-of-period values.""" + # 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, @@ -429,26 +723,22 @@ 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_human_wealth, + ShockDstn, + args=(PermGroFac, ShareLimit, Rfree, CRRA, solution_next.hNrm), + ) + / R_adj + ) # This basic equation works if there's no correlation among shocks # hNrmNow = (PermGroFac/Rfree)*(1 + solution_next.hNrm) @@ -481,148 +771,68 @@ def calc_hNrm(S): 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 @@ -630,22 +840,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) @@ -662,76 +861,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) @@ -739,7 +886,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) @@ -869,17 +1027,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 @@ -937,13 +1084,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"], ) solution_now.hNrm = hNrmNow solution_now.MPCmin = MPCminNow diff --git a/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py b/HARK/ConsumptionSaving/ConsWealthPortfolioModel.py index 27a8c7488..bb4bfbb16 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 ( @@ -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