Coverage for HARK/ConsumptionSaving/ConsIndShockModelFast.py: 99%
160 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +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"""
16from copy import deepcopy
18import numpy as np
19from interpolation import interp
20from numba import njit
21from quantecon.optimize import newton_secant
23from HARK import make_one_period_oo_solver
24from HARK.ConsumptionSaving.ConsIndShockModel import (
25 ConsumerSolution,
26 IndShockConsumerType,
27 PerfForesightConsumerType,
28 init_perfect_foresight,
29 init_idiosyncratic_shocks,
30)
31from HARK.ConsumptionSaving.LegacyOOsolvers import (
32 ConsIndShockSolverBasic,
33 ConsPerfForesightSolver,
34)
35from HARK.interpolation import (
36 CubicInterp,
37 LinearInterp,
38 LowerEnvelope,
39 MargValueFuncCRRA,
40 MargMargValueFuncCRRA,
41 ValueFuncCRRA,
42)
43from HARK.metric import MetricObject
44from HARK.numba_tools import (
45 CRRAutility,
46 CRRAutility_inv,
47 CRRAutility_invP,
48 CRRAutilityP,
49 CRRAutilityP_inv,
50 CRRAutilityP_invP,
51 CRRAutilityPP,
52 cubic_interp_fast,
53 linear_interp_deriv_fast,
54 linear_interp_fast,
55)
56from HARK.utilities import NullFunc
58__all__ = [
59 "PerfForesightSolution",
60 "IndShockSolution",
61 "PerfForesightConsumerTypeFast",
62 "IndShockConsumerTypeFast",
63]
65utility = CRRAutility
66utilityP = CRRAutilityP
67utilityPP = CRRAutilityPP
68utilityP_inv = CRRAutilityP_inv
69utility_invP = CRRAutility_invP
70utility_inv = CRRAutility_inv
71utilityP_invP = CRRAutilityP_invP
74# =====================================================================
75# === Classes that help solve consumption-saving models ===
76# =====================================================================
79class PerfForesightSolution(MetricObject):
80 r"""
81 A class representing the solution of a single period of a consumption-saving
82 perfect foresight problem.
84 Here and elsewhere in the code, Nrm indicates that variables are normalized
85 by permanent income.
87 Parameters
88 ----------
89 mNrm: np.array
90 (Normalized) corresponding market resource points for interpolation.
91 cNrm : np.array
92 (Normalized) consumption points for interpolation.
93 vFuncNvrsSlope: float
94 Constant slope of inverse value vFuncNvrs
95 mNrmMin : float
96 The minimum allowable market resources for this period; the consump-
97 tion function (etc) are undefined for m < mNrmMin.
98 hNrm : float
99 Human wealth after receiving income this period: PDV of all future
100 income, ignoring mortality.
101 MPCmin : float
102 Infimum of the marginal propensity to consume this period.
103 MPC --> MPCmin as m --> infinity.
104 MPCmax : float
105 Supremum of the marginal propensity to consume this period.
106 MPC --> MPCmax as m --> mNrmMin.
107 """
109 distance_criteria = ["cNrm", "mNrm"]
111 def __init__(
112 self,
113 mNrm=np.array([0.0, 1.0]),
114 cNrm=np.array([0.0, 1.0]),
115 vFuncNvrsSlope=0.0,
116 mNrmMin=0.0,
117 hNrm=0.0,
118 MPCmin=1.0,
119 MPCmax=1.0,
120 ):
121 self.mNrm = mNrm
122 self.cNrm = cNrm
123 self.vFuncNvrsSlope = vFuncNvrsSlope
124 self.mNrmMin = mNrmMin
125 self.hNrm = hNrm
126 self.MPCmin = MPCmin
127 self.MPCmax = MPCmax
130class IndShockSolution(MetricObject):
131 """
132 A class representing the solution of a single period of a consumption-saving
133 idiosyncratic shocks to permanent and transitory income problem.
135 Parameters
136 ----------
137 mNrm: np.array
138 (Normalized) corresponding market resource points for interpolation.
139 cNrm : np.array
140 (Normalized) consumption points for interpolation.
141 vFuncNvrsSlope: float
142 Constant slope of inverse value ``vFuncNvrs``
143 mNrmMin : float
144 The minimum allowable market resources for this period; the consump-
145 tion function (etc) are undefined for m < mNrmMin.
146 hNrm : float
147 Human wealth after receiving income this period: PDV of all future
148 income, ignoring mortality.
149 MPCmin : float
150 Infimum of the marginal propensity to consume this period.
151 MPC --> MPCmin as m --> infinity.
152 MPCmax : float
153 Supremum of the marginal propensity to consume this period.
154 MPC --> MPCmax as m --> mNrmMin.
155 """
157 distance_criteria = ["cNrm", "mNrm", "mNrmMin"]
159 def __init__(
160 self,
161 mNrm=np.linspace(0, 1),
162 cNrm=np.linspace(0, 1),
163 cFuncLimitIntercept=None,
164 cFuncLimitSlope=None,
165 mNrmMin=0.0,
166 hNrm=0.0,
167 MPCmin=1.0,
168 MPCmax=1.0,
169 Ex_IncNext=0.0,
170 MPC=None,
171 mNrmGrid=None,
172 vNvrs=None,
173 vNvrsP=None,
174 MPCminNvrs=None,
175 ):
176 self.mNrm = mNrm
177 self.cNrm = cNrm
178 self.cFuncLimitIntercept = cFuncLimitIntercept
179 self.cFuncLimitSlope = cFuncLimitSlope
180 self.mNrmMin = mNrmMin
181 self.hNrm = hNrm
182 self.MPCmin = MPCmin
183 self.MPCmax = MPCmax
184 self.Ex_IncNext = Ex_IncNext
185 self.mNrmGrid = mNrmGrid
186 self.vNvrs = vNvrs
187 self.MPCminNvrs = MPCminNvrs
188 self.MPC = MPC
189 self.vNvrsP = vNvrsP
192# =====================================================================
193# === Classes and functions that solve consumption-saving models ===
194# =====================================================================
197def make_solution_terminal_fast(solution_terminal_class, CRRA):
198 solution_terminal = solution_terminal_class()
199 cFunc_terminal = LinearInterp([0.0, 1.0], [0.0, 1.0])
200 solution_terminal.cFunc = cFunc_terminal # c=m at t=T
201 solution_terminal.vFunc = ValueFuncCRRA(cFunc_terminal, CRRA)
202 solution_terminal.vPfunc = MargValueFuncCRRA(cFunc_terminal, CRRA)
203 solution_terminal.vPPfunc = MargMargValueFuncCRRA(cFunc_terminal, CRRA)
204 solution_terminal.MPC = np.array([1.0, 1.0])
205 solution_terminal.MPCminNvrs = 0.0
206 solution_terminal.vNvrs = utility(np.linspace(0.0, 1.0), CRRA)
207 solution_terminal.vNvrsP = utilityP(np.linspace(0.0, 1.0), CRRA)
208 solution_terminal.mNrmGrid = np.linspace(0.0, 1.0)
209 return solution_terminal
212@njit(cache=True)
213def _find_mNrmStE(m, Rfree, PermGroFac, mNrm, cNrm, Ex_IncNext): # pragma: nocover
214 # Make a linear function of all combinations of c and m that yield mNext = mNow
215 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
217 # Find the steady state level of market resources
218 res = interp(mNrm, cNrm, m) - mZeroChange
219 # A zero of this is SS market resources
220 return res
223# @njit(cache=True) can't cache because of use of globals, perhaps newton_secant?
224@njit
225def _add_mNrmStENumba(
226 Rfree, PermGroFac, mNrm, cNrm, mNrmMin, Ex_IncNext, _find_mNrmStE
227): # pragma: nocover
228 """
229 Finds steady state (normalized) market resources and adds it to the
230 solution. This is the level of market resources such that the expectation
231 of market resources in the next period is unchanged. This value doesn't
232 necessarily exist.
233 """
235 # Minimum market resources plus next income is okay starting guess
236 m_init_guess = mNrmMin + Ex_IncNext
238 mNrmStE = newton_secant(
239 _find_mNrmStE,
240 m_init_guess,
241 args=(Rfree, PermGroFac, mNrm, cNrm, Ex_IncNext),
242 disp=False,
243 )
245 if mNrmStE.converged:
246 return mNrmStE.root
247 else:
248 return None
251@njit(cache=True, parallel=True)
252def _solveConsPerfForesightNumba(
253 DiscFac,
254 LivPrb,
255 CRRA,
256 Rfree,
257 PermGroFac,
258 BoroCnstArt,
259 MaxKinks,
260 mNrmNext,
261 cNrmNext,
262 hNrmNext,
263 MPCminNext,
264): # pragma: nocover
265 """
266 Makes the (linear) consumption function for this period.
267 """
269 DiscFacEff = DiscFac * LivPrb
271 # Calculate human wealth this period
272 hNrmNow = (PermGroFac / Rfree) * (hNrmNext + 1.0)
274 # Calculate the lower bound of the marginal propensity to consume
275 APF = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
276 MPCmin = 1.0 / (1.0 + APF / MPCminNext)
278 # Extract the discrete kink points in next period's consumption function;
279 # don't take the last one, as it only defines the extrapolation and is not a kink.
280 mNrmNext = mNrmNext[:-1]
281 cNrmNext = cNrmNext[:-1]
283 # Calculate the end-of-period asset values that would reach those kink points
284 # next period, then invert the first order condition to get consumption. Then
285 # find the endogenous gridpoint (kink point) today that corresponds to each kink
286 aNrmNow = (PermGroFac / Rfree) * (mNrmNext - 1.0)
287 cNrmNow = (DiscFacEff * Rfree) ** (-1.0 / CRRA) * (PermGroFac * cNrmNext)
288 mNrmNow = aNrmNow + cNrmNow
290 # Add an additional point to the list of gridpoints for the extrapolation,
291 # using the new value of the lower bound of the MPC.
292 mNrmNow = np.append(mNrmNow, mNrmNow[-1] + 1.0)
293 cNrmNow = np.append(cNrmNow, cNrmNow[-1] + MPCmin)
295 # If the artificial borrowing constraint binds, combine the constrained and
296 # unconstrained consumption functions.
297 if BoroCnstArt > mNrmNow[0]:
298 # Find the highest index where constraint binds
299 cNrmCnst = mNrmNow - BoroCnstArt
300 CnstBinds = cNrmCnst < cNrmNow
301 idx = np.where(CnstBinds)[0][-1]
303 if idx < (mNrmNow.size - 1):
304 # If it is not the *very last* index, find the the critical level
305 # of mNrm where the artificial borrowing contraint begins to bind.
306 d0 = cNrmNow[idx] - cNrmCnst[idx]
307 d1 = cNrmCnst[idx + 1] - cNrmNow[idx + 1]
308 m0 = mNrmNow[idx]
309 m1 = mNrmNow[idx + 1]
310 alpha = d0 / (d0 + d1)
311 mCrit = m0 + alpha * (m1 - m0)
313 # Adjust the grids of mNrm and cNrm to account for the borrowing constraint.
314 cCrit = mCrit - BoroCnstArt
315 mNrmNow = np.concatenate(
316 (np.array([BoroCnstArt, mCrit]), mNrmNow[(idx + 1) :])
317 )
318 cNrmNow = np.concatenate((np.array([0.0, cCrit]), cNrmNow[(idx + 1) :]))
320 else:
321 # If it *is* the very last index, then there are only three points
322 # that characterize the consumption function: the artificial borrowing
323 # constraint, the constraint kink, and the extrapolation point.
324 mXtra = (cNrmNow[-1] - cNrmCnst[-1]) / (1.0 - MPCmin)
325 mCrit = mNrmNow[-1] + mXtra
326 cCrit = mCrit - BoroCnstArt
327 mNrmNow = np.array([BoroCnstArt, mCrit, mCrit + 1.0])
328 cNrmNow = np.array([0.0, cCrit, cCrit + MPCmin])
330 # If the mNrm and cNrm grids have become too large, throw out the last
331 # kink point, being sure to adjust the extrapolation.
332 if mNrmNow.size > MaxKinks:
333 mNrmNow = np.concatenate((mNrmNow[:-2], np.array([mNrmNow[-3] + 1.0])))
334 cNrmNow = np.concatenate((cNrmNow[:-2], np.array([cNrmNow[-3] + MPCmin])))
336 # Calculate the upper bound of the MPC as the slope of the bottom segment.
337 MPCmax = (cNrmNow[1] - cNrmNow[0]) / (mNrmNow[1] - mNrmNow[0])
339 # Add attributes to enable calculation of steady state market resources.
340 # Relabeling for compatibility with add_mNrmStE
341 mNrmMinNow = mNrmNow[0]
343 # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations
344 vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
346 return (
347 mNrmNow,
348 cNrmNow,
349 vFuncNvrsSlope,
350 mNrmMinNow,
351 hNrmNow,
352 MPCmin,
353 MPCmax,
354 )
357class ConsPerfForesightSolverFast(ConsPerfForesightSolver):
358 """
359 A class for solving a one period perfect foresight consumption-saving problem.
360 An instance of this class is created by the function solvePerfForesight in each period.
361 """
363 def solve(self):
364 """
365 Solves the one period perfect foresight consumption-saving problem.
367 Parameters
368 ----------
369 None
371 Returns
372 -------
373 solution : PerfForesightSolution
374 The solution to this period's problem.
375 """
377 # Use a local value of BoroCnstArt to prevent comparing None and float below.
378 if self.BoroCnstArt is None:
379 BoroCnstArt = -np.inf
380 else:
381 BoroCnstArt = self.BoroCnstArt
383 (
384 self.mNrmNow,
385 self.cNrmNow,
386 self.vFuncNvrsSlope,
387 self.mNrmMinNow,
388 self.hNrmNow,
389 self.MPCmin,
390 self.MPCmax,
391 ) = _solveConsPerfForesightNumba(
392 self.DiscFac,
393 self.LivPrb,
394 self.CRRA,
395 self.Rfree,
396 self.PermGroFac,
397 BoroCnstArt,
398 self.MaxKinks,
399 self.solution_next.mNrm,
400 self.solution_next.cNrm,
401 self.solution_next.hNrm,
402 self.solution_next.MPCmin,
403 )
405 solution = PerfForesightSolution(
406 self.mNrmNow,
407 self.cNrmNow,
408 self.vFuncNvrsSlope,
409 self.mNrmMinNow,
410 self.hNrmNow,
411 self.MPCmin,
412 self.MPCmax,
413 )
414 return solution
417@njit(cache=True)
418def _np_tile(A, reps): # pragma: nocover
419 return A.repeat(reps[0]).reshape(A.size, -1).transpose()
422@njit(cache=True)
423def _np_insert(arr, obj, values, axis=-1): # pragma: nocover
424 return np.append(np.array(values), arr)
427@njit(cache=True, parallel=True)
428def _prepare_to_solveConsIndShockNumba(
429 DiscFac,
430 LivPrb,
431 CRRA,
432 Rfree,
433 PermGroFac,
434 BoroCnstArt,
435 aXtraGrid,
436 hNrmNext,
437 mNrmMinNext,
438 MPCminNext,
439 MPCmaxNext,
440 PermShkValsNext,
441 TranShkValsNext,
442 ShkPrbsNext,
443): # pragma: nocover
444 """
445 Unpacks some of the inputs (and calculates simple objects based on them),
446 storing the results in self for use by other methods. These include:
447 income shocks and probabilities, next period's marginal value function
448 (etc), the probability of getting the worst income shock next period,
449 the patience factor, human wealth, and the bounding MPCs.
451 Defines the constrained portion of the consumption function as cFuncNowCnst,
452 an attribute of self. Uses the artificial and natural borrowing constraints.
454 """
456 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
457 PermShkMinNext = np.min(PermShkValsNext)
458 TranShkMinNext = np.min(TranShkValsNext)
459 WorstIncPrb = np.sum(
460 ShkPrbsNext[
461 (PermShkValsNext * TranShkValsNext) == (PermShkMinNext * TranShkMinNext)
462 ]
463 )
465 # Update the bounding MPCs and PDV of human wealth:
466 APF = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
467 MPCminNow = 1.0 / (1.0 + APF / MPCminNext)
468 Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
469 hNrmNow = PermGroFac / Rfree * (Ex_IncNext + hNrmNext)
470 MPCmaxNow = 1.0 / (1.0 + (WorstIncPrb ** (1.0 / CRRA)) * APF / MPCmaxNext)
472 cFuncLimitIntercept = MPCminNow * hNrmNow
473 cFuncLimitSlope = MPCminNow
475 # Calculate the minimum allowable value of money resources in this period
476 BoroCnstNat = (mNrmMinNext - TranShkMinNext) * (PermGroFac * PermShkMinNext) / Rfree
478 # Note: need to be sure to handle BoroCnstArt==None appropriately.
479 # In Py2, this would evaluate to 5.0: np.max([None, 5.0]).
480 # However in Py3, this raises a TypeError. Thus here we need to directly
481 # address the situation in which BoroCnstArt == None:
482 if BoroCnstArt is None:
483 mNrmMinNow = BoroCnstNat
484 else:
485 mNrmMinNow = np.max(np.array([BoroCnstNat, BoroCnstArt]))
486 if BoroCnstNat < mNrmMinNow:
487 MPCmaxEff = 1.0 # If actually constrained, MPC near limit is 1
488 else:
489 MPCmaxEff = MPCmaxNow
491 """
492 Prepare to calculate end-of-period marginal value by creating an array
493 of market resources that the agent could have next period, considering
494 the grid of end-of-period assets and the distribution of shocks he might
495 experience next period.
496 """
498 # We define aNrmNow all the way from BoroCnstNat up to max(self.aXtraGrid)
499 # even if BoroCnstNat < BoroCnstArt, so we can construct the consumption
500 # function as the lower envelope of the (by the artificial borrowing con-
501 # straint) uconstrained consumption function, and the artificially con-
502 # strained consumption function.
503 aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat
504 ShkCount = TranShkValsNext.size
505 aNrm_temp = _np_tile(aNrmNow, (ShkCount, 1))
507 # Tile arrays of the income shocks and put them into useful shapes
508 aNrmCount = aNrmNow.shape[0]
509 PermShkVals_temp = (_np_tile(PermShkValsNext, (aNrmCount, 1))).transpose()
510 TranShkVals_temp = (_np_tile(TranShkValsNext, (aNrmCount, 1))).transpose()
511 ShkPrbs_temp = (_np_tile(ShkPrbsNext, (aNrmCount, 1))).transpose()
513 # Get cash on hand next period
514 mNrmNext = Rfree / (PermGroFac * PermShkVals_temp) * aNrm_temp + TranShkVals_temp
515 # CDC 20191205: This should be divided by LivPrb[0] for Blanchard insurance
517 return (
518 DiscFacEff,
519 BoroCnstNat,
520 cFuncLimitIntercept,
521 cFuncLimitSlope,
522 mNrmMinNow,
523 hNrmNow,
524 MPCminNow,
525 MPCmaxNow,
526 MPCmaxEff,
527 Ex_IncNext,
528 mNrmNext,
529 PermShkVals_temp,
530 ShkPrbs_temp,
531 aNrmNow,
532 )
535@njit(cache=True, parallel=True)
536def _solveConsIndShockLinearNumba(
537 mNrmMinNext,
538 mNrmNext,
539 CRRA,
540 mNrmUnc,
541 cNrmUnc,
542 DiscFacEff,
543 Rfree,
544 PermGroFac,
545 PermShkVals_temp,
546 ShkPrbs_temp,
547 aNrmNow,
548 BoroCnstNat,
549 cFuncInterceptNext,
550 cFuncSlopeNext,
551): # pragma: nocover
552 """
553 Calculate end-of-period marginal value of assets at each point in aNrmNow.
554 Does so by taking a weighted sum of next period marginal values across
555 income shocks (in a preconstructed grid self.mNrmNext).
556 """
558 mNrmCnst = np.array([mNrmMinNext, mNrmMinNext + 1])
559 cNrmCnst = np.array([0.0, 1.0])
560 cFuncNextCnst = linear_interp_fast(mNrmNext.flatten(), mNrmCnst, cNrmCnst)
561 cFuncNextUnc = linear_interp_fast(
562 mNrmNext.flatten(), mNrmUnc, cNrmUnc, cFuncInterceptNext, cFuncSlopeNext
563 )
564 cFuncNext = np.minimum(cFuncNextCnst, cFuncNextUnc)
565 vPfuncNext = utilityP(cFuncNext, CRRA).reshape(mNrmNext.shape)
567 EndOfPrdvP = (
568 DiscFacEff
569 * Rfree
570 * PermGroFac ** (-CRRA)
571 * np.sum(PermShkVals_temp ** (-CRRA) * vPfuncNext * ShkPrbs_temp, axis=0)
572 )
574 # Finds interpolation points (c,m) for the consumption function.
576 cNrmNow = utilityP_inv(EndOfPrdvP, CRRA)
577 mNrmNow = cNrmNow + aNrmNow
579 # Limiting consumption is zero as m approaches mNrmMin
580 cNrm = _np_insert(cNrmNow, 0, 0.0, axis=-1)
581 mNrm = _np_insert(mNrmNow, 0, BoroCnstNat, axis=-1)
583 return (cNrm, mNrm, EndOfPrdvP)
586class ConsIndShockSolverBasicFast(ConsIndShockSolverBasic):
587 """
588 This class solves a single period of a standard consumption-saving problem,
589 using linear interpolation and without the ability to calculate the value
590 function. ConsIndShockSolver inherits from this class and adds the ability
591 to perform cubic interpolation and to calculate the value function.
593 Note that this class does not have its own initializing method. It initial-
594 izes the same problem in the same way as ConsIndShockSetup, from which it
595 inherits.
596 """
598 def prepare_to_solve(self):
599 """
600 Perform preparatory work before calculating the unconstrained consumption
601 function.
602 Parameters
603 ----------
604 none
605 Returns
606 -------
607 none
608 """
610 self.ShkPrbsNext = self.IncShkDstn.pmv
611 self.PermShkValsNext = self.IncShkDstn.atoms[0]
612 self.TranShkValsNext = self.IncShkDstn.atoms[1]
614 (
615 self.DiscFacEff,
616 self.BoroCnstNat,
617 self.cFuncLimitIntercept,
618 self.cFuncLimitSlope,
619 self.mNrmMinNow,
620 self.hNrmNow,
621 self.MPCminNow,
622 self.MPCmaxNow,
623 self.MPCmaxEff,
624 self.Ex_IncNext,
625 self.mNrmNext,
626 self.PermShkVals_temp,
627 self.ShkPrbs_temp,
628 self.aNrmNow,
629 ) = _prepare_to_solveConsIndShockNumba(
630 self.DiscFac,
631 self.LivPrb,
632 self.CRRA,
633 self.Rfree,
634 self.PermGroFac,
635 self.BoroCnstArt,
636 self.aXtraGrid,
637 self.solution_next.hNrm,
638 self.solution_next.mNrmMin,
639 self.solution_next.MPCmin,
640 self.solution_next.MPCmax,
641 self.PermShkValsNext,
642 self.TranShkValsNext,
643 self.ShkPrbsNext,
644 )
646 def solve(self):
647 """
648 Solves a one period consumption saving problem with risky income.
649 Parameters
650 ----------
651 None
652 Returns
653 -------
654 solution : ConsumerSolution
655 The solution to the one period problem.
656 """
658 self.cNrm, self.mNrm, self.EndOfPrdvP = _solveConsIndShockLinearNumba(
659 self.solution_next.mNrmMin,
660 self.mNrmNext,
661 self.CRRA,
662 self.solution_next.mNrm,
663 self.solution_next.cNrm,
664 self.DiscFacEff,
665 self.Rfree,
666 self.PermGroFac,
667 self.PermShkVals_temp,
668 self.ShkPrbs_temp,
669 self.aNrmNow,
670 self.BoroCnstNat,
671 self.solution_next.cFuncLimitIntercept,
672 self.solution_next.cFuncLimitSlope,
673 )
675 # Pack up the solution and return it
676 solution = IndShockSolution(
677 self.mNrm,
678 self.cNrm,
679 self.cFuncLimitIntercept,
680 self.cFuncLimitSlope,
681 self.mNrmMinNow,
682 self.hNrmNow,
683 self.MPCminNow,
684 self.MPCmaxEff,
685 self.Ex_IncNext,
686 )
688 return solution
691@njit(cache=True, parallel=True)
692def _solveConsIndShockCubicNumba(
693 mNrmMinNext,
694 mNrmNext,
695 mNrmUnc,
696 cNrmUnc,
697 MPCNext,
698 cFuncInterceptNext,
699 cFuncSlopeNext,
700 CRRA,
701 DiscFacEff,
702 Rfree,
703 PermGroFac,
704 PermShkVals_temp,
705 ShkPrbs_temp,
706 aNrmNow,
707 BoroCnstNat,
708 MPCmaxNow,
709): # pragma: nocover
710 mNrmCnst = np.array([mNrmMinNext, mNrmMinNext + 1])
711 cNrmCnst = np.array([0.0, 1.0])
712 cFuncNextCnst, MPCNextCnst = linear_interp_deriv_fast(
713 mNrmNext.flatten(), mNrmCnst, cNrmCnst
714 )
715 cFuncNextUnc, MPCNextUnc = cubic_interp_fast(
716 mNrmNext.flatten(),
717 mNrmUnc,
718 cNrmUnc,
719 MPCNext,
720 cFuncInterceptNext,
721 cFuncSlopeNext,
722 )
724 cFuncNext = np.where(cFuncNextCnst <= cFuncNextUnc, cFuncNextCnst, cFuncNextUnc)
726 vPfuncNext = utilityP(cFuncNext, CRRA).reshape(mNrmNext.shape)
728 EndOfPrdvP = (
729 DiscFacEff
730 * Rfree
731 * PermGroFac ** (-CRRA)
732 * np.sum(PermShkVals_temp ** (-CRRA) * vPfuncNext * ShkPrbs_temp, axis=0)
733 )
734 # Finds interpolation points (c,m) for the consumption function.
736 cNrmNow = EndOfPrdvP ** (-1.0 / CRRA)
737 mNrmNow = cNrmNow + aNrmNow
739 # Limiting consumption is zero as m approaches mNrmMin
740 cNrm = _np_insert(cNrmNow, 0, 0.0, axis=-1)
741 mNrm = _np_insert(mNrmNow, 0, BoroCnstNat, axis=-1)
743 """
744 Makes a cubic spline interpolation of the unconstrained consumption
745 function for this period.
746 """
748 MPCinterpNext = np.where(cFuncNextCnst <= cFuncNextUnc, MPCNextCnst, MPCNextUnc)
749 vPPfuncNext = (MPCinterpNext * utilityPP(cFuncNext, CRRA)).reshape(mNrmNext.shape)
751 EndOfPrdvPP = (
752 DiscFacEff
753 * Rfree
754 * Rfree
755 * PermGroFac ** (-CRRA - 1.0)
756 * np.sum(PermShkVals_temp ** (-CRRA - 1.0) * vPPfuncNext * ShkPrbs_temp, axis=0)
757 )
758 dcda = EndOfPrdvPP / utilityPP(cNrm[1:], CRRA)
759 MPC = dcda / (dcda + 1.0)
760 MPC = _np_insert(MPC, 0, MPCmaxNow)
762 return cNrm, mNrm, MPC, EndOfPrdvP
765@njit(cache=True)
766def _cFuncCubic(
767 aXtraGrid, mNrmMinNow, mNrmNow, cNrmNow, MPCNow, MPCminNow, hNrmNow
768): # pragma: nocover
769 mNrmGrid = mNrmMinNow + aXtraGrid
770 mNrmCnst = np.array([mNrmMinNow, mNrmMinNow + 1])
771 cNrmCnst = np.array([0.0, 1.0])
772 cFuncNowCnst = linear_interp_fast(mNrmGrid.flatten(), mNrmCnst, cNrmCnst)
773 cFuncNowUnc, MPCNowUnc = cubic_interp_fast(
774 mNrmGrid.flatten(), mNrmNow, cNrmNow, MPCNow, MPCminNow * hNrmNow, MPCminNow
775 )
777 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
779 return cNrmNow, mNrmGrid
782@njit(cache=True)
783def _cFuncLinear(
784 aXtraGrid, mNrmMinNow, mNrmNow, cNrmNow, MPCminNow, hNrmNow
785): # pragma: nocover
786 mNrmGrid = mNrmMinNow + aXtraGrid
787 mNrmCnst = np.array([mNrmMinNow, mNrmMinNow + 1])
788 cNrmCnst = np.array([0.0, 1.0])
789 cFuncNowCnst = linear_interp_fast(mNrmGrid.flatten(), mNrmCnst, cNrmCnst)
790 cFuncNowUnc = linear_interp_fast(
791 mNrmGrid.flatten(), mNrmNow, cNrmNow, MPCminNow * hNrmNow, MPCminNow
792 )
794 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
796 return cNrmNow, mNrmGrid
799@njit(cache=True)
800def _add_vFuncNumba(
801 mNrmNext,
802 mNrmGridNext,
803 vNvrsNext,
804 vNvrsPNext,
805 MPCminNvrsNext,
806 hNrmNext,
807 CRRA,
808 PermShkVals_temp,
809 PermGroFac,
810 DiscFacEff,
811 ShkPrbs_temp,
812 EndOfPrdvP,
813 aNrmNow,
814 BoroCnstNat,
815 mNrmGrid,
816 cFuncNow,
817 mNrmMinNow,
818 MPCmaxEff,
819 MPCminNow,
820): # pragma: nocover
821 """
822 Construct the end-of-period value function for this period, storing it
823 as an attribute of self for use by other methods.
824 """
826 # vFunc always cubic
828 vNvrsFuncNow, _ = cubic_interp_fast(
829 mNrmNext.flatten(),
830 mNrmGridNext,
831 vNvrsNext,
832 vNvrsPNext,
833 MPCminNvrsNext * hNrmNext,
834 MPCminNvrsNext,
835 )
837 vFuncNext = utility(vNvrsFuncNow, CRRA).reshape(mNrmNext.shape)
839 VLvlNext = (
840 PermShkVals_temp ** (1.0 - CRRA) * PermGroFac ** (1.0 - CRRA)
841 ) * vFuncNext
842 EndOfPrdv = DiscFacEff * np.sum(VLvlNext * ShkPrbs_temp, axis=0)
844 # value transformed through inverse utility
845 EndOfPrdvNvrs = utility_inv(EndOfPrdv, CRRA)
846 EndOfPrdvNvrsP = EndOfPrdvP * utility_invP(EndOfPrdv, CRRA)
847 EndOfPrdvNvrs = _np_insert(EndOfPrdvNvrs, 0, 0.0)
849 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
850 EndOfPrdvNvrsP = _np_insert(EndOfPrdvNvrsP, 0, EndOfPrdvNvrsP[0])
851 aNrm_temp = _np_insert(aNrmNow, 0, BoroCnstNat)
853 """
854 Creates the value function for this period, defined over market resources m.
855 self must have the attribute EndOfPrdvFunc in order to execute.
856 """
858 # Compute expected value and marginal value on a grid of market resources
860 aNrmNow = mNrmGrid - cFuncNow
862 EndOfPrdvNvrsFunc, _ = cubic_interp_fast(
863 aNrmNow, aNrm_temp, EndOfPrdvNvrs, EndOfPrdvNvrsP
864 )
865 EndOfPrdvFunc = utility(EndOfPrdvNvrsFunc, CRRA)
867 vNrmNow = utility(cFuncNow, CRRA) + EndOfPrdvFunc
868 vPnow = utilityP(cFuncNow, CRRA)
870 # Construct the beginning-of-period value function
871 vNvrs = utility_inv(vNrmNow, CRRA) # value transformed through inverse utility
872 vNvrsP = vPnow * utility_invP(vNrmNow, CRRA)
873 mNrmGrid = _np_insert(mNrmGrid, 0, mNrmMinNow)
874 vNvrs = _np_insert(vNvrs, 0, 0.0)
875 vNvrsP = _np_insert(vNvrsP, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
876 MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
878 return (
879 mNrmGrid,
880 vNvrs,
881 vNvrsP,
882 MPCminNvrs,
883 )
886@njit
887def _add_mNrmStEIndNumba(
888 PermGroFac,
889 Rfree,
890 Ex_IncNext,
891 mNrmMin,
892 mNrm,
893 cNrm,
894 MPC,
895 MPCmin,
896 hNrm,
897 _searchfunc,
898): # pragma: nocover
899 """
900 Finds steady state (normalized) market resources and adds it to the
901 solution. This is the level of market resources such that the expectation
902 of market resources in the next period is unchanged. This value doesn't
903 necessarily exist.
904 """
906 # Minimum market resources plus next income is okay starting guess
907 m_init_guess = mNrmMin + Ex_IncNext
909 mNrmStE = newton_secant(
910 _searchfunc,
911 m_init_guess,
912 args=(PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm),
913 disp=False,
914 )
916 if mNrmStE.converged:
917 return mNrmStE.root
918 else:
919 return None
922@njit(cache=True)
923def _find_mNrmStELinear(
924 m, PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm
925): # pragma: nocover
926 # Make a linear function of all combinations of c and m that yield mNext = mNow
927 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
929 mNrmCnst = np.array([mNrmMin, mNrmMin + 1])
930 cNrmCnst = np.array([0.0, 1.0])
931 cFuncNowCnst = linear_interp_fast(np.array([m]), mNrmCnst, cNrmCnst)
932 cFuncNowUnc = linear_interp_fast(np.array([m]), mNrm, cNrm, MPCmin * hNrm, MPCmin)
934 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
936 # Find the steady state level of market resources
937 res = cNrmNow[0] - mZeroChange
938 # A zero of this is SS market resources
939 return res
942@njit(cache=True)
943def _find_mNrmStECubic(
944 m, PermGroFac, Rfree, Ex_IncNext, mNrmMin, mNrm, cNrm, MPC, MPCmin, hNrm
945): # pragma: nocover
946 # Make a linear function of all combinations of c and m that yield mNext = mNow
947 mZeroChange = (1.0 - PermGroFac / Rfree) * m + (PermGroFac / Rfree) * Ex_IncNext
949 mNrmCnst = np.array([mNrmMin, mNrmMin + 1])
950 cNrmCnst = np.array([0.0, 1.0])
951 cFuncNowCnst = linear_interp_fast(np.array([m]), mNrmCnst, cNrmCnst)
952 cFuncNowUnc, MPCNowUnc = cubic_interp_fast(
953 np.array([m]), mNrm, cNrm, MPC, MPCmin * hNrm, MPCmin
954 )
956 cNrmNow = np.where(cFuncNowCnst <= cFuncNowUnc, cFuncNowCnst, cFuncNowUnc)
958 # Find the steady state level of market resources
959 res = cNrmNow[0] - mZeroChange
960 # A zero of this is SS market resources
961 return res
964class ConsIndShockSolverFast(ConsIndShockSolverBasicFast):
965 r"""
966 This class solves a single period of a standard consumption-saving problem.
967 It inherits from ConsIndShockSolverBasic, adding the ability to perform cubic
968 interpolation and to calculate the value function.
969 """
971 def solve(self):
972 """
973 Solves a one period consumption saving problem with risky income.
974 Parameters
975 ----------
976 None
977 Returns
978 -------
979 solution : ConsumerSolution
980 The solution to the one period problem.
981 """
983 if self.CubicBool:
984 (
985 self.cNrm,
986 self.mNrm,
987 self.MPC,
988 self.EndOfPrdvP,
989 ) = _solveConsIndShockCubicNumba(
990 self.solution_next.mNrmMin,
991 self.mNrmNext,
992 self.solution_next.mNrm,
993 self.solution_next.cNrm,
994 self.solution_next.MPC,
995 self.solution_next.cFuncLimitIntercept,
996 self.solution_next.cFuncLimitSlope,
997 self.CRRA,
998 self.DiscFacEff,
999 self.Rfree,
1000 self.PermGroFac,
1001 self.PermShkVals_temp,
1002 self.ShkPrbs_temp,
1003 self.aNrmNow,
1004 self.BoroCnstNat,
1005 self.MPCmaxNow,
1006 )
1007 # Pack up the solution and return it
1008 solution = IndShockSolution(
1009 self.mNrm,
1010 self.cNrm,
1011 self.cFuncLimitIntercept,
1012 self.cFuncLimitSlope,
1013 self.mNrmMinNow,
1014 self.hNrmNow,
1015 self.MPCminNow,
1016 self.MPCmaxEff,
1017 self.Ex_IncNext,
1018 self.MPC,
1019 )
1020 else:
1021 self.cNrm, self.mNrm, self.EndOfPrdvP = _solveConsIndShockLinearNumba(
1022 self.solution_next.mNrmMin,
1023 self.mNrmNext,
1024 self.CRRA,
1025 self.solution_next.mNrm,
1026 self.solution_next.cNrm,
1027 self.DiscFacEff,
1028 self.Rfree,
1029 self.PermGroFac,
1030 self.PermShkVals_temp,
1031 self.ShkPrbs_temp,
1032 self.aNrmNow,
1033 self.BoroCnstNat,
1034 self.solution_next.cFuncLimitIntercept,
1035 self.solution_next.cFuncLimitSlope,
1036 )
1038 # Pack up the solution and return it
1039 solution = IndShockSolution(
1040 self.mNrm,
1041 self.cNrm,
1042 self.cFuncLimitIntercept,
1043 self.cFuncLimitSlope,
1044 self.mNrmMinNow,
1045 self.hNrmNow,
1046 self.MPCminNow,
1047 self.MPCmaxEff,
1048 self.Ex_IncNext,
1049 )
1051 if self.vFuncBool:
1052 if self.CubicBool:
1053 self.cFuncNow, self.mNrmGrid = _cFuncCubic(
1054 self.aXtraGrid,
1055 self.mNrmMinNow,
1056 self.mNrm,
1057 self.cNrm,
1058 self.MPC,
1059 self.MPCminNow,
1060 self.hNrmNow,
1061 )
1062 else:
1063 self.cFuncNow, self.mNrmGrid = _cFuncLinear(
1064 self.aXtraGrid,
1065 self.mNrmMinNow,
1066 self.mNrm,
1067 self.cNrm,
1068 self.MPCminNow,
1069 self.hNrmNow,
1070 )
1072 self.mNrmGrid, self.vNvrs, self.vNvrsP, self.MPCminNvrs = _add_vFuncNumba(
1073 self.mNrmNext,
1074 self.solution_next.mNrmGrid,
1075 self.solution_next.vNvrs,
1076 self.solution_next.vNvrsP,
1077 self.solution_next.MPCminNvrs,
1078 self.solution_next.hNrm,
1079 self.CRRA,
1080 self.PermShkVals_temp,
1081 self.PermGroFac,
1082 self.DiscFacEff,
1083 self.ShkPrbs_temp,
1084 self.EndOfPrdvP,
1085 self.aNrmNow,
1086 self.BoroCnstNat,
1087 self.mNrmGrid,
1088 self.cFuncNow,
1089 self.mNrmMinNow,
1090 self.MPCmaxEff,
1091 self.MPCminNow,
1092 )
1094 # Pack up the solution and return it
1096 solution.mNrmGrid = self.mNrmGrid
1097 solution.vNvrs = self.vNvrs
1098 solution.vNvrsP = self.vNvrsP
1099 solution.MPCminNvrs = self.MPCminNvrs
1101 return solution
1104# ============================================================================
1105# == Classes for representing types of consumer agents (and things they do) ==
1106# ============================================================================
1109init_perfect_foresight_fast = init_perfect_foresight.copy()
1110perf_foresight_constructor_dict = init_perfect_foresight["constructors"].copy()
1111perf_foresight_constructor_dict["solution_terminal"] = make_solution_terminal_fast
1112init_perfect_foresight_fast["constructors"] = perf_foresight_constructor_dict
1115class PerfForesightConsumerTypeFast(PerfForesightConsumerType):
1116 r"""
1117 A version of the perfect foresight consumer type speed up by numba.
1118 """
1120 solution_terminal_class = PerfForesightSolution
1121 default_ = {
1122 "params": init_perfect_foresight_fast,
1123 "solver": make_one_period_oo_solver(ConsPerfForesightSolverFast),
1124 "model": "ConsPerfForesight.yaml",
1125 }
1127 def post_solve(self):
1128 self.solution_fast = deepcopy(self.solution)
1130 if self.cycles == 0:
1131 terminal = 1
1132 else:
1133 terminal = self.cycles
1134 self.solution[terminal] = self.solution_terminal
1136 for i in range(terminal):
1137 solution = self.solution[i]
1139 # Construct the consumption function as a linear interpolation.
1140 cFunc = LinearInterp(solution.mNrm, solution.cNrm)
1142 """
1143 Defines the value and marginal value functions for this period.
1144 Uses the fact that for a perfect foresight CRRA utility problem,
1145 if the MPC in period t is :math:`\\kappa_{t}`, and relative risk
1146 aversion :math:`\rho`, then the inverse value vFuncNvrs has a
1147 constant slope of :math:`\\kappa_{t}^{-\rho/(1-\rho)}` and
1148 vFuncNvrs has value of zero at the lower bound of market resources
1149 mNrmMin. See PerfForesightConsumerType.ipynb documentation notebook
1150 for a brief explanation and the links below for a fuller treatment.
1152 https://www.econ2.jhu.edu/people/ccarroll/public/lecturenotes/consumption/PerfForesightCRRA/#vFuncAnalytical
1153 https://www.econ2.jhu.edu/people/ccarroll/SolvingMicroDSOPs/#vFuncPF
1154 """
1156 vFuncNvrs = LinearInterp(
1157 np.array([solution.mNrmMin, solution.mNrmMin + 1.0]),
1158 np.array([0.0, solution.vFuncNvrsSlope]),
1159 )
1160 vFunc = ValueFuncCRRA(vFuncNvrs, self.CRRA)
1161 vPfunc = MargValueFuncCRRA(cFunc, self.CRRA)
1163 consumer_solution = ConsumerSolution(
1164 cFunc=cFunc,
1165 vFunc=vFunc,
1166 vPfunc=vPfunc,
1167 mNrmMin=solution.mNrmMin,
1168 hNrm=solution.hNrm,
1169 MPCmin=solution.MPCmin,
1170 MPCmax=solution.MPCmax,
1171 )
1173 Ex_IncNext = 1.0 # Perfect foresight income of 1
1175 # Add mNrmStE to the solution and return it
1176 consumer_solution.mNrmStE = _add_mNrmStENumba(
1177 self.Rfree[i],
1178 self.PermGroFac[i],
1179 solution.mNrm,
1180 solution.cNrm,
1181 solution.mNrmMin,
1182 Ex_IncNext,
1183 _find_mNrmStE,
1184 )
1186 self.solution[i] = consumer_solution
1189###############################################################################
1192def select_fast_solver(CubicBool, vFuncBool):
1193 if (not CubicBool) and (not vFuncBool):
1194 solver = ConsIndShockSolverBasicFast
1195 else: # Use the "advanced" solver if either is requested
1196 solver = ConsIndShockSolverFast
1197 solve_one_period = make_one_period_oo_solver(solver)
1198 return solve_one_period
1201init_idiosyncratic_shocks_fast = init_idiosyncratic_shocks.copy()
1202ind_shock_fast_constructor_dict = init_idiosyncratic_shocks["constructors"].copy()
1203ind_shock_fast_constructor_dict["solution_terminal"] = make_solution_terminal_fast
1204ind_shock_fast_constructor_dict["solve_one_period"] = select_fast_solver
1205init_idiosyncratic_shocks_fast["constructors"] = ind_shock_fast_constructor_dict
1208class IndShockConsumerTypeFast(IndShockConsumerType, PerfForesightConsumerTypeFast):
1209 r"""
1210 A version of the idiosyncratic shock consumer type sped up by numba.
1212 If CubicBool and vFuncBool are both set to false it's further optimized.
1213 """
1215 solution_terminal_class = IndShockSolution
1216 default_ = {
1217 "params": init_idiosyncratic_shocks_fast,
1218 "solver": NullFunc(),
1219 "model": "ConsIndShock.yaml",
1220 }
1222 def post_solve(self):
1223 self.solution_fast = deepcopy(self.solution)
1225 if self.cycles == 0:
1226 cycles = 1
1227 else:
1228 cycles = self.cycles
1229 self.solution[-1] = init_idiosyncratic_shocks["constructors"][
1230 "solution_terminal"
1231 ](self.CRRA)
1233 for i in range(cycles):
1234 for j in range(self.T_cycle):
1235 solution = self.solution[i * self.T_cycle + j]
1237 # Define the borrowing constraint (limiting consumption function)
1238 cFuncNowCnst = LinearInterp(
1239 np.array([solution.mNrmMin, solution.mNrmMin + 1]),
1240 np.array([0.0, 1.0]),
1241 )
1243 """
1244 Constructs a basic solution for this period, including the consumption
1245 function and marginal value function.
1246 """
1248 if self.CubicBool:
1249 # Makes a cubic spline interpolation of the unconstrained consumption
1250 # function for this period.
1251 cFuncNowUnc = CubicInterp(
1252 solution.mNrm,
1253 solution.cNrm,
1254 solution.MPC,
1255 solution.cFuncLimitIntercept,
1256 solution.cFuncLimitSlope,
1257 )
1258 else:
1259 # Makes a linear interpolation to represent the (unconstrained) consumption function.
1260 # Construct the unconstrained consumption function
1261 cFuncNowUnc = LinearInterp(
1262 solution.mNrm,
1263 solution.cNrm,
1264 solution.cFuncLimitIntercept,
1265 solution.cFuncLimitSlope,
1266 )
1268 # Combine the constrained and unconstrained functions into the true consumption function
1269 cFuncNow = LowerEnvelope(cFuncNowUnc, cFuncNowCnst)
1271 # Make the marginal value function and the marginal marginal value function
1272 vPfuncNow = MargValueFuncCRRA(cFuncNow, self.CRRA)
1274 # Pack up the solution and return it
1275 consumer_solution = ConsumerSolution(
1276 cFunc=cFuncNow,
1277 vPfunc=vPfuncNow,
1278 mNrmMin=solution.mNrmMin,
1279 hNrm=solution.hNrm,
1280 MPCmin=solution.MPCmin,
1281 MPCmax=solution.MPCmax,
1282 )
1284 if self.vFuncBool:
1285 vNvrsFuncNow = CubicInterp(
1286 solution.mNrmGrid,
1287 solution.vNvrs,
1288 solution.vNvrsP,
1289 solution.MPCminNvrs * solution.hNrm,
1290 solution.MPCminNvrs,
1291 )
1292 vFuncNow = ValueFuncCRRA(vNvrsFuncNow, self.CRRA)
1294 consumer_solution.vFunc = vFuncNow
1296 if self.CubicBool or self.vFuncBool:
1297 _searchFunc = (
1298 _find_mNrmStECubic if self.CubicBool else _find_mNrmStELinear
1299 )
1300 # Add mNrmStE to the solution and return it
1301 consumer_solution.mNrmStE = _add_mNrmStEIndNumba(
1302 self.PermGroFac[j],
1303 self.Rfree[j],
1304 solution.Ex_IncNext,
1305 solution.mNrmMin,
1306 solution.mNrm,
1307 solution.cNrm,
1308 solution.MPC,
1309 solution.MPCmin,
1310 solution.hNrm,
1311 _searchFunc,
1312 )
1314 self.solution[i * self.T_cycle + j] = consumer_solution
1316 if (self.cycles == 0) and (self.T_cycle == 1):
1317 self.calc_stable_points(force=True)