diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 37ceaff43..d9d2a8e0f 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -397,158 +397,8 @@ def make_cFunc_PF(self): # Add two attributes to enable calculation of steady state market resources. self.Ex_IncNext = 1.0 # Perfect foresight income of 1 - # Relabeling for compatibility with add_mNrmStE self.mNrmMinNow = mNrmNow[0] - def add_mNrmTrg(self, solution): - """ - Finds value of (normalized) market resources m at which individual consumer - expects m not to change. - This will exist if the GICNrm holds. - - https://econ-ark.github.io/BufferStockTheory#Unique-Stable-Points - - Parameters - ---------- - solution : ConsumerSolution - Solution to this period's problem, which must have attribute cFunc. - Returns - ------- - solution : ConsumerSolution - Same solution that was passed, but now with the attribute mNrmStE. - """ - - # If no uncertainty, return the degenerate targets for the PF model - if hasattr(self, "TranShkMinNext"): # Then it has transitory shocks - # Handle the degenerate case where shocks are of size zero - if (self.TranShkMinNext == 1.0) and (self.PermShkMinNext == 1.0): - # but they are of zero size (and also permanent are zero) - if self.GICRaw: # max of nat and art boro cnst - if type(self.BoroCnstArt) == type(None): - solution.mNrmStE = -self.hNrmNow - solution.mNrmTrg = -self.hNrmNow - else: - bNrmNxt = -self.BoroCnstArt * self.Rfree / self.PermGroFac - solution.mNrmStE = bNrmNxt + 1.0 - solution.mNrmTrg = bNrmNxt + 1.0 - else: # infinity - solution.mNrmStE = float("inf") - solution.mNrmTrg = float("inf") - return solution - - # First find - # \bar{\mathcal{R}} = E_t[R/Gamma_{t+1}] = R/Gamma E_t[1/psi_{t+1}] - if type(self) == ConsPerfForesightSolver: - Ex_PermShkInv = 1.0 - else: - Ex_PermShkInv = np.dot(1 / self.PermShkValsNext, self.ShkPrbsNext) - - Ex_RNrmFac = (self.Rfree / self.PermGroFac) * Ex_PermShkInv - - # mNrmTrg solves Rcalbar*(m - c(m)) + E[inc_next] = m. Define a - # rearranged version. - def Ex_m_tp1_minus_m_t(m): - return Ex_RNrmFac * (m - solution.cFunc(m)) + self.Ex_IncNext - m - - # Minimum market resources plus next income is okay starting guess - m_init_guess = self.mNrmMinNow + self.Ex_IncNext - try: - mNrmTrg = newton(Ex_m_tp1_minus_m_t, m_init_guess) - except: - mNrmTrg = None - - # Add mNrmTrg to the solution and return it - solution.mNrmTrg = mNrmTrg - return solution - - def add_mNrmStE(self, solution): - """ - Finds market resources ratio at which 'balanced growth' is expected. - This is the m ratio such that the expected growth rate of the M level - matches the expected growth rate of permanent income. This value does - not exist if the Growth Impatience Condition does not hold. - - https://econ-ark.github.io/BufferStockTheory#Unique-Stable-Points - - Parameters - ---------- - solution : ConsumerSolution - Solution to this period's problem, which must have attribute cFunc. - Returns - ------- - solution : ConsumerSolution - Same solution that was passed, but now with the attribute mNrmStE - """ - # Probably should test whether GICRaw holds and log error if it does not - # using check_conditions - # All combinations of c and m that yield E[PermGroFac PermShkVal mNext] = mNow - # https://econ-ark.github.io/BufferStockTheory/#The-Individual-Steady-State - - PF_RNrm = self.Rfree / self.PermGroFac - # If we are working with a model that permits uncertainty but that - # uncertainty has been set to zero, return the correct answer - # by hand because in this degenerate case numerical search may - # have trouble - if hasattr(self, "TranShkMinNext"): # Then it has transitory shocks - if (self.TranShkMinNext == 1.0) and (self.PermShkMinNext == 1.0): - # but they are of zero size (and permanent shocks also not there) - if self.GICRaw: # max of nat and art boro cnst - # breakpoint() - if type(self.BoroCnstArt) == type(None): - solution.mNrmStE = -self.hNrmNow - solution.mNrmTrg = -self.hNrmNow - else: - bNrmNxt = -self.BoroCnstArt * self.Rfree / self.PermGroFac - solution.mNrmStE = bNrmNxt + 1.0 - solution.mNrmTrg = bNrmNxt + 1.0 - else: # infinity - solution.mNrmStE = float("inf") - solution.mNrmTrg = float("inf") - return solution - - def Ex_PermShk_tp1_times_m_tp1_minus_m_t(mStE): - return PF_RNrm * (mStE - solution.cFunc(mStE)) + 1.0 - mStE - - # Minimum market resources plus next income is okay starting guess - m_init_guess = self.mNrmMinNow + self.Ex_IncNext - try: - mNrmStE = newton(Ex_PermShk_tp1_times_m_tp1_minus_m_t, m_init_guess) - except: - mNrmStE = None - - solution.mNrmStE = mNrmStE - return solution - - def add_stable_points(self, solution): - """ - Checks necessary conditions for the existence of the individual steady - state and target levels of market resources (see above). - If the conditions are satisfied, computes and adds the stable points - to the solution. - - Parameters - ---------- - solution : ConsumerSolution - Solution to this period's problem, which must have attribute cFunc. - - Returns - ------- - solution : ConsumerSolution - Same solution that was provided, augmented with attributes mNrmStE and - mNrmTrg, if they exist. - - """ - - # 0. There is no non-degenerate steady state for any unconstrained PF model. - # 1. There is a non-degenerate SS for constrained PF model if GICRaw holds. - # Therefore - # Check if (GICRaw and BoroCnstArt) and if so compute them both - thorn = (self.Rfree * self.DiscFacEff) ** (1 / self.CRRA) - GICRaw = 1 > thorn / self.PermGroFac - if self.BoroCnstArt is not None and GICRaw: - solution = self.add_mNrmStE(solution) - solution = self.add_mNrmTrg(solution) - return solution def solve(self): """ @@ -993,51 +843,7 @@ def add_MPC_and_human_wealth(self, solution): solution.MPCmin = self.MPCminNow solution.MPCmax = self.MPCmaxEff return solution - - def add_stable_points(self, solution): - """ - Checks necessary conditions for the existence of the individual steady - state and target levels of market resources (see above). - If the conditions are satisfied, computes and adds the stable points - to the solution. - - Parameters - ---------- - solution : ConsumerSolution - Solution to this period's problem, which must have attribute cFunc. - - Returns - ------- - solution : ConsumerSolution - Same solution that was passed, but now with attributes mNrmStE and - mNrmTrg, if they exist. - - """ - - # 0. Check if GICRaw holds. If so, then mNrmStE will exist. So, compute it. - # 1. Check if GICNrm holds. If so, then mNrmTrg will exist. So, compute it. - - thorn = (self.Rfree * self.DiscFacEff) ** (1 / self.CRRA) - - GPFRaw = thorn / self.PermGroFac - self.GPFRaw = GPFRaw - GPFNrm = ( - thorn / self.PermGroFac / np.dot(1 / self.PermShkValsNext, self.ShkPrbsNext) - ) - self.GPFNrm = GPFNrm - GICRaw = 1 > thorn / self.PermGroFac - self.GICRaw = GICRaw - GICNrm = 1 > GPFNrm - self.GICNrm = GICNrm - - if GICRaw: - # find steady state m, if it exists - solution = self.add_mNrmStE(solution) - if GICNrm: - # find target m, if it exists - solution = self.add_mNrmTrg(solution) - - return solution + def make_linear_cFunc(self, mNrm, cNrm): """ @@ -1077,7 +883,6 @@ def solve(self): EndOfPrdvP = self.calc_EndOfPrdvP() solution = self.make_basic_solution(EndOfPrdvP, aNrmNow, self.make_linear_cFunc) solution = self.add_MPC_and_human_wealth(solution) - solution = self.add_stable_points(solution) return solution @@ -1286,7 +1091,6 @@ def solve(self): ) solution = self.add_MPC_and_human_wealth(solution) # add a few things - solution = self.add_stable_points(solution) # Add the value function if requested, as well as the marginal marginal # value function if cubic splines were used (to prepare for next period) @@ -2049,7 +1853,7 @@ def calc_limiting_values(self): solution to an infinite horizon problem. This method should only be called when T_cycle=1 and cycles=0, otherwise the values generated are meaningless. This method adds the following values to the instance in the dictionary - attribute auxiliary. + attribute called bilt. APFac : Absolute Patience Factor GPFacRaw : Growth Patience Factor @@ -3261,7 +3065,7 @@ def calc_limiting_values(self): solution to an infinite horizon problem. This method should only be called when T_cycle=1 and cycles=0, otherwise the values generated are meaningless. This method adds the following values to this instance in the dictionary - attribute auxiliary. + attribute called bilt. APFac : Absolute Patience Factor GPFacRaw : Growth Patience Factor @@ -3280,6 +3084,8 @@ def calc_limiting_values(self): hNrm : Human wealth divided by permanent income. ELogPermShk : Expected log permanent income shock WorstPrb : Probability of worst income shock realization + Delta_mNrm_ZeroFunc : Linear consumption function where expected change in market resource ratio is zero + BalGroFunc : Linear consumption function where the level of market resources grows at the same rate as permanent income Returns ------- @@ -3291,7 +3097,8 @@ def calc_limiting_values(self): # Calculate the risk-modified growth impatience factor PermShkDstn = self.PermShkDstn[0] inv_func = lambda x : x**(-1.) - GroCompPermShk = expected(inv_func, PermShkDstn)[0]**(-1.) + Ex_PermShkInv = expected(inv_func, PermShkDstn)[0] + GroCompPermShk = Ex_PermShkInv**(-1.) aux_dict['GPFacMod'] = aux_dict['APFac'] / (self.PermGroFac[0] * GroCompPermShk) # Calculate the mortality-adjusted growth impatience factor (and version @@ -3353,6 +3160,14 @@ def calc_limiting_values(self): aux_dict['hNrm'] = hNrm aux_dict['MPCmax'] = MPCmax + # Generate the "Delta m = 0" function, which is used to find target market resources + Ex_Rnrm = self.Rfree / self.PermGroFac[0] * Ex_PermShkInv + aux_dict['Delta_mNrm_ZeroFunc'] = lambda m : (1. - 1./Ex_Rnrm) * m + 1./Ex_Rnrm + + # Generate the "E[M_tp1 / M_t] = G" function, which is used to find balanced growth market resources + PF_Rnrm = self.Rfree / self.PermGroFac[0] + aux_dict['BalGroFunc'] = lambda m : (1. - 1./PF_Rnrm) * m + 1./PF_Rnrm + self.bilt = aux_dict @@ -3591,11 +3406,49 @@ def calc_stable_points(self): None """ infinite_horizon = self.cycles == 0 + single_period = self.T_cycle = 1 if not infinite_horizon: - _log.warning( - "The calc_stable_points method works only for infinite horizon models." - ) + _log.warning("The calc_stable_points method works only for infinite horizon models.") + return + if not single_period: + _log.warning("The calc_stable_points method works only with a single infinitely repeated period.") + return + if not hasattr(self, 'conditions'): + _log.warning("The calc_limiting_values method must be run before the calc_stable_points method.") return + if not hasattr(self, 'solution'): + _log.warning("The solve method must be run before the calc_stable_points method.") + return + + # Extract balanced growth and delta m_t+1 = 0 functions + BalGroFunc = self.bilt['BalGroFunc'] + Delta_mNrm_ZeroFunc = self.bilt['Delta_mNrm_ZeroFunc'] + + # If the GICRaw holds, then there is a balanced growth market resources ratio + if self.conditions['GICRaw']: + cFunc = self.solution[0].cFunc + func_to_zero = lambda m : BalGroFunc(m) - cFunc(m) + m0 = 1.0 + try: + mNrmStE = newton(func_to_zero, m0) + except: + mNrmStE = np.nan + + # A target level of assets *might* exist even if the GICMod fails, so check no matter what + func_to_zero = lambda m : Delta_mNrm_ZeroFunc(m) - cFunc(m) + m0 = 1.0 if np.isnan(mNrmStE) else mNrmStE + try: + mNrmTrg = newton(func_to_zero, m0, maxiter=200) + except: + mNrmTrg = np.nan + else: + mNrmStE = np.nan + mNrmTrg = np.nan + + self.solution[0].mNrmStE = mNrmStE + self.solution[0].mNrmTrg = mNrmTrg + self.bilt['mNrmStE'] = mNrmStE + self.bilt['mNrmTrg'] = mNrmTrg # = Functions for generating discrete income processes and