Coverage for HARK / ConsumptionSaving / ConsMedModel.py: 98%
559 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +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_constant_mean,
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):
279 r"""
280 Constructs discretized lognormal distributions of medical preference shocks
281 for each period in the cycle.
283 .. math::
284 \text{ medShk}_t \sim \exp(\mathcal{N}(\textbf{MedShkStd}^2)) \\
285 \mathbb{E}[\text{medShk}_t]=\textbf{MedShkAvg}
288 Parameters
289 ----------
290 T_cycle : int
291 Number of non-terminal periods in the agent's cycle.
292 MedShkAvg : [float]
293 Mean of medical needs shock in each period of the problem.
294 MedShkStd : [float]
295 Standard deviation of log medical needs shock in each period of the problem.
296 MedShkCount : int
297 Number of equiprobable nodes in the "body" of the discretization.
298 MedShkCountTail : int
299 Number of nodes in each "tail" of the discretization.
300 RNG : RandomState
301 The AgentType's internal random number generator.
302 MedShkTailBound : [float,float]
303 CDF bounds for the tail of the discretization.
305 Returns
306 -------
307 MedShkDstn : [DiscreteDistribuion]
308 """
309 MedShkDstn = [] # empty list for medical shock distribution each period
310 for t in range(T_cycle):
311 # get shock distribution parameters
312 MedShkAvg_t = MedShkAvg[t]
313 MedShkStd_t = MedShkStd[t]
314 MedShkDstn_t = Lognormal(
315 mu=np.log(MedShkAvg_t) - 0.5 * MedShkStd_t**2,
316 sigma=MedShkStd_t,
317 seed=RNG.integers(0, 2**31 - 1),
318 ).discretize(
319 N=MedShkCount,
320 method="equiprobable",
321 tail_N=MedShkCountTail,
322 tail_bound=MedShkTailBound,
323 )
324 MedShkDstn_t = add_discrete_outcome_constant_mean(
325 MedShkDstn_t, 0.0, 0.0, sort=True
326 ) # add point at zero with no probability
327 MedShkDstn.append(MedShkDstn_t)
328 return MedShkDstn
331def make_continuous_MedShockDstn(
332 MedShkLogMean, MedShkLogStd, MedCostLogMean, MedCostLogStd, MedCorr, T_cycle, RNG
333):
334 """
335 Construct a time-varying list of bivariate lognormals for the medical shocks
336 distribution. This representation uses fully continuous distributions, with
337 no discretization in either dimension.
339 Parameters
340 ----------
341 MedShkLogMean : [float]
342 Age-varying list of means of log medical needs (utility) shocks.
343 MedShkLogStd : [float]
344 Age-varying list of standard deviations of log medical needs (utility) shocks.
345 MedCostLogMean : [float]
346 Age-varying list of means of log medical expense shocks.
347 MedCostLogStd : [float]
348 Age-varying list of standard deviations of log medical expense shocks..
349 MedCorr : [float]
350 Age-varying correlation coefficient between log medical expenses and utility shocks.
351 T_cycle : int
352 Number of periods in the agent's sequence.
353 RNG : RandomState
354 Random number generator for this type.
356 Returns
357 -------
358 MedShockDstn : [MultivariateLognormal]
359 Age-varying list of bivariate lognormal distributions, ordered as (MedCost,MedShk).
360 """
361 MedShockDstn = []
362 for t in range(T_cycle):
363 s1 = MedCostLogStd[t]
364 s2 = MedShkLogStd[t]
365 diag = MedCorr[t] * s1 * s2
366 S = np.array([[s1**2, diag], [diag, s2**2]])
367 M = np.array([MedCostLogMean[t], MedShkLogMean[t]])
368 seed_t = RNG.integers(0, 2**31 - 1)
369 dstn_t = MultivariateLogNormal(mu=M, Sigma=S, seed=seed_t)
370 MedShockDstn.append(dstn_t)
371 return MedShockDstn
374def make_MedShock_solution_terminal(
375 CRRA,
376 CRRAmed,
377 MedShkDstn,
378 MedPrice,
379 MedShift,
380 aXtraGrid,
381 pLvlGrid,
382 qFunc,
383 CubicBool,
384):
385 """
386 Construct the terminal period solution for this type. Similar to other models,
387 optimal behavior involves spending all available market resources; however,
388 the agent must split his resources between consumption and medical care.
389 """
390 # Take last period data
391 MedPrice_T = MedPrice[-1]
392 MedShkDstn_T = MedShkDstn[-1]
393 pLvlGrid_T = pLvlGrid[-1]
394 MedShift_T = MedShift[-1]
396 # Make the policy functions for the terminal period
397 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
398 xFunc_terminal = TrilinearInterp(
399 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
400 trivial_grid,
401 trivial_grid,
402 trivial_grid,
403 )
404 PolicyFunc_terminal = cAndMedFunc(
405 xFunc_terminal,
406 qFunc,
407 MedShift_T,
408 MedPrice_T,
409 )
411 # Make state space grids for the terminal period
412 xLvlMin = np.min(aXtraGrid) * np.min(pLvlGrid_T)
413 xLvlMax = np.max(aXtraGrid) * np.max(pLvlGrid_T)
414 xLvlGrid = make_grid_exp_mult(xLvlMin, xLvlMax, 3 * aXtraGrid.size, 4)
415 MedShkGrid = MedShkDstn_T.atoms[0]
416 MedShkPrbs = MedShkDstn_T.pmv
418 # Calculate optimal consumption on a grid of market resources and medical shocks
419 mLvlGrid = xLvlGrid
420 mLvlGrid_tiled = np.tile(
421 np.reshape(mLvlGrid, (mLvlGrid.size, 1)), (1, MedShkGrid.size)
422 )
423 # Permanent income irrelevant in terminal period
424 pLvlGrid_tiled = np.ones_like(mLvlGrid_tiled)
425 MedShkGrid_tiled = np.tile(
426 np.reshape(MedShkGrid, (1, MedShkGrid.size)), (mLvlGrid.size, 1)
427 )
428 cLvlGrid, MedGrid, xTrash = PolicyFunc_terminal(
429 mLvlGrid_tiled, pLvlGrid_tiled, MedShkGrid_tiled
430 )
432 # Integrate marginal value across shocks to get expected marginal value
433 vPgrid = cLvlGrid ** (-CRRA)
434 vPgrid[np.isinf(vPgrid)] = 0.0 # correct for issue at bottom edges
435 PrbGrid = np.tile(np.reshape(MedShkPrbs, (1, MedShkGrid.size)), (mLvlGrid.size, 1))
436 vP_expected = np.sum(vPgrid * PrbGrid, axis=1)
438 # Construct the marginal (marginal) value function for the terminal period
439 vPnvrs = np.insert(vP_expected ** (-1.0 / CRRA), 0, 0.0)
440 vPnvrsFunc = BilinearInterp(
441 np.tile(np.reshape(vPnvrs, (vPnvrs.size, 1)), (1, trivial_grid.size)),
442 np.insert(mLvlGrid, 0, 0.0),
443 trivial_grid,
444 )
445 vPfunc_terminal = MargValueFuncCRRA(vPnvrsFunc, CRRA)
446 vPPfunc_terminal = MargMargValueFuncCRRA(vPnvrsFunc, CRRA)
448 # Integrate value across shocks to get expected value
449 vGrid = utility(cLvlGrid, rho=CRRA) + utility(
450 (MedGrid + MedShift_T) / MedShkGrid_tiled, rho=CRRAmed
451 )
452 vGrid[np.isinf(vGrid)] = 0.0 # correct for issue at bottom edges
453 v_expected = np.sum(vGrid * PrbGrid, axis=1)
455 # Construct the value function for the terminal period
456 vNvrs = np.insert(utility_inv(v_expected, rho=CRRA), 0, 0.0)
457 vNvrsP = vP_expected * utility_invP(v_expected, rho=CRRA)
458 vNvrsP = np.insert(vNvrsP, 0, 0.0)
459 # tempFunc = CubicInterp(np.insert(mLvlGrid, 0, 0.0), vNvrs, vNvrsP)
460 tempFunc = LinearInterp(np.insert(mLvlGrid, 0, 0.0), vNvrs)
461 vNvrsFunc = LinearInterpOnInterp1D([tempFunc, tempFunc], trivial_grid)
462 vFunc_terminal = ValueFuncCRRA(vNvrsFunc, CRRA)
464 # Make and return the terminal period solution
465 solution_terminal = {
466 "distance_criteria": ["vPfunc"],
467 "PolicyFunc": PolicyFunc_terminal,
468 "vFunc": vFunc_terminal,
469 "vPfunc": vPfunc_terminal,
470 "vPPfunc": vPPfunc_terminal,
471 "hLvl": ConstantFunction(0.0),
472 "mLvlMin": ConstantFunction(0.0),
473 }
474 return solution_terminal
477###############################################################################
480def solve_one_period_ConsMedShock(
481 solution_next,
482 IncShkDstn,
483 MedShkDstn,
484 LivPrb,
485 DiscFac,
486 CRRA,
487 CRRAmed,
488 MedShift,
489 Rfree,
490 MedPrice,
491 pLvlNextFunc,
492 BoroCnstArt,
493 aXtraGrid,
494 pLvlGrid,
495 vFuncBool,
496 CubicBool,
497 qFunc,
498):
499 """
500 Class for solving the one period problem for the "medical shocks" model, in
501 which consumers receive shocks to permanent and transitory income as well as
502 shocks to "medical need"-- utility shocks for a second good.
504 Parameters
505 ----------
506 solution_next : dict
507 The solution to next period's one period problem, represented as a dictionary.
508 IncShkDstn : distribution.Distribution
509 A discrete approximation to the income process between the period being
510 solved and the one immediately following (in solution_next).
511 MedShkDstn : distribution.Distribution
512 Discrete distribution of the need shifter for medical care.
513 LivPrb : float
514 Survival probability; likelihood of being alive at the beginning of
515 the succeeding period.
516 DiscFac : float
517 Intertemporal discount factor for future utility.
518 CRRA : float
519 Coefficient of relative risk aversion for composite consumption.
520 CRRAmed : float
521 Coefficient of relative risk aversion for medical care.
522 MedShift : float
523 Additive shifter on medical care in the utility function. Could represent
524 self-care for small medical needs that would not be observed / reported.
525 Rfree : float
526 Risk free interest factor on end-of-period assets.
527 MedPrice : float
528 Price of unit of medical care relative to unit of consumption.
529 pLvlNextFunc : float
530 Expected permanent income next period as a function of current pLvl.
531 BoroCnstArt: float or None
532 Borrowing constraint for the minimum allowable assets to end the
533 period with.
534 aXtraGrid: np.array
535 Array of "extra" end-of-period (normalized) asset values-- assets
536 above the absolute minimum acceptable level.
537 pLvlGrid: np.array
538 Array of permanent income levels at which to solve the problem.
539 vFuncBool: boolean
540 An indicator for whether the value function should be computed and
541 included in the reported solution.
542 CubicBool: boolean
543 An indicator for whether the solver should use cubic or linear interpolation.
544 qFunc : TransConShareFunc
545 Function that maps spending xLvl and "effective" medical needs shock MedShkEff
546 to a transformation of the consumption share.
548 Returns
549 -------
550 solution_now : dict
551 Solution to this period's consumption-saving problem, as a dictionary.
552 """
553 # Define the utility functions for this period
554 uCon = UtilityFuncCRRA(CRRA)
555 uMed = UtilityFuncCRRA(CRRAmed) # Utility function for normalized medical care
556 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
558 # Unpack next period's income shock distribution
559 ShkPrbsNext = IncShkDstn.pmv
560 PermShkValsNext = IncShkDstn.atoms[0]
561 TranShkValsNext = IncShkDstn.atoms[1]
562 PermShkMinNext = np.min(PermShkValsNext)
563 TranShkMinNext = np.min(TranShkValsNext)
564 MedShkPrbs = MedShkDstn.pmv
565 MedShkVals = MedShkDstn.atoms.flatten()
567 # Calculate the probability that we get the worst possible income draw
568 IncNext = PermShkValsNext * TranShkValsNext
569 WorstIncNext = PermShkMinNext * TranShkMinNext
570 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
571 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
573 # Unpack next period's (marginal) value function
574 vFuncNext = solution_next["vFunc"] # This is a NullFunc when vFuncBool is False
575 vPfuncNext = solution_next["vPfunc"]
576 vPPfuncNext = solution_next["vPPfunc"] # This is a NullFunc when CubicBool is False
578 # Update the bounding MPCs and PDV of human wealth:
579 PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
580 try:
581 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
582 except:
583 MPCminNow = 0.0
584 mLvlMinNext = solution_next["mLvlMin"]
586 # TODO: Deal with this unused code for the upper bound of MPC (should be a function now)
587 # Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
588 # hNrmNow = 0.0
589 # temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
590 # MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
592 # Define some functions for calculating future expectations
593 def calc_pLvl_next(S, p):
594 return pLvlNextFunc(p) * S["PermShk"]
596 def calc_mLvl_next(S, a, p_next):
597 return Rfree * a + p_next * S["TranShk"]
599 def calc_hLvl(S, p):
600 pLvl_next = calc_pLvl_next(S, p)
601 hLvl = S["TranShk"] * pLvl_next + solution_next["hLvl"](pLvl_next)
602 return hLvl
604 def calc_v_next(S, a, p):
605 pLvl_next = calc_pLvl_next(S, p)
606 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
607 v_next = vFuncNext(mLvl_next, pLvl_next)
608 return v_next
610 def calc_vP_next(S, a, p):
611 pLvl_next = calc_pLvl_next(S, p)
612 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
613 vP_next = vPfuncNext(mLvl_next, pLvl_next)
614 return vP_next
616 def calc_vPP_next(S, a, p):
617 pLvl_next = calc_pLvl_next(S, p)
618 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
619 vPP_next = vPPfuncNext(mLvl_next, pLvl_next)
620 return vPP_next
622 # Construct human wealth level as a function of productivity pLvl
623 hLvlGrid = 1.0 / Rfree * expected(calc_hLvl, IncShkDstn, args=(pLvlGrid))
624 hLvlNow = LinearInterp(np.insert(pLvlGrid, 0, 0.0), np.insert(hLvlGrid, 0, 0.0))
626 # Make temporary grids of income shocks and next period income values
627 ShkCount = TranShkValsNext.size
628 pLvlCount = pLvlGrid.size
629 PermShkVals_temp = np.tile(
630 np.reshape(PermShkValsNext, (1, ShkCount)), (pLvlCount, 1)
631 )
632 TranShkVals_temp = np.tile(
633 np.reshape(TranShkValsNext, (1, ShkCount)), (pLvlCount, 1)
634 )
635 pLvlNext_temp = (
636 np.tile(
637 np.reshape(pLvlNextFunc(pLvlGrid), (pLvlCount, 1)),
638 (1, ShkCount),
639 )
640 * PermShkVals_temp
641 )
643 # Find the natural borrowing constraint for each persistent income level
644 aLvlMin_candidates = (
645 mLvlMinNext(pLvlNext_temp) - TranShkVals_temp * pLvlNext_temp
646 ) / Rfree
647 aLvlMinNow = np.max(aLvlMin_candidates, axis=1)
648 BoroCnstNat = LinearInterp(
649 np.insert(pLvlGrid, 0, 0.0), np.insert(aLvlMinNow, 0, 0.0)
650 )
652 # Define the minimum allowable mLvl by pLvl as the greater of the natural and artificial borrowing constraints
653 if BoroCnstArt is not None:
654 BoroCnstArt = LinearInterp(np.array([0.0, 1.0]), np.array([0.0, BoroCnstArt]))
655 mLvlMinNow = UpperEnvelope(BoroCnstArt, BoroCnstNat)
656 else:
657 mLvlMinNow = BoroCnstNat
659 # Make the constrained total spending function: spend all market resources
660 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
661 SpendAllFunc = TrilinearInterp(
662 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
663 trivial_grid,
664 trivial_grid,
665 trivial_grid,
666 )
667 xFuncNowCnst = VariableLowerBoundFunc3D(SpendAllFunc, mLvlMinNow)
669 # Define grids of pLvl and aLvl on which to compute future expectations
670 pLvlCount = pLvlGrid.size
671 aNrmCount = aXtraGrid.size
672 MedCount = MedShkVals.size
673 pLvlNow = np.tile(np.reshape(pLvlGrid, (1, pLvlCount)), (aNrmCount, 1))
674 aLvlNow = np.reshape(aXtraGrid, (aNrmCount, 1)) * pLvlNow + BoroCnstNat(pLvlNow)
675 if pLvlGrid[0] == 0.0: # aLvl turns out badly if pLvl is 0 at bottom
676 aLvlNow[:, 0] = aXtraGrid
678 # Calculate end-of-period marginal value of assets
679 EndOfPrd_vP = (
680 DiscFacEff * Rfree * expected(calc_vP_next, IncShkDstn, args=(aLvlNow, pLvlNow))
681 )
683 # If the value function has been requested, construct the end-of-period vFunc
684 if vFuncBool:
685 # Compute expected value from end-of-period states
686 EndOfPrd_v = expected(calc_v_next, IncShkDstn, args=(aLvlNow, pLvlNow))
687 EndOfPrd_v *= DiscFacEff
689 # Transformed value through inverse utility function to "decurve" it
690 EndOfPrd_vNvrs = uCon.inv(EndOfPrd_v)
692 # Add points at mLvl=zero
693 EndOfPrd_vNvrs = np.concatenate(
694 (np.zeros((1, pLvlCount)), EndOfPrd_vNvrs), axis=0
695 )
697 # Make a temporary aLvl grid for interpolating the end-of-period value function
698 aLvl_temp = np.concatenate(
699 (
700 np.reshape(BoroCnstNat(pLvlGrid), (1, pLvlGrid.size)),
701 aLvlNow,
702 ),
703 axis=0,
704 )
706 # Make an end-of-period value function for each persistent income level in the grid
707 EndOfPrd_vNvrsFunc_list = []
708 for p in range(pLvlCount):
709 EndOfPrd_vNvrsFunc_list.append(
710 LinearInterp(
711 aLvl_temp[:, p] - BoroCnstNat(pLvlGrid[p]),
712 EndOfPrd_vNvrs[:, p],
713 )
714 )
715 EndOfPrd_vNvrsFuncBase = LinearInterpOnInterp1D(
716 EndOfPrd_vNvrsFunc_list, pLvlGrid
717 )
719 # Re-adjust the combined end-of-period value function to account for the
720 # natural borrowing constraint shifter and "re-curve" it
721 EndOfPrd_vNvrsFunc = VariableLowerBoundFunc2D(
722 EndOfPrd_vNvrsFuncBase, BoroCnstNat
723 )
724 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
726 # Solve the first order condition to get optimal consumption and medical
727 # spending, then find the endogenous mLvl gridpoints
728 cLvlNow = np.tile(
729 np.reshape(uCon.derinv(EndOfPrd_vP, order=(1, 0)), (aNrmCount, pLvlCount, 1)),
730 (1, 1, MedCount),
731 )
732 MedShkVals_adj = np.reshape(MedShkVals ** (1.0 - 1.0 / CRRAmed), (1, 1, MedCount))
733 MedBase = MedPrice ** (-1.0 / CRRAmed) * uMed.derinv(EndOfPrd_vP, order=(1, 0))
734 MedLvlNow = np.reshape(MedBase, (aNrmCount, pLvlCount, 1)) * MedShkVals_adj
735 MedLvlNow = np.maximum(MedLvlNow - MedShift, 0.0)
736 aLvlNow_tiled = np.tile(
737 np.reshape(aLvlNow, (aNrmCount, pLvlCount, 1)), (1, 1, MedCount)
738 )
739 xLvlNow = cLvlNow + MedPrice * MedLvlNow
740 mLvlNow = xLvlNow + aLvlNow_tiled
742 # Limiting consumption is zero as m approaches the natural borrowing constraint
743 x_for_interpolation = np.concatenate(
744 (np.zeros((1, pLvlCount, MedCount)), xLvlNow), axis=0
745 )
746 temp = np.tile(
747 BoroCnstNat(np.reshape(pLvlGrid, (1, pLvlCount, 1))),
748 (1, 1, MedCount),
749 )
750 m_for_interpolation = np.concatenate((temp, mLvlNow), axis=0)
752 # Make a 3D array of permanent income for interpolation
753 p_for_interpolation = np.tile(
754 np.reshape(pLvlGrid, (1, pLvlCount, 1)), (aNrmCount + 1, 1, MedCount)
755 )
756 MedShkVals_tiled = np.tile( # This does *not* have the CRRA adjustment
757 np.reshape(MedShkVals, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
758 )
760 # Build the set of cFuncs by pLvl, gathered in a list
761 xFunc_by_pLvl_and_MedShk = [] # Initialize the empty list of lists of 1D xFuncs
762 if CubicBool:
763 raise NotImplementedError(
764 "Cubic spline interpolation has not been implemented yet!"
765 )
767 # Basic version: use linear interpolation within a pLvl and MedShk
768 else:
769 # Loop over pLvl and then MedShk within that
770 for i in range(pLvlCount):
771 temp_list = []
772 pLvl_i = p_for_interpolation[0, i, 0]
773 mLvlMin_i = BoroCnstNat(pLvl_i)
774 for j in range(MedCount):
775 m_temp = m_for_interpolation[:, i, j] - mLvlMin_i
776 x_temp = x_for_interpolation[:, i, j]
777 temp_list.append(LinearInterp(m_temp, x_temp))
778 xFunc_by_pLvl_and_MedShk.append(deepcopy(temp_list))
780 # Combine the nested list of linear xFuncs into a single function
781 pLvl_temp = p_for_interpolation[0, :, 0]
782 MedShk_temp = MedShkVals_tiled[0, 0, :]
783 xFuncUncBase = BilinearInterpOnInterp1D(
784 xFunc_by_pLvl_and_MedShk, pLvl_temp, MedShk_temp
785 )
786 xFuncNowUnc = VariableLowerBoundFunc3D(xFuncUncBase, BoroCnstNat)
787 # Re-adjust for lower bound of natural borrowing constraint
789 # Combine the constrained and unconstrained functions into the true consumption function
790 xFuncNow = LowerEnvelope3D(xFuncNowUnc, xFuncNowCnst)
792 # Build the policy function for this period as a single object
793 PolicyFuncNow = cAndMedFunc(xFuncNow, qFunc, MedShift, MedPrice)
795 # Make the marginal value function by integrating over medical shocks
796 # Make temporary grids to evaluate the consumption function
797 temp_grid = np.tile(
798 np.reshape(aXtraGrid, (aNrmCount, 1, 1)), (1, pLvlCount, MedCount)
799 )
800 aMinGrid = np.tile(
801 np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount, 1)),
802 (aNrmCount, 1, MedCount),
803 )
804 pGrid = np.tile(np.reshape(pLvlGrid, (1, pLvlCount, 1)), (aNrmCount, 1, MedCount))
805 mGrid = temp_grid * pGrid + aMinGrid
806 if pLvlGrid[0] == 0:
807 mGrid[:, 0, :] = np.tile(np.reshape(aXtraGrid, (aNrmCount, 1)), (1, MedCount))
808 MedShkGrid = np.tile(
809 np.reshape(MedShkVals, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
810 )
811 ProbsGrid = np.reshape(MedShkPrbs, (1, 1, MedCount))
813 # Get optimal consumption (and medical care) for each state
814 cGrid, MedGrid, xTrash = PolicyFuncNow(mGrid, pGrid, MedShkGrid)
816 # Calculate expected marginal value by "integrating" across medical shocks
817 vPgrid = uCon.der(cGrid)
818 vPnow = np.sum(vPgrid * ProbsGrid, axis=2)
820 # Add vPnvrs=0 at m=mLvlMin to close it off at the bottom (and vNvrs=0)
821 mGrid_small = np.concatenate(
822 (np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount)), mGrid[:, :, 0])
823 )
824 vPnvrsNow = np.concatenate(
825 (np.zeros((1, pLvlCount)), uCon.derinv(vPnow, order=(1, 0)))
826 )
828 # Calculate expected value by "integrating" across medical shocks
829 if vFuncBool:
830 # interpolation error sometimes makes Med < 0 (barely), so fix that
831 MedGrid = np.maximum(MedGrid, 1e-100)
832 # interpolation error sometimes makes tiny violations, so fix that
833 aGrid = np.maximum(mGrid - cGrid - MedPrice * MedGrid, aMinGrid)
834 MedEff = (MedGrid + MedShift) / MedShkGrid
835 vGrid = uCon(cGrid) + uMed(MedEff) + EndOfPrd_vFunc(aGrid, pGrid)
836 vNow = np.sum(vGrid * ProbsGrid, axis=2)
838 # Switch to pseudo-inverse value and add a point at bottom
839 vNvrsNow = np.concatenate((np.zeros((1, pLvlCount)), uCon.inv(vNow)), axis=0)
841 # Construct the pseudo-inverse value and marginal value functions over mLvl,pLvl
842 vPnvrsFunc_by_pLvl = []
843 vNvrsFunc_by_pLvl = []
844 # Make a pseudo inverse marginal value function for each pLvl
845 for j in range(pLvlCount):
846 pLvl = pLvlGrid[j]
847 m_temp = mGrid_small[:, j] - mLvlMinNow(pLvl)
848 vPnvrs_temp = vPnvrsNow[:, j]
849 vPnvrsFunc_by_pLvl.append(LinearInterp(m_temp, vPnvrs_temp))
850 if vFuncBool:
851 vNvrs_temp = vNvrsNow[:, j]
852 vNvrsFunc_by_pLvl.append(LinearInterp(m_temp, vNvrs_temp))
854 # Combine those functions across pLvls, and adjust for the lower bound of mLvl
855 vPnvrsFuncBase = LinearInterpOnInterp1D(vPnvrsFunc_by_pLvl, pLvlGrid)
856 vPnvrsFunc = VariableLowerBoundFunc2D(vPnvrsFuncBase, mLvlMinNow)
857 if vFuncBool:
858 vNvrsFuncBase = LinearInterpOnInterp1D(vNvrsFunc_by_pLvl, pLvlGrid)
859 vNvrsFunc = VariableLowerBoundFunc2D(vNvrsFuncBase, mLvlMinNow)
861 # "Re-curve" the (marginal) value function
862 vPfuncNow = MargValueFuncCRRA(vPnvrsFunc, CRRA)
863 if vFuncBool:
864 vFuncNow = ValueFuncCRRA(vNvrsFunc, CRRA)
865 else:
866 vFuncNow = NullFunc()
868 # If using cubic spline interpolation, construct the marginal marginal value function
869 if CubicBool:
870 vPPfuncNow = MargMargValueFuncCRRA(vPfuncNow.cFunc, CRRA)
871 else:
872 vPPfuncNow = NullFunc()
874 # Package and return the solution as a dictionary
875 solution_now = {
876 "distance_criteria": ["vPfunc"],
877 "PolicyFunc": PolicyFuncNow,
878 "vFunc": vFuncNow,
879 "vPfunc": vPfuncNow,
880 "vPPfunc": vPPfuncNow,
881 "hLvl": hLvlNow,
882 "mLvlMin": mLvlMinNow,
883 }
884 return solution_now
887###############################################################################
889# Make a constructor dictionary for the general income process consumer type
890medshock_constructor_dict = {
891 "IncShkDstn": construct_lognormal_income_process_unemployment,
892 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
893 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
894 "aXtraGrid": make_assets_grid,
895 "pLvlPctiles": make_basic_pLvlPctiles,
896 "pLvlGrid": make_pLvlGrid_by_simulation,
897 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
898 "MedShkDstn": make_lognormal_MedShkDstn,
899 "qFunc": make_qFunc,
900 "solution_terminal": make_MedShock_solution_terminal,
901 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
902 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
903}
905# Make a dictionary with parameters for the default constructor for kNrmInitDstn
906default_kNrmInitDstn_params = {
907 "kLogInitMean": -6.0, # Mean of log initial capital
908 "kLogInitStd": 1.0, # Stdev of log initial capital
909 "kNrmInitCount": 15, # Number of points in initial capital discretization
910}
912# Make a dictionary with parameters for the default constructor for pLvlInitDstn
913default_pLvlInitDstn_params = {
914 "pLogInitMean": 0.0, # Mean of log permanent income
915 "pLogInitStd": 0.4, # Stdev of log permanent income
916 "pLvlInitCount": 15, # Number of points in initial capital discretization
917}
919# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
920default_IncShkDstn_params = {
921 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
922 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
923 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
924 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
925 "UnempPrb": 0.05, # Probability of unemployment while working
926 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
927 "T_retire": 0, # Period of retirement (0 --> no retirement)
928 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
929 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
930}
932# Default parameters to make aXtraGrid using make_assets_grid
933default_aXtraGrid_params = {
934 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
935 "aXtraMax": 30, # Maximum end-of-period "assets above minimum" value
936 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
937 "aXtraCount": 32, # Number of points in the grid of "assets above minimum"
938 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
939}
941# Default parameters to make pLvlGrid using make_basic_pLvlPctiles
942default_pLvlPctiles_params = {
943 "pLvlPctiles_count": 19, # Number of points in the "body" of the grid
944 "pLvlPctiles_bound": [0.05, 0.95], # Percentile bounds of the "body"
945 "pLvlPctiles_tail_count": 4, # Number of points in each tail of the grid
946 "pLvlPctiles_tail_order": np.e, # Scaling factor for points in each tail
947}
949# Default parameters to make pLvlGrid using make_trivial_pLvlNextFunc
950default_pLvlGrid_params = {
951 "pLvlInitMean": 0.0, # Mean of log initial permanent income
952 "pLvlInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
953 # "pLvlPctiles": pLvlPctiles, # Percentiles of permanent income to use for the grid
954 "pLvlExtra": [
955 0.0001
956 ], # Additional permanent income points to automatically add to the grid, optional
957}
959# Default parameters to make pLvlNextFunc using make_AR1_style_pLvlNextFunc
960default_pLvlNextFunc_params = {
961 "PermGroFac": [1.0], # Permanent income growth factor
962 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
963}
965# Default parameters to make MedShkDstn using make_lognormal_MedShkDstn
966default_MedShkDstn_params = {
967 "MedShkAvg": [0.1], # Average of medical need shocks
968 "MedShkStd": [1.5], # Standard deviation of (log) medical need shocks
969 "MedShkCount": 5, # Number of medical shock points in "body"
970 "MedShkCountTail": 15, # Number of medical shock points in "tail" (upper only)
971 "MedPrice": [1.5], # Relative price of a unit of medical care
972}
974# Make a dictionary to specify a medical shocks consumer type
975init_medical_shocks = {
976 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
977 "cycles": 1, # Finite, non-cyclic model
978 "T_cycle": 1, # Number of periods in the cycle for this agent type
979 "pseudo_terminal": False, # Terminal period really does exist
980 "constructors": medshock_constructor_dict, # See dictionary above
981 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
982 "CRRA": 2.0, # Coefficient of relative risk aversion on consumption
983 "CRRAmed": 5.0, # Coefficient of relative risk aversion on medical care
984 "MedShift": [1e-8], # Amount of medical care the agent can self-provide
985 "Rfree": [1.03], # Interest factor on retained assets
986 "DiscFac": 0.96, # Intertemporal discount factor
987 "LivPrb": [0.99], # Survival probability after each period
988 "BoroCnstArt": 0.0, # Artificial borrowing constraint
989 "vFuncBool": False, # Whether to calculate the value function during solution
990 "CubicBool": False, # Whether to use cubic spline interpolation when True
991 # (Uses linear spline interpolation for cFunc when False)
992 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
993 "AgentCount": 10000, # Number of agents of this type
994 "T_age": None, # Age after which simulated agents are automatically killed
995 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
996 # (The portion of PermGroFac attributable to aggregate productivity growth)
997 "NewbornTransShk": False, # Whether Newborns have transitory shock
998 # ADDITIONAL OPTIONAL PARAMETERS
999 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1000 # (Forces Newborns to follow solution path of the agent they replaced if True)
1001 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1002}
1003init_medical_shocks.update(default_IncShkDstn_params)
1004init_medical_shocks.update(default_aXtraGrid_params)
1005init_medical_shocks.update(default_pLvlPctiles_params)
1006init_medical_shocks.update(default_pLvlGrid_params)
1007init_medical_shocks.update(default_MedShkDstn_params)
1008init_medical_shocks.update(default_pLvlNextFunc_params)
1009init_medical_shocks.update(default_pLvlInitDstn_params)
1010init_medical_shocks.update(default_kNrmInitDstn_params)
1013class MedShockConsumerType(PersistentShockConsumerType):
1014 r"""
1015 A consumer type based on GenIncShockConsumerType, with two types of consumption goods (medical and nonmedical) and random shocks to medical utility.
1017 .. math::
1018 \begin{eqnarray*}
1019 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})], \\
1020 A_t &=& M_t - X_t, \\
1021 X_t &=& C_t +med_t \textbf{ medPrice}_t,\\
1022 A_t/P_t &\geq& \underline{a}, \\
1023 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
1024 Y_{t+1} &=& P_{t+1} \theta_{t+1}
1025 M_{t+1} &=& R A_t + Y_{t+1}, \\
1026 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
1027 \eta_t &~\sim& G_t,\\
1028 U_t(C, med; \eta) &=& \frac{C^{1-\rho}}{1-\rho} +\eta \frac{med^{1-\nu}}{1-\nu}.
1029 \end{eqnarray*}
1032 Constructors
1033 ------------
1034 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
1035 The agent's income shock distributions.
1036 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
1037 aXtraGrid: Constructor
1038 The agent's asset grid.
1039 Its default constructor is :func:`HARK.utilities.make_assets_grid`
1040 pLvlNextFunc: Constructor
1041 An arbitrary function used to evolve the GenIncShockConsumerType's permanent income
1042 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_trivial_pLvlNextFunc`
1043 pLvlGrid: Constructor
1044 The agent's pLvl grid
1045 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_pLvlGrid_by_simulation`
1046 pLvlPctiles: Constructor
1047 The agents income level percentile grid
1048 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_basic_pLvlPctiles`
1049 MedShkDstn: Constructor, :math:`\text{medShk}`
1050 The agent's Medical utility shock distribution.
1051 Its default constructor is :func:`HARK.ConsumptionSaving.ConsMedModel.make_lognormal_MedShkDstn`
1053 Solving Parameters
1054 ------------------
1055 cycles: int
1056 0 specifies an infinite horizon model, 1 specifies a finite model.
1057 T_cycle: int
1058 Number of periods in the cycle for this agent type.
1059 CRRA: float, :math:`\rho`
1060 Coefficient of Relative Risk Aversion for consumption.
1061 CRRAmed: float, :math:`\nu`
1062 Coefficient of Relative Risk Aversion for medical care.
1063 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1064 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1065 DiscFac: float, :math:`\beta`
1066 Intertemporal discount factor.
1067 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1068 Survival probability after each period.
1069 PermGroFac: list[float], time varying, :math:`\Gamma`
1070 Permanent income growth factor.
1071 BoroCnstArt: float, :math:`\underline{a}`
1072 The minimum Asset/Perminant Income ratio, None to ignore.
1073 vFuncBool: bool
1074 Whether to calculate the value function during solution.
1075 CubicBool: bool
1076 Whether to use cubic spline interpoliation.
1078 Simulation Parameters
1079 ---------------------
1080 AgentCount: int
1081 Number of agents of this kind that are created during simulations.
1082 T_age: int
1083 Age after which to automatically kill agents, None to ignore.
1084 T_sim: int, required for simulation
1085 Number of periods to simulate.
1086 track_vars: list[strings]
1087 List of variables that should be tracked when running the simulation.
1088 For this agent, the options are 'Med', 'MedShk', 'PermShk', 'TranShk', 'aLvl', 'cLvl', 'mLvl', 'pLvl', and 'who_dies'.
1090 PermShk is the agent's permanent income shock
1092 MedShk is the agent's medical utility shock
1094 TranShk is the agent's transitory income shock
1096 aLvl is the nominal asset level
1098 cLvl is the nominal consumption level
1100 Med is the nominal medical spending level
1102 mLvl is the nominal market resources
1104 pLvl is the permanent income level
1106 who_dies is the array of which agents died
1107 aNrmInitMean: float
1108 Mean of Log initial Normalized Assets.
1109 aNrmInitStd: float
1110 Std of Log initial Normalized Assets.
1111 pLvlInitMean: float
1112 Mean of Log initial permanent income.
1113 pLvlInitStd: float
1114 Std of Log initial permanent income.
1115 PermGroFacAgg: float
1116 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1117 PerfMITShk: boolean
1118 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1119 NewbornTransShk: boolean
1120 Whether Newborns have transitory shock.
1122 Attributes
1123 ----------
1124 solution: list[Consumer solution object]
1125 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1126 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1128 Unlike other models with this solution type, this model's variables are NOT normalized.
1129 The solution functions additionally depend on the permanent income level and the medical shock.
1130 For example, :math:`C=\text{cFunc}(M,P,MedShk)`.
1131 hNrm has been replaced by hLvl which is a function of permanent income.
1132 MPC max has not yet been implemented for this class. It will be a function of permanent income.
1134 This solution has two additional functions
1135 :math:`\text{Med}=\text{MedFunc}(M,P,\text{MedShk})`: returns the agent's spending on Medical care
1137 :math:`[C,Med]=\text{policyFunc}(M,P,\text{MedShk})`: returns the agent's spending on consumption and Medical care as numpy arrays
1139 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1140 history: Dict[Array]
1141 Created by running the :func:`.simulate()` method.
1142 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1143 Visit :class:`HARK.core.AgentType.simulate` for more information.
1144 """
1146 default_ = {
1147 "params": init_medical_shocks,
1148 "solver": solve_one_period_ConsMedShock,
1149 "model": "ConsMedShock.yaml",
1150 "track_vars": ["aLvl", "cLvl", "Med", "mLvl", "pLvl"],
1151 }
1153 time_vary_ = PersistentShockConsumerType.time_vary_ + [
1154 "MedPrice",
1155 "MedShkDstn",
1156 "MedShift",
1157 ]
1158 time_inv_ = PersistentShockConsumerType.time_inv_ + ["CRRAmed", "qFunc"]
1159 shock_vars_ = PersistentShockConsumerType.shock_vars_ + ["MedShk"]
1160 state_vars = PersistentShockConsumerType.state_vars + ["mLvl"]
1161 distributions = [
1162 "IncShkDstn",
1163 "PermShkDstn",
1164 "TranShkDstn",
1165 "kNrmInitDstn",
1166 "pLvlInitDstn",
1167 "MedShkDstn",
1168 ]
1170 def pre_solve(self):
1171 self.construct("solution_terminal")
1173 def get_shocks(self):
1174 """
1175 Gets permanent and transitory income shocks for this period as well as medical need shocks
1176 and the price of medical care.
1178 Parameters
1179 ----------
1180 None
1182 Returns
1183 -------
1184 None
1185 """
1186 # Get permanent and transitory income shocks
1187 PersistentShockConsumerType.get_shocks(self)
1189 # Initialize medical shock array and relative price array
1190 MedShkNow = np.zeros(self.AgentCount)
1191 MedPriceNow = np.zeros(self.AgentCount)
1193 # Get shocks for each period of the cycle
1194 for t in range(self.T_cycle):
1195 these = t == self.t_cycle
1196 N = np.sum(these)
1197 if N > 0:
1198 MedShkNow[these] = self.MedShkDstn[t].draw(N)
1199 MedPriceNow[these] = self.MedPrice[t]
1200 self.shocks["MedShk"] = MedShkNow
1201 self.shocks["MedPrice"] = MedPriceNow
1203 def get_controls(self):
1204 """
1205 Calculates consumption and medical care for each consumer of this type
1206 using the consumption and medical care functions.
1208 Parameters
1209 ----------
1210 None
1212 Returns
1213 -------
1214 None
1215 """
1216 cLvlNow = np.empty(self.AgentCount)
1217 MedNow = np.empty(self.AgentCount)
1218 for t in range(self.T_cycle):
1219 these = t == self.t_cycle
1220 cLvlNow[these], MedNow[these], x = self.solution[t]["PolicyFunc"](
1221 self.state_now["mLvl"][these],
1222 self.state_now["pLvl"][these],
1223 self.shocks["MedShk"][these],
1224 )
1225 self.controls["cLvl"] = cLvlNow
1226 self.controls["Med"] = MedNow
1228 def get_poststates(self):
1229 """
1230 Calculates end-of-period assets for each consumer of this type.
1232 Parameters
1233 ----------
1234 None
1236 Returns
1237 -------
1238 None
1239 """
1240 self.state_now["aLvl"] = (
1241 self.state_now["mLvl"]
1242 - self.controls["cLvl"]
1243 - self.shocks["MedPrice"] * self.controls["Med"]
1244 )
1247###############################################################################
1250class ConsMedExtMargSolution(MetricObject):
1251 """
1252 Representation of the solution to one period's problem in the extensive margin
1253 medical expense model. If no inputs are passed, a trivial object is constructed,
1254 which can be used as the pseudo-terminal solution.
1256 Parameters
1257 ----------
1258 vFunc_by_pLvl : [function]
1259 List of beginning-of-period value functions over kLvl, by pLvl.
1260 vPfunc_by_pLvl : [function]
1261 List of beginning-of-period marginal functions over kLvl, by pLvl.
1262 cFunc_by_pLvl : [function]
1263 List of consumption functions over bLvl, by pLvl.
1264 vNvrsFuncMid_by_pLvl : [function]
1265 List of pseudo-inverse value function for consumption phase over bLvl, by pLvl.
1266 ExpMedFunc : function
1267 Expected medical care as a function of mLvl and pLvl, just before medical
1268 shock is realized.
1269 CareProbFunc : function
1270 Probability of getting medical treatment as a function of mLvl and pLvl,
1271 just before medical shock is realized.
1272 pLvl : np.array
1273 Grid of permanent income levels during the period (after shocks).
1274 CRRA : float
1275 Coefficient of relative risk aversion
1276 """
1278 distance_criteria = ["cFunc"]
1280 def __init__(
1281 self,
1282 vFunc_by_pLvl=None,
1283 vPfunc_by_pLvl=None,
1284 cFunc_by_pLvl=None,
1285 vNvrsFuncMid_by_pLvl=None,
1286 ExpMedFunc=None,
1287 CareProbFunc=None,
1288 pLvl=None,
1289 CRRA=None,
1290 ):
1291 self.pLvl = pLvl
1292 self.CRRA = CRRA
1293 if vFunc_by_pLvl is None:
1294 self.vFunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1295 else:
1296 self.vFunc_by_pLvl = vFunc_by_pLvl
1297 if vPfunc_by_pLvl is None:
1298 self.vPfunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1299 else:
1300 self.vPfunc_by_pLvl = vPfunc_by_pLvl
1301 if cFunc_by_pLvl is not None:
1302 self.cFunc = LinearInterpOnInterp1D(cFunc_by_pLvl, pLvl)
1303 else:
1304 self.cFunc = None
1305 if vNvrsFuncMid_by_pLvl is not None:
1306 vNvrsFuncMid = LinearInterpOnInterp1D(vNvrsFuncMid_by_pLvl, pLvl)
1307 self.vFuncMid = ValueFuncCRRA(vNvrsFuncMid, CRRA, illegal_value=-np.inf)
1308 if ExpMedFunc is not None:
1309 self.ExpMedFunc = ExpMedFunc
1310 if CareProbFunc is not None:
1311 self.CareProbFunc = CareProbFunc
1314def make_MedExtMarg_solution_terminal(pLogCount):
1315 """
1316 Construct a trivial pseudo-terminal solution for the extensive margin medical
1317 spending model: a list of constant zero functions for (marginal) value. The
1318 only piece of information needed for this is how many such functions to include.
1319 """
1320 pLvl_terminal = np.arange(pLogCount)
1321 solution_terminal = ConsMedExtMargSolution(pLvl=pLvl_terminal)
1322 return solution_terminal
1325###############################################################################
1328def solve_one_period_ConsMedExtMarg(
1329 solution_next,
1330 DiscFac,
1331 CRRA,
1332 BeqFac,
1333 BeqShift,
1334 Rfree,
1335 LivPrb,
1336 MedShkLogMean,
1337 MedShkLogStd,
1338 MedCostLogMean,
1339 MedCostLogStd,
1340 MedCorr,
1341 MedCostBot,
1342 MedCostTop,
1343 MedCostCount,
1344 aNrmGrid,
1345 pLogGrid,
1346 pLvlMean,
1347 TranShkDstn,
1348 pLogMrkvArray,
1349 mNrmGrid,
1350 kLvlGrid,
1351):
1352 """
1353 Solve one period of the "extensive margin medical care" model. Each period, the
1354 agent receives a persistent and transitory shock to income, and then a medical
1355 shock with two components: utility and cost. He makes a binary choice between
1356 paying the cost in medical expenses or suffering the utility loss, then makes
1357 his ordinary consumption-saving decision (technically made simultaneously, but
1358 solved as if sequential). This version has one health state and no insurance choice
1359 and hardcodes a liquidity constraint.
1361 Parameters
1362 ----------
1363 solution_next : ConsMedExtMargSolution
1364 Solution to the succeeding period's problem.
1365 DiscFac : float
1366 Intertemporal discount factor.
1367 CRRA : float
1368 Coefficient of relative risk aversion.
1369 BeqFac : float
1370 Scaling factor for bequest motive.
1371 BeqShift : float
1372 Shifter for bequest motive.
1373 Rfree : float
1374 Risk free return factor on saving.
1375 LivPrb : float
1376 Survival probability from this period to the next one.
1377 MedShkLogMean : float
1378 Mean of log utility shocks, assumed to be lognormally distributed.
1379 MedShkLogStd : float
1380 Stdev of log utility shocks, assumed to be lognormally distributed.
1381 MedCostLogMean : float
1382 Mean of log medical expense shocks, assumed to be lognormally distributed.
1383 MedCostLogStd : float
1384 Stdev of log medical expense shocks, assumed to be lognormally distributed.
1385 MedCorr : float
1386 Correlation coefficient betwen log utility shocks and log medical expense
1387 shocks, assumed to be joint normal (in logs).
1388 MedCostBot : float
1389 Lower bound of medical costs to consider, as standard deviations of log
1390 expenses away from the mean.
1391 MedCostTop : float
1392 Upper bound of medical costs to consider, as standard deviations of log
1393 expenses away from the mean.
1394 MedCostCount : int
1395 Number of points to use when discretizing MedCost
1396 aNrmGrid : np.array
1397 Exogenous grid of end-of-period assets, normalized by income level.
1398 pLogGrid : np.array
1399 Exogenous grid of *deviations from mean* log income level.
1400 pLvlMean : float
1401 Mean income level at this age, for generating actual income levels from
1402 pLogGrid as pLvl = pLvlMean * np.exp(pLogGrid).
1403 TranShkDstn : DiscreteDistribution
1404 Discretized transitory income shock distribution.
1405 pLogMrkvArray : np.array
1406 Markov transition array from beginning-of-period (prior) income levels
1407 to this period's levels. Pre-computed by (e.g.) Tauchen's method.
1408 mNrmGrid : np.array
1409 Exogenous grid of decision-time normalized market resources,
1410 kLvlGrid : np.array
1411 Beginning-of-period capital grid (spanning zero to very high wealth).
1413 Returns
1414 -------
1415 solution_now : ConsMedExtMargSolution
1416 Representation of the solution to this period's problem, including the
1417 beginning-of-period (marginal) value function by pLvl, the consumption
1418 function by pLvl, and the (pseudo-inverse) value function for the consumption
1419 phase (as a list by pLvl).
1420 """
1421 # Define (marginal) utility and bequest motive functions
1422 u = lambda x: CRRAutility(x, rho=CRRA)
1423 uP = lambda x: CRRAutilityP(x, rho=CRRA)
1424 W = lambda x: BeqFac * u(x + BeqShift)
1425 Wp = lambda x: BeqFac * uP(x + BeqShift)
1426 n = lambda x: CRRAutility_inv(x, rho=CRRA)
1428 # Make grids of pLvl x aLvl
1429 pLvl = np.exp(pLogGrid) * pLvlMean
1430 aLvl = np.dot(
1431 np.reshape(aNrmGrid, (aNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvl.size))
1432 )
1433 aLvl = np.concatenate([np.zeros((1, pLvl.size)), aLvl]) # add zero entries
1435 # Evaluate end-of-period marginal value at each combination of pLvl x aLvl
1436 pLvlCount = pLvl.size
1437 EndOfPrd_vP = np.empty_like(aLvl)
1438 EndOfPrd_v = np.empty_like(aLvl)
1439 for j in range(pLvlCount):
1440 EndOfPrdvFunc_this_pLvl = solution_next.vFunc_by_pLvl[j]
1441 EndOfPrdvPfunc_this_pLvl = solution_next.vPfunc_by_pLvl[j]
1442 EndOfPrd_v[:, j] = DiscFac * LivPrb * EndOfPrdvFunc_this_pLvl(aLvl[:, j])
1443 EndOfPrd_vP[:, j] = DiscFac * LivPrb * EndOfPrdvPfunc_this_pLvl(aLvl[:, j])
1444 EndOfPrd_v += (1.0 - LivPrb) * W(aLvl)
1445 EndOfPrd_vP += (1.0 - LivPrb) * Wp(aLvl)
1447 # Calculate optimal consumption for each (aLvl,pLvl) gridpoint, roll back to bLvl
1448 cLvl = CRRAutilityP_inv(EndOfPrd_vP, CRRA)
1449 bLvl = aLvl + cLvl
1451 # Construct consumption functions over bLvl for each pLvl
1452 cFunc_by_pLvl = []
1453 for j in range(pLvlCount):
1454 cFunc_j = LinearInterp(
1455 np.insert(bLvl[:, j], 0, 0.0), np.insert(cLvl[:, j], 0, 0.0)
1456 )
1457 cFunc_by_pLvl.append(cFunc_j)
1459 # Construct pseudo-inverse value functions over bLvl for each pLvl
1460 v_mid = u(cLvl) + EndOfPrd_v # value of reaching consumption phase
1461 vNvrsFuncMid_by_pLvl = []
1462 for j in range(pLvlCount):
1463 b_cnst = np.linspace(0.001, 0.95, 10) * bLvl[0, j] # constrained wealth levels
1464 c_cnst = b_cnst
1465 v_cnst = u(c_cnst) + EndOfPrd_v[0, j]
1466 b_temp = np.concatenate([b_cnst, bLvl[:, j]])
1467 v_temp = np.concatenate([v_cnst, v_mid[:, j]])
1468 vNvrs_temp = n(v_temp)
1469 vNvrsFunc_j = LinearInterp(
1470 np.insert(b_temp, 0, 0.0), np.insert(vNvrs_temp, 0, 0.0)
1471 )
1472 vNvrsFuncMid_by_pLvl.append(vNvrsFunc_j)
1474 # Make a grid of (log) medical expenses (and probs), cross it with (mLvl,pLvl)
1475 MedCostBaseGrid = np.linspace(MedCostBot, MedCostTop, MedCostCount)
1476 MedCostLogGrid = MedCostLogMean + MedCostBaseGrid * MedCostLogStd
1477 MedCostGrid = np.exp(MedCostLogGrid)
1478 mLvl_base = np.dot(
1479 np.reshape(mNrmGrid, (mNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvlCount))
1480 )
1481 mLvl = np.reshape(mLvl_base, (mNrmGrid.size, pLvlCount, 1))
1482 bLvl_if_care = mLvl - np.reshape(MedCostGrid, (1, 1, MedCostCount))
1483 bLvl_if_not = mLvl_base
1485 # Calculate mean (log) utility shock for each MedCost gridpoint, and conditional stdev
1486 MedShkLog_cond_mean = MedShkLogMean + MedCorr * MedShkLogStd * MedCostBaseGrid
1487 MedShkLog_cond_mean = np.reshape(MedShkLog_cond_mean, (1, MedCostCount))
1488 MedShkLog_cond_std = np.sqrt(MedShkLogStd**2 * (1.0 - MedCorr**2))
1489 MedShk_cond_mean = np.exp(MedShkLog_cond_mean + 0.5 * MedShkLog_cond_std**2)
1491 # Initialize (marginal) value function arrays over (mLvl,pLvl,MedCost)
1492 v_at_Dcsn = np.empty_like(bLvl_if_care)
1493 vP_at_Dcsn = np.empty_like(bLvl_if_care)
1494 care_prob_array = np.empty_like(bLvl_if_care)
1495 for j in range(pLvlCount):
1496 # Evaluate value function for (bLvl,pLvl_j), including MedCost=0
1497 v_if_care = u(vNvrsFuncMid_by_pLvl[j](bLvl_if_care[:, j, :]))
1498 v_if_not = np.reshape(
1499 u(vNvrsFuncMid_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1500 )
1501 cant_pay = bLvl_if_care[:, j, :] <= 0.0
1502 v_if_care[cant_pay] = -np.inf
1504 # Find value difference at each gridpoint, convert to MedShk stdev; find prob of care
1505 v_diff = v_if_not - v_if_care
1506 log_v_diff = np.log(v_diff)
1507 crit_stdev = (log_v_diff - MedShkLog_cond_mean) / MedShkLog_cond_std
1508 prob_no_care = norm.cdf(crit_stdev)
1509 prob_get_care = 1.0 - prob_no_care
1510 care_prob_array[:, j, :] = prob_get_care
1512 # Calculate expected MedShk conditional on not getting medical care
1513 crit_z = crit_stdev - MedShkLog_cond_std
1514 MedShk_no_care_cond_mean = 0.5 * MedShk_cond_mean * erfc(crit_z) / prob_no_care
1516 # Compute expected (marginal) value over MedShk for each (mLvl,pLvl_j,MedCost)
1517 v_if_care[cant_pay] = 0.0
1518 v_at_Dcsn[:, j, :] = (
1519 prob_no_care * (v_if_not - MedShk_no_care_cond_mean)
1520 + prob_get_care * v_if_care
1521 )
1522 vP_if_care = uP(cFunc_by_pLvl[j](bLvl_if_care[:, j, :]))
1523 vP_if_not = np.reshape(
1524 uP(cFunc_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1525 )
1526 vP_if_care[cant_pay] = 0.0
1527 MedShk_rate_of_change = (
1528 norm.pdf(crit_stdev) * (vP_if_care - vP_if_not) * MedShk_no_care_cond_mean
1529 )
1530 vP_at_Dcsn[:, j, :] = (
1531 prob_no_care * vP_if_not
1532 + prob_get_care * vP_if_care
1533 + MedShk_rate_of_change
1534 )
1536 # Compute expected (marginal) value over MedCost for each (mLvl,pLvl)
1537 temp_grid = np.linspace(MedCostBot, MedCostTop, MedCostCount)
1538 MedCost_pmv = norm.pdf(temp_grid)
1539 MedCost_pmv /= np.sum(MedCost_pmv)
1540 MedCost_probs = np.reshape(MedCost_pmv, (1, 1, MedCostCount))
1541 v_before_shk = np.sum(v_at_Dcsn * MedCost_probs, axis=2)
1542 vP_before_shk = np.sum(vP_at_Dcsn * MedCost_probs, axis=2)
1543 vNvrs_before_shk = n(v_before_shk)
1544 vPnvrs_before_shk = CRRAutilityP_inv(vP_before_shk, CRRA)
1546 # Compute expected medical expenses at each state space point
1547 ExpCare_all = care_prob_array * np.reshape(MedCostGrid, (1, 1, MedCostCount))
1548 ExpCare = np.sum(ExpCare_all * MedCost_probs, axis=2)
1549 ProbCare = np.sum(care_prob_array * MedCost_probs, axis=2)
1550 ExpCareFunc_by_pLvl = []
1551 CareProbFunc_by_pLvl = []
1552 for j in range(pLvlCount):
1553 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1554 EC_temp = np.insert(ExpCare[:, j], 0, 0.0)
1555 prob_temp = np.insert(ProbCare[:, j], 0, 0.0)
1556 ExpCareFunc_by_pLvl.append(LinearInterp(m_temp, EC_temp))
1557 CareProbFunc_by_pLvl.append(LinearInterp(m_temp, prob_temp))
1558 ExpCareFunc = LinearInterpOnInterp1D(ExpCareFunc_by_pLvl, pLvl)
1559 ProbCareFunc = LinearInterpOnInterp1D(CareProbFunc_by_pLvl, pLvl)
1561 # Fixing kLvlGrid, compute expected (marginal) value over TranShk for each (kLvl,pLvl)
1562 v_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1563 vP_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1564 for j in range(pLvlCount):
1565 p = pLvl[j]
1567 # Make (marginal) value functions over mLvl for this pLvl
1568 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1569 vNvrs_temp = np.insert(vNvrs_before_shk[:, j], 0, 0.0)
1570 vPnvrs_temp = np.insert(vPnvrs_before_shk[:, j], 0, 0.0)
1571 vNvrsFunc_temp = LinearInterp(m_temp, vNvrs_temp)
1572 vPnvrsFunc_temp = LinearInterp(m_temp, vPnvrs_temp)
1573 vFunc_temp = lambda x: u(vNvrsFunc_temp(x))
1574 vPfunc_temp = lambda x: uP(vPnvrsFunc_temp(x))
1576 # Compute expectation over TranShkDstn
1577 v = lambda TranShk, kLvl: vFunc_temp(kLvl + TranShk * p)
1578 vP = lambda TranShk, kLvl: vPfunc_temp(kLvl + TranShk * p)
1579 v_by_kLvl_and_pLvl[:, j] = expected(v, TranShkDstn, args=(kLvlGrid,))
1580 vP_by_kLvl_and_pLvl[:, j] = expected(vP, TranShkDstn, args=(kLvlGrid,))
1582 # Compute expectation over persistent shocks by using pLvlMrkvArray
1583 v_arvl = np.dot(v_by_kLvl_and_pLvl, pLogMrkvArray.T)
1584 vP_arvl = np.dot(vP_by_kLvl_and_pLvl, pLogMrkvArray.T)
1585 vNvrs_arvl = n(v_arvl)
1586 vPnvrs_arvl = CRRAutilityP_inv(vP_arvl, CRRA)
1588 # Construct "arrival" (marginal) value function by pLvl
1589 vFuncArvl_by_pLvl = []
1590 vPfuncArvl_by_pLvl = []
1591 for j in range(pLvlCount):
1592 vNvrsFunc_temp = LinearInterp(kLvlGrid, vNvrs_arvl[:, j])
1593 vPnvrsFunc_temp = LinearInterp(kLvlGrid, vPnvrs_arvl[:, j])
1594 vFuncArvl_by_pLvl.append(ValueFuncCRRA(vNvrsFunc_temp, CRRA))
1595 vPfuncArvl_by_pLvl.append(MargValueFuncCRRA(vPnvrsFunc_temp, CRRA))
1597 # Gather elements and return the solution object
1598 solution_now = ConsMedExtMargSolution(
1599 vFunc_by_pLvl=vFuncArvl_by_pLvl,
1600 vPfunc_by_pLvl=vPfuncArvl_by_pLvl,
1601 cFunc_by_pLvl=cFunc_by_pLvl,
1602 vNvrsFuncMid_by_pLvl=vNvrsFuncMid_by_pLvl,
1603 pLvl=pLvl,
1604 CRRA=CRRA,
1605 ExpMedFunc=ExpCareFunc,
1606 CareProbFunc=ProbCareFunc,
1607 )
1608 return solution_now
1611###############################################################################
1614# Define a dictionary of constructors for the extensive margin medical spending model
1615med_ext_marg_constructors = {
1616 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
1617 "IncomeProcessDict": make_persistent_income_process_dict,
1618 "pLogGrid": get_it_from("IncomeProcessDict"),
1619 "pLvlMean": get_it_from("IncomeProcessDict"),
1620 "pLogMrkvArray": get_it_from("IncomeProcessDict"),
1621 "IncShkDstn": construct_lognormal_income_process_unemployment,
1622 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1623 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1624 "BeqFac": get_it_from("BeqParamDict"),
1625 "BeqShift": get_it_from("BeqParamDict"),
1626 "BeqParamDict": reformat_bequest_motive,
1627 "aNrmGrid": make_assets_grid,
1628 "mNrmGrid": make_market_resources_grid,
1629 "kLvlGrid": make_capital_grid,
1630 "solution_terminal": make_MedExtMarg_solution_terminal,
1631 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1632 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1633 "MedShockDstn": make_continuous_MedShockDstn,
1634}
1636# Make a dictionary with default bequest motive parameters
1637default_BeqParam_dict = {
1638 "BeqMPC": 0.1, # Hypothetical "MPC at death"
1639 "BeqInt": 1.0, # Intercept term for hypothetical "consumption function at death"
1640}
1642# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1643default_kNrmInitDstn_params_ExtMarg = {
1644 "kLogInitMean": -6.0, # Mean of log initial capital
1645 "kLogInitStd": 1.0, # Stdev of log initial capital
1646 "kNrmInitCount": 15, # Number of points in initial capital discretization
1647}
1649# Default parameters to make IncomeProcessDict using make_persistent_income_process_dict;
1650# some of these are used by construct_lognormal_income_process_unemployment as well
1651default_IncomeProcess_params = {
1652 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1653 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1654 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1655 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1656 "UnempPrb": 0.05, # Probability of unemployment while working
1657 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1658 "T_retire": 0, # Period of retirement (0 --> no retirement)
1659 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1660 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1661 "pLogInitMean": 0.0, # Mean of log initial permanent income
1662 "pLogInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
1663 "pLvlInitCount": 25, # Number of discrete nodes in initial permanent income level dstn
1664 "PermGroFac": [1.0], # Permanent income growth factor
1665 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
1666 "pLogCount": 45, # Number of points in persistent income grid each period
1667 "pLogRange": 3.5, # Upper/lower bound of persistent income, in unconditional standard deviations
1668}
1670# Default parameters to make aNrmGrid using make_assets_grid
1671default_aNrmGrid_params = {
1672 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1673 "aXtraMax": 40.0, # Maximum end-of-period "assets above minimum" value
1674 "aXtraNestFac": 2, # Exponential nesting factor for aXtraGrid
1675 "aXtraCount": 96, # Number of points in the grid of "assets above minimum"
1676 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
1677}
1679# Default parameters to make mLvlGrid using make_market_resources_grid
1680default_mNrmGrid_params = {
1681 "mNrmMin": 0.001,
1682 "mNrmMax": 40.0,
1683 "mNrmNestFac": 2,
1684 "mNrmCount": 72,
1685 "mNrmExtra": None,
1686}
1688# Default parameters to make kLvlGrid using make_capital_grid
1689default_kLvlGrid_params = {
1690 "kLvlMin": 0.0,
1691 "kLvlMax": 200,
1692 "kLvlCount": 250,
1693 "kLvlOrder": 2.0,
1694}
1696# Default "basic" parameters
1697med_ext_marg_basic_params = {
1698 "constructors": med_ext_marg_constructors,
1699 "cycles": 1, # Lifecycle by default
1700 "T_cycle": 1, # Number of periods in the default sequence
1701 "T_age": None,
1702 "AgentCount": 10000, # Number of agents to simulate
1703 "DiscFac": 0.96, # intertemporal discount factor
1704 "CRRA": 1.5, # coefficient of relative risk aversion
1705 "Rfree": [1.02], # risk free interest factor
1706 "LivPrb": [0.99], # survival probability
1707 "MedCostBot": -3.1, # lower bound of medical cost distribution, in stdevs
1708 "MedCostTop": 5.2, # upper bound of medical cost distribution, in stdevs
1709 "MedCostCount": 84, # number of nodes in medical cost discretization
1710 "MedShkLogMean": [-2.0], # mean of log utility shocks
1711 "MedShkLogStd": [1.5], # standard deviation of log utility shocks
1712 "MedCostLogMean": [-1.0], # mean of log medical expenses
1713 "MedCostLogStd": [1.0], # standard deviation of log medical expenses
1714 "MedCorr": [0.3], # correlation coefficient between utility shock and expenses
1715 "PermGroFacAgg": 1.0, # Aggregate productivity growth rate (leave at 1)
1716 "NewbornTransShk": True, # Whether "newborns" have a transitory income shock
1717}
1719# Combine the dictionaries into a single default dictionary
1720init_med_ext_marg = med_ext_marg_basic_params.copy()
1721init_med_ext_marg.update(default_IncomeProcess_params)
1722init_med_ext_marg.update(default_aNrmGrid_params)
1723init_med_ext_marg.update(default_mNrmGrid_params)
1724init_med_ext_marg.update(default_kLvlGrid_params)
1725init_med_ext_marg.update(default_kNrmInitDstn_params_ExtMarg)
1726init_med_ext_marg.update(default_BeqParam_dict)
1729class MedExtMargConsumerType(PersistentShockConsumerType):
1730 r"""
1731 Class for representing agents in the extensive margin medical expense model.
1732 Such agents have labor income dynamics identical to the "general income process"
1733 model (permanent income is not normalized out), and also experience a medical
1734 shock with two components: medical cost and utility loss. They face a binary
1735 choice of whether to pay the cost or suffer the loss, then make a consumption-
1736 saving decision as normal. To simplify the computation, the joint distribution
1737 of medical shocks is specified as bivariate lognormal. This can be loosened to
1738 accommodate insurance contracts as mappings from total to out-of-pocket expenses.
1739 Can also be extended to include a health process.
1741 .. math::
1742 \begin{eqnarray*}
1743 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}], \\
1744 A_t &=& M_t - C_t - D_t med_t, \\
1745 A_t/ &\geq& 0, \\
1746 D_t &\in& \{0,1\}, \\
1747 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
1748 Y_{t+1} &=& P_{t+1} \theta_{t+1}
1749 M_{t+1} &=& R A_t + Y_{t+1}, \\
1750 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
1751 (med_t,\eta_t) &~\sim& G_t,\\
1752 U_t(C) &=& \frac{C^{1-\rho}}{1-\rho}.
1753 \end{eqnarray*}
1754 """
1756 default_ = {
1757 "params": init_med_ext_marg,
1758 "solver": solve_one_period_ConsMedExtMarg,
1759 "model": "ConsExtMargMed.yaml",
1760 "track_vars": ["aLvl", "cLvl", "Med", "mLvl", "pLvl"],
1761 }
1763 time_vary_ = [
1764 "Rfree",
1765 "LivPrb",
1766 "MedShkLogMean",
1767 "MedShkLogStd",
1768 "MedCostLogMean",
1769 "MedCostLogStd",
1770 "MedCorr",
1771 "pLogGrid",
1772 "pLvlMean",
1773 "TranShkDstn",
1774 "pLogMrkvArray",
1775 "pLvlNextFunc",
1776 "IncShkDstn",
1777 "MedShockDstn",
1778 ]
1779 time_inv_ = [
1780 "DiscFac",
1781 "CRRA",
1782 "BeqFac",
1783 "BeqShift",
1784 "MedCostBot",
1785 "MedCostTop",
1786 "MedCostCount",
1787 "aNrmGrid",
1788 "mNrmGrid",
1789 "kLvlGrid",
1790 ]
1791 shock_vars = ["PermShk", "TranShk", "MedShk", "MedCost"]
1793 def get_shocks(self):
1794 """
1795 Gets permanent and transitory income shocks for this period as well as
1796 medical need and cost shocks.
1797 """
1798 # Get permanent and transitory income shocks
1799 PersistentShockConsumerType.get_shocks(self)
1801 # Initialize medical shock array and cost of care array
1802 MedShkNow = np.zeros(self.AgentCount)
1803 MedCostNow = np.zeros(self.AgentCount)
1805 # Get shocks for each period of the cycle
1806 for t in range(self.T_cycle):
1807 these = t == self.t_cycle
1808 if np.any(these):
1809 N = np.sum(these)
1810 dstn_t = self.MedShockDstn[t]
1811 draws_t = dstn_t.draw(N)
1812 MedCostNow[these] = draws_t[0, :]
1813 MedShkNow[these] = draws_t[1, :]
1814 self.shocks["MedShk"] = MedShkNow
1815 self.shocks["MedCost"] = MedCostNow
1817 def get_controls(self):
1818 """
1819 Finds consumption for each agent, along with whether or not they get care.
1820 """
1821 # Initialize output
1822 cLvlNow = np.empty(self.AgentCount)
1823 CareNow = np.zeros(self.AgentCount, dtype=bool)
1825 # Get states and shocks
1826 mLvl = self.state_now["mLvl"]
1827 pLvl = self.state_now["pLvl"]
1828 MedCost = self.shocks["MedCost"]
1829 MedShk = self.shocks["MedShk"]
1831 # Find remaining resources with and without care
1832 bLvl_no_care = mLvl
1833 bLvl_with_care = mLvl - MedCost
1835 # Get controls for each period of the cycle
1836 for t in range(self.T_cycle):
1837 these = t == self.t_cycle
1838 if np.any(these):
1839 vFunc_t = self.solution[t].vFuncMid
1840 cFunc_t = self.solution[t].cFunc
1842 v_no_care = vFunc_t(bLvl_no_care[these], pLvl[these]) - MedShk[these]
1843 v_if_care = vFunc_t(bLvl_with_care[these], pLvl[these])
1844 get_care = v_if_care > v_no_care
1846 b_temp = bLvl_no_care[these]
1847 b_temp[get_care] = bLvl_with_care[get_care]
1848 cLvlNow[these] = cFunc_t(b_temp, pLvl[these])
1849 CareNow[these] = get_care
1851 # Store the results
1852 self.controls["cLvl"] = cLvlNow
1853 self.controls["Care"] = CareNow
1855 def get_poststates(self):
1856 """
1857 Calculates end-of-period assets for each consumer of this type.
1858 """
1859 self.state_now["Med"] = self.shocks["MedCost"] * self.controls["Care"]
1860 self.state_now["aLvl"] = (
1861 self.state_now["mLvl"] - self.controls["cLvl"] - self.state_now["Med"]
1862 )
1863 # Move now to prev
1864 AgentType.get_poststates(self)