Coverage for HARK / ConsumptionSaving / ConsIndShockModelFast.py: 99%
178 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +0000
1"""
2Classes to solve canonical consumption-savings models with idiosyncratic shocks
3to income. All models here assume CRRA utility with geometric discounting, no
4bequest motive, and income shocks are fully transitory or fully permanent.
6It currently solves three types of models:
7 1) A very basic "perfect foresight" consumption-savings model with no uncertainty.
8 2) A consumption-savings model with risk over transitory and permanent income shocks.
9 3) The model described in (2), with an interest rate for debt that differs
10 from the interest rate for savings. #todo
12See NARK https://github.com/econ-ark/HARK/blob/master/docs/NARK/NARK.pdf for information on variable naming conventions.
13See HARK documentation for mathematical descriptions of the models being solved.
14"""
16import warnings
17from copy import deepcopy
19import numpy as np
20from interpolation import interp
21from numba import njit
22from quantecon.optimize import newton_secant
24from HARK import make_one_period_oo_solver
25from HARK.ConsumptionSaving.ConsIndShockModel import (
26 ConsumerSolution,
27 IndShockConsumerType,
28 PerfForesightConsumerType,
29 init_perfect_foresight,
30 init_idiosyncratic_shocks,
31)
32from HARK.ConsumptionSaving.LegacyOOsolvers import (
33 ConsIndShockSolverBasic,
34 ConsPerfForesightSolver,
35)
36from HARK.interpolation import (
37 CubicInterp,
38 LinearInterp,
39 LowerEnvelope,
40 MargValueFuncCRRA,
41 MargMargValueFuncCRRA,
42 ValueFuncCRRA,
43)
44from HARK.metric import MetricObject
45from HARK.numba_tools import (
46 CRRAutility,
47 CRRAutility_inv,
48 CRRAutility_invP,
49 CRRAutilityP,
50 CRRAutilityP_inv,
51 CRRAutilityP_invP,
52 CRRAutilityPP,
53 cubic_interp_fast,
54 linear_interp_deriv_fast,
55 linear_interp_fast,
56)
57from HARK.utilities import NullFunc
59__all__ = [
60 "PerfForesightSolution",
61 "IndShockSolution",
62 "PerfForesightConsumerTypeFast",
63 "IndShockConsumerTypeFast",
64]
66utility = CRRAutility
67utilityP = CRRAutilityP
68utilityPP = CRRAutilityPP
69utilityP_inv = CRRAutilityP_inv
70utility_invP = CRRAutility_invP
71utility_inv = CRRAutility_inv
72utilityP_invP = CRRAutilityP_invP
74# =====================================================================
75# === Terminal solution grid parameters for numba compatibility ===
76# =====================================================================
77# These parameters define the grid used to initialize vNvrs and vNvrsP arrays
78# in the terminal period solution. The grid needs to cover the range of mNrmNext
79# values that may be encountered during backward induction.
81# Minimum value for terminal grid (near-zero to avoid division issues)
82TERMINAL_GRID_MIN = 1e-6
84# Maximum value for terminal grid (should be large enough to cover typical
85# mNrmNext values during backward induction; 100 covers most standard cases)
86TERMINAL_GRID_MAX = 100.0
88# Number of points in the terminal grid (more points = better interpolation
89# accuracy but slightly higher memory usage)
90TERMINAL_GRID_SIZE = 200
93# =====================================================================
94# === Classes that help solve consumption-saving models ===
95# =====================================================================
98class PerfForesightSolution(MetricObject):
99 r"""
100 A class representing the solution of a single period of a consumption-saving
101 perfect foresight problem.
103 Here and elsewhere in the code, Nrm indicates that variables are normalized
104 by permanent income.
106 Parameters
107 ----------
108 mNrm: np.array
109 (Normalized) corresponding market resource points for interpolation.
110 cNrm : np.array
111 (Normalized) consumption points for interpolation.
112 vFuncNvrsSlope: float
113 Constant slope of inverse value vFuncNvrs
114 mNrmMin : float
115 The minimum allowable market resources for this period; the consump-
116 tion function (etc) are undefined for m < mNrmMin.
117 hNrm : float
118 Human wealth after receiving income this period: PDV of all future
119 income, ignoring mortality.
120 MPCmin : float
121 Infimum of the marginal propensity to consume this period.
122 MPC --> MPCmin as m --> infinity.
123 MPCmax : float
124 Supremum of the marginal propensity to consume this period.
125 MPC --> MPCmax as m --> mNrmMin.
126 """
128 distance_criteria = ["cNrm", "mNrm"]
130 def __init__(
131 self,
132 mNrm=np.array([0.0, 1.0]),
133 cNrm=np.array([0.0, 1.0]),
134 vFuncNvrsSlope=0.0,
135 mNrmMin=0.0,
136 hNrm=0.0,
137 MPCmin=1.0,
138 MPCmax=1.0,
139 ):
140 self.mNrm = mNrm
141 self.cNrm = cNrm
142 self.vFuncNvrsSlope = vFuncNvrsSlope
143 self.mNrmMin = mNrmMin
144 self.hNrm = hNrm
145 self.MPCmin = MPCmin
146 self.MPCmax = MPCmax
149class IndShockSolution(MetricObject):
150 """
151 A class representing the solution of a single period of a consumption-saving
152 idiosyncratic shocks to permanent and transitory income problem.
154 Parameters
155 ----------
156 mNrm: np.array
157 (Normalized) corresponding market resource points for interpolation.
158 cNrm : np.array
159 (Normalized) consumption points for interpolation.
160 vFuncNvrsSlope: float
161 Constant slope of inverse value ``vFuncNvrs``
162 mNrmMin : float
163 The minimum allowable market resources for this period; the consump-
164 tion function (etc) are undefined for m < mNrmMin.
165 hNrm : float
166 Human wealth after receiving income this period: PDV of all future
167 income, ignoring mortality.
168 MPCmin : float
169 Infimum of the marginal propensity to consume this period.
170 MPC --> MPCmin as m --> infinity.
171 MPCmax : float
172 Supremum of the marginal propensity to consume this period.
173 MPC --> MPCmax as m --> mNrmMin.
174 """
176 distance_criteria = ["cNrm", "mNrm", "mNrmMin"]
178 def __init__(
179 self,
180 mNrm=None,
181 cNrm=None,
182 cFuncLimitIntercept=None,
183 cFuncLimitSlope=None,
184 mNrmMin=0.0,
185 hNrm=0.0,
186 MPCmin=1.0,
187 MPCmax=1.0,
188 Ex_IncNext=0.0,
189 MPC=None,
190 mNrmGrid=None,
191 vNvrs=None,
192 vNvrsP=None,
193 MPCminNvrs=None,
194 ):
195 self.mNrm = mNrm if mNrm is not None else np.linspace(0, 1)
196 self.cNrm = cNrm if cNrm is not None else np.linspace(0, 1)
197 self.cFuncLimitIntercept = cFuncLimitIntercept
198 self.cFuncLimitSlope = cFuncLimitSlope
199 self.mNrmMin = mNrmMin
200 self.hNrm = hNrm
201 self.MPCmin = MPCmin
202 self.MPCmax = MPCmax
203 self.Ex_IncNext = Ex_IncNext
204 self.mNrmGrid = mNrmGrid
205 self.vNvrs = vNvrs
206 self.MPCminNvrs = MPCminNvrs
207 self.MPC = MPC
208 self.vNvrsP = vNvrsP
211# =====================================================================
212# === Classes and functions that solve consumption-saving models ===
213# =====================================================================
216def make_solution_terminal_fast(solution_terminal_class, CRRA):
217 """
218 Construct the terminal period solution for the fast solver.
220 At terminal period, consumer consumes everything: c = m.
221 Therefore (for CRRA != 1):
222 - v(m) = u(m)
223 - vNvrs(m) = u_inv(v(m)) = u_inv(u(m)) = m
224 - vNvrsP = d(vNvrs)/dm = 1
225 - MPCmin = 1 (consume everything)
226 - MPCminNvrs = 1 (since MPCmin = 1)
228 Note: This function requires CRRA != 1 because the vNvrs transformation
229 vNvrs(m) = u_inv(u(m)) = m only holds for CRRA utility. For log utility
230 (CRRA = 1), u(c) = log(c) and the inverse differs fundamentally.
232 Parameters
233 ----------
234 solution_terminal_class : class
235 The solution class to instantiate (PerfForesightSolution or IndShockSolution).
236 CRRA : float
237 Coefficient of relative risk aversion.
239 Returns
240 -------
241 solution_terminal : solution_terminal_class
242 The terminal period solution with properly initialized arrays for numba.
243 """
244 solution_terminal = solution_terminal_class()
246 # Terminal consumption function: c = m
247 cFunc_terminal = LinearInterp([0.0, 1.0], [0.0, 1.0])
248 solution_terminal.cFunc = cFunc_terminal
249 solution_terminal.vFunc = ValueFuncCRRA(cFunc_terminal, CRRA)
250 solution_terminal.vPfunc = MargValueFuncCRRA(cFunc_terminal, CRRA)
251 solution_terminal.vPPfunc = MargMargValueFuncCRRA(cFunc_terminal, CRRA)
253 # MPC is 1 everywhere at terminal (consume everything)
254 solution_terminal.MPC = np.array([1.0, 1.0])
256 # At terminal, MPCmin = 1 (consume everything), so MPCminNvrs = 1
257 solution_terminal.MPCminNvrs = 1.0
259 # Create grid that covers typical mNrmNext range during backward induction
260 # Uses module-level constants for configurability
261 mNrmGrid = np.linspace(TERMINAL_GRID_MIN, TERMINAL_GRID_MAX, TERMINAL_GRID_SIZE)
262 solution_terminal.mNrmGrid = mNrmGrid
264 # At terminal: vNvrs(m) = u_inv(u(m)) = m (since c = m, for CRRA != 1)
265 solution_terminal.vNvrs = mNrmGrid.copy()
267 # vNvrsP = d(vNvrs)/dm = d(m)/dm = 1 everywhere
268 solution_terminal.vNvrsP = np.ones_like(mNrmGrid)
270 # hNrm = 0 at terminal (no future income)
271 solution_terminal.hNrm = 0.0
273 return solution_terminal
276@njit(cache=True)
277def _find_mNrmStE(m, Rfree, PermGroFac, mNrm, cNrm, Ex_IncNext): # pragma: nocover
278 # Make a linear function of all combinations of c and m that yield mNext = mNow
279 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
281 # Find the steady state level of market resources
282 res = interp(mNrm, cNrm, m) - mZeroChange
283 # A zero of this is SS market resources
284 return res
287# @njit(cache=True) can't cache because of use of globals, perhaps newton_secant?
288@njit
289def _add_mNrmStENumba(
290 Rfree, PermGroFac, mNrm, cNrm, mNrmMin, Ex_IncNext, _find_mNrmStE
291): # pragma: nocover
292 """
293 Finds steady state (normalized) market resources and adds it to the
294 solution. This is the level of market resources such that the expectation
295 of market resources in the next period is unchanged. This value doesn't
296 necessarily exist.
297 """
299 # Minimum market resources plus next income is okay starting guess
300 m_init_guess = mNrmMin + Ex_IncNext
302 mNrmStE = newton_secant(
303 _find_mNrmStE,
304 m_init_guess,
305 args=(Rfree, PermGroFac, mNrm, cNrm, Ex_IncNext),
306 disp=False,
307 )
309 if mNrmStE.converged:
310 return mNrmStE.root
311 else:
312 return None
315@njit(cache=True, parallel=True)
316def _solveConsPerfForesightNumba(
317 DiscFac,
318 LivPrb,
319 CRRA,
320 Rfree,
321 PermGroFac,
322 BoroCnstArt,
323 MaxKinks,
324 mNrmNext,
325 cNrmNext,
326 hNrmNext,
327 MPCminNext,
328): # pragma: nocover
329 """
330 Makes the (linear) consumption function for this period.
331 """
333 DiscFacEff = DiscFac * LivPrb
335 # Calculate human wealth this period
336 hNrmNow = (PermGroFac / Rfree) * (hNrmNext + 1.0)
338 # Calculate the lower bound of the marginal propensity to consume
339 APF = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
340 MPCmin = 1.0 / (1.0 + APF / MPCminNext)
342 # Extract the discrete kink points in next period's consumption function;
343 # don't take the last one, as it only defines the extrapolation and is not a kink.
344 mNrmNext = mNrmNext[:-1]
345 cNrmNext = cNrmNext[:-1]
347 # Calculate the end-of-period asset values that would reach those kink points
348 # next period, then invert the first order condition to get consumption. Then
349 # find the endogenous gridpoint (kink point) today that corresponds to each kink
350 aNrmNow = (PermGroFac / Rfree) * (mNrmNext - 1.0)
351 cNrmNow = (DiscFacEff * Rfree) ** (-1.0 / CRRA) * (PermGroFac * cNrmNext)
352 mNrmNow = aNrmNow + cNrmNow
354 # Add an additional point to the list of gridpoints for the extrapolation,
355 # using the new value of the lower bound of the MPC.
356 mNrmNow = np.append(mNrmNow, mNrmNow[-1] + 1.0)
357 cNrmNow = np.append(cNrmNow, cNrmNow[-1] + MPCmin)
359 # If the artificial borrowing constraint binds, combine the constrained and
360 # unconstrained consumption functions.
361 if BoroCnstArt > mNrmNow[0]:
362 # Find the highest index where constraint binds
363 cNrmCnst = mNrmNow - BoroCnstArt
364 CnstBinds = cNrmCnst < cNrmNow
365 idx = np.where(CnstBinds)[0][-1]
367 if idx < (mNrmNow.size - 1):
368 # If it is not the *very last* index, find the the critical level
369 # of mNrm where the artificial borrowing contraint begins to bind.
370 d0 = cNrmNow[idx] - cNrmCnst[idx]
371 d1 = cNrmCnst[idx + 1] - cNrmNow[idx + 1]
372 m0 = mNrmNow[idx]
373 m1 = mNrmNow[idx + 1]
374 alpha = d0 / (d0 + d1)
375 mCrit = m0 + alpha * (m1 - m0)
377 # Adjust the grids of mNrm and cNrm to account for the borrowing constraint.
378 cCrit = mCrit - BoroCnstArt
379 mNrmNow = np.concatenate(
380 (np.array([BoroCnstArt, mCrit]), mNrmNow[(idx + 1) :])
381 )
382 cNrmNow = np.concatenate((np.array([0.0, cCrit]), cNrmNow[(idx + 1) :]))
384 else:
385 # If it *is* the very last index, then there are only three points
386 # that characterize the consumption function: the artificial borrowing
387 # constraint, the constraint kink, and the extrapolation point.
388 mXtra = (cNrmNow[-1] - cNrmCnst[-1]) / (1.0 - MPCmin)
389 mCrit = mNrmNow[-1] + mXtra
390 cCrit = mCrit - BoroCnstArt
391 mNrmNow = np.array([BoroCnstArt, mCrit, mCrit + 1.0])
392 cNrmNow = np.array([0.0, cCrit, cCrit + MPCmin])
394 # If the mNrm and cNrm grids have become too large, throw out the last
395 # kink point, being sure to adjust the extrapolation.
396 if mNrmNow.size > MaxKinks:
397 mNrmNow = np.concatenate((mNrmNow[:-2], np.array([mNrmNow[-3] + 1.0])))
398 cNrmNow = np.concatenate((cNrmNow[:-2], np.array([cNrmNow[-3] + MPCmin])))
400 # Calculate the upper bound of the MPC as the slope of the bottom segment.
401 MPCmax = (cNrmNow[1] - cNrmNow[0]) / (mNrmNow[1] - mNrmNow[0])
403 # Add attributes to enable calculation of steady state market resources.
404 # Relabeling for compatibility with add_mNrmStE
405 mNrmMinNow = mNrmNow[0]
407 # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations
408 vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
410 return (
411 mNrmNow,
412 cNrmNow,
413 vFuncNvrsSlope,
414 mNrmMinNow,
415 hNrmNow,
416 MPCmin,
417 MPCmax,
418 )
421class ConsPerfForesightSolverFast(ConsPerfForesightSolver):
422 """
423 A class for solving a one period perfect foresight consumption-saving problem.
424 An instance of this class is created by the function solvePerfForesight in each period.
425 """
427 def solve(self):
428 """
429 Solves the one period perfect foresight consumption-saving problem.
431 Parameters
432 ----------
433 None
435 Returns
436 -------
437 solution : PerfForesightSolution
438 The solution to this period's problem.
439 """
441 # Use a local value of BoroCnstArt to prevent comparing None and float below.
442 if self.BoroCnstArt is None:
443 BoroCnstArt = -np.inf
444 else:
445 BoroCnstArt = self.BoroCnstArt
447 (
448 self.mNrmNow,
449 self.cNrmNow,
450 self.vFuncNvrsSlope,
451 self.mNrmMinNow,
452 self.hNrmNow,
453 self.MPCmin,
454 self.MPCmax,
455 ) = _solveConsPerfForesightNumba(
456 self.DiscFac,
457 self.LivPrb,
458 self.CRRA,
459 self.Rfree,
460 self.PermGroFac,
461 BoroCnstArt,
462 self.MaxKinks,
463 self.solution_next.mNrm,
464 self.solution_next.cNrm,
465 self.solution_next.hNrm,
466 self.solution_next.MPCmin,
467 )
469 solution = PerfForesightSolution(
470 self.mNrmNow,
471 self.cNrmNow,
472 self.vFuncNvrsSlope,
473 self.mNrmMinNow,
474 self.hNrmNow,
475 self.MPCmin,
476 self.MPCmax,
477 )
478 return solution
481@njit(cache=True)
482def _np_tile(A, reps): # pragma: nocover
483 return A.repeat(reps[0]).reshape(A.size, -1).transpose()
486@njit(cache=True)
487def _np_insert(arr, obj, values, axis=-1): # pragma: nocover
488 return np.append(np.array(values), arr)
491@njit(cache=True, parallel=True)
492def _prepare_to_solveConsIndShockNumba(
493 DiscFac,
494 LivPrb,
495 CRRA,
496 Rfree,
497 PermGroFac,
498 BoroCnstArt,
499 aXtraGrid,
500 hNrmNext,
501 mNrmMinNext,
502 MPCminNext,
503 MPCmaxNext,
504 PermShkValsNext,
505 TranShkValsNext,
506 ShkPrbsNext,
507): # pragma: nocover
508 """
509 Unpacks some of the inputs (and calculates simple objects based on them),
510 storing the results in self for use by other methods. These include:
511 income shocks and probabilities, next period's marginal value function
512 (etc), the probability of getting the worst income shock next period,
513 the patience factor, human wealth, and the bounding MPCs.
515 Defines the constrained portion of the consumption function as cFuncNowCnst,
516 an attribute of self. Uses the artificial and natural borrowing constraints.
518 """
520 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
521 PermShkMinNext = np.min(PermShkValsNext)
522 TranShkMinNext = np.min(TranShkValsNext)
523 WorstIncPrb = np.sum(
524 ShkPrbsNext[
525 (PermShkValsNext * TranShkValsNext) == (PermShkMinNext * TranShkMinNext)
526 ]
527 )
529 # Update the bounding MPCs and PDV of human wealth:
530 APF = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
531 MPCminNow = 1.0 / (1.0 + APF / MPCminNext)
532 Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
533 hNrmNow = PermGroFac / Rfree * (Ex_IncNext + hNrmNext)
534 MPCmaxNow = 1.0 / (1.0 + (WorstIncPrb ** (1.0 / CRRA)) * APF / MPCmaxNext)
536 cFuncLimitIntercept = MPCminNow * hNrmNow
537 cFuncLimitSlope = MPCminNow
539 # Calculate the minimum allowable value of money resources in this period
540 BoroCnstNat = (mNrmMinNext - TranShkMinNext) * (PermGroFac * PermShkMinNext) / Rfree
542 # Note: need to be sure to handle BoroCnstArt==None appropriately.
543 # In Py2, this would evaluate to 5.0: np.max([None, 5.0]).
544 # However in Py3, this raises a TypeError. Thus here we need to directly
545 # address the situation in which BoroCnstArt == None:
546 if BoroCnstArt is None:
547 mNrmMinNow = BoroCnstNat
548 else:
549 mNrmMinNow = np.max(np.array([BoroCnstNat, BoroCnstArt]))
550 if BoroCnstNat < mNrmMinNow:
551 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
552 else:
553 MPCmaxEff = MPCmaxNow
555 """
556 Prepare to calculate end-of-period marginal value by creating an array
557 of market resources that the agent could have next period, considering
558 the grid of end-of-period assets and the distribution of shocks he might
559 experience next period.
560 """
562 # We define aNrmNow all the way from BoroCnstNat up to max(self.aXtraGrid)
563 # even if BoroCnstNat < BoroCnstArt, so we can construct the consumption
564 # function as the lower envelope of the (by the artificial borrowing con-
565 # straint) uconstrained consumption function, and the artificially con-
566 # strained consumption function.
567 aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
568 ShkCount = TranShkValsNext.size
569 aNrm_temp = _np_tile(aNrmNow, (ShkCount, 1))
571 # Tile arrays of the income shocks and put them into useful shapes
572 aNrmCount = aNrmNow.shape[0]
573 PermShkVals_temp = (_np_tile(PermShkValsNext, (aNrmCount, 1))).transpose()
574 TranShkVals_temp = (_np_tile(TranShkValsNext, (aNrmCount, 1))).transpose()
575 ShkPrbs_temp = (_np_tile(ShkPrbsNext, (aNrmCount, 1))).transpose()
577 # Get cash on hand next period
578 mNrmNext = Rfree / (PermGroFac * PermShkVals_temp) * aNrm_temp + TranShkVals_temp
579 # CDC 20191205: This should be divided by LivPrb[0] for Blanchard insurance
581 return (
582 DiscFacEff,
583 BoroCnstNat,
584 cFuncLimitIntercept,
585 cFuncLimitSlope,
586 mNrmMinNow,
587 hNrmNow,
588 MPCminNow,
589 MPCmaxNow,
590 MPCmaxEff,
591 Ex_IncNext,
592 mNrmNext,
593 PermShkVals_temp,
594 ShkPrbs_temp,
595 aNrmNow,
596 )
599@njit(cache=True, parallel=True)
600def _solveConsIndShockLinearNumba(
601 mNrmMinNext,
602 mNrmNext,
603 CRRA,
604 mNrmUnc,
605 cNrmUnc,
606 DiscFacEff,
607 Rfree,
608 PermGroFac,
609 PermShkVals_temp,
610 ShkPrbs_temp,
611 aNrmNow,
612 BoroCnstNat,
613 cFuncInterceptNext,
614 cFuncSlopeNext,
615): # pragma: nocover
616 """
617 Calculate end-of-period marginal value of assets at each point in aNrmNow.
618 Does so by taking a weighted sum of next period marginal values across
619 income shocks (in a preconstructed grid self.mNrmNext).
620 """
622 mNrmCnst = np.array([mNrmMinNext, mNrmMinNext + 1])
623 cNrmCnst = np.array([0.0, 1.0])
624 cFuncNextCnst = linear_interp_fast(mNrmNext.flatten(), mNrmCnst, cNrmCnst)
625 cFuncNextUnc = linear_interp_fast(
626 mNrmNext.flatten(), mNrmUnc, cNrmUnc, cFuncInterceptNext, cFuncSlopeNext
627 )
628 cFuncNext = np.minimum(cFuncNextCnst, cFuncNextUnc)
629 vPfuncNext = utilityP(cFuncNext, CRRA).reshape(mNrmNext.shape)
631 EndOfPrdvP = (
632 DiscFacEff
633 * Rfree
634 * PermGroFac ** (-CRRA)
635 * np.sum(PermShkVals_temp ** (-CRRA) * vPfuncNext * ShkPrbs_temp, axis=0)
636 )
638 # Finds interpolation points (c,m) for the consumption function.
640 cNrmNow = utilityP_inv(EndOfPrdvP, CRRA)
641 mNrmNow = cNrmNow + aNrmNow
643 # Limiting consumption is zero as m approaches mNrmMin
644 cNrm = _np_insert(cNrmNow, 0, 0.0, axis=-1)
645 mNrm = _np_insert(mNrmNow, 0, BoroCnstNat, axis=-1)
647 return (cNrm, mNrm, EndOfPrdvP)
650class ConsIndShockSolverBasicFast(ConsIndShockSolverBasic):
651 """
652 This class solves a single period of a standard consumption-saving problem,
653 using linear interpolation and without the ability to calculate the value
654 function. ConsIndShockSolver inherits from this class and adds the ability
655 to perform cubic interpolation and to calculate the value function.
657 Note that this class does not have its own initializing method. It initial-
658 izes the same problem in the same way as ConsIndShockSetup, from which it
659 inherits.
660 """
662 def prepare_to_solve(self):
663 """
664 Perform preparatory work before calculating the unconstrained consumption
665 function.
666 Parameters
667 ----------
668 none
669 Returns
670 -------
671 none
672 """
674 self.ShkPrbsNext = self.IncShkDstn.pmv
675 self.PermShkValsNext = self.IncShkDstn.atoms[0]
676 self.TranShkValsNext = self.IncShkDstn.atoms[1]
678 (
679 self.DiscFacEff,
680 self.BoroCnstNat,
681 self.cFuncLimitIntercept,
682 self.cFuncLimitSlope,
683 self.mNrmMinNow,
684 self.hNrmNow,
685 self.MPCminNow,
686 self.MPCmaxNow,
687 self.MPCmaxEff,
688 self.Ex_IncNext,
689 self.mNrmNext,
690 self.PermShkVals_temp,
691 self.ShkPrbs_temp,
692 self.aNrmNow,
693 ) = _prepare_to_solveConsIndShockNumba(
694 self.DiscFac,
695 self.LivPrb,
696 self.CRRA,
697 self.Rfree,
698 self.PermGroFac,
699 self.BoroCnstArt,
700 self.aXtraGrid,
701 self.solution_next.hNrm,
702 self.solution_next.mNrmMin,
703 self.solution_next.MPCmin,
704 self.solution_next.MPCmax,
705 self.PermShkValsNext,
706 self.TranShkValsNext,
707 self.ShkPrbsNext,
708 )
710 def solve(self):
711 """
712 Solves a one period consumption saving problem with risky income.
713 Parameters
714 ----------
715 None
716 Returns
717 -------
718 solution : ConsumerSolution
719 The solution to the one period problem.
720 """
722 self.cNrm, self.mNrm, self.EndOfPrdvP = _solveConsIndShockLinearNumba(
723 self.solution_next.mNrmMin,
724 self.mNrmNext,
725 self.CRRA,
726 self.solution_next.mNrm,
727 self.solution_next.cNrm,
728 self.DiscFacEff,
729 self.Rfree,
730 self.PermGroFac,
731 self.PermShkVals_temp,
732 self.ShkPrbs_temp,
733 self.aNrmNow,
734 self.BoroCnstNat,
735 self.solution_next.cFuncLimitIntercept,
736 self.solution_next.cFuncLimitSlope,
737 )
739 # Pack up the solution and return it
740 solution = IndShockSolution(
741 self.mNrm,
742 self.cNrm,
743 self.cFuncLimitIntercept,
744 self.cFuncLimitSlope,
745 self.mNrmMinNow,
746 self.hNrmNow,
747 self.MPCminNow,
748 self.MPCmaxEff,
749 self.Ex_IncNext,
750 )
752 return solution
755@njit(cache=True, parallel=True)
756def _solveConsIndShockCubicNumba(
757 mNrmMinNext,
758 mNrmNext,
759 mNrmUnc,
760 cNrmUnc,
761 MPCNext,
762 cFuncInterceptNext,
763 cFuncSlopeNext,
764 CRRA,
765 DiscFacEff,
766 Rfree,
767 PermGroFac,
768 PermShkVals_temp,
769 ShkPrbs_temp,
770 aNrmNow,
771 BoroCnstNat,
772 MPCmaxNow,
773): # pragma: nocover
774 mNrmCnst = np.array([mNrmMinNext, mNrmMinNext + 1])
775 cNrmCnst = np.array([0.0, 1.0])
776 cFuncNextCnst, MPCNextCnst = linear_interp_deriv_fast(
777 mNrmNext.flatten(), mNrmCnst, cNrmCnst
778 )
779 cFuncNextUnc, MPCNextUnc = cubic_interp_fast(
780 mNrmNext.flatten(),
781 mNrmUnc,
782 cNrmUnc,
783 MPCNext,
784 cFuncInterceptNext,
785 cFuncSlopeNext,
786 )
788 cFuncNext = np.where(cFuncNextCnst <= cFuncNextUnc, cFuncNextCnst, cFuncNextUnc)
790 vPfuncNext = utilityP(cFuncNext, CRRA).reshape(mNrmNext.shape)
792 EndOfPrdvP = (
793 DiscFacEff
794 * Rfree
795 * PermGroFac ** (-CRRA)
796 * np.sum(PermShkVals_temp ** (-CRRA) * vPfuncNext * ShkPrbs_temp, axis=0)
797 )
798 # Finds interpolation points (c,m) for the consumption function.
800 cNrmNow = EndOfPrdvP ** (-1.0 / CRRA)
801 mNrmNow = cNrmNow + aNrmNow
803 # Limiting consumption is zero as m approaches mNrmMin
804 cNrm = _np_insert(cNrmNow, 0, 0.0, axis=-1)
805 mNrm = _np_insert(mNrmNow, 0, BoroCnstNat, axis=-1)
807 """
808 Makes a cubic spline interpolation of the unconstrained consumption
809 function for this period.
810 """
812 MPCinterpNext = np.where(cFuncNextCnst <= cFuncNextUnc, MPCNextCnst, MPCNextUnc)
813 vPPfuncNext = (MPCinterpNext * utilityPP(cFuncNext, CRRA)).reshape(mNrmNext.shape)
815 EndOfPrdvPP = (
816 DiscFacEff
817 * Rfree
818 * Rfree
819 * PermGroFac ** (-CRRA - 1.0)
820 * np.sum(PermShkVals_temp ** (-CRRA - 1.0) * vPPfuncNext * ShkPrbs_temp, axis=0)
821 )
822 dcda = EndOfPrdvPP / utilityPP(cNrm[1:], CRRA)
823 MPC = dcda / (dcda + 1.0)
824 MPC = _np_insert(MPC, 0, MPCmaxNow)
826 return cNrm, mNrm, MPC, EndOfPrdvP
829@njit(cache=True)
830def _cFuncCubic(
831 aXtraGrid, mNrmMinNow, mNrmNow, cNrmNow, MPCNow, MPCminNow, hNrmNow
832): # pragma: nocover
833 mNrmGrid = mNrmMinNow + aXtraGrid
834 mNrmCnst = np.array([mNrmMinNow, mNrmMinNow + 1])
835 cNrmCnst = np.array([0.0, 1.0])
836 cFuncNowCnst = linear_interp_fast(mNrmGrid.flatten(), mNrmCnst, cNrmCnst)
837 cFuncNowUnc, MPCNowUnc = cubic_interp_fast(
838 mNrmGrid.flatten(), mNrmNow, cNrmNow, MPCNow, MPCminNow * hNrmNow, MPCminNow
839 )
841 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
843 return cNrmNow, mNrmGrid
846@njit(cache=True)
847def _cFuncLinear(
848 aXtraGrid, mNrmMinNow, mNrmNow, cNrmNow, MPCminNow, hNrmNow
849): # pragma: nocover
850 mNrmGrid = mNrmMinNow + aXtraGrid
851 mNrmCnst = np.array([mNrmMinNow, mNrmMinNow + 1])
852 cNrmCnst = np.array([0.0, 1.0])
853 cFuncNowCnst = linear_interp_fast(mNrmGrid.flatten(), mNrmCnst, cNrmCnst)
854 cFuncNowUnc = linear_interp_fast(
855 mNrmGrid.flatten(), mNrmNow, cNrmNow, MPCminNow * hNrmNow, MPCminNow
856 )
858 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
860 return cNrmNow, mNrmGrid
863@njit(cache=True)
864def _add_vFuncNumba(
865 mNrmNext,
866 mNrmGridNext,
867 vNvrsNext,
868 vNvrsPNext,
869 MPCminNvrsNext,
870 hNrmNext,
871 CRRA,
872 PermShkVals_temp,
873 PermGroFac,
874 DiscFacEff,
875 ShkPrbs_temp,
876 EndOfPrdvP,
877 aNrmNow,
878 BoroCnstNat,
879 mNrmGrid,
880 cFuncNow,
881 mNrmMinNow,
882 MPCmaxEff,
883 MPCminNow,
884): # pragma: nocover
885 """
886 Construct the end-of-period value function for this period, storing it
887 as an attribute of self for use by other methods.
888 """
890 # vFunc always cubic
892 vNvrsFuncNow, _ = cubic_interp_fast(
893 mNrmNext.flatten(),
894 mNrmGridNext,
895 vNvrsNext,
896 vNvrsPNext,
897 MPCminNvrsNext * hNrmNext,
898 MPCminNvrsNext,
899 )
901 vFuncNext = utility(vNvrsFuncNow, CRRA).reshape(mNrmNext.shape)
903 VLvlNext = (
904 PermShkVals_temp ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)
905 ) * vFuncNext
906 EndOfPrdv = DiscFacEff * np.sum(VLvlNext * ShkPrbs_temp, axis=0)
908 # value transformed through inverse utility
909 EndOfPrdvNvrs = utility_inv(EndOfPrdv, CRRA)
910 EndOfPrdvNvrsP = EndOfPrdvP * utility_invP(EndOfPrdv, CRRA)
911 EndOfPrdvNvrs = _np_insert(EndOfPrdvNvrs, 0, 0.0)
913 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
914 EndOfPrdvNvrsP = _np_insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
915 aNrm_temp = _np_insert(aNrmNow, 0, BoroCnstNat)
917 """
918 Creates the value function for this period, defined over market resources m.
919 self must have the attribute EndOfPrdvFunc in order to execute.
920 """
922 # Compute expected value and marginal value on a grid of market resources
924 aNrmNow = mNrmGrid - cFuncNow
926 EndOfPrdvNvrsFunc, _ = cubic_interp_fast(
927 aNrmNow, aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP
928 )
929 EndOfPrdvFunc = utility(EndOfPrdvNvrsFunc, CRRA)
931 vNrmNow = utility(cFuncNow, CRRA) + EndOfPrdvFunc
932 vPnow = utilityP(cFuncNow, CRRA)
934 # Construct the beginning-of-period value function
935 vNvrs = utility_inv(vNrmNow, CRRA) # value transformed through inverse utility
936 vNvrsP = vPnow * utility_invP(vNrmNow, CRRA)
937 mNrmGrid = _np_insert(mNrmGrid, 0, mNrmMinNow)
938 vNvrs = _np_insert(vNvrs, 0, 0.0)
939 vNvrsP = _np_insert(vNvrsP, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
940 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
942 return (
943 mNrmGrid,
944 vNvrs,
945 vNvrsP,
946 MPCminNvrs,
947 )
950@njit
951def _add_mNrmStEIndNumba(
952 PermGroFac,
953 Rfree,
954 Ex_IncNext,
955 mNrmMin,
956 mNrm,
957 cNrm,
958 MPC,
959 MPCmin,
960 hNrm,
961 _searchfunc,
962): # pragma: nocover
963 """
964 Finds steady state (normalized) market resources and adds it to the
965 solution. This is the level of market resources such that the expectation
966 of market resources in the next period is unchanged. This value doesn't
967 necessarily exist.
968 """
970 # Minimum market resources plus next income is okay starting guess
971 m_init_guess = mNrmMin + Ex_IncNext
973 mNrmStE = newton_secant(
974 _searchfunc,
975 m_init_guess,
976 args=(PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm),
977 disp=False,
978 )
980 if mNrmStE.converged:
981 return mNrmStE.root
982 else:
983 return None
986@njit(cache=True)
987def _find_mNrmStELinear(
988 m, PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm
989): # pragma: nocover
990 # Make a linear function of all combinations of c and m that yield mNext = mNow
991 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
993 mNrmCnst = np.array([mNrmMin, mNrmMin + 1])
994 cNrmCnst = np.array([0.0, 1.0])
995 cFuncNowCnst = linear_interp_fast(np.array([m]), mNrmCnst, cNrmCnst)
996 cFuncNowUnc = linear_interp_fast(np.array([m]), mNrm, cNrm, MPCmin * hNrm, MPCmin)
998 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
1000 # Find the steady state level of market resources
1001 res = cNrmNow[0] - mZeroChange
1002 # A zero of this is SS market resources
1003 return res
1006@njit(cache=True)
1007def _find_mNrmStECubic(
1008 m, PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm
1009): # pragma: nocover
1010 # Make a linear function of all combinations of c and m that yield mNext = mNow
1011 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
1013 mNrmCnst = np.array([mNrmMin, mNrmMin + 1])
1014 cNrmCnst = np.array([0.0, 1.0])
1015 cFuncNowCnst = linear_interp_fast(np.array([m]), mNrmCnst, cNrmCnst)
1016 cFuncNowUnc, MPCNowUnc = cubic_interp_fast(
1017 np.array([m]), mNrm, cNrm, MPC, MPCmin * hNrm, MPCmin
1018 )
1020 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
1022 # Find the steady state level of market resources
1023 res = cNrmNow[0] - mZeroChange
1024 # A zero of this is SS market resources
1025 return res
1028class ConsIndShockSolverFast(ConsIndShockSolverBasicFast):
1029 r"""
1030 This class solves a single period of a standard consumption-saving problem.
1031 It inherits from ConsIndShockSolverBasic, adding the ability to perform cubic
1032 interpolation and to calculate the value function.
1033 """
1035 def solve(self):
1036 """
1037 Solves a one period consumption saving problem with risky income.
1038 Parameters
1039 ----------
1040 None
1041 Returns
1042 -------
1043 solution : ConsumerSolution
1044 The solution to the one period problem.
1045 """
1047 if self.CubicBool:
1048 (
1049 self.cNrm,
1050 self.mNrm,
1051 self.MPC,
1052 self.EndOfPrdvP,
1053 ) = _solveConsIndShockCubicNumba(
1054 self.solution_next.mNrmMin,
1055 self.mNrmNext,
1056 self.solution_next.mNrm,
1057 self.solution_next.cNrm,
1058 self.solution_next.MPC,
1059 self.solution_next.cFuncLimitIntercept,
1060 self.solution_next.cFuncLimitSlope,
1061 self.CRRA,
1062 self.DiscFacEff,
1063 self.Rfree,
1064 self.PermGroFac,
1065 self.PermShkVals_temp,
1066 self.ShkPrbs_temp,
1067 self.aNrmNow,
1068 self.BoroCnstNat,
1069 self.MPCmaxNow,
1070 )
1071 # Pack up the solution and return it
1072 solution = IndShockSolution(
1073 self.mNrm,
1074 self.cNrm,
1075 self.cFuncLimitIntercept,
1076 self.cFuncLimitSlope,
1077 self.mNrmMinNow,
1078 self.hNrmNow,
1079 self.MPCminNow,
1080 self.MPCmaxEff,
1081 self.Ex_IncNext,
1082 self.MPC,
1083 )
1084 else:
1085 self.cNrm, self.mNrm, self.EndOfPrdvP = _solveConsIndShockLinearNumba(
1086 self.solution_next.mNrmMin,
1087 self.mNrmNext,
1088 self.CRRA,
1089 self.solution_next.mNrm,
1090 self.solution_next.cNrm,
1091 self.DiscFacEff,
1092 self.Rfree,
1093 self.PermGroFac,
1094 self.PermShkVals_temp,
1095 self.ShkPrbs_temp,
1096 self.aNrmNow,
1097 self.BoroCnstNat,
1098 self.solution_next.cFuncLimitIntercept,
1099 self.solution_next.cFuncLimitSlope,
1100 )
1102 # Pack up the solution and return it
1103 solution = IndShockSolution(
1104 self.mNrm,
1105 self.cNrm,
1106 self.cFuncLimitIntercept,
1107 self.cFuncLimitSlope,
1108 self.mNrmMinNow,
1109 self.hNrmNow,
1110 self.MPCminNow,
1111 self.MPCmaxEff,
1112 self.Ex_IncNext,
1113 )
1115 if self.vFuncBool:
1116 if self.CubicBool:
1117 self.cFuncNow, self.mNrmGrid = _cFuncCubic(
1118 self.aXtraGrid,
1119 self.mNrmMinNow,
1120 self.mNrm,
1121 self.cNrm,
1122 self.MPC,
1123 self.MPCminNow,
1124 self.hNrmNow,
1125 )
1126 else:
1127 self.cFuncNow, self.mNrmGrid = _cFuncLinear(
1128 self.aXtraGrid,
1129 self.mNrmMinNow,
1130 self.mNrm,
1131 self.cNrm,
1132 self.MPCminNow,
1133 self.hNrmNow,
1134 )
1136 self.mNrmGrid, self.vNvrs, self.vNvrsP, self.MPCminNvrs = _add_vFuncNumba(
1137 self.mNrmNext,
1138 self.solution_next.mNrmGrid,
1139 self.solution_next.vNvrs,
1140 self.solution_next.vNvrsP,
1141 self.solution_next.MPCminNvrs,
1142 self.solution_next.hNrm,
1143 self.CRRA,
1144 self.PermShkVals_temp,
1145 self.PermGroFac,
1146 self.DiscFacEff,
1147 self.ShkPrbs_temp,
1148 self.EndOfPrdvP,
1149 self.aNrmNow,
1150 self.BoroCnstNat,
1151 self.mNrmGrid,
1152 self.cFuncNow,
1153 self.mNrmMinNow,
1154 self.MPCmaxEff,
1155 self.MPCminNow,
1156 )
1158 # Pack up the solution and return it
1160 solution.mNrmGrid = self.mNrmGrid
1161 solution.vNvrs = self.vNvrs
1162 solution.vNvrsP = self.vNvrsP
1163 solution.MPCminNvrs = self.MPCminNvrs
1165 return solution
1168# ============================================================================
1169# == Classes for representing types of consumer agents (and things they do) ==
1170# ============================================================================
1173init_perfect_foresight_fast = init_perfect_foresight.copy()
1174perf_foresight_constructor_dict = init_perfect_foresight["constructors"].copy()
1175perf_foresight_constructor_dict["solution_terminal"] = make_solution_terminal_fast
1176init_perfect_foresight_fast["constructors"] = perf_foresight_constructor_dict
1179class PerfForesightConsumerTypeFast(PerfForesightConsumerType):
1180 r"""
1181 A version of the perfect foresight consumer type speed up by numba.
1183 Note: This fast solver does not support CRRA=1 (log utility) due to the
1184 mathematical singularity in the inverse value function transformation.
1185 Use the standard PerfForesightConsumerType for log utility.
1186 """
1188 solution_terminal_class = PerfForesightSolution
1189 default_ = {
1190 "params": init_perfect_foresight_fast,
1191 "solver": make_one_period_oo_solver(ConsPerfForesightSolverFast),
1192 "model": "ConsPerfForesight.yaml",
1193 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
1194 }
1196 def pre_solve(self):
1197 """
1198 Perform pre-solve checks and setup.
1200 Raises
1201 ------
1202 ValueError
1203 If CRRA equals 1 (log utility), which is not supported by the fast solver.
1205 Warns
1206 -----
1207 UserWarning
1208 If CRRA is very close to 1 (between 0.99 and 1.01), which may cause
1209 numerical instability.
1210 """
1211 if np.isclose(self.CRRA, 1.0):
1212 raise ValueError(
1213 "PerfForesightConsumerTypeFast does not support CRRA=1 (log utility) "
1214 "due to mathematical singularities in the numba-optimized solver. "
1215 "Please use PerfForesightConsumerType instead for log utility preferences."
1216 )
1217 # Warn for CRRA values that are close to 1 but not caught by np.isclose
1218 if 0.99 < self.CRRA < 1.01 and not np.isclose(self.CRRA, 1.0):
1219 warnings.warn(
1220 f"CRRA={self.CRRA} is very close to 1, which may cause numerical "
1221 "instability. Consider using the standard solver or a CRRA value "
1222 "further from 1.",
1223 UserWarning,
1224 stacklevel=2,
1225 )
1226 # Call parent's pre_solve
1227 super().pre_solve()
1229 def post_solve(self):
1230 self.solution_fast = deepcopy(self.solution)
1232 if self.cycles == 0:
1233 terminal = 1
1234 else:
1235 terminal = self.cycles
1236 self.solution[terminal] = self.solution_terminal
1238 for i in range(terminal):
1239 solution = self.solution[i]
1241 # Construct the consumption function as a linear interpolation.
1242 cFunc = LinearInterp(solution.mNrm, solution.cNrm)
1244 """
1245 Defines the value and marginal value functions for this period.
1246 Uses the fact that for a perfect foresight CRRA utility problem,
1247 if the MPC in period t is :math:`\\kappa_{t}`, and relative risk
1248 aversion :math:`\rho`, then the inverse value vFuncNvrs has a
1249 constant slope of :math:`\\kappa_{t}^{-\rho/(1-\rho)}` and
1250 vFuncNvrs has value of zero at the lower bound of market resources
1251 mNrmMin. See PerfForesightConsumerType.ipynb documentation notebook
1252 for a brief explanation and the links below for a fuller treatment.
1254 https://www.econ2.jhu.edu/people/ccarroll/public/lecturenotes/consumption/PerfForesightCRRA/#vFuncAnalytical
1255 https://www.econ2.jhu.edu/people/ccarroll/SolvingMicroDSOPs/#vFuncPF
1256 """
1258 vFuncNvrs = LinearInterp(
1259 np.array([solution.mNrmMin, solution.mNrmMin + 1.0]),
1260 np.array([0.0, solution.vFuncNvrsSlope]),
1261 )
1262 vFunc = ValueFuncCRRA(vFuncNvrs, self.CRRA)
1263 vPfunc = MargValueFuncCRRA(cFunc, self.CRRA)
1265 consumer_solution = ConsumerSolution(
1266 cFunc=cFunc,
1267 vFunc=vFunc,
1268 vPfunc=vPfunc,
1269 mNrmMin=solution.mNrmMin,
1270 hNrm=solution.hNrm,
1271 MPCmin=solution.MPCmin,
1272 MPCmax=solution.MPCmax,
1273 )
1275 Ex_IncNext = 1.0 # Perfect foresight income of 1
1277 # Add mNrmStE to the solution and return it
1278 consumer_solution.mNrmStE = _add_mNrmStENumba(
1279 self.Rfree[i],
1280 self.PermGroFac[i],
1281 solution.mNrm,
1282 solution.cNrm,
1283 solution.mNrmMin,
1284 Ex_IncNext,
1285 _find_mNrmStE,
1286 )
1288 self.solution[i] = consumer_solution
1291###############################################################################
1294def select_fast_solver(CubicBool, vFuncBool):
1295 if (not CubicBool) and (not vFuncBool):
1296 solver = ConsIndShockSolverBasicFast
1297 else: # Use the "advanced" solver if either is requested
1298 solver = ConsIndShockSolverFast
1299 solve_one_period = make_one_period_oo_solver(solver)
1300 return solve_one_period
1303init_idiosyncratic_shocks_fast = init_idiosyncratic_shocks.copy()
1304ind_shock_fast_constructor_dict = init_idiosyncratic_shocks["constructors"].copy()
1305ind_shock_fast_constructor_dict["solution_terminal"] = make_solution_terminal_fast
1306ind_shock_fast_constructor_dict["solve_one_period"] = select_fast_solver
1307init_idiosyncratic_shocks_fast["constructors"] = ind_shock_fast_constructor_dict
1310class IndShockConsumerTypeFast(IndShockConsumerType, PerfForesightConsumerTypeFast):
1311 r"""
1312 A version of the idiosyncratic shock consumer type sped up by numba.
1314 If CubicBool and vFuncBool are both set to false it's further optimized.
1316 Note: This fast solver does not support CRRA=1 (log utility) due to the
1317 mathematical singularity in the inverse value function transformation.
1318 Use the standard IndShockConsumerType for log utility.
1319 """
1321 solution_terminal_class = IndShockSolution
1322 default_ = {
1323 "params": init_idiosyncratic_shocks_fast,
1324 "solver": NullFunc(),
1325 "model": "ConsIndShock.yaml",
1326 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
1327 }
1329 def pre_solve(self):
1330 """
1331 Perform pre-solve checks and setup.
1333 Raises
1334 ------
1335 ValueError
1336 If CRRA equals 1 (log utility), which is not supported by the fast solver.
1338 Warns
1339 -----
1340 UserWarning
1341 If CRRA is very close to 1 (between 0.99 and 1.01), which may cause
1342 numerical instability.
1343 """
1344 if np.isclose(self.CRRA, 1.0):
1345 raise ValueError(
1346 "IndShockConsumerTypeFast does not support CRRA=1 (log utility) "
1347 "due to mathematical singularities in the numba-optimized solver. "
1348 "Please use IndShockConsumerType instead for log utility preferences."
1349 )
1350 # Warn for CRRA values that are close to 1 but not caught by np.isclose
1351 if 0.99 < self.CRRA < 1.01 and not np.isclose(self.CRRA, 1.0):
1352 warnings.warn(
1353 f"CRRA={self.CRRA} is very close to 1, which may cause numerical "
1354 "instability. Consider using the standard solver or a CRRA value "
1355 "further from 1.",
1356 UserWarning,
1357 stacklevel=2,
1358 )
1359 # Call parent's pre_solve
1360 super().pre_solve()
1362 def post_solve(self):
1363 self.solution_fast = deepcopy(self.solution)
1365 if self.cycles == 0:
1366 cycles = 1
1367 else:
1368 cycles = self.cycles
1369 self.solution[-1] = init_idiosyncratic_shocks["constructors"][
1370 "solution_terminal"
1371 ](self.CRRA)
1373 for i in range(cycles):
1374 for j in range(self.T_cycle):
1375 solution = self.solution[i * self.T_cycle + j]
1377 # Define the borrowing constraint (limiting consumption function)
1378 cFuncNowCnst = LinearInterp(
1379 np.array([solution.mNrmMin, solution.mNrmMin + 1]),
1380 np.array([0.0, 1.0]),
1381 )
1383 """
1384 Constructs a basic solution for this period, including the consumption
1385 function and marginal value function.
1386 """
1388 if self.CubicBool:
1389 # Makes a cubic spline interpolation of the unconstrained consumption
1390 # function for this period.
1391 cFuncNowUnc = CubicInterp(
1392 solution.mNrm,
1393 solution.cNrm,
1394 solution.MPC,
1395 solution.cFuncLimitIntercept,
1396 solution.cFuncLimitSlope,
1397 )
1398 else:
1399 # Makes a linear interpolation to represent the (unconstrained) consumption function.
1400 # Construct the unconstrained consumption function
1401 cFuncNowUnc = LinearInterp(
1402 solution.mNrm,
1403 solution.cNrm,
1404 solution.cFuncLimitIntercept,
1405 solution.cFuncLimitSlope,
1406 )
1408 # Combine the constrained and unconstrained functions into the true consumption function
1409 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst)
1411 # Make the marginal value function and the marginal marginal value function
1412 vPfuncNow = MargValueFuncCRRA(cFuncNow, self.CRRA)
1414 # Pack up the solution and return it
1415 consumer_solution = ConsumerSolution(
1416 cFunc=cFuncNow,
1417 vPfunc=vPfuncNow,
1418 mNrmMin=solution.mNrmMin,
1419 hNrm=solution.hNrm,
1420 MPCmin=solution.MPCmin,
1421 MPCmax=solution.MPCmax,
1422 )
1424 if self.vFuncBool:
1425 vNvrsFuncNow = CubicInterp(
1426 solution.mNrmGrid,
1427 solution.vNvrs,
1428 solution.vNvrsP,
1429 solution.MPCminNvrs * solution.hNrm,
1430 solution.MPCminNvrs,
1431 )
1432 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, self.CRRA)
1434 consumer_solution.vFunc = vFuncNow
1436 if self.CubicBool or self.vFuncBool:
1437 _searchFunc = (
1438 _find_mNrmStECubic if self.CubicBool else _find_mNrmStELinear
1439 )
1440 # Add mNrmStE to the solution and return it
1441 consumer_solution.mNrmStE = _add_mNrmStEIndNumba(
1442 self.PermGroFac[j],
1443 self.Rfree[j],
1444 solution.Ex_IncNext,
1445 solution.mNrmMin,
1446 solution.mNrm,
1447 solution.cNrm,
1448 solution.MPC,
1449 solution.MPCmin,
1450 solution.hNrm,
1451 _searchFunc,
1452 )
1454 self.solution[i * self.T_cycle + j] = consumer_solution
1456 if (self.cycles == 0) and (self.T_cycle == 1):
1457 self.calc_stable_points(force=True)