Coverage for HARK / ConsumptionSaving / ConsHabitModel.py: 99%
227 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 06:00 +0000
1"""
2This module has consumption-models with habit formation. Right now, it only has
3a single basic model with permanent and transitory income shocks, one risk-free
4asset, and a habit stock that evolves as a weighted average of current consumption
5and the prior habit stock.
6"""
8import numpy as np
9from HARK.utilities import make_exponential_grid
10from HARK.interpolation import (
11 LinearInterp,
12 LinearInterpOnInterp1D,
13 ConstantFunction,
14 Curvilinear2DInterp,
15 LowerEnvelope2D,
16 IdentityFunction,
17 BilinearInterp,
18 MargValueFuncCRRA,
19 ValueFuncCRRA,
20)
21from HARK.distributions import expected, Lognormal
22from HARK.core import AgentType, get_it_from
23from HARK.Calibration.Income.IncomeProcesses import (
24 construct_lognormal_income_process_unemployment,
25 get_PermShkDstn_from_IncShkDstn,
26 get_TranShkDstn_from_IncShkDstn,
27)
28from HARK.utilities import make_assets_grid
29from HARK.rewards import UtilityFuncCRRA
30from HARK.ConsumptionSaving.ConsIndShockModel import (
31 make_lognormal_kNrm_init_dstn,
32 make_lognormal_pLvl_init_dstn,
33)
34from HARK.Calibration.Assets.AssetProcesses import (
35 make_lognormal_RiskyDstn,
36 calc_ShareLimit_for_CRRA,
37)
38from HARK.ConsumptionSaving.ConsRiskyAssetModel import make_simple_ShareGrid
41def make_lognormal_habit_init_dstn(hLogInitMean, hLogInitStd, HabitInitCount, RNG):
42 """
43 Construct a lognormal distribution for (normalized) initial habit stock
44 of newborns, hNrm. This is the default constructor for HabitInitDstn.
46 Parameters
47 ----------
48 hLogInitMean : float
49 Mean of log habit stock for newborns.
50 hLogInitStd : float
51 Stdev of log habit stock for newborns.
52 HabitInitCount : int
53 Number of points in the discretization.
54 RNG : np.random.RandomState
55 Agent's internal RNG.
57 Returns
58 -------
59 HabitInitDstn : DiscreteDistribution
60 Discretized distribution of initial habit stock for newborns.
61 """
62 dstn = Lognormal(
63 mu=hLogInitMean,
64 sigma=hLogInitStd,
65 seed=RNG.integers(0, 2**31 - 1),
66 )
67 HabitInitDstn = dstn.discretize(HabitInitCount)
68 return HabitInitDstn
71class HabitFormationInverter:
72 """
73 A class for solving the first order conditions of a consumption-saving model
74 with habit formation. In this notation, HabitRte is a parameter on the unit
75 interval representing how fast the habit stock evolves; a value of zero means
76 no habit dynamics and a value of one means that H_t = c_t, complete updating.
77 HabitWgt is also on the unit interval and represents the exponent on the
78 habit stock, which is used as a divisor on consumption in the utility function.
80 Instances of this class take two arguments when called as a function: end-of-
81 period habit stock H and transformed end-of-period marginal value chi.
83 chi = (W_a(a,H) - lambda * W_H(a,H)) ** (-1/rho)
85 a = m - c
86 H = lambda * c + (1-lambda) * h
87 m' = R a / psi + theta
88 h' = H / psi
90 Parameters
91 ----------
92 CRRA : float
93 Coefficient of relative risk aversion, rho.
94 HabitRte : float
95 Rate of habit stock updating with new consumption, lambda. Must be greater
96 than zero but less than one.
97 HabitWgt : float
98 Weight of habit stock in preferences, alpha; exponent on habits as a divisor
99 in utility function. Must be greater than zero but less than one.
100 ChiMax : float
101 Largest value of "transformed marginal value" to consider in the grid.
102 The minimum value is always zero. These chi values are "consumption-like".
103 ChiCount : int
104 Number of gridpoints in the "transformed marginal value" grid.
105 ChiOrder : float
106 Strictly positive exponential order for the "transformed marginal value" grid.
107 HabitMax : float
108 Largest value in the habit grid to consider; minimum is always zero.
109 HabitCount : int
110 Number of gridpoints in the habit stock grid.
111 HabitOrder : float
112 Strictly positive exponential order for the habit stock grid.
113 """
115 def __init__(
116 self,
117 CRRA,
118 HabitRte,
119 HabitWgt,
120 ChiMax,
121 ChiCount,
122 ChiOrder,
123 HabitMax,
124 HabitCount,
125 HabitOrder,
126 ):
127 # Validation
128 if HabitRte > 1.0:
129 raise ValueError("HabitRte must be no greater than 1!")
130 if HabitRte <= 0.0:
131 raise ValueError("HabitRte must be strictly positive!")
132 if HabitWgt > 1.0:
133 raise ValueError("HabitWgt must be no greater than 1!")
134 if HabitWgt <= 0.0:
135 raise ValueError("HabitWgt must be strictly positive!")
137 xGrid = make_exponential_grid(0.0, ChiMax, ChiCount, ChiOrder)
138 hGrid = make_exponential_grid(0.0, HabitMax, HabitCount, HabitOrder)
139 hMesh, xMesh = np.meshgrid(hGrid, xGrid, indexing="ij")
141 PostHabit = (
142 HabitRte * hMesh ** (HabitWgt * (1.0 - 1.0 / CRRA)) * xMesh
143 + (1.0 - HabitRte) * hMesh
144 )
146 func_by_chi = [LinearInterp(PostHabit[:, j], hGrid) for j in range(ChiCount)]
147 func = LinearInterpOnInterp1D(func_by_chi, xGrid)
149 self.hFunc = func
150 self.rate = HabitRte
152 def __call__(self, H, chi):
153 h = self.hFunc(H, chi)
154 rate = self.rate
155 c = (H - (1 - rate) * h) / rate
156 return c, h
158 def cFunc(self, H, chi):
159 return self(H, chi)[0] # just return consumption
162def make_habit_inverter(
163 CRRA,
164 HabitRte,
165 HabitWgt,
166 ChiMax,
167 ChiCount,
168 ChiOrder,
169 HabitMax,
170 HabitCount,
171 HabitOrder,
172):
173 return HabitFormationInverter(
174 CRRA,
175 HabitRte,
176 HabitWgt,
177 ChiMax,
178 ChiCount,
179 ChiOrder,
180 HabitMax,
181 HabitCount,
182 HabitOrder,
183 )
186def make_habit_grid(HabitMin, HabitMax, HabitCount, HabitOrder):
187 return make_exponential_grid(HabitMin, HabitMax, HabitCount, HabitOrder)
190def make_dense_grids(
191 HabitMin,
192 HabitMax,
193 HabitCount,
194 HabitOrder,
195 aXtraMin,
196 aXtraMax,
197 aXtraCount,
198 aXtraNestFac,
199 aXtraExtra,
200 DenseFactor,
201):
202 HabitGridDense = make_exponential_grid(
203 HabitMin, HabitMax, HabitCount * DenseFactor, HabitOrder
204 )
205 mXtraGridDense = make_assets_grid(
206 aXtraMin, aXtraMax, aXtraCount * DenseFactor, aXtraExtra, aXtraNestFac
207 )
208 out = {"hGridDense": HabitGridDense, "mXtraGrid": mXtraGridDense}
209 return out
212def make_habit_solution_terminal():
213 """
214 Make a pseudo-terminal solution for the habit formation model, which has zero
215 functions for (marginal) value.
216 """
217 dvdkFunc_terminal = ConstantFunction(0.0)
218 dvdhFunc_terminal = ConstantFunction(0.0)
219 solution_terminal = {
220 "dvdkFunc": dvdkFunc_terminal,
221 "dvdhFunc": dvdhFunc_terminal,
222 "kNrmMin": 0.0,
223 }
224 return solution_terminal
227def calc_marg_values(S, k, hpre, rho, R, Gamma, alpha, lamda, beta, C, Vp):
228 """
229 Helper function for computing expected marginal value with respect to market
230 resources and habit stock. Used internally by solve_one_period_ConsHabit.
232 The code here uses "math notation" for quick programming. See the only place
233 in the code where this function is used for a translation of the symbols.
234 """
235 psi = S["PermShk"]
236 theta = S["TranShk"]
237 G = psi * Gamma
238 m = R * k / G + theta
239 h = hpre / G
240 c = C(m, h)
241 a = m - c
242 H = lamda * c + (1 - lamda) * h
243 dvdH = beta * Vp(a, H)
244 temp = h ** (-alpha * (1 - rho))
245 dudc = temp * c ** (-rho)
246 dudh = c ** (1 - rho) * (-alpha) * temp / h
247 dvdm = dudc + lamda * dvdH
248 dvdh = dudh + (1 - lamda) * dvdH
249 G_adj = G ** ((1 - rho) * (1 - alpha) - 1.0)
250 dvdk = R * G_adj * dvdm
251 dvdh = G_adj * dvdh
252 return dvdk, dvdh
255def calc_mid_dvds(shocks, w_nrm, H_nrm, share, rfree, dvdkFunc_next):
256 """
257 Compute middle-of-period marginal value of risky asset share by taking expec-
258 tations over risky return shocks. Uses the next period's marginal value func-
259 tions directly (they already integrate over income shocks).
261 Parameters
262 ----------
263 shocks : float or np.array
264 Risky return realizations.
265 w_nrm : np.array
266 Pre-return savings (w = m - c).
267 H_nrm : np.array
268 End-of-period habit stock.
269 share : np.array
270 Risky share.
271 rfree : float
272 Risk-free return factor.
273 dvdkFunc_next : callable
274 Next period's marginal value of capital, dvdk(k, hPre).
276 Returns
277 -------
278 dvds : np.array
279 Marginal value of risky share (zero at optimum).
280 """
281 ex_ret = shocks - rfree
282 rport = rfree + share * ex_ret
283 a_nrm = rport * w_nrm # post-return assets = next period's kNrm
285 dvdk = dvdkFunc_next(a_nrm, H_nrm)
286 dvds = ex_ret * w_nrm * dvdk
287 return dvds
290def calc_mid_dvdx(shocks, w_nrm, H_nrm, share, rfree, dvdkFunc_next, dvdhFunc_next):
291 """
292 Compute middle-of-period marginal value of wealth and habit stock by taking
293 expectations over risky return shocks. Uses the next period's marginal value
294 function directly (they already integrate over income shocks).
296 Parameters
297 ----------
298 shocks : float or np.array
299 Risky return realizations.
300 w_nrm : np.array
301 Pre-return savings (w = m - c).
302 H_nrm : np.array
303 End-of-period habit stock.
304 share : np.array
305 Risky share.
306 rfree : float
307 Risk-free return factor.
308 dvdkFunc_next : callable
309 Next period's marginal value of capital, dvdk(k, hPre).
310 dvdhFunc_next : callable
311 Next period's marginal value of habit stock, dvdh(k, hPre).
313 Returns
314 -------
315 dvdw : np.array
316 Marginal value of pre-return savings.
317 dvdH : np.array
318 Marginal value of end-of-period habit stock.
319 """
320 ex_ret = shocks - rfree
321 rport = rfree + share * ex_ret
322 a_nrm = rport * w_nrm # post-return assets = next period's kNrm
324 dvdk = dvdkFunc_next(a_nrm, H_nrm)
325 dvdw = rport * dvdk
326 dvdH = dvdhFunc_next(a_nrm, H_nrm)
327 return dvdw, dvdH
330###############################################################################
333def solve_one_period_ConsHabit(
334 solution_next,
335 IncShkDstn,
336 LivPrb,
337 DiscFac,
338 CRRA,
339 Rfree,
340 PermGroFac,
341 BoroCnstArt,
342 aXtraGrid,
343 HabitGrid,
344 FOCinverter,
345 HabitWgt,
346 HabitRte,
347 mXtraGrid,
348 hGridDense,
349):
350 """
351 Solve one period of the consumption-saving model with habit formation.
353 Parameters
354 ----------
355 solution_next : dict
356 Dictionary with next period's solution.
357 IncShkDstn : DiscreteDistribution
358 Discretized permanent and transitory income shock distribution this period.
359 LivPrb : float
360 Survival probability at the end of this period.
361 DiscFac : float
362 Intertemporal discount factor.
363 CRRA : float
364 Coefficient of relative risk aversion.
365 Rfree : float
366 Interest factor on capital at the start of this period.
367 PermGroFac : float
368 Permanent income growth factor at the start of this period.
369 BoroCnstArt : float or None
370 Artificial borrowing constraint on assets at the end of this period,
371 as a fraction of permanent income.
372 aXtraGrid : np.array
373 Grid of "assets above minimum".
374 HabitGrid : np.array
375 Grid of consumption habit stocks on which to solve the problem.
376 FOCinverter : HabitFormationInverter
377 Function that inverts the first order conditions to yield optimal consumption
378 and the decision-time habit stock from which it was chosen.
379 HabitWgt : float
380 Exponent on habit stock, which is used as a divisor on consumption in
381 the utility function: U(c,h) = u(c / h**alpha). Should be on unit interval.
382 HabitRte : float
383 Rate at which habit stock is updated by new consumption: H = lambda*c + (1-lambda)*h.
384 Should be on the unit interval.
385 mXtraGrid : np.array
386 Dense grid of market resources, used to "re-interpolate" the curvilinear
387 consumption function onto a rectilinear grid.
388 hGridDense : np.array
389 Dense grid of habit stocks, used to "re-interpolate" the curvilinear
390 consumption function onto a rectilinear grid.
392 Returns
393 -------
394 solution_now : dict
395 Solution to this period's problem, with the following keys:
396 cFunc : Consumption function over (mNrm, hNrm).
397 dvdkFunc : Marginal value of beginning-of-period capital, defined on (kNrm, hPre).
398 dvdhFunc : Marginal value of beginning-of-period habit stock, defined on (kNrm, hPre).
399 kNrmMin : Minimum allowable beginning-of-period capital.
400 """
401 U = UtilityFuncCRRA(CRRA)
402 DiscFacEff = DiscFac * LivPrb
404 # Make end-of-period state grids
405 aNrmMin = solution_next["kNrmMin"]
406 aGrid = aXtraGrid + aNrmMin
407 aNrm, HNrm = np.meshgrid(aGrid, HabitGrid, indexing="ij")
409 # Construct the consumption function
410 if type(solution_next["dvdkFunc"]) is ConstantFunction:
411 # This is the terminal period, and the consumption function is to consume all
412 cFunc = IdentityFunction(i_dim=0, n_dims=2)
413 mNrmMin = 0.0
415 else:
416 # Evaluate end-of-period marginal value on those grids, then calculate chi
417 EndOfPrd_dvda = DiscFacEff * solution_next["dvdkFunc"](aNrm, HNrm)
418 EndOfPrd_dvdH = DiscFacEff * solution_next["dvdhFunc"](aNrm, HNrm)
419 chi = U.derinv(EndOfPrd_dvda - HabitRte * EndOfPrd_dvdH)
421 # Recover c and h using the FOC inverter, then find endogenous m gridpoints
422 cNrm, hNrm = FOCinverter(HNrm, chi)
423 mNrm = aNrm + cNrm
425 # Construct the unconstrained consumption as a Curvilinear2Dinterp
426 cNrm = np.concatenate((np.zeros((1, HabitGrid.size)), cNrm), axis=0)
427 mNrm = np.concatenate((aNrmMin * np.ones((1, HabitGrid.size)), mNrm), axis=0)
428 hBot = (
429 np.reshape(hNrm[0, :], (1, HabitGrid.size))
430 if HabitRte == 1.0
431 else np.reshape(HabitGrid / (1.0 - HabitRte), (1, HabitGrid.size))
432 )
433 hNrm = np.concatenate((hBot, hNrm), axis=0)
434 cFuncUnc_base = Curvilinear2DInterp(cNrm, mNrm, hNrm)
436 # Re-interpolate the curvilinear consumption function onto an ordinary grid
437 mGridDense = mXtraGrid + aNrmMin
438 mMesh, hMesh = np.meshgrid(mGridDense, hGridDense, indexing="ij")
439 cMesh = cFuncUnc_base(mMesh, hMesh)
440 cMesh = np.concatenate((np.zeros((mGridDense.size, 1)), cMesh), axis=1)
441 cMesh = np.concatenate((np.zeros((1, hGridDense.size + 1)), cMesh), axis=0)
442 cFuncUnc = BilinearInterp(
443 cMesh, np.insert(mGridDense, 0, aNrmMin), np.insert(hGridDense, 0, 0.0)
444 )
446 # Add the constrained consumption function to that
447 if (BoroCnstArt is not None) and (BoroCnstArt > -np.inf):
448 cFuncCnst_temp = LinearInterp([BoroCnstArt, BoroCnstArt + 1.0], [0.0, 1.0])
449 cFuncCnst = LinearInterpOnInterp1D(
450 [cFuncCnst_temp, cFuncCnst_temp], np.array([0.0, 1.0])
451 )
452 cFunc = LowerEnvelope2D(cFuncUnc, cFuncCnst)
453 mNrmMin = np.maximum(aNrmMin, BoroCnstArt)
454 else:
455 cFunc = cFuncUnc
456 mNrmMin = aNrmMin
458 # Calculate the natural borrowing constraint (lowest allowable beginning-of-period capital)
459 PermShkVals = IncShkDstn.atoms[0, :]
460 TranShkVals = IncShkDstn.atoms[1, :]
461 kNrmMin_cand = (mNrmMin - TranShkVals) / Rfree * (PermShkVals * PermGroFac)
462 kNrmMin = np.max(kNrmMin_cand)
464 # Make beginning-of-period state grids
465 kGrid = kNrmMin + aXtraGrid
466 kNrm, hPre = np.meshgrid(kGrid, HabitGrid, indexing="ij")
468 # Compute expected marginal value over income shocks from beginning-of-period states
469 dvdk, dvdh = expected(
470 calc_marg_values,
471 IncShkDstn,
472 args=(
473 kNrm,
474 hPre,
475 CRRA,
476 Rfree,
477 PermGroFac,
478 HabitWgt,
479 HabitRte,
480 DiscFacEff,
481 cFunc,
482 solution_next["dvdhFunc"],
483 ),
484 )
486 # Transform and package the marginal value functions
487 dvdkNvrs = np.concatenate((np.zeros((1, HabitGrid.size)), U.derinv(dvdk)), axis=0)
488 dvdkNvrsFunc = BilinearInterp(dvdkNvrs, np.insert(kGrid, 0, kNrmMin), HabitGrid)
489 dvdkFunc = MargValueFuncCRRA(dvdkNvrsFunc, CRRA)
490 dvdhNvrs = U.inv(dvdh)
491 dvdhNvrsFunc = BilinearInterp(dvdhNvrs, kGrid, HabitGrid)
492 dvdhFunc = ValueFuncCRRA(dvdhNvrsFunc, CRRA)
494 # Package the solution as a dictionary and return it
495 solution_now = {
496 "cFunc": cFunc,
497 "dvdkFunc": dvdkFunc,
498 "dvdhFunc": dvdhFunc,
499 "kNrmMin": kNrmMin,
500 "distance_criteria": ["cFunc"],
501 }
502 return solution_now
505###############################################################################
508def solve_optimal_share_habit(
509 solution_next,
510 RiskyDstn,
511 LivPrb,
512 DiscFac,
513 CRRA,
514 Rfree,
515 BoroCnstArt,
516 aXtraGrid,
517 HabitGrid,
518 ShareGrid,
519 ShareLimit,
520):
521 """
522 Solve only the portfolio allocation problem of the consumption-saving model
523 with habit formation.
525 Parameters
526 ----------
527 solution_next : dict
528 Dictionary with next period's solution.
529 RiskyDstn : DiscreteDistribution
530 Discretized risky asset return distribution.
531 LivPrb : float
532 Survival probability at the end of this period.
533 DiscFac : float
534 Intertemporal discount factor.
535 CRRA : float
536 Coefficient of relative risk aversion.
537 Rfree : float
538 Interest factor on risk-free asset.
539 BoroCnstArt : float or None
540 Artificial borrowing constraint on end-of-period assets.
541 aXtraGrid : np.array
542 Grid of "assets above minimum".
543 HabitGrid : np.array
544 Grid of habit stock values.
545 ShareGrid : np.array
546 Grid of risky share values on [0,1].
547 ShareLimit : float
548 Merton-Samuelson limiting share as wealth -> infinity.
550 Returns
551 -------
552 solution_mid : dict
553 Solution to the portfolio allocation problem, which has a similar form
554 as the solution to the period overall. The dictionary key dvdkFunc refers
555 to dvdwFunc, marginal value of post-consumption wealth. The key dvdHfunc
556 is the *mid-period* marginal value of post-consumption habit stock, just
557 before the portfolio allocation decision is made.
558 """
559 U = UtilityFuncCRRA(CRRA)
560 DiscFacEff = DiscFac * LivPrb
562 # Unpack next period's solution
563 dvdkFunc_next = solution_next["dvdkFunc"]
564 dvdhFunc_next = solution_next["dvdhFunc"]
566 # This solver assumes the post-consumption wealth state satisfies wNrm >= 0.
567 # Reject unsupported borrowing configurations rather than silently producing
568 # share policies on an inconsistent state space.
569 if not np.isclose(BoroCnstArt, 0.0):
570 raise NotImplementedError(
571 "solve_optimal_share_habit only supports BoroCnstArt == 0.0; "
572 "negative borrowing constraints are not incorporated into the "
573 "post-consumption wealth/share state space."
574 )
576 # Minimum post-consumption wealth in the supported no-borrowing case
577 wNrmMin = 0.0
579 # Handle the terminal period: share doesn't matter, marginal value still zero
580 if isinstance(dvdkFunc_next, ConstantFunction):
581 ShareFunc_mid = ConstantFunction(ShareLimit)
582 dvdwFunc = ConstantFunction(0.0)
583 dvdHmidFunc = ConstantFunction(0.0)
584 solution_mid = {
585 "ShareFunc": ShareFunc_mid,
586 "dvdkFunc": dvdwFunc,
587 "dvdhFunc": dvdHmidFunc,
588 "kNrmMin": wNrmMin,
589 }
590 return solution_mid
592 # For each (w, H, s), compute E_R[Rport * dvdkFunc_next(Rport*w, H)]
593 # and E_R[dvdhFunc_next(Rport*w, H)] and
594 # E_R[(Risky-Rfree)*w * dvdkFunc_next(Rport*w, H)].
596 # Build 3D meshes (w, H, s). These are for mid-period state space points
597 # combined with candidate risky share values
598 wGrid = aXtraGrid + wNrmMin
599 wCount = wGrid.size
600 HabitCount = HabitGrid.size
601 ShareCount = ShareGrid.size
602 wMesh, Hmesh, sMesh = np.meshgrid(wGrid, HabitGrid, ShareGrid, indexing="ij")
604 # Compute expected marginal value of wealth, habit stock, and risky share for each (w,H,S)
605 dvds_mid = DiscFacEff * expected(
606 calc_mid_dvds,
607 RiskyDstn,
608 args=(wMesh, Hmesh, sMesh, Rfree, dvdkFunc_next),
609 )
611 # For each (w, H), find optimal share where dvds == 0 by looking for a
612 # sign change: dvds goes from positive to negative
613 focs = dvds_mid
614 crossing = np.logical_and(focs[:, :, 1:] <= 0.0, focs[:, :, :-1] >= 0.0)
615 share_idx = np.argmax(crossing, axis=2)
617 # Find the optimal risky share by solving a linear equation for the FOC
618 # crossing point, given that we know upper and lower bounding points for it
619 w_idx, h_idx = np.meshgrid(np.arange(wCount), np.arange(HabitCount), indexing="ij")
620 bot_s = ShareGrid[share_idx]
621 top_s = ShareGrid[np.minimum(share_idx + 1, ShareCount - 1)]
622 bot_f = focs[w_idx, h_idx, share_idx]
623 top_f = focs[w_idx, h_idx, np.minimum(share_idx + 1, ShareCount - 1)]
624 alpha_interp = np.where(
625 (top_f - bot_f) != 0.0,
626 1.0 - top_f / (top_f - bot_f),
627 0.5,
628 )
629 Share_opt = (1.0 - alpha_interp) * bot_s + alpha_interp * top_s
631 # Handle corner solutions
632 constrained_top = focs[:, :, -1] > 0.0
633 constrained_bot = focs[:, :, 0] < 0.0
634 Share_opt[constrained_top] = 1.0
635 Share_opt[constrained_bot] = 0.0
637 # Construct the share function over (w,H)
638 ShareFunc_by_HNrm = [
639 LinearInterp(wGrid, Share_opt[:, j], ShareLimit, 0.0, lower_extrap=True)
640 for j in range(HabitCount)
641 ]
642 ShareFunc_mid = LinearInterpOnInterp1D(ShareFunc_by_HNrm, HabitGrid)
644 # Extract optimized end-of-period marginal values at optimal share
645 wNrm, HNrm = np.meshgrid(wGrid, HabitGrid, indexing="ij")
646 dvdw_opt, dvdH_opt = DiscFacEff * expected(
647 calc_mid_dvdx,
648 RiskyDstn,
649 args=(wNrm, HNrm, Share_opt, Rfree, dvdkFunc_next, dvdhFunc_next),
650 )
652 # Build interpolant for mid-period marginal value of habit stock on (w, H) grid;
653 # dvdH_opt already includes DiscFacEff.
654 dvdHnvrs = U.inv(dvdH_opt)
655 dvdHNvrsFunc = BilinearInterp(dvdHnvrs, wGrid, HabitGrid)
656 dvdHmidFunc = ValueFuncCRRA(dvdHNvrsFunc, CRRA)
658 # Build interpolant for mid-period marginal value of wealth (after consumption
659 # but before portfolio allocation); incorporates DiscFacEff
660 dvdwNvrs = np.concatenate(
661 (np.zeros((1, HabitGrid.size)), U.derinv(dvdw_opt)), axis=0
662 )
663 dvdwNvrsFunc = BilinearInterp(dvdwNvrs, np.insert(wGrid, 0, wNrmMin), HabitGrid)
664 dvdwFunc = MargValueFuncCRRA(dvdwNvrsFunc, CRRA)
666 # Package and return the mid-period solution as a dictionary
667 solution_mid = {
668 "ShareFunc": ShareFunc_mid,
669 "dvdkFunc": dvdwFunc,
670 "dvdhFunc": dvdHmidFunc,
671 "kNrmMin": wNrmMin,
672 }
673 return solution_mid
676###############################################################################
679def solve_one_period_HabitPortfolio(
680 solution_next,
681 IncShkDstn,
682 RiskyDstn,
683 LivPrb,
684 DiscFac,
685 CRRA,
686 Rfree,
687 PermGroFac,
688 BoroCnstArt,
689 aXtraGrid,
690 HabitGrid,
691 ShareGrid,
692 ShareLimit,
693 FOCinverter,
694 HabitWgt,
695 HabitRte,
696 mXtraGrid,
697 hGridDense,
698):
699 """
700 Solve one period of the consumption-saving model with habit formation and
701 portfolio choice. This version uses a "modular" structure, solving the portfolio
702 problem and consumption problem using separate sub-functions. The latter sub-
703 function is literally the basic habit formation solver with some parameters
704 turned off (i.e. set to 1.0) because they were handled in the portfolio block.
706 Parameters
707 ----------
708 solution_next : dict
709 Dictionary with next period's solution.
710 IncShkDstn : DiscreteDistribution
711 Discretized permanent and transitory income shock distribution this period.
712 RiskyDstn : DiscreteDistribution
713 Discretized risky asset return distribution.
714 LivPrb : float
715 Survival probability at the end of this period.
716 DiscFac : float
717 Intertemporal discount factor.
718 CRRA : float
719 Coefficient of relative risk aversion.
720 Rfree : float
721 Interest factor on capital at the start of this period.
722 PermGroFac : float
723 Permanent income growth factor at the start of this period.
724 BoroCnstArt : float or None
725 Artificial borrowing constraint on assets at the end of this period,
726 as a fraction of permanent income.
727 aXtraGrid : np.array
728 Grid of "assets above minimum".
729 HabitGrid : np.array
730 Grid of consumption habit stocks on which to solve the problem.
731 ShareGrid : np.array
732 Grid of risky share values on [0,1].
733 ShareLimit : float
734 Merton-Samuelson limiting share as wealth -> infinity.
735 FOCinverter : HabitFormationInverter
736 Function that inverts the first order conditions to yield optimal consumption
737 and the decision-time habit stock from which it was chosen.
738 HabitWgt : float
739 Exponent on habit stock, which is used as a divisor on consumption in
740 the utility function: U(c,h) = u(c / h**alpha). Should be on unit interval.
741 HabitRte : float
742 Rate at which habit stock is updated by new consumption: H = lambda*c + (1-lambda)*h.
743 Should be on the unit interval.
744 mXtraGrid : np.array
745 Dense grid of market resources, used to "re-interpolate" the curvilinear
746 consumption function onto a rectilinear grid.
747 hGridDense : np.array
748 Dense grid of habit stocks, used to "re-interpolate" the curvilinear
749 consumption function onto a rectilinear grid.
751 Returns
752 -------
753 solution_now : dict
754 Solution to this period's problem, with the following keys:
755 cFunc : Consumption function over (mNrm, hNrm).
756 ShareFunc : Risky share function over (wNrm, HNrm).
757 dvdkFunc : Marginal value of beginning-of-period capital, defined on (kNrm, hPre).
758 dvdhFunc : Marginal value of beginning-of-period habit stock, defined on (kNrm, hPre).
759 kNrmMin : Minimum allowable beginning-of-period capital.
760 """
761 # Solve the portfolio allocation problem, yielding the "mid period solution"
762 solution_mid = solve_optimal_share_habit(
763 solution_next,
764 RiskyDstn,
765 LivPrb,
766 DiscFac,
767 CRRA,
768 Rfree,
769 BoroCnstArt,
770 aXtraGrid,
771 HabitGrid,
772 ShareGrid,
773 ShareLimit,
774 )
776 # Solve the consumption-saving problem, yielding the solution to this period
777 solution_now = solve_one_period_ConsHabit(
778 solution_mid,
779 IncShkDstn,
780 1.0, # LivPrb accounted for above, turn off
781 1.0, # DiscFac accounted for above, turn off
782 CRRA,
783 1.0, # Rfree accounted for above, turn off
784 PermGroFac,
785 BoroCnstArt,
786 aXtraGrid,
787 HabitGrid,
788 FOCinverter,
789 HabitWgt,
790 HabitRte,
791 mXtraGrid,
792 hGridDense,
793 )
795 # Add the risky share function to the period's solution and return it
796 solution_now["ShareFunc"] = solution_mid["ShareFunc"]
797 return solution_now
800###############################################################################
802# Make a dictionary of constructors for the habit formation model
803HabitConsumerType_constructors_default = {
804 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
805 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
806 "HabitInitDstn": make_lognormal_habit_init_dstn,
807 "IncShkDstn": construct_lognormal_income_process_unemployment,
808 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
809 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
810 "aXtraGrid": make_assets_grid,
811 "FOCinverter": make_habit_inverter,
812 "HabitGrid": make_habit_grid,
813 "DenseGrids": make_dense_grids,
814 "mXtraGrid": get_it_from("DenseGrids"),
815 "hGridDense": get_it_from("DenseGrids"),
816 "solution_terminal": make_habit_solution_terminal,
817}
819# Make a dictionary with parameters for the default constructor for kNrmInitDstn
820HabitConsumerType_kNrmInitDstn_default = {
821 "kLogInitMean": -12.0, # Mean of log initial capital
822 "kLogInitStd": 0.0, # Stdev of log initial capital
823 "kNrmInitCount": 15, # Number of points in initial capital discretization
824}
826# Make a dictionary with parameters for the default constructor for pLvlInitDstn
827HabitConsumerType_pLvlInitDstn_default = {
828 "pLogInitMean": 0.0, # Mean of log permanent income
829 "pLogInitStd": 0.0, # Stdev of log permanent income
830 "pLvlInitCount": 15, # Number of points in initial capital discretization
831}
833# Make a dictionary with parameters for the default constructor for HabitInitDstn
834HabitConsumerType_HabitInitDstn_default = {
835 "hLogInitMean": -0.5, # Mean of log habit stock
836 "hLogInitStd": 0.2, # Stdev of log initial habit stock
837 "HabitInitCount": 15, # Number of points in initial habit stock discretization
838}
840# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
841HabitConsumerType_IncShkDstn_default = {
842 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
843 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
844 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
845 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
846 "UnempPrb": 0.05, # Probability of unemployment while working
847 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
848 "T_retire": 0, # Period of retirement (0 --> no retirement)
849 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
850 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
851}
853# Default parameters to make aXtraGrid using make_assets_grid
854HabitConsumerType_aXtraGrid_default = {
855 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
856 "aXtraMax": 30, # Maximum end-of-period "assets above minimum" value
857 "aXtraNestFac": 2, # Exponential nesting factor for aXtraGrid
858 "aXtraCount": 72, # Number of points in the grid of "assets above minimum"
859 "aXtraExtra": None, # Additional other values to add in grid (optional)
860}
862# Default parameters to make HabitGrid using make_habit_grid
863HabitConsumerType_HabitGrid_default = {
864 "HabitMin": 0.2,
865 "HabitMax": 5.0,
866 "HabitCount": 41,
867 "HabitOrder": 2.0,
868}
870# Default parameters to make the FOC inverter
871HabitConsumerType_inverter_default = {
872 "ChiMax": 50.0,
873 "ChiCount": 251,
874 "ChiOrder": 1.5,
875}
877# Make a dictionary to specify an habit formation consumer type
878HabitConsumerType_solving_default = {
879 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
880 "cycles": 1, # Finite, non-cyclic model
881 "T_cycle": 1, # Number of periods in the cycle for this agent type
882 "pseudo_terminal": True, # It's a fake stub
883 "constructors": HabitConsumerType_constructors_default, # See dictionary above
884 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
885 "CRRA": 2.0, # Coefficient of relative risk aversion
886 "Rfree": [1.03], # Interest factor on retained assets
887 "DiscFac": 0.96, # Intertemporal discount factor
888 "LivPrb": [0.98], # Survival probability after each period
889 "PermGroFac": [1.01], # Permanent income growth factor
890 "BoroCnstArt": 0.0, # Artificial borrowing constraint
891 "HabitWgt": 0.5, # Weight on consumption habit; exponent on habit divisor in utility
892 "HabitRte": 0.2, # Speed of consumption habit updating
893 "DenseFactor": 3, # Density factor for re-interpolation
894}
895HabitConsumerType_simulation_default = {
896 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
897 "AgentCount": 10000, # Number of agents of this type
898 "T_age": None, # Age after which simulated agents are automatically killed
899}
901HabitConsumerType_defaults = {}
902HabitConsumerType_defaults.update(HabitConsumerType_IncShkDstn_default)
903HabitConsumerType_defaults.update(HabitConsumerType_kNrmInitDstn_default)
904HabitConsumerType_defaults.update(HabitConsumerType_pLvlInitDstn_default)
905HabitConsumerType_defaults.update(HabitConsumerType_HabitInitDstn_default)
906HabitConsumerType_defaults.update(HabitConsumerType_aXtraGrid_default)
907HabitConsumerType_defaults.update(HabitConsumerType_HabitGrid_default)
908HabitConsumerType_defaults.update(HabitConsumerType_inverter_default)
909HabitConsumerType_defaults.update(HabitConsumerType_solving_default)
910HabitConsumerType_defaults.update(HabitConsumerType_simulation_default)
913class HabitConsumerType(AgentType):
914 r"""
915 A class for representing consumers who form consumption habits. Agents get
916 flow utility according to a CRRA felicity function that depends on both current
917 consumption and the habit stock h_t. The habit stock evolves as a weighted
918 average of current consumption and prior habit stock. Consumers can save in
919 a single risk-free asset, so this model is an extension of the workhorse
920 IndShockConsumerType to include a habit stock.
922 .. math::
923 \newcommand{\CRRA}{\rho}
924 \newcommand{\LivPrb}{\mathsf{S}}
925 \newcommand{\PermGroFac}{\Gamma}
926 \newcommand{\Rfree}{\mathsf{R}}
927 \newcommand{\DiscFac}{\beta}
928 \newcommand{\HabitWgt}{\alpha}
929 \newcommand{\HabitRte}{\lambda}
931 \begin{align*}
932 v_t(m_t,h_t) &= \max_{c_t}u(c_t,h_t) + \DiscFac \LivPrb_t \mathbb{E}_{t} \left[ (\PermGroFac_{t+1} \psi_{t+1})^{(1-\HabitWgt)(1-\CRRA)} v_{t+1}(m_{t+1}, h_{t+1}) \right] \\
933 & \text{s.t.} \\
934 a_t &= m_t - c_t, \\
935 H_t &= \HabitRte c_t + (1-\HabitRte) h_t, \\
936 a_t &\geq \underline{a}, \\
937 m_{t+1} &= a_t \Rfree_{t+1}/(\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
938 h_{t+1} &= H_t / (\PermGroFac_{t+1} \psi_{t+1}), \\
939 (\psi_{t+1},\theta_{t+1}) &\sim F_{t+1}, \\
940 \mathbb{E}[\psi] &= 1, \\
941 u(c,h) &= \frac{(c/h^\HabitWgt)^{1-\CRRA}}{1-\CRRA}.
942 \end{align*}
943 """
945 default_ = {
946 "params": HabitConsumerType_defaults,
947 "solver": solve_one_period_ConsHabit,
948 "model": "ConsHabit.yaml",
949 "track_vars": ["aNrm", "cNrm", "mNrm", "hNrm", "pLvl"],
950 }
952 time_inv_ = [
953 "DiscFac",
954 "CRRA",
955 "BoroCnstArt",
956 "aXtraGrid",
957 "HabitGrid",
958 "FOCinverter",
959 "HabitWgt",
960 "HabitRte",
961 "mXtraGrid",
962 "hGridDense",
963 ]
964 time_vary_ = ["IncShkDstn", "Rfree", "PermGroFac", "LivPrb"]
966 shock_vars_ = ["PermShk", "TranShk"]
967 distributions = [
968 "IncShkDstn",
969 "PermShkDstn",
970 "TranShkDstn",
971 "kNrmInitDstn",
972 "pLvlInitDstn",
973 "HabitInitDstn",
974 ]
977###############################################################################
979# Make a dictionary of constructors for the habit-portfolio model
980HabitPortfolio_constructors_default = HabitConsumerType_constructors_default.copy()
981HabitPortfolio_additional_constructors = {
982 "RiskyDstn": make_lognormal_RiskyDstn,
983 "ShareGrid": make_simple_ShareGrid,
984 "ShareLimit": calc_ShareLimit_for_CRRA,
985}
986HabitPortfolio_constructors_default.update(HabitPortfolio_additional_constructors)
988HabitPortfolio_RiskyDstn_default = {
989 "RiskyAvg": 1.08,
990 "RiskyStd": 0.18,
991 "RiskyCount": 5,
992}
994HabitPortfolio_ShareGrid_default = {
995 "ShareCount": 26,
996}
998HabitPortfolioConsumerType_defaults = HabitConsumerType_defaults.copy()
999HabitPortfolioConsumerType_defaults["constructors"] = (
1000 HabitPortfolio_constructors_default
1001)
1002HabitPortfolioConsumerType_defaults.update(HabitPortfolio_RiskyDstn_default)
1003HabitPortfolioConsumerType_defaults.update(HabitPortfolio_ShareGrid_default)
1006class HabitPortfolioConsumerType(HabitConsumerType):
1007 r"""
1008 A class for representing consumers who form consumption habits and can allocate
1009 their wealth between a risky and riskless asset. Agents get flow utility according
1010 to a CRRA felicity function that depends on both current consumption and the habit
1011 stock h_t. The habit stock evolves as a weighted average of current consumption
1012 and prior habit stock.
1014 This type's solver uses a "modular" approach that separates the portfolio allo-
1015 cation and consumption-saving problems into two functions. The latter function
1016 is just the solver for HabitConsumerType. Consequently, the consumption function
1017 is defined over (m,h) but the risky share function is defined over (w,H).
1019 .. math::
1020 \newcommand{\CRRA}{\rho}
1021 \newcommand{\LivPrb}{\mathsf{S}}
1022 \newcommand{\PermGroFac}{\Gamma}
1023 \newcommand{\Rfree}{\mathsf{R}}
1024 \newcommand{\Risky}{\mathfrak{R}}
1025 \newcommand{\DiscFac}{\beta}
1026 \newcommand{\HabitWgt}{\alpha}
1027 \newcommand{\HabitRte}{\lambda}
1029 \begin{align*}
1030 v_t(m_t,h_t) &= \max_{c_t, s_t} u(c_t,h_t) + \DiscFac \LivPrb_t
1031 \mathbb{E}_{t} \left[ (\PermGroFac_{t+1} \psi_{t+1})^{(1-\HabitWgt)(1-\CRRA)}
1032 v_{t+1}(m_{t+1}, h_{t+1}) \right] \\
1033 & \text{s.t.} \\
1034 w_t &= m_t - c_t, \\
1035 H_t &= \HabitRte c_t + (1-\HabitRte) h_t, \\
1036 w_t &\geq 0, \\
1037 s_t &\in [0,1], \\
1038 a_t &= R_t w_t, \\
1039 R_{t} &= s_t \Risky_{t} + (1-s_t) \Rfree_{t}, \\
1040 m_{t+1} &= a_t / (\PermGroFac_{t+1} \psi_{t+1}) + \theta_{t+1}, \\
1041 h_{t+1} &= H_t / (\PermGroFac_{t+1} \psi_{t+1}), \\
1042 (\psi_{t+1}, \theta_{t+1}) &\sim F_{t+1}, \\
1043 \Risky_{t} & \sim G, \\
1044 u(c,h) &= \frac{(c/h^\HabitWgt)^{1-\CRRA}}{1-\CRRA}.
1045 \end{align*}
1046 """
1048 default_ = {
1049 "params": HabitPortfolioConsumerType_defaults,
1050 "solver": solve_one_period_HabitPortfolio,
1051 "model": "ConsHabitPortfolio.yaml",
1052 "track_vars": ["aNrm", "cNrm", "mNrm", "hNrm", "Share", "pLvl"],
1053 }
1055 time_inv_ = HabitConsumerType.time_inv_ + ["RiskyDstn", "ShareGrid"]
1056 time_vary_ = HabitConsumerType.time_vary_ + ["ShareLimit"]
1057 shock_vars_ = HabitConsumerType.shock_vars_ + ["Risky"]
1058 distributions = HabitConsumerType.distributions + ["RiskyDstn"]