Skip to content

Commit

Permalink
Move addStablePoints from solver to post-solve
Browse files Browse the repository at this point in the history 
Calculating the "stable points" mNrmTrg and mNrmStE is a costly operation that is not useful in most contexts. This commit moves that code out of the solver and into a method on IndShockConsumerType, which is called as part of post_solve only if it is appropriate. The results get put into both self.bilt and self.solution (for compatibility with tests).
  • Loading branch information
@mnwhite
mnwhite committed Sep 15, 2023
1 parent dce9901 commit 7fa0abf
Showing 1 changed file with 56 additions and 203 deletions.
259 changes: 56 additions & 203 deletions HARK/ConsumptionSaving/ConsIndShockModel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
-------
Expand All @@ -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
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 7fa0abf

Please sign in to comment.