Coverage for HARK / ConsumptionSaving / ConsMedModel.py: 98%
560 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-08 05:31 +0000
1"""
2Consumption-saving models that also include medical spending.
3"""
5from copy import deepcopy
7import numpy as np
8from scipy.stats import norm
9from scipy.special import erfc
11from HARK import AgentType
12from HARK.Calibration.Income.IncomeProcesses import (
13 construct_lognormal_income_process_unemployment,
14 get_PermShkDstn_from_IncShkDstn,
15 get_TranShkDstn_from_IncShkDstn,
16 make_AR1_style_pLvlNextFunc,
17 make_pLvlGrid_by_simulation,
18 make_basic_pLvlPctiles,
19 make_persistent_income_process_dict,
20)
21from HARK.ConsumptionSaving.ConsIndShockModel import (
22 make_lognormal_kNrm_init_dstn,
23 make_lognormal_pLvl_init_dstn,
24)
25from HARK.ConsumptionSaving.ConsGenIncProcessModel import (
26 PersistentShockConsumerType,
27 VariableLowerBoundFunc2D,
28)
29from HARK.distributions import (
30 Lognormal,
31 MultivariateLogNormal,
32 add_discrete_outcome,
33 expected,
34)
35from HARK.interpolation import (
36 BilinearInterp,
37 BilinearInterpOnInterp1D,
38 ConstantFunction,
39 CubicInterp,
40 LinearInterp,
41 LinearInterpOnInterp1D,
42 LowerEnvelope3D,
43 MargMargValueFuncCRRA,
44 MargValueFuncCRRA,
45 TrilinearInterp,
46 UpperEnvelope,
47 ValueFuncCRRA,
48 VariableLowerBoundFunc3D,
49)
50from HARK.metric import MetricObject
51from HARK.rewards import (
52 CRRAutility,
53 CRRAutilityP,
54 CRRAutility_inv,
55 CRRAutility_invP,
56 CRRAutilityP_inv,
57 CRRAutilityPP,
58 UtilityFuncCRRA,
59)
60from HARK.utilities import NullFunc, make_grid_exp_mult, make_assets_grid, get_it_from
62__all__ = [
63 "MedShockConsumerType",
64 "MedExtMargConsumerType",
65 "make_lognormal_MedShkDstn",
66 "make_continuous_MedShockDstn",
67]
69utility_inv = CRRAutility_inv
70utilityP_inv = CRRAutilityP_inv
71utility = CRRAutility
72utility_invP = CRRAutility_invP
73utilityPP = CRRAutilityPP
76class TransConShareFunc(MetricObject):
77 """
78 A class for representing the "transformed consumption share" function with
79 "double CRRA" utility. Instances of this class take as inputs xLvl and MedShkEff
80 and return a transformed consumption share q. By definition, optimal consumption
81 cLvl = xLvl/(1 + exp(-q)). MedShkEff = MedShk*MedPrice.
82 """
84 distance_criteria = ["CRRA", "CRRAmed"]
86 def __init__(self, CRRA, CRRAmed):
87 Qhalf = make_grid_exp_mult(0.0, 30.0, 40, 3)
88 Q = np.unique(np.concatenate([-Qhalf, Qhalf]))
89 f = lambda q: q + (1.0 - CRRA / CRRAmed) * np.log(1.0 + np.exp(-q))
90 fp = lambda q: 1.0 - (1.0 - CRRA / CRRAmed) * np.exp(-q) / (1.0 + np.exp(-q))
91 G = f(Q)
92 D = 1.0 / fp(Q)
93 f_inv = CubicInterp(G, Q, D, lower_extrap=True)
95 self.CRRA = CRRA
96 self.CRRAmed = CRRAmed
97 self.f_inv = f_inv
98 self.coeff_x = 1.0 - CRRA / CRRAmed
99 self.coeff_Shk = 1.0 - 1.0 / CRRAmed
101 def __call__(self, xLvl, MedShkEff):
102 """
103 Evaluate the "transformed consumption share" function.
105 Parameters
106 ----------
107 xLvl : np.array
108 Array of expenditure levels.
109 MedShkEff : np.array
110 Identically shaped array of effective medical need shocks: MedShk*MedPrice.
112 Returns
113 -------
114 q : np.array
115 Identically shaped array of transformed consumption shares.
116 """
117 q = self.f_inv(self.coeff_x * np.log(xLvl) - self.coeff_Shk * np.log(MedShkEff))
118 return q
121def make_qFunc(CRRA, CRRAmed):
122 """
123 Basic constructor that makes the transformed-consumption-share function
124 from primitive parameters CRRA and CRRAmed.
125 """
126 return TransConShareFunc(CRRA, CRRAmed)
129class cAndMedFunc(MetricObject):
130 """
131 A class representing the consumption and medical care function based on the total
132 expenditure function. Its call function returns cLvl, MedLvl, and xLvl.
133 Also has functions for the two main controls individually.
134 """
136 def __init__(self, xFunc, qFunc, MedShift, MedPrice):
137 """
138 Constructor method for a new instance of cAndMedFunc.
140 Parameters
141 ----------
142 xFunc : function
143 Expenditure function (cLvl & MedLvl), defined over (mLvl,pLvl,MedShk).
144 qFunc : function
145 Transformed consumption share function, defined over (xLvl, MedShkEff).
146 MedShift : float
147 Shifter term for medical care, representing the quantity of care that the
148 agent "gets for free" automatically-- self care, perhaps.
149 MedPrice : float
150 Relative price of a unit of medical care.
152 Returns
153 -------
154 None
155 """
156 # Store the data
157 self.xFunc = xFunc
158 self.qFunc = qFunc
159 self.MedShift = MedShift
160 self.MedPrice = MedPrice
161 self.update()
163 def update(self):
164 # Update "constructed" attributes
165 self.CRRA = self.qFunc.CRRA
166 self.CRRAmed = self.qFunc.CRRAmed
167 self.factor = self.MedPrice ** (1.0 / self.CRRA) * self.MedShift ** (
168 self.CRRAmed / self.CRRA
169 )
171 def __call__(self, mLvl, pLvl, MedShk):
172 """
173 Evaluates the policy function and returns cLvl and MedLvl.
175 Parameters
176 ----------
177 mLvl : np.array
178 Array of market resources values.
179 pLvl : np.array
180 Array of permanent income levels.
181 MedShk : np.array
182 Array of medical need shocks.
184 Returns
185 -------
186 cLvl : np.array
187 Array of consumption levels.
188 MedLvl : np.array
189 Array of medical care levels.
190 xLvl : np.array
191 Array of total expenditure levels.
192 """
193 if type(mLvl) is float:
194 mLvl = np.array([mLvl])
195 pLvl = np.array([pLvl])
196 MedShk = np.array([MedShk])
197 singleton = True
198 else:
199 singleton = False
201 # Initialize the output arrays
202 cLvl = np.zeros_like(mLvl)
203 MedLvl = np.zeros_like(mLvl)
205 # Evaluate total expenditure
206 xLvl = self.xFunc(mLvl, pLvl, MedShk)
208 # Determine which inputs should buy zero medical care
209 xLvlCrit = self.factor * MedShk ** ((1.0 - self.CRRAmed) / self.CRRA)
210 ZeroShk = MedShk == 0.0
211 ZeroMed = np.logical_or(xLvl <= xLvlCrit, ZeroShk)
212 SomeMed = np.logical_not(ZeroMed)
214 # Fill in consumption for those who buy zero medical care
215 cLvl[ZeroMed] = xLvl[ZeroMed]
217 # Calculate consumption and medical care for those who buy some care
218 xLvl_temp = xLvl[SomeMed] + self.MedShift * self.MedPrice
219 q = self.qFunc(xLvl_temp, MedShk[SomeMed] * self.MedPrice)
220 z = np.exp(-q)
221 cLvl[SomeMed] = xLvl_temp / (1.0 + z)
222 MedLvl[SomeMed] = xLvl_temp / self.MedPrice * (z / (1.0 + z)) - self.MedShift
224 # Make sure MedLvl is non-negative (can clip over by 1e-14 sometimes)
225 MedLvl = np.maximum(MedLvl, 0.0)
227 # If the inputs were singletons, extract them before returning output
228 if singleton:
229 return cLvl[0], MedLvl[0], xLvl[0]
230 else:
231 return cLvl, MedLvl, xLvl
233 def cFunc(self, mLvl, pLvl, MedShk):
234 cLvl, MedLvl, xLvl = self(mLvl, pLvl, MedShk)
235 return cLvl
237 def MedFunc(self, mLvl, pLvl, MedShk):
238 cLvl, MedLvl, xLvl = self(mLvl, pLvl, MedShk)
239 return MedLvl
242def make_market_resources_grid(mNrmMin, mNrmMax, mNrmNestFac, mNrmCount, mNrmExtra):
243 """
244 Constructor for mNrmGrid that aliases make_assets_grid.
245 """
246 return make_assets_grid(mNrmMin, mNrmMax, mNrmCount, mNrmExtra, mNrmNestFac)
249def make_capital_grid(kLvlMin, kLvlMax, kLvlCount, kLvlOrder):
250 """
251 Constructor for kLvlGrid, using a simple "invertible" format.
252 """
253 base_grid = np.linspace(0.0, 1.0, kLvlCount) ** kLvlOrder
254 kLvlGrid = (kLvlMax - kLvlMin) * base_grid + kLvlMin
255 return kLvlGrid
258def reformat_bequest_motive(BeqMPC, BeqInt, CRRA):
259 """
260 Reformats interpretable bequest motive parameters (terminal intercept and MPC)
261 into parameters that are easily useable in math (shifter and scaler).
262 """
263 BeqParamDict = {
264 "BeqFac": BeqMPC ** (-CRRA),
265 "BeqShift": BeqInt / BeqMPC,
266 }
267 return BeqParamDict
270def make_lognormal_MedShkDstn(
271 T_cycle,
272 MedShkAvg,
273 MedShkStd,
274 MedShkCount,
275 MedShkCountTail,
276 RNG,
277 MedShkTailBound=[0.0, 0.9],
278 MedShkZeroPrb=[0.0],
279):
280 r"""
281 Constructs discretized lognormal distributions of medical preference shocks
282 for each period in the cycle.
284 .. math::
285 \text{ medShk}_t \sim \exp(\mathcal{N}(\textbf{MedShkStd}^2)) \\
286 \mathbb{E}[\text{medShk}_t]=\textbf{MedShkAvg}
289 Parameters
290 ----------
291 T_cycle : int
292 Number of non-terminal periods in the agent's cycle.
293 MedShkAvg : [float]
294 Mean of non-zero medical needs shock in each period of the problem.
295 MedShkStd : [float]
296 Standard deviation of log (non-zero) medical needs shock in each period of the problem.
297 MedShkCount : int
298 Number of equiprobable nodes in the "body" of the discretization.
299 MedShkCountTail : int
300 Number of nodes in each "tail" of the discretization.
301 RNG : RandomState
302 The AgentType's internal random number generator.
303 MedShkTailBound : [float,float]
304 CDF bounds for the tail of the discretization.
305 MedShkZeroPrb : [float]
306 Probability of getting a zero medical need shock in each period (default zero).
308 Returns
309 -------
310 MedShkDstn : [DiscreteDistribuion]
311 """
312 MedShkDstn = [] # empty list for medical shock distribution each period
313 for t in range(T_cycle):
314 # get shock distribution parameters
315 MedShkAvg_t = MedShkAvg[t]
316 MedShkStd_t = MedShkStd[t]
317 MedShkDstn_t = Lognormal(
318 mu=np.log(MedShkAvg_t) - 0.5 * MedShkStd_t**2,
319 sigma=MedShkStd_t,
320 seed=RNG.integers(0, 2**31 - 1),
321 ).discretize(
322 N=MedShkCount,
323 method="equiprobable",
324 tail_N=MedShkCountTail,
325 tail_bound=MedShkTailBound,
326 )
327 MedShkDstn_t = add_discrete_outcome(
328 MedShkDstn_t, 0.0, MedShkZeroPrb[t], sort=True
329 ) # add point at zero
330 MedShkDstn.append(MedShkDstn_t)
331 return MedShkDstn
334def make_continuous_MedShockDstn(
335 MedShkLogMean, MedShkLogStd, MedCostLogMean, MedCostLogStd, MedCorr, T_cycle, RNG
336):
337 """
338 Construct a time-varying list of bivariate lognormals for the medical shocks
339 distribution. This representation uses fully continuous distributions, with
340 no discretization in either dimension.
342 Parameters
343 ----------
344 MedShkLogMean : [float]
345 Age-varying list of means of log medical needs (utility) shocks.
346 MedShkLogStd : [float]
347 Age-varying list of standard deviations of log medical needs (utility) shocks.
348 MedCostLogMean : [float]
349 Age-varying list of means of log medical expense shocks.
350 MedCostLogStd : [float]
351 Age-varying list of standard deviations of log medical expense shocks..
352 MedCorr : [float]
353 Age-varying correlation coefficient between log medical expenses and utility shocks.
354 T_cycle : int
355 Number of periods in the agent's sequence.
356 RNG : RandomState
357 Random number generator for this type.
359 Returns
360 -------
361 MedShockDstn : [MultivariateLognormal]
362 Age-varying list of bivariate lognormal distributions, ordered as (MedCost,MedShk).
363 """
364 MedShockDstn = []
365 for t in range(T_cycle):
366 s1 = MedCostLogStd[t]
367 s2 = MedShkLogStd[t]
368 diag = MedCorr[t] * s1 * s2
369 S = np.array([[s1**2, diag], [diag, s2**2]])
370 M = np.array([MedCostLogMean[t], MedShkLogMean[t]])
371 seed_t = RNG.integers(0, 2**31 - 1)
372 dstn_t = MultivariateLogNormal(mu=M, Sigma=S, seed=seed_t)
373 MedShockDstn.append(dstn_t)
374 return MedShockDstn
377def make_MedShock_solution_terminal(
378 CRRA,
379 CRRAmed,
380 MedShkDstn,
381 MedPrice,
382 MedShift,
383 aXtraGrid,
384 pLvlGrid,
385 qFunc,
386 CubicBool,
387):
388 """
389 Construct the terminal period solution for this type. Similar to other models,
390 optimal behavior involves spending all available market resources; however,
391 the agent must split his resources between consumption and medical care.
392 """
393 # Take last period data
394 MedPrice_T = MedPrice[-1]
395 MedShkDstn_T = MedShkDstn[-1]
396 pLvlGrid_T = pLvlGrid[-1]
397 MedShift_T = MedShift[-1]
399 # Make the policy functions for the terminal period
400 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
401 xFunc_terminal = TrilinearInterp(
402 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
403 trivial_grid,
404 trivial_grid,
405 trivial_grid,
406 )
407 PolicyFunc_terminal = cAndMedFunc(
408 xFunc_terminal,
409 qFunc,
410 MedShift_T,
411 MedPrice_T,
412 )
414 # Make state space grids for the terminal period
415 xLvlMin = np.min(aXtraGrid) * np.min(pLvlGrid_T)
416 xLvlMax = np.max(aXtraGrid) * np.max(pLvlGrid_T)
417 xLvlGrid = make_grid_exp_mult(xLvlMin, xLvlMax, 3 * aXtraGrid.size, 4)
418 MedShkGrid = MedShkDstn_T.atoms[0]
419 MedShkPrbs = MedShkDstn_T.pmv
421 # Calculate optimal consumption on a grid of market resources and medical shocks
422 mLvlGrid = xLvlGrid
423 mLvlGrid_tiled = np.tile(
424 np.reshape(mLvlGrid, (mLvlGrid.size, 1)), (1, MedShkGrid.size)
425 )
426 # Permanent income irrelevant in terminal period
427 pLvlGrid_tiled = np.ones_like(mLvlGrid_tiled)
428 MedShkGrid_tiled = np.tile(
429 np.reshape(MedShkGrid, (1, MedShkGrid.size)), (mLvlGrid.size, 1)
430 )
431 cLvlGrid, MedGrid, xTrash = PolicyFunc_terminal(
432 mLvlGrid_tiled, pLvlGrid_tiled, MedShkGrid_tiled
433 )
435 # Integrate marginal value across shocks to get expected marginal value
436 vPgrid = cLvlGrid ** (-CRRA)
437 vPgrid[np.isinf(vPgrid)] = 0.0 # correct for issue at bottom edges
438 PrbGrid = np.tile(np.reshape(MedShkPrbs, (1, MedShkGrid.size)), (mLvlGrid.size, 1))
439 vP_expected = np.sum(vPgrid * PrbGrid, axis=1)
441 # Construct the marginal (marginal) value function for the terminal period
442 vPnvrs = np.insert(vP_expected ** (-1.0 / CRRA), 0, 0.0)
443 vPnvrsFunc = BilinearInterp(
444 np.tile(np.reshape(vPnvrs, (vPnvrs.size, 1)), (1, trivial_grid.size)),
445 np.insert(mLvlGrid, 0, 0.0),
446 trivial_grid,
447 )
448 vPfunc_terminal = MargValueFuncCRRA(vPnvrsFunc, CRRA)
449 vPPfunc_terminal = MargMargValueFuncCRRA(vPnvrsFunc, CRRA)
451 # Integrate value across shocks to get expected value
452 vGrid = utility(cLvlGrid, rho=CRRA) + utility(
453 (MedGrid + MedShift_T) / MedShkGrid_tiled, rho=CRRAmed
454 )
455 vGrid[np.isinf(vGrid)] = 0.0 # correct for issue at bottom edges
456 v_expected = np.sum(vGrid * PrbGrid, axis=1)
458 # Construct the value function for the terminal period
459 vNvrs = np.insert(utility_inv(v_expected, rho=CRRA), 0, 0.0)
460 vNvrsP = vP_expected * utility_invP(v_expected, rho=CRRA)
461 vNvrsP = np.insert(vNvrsP, 0, 0.0)
462 # tempFunc = CubicInterp(np.insert(mLvlGrid, 0, 0.0), vNvrs, vNvrsP)
463 tempFunc = LinearInterp(np.insert(mLvlGrid, 0, 0.0), vNvrs)
464 vNvrsFunc = LinearInterpOnInterp1D([tempFunc, tempFunc], trivial_grid)
465 vFunc_terminal = ValueFuncCRRA(vNvrsFunc, CRRA)
467 # Make and return the terminal period solution
468 solution_terminal = {
469 "distance_criteria": ["vPfunc"],
470 "PolicyFunc": PolicyFunc_terminal,
471 "vFunc": vFunc_terminal,
472 "vPfunc": vPfunc_terminal,
473 "vPPfunc": vPPfunc_terminal,
474 "hLvl": ConstantFunction(0.0),
475 "mLvlMin": ConstantFunction(0.0),
476 }
477 return solution_terminal
480###############################################################################
483def solve_one_period_ConsMedShock(
484 solution_next,
485 IncShkDstn,
486 MedShkDstn,
487 LivPrb,
488 DiscFac,
489 CRRA,
490 CRRAmed,
491 MedShift,
492 Rfree,
493 MedPrice,
494 pLvlNextFunc,
495 BoroCnstArt,
496 aXtraGrid,
497 pLvlGrid,
498 vFuncBool,
499 CubicBool,
500 qFunc,
501):
502 """
503 Class for solving the one period problem for the "medical shocks" model, in
504 which consumers receive shocks to permanent and transitory income as well as
505 shocks to "medical need"-- utility shocks for a second good.
507 Parameters
508 ----------
509 solution_next : dict
510 The solution to next period's one period problem, represented as a dictionary.
511 IncShkDstn : distribution.Distribution
512 A discrete approximation to the income process between the period being
513 solved and the one immediately following (in solution_next).
514 MedShkDstn : distribution.Distribution
515 Discrete distribution of the need shifter for medical care.
516 LivPrb : float
517 Survival probability; likelihood of being alive at the beginning of
518 the succeeding period.
519 DiscFac : float
520 Intertemporal discount factor for future utility.
521 CRRA : float
522 Coefficient of relative risk aversion for composite consumption.
523 CRRAmed : float
524 Coefficient of relative risk aversion for medical care.
525 MedShift : float
526 Additive shifter on medical care in the utility function. Could represent
527 self-care for small medical needs that would not be observed / reported.
528 Rfree : float
529 Risk free interest factor on end-of-period assets.
530 MedPrice : float
531 Price of unit of medical care relative to unit of consumption.
532 pLvlNextFunc : float
533 Expected permanent income next period as a function of current pLvl.
534 BoroCnstArt: float or None
535 Borrowing constraint for the minimum allowable assets to end the
536 period with.
537 aXtraGrid: np.array
538 Array of "extra" end-of-period (normalized) asset values-- assets
539 above the absolute minimum acceptable level.
540 pLvlGrid: np.array
541 Array of permanent income levels at which to solve the problem.
542 vFuncBool: boolean
543 An indicator for whether the value function should be computed and
544 included in the reported solution.
545 CubicBool: boolean
546 An indicator for whether the solver should use cubic or linear interpolation.
547 qFunc : TransConShareFunc
548 Function that maps spending xLvl and "effective" medical needs shock MedShkEff
549 to a transformation of the consumption share.
551 Returns
552 -------
553 solution_now : dict
554 Solution to this period's consumption-saving problem, as a dictionary.
555 """
556 # Define the utility functions for this period
557 uCon = UtilityFuncCRRA(CRRA)
558 uMed = UtilityFuncCRRA(CRRAmed) # Utility function for normalized medical care
559 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
561 # Unpack next period's income shock distribution
562 ShkPrbsNext = IncShkDstn.pmv
563 PermShkValsNext = IncShkDstn.atoms[0]
564 TranShkValsNext = IncShkDstn.atoms[1]
565 PermShkMinNext = np.min(PermShkValsNext)
566 TranShkMinNext = np.min(TranShkValsNext)
567 MedShkPrbs = MedShkDstn.pmv
568 MedShkVals = MedShkDstn.atoms.flatten()
570 # Calculate the probability that we get the worst possible income draw
571 IncNext = PermShkValsNext * TranShkValsNext
572 WorstIncNext = PermShkMinNext * TranShkMinNext
573 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
574 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
576 # Unpack next period's (marginal) value function
577 vFuncNext = solution_next["vFunc"] # This is a NullFunc when vFuncBool is False
578 vPfuncNext = solution_next["vPfunc"]
579 vPPfuncNext = solution_next["vPPfunc"] # This is a NullFunc when CubicBool is False
581 # Update the bounding MPCs and PDV of human wealth:
582 PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
583 try:
584 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
585 except:
586 MPCminNow = 0.0
587 mLvlMinNext = solution_next["mLvlMin"]
589 # TODO: Deal with this unused code for the upper bound of MPC (should be a function now)
590 # Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
591 # hNrmNow = 0.0
592 # temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
593 # MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
595 # Define some functions for calculating future expectations
596 def calc_pLvl_next(S, p):
597 return pLvlNextFunc(p) * S["PermShk"]
599 def calc_mLvl_next(S, a, p_next):
600 return Rfree * a + p_next * S["TranShk"]
602 def calc_hLvl(S, p):
603 pLvl_next = calc_pLvl_next(S, p)
604 hLvl = S["TranShk"] * pLvl_next + solution_next["hLvl"](pLvl_next)
605 return hLvl
607 def calc_v_next(S, a, p):
608 pLvl_next = calc_pLvl_next(S, p)
609 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
610 v_next = vFuncNext(mLvl_next, pLvl_next)
611 return v_next
613 def calc_vP_next(S, a, p):
614 pLvl_next = calc_pLvl_next(S, p)
615 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
616 vP_next = vPfuncNext(mLvl_next, pLvl_next)
617 return vP_next
619 def calc_vPP_next(S, a, p):
620 pLvl_next = calc_pLvl_next(S, p)
621 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
622 vPP_next = vPPfuncNext(mLvl_next, pLvl_next)
623 return vPP_next
625 # Construct human wealth level as a function of productivity pLvl
626 hLvlGrid = 1.0 / Rfree * expected(calc_hLvl, IncShkDstn, args=(pLvlGrid))
627 hLvlNow = LinearInterp(np.insert(pLvlGrid, 0, 0.0), np.insert(hLvlGrid, 0, 0.0))
629 # Make temporary grids of income shocks and next period income values
630 ShkCount = TranShkValsNext.size
631 pLvlCount = pLvlGrid.size
632 PermShkVals_temp = np.tile(
633 np.reshape(PermShkValsNext, (1, ShkCount)), (pLvlCount, 1)
634 )
635 TranShkVals_temp = np.tile(
636 np.reshape(TranShkValsNext, (1, ShkCount)), (pLvlCount, 1)
637 )
638 pLvlNext_temp = (
639 np.tile(
640 np.reshape(pLvlNextFunc(pLvlGrid), (pLvlCount, 1)),
641 (1, ShkCount),
642 )
643 * PermShkVals_temp
644 )
646 # Find the natural borrowing constraint for each persistent income level
647 aLvlMin_candidates = (
648 mLvlMinNext(pLvlNext_temp) - TranShkVals_temp * pLvlNext_temp
649 ) / Rfree
650 aLvlMinNow = np.max(aLvlMin_candidates, axis=1)
651 BoroCnstNat = LinearInterp(
652 np.insert(pLvlGrid, 0, 0.0), np.insert(aLvlMinNow, 0, 0.0)
653 )
655 # Define the minimum allowable mLvl by pLvl as the greater of the natural and artificial borrowing constraints
656 if BoroCnstArt is not None:
657 BoroCnstArt = LinearInterp(np.array([0.0, 1.0]), np.array([0.0, BoroCnstArt]))
658 mLvlMinNow = UpperEnvelope(BoroCnstArt, BoroCnstNat)
659 else:
660 mLvlMinNow = BoroCnstNat
662 # Make the constrained total spending function: spend all market resources
663 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
664 SpendAllFunc = TrilinearInterp(
665 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
666 trivial_grid,
667 trivial_grid,
668 trivial_grid,
669 )
670 xFuncNowCnst = VariableLowerBoundFunc3D(SpendAllFunc, mLvlMinNow)
672 # Define grids of pLvl and aLvl on which to compute future expectations
673 pLvlCount = pLvlGrid.size
674 aNrmCount = aXtraGrid.size
675 MedCount = MedShkVals.size
676 pLvlNow = np.tile(np.reshape(pLvlGrid, (1, pLvlCount)), (aNrmCount, 1))
677 aLvlNow = np.reshape(aXtraGrid, (aNrmCount, 1)) * pLvlNow + BoroCnstNat(pLvlNow)
678 if pLvlGrid[0] == 0.0: # aLvl turns out badly if pLvl is 0 at bottom
679 aLvlNow[:, 0] = aXtraGrid
681 # Calculate end-of-period marginal value of assets
682 EndOfPrd_vP = (
683 DiscFacEff * Rfree * expected(calc_vP_next, IncShkDstn, args=(aLvlNow, pLvlNow))
684 )
686 # If the value function has been requested, construct the end-of-period vFunc
687 if vFuncBool:
688 # Compute expected value from end-of-period states
689 EndOfPrd_v = expected(calc_v_next, IncShkDstn, args=(aLvlNow, pLvlNow))
690 EndOfPrd_v *= DiscFacEff
692 # Transformed value through inverse utility function to "decurve" it
693 EndOfPrd_vNvrs = uCon.inv(EndOfPrd_v)
695 # Add points at mLvl=zero
696 EndOfPrd_vNvrs = np.concatenate(
697 (np.zeros((1, pLvlCount)), EndOfPrd_vNvrs), axis=0
698 )
700 # Make a temporary aLvl grid for interpolating the end-of-period value function
701 aLvl_temp = np.concatenate(
702 (
703 np.reshape(BoroCnstNat(pLvlGrid), (1, pLvlGrid.size)),
704 aLvlNow,
705 ),
706 axis=0,
707 )
709 # Make an end-of-period value function for each persistent income level in the grid
710 EndOfPrd_vNvrsFunc_list = []
711 for p in range(pLvlCount):
712 EndOfPrd_vNvrsFunc_list.append(
713 LinearInterp(
714 aLvl_temp[:, p] - BoroCnstNat(pLvlGrid[p]),
715 EndOfPrd_vNvrs[:, p],
716 )
717 )
718 EndOfPrd_vNvrsFuncBase = LinearInterpOnInterp1D(
719 EndOfPrd_vNvrsFunc_list, pLvlGrid
720 )
722 # Re-adjust the combined end-of-period value function to account for the
723 # natural borrowing constraint shifter and "re-curve" it
724 EndOfPrd_vNvrsFunc = VariableLowerBoundFunc2D(
725 EndOfPrd_vNvrsFuncBase, BoroCnstNat
726 )
727 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
729 # Solve the first order condition to get optimal consumption and medical
730 # spending, then find the endogenous mLvl gridpoints
731 cLvlNow = np.tile(
732 np.reshape(uCon.derinv(EndOfPrd_vP, order=(1, 0)), (aNrmCount, pLvlCount, 1)),
733 (1, 1, MedCount),
734 )
735 MedShkVals_adj = np.reshape(MedShkVals ** (1.0 - 1.0 / CRRAmed), (1, 1, MedCount))
736 MedBase = MedPrice ** (-1.0 / CRRAmed) * uMed.derinv(EndOfPrd_vP, order=(1, 0))
737 MedLvlNow = np.reshape(MedBase, (aNrmCount, pLvlCount, 1)) * MedShkVals_adj
738 MedLvlNow = np.maximum(MedLvlNow - MedShift, 0.0)
739 aLvlNow_tiled = np.tile(
740 np.reshape(aLvlNow, (aNrmCount, pLvlCount, 1)), (1, 1, MedCount)
741 )
742 xLvlNow = cLvlNow + MedPrice * MedLvlNow
743 mLvlNow = xLvlNow + aLvlNow_tiled
745 # Limiting consumption is zero as m approaches the natural borrowing constraint
746 x_for_interpolation = np.concatenate(
747 (np.zeros((1, pLvlCount, MedCount)), xLvlNow), axis=0
748 )
749 temp = np.tile(
750 BoroCnstNat(np.reshape(pLvlGrid, (1, pLvlCount, 1))),
751 (1, 1, MedCount),
752 )
753 m_for_interpolation = np.concatenate((temp, mLvlNow), axis=0)
755 # Make a 3D array of permanent income for interpolation
756 p_for_interpolation = np.tile(
757 np.reshape(pLvlGrid, (1, pLvlCount, 1)), (aNrmCount + 1, 1, MedCount)
758 )
759 MedShkVals_tiled = np.tile( # This does *not* have the CRRA adjustment
760 np.reshape(MedShkVals, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
761 )
763 # Build the set of cFuncs by pLvl, gathered in a list
764 xFunc_by_pLvl_and_MedShk = [] # Initialize the empty list of lists of 1D xFuncs
765 if CubicBool:
766 raise NotImplementedError(
767 "Cubic spline interpolation has not been implemented yet!"
768 )
770 # Basic version: use linear interpolation within a pLvl and MedShk
771 else:
772 # Loop over pLvl and then MedShk within that
773 for i in range(pLvlCount):
774 temp_list = []
775 pLvl_i = p_for_interpolation[0, i, 0]
776 mLvlMin_i = BoroCnstNat(pLvl_i)
777 for j in range(MedCount):
778 m_temp = m_for_interpolation[:, i, j] - mLvlMin_i
779 x_temp = x_for_interpolation[:, i, j]
780 temp_list.append(LinearInterp(m_temp, x_temp))
781 xFunc_by_pLvl_and_MedShk.append(deepcopy(temp_list))
783 # Combine the nested list of linear xFuncs into a single function
784 pLvl_temp = p_for_interpolation[0, :, 0]
785 MedShk_temp = MedShkVals_tiled[0, 0, :]
786 xFuncUncBase = BilinearInterpOnInterp1D(
787 xFunc_by_pLvl_and_MedShk, pLvl_temp, MedShk_temp
788 )
789 xFuncNowUnc = VariableLowerBoundFunc3D(xFuncUncBase, BoroCnstNat)
790 # Re-adjust for lower bound of natural borrowing constraint
792 # Combine the constrained and unconstrained functions into the true consumption function
793 xFuncNow = LowerEnvelope3D(xFuncNowUnc, xFuncNowCnst)
795 # Build the policy function for this period as a single object
796 PolicyFuncNow = cAndMedFunc(xFuncNow, qFunc, MedShift, MedPrice)
798 # Make the marginal value function by integrating over medical shocks
799 # Make temporary grids to evaluate the consumption function
800 temp_grid = np.tile(
801 np.reshape(aXtraGrid, (aNrmCount, 1, 1)), (1, pLvlCount, MedCount)
802 )
803 aMinGrid = np.tile(
804 np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount, 1)),
805 (aNrmCount, 1, MedCount),
806 )
807 pGrid = np.tile(np.reshape(pLvlGrid, (1, pLvlCount, 1)), (aNrmCount, 1, MedCount))
808 mGrid = temp_grid * pGrid + aMinGrid
809 if pLvlGrid[0] == 0:
810 mGrid[:, 0, :] = np.tile(np.reshape(aXtraGrid, (aNrmCount, 1)), (1, MedCount))
811 MedShkGrid = np.tile(
812 np.reshape(MedShkVals, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
813 )
814 ProbsGrid = np.reshape(MedShkPrbs, (1, 1, MedCount))
816 # Get optimal consumption (and medical care) for each state
817 cGrid, MedGrid, xTrash = PolicyFuncNow(mGrid, pGrid, MedShkGrid)
819 # Calculate expected marginal value by "integrating" across medical shocks
820 vPgrid = uCon.der(cGrid)
821 vPnow = np.sum(vPgrid * ProbsGrid, axis=2)
823 # Add vPnvrs=0 at m=mLvlMin to close it off at the bottom (and vNvrs=0)
824 mGrid_small = np.concatenate(
825 (np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount)), mGrid[:, :, 0])
826 )
827 vPnvrsNow = np.concatenate(
828 (np.zeros((1, pLvlCount)), uCon.derinv(vPnow, order=(1, 0)))
829 )
831 # Calculate expected value by "integrating" across medical shocks
832 if vFuncBool:
833 # interpolation error sometimes makes Med < 0 (barely), so fix that
834 MedGrid = np.maximum(MedGrid, 1e-100)
835 # interpolation error sometimes makes tiny violations, so fix that
836 aGrid = np.maximum(mGrid - cGrid - MedPrice * MedGrid, aMinGrid)
837 MedEff = (MedGrid + MedShift) / MedShkGrid
838 vGrid = uCon(cGrid) + uMed(MedEff) + EndOfPrd_vFunc(aGrid, pGrid)
839 vNow = np.sum(vGrid * ProbsGrid, axis=2)
841 # Switch to pseudo-inverse value and add a point at bottom
842 vNvrsNow = np.concatenate((np.zeros((1, pLvlCount)), uCon.inv(vNow)), axis=0)
844 # Construct the pseudo-inverse value and marginal value functions over mLvl,pLvl
845 vPnvrsFunc_by_pLvl = []
846 vNvrsFunc_by_pLvl = []
847 # Make a pseudo inverse marginal value function for each pLvl
848 for j in range(pLvlCount):
849 pLvl = pLvlGrid[j]
850 m_temp = mGrid_small[:, j] - mLvlMinNow(pLvl)
851 vPnvrs_temp = vPnvrsNow[:, j]
852 vPnvrsFunc_by_pLvl.append(LinearInterp(m_temp, vPnvrs_temp))
853 if vFuncBool:
854 vNvrs_temp = vNvrsNow[:, j]
855 vNvrsFunc_by_pLvl.append(LinearInterp(m_temp, vNvrs_temp))
857 # Combine those functions across pLvls, and adjust for the lower bound of mLvl
858 vPnvrsFuncBase = LinearInterpOnInterp1D(vPnvrsFunc_by_pLvl, pLvlGrid)
859 vPnvrsFunc = VariableLowerBoundFunc2D(vPnvrsFuncBase, mLvlMinNow)
860 if vFuncBool:
861 vNvrsFuncBase = LinearInterpOnInterp1D(vNvrsFunc_by_pLvl, pLvlGrid)
862 vNvrsFunc = VariableLowerBoundFunc2D(vNvrsFuncBase, mLvlMinNow)
864 # "Re-curve" the (marginal) value function
865 vPfuncNow = MargValueFuncCRRA(vPnvrsFunc, CRRA)
866 if vFuncBool:
867 vFuncNow = ValueFuncCRRA(vNvrsFunc, CRRA)
868 else:
869 vFuncNow = NullFunc()
871 # If using cubic spline interpolation, construct the marginal marginal value function
872 if CubicBool:
873 vPPfuncNow = MargMargValueFuncCRRA(vPfuncNow.cFunc, CRRA)
874 else:
875 vPPfuncNow = NullFunc()
877 # Package and return the solution as a dictionary
878 solution_now = {
879 "distance_criteria": ["vPfunc"],
880 "PolicyFunc": PolicyFuncNow,
881 "vFunc": vFuncNow,
882 "vPfunc": vPfuncNow,
883 "vPPfunc": vPPfuncNow,
884 "hLvl": hLvlNow,
885 "mLvlMin": mLvlMinNow,
886 }
887 return solution_now
890###############################################################################
892# Make a constructor dictionary for the general income process consumer type
893medshock_constructor_dict = {
894 "IncShkDstn": construct_lognormal_income_process_unemployment,
895 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
896 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
897 "aXtraGrid": make_assets_grid,
898 "pLvlPctiles": make_basic_pLvlPctiles,
899 "pLvlGrid": make_pLvlGrid_by_simulation,
900 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
901 "MedShkDstn": make_lognormal_MedShkDstn,
902 "qFunc": make_qFunc,
903 "solution_terminal": make_MedShock_solution_terminal,
904 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
905 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
906}
908# Make a dictionary with parameters for the default constructor for kNrmInitDstn
909default_kNrmInitDstn_params = {
910 "kLogInitMean": -6.0, # Mean of log initial capital
911 "kLogInitStd": 1.0, # Stdev of log initial capital
912 "kNrmInitCount": 15, # Number of points in initial capital discretization
913}
915# Make a dictionary with parameters for the default constructor for pLvlInitDstn
916default_pLvlInitDstn_params = {
917 "pLogInitMean": 0.0, # Mean of log permanent income
918 "pLogInitStd": 0.4, # Stdev of log permanent income
919 "pLvlInitCount": 15, # Number of points in initial capital discretization
920}
922# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
923default_IncShkDstn_params = {
924 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
925 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
926 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
927 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
928 "UnempPrb": 0.05, # Probability of unemployment while working
929 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
930 "T_retire": 0, # Period of retirement (0 --> no retirement)
931 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
932 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
933}
935# Default parameters to make aXtraGrid using make_assets_grid
936default_aXtraGrid_params = {
937 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
938 "aXtraMax": 30, # Maximum end-of-period "assets above minimum" value
939 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
940 "aXtraCount": 32, # Number of points in the grid of "assets above minimum"
941 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
942}
944# Default parameters to make pLvlGrid using make_basic_pLvlPctiles
945default_pLvlPctiles_params = {
946 "pLvlPctiles_count": 19, # Number of points in the "body" of the grid
947 "pLvlPctiles_bound": [0.05, 0.95], # Percentile bounds of the "body"
948 "pLvlPctiles_tail_count": 4, # Number of points in each tail of the grid
949 "pLvlPctiles_tail_order": np.e, # Scaling factor for points in each tail
950}
952# Default parameters to make pLvlGrid using make_trivial_pLvlNextFunc
953default_pLvlGrid_params = {
954 "pLvlInitMean": 0.0, # Mean of log initial permanent income
955 "pLvlInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
956 "pLvlExtra": [0.0001],
957 # Additional permanent income points to automatically add to the grid, optional
958}
960# Default parameters to make pLvlNextFunc using make_AR1_style_pLvlNextFunc
961default_pLvlNextFunc_params = {
962 "PermGroFac": [1.0], # Permanent income growth factor
963 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
964}
966# Default parameters to make MedShkDstn using make_lognormal_MedShkDstn
967default_MedShkDstn_params = {
968 "MedShkAvg": [0.1], # Average of medical need shocks
969 "MedShkStd": [1.5], # Standard deviation of (log) medical need shocks
970 "MedShkCount": 5, # Number of medical shock points in "body"
971 "MedShkCountTail": 15, # Number of medical shock points in "tail" (upper only)
972 "MedPrice": [1.5], # Relative price of a unit of medical care
973}
975# Make a dictionary to specify a medical shocks consumer type
976init_medical_shocks = {
977 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
978 "cycles": 1, # Finite, non-cyclic model
979 "T_cycle": 1, # Number of periods in the cycle for this agent type
980 "pseudo_terminal": False, # Terminal period really does exist
981 "constructors": medshock_constructor_dict, # See dictionary above
982 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
983 "CRRA": 2.0, # Coefficient of relative risk aversion on consumption
984 "CRRAmed": 5.0, # Coefficient of relative risk aversion on medical care
985 "MedShift": [1e-8], # Amount of medical care the agent can self-provide
986 "Rfree": [1.03], # Interest factor on retained assets
987 "DiscFac": 0.96, # Intertemporal discount factor
988 "LivPrb": [0.99], # Survival probability after each period
989 "BoroCnstArt": 0.0, # Artificial borrowing constraint
990 "vFuncBool": False, # Whether to calculate the value function during solution
991 "CubicBool": False, # Whether to use cubic spline interpolation when True
992 # (Uses linear spline interpolation for cFunc when False)
993 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
994 "AgentCount": 10000, # Number of agents of this type
995 "T_age": None, # Age after which simulated agents are automatically killed
996 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
997 # (The portion of PermGroFac attributable to aggregate productivity growth)
998 "NewbornTransShk": False, # Whether Newborns have transitory shock
999 # ADDITIONAL OPTIONAL PARAMETERS
1000 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1001 # (Forces Newborns to follow solution path of the agent they replaced if True)
1002 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1003}
1004init_medical_shocks.update(default_IncShkDstn_params)
1005init_medical_shocks.update(default_aXtraGrid_params)
1006init_medical_shocks.update(default_pLvlPctiles_params)
1007init_medical_shocks.update(default_pLvlGrid_params)
1008init_medical_shocks.update(default_MedShkDstn_params)
1009init_medical_shocks.update(default_pLvlNextFunc_params)
1010init_medical_shocks.update(default_pLvlInitDstn_params)
1011init_medical_shocks.update(default_kNrmInitDstn_params)
1014# This dictionary is based on the main specification results in Fulford and Low's
1015# "Expense Shocks Matter". These expenses represent *all* unexpected spending, not
1016# just medical expenses. It is calibrated at an annual frequency. The specification
1017# in their paper has serially correlated expense shocks (with a low correlation
1018# coefficient of about 0.086) and serially correlated unemployment ("crisis income"),
1019# which are not present for MedShockConsumerType. The unemployment probability here
1020# is thus the (approximate) fraction of the time a consumer will spend in the crisis
1021# state. Moreover, the medical shock parameters here differ from those in Fulford &
1022# Low's January 2026 paper version. In private correspondence with them, MNW found
1023# that their results hinged critically on the specific method they used to discretize
1024# the lognormal shock distribution. The parameters below match the following three
1025# key statistics from their exercise: mean of (non-zero) expense ratio is about 16.2%,
1026# standard deviation of log non-zero expense ratio is about 1.1, and mean wealth is
1027# about 0.6 times income.
1028Fulford_and_Low_params = {
1029 "cycles": 0,
1030 "DiscFac": 0.85,
1031 "LivPrb": [1.0],
1032 "CRRA": 2.0,
1033 "CRRAmed": 4.0,
1034 "Rfree": [1.01],
1035 "TranShkStd": [0.2],
1036 "PermShkStd": [0.117],
1037 "PrstIncCorr": 0.97,
1038 "BoroCnstArt": -0.185,
1039 "MedShkAvg": [0.121],
1040 "MedShkStd": [1.55],
1041 "MedShkCountTail": 5,
1042 "MedShkTailBound": [0.0, 0.9],
1043 "MedShkZeroPrb": [0.31],
1044 "MedPrice": [1.0],
1045 "IncUnemp": 0.195,
1046 "UnempPrb": 0.018,
1047}
1050class MedShockConsumerType(PersistentShockConsumerType):
1051 r"""
1052 A consumer type based on GenIncShockConsumerType, with two types of consumption goods (medical and nonmedical) and random shocks to medical utility.
1054 .. math::
1055 \begin{eqnarray*}
1056 V_t(M_t,P_t,\eta_t) &=& \max_{C_t, med_t} U_t(C_t, med_t; \eta_t) + \beta (1-\mathsf{D}_{t+1}) \mathbb{E} [V_{t+1}(M_{t+1}, P_{t+1}, \text{medShk}_{t+1})], \\
1057 A_t &=& M_t - X_t, \\
1058 X_t &=& C_t +med_t \textbf{ medPrice}_t,\\
1059 A_t/P_t &\geq& \underline{a}, \\
1060 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
1061 Y_{t+1} &=& P_{t+1} \theta_{t+1}
1062 M_{t+1} &=& R A_t + Y_{t+1}, \\
1063 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
1064 \eta_t &~\sim& G_t,\\
1065 U_t(C, med; \eta) &=& \frac{C^{1-\rho}}{1-\rho} +\eta \frac{med^{1-\nu}}{1-\nu}.
1066 \end{eqnarray*}
1069 Constructors
1070 ------------
1071 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
1072 The agent's income shock distributions.
1073 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
1074 aXtraGrid: Constructor
1075 The agent's asset grid.
1076 Its default constructor is :func:`HARK.utilities.make_assets_grid`
1077 pLvlNextFunc: Constructor
1078 An arbitrary function used to evolve the GenIncShockConsumerType's permanent income
1079 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_trivial_pLvlNextFunc`
1080 pLvlGrid: Constructor
1081 The agent's pLvl grid
1082 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_pLvlGrid_by_simulation`
1083 pLvlPctiles: Constructor
1084 The agents income level percentile grid
1085 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_basic_pLvlPctiles`
1086 MedShkDstn: Constructor, :math:`\text{medShk}`
1087 The agent's Medical utility shock distribution.
1088 Its default constructor is :func:`HARK.ConsumptionSaving.ConsMedModel.make_lognormal_MedShkDstn`
1090 Solving Parameters
1091 ------------------
1092 cycles: int
1093 0 specifies an infinite horizon model, 1 specifies a finite model.
1094 T_cycle: int
1095 Number of periods in the cycle for this agent type.
1096 CRRA: float, :math:`\rho`
1097 Coefficient of Relative Risk Aversion for consumption.
1098 CRRAmed: float, :math:`\nu`
1099 Coefficient of Relative Risk Aversion for medical care.
1100 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1101 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1102 DiscFac: float, :math:`\beta`
1103 Intertemporal discount factor.
1104 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1105 Survival probability after each period.
1106 PermGroFac: list[float], time varying, :math:`\Gamma`
1107 Permanent income growth factor.
1108 BoroCnstArt: float, :math:`\underline{a}`
1109 The minimum Asset/Perminant Income ratio, None to ignore.
1110 vFuncBool: bool
1111 Whether to calculate the value function during solution.
1112 CubicBool: bool
1113 Whether to use cubic spline interpoliation.
1115 Simulation Parameters
1116 ---------------------
1117 AgentCount: int
1118 Number of agents of this kind that are created during simulations.
1119 T_age: int
1120 Age after which to automatically kill agents, None to ignore.
1121 T_sim: int, required for simulation
1122 Number of periods to simulate.
1123 track_vars: list[strings]
1124 List of variables that should be tracked when running the simulation.
1125 For this agent, the options are 'Med', 'MedShk', 'PermShk', 'TranShk', 'aLvl', 'cLvl', 'mLvl', 'pLvl', and 'who_dies'.
1127 PermShk is the agent's permanent income shock
1129 MedShk is the agent's medical utility shock
1131 TranShk is the agent's transitory income shock
1133 aLvl is the nominal asset level
1135 cLvl is the nominal consumption level
1137 MedLvl is the nominal medical spending level
1139 mLvl is the nominal market resources
1141 pLvl is the permanent income level
1143 who_dies is the array of which agents died
1144 aNrmInitMean: float
1145 Mean of Log initial Normalized Assets.
1146 aNrmInitStd: float
1147 Std of Log initial Normalized Assets.
1148 pLvlInitMean: float
1149 Mean of Log initial permanent income.
1150 pLvlInitStd: float
1151 Std of Log initial permanent income.
1152 PermGroFacAgg: float
1153 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1154 PerfMITShk: boolean
1155 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1156 NewbornTransShk: boolean
1157 Whether Newborns have transitory shock.
1159 Attributes
1160 ----------
1161 solution: list[Consumer solution object]
1162 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1163 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1165 Unlike other models with this solution type, this model's variables are NOT normalized.
1166 The solution functions additionally depend on the permanent income level and the medical shock.
1167 For example, :math:`C=\text{cFunc}(M,P,MedShk)`.
1168 hNrm has been replaced by hLvl which is a function of permanent income.
1169 MPC max has not yet been implemented for this class. It will be a function of permanent income.
1171 This solution has two additional functions
1172 :math:`\text{Med}=\text{MedFunc}(M,P,\text{MedShk})`: returns the agent's spending on Medical care
1174 :math:`[C,Med]=\text{policyFunc}(M,P,\text{MedShk})`: returns the agent's spending on consumption and Medical care as numpy arrays
1176 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1177 history: Dict[Array]
1178 Created by running the :func:`.simulate()` method.
1179 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1180 Visit :class:`HARK.core.AgentType.simulate` for more information.
1181 """
1183 default_ = {
1184 "params": init_medical_shocks,
1185 "solver": solve_one_period_ConsMedShock,
1186 "model": "ConsMedShock.yaml",
1187 "track_vars": ["aLvl", "cLvl", "MedLvl", "mLvl", "pLvl"],
1188 }
1190 time_vary_ = PersistentShockConsumerType.time_vary_ + [
1191 "MedPrice",
1192 "MedShkDstn",
1193 "MedShift",
1194 ]
1195 time_inv_ = PersistentShockConsumerType.time_inv_ + ["CRRAmed", "qFunc"]
1196 shock_vars_ = PersistentShockConsumerType.shock_vars_ + ["MedShk"]
1197 state_vars = PersistentShockConsumerType.state_vars + ["mLvl"]
1198 distributions = [
1199 "IncShkDstn",
1200 "PermShkDstn",
1201 "TranShkDstn",
1202 "kNrmInitDstn",
1203 "pLvlInitDstn",
1204 "MedShkDstn",
1205 ]
1207 def pre_solve(self):
1208 self.construct("solution_terminal")
1210 def get_shocks(self):
1211 """
1212 Gets permanent and transitory income shocks for this period as well as medical need shocks
1213 and the price of medical care.
1215 Parameters
1216 ----------
1217 None
1219 Returns
1220 -------
1221 None
1222 """
1223 # Get permanent and transitory income shocks
1224 PersistentShockConsumerType.get_shocks(self)
1226 # Initialize medical shock array and relative price array
1227 MedShkNow = np.zeros(self.AgentCount)
1228 MedPriceNow = np.zeros(self.AgentCount)
1230 # Get shocks for each period of the cycle
1231 for t in range(self.T_cycle):
1232 these = t == self.t_cycle
1233 N = np.sum(these)
1234 if N > 0:
1235 MedShkNow[these] = self.MedShkDstn[t].draw(N)
1236 MedPriceNow[these] = self.MedPrice[t]
1237 self.shocks["MedShk"] = MedShkNow
1238 self.shocks["MedPrice"] = MedPriceNow
1240 def get_controls(self):
1241 """
1242 Calculates consumption and medical care for each consumer of this type
1243 using the consumption and medical care functions.
1245 Parameters
1246 ----------
1247 None
1249 Returns
1250 -------
1251 None
1252 """
1253 cLvlNow = np.empty(self.AgentCount)
1254 MedNow = np.empty(self.AgentCount)
1255 for t in range(self.T_cycle):
1256 these = t == self.t_cycle
1257 cLvlNow[these], MedNow[these], x = self.solution[t]["PolicyFunc"](
1258 self.state_now["mLvl"][these],
1259 self.state_now["pLvl"][these],
1260 self.shocks["MedShk"][these],
1261 )
1262 self.controls["cLvl"] = cLvlNow
1263 self.controls["MedLvl"] = MedNow
1265 def get_poststates(self):
1266 """
1267 Calculates end-of-period assets for each consumer of this type.
1269 Parameters
1270 ----------
1271 None
1273 Returns
1274 -------
1275 None
1276 """
1277 self.state_now["aLvl"] = (
1278 self.state_now["mLvl"]
1279 - self.controls["cLvl"]
1280 - self.shocks["MedPrice"] * self.controls["MedLvl"]
1281 )
1284###############################################################################
1287class ConsMedExtMargSolution(MetricObject):
1288 """
1289 Representation of the solution to one period's problem in the extensive margin
1290 medical expense model. If no inputs are passed, a trivial object is constructed,
1291 which can be used as the pseudo-terminal solution.
1293 Parameters
1294 ----------
1295 vFunc_by_pLvl : [function]
1296 List of beginning-of-period value functions over kLvl, by pLvl.
1297 vPfunc_by_pLvl : [function]
1298 List of beginning-of-period marginal functions over kLvl, by pLvl.
1299 cFunc_by_pLvl : [function]
1300 List of consumption functions over bLvl, by pLvl.
1301 vNvrsFuncMid_by_pLvl : [function]
1302 List of pseudo-inverse value function for consumption phase over bLvl, by pLvl.
1303 ExpMedFunc : function
1304 Expected medical care as a function of mLvl and pLvl, just before medical
1305 shock is realized.
1306 CareProbFunc : function
1307 Probability of getting medical treatment as a function of mLvl and pLvl,
1308 just before medical shock is realized.
1309 pLvl : np.array
1310 Grid of permanent income levels during the period (after shocks).
1311 CRRA : float
1312 Coefficient of relative risk aversion
1313 """
1315 distance_criteria = ["cFunc"]
1317 def __init__(
1318 self,
1319 vFunc_by_pLvl=None,
1320 vPfunc_by_pLvl=None,
1321 cFunc_by_pLvl=None,
1322 vNvrsFuncMid_by_pLvl=None,
1323 ExpMedFunc=None,
1324 CareProbFunc=None,
1325 pLvl=None,
1326 CRRA=None,
1327 ):
1328 self.pLvl = pLvl
1329 self.CRRA = CRRA
1330 if vFunc_by_pLvl is None:
1331 self.vFunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1332 else:
1333 self.vFunc_by_pLvl = vFunc_by_pLvl
1334 if vPfunc_by_pLvl is None:
1335 self.vPfunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1336 else:
1337 self.vPfunc_by_pLvl = vPfunc_by_pLvl
1338 if cFunc_by_pLvl is not None:
1339 self.cFunc = LinearInterpOnInterp1D(cFunc_by_pLvl, pLvl)
1340 else:
1341 self.cFunc = None
1342 if vNvrsFuncMid_by_pLvl is not None:
1343 vNvrsFuncMid = LinearInterpOnInterp1D(vNvrsFuncMid_by_pLvl, pLvl)
1344 self.vFuncMid = ValueFuncCRRA(vNvrsFuncMid, CRRA, illegal_value=-np.inf)
1345 if ExpMedFunc is not None:
1346 self.ExpMedFunc = ExpMedFunc
1347 if CareProbFunc is not None:
1348 self.CareProbFunc = CareProbFunc
1351def make_MedExtMarg_solution_terminal(pLogCount):
1352 """
1353 Construct a trivial pseudo-terminal solution for the extensive margin medical
1354 spending model: a list of constant zero functions for (marginal) value. The
1355 only piece of information needed for this is how many such functions to include.
1356 """
1357 pLvl_terminal = np.arange(pLogCount)
1358 solution_terminal = ConsMedExtMargSolution(pLvl=pLvl_terminal)
1359 return solution_terminal
1362###############################################################################
1365def solve_one_period_ConsMedExtMarg(
1366 solution_next,
1367 DiscFac,
1368 CRRA,
1369 BeqFac,
1370 BeqShift,
1371 Rfree,
1372 LivPrb,
1373 MedShkLogMean,
1374 MedShkLogStd,
1375 MedCostLogMean,
1376 MedCostLogStd,
1377 MedCorr,
1378 MedCostBot,
1379 MedCostTop,
1380 MedCostCount,
1381 aNrmGrid,
1382 pLogGrid,
1383 pLvlMean,
1384 TranShkDstn,
1385 pLogMrkvArray,
1386 mNrmGrid,
1387 kLvlGrid,
1388):
1389 """
1390 Solve one period of the "extensive margin medical care" model. Each period, the
1391 agent receives a persistent and transitory shock to income, and then a medical
1392 shock with two components: utility and cost. He makes a binary choice between
1393 paying the cost in medical expenses or suffering the utility loss, then makes
1394 his ordinary consumption-saving decision (technically made simultaneously, but
1395 solved as if sequential). This version has one health state and no insurance choice
1396 and hardcodes a liquidity constraint.
1398 Parameters
1399 ----------
1400 solution_next : ConsMedExtMargSolution
1401 Solution to the succeeding period's problem.
1402 DiscFac : float
1403 Intertemporal discount factor.
1404 CRRA : float
1405 Coefficient of relative risk aversion.
1406 BeqFac : float
1407 Scaling factor for bequest motive.
1408 BeqShift : float
1409 Shifter for bequest motive.
1410 Rfree : float
1411 Risk free return factor on saving.
1412 LivPrb : float
1413 Survival probability from this period to the next one.
1414 MedShkLogMean : float
1415 Mean of log utility shocks, assumed to be lognormally distributed.
1416 MedShkLogStd : float
1417 Stdev of log utility shocks, assumed to be lognormally distributed.
1418 MedCostLogMean : float
1419 Mean of log medical expense shocks, assumed to be lognormally distributed.
1420 MedCostLogStd : float
1421 Stdev of log medical expense shocks, assumed to be lognormally distributed.
1422 MedCorr : float
1423 Correlation coefficient betwen log utility shocks and log medical expense
1424 shocks, assumed to be joint normal (in logs).
1425 MedCostBot : float
1426 Lower bound of medical costs to consider, as standard deviations of log
1427 expenses away from the mean.
1428 MedCostTop : float
1429 Upper bound of medical costs to consider, as standard deviations of log
1430 expenses away from the mean.
1431 MedCostCount : int
1432 Number of points to use when discretizing MedCost
1433 aNrmGrid : np.array
1434 Exogenous grid of end-of-period assets, normalized by income level.
1435 pLogGrid : np.array
1436 Exogenous grid of *deviations from mean* log income level.
1437 pLvlMean : float
1438 Mean income level at this age, for generating actual income levels from
1439 pLogGrid as pLvl = pLvlMean * np.exp(pLogGrid).
1440 TranShkDstn : DiscreteDistribution
1441 Discretized transitory income shock distribution.
1442 pLogMrkvArray : np.array
1443 Markov transition array from beginning-of-period (prior) income levels
1444 to this period's levels. Pre-computed by (e.g.) Tauchen's method.
1445 mNrmGrid : np.array
1446 Exogenous grid of decision-time normalized market resources,
1447 kLvlGrid : np.array
1448 Beginning-of-period capital grid (spanning zero to very high wealth).
1450 Returns
1451 -------
1452 solution_now : ConsMedExtMargSolution
1453 Representation of the solution to this period's problem, including the
1454 beginning-of-period (marginal) value function by pLvl, the consumption
1455 function by pLvl, and the (pseudo-inverse) value function for the consumption
1456 phase (as a list by pLvl).
1457 """
1458 # Define (marginal) utility and bequest motive functions
1459 u = lambda x: CRRAutility(x, rho=CRRA)
1460 uP = lambda x: CRRAutilityP(x, rho=CRRA)
1461 W = lambda x: BeqFac * u(x + BeqShift)
1462 Wp = lambda x: BeqFac * uP(x + BeqShift)
1463 n = lambda x: CRRAutility_inv(x, rho=CRRA)
1465 # Make grids of pLvl x aLvl
1466 pLvl = np.exp(pLogGrid) * pLvlMean
1467 aLvl = np.dot(
1468 np.reshape(aNrmGrid, (aNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvl.size))
1469 )
1470 aLvl = np.concatenate([np.zeros((1, pLvl.size)), aLvl]) # add zero entries
1472 # Evaluate end-of-period marginal value at each combination of pLvl x aLvl
1473 pLvlCount = pLvl.size
1474 EndOfPrd_vP = np.empty_like(aLvl)
1475 EndOfPrd_v = np.empty_like(aLvl)
1476 for j in range(pLvlCount):
1477 EndOfPrdvFunc_this_pLvl = solution_next.vFunc_by_pLvl[j]
1478 EndOfPrdvPfunc_this_pLvl = solution_next.vPfunc_by_pLvl[j]
1479 EndOfPrd_v[:, j] = DiscFac * LivPrb * EndOfPrdvFunc_this_pLvl(aLvl[:, j])
1480 EndOfPrd_vP[:, j] = DiscFac * LivPrb * EndOfPrdvPfunc_this_pLvl(aLvl[:, j])
1481 EndOfPrd_v += (1.0 - LivPrb) * W(aLvl)
1482 EndOfPrd_vP += (1.0 - LivPrb) * Wp(aLvl)
1484 # Calculate optimal consumption for each (aLvl,pLvl) gridpoint, roll back to bLvl
1485 cLvl = CRRAutilityP_inv(EndOfPrd_vP, CRRA)
1486 bLvl = aLvl + cLvl
1488 # Construct consumption functions over bLvl for each pLvl
1489 cFunc_by_pLvl = []
1490 for j in range(pLvlCount):
1491 cFunc_j = LinearInterp(
1492 np.insert(bLvl[:, j], 0, 0.0), np.insert(cLvl[:, j], 0, 0.0)
1493 )
1494 cFunc_by_pLvl.append(cFunc_j)
1496 # Construct pseudo-inverse value functions over bLvl for each pLvl
1497 v_mid = u(cLvl) + EndOfPrd_v # value of reaching consumption phase
1498 vNvrsFuncMid_by_pLvl = []
1499 for j in range(pLvlCount):
1500 b_cnst = np.linspace(0.001, 0.95, 10) * bLvl[0, j] # constrained wealth levels
1501 c_cnst = b_cnst
1502 v_cnst = u(c_cnst) + EndOfPrd_v[0, j]
1503 b_temp = np.concatenate([b_cnst, bLvl[:, j]])
1504 v_temp = np.concatenate([v_cnst, v_mid[:, j]])
1505 vNvrs_temp = n(v_temp)
1506 vNvrsFunc_j = LinearInterp(
1507 np.insert(b_temp, 0, 0.0), np.insert(vNvrs_temp, 0, 0.0)
1508 )
1509 vNvrsFuncMid_by_pLvl.append(vNvrsFunc_j)
1511 # Make a grid of (log) medical expenses (and probs), cross it with (mLvl,pLvl)
1512 MedCostBaseGrid = np.linspace(MedCostBot, MedCostTop, MedCostCount)
1513 MedCostLogGrid = MedCostLogMean + MedCostBaseGrid * MedCostLogStd
1514 MedCostGrid = np.exp(MedCostLogGrid)
1515 mLvl_base = np.dot(
1516 np.reshape(mNrmGrid, (mNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvlCount))
1517 )
1518 mLvl = np.reshape(mLvl_base, (mNrmGrid.size, pLvlCount, 1))
1519 bLvl_if_care = mLvl - np.reshape(MedCostGrid, (1, 1, MedCostCount))
1520 bLvl_if_not = mLvl_base
1522 # Calculate mean (log) utility shock for each MedCost gridpoint, and conditional stdev
1523 MedShkLog_cond_mean = MedShkLogMean + MedCorr * MedShkLogStd * MedCostBaseGrid
1524 MedShkLog_cond_mean = np.reshape(MedShkLog_cond_mean, (1, MedCostCount))
1525 MedShkLog_cond_std = np.sqrt(MedShkLogStd**2 * (1.0 - MedCorr**2))
1526 MedShk_cond_mean = np.exp(MedShkLog_cond_mean + 0.5 * MedShkLog_cond_std**2)
1528 # Initialize (marginal) value function arrays over (mLvl,pLvl,MedCost)
1529 v_at_Dcsn = np.empty_like(bLvl_if_care)
1530 vP_at_Dcsn = np.empty_like(bLvl_if_care)
1531 care_prob_array = np.empty_like(bLvl_if_care)
1532 for j in range(pLvlCount):
1533 # Evaluate value function for (bLvl,pLvl_j), including MedCost=0
1534 v_if_care = u(vNvrsFuncMid_by_pLvl[j](bLvl_if_care[:, j, :]))
1535 v_if_not = np.reshape(
1536 u(vNvrsFuncMid_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1537 )
1538 cant_pay = bLvl_if_care[:, j, :] <= 0.0
1539 v_if_care[cant_pay] = -np.inf
1541 # Find value difference at each gridpoint, convert to MedShk stdev; find prob of care
1542 v_diff = v_if_not - v_if_care
1543 log_v_diff = np.log(v_diff)
1544 crit_stdev = (log_v_diff - MedShkLog_cond_mean) / MedShkLog_cond_std
1545 prob_no_care = norm.cdf(crit_stdev)
1546 prob_get_care = 1.0 - prob_no_care
1547 care_prob_array[:, j, :] = prob_get_care
1549 # Calculate expected MedShk conditional on not getting medical care
1550 crit_z = crit_stdev - MedShkLog_cond_std
1551 MedShk_no_care_cond_mean = 0.5 * MedShk_cond_mean * erfc(crit_z) / prob_no_care
1553 # Compute expected (marginal) value over MedShk for each (mLvl,pLvl_j,MedCost)
1554 v_if_care[cant_pay] = 0.0
1555 v_at_Dcsn[:, j, :] = (
1556 prob_no_care * (v_if_not - MedShk_no_care_cond_mean)
1557 + prob_get_care * v_if_care
1558 )
1559 vP_if_care = uP(cFunc_by_pLvl[j](bLvl_if_care[:, j, :]))
1560 vP_if_not = np.reshape(
1561 uP(cFunc_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1562 )
1563 vP_if_care[cant_pay] = 0.0
1564 MedShk_rate_of_change = (
1565 norm.pdf(crit_stdev) * (vP_if_care - vP_if_not) * MedShk_no_care_cond_mean
1566 )
1567 vP_at_Dcsn[:, j, :] = (
1568 prob_no_care * vP_if_not
1569 + prob_get_care * vP_if_care
1570 + MedShk_rate_of_change
1571 )
1573 # Compute expected (marginal) value over MedCost for each (mLvl,pLvl)
1574 temp_grid = np.linspace(MedCostBot, MedCostTop, MedCostCount)
1575 MedCost_pmv = norm.pdf(temp_grid)
1576 MedCost_pmv /= np.sum(MedCost_pmv)
1577 MedCost_probs = np.reshape(MedCost_pmv, (1, 1, MedCostCount))
1578 v_before_shk = np.sum(v_at_Dcsn * MedCost_probs, axis=2)
1579 vP_before_shk = np.sum(vP_at_Dcsn * MedCost_probs, axis=2)
1580 vNvrs_before_shk = n(v_before_shk)
1581 vPnvrs_before_shk = CRRAutilityP_inv(vP_before_shk, CRRA)
1583 # Compute expected medical expenses at each state space point
1584 ExpCare_all = care_prob_array * np.reshape(MedCostGrid, (1, 1, MedCostCount))
1585 ExpCare = np.sum(ExpCare_all * MedCost_probs, axis=2)
1586 ProbCare = np.sum(care_prob_array * MedCost_probs, axis=2)
1587 ExpCareFunc_by_pLvl = []
1588 CareProbFunc_by_pLvl = []
1589 for j in range(pLvlCount):
1590 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1591 EC_temp = np.insert(ExpCare[:, j], 0, 0.0)
1592 prob_temp = np.insert(ProbCare[:, j], 0, 0.0)
1593 ExpCareFunc_by_pLvl.append(LinearInterp(m_temp, EC_temp))
1594 CareProbFunc_by_pLvl.append(LinearInterp(m_temp, prob_temp))
1595 ExpCareFunc = LinearInterpOnInterp1D(ExpCareFunc_by_pLvl, pLvl)
1596 ProbCareFunc = LinearInterpOnInterp1D(CareProbFunc_by_pLvl, pLvl)
1598 # Fixing kLvlGrid, compute expected (marginal) value over TranShk for each (kLvl,pLvl)
1599 v_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1600 vP_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1601 for j in range(pLvlCount):
1602 p = pLvl[j]
1604 # Make (marginal) value functions over mLvl for this pLvl
1605 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1606 vNvrs_temp = np.insert(vNvrs_before_shk[:, j], 0, 0.0)
1607 vPnvrs_temp = np.insert(vPnvrs_before_shk[:, j], 0, 0.0)
1608 vNvrsFunc_temp = LinearInterp(m_temp, vNvrs_temp)
1609 vPnvrsFunc_temp = LinearInterp(m_temp, vPnvrs_temp)
1610 vFunc_temp = lambda x: u(vNvrsFunc_temp(x))
1611 vPfunc_temp = lambda x: uP(vPnvrsFunc_temp(x))
1613 # Compute expectation over TranShkDstn
1614 v = lambda TranShk, kLvl: vFunc_temp(kLvl + TranShk * p)
1615 vP = lambda TranShk, kLvl: vPfunc_temp(kLvl + TranShk * p)
1616 v_by_kLvl_and_pLvl[:, j] = expected(v, TranShkDstn, args=(kLvlGrid,))
1617 vP_by_kLvl_and_pLvl[:, j] = expected(vP, TranShkDstn, args=(kLvlGrid,))
1619 # Compute expectation over persistent shocks by using pLvlMrkvArray
1620 v_arvl = np.dot(v_by_kLvl_and_pLvl, pLogMrkvArray.T)
1621 vP_arvl = np.dot(vP_by_kLvl_and_pLvl, pLogMrkvArray.T)
1622 vNvrs_arvl = n(v_arvl)
1623 vPnvrs_arvl = CRRAutilityP_inv(vP_arvl, CRRA)
1625 # Construct "arrival" (marginal) value function by pLvl
1626 vFuncArvl_by_pLvl = []
1627 vPfuncArvl_by_pLvl = []
1628 for j in range(pLvlCount):
1629 vNvrsFunc_temp = LinearInterp(kLvlGrid, vNvrs_arvl[:, j])
1630 vPnvrsFunc_temp = LinearInterp(kLvlGrid, vPnvrs_arvl[:, j])
1631 vFuncArvl_by_pLvl.append(ValueFuncCRRA(vNvrsFunc_temp, CRRA))
1632 vPfuncArvl_by_pLvl.append(MargValueFuncCRRA(vPnvrsFunc_temp, CRRA))
1634 # Gather elements and return the solution object
1635 solution_now = ConsMedExtMargSolution(
1636 vFunc_by_pLvl=vFuncArvl_by_pLvl,
1637 vPfunc_by_pLvl=vPfuncArvl_by_pLvl,
1638 cFunc_by_pLvl=cFunc_by_pLvl,
1639 vNvrsFuncMid_by_pLvl=vNvrsFuncMid_by_pLvl,
1640 pLvl=pLvl,
1641 CRRA=CRRA,
1642 ExpMedFunc=ExpCareFunc,
1643 CareProbFunc=ProbCareFunc,
1644 )
1645 return solution_now
1648###############################################################################
1651# Define a dictionary of constructors for the extensive margin medical spending model
1652med_ext_marg_constructors = {
1653 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
1654 "IncomeProcessDict": make_persistent_income_process_dict,
1655 "pLogGrid": get_it_from("IncomeProcessDict"),
1656 "pLvlMean": get_it_from("IncomeProcessDict"),
1657 "pLogMrkvArray": get_it_from("IncomeProcessDict"),
1658 "IncShkDstn": construct_lognormal_income_process_unemployment,
1659 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1660 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1661 "BeqFac": get_it_from("BeqParamDict"),
1662 "BeqShift": get_it_from("BeqParamDict"),
1663 "BeqParamDict": reformat_bequest_motive,
1664 "aNrmGrid": make_assets_grid,
1665 "mNrmGrid": make_market_resources_grid,
1666 "kLvlGrid": make_capital_grid,
1667 "solution_terminal": make_MedExtMarg_solution_terminal,
1668 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1669 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1670 "MedShockDstn": make_continuous_MedShockDstn,
1671}
1673# Make a dictionary with default bequest motive parameters
1674default_BeqParam_dict = {
1675 "BeqMPC": 0.1, # Hypothetical "MPC at death"
1676 "BeqInt": 1.0, # Intercept term for hypothetical "consumption function at death"
1677}
1679# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1680default_kNrmInitDstn_params_ExtMarg = {
1681 "kLogInitMean": -6.0, # Mean of log initial capital
1682 "kLogInitStd": 1.0, # Stdev of log initial capital
1683 "kNrmInitCount": 15, # Number of points in initial capital discretization
1684}
1686# Default parameters to make IncomeProcessDict using make_persistent_income_process_dict;
1687# some of these are used by construct_lognormal_income_process_unemployment as well
1688default_IncomeProcess_params = {
1689 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1690 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1691 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1692 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1693 "UnempPrb": 0.05, # Probability of unemployment while working
1694 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1695 "T_retire": 0, # Period of retirement (0 --> no retirement)
1696 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1697 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1698 "pLogInitMean": 0.0, # Mean of log initial permanent income
1699 "pLogInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
1700 "pLvlInitCount": 25, # Number of discrete nodes in initial permanent income level dstn
1701 "PermGroFac": [1.0], # Permanent income growth factor
1702 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
1703 "pLogCount": 45, # Number of points in persistent income grid each period
1704 "pLogRange": 3.5, # Upper/lower bound of persistent income, in unconditional standard deviations
1705}
1707# Default parameters to make aNrmGrid using make_assets_grid
1708default_aNrmGrid_params = {
1709 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1710 "aXtraMax": 40.0, # Maximum end-of-period "assets above minimum" value
1711 "aXtraNestFac": 2, # Exponential nesting factor for aXtraGrid
1712 "aXtraCount": 96, # Number of points in the grid of "assets above minimum"
1713 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
1714}
1716# Default parameters to make mLvlGrid using make_market_resources_grid
1717default_mNrmGrid_params = {
1718 "mNrmMin": 0.001,
1719 "mNrmMax": 40.0,
1720 "mNrmNestFac": 2,
1721 "mNrmCount": 72,
1722 "mNrmExtra": None,
1723}
1725# Default parameters to make kLvlGrid using make_capital_grid
1726default_kLvlGrid_params = {
1727 "kLvlMin": 0.0,
1728 "kLvlMax": 200,
1729 "kLvlCount": 250,
1730 "kLvlOrder": 2.0,
1731}
1733# Default "basic" parameters
1734med_ext_marg_basic_params = {
1735 "constructors": med_ext_marg_constructors,
1736 "cycles": 1, # Lifecycle by default
1737 "T_cycle": 1, # Number of periods in the default sequence
1738 "T_age": None,
1739 "AgentCount": 10000, # Number of agents to simulate
1740 "DiscFac": 0.96, # intertemporal discount factor
1741 "CRRA": 1.5, # coefficient of relative risk aversion
1742 "Rfree": [1.02], # risk free interest factor
1743 "LivPrb": [0.99], # survival probability
1744 "MedCostBot": -3.1, # lower bound of medical cost distribution, in stdevs
1745 "MedCostTop": 5.2, # upper bound of medical cost distribution, in stdevs
1746 "MedCostCount": 84, # number of nodes in medical cost discretization
1747 "MedShkLogMean": [-2.0], # mean of log utility shocks
1748 "MedShkLogStd": [1.5], # standard deviation of log utility shocks
1749 "MedCostLogMean": [-1.0], # mean of log medical expenses
1750 "MedCostLogStd": [1.0], # standard deviation of log medical expenses
1751 "MedCorr": [0.3], # correlation coefficient between utility shock and expenses
1752 "PermGroFacAgg": 1.0, # Aggregate productivity growth rate (leave at 1)
1753 "NewbornTransShk": True, # Whether "newborns" have a transitory income shock
1754}
1756# Combine the dictionaries into a single default dictionary
1757init_med_ext_marg = med_ext_marg_basic_params.copy()
1758init_med_ext_marg.update(default_IncomeProcess_params)
1759init_med_ext_marg.update(default_aNrmGrid_params)
1760init_med_ext_marg.update(default_mNrmGrid_params)
1761init_med_ext_marg.update(default_kLvlGrid_params)
1762init_med_ext_marg.update(default_kNrmInitDstn_params_ExtMarg)
1763init_med_ext_marg.update(default_BeqParam_dict)
1766class MedExtMargConsumerType(PersistentShockConsumerType):
1767 r"""
1768 Class for representing agents in the extensive margin medical expense model.
1769 Such agents have labor income dynamics identical to the "general income process"
1770 model (permanent income is not normalized out), and also experience a medical
1771 shock with two components: medical cost and utility loss. They face a binary
1772 choice of whether to pay the cost or suffer the loss, then make a consumption-
1773 saving decision as normal. To simplify the computation, the joint distribution
1774 of medical shocks is specified as bivariate lognormal. This can be loosened to
1775 accommodate insurance contracts as mappings from total to out-of-pocket expenses.
1776 Can also be extended to include a health process.
1778 .. math::
1779 \begin{eqnarray*}
1780 V_t(M_t,P_t) &=& \max_{C_t, D_t} U_t(C_t) - (1-D_t) \eta_t + \beta (1-\mathsf{D}_{t+1}) \mathbb{E} [V_{t+1}(M_{t+1}, P_{t+1}], \\
1781 A_t &=& M_t - C_t - D_t med_t, \\
1782 A_t/ &\geq& 0, \\
1783 D_t &\in& \{0,1\}, \\
1784 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
1785 Y_{t+1} &=& P_{t+1} \theta_{t+1}
1786 M_{t+1} &=& R A_t + Y_{t+1}, \\
1787 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
1788 (med_t,\eta_t) &~\sim& G_t,\\
1789 U_t(C) &=& \frac{C^{1-\rho}}{1-\rho}.
1790 \end{eqnarray*}
1791 """
1793 default_ = {
1794 "params": init_med_ext_marg,
1795 "solver": solve_one_period_ConsMedExtMarg,
1796 "model": "ConsExtMargMed.yaml",
1797 "track_vars": ["aLvl", "cLvl", "Med", "mLvl", "pLvl"],
1798 }
1800 time_vary_ = [
1801 "Rfree",
1802 "LivPrb",
1803 "MedShkLogMean",
1804 "MedShkLogStd",
1805 "MedCostLogMean",
1806 "MedCostLogStd",
1807 "MedCorr",
1808 "pLogGrid",
1809 "pLvlMean",
1810 "TranShkDstn",
1811 "pLogMrkvArray",
1812 "pLvlNextFunc",
1813 "IncShkDstn",
1814 "MedShockDstn",
1815 ]
1816 time_inv_ = [
1817 "DiscFac",
1818 "CRRA",
1819 "BeqFac",
1820 "BeqShift",
1821 "MedCostBot",
1822 "MedCostTop",
1823 "MedCostCount",
1824 "aNrmGrid",
1825 "mNrmGrid",
1826 "kLvlGrid",
1827 ]
1828 shock_vars = ["PermShk", "TranShk", "MedShk", "MedCost"]
1830 def get_shocks(self):
1831 """
1832 Gets permanent and transitory income shocks for this period as well as
1833 medical need and cost shocks.
1834 """
1835 # Get permanent and transitory income shocks
1836 PersistentShockConsumerType.get_shocks(self)
1838 # Initialize medical shock array and cost of care array
1839 MedShkNow = np.zeros(self.AgentCount)
1840 MedCostNow = np.zeros(self.AgentCount)
1842 # Get shocks for each period of the cycle
1843 for t in range(self.T_cycle):
1844 these = t == self.t_cycle
1845 if np.any(these):
1846 N = np.sum(these)
1847 dstn_t = self.MedShockDstn[t]
1848 draws_t = dstn_t.draw(N)
1849 MedCostNow[these] = draws_t[0, :]
1850 MedShkNow[these] = draws_t[1, :]
1851 self.shocks["MedShk"] = MedShkNow
1852 self.shocks["MedCost"] = MedCostNow
1854 def get_controls(self):
1855 """
1856 Finds consumption for each agent, along with whether or not they get care.
1857 """
1858 # Initialize output
1859 cLvlNow = np.empty(self.AgentCount)
1860 CareNow = np.zeros(self.AgentCount, dtype=bool)
1862 # Get states and shocks
1863 mLvl = self.state_now["mLvl"]
1864 pLvl = self.state_now["pLvl"]
1865 MedCost = self.shocks["MedCost"]
1866 MedShk = self.shocks["MedShk"]
1868 # Find remaining resources with and without care
1869 bLvl_no_care = mLvl
1870 bLvl_with_care = mLvl - MedCost
1872 # Get controls for each period of the cycle
1873 for t in range(self.T_cycle):
1874 these = t == self.t_cycle
1875 if np.any(these):
1876 vFunc_t = self.solution[t].vFuncMid
1877 cFunc_t = self.solution[t].cFunc
1879 v_no_care = vFunc_t(bLvl_no_care[these], pLvl[these]) - MedShk[these]
1880 v_if_care = vFunc_t(bLvl_with_care[these], pLvl[these])
1881 get_care = v_if_care > v_no_care
1883 b_temp = bLvl_no_care[these]
1884 b_temp[get_care] = bLvl_with_care[get_care]
1885 cLvlNow[these] = cFunc_t(b_temp, pLvl[these])
1886 CareNow[these] = get_care
1888 # Store the results
1889 self.controls["cLvl"] = cLvlNow
1890 self.controls["Care"] = CareNow
1892 def get_poststates(self):
1893 """
1894 Calculates end-of-period assets for each consumer of this type.
1895 """
1896 self.state_now["Med"] = self.shocks["MedCost"] * self.controls["Care"]
1897 self.state_now["aLvl"] = (
1898 self.state_now["mLvl"] - self.controls["cLvl"] - self.state_now["Med"]
1899 )
1900 # Move now to prev
1901 AgentType.get_poststates(self)