Coverage for HARK/ConsumptionSaving/ConsMedModel.py: 92%
658 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +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
10from scipy.optimize import brentq
12from HARK import AgentType
13from HARK.Calibration.Income.IncomeProcesses import (
14 construct_lognormal_income_process_unemployment,
15 get_PermShkDstn_from_IncShkDstn,
16 get_TranShkDstn_from_IncShkDstn,
17 make_AR1_style_pLvlNextFunc,
18 make_pLvlGrid_by_simulation,
19 make_basic_pLvlPctiles,
20 make_persistent_income_process_dict,
21)
22from HARK.ConsumptionSaving.ConsIndShockModel import (
23 make_lognormal_kNrm_init_dstn,
24 make_lognormal_pLvl_init_dstn,
25)
26from HARK.ConsumptionSaving.ConsGenIncProcessModel import (
27 PersistentShockConsumerType,
28 VariableLowerBoundFunc2D,
29)
30from HARK.ConsumptionSaving.ConsIndShockModel import ConsumerSolution
31from HARK.distributions import (
32 Lognormal,
33 MultivariateLognormal,
34 add_discrete_outcome_constant_mean,
35 expected,
36)
37from HARK.interpolation import (
38 BilinearInterp,
39 BilinearInterpOnInterp1D,
40 ConstantFunction,
41 CubicInterp,
42 LinearInterp,
43 LinearInterpOnInterp1D,
44 LowerEnvelope3D,
45 MargMargValueFuncCRRA,
46 MargValueFuncCRRA,
47 TrilinearInterp,
48 UpperEnvelope,
49 ValueFuncCRRA,
50 VariableLowerBoundFunc3D,
51)
52from HARK.metric import MetricObject
53from HARK.rewards import (
54 CRRAutility,
55 CRRAutilityP,
56 CRRAutility_inv,
57 CRRAutility_invP,
58 CRRAutilityP_inv,
59 CRRAutilityPP,
60 UtilityFuncCRRA,
61)
62from HARK.utilities import NullFunc, make_grid_exp_mult, make_assets_grid, get_it_from
64__all__ = [
65 "MedShockPolicyFunc",
66 "cThruXfunc",
67 "MedThruXfunc",
68 "MedShockConsumerType",
69 "make_lognormal_MedShkDstn",
70]
72utility_inv = CRRAutility_inv
73utilityP_inv = CRRAutilityP_inv
74utility = CRRAutility
75utility_invP = CRRAutility_invP
76utilityPP = CRRAutilityPP
79class MedShockPolicyFunc(MetricObject):
80 """
81 Class for representing the policy function in the medical shocks model: opt-
82 imal consumption and medical care for given market resources, permanent income,
83 and medical need shock. Always obeys Con + MedPrice*Med = optimal spending.
85 Parameters
86 ----------
87 xFunc : function
88 Optimal total spending as a function of market resources, permanent
89 income, and the medical need shock.
90 xLvlGrid : np.array
91 1D array of total expenditure levels.
92 MedShkGrid : np.array
93 1D array of medical shocks.
94 MedPrice : float
95 Relative price of a unit of medical care.
96 CRRAcon : float
97 Coefficient of relative risk aversion for consumption.
98 CRRAmed : float
99 Coefficient of relative risk aversion for medical care.
100 xLvlCubicBool : boolean
101 Indicator for whether cubic spline interpolation (rather than linear)
102 should be used in the xLvl dimension.
103 MedShkCubicBool : boolean
104 Indicator for whether bicubic interpolation should be used; only
105 operative when xLvlCubicBool=True.
106 """
108 distance_criteria = ["xFunc", "cFunc", "MedPrice"]
110 def __init__(
111 self,
112 xFunc,
113 xLvlGrid,
114 MedShkGrid,
115 MedPrice,
116 CRRAcon,
117 CRRAmed,
118 xLvlCubicBool=False,
119 MedShkCubicBool=False,
120 ):
121 # Store some of the inputs in self
122 self.MedPrice = MedPrice
123 self.xFunc = xFunc
125 # Calculate optimal consumption at each combination of mLvl and MedShk.
126 cLvlGrid = np.zeros(
127 (xLvlGrid.size, MedShkGrid.size)
128 ) # Initialize consumption grid
129 for i in range(xLvlGrid.size):
130 xLvl = xLvlGrid[i]
131 for j in range(MedShkGrid.size):
132 MedShk = MedShkGrid[j]
133 if xLvl == 0: # Zero consumption when mLvl = 0
134 cLvl = 0.0
135 elif MedShk == 0: # All consumption when MedShk = 0
136 cLvl = xLvl
137 else:
139 def optMedZeroFunc(c):
140 return (MedShk / MedPrice) ** (-1.0 / CRRAcon) * (
141 (xLvl - c) / MedPrice
142 ) ** (CRRAmed / CRRAcon) - c
144 # Find solution to FOC
145 cLvl = brentq(optMedZeroFunc, 0.0, xLvl)
146 cLvlGrid[i, j] = cLvl
148 # Construct the consumption function and medical care function
149 if xLvlCubicBool:
150 if MedShkCubicBool:
151 raise NotImplementedError()("Bicubic interpolation not yet implemented")
152 else:
153 xLvlGrid_tiled = np.tile(
154 np.reshape(xLvlGrid, (xLvlGrid.size, 1)), (1, MedShkGrid.size)
155 )
156 MedShkGrid_tiled = np.tile(
157 np.reshape(MedShkGrid, (1, MedShkGrid.size)), (xLvlGrid.size, 1)
158 )
159 dfdx = (
160 (CRRAmed / (CRRAcon * MedPrice))
161 * (MedShkGrid_tiled / MedPrice) ** (-1.0 / CRRAcon)
162 * ((xLvlGrid_tiled - cLvlGrid) / MedPrice)
163 ** (CRRAmed / CRRAcon - 1.0)
164 )
165 dcdx = dfdx / (dfdx + 1.0)
166 # approximation; function goes crazy otherwise
167 dcdx[0, :] = dcdx[1, :]
168 dcdx[:, 0] = 1.0 # no Med when MedShk=0, so all x is c
169 cFromxFunc_by_MedShk = []
170 for j in range(MedShkGrid.size):
171 cFromxFunc_by_MedShk.append(
172 CubicInterp(xLvlGrid, cLvlGrid[:, j], dcdx[:, j])
173 )
174 cFunc = LinearInterpOnInterp1D(cFromxFunc_by_MedShk, MedShkGrid)
175 else:
176 cFunc = BilinearInterp(cLvlGrid, xLvlGrid, MedShkGrid)
177 self.cFunc = cFunc
179 def __call__(self, mLvl, pLvl, MedShk):
180 """
181 Evaluate optimal consumption and medical care at given levels of market
182 resources, permanent income, and medical need shocks.
184 Parameters
185 ----------
186 mLvl : np.array
187 Market resource levels.
188 pLvl : np.array
189 Permanent income levels; should be same size as mLvl.
190 MedShk : np.array
191 Medical need shocks; should be same size as mLvl.
193 Returns
194 -------
195 cLvl : np.array
196 Optimal consumption for each point in (xLvl,MedShk).
197 Med : np.array
198 Optimal medical care for each point in (xLvl,MedShk).
199 """
200 xLvl = self.xFunc(mLvl, pLvl, MedShk)
201 cLvl = self.cFunc(xLvl, MedShk)
202 Med = (xLvl - cLvl) / self.MedPrice
203 return cLvl, Med
205 def derivativeX(self, mLvl, pLvl, MedShk):
206 """
207 Evaluate the derivative of consumption and medical care with respect to
208 market resources at given levels of market resources, permanent income,
209 and medical need shocks.
211 Parameters
212 ----------
213 mLvl : np.array
214 Market resource levels.
215 pLvl : np.array
216 Permanent income levels; should be same size as mLvl.
217 MedShk : np.array
218 Medical need shocks; should be same size as mLvl.
220 Returns
221 -------
222 dcdm : np.array
223 Derivative of consumption with respect to market resources for each
224 point in (xLvl,MedShk).
225 dMeddm : np.array
226 Derivative of medical care with respect to market resources for each
227 point in (xLvl,MedShk).
228 """
229 xLvl = self.xFunc(mLvl, pLvl, MedShk)
230 dxdm = self.xFunc.derivativeX(mLvl, pLvl, MedShk)
231 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
232 dcdm = dxdm * dcdx
233 dMeddm = (dxdm - dcdm) / self.MedPrice
234 return dcdm, dMeddm
236 def derivativeY(self, mLvl, pLvl, MedShk):
237 """
238 Evaluate the derivative of consumption and medical care with respect to
239 permanent income at given levels of market resources, permanent income,
240 and medical need shocks.
242 Parameters
243 ----------
244 mLvl : np.array
245 Market resource levels.
246 pLvl : np.array
247 Permanent income levels; should be same size as mLvl.
248 MedShk : np.array
249 Medical need shocks; should be same size as mLvl.
251 Returns
252 -------
253 dcdp : np.array
254 Derivative of consumption with respect to permanent income for each
255 point in (xLvl,MedShk).
256 dMeddp : np.array
257 Derivative of medical care with respect to permanent income for each
258 point in (xLvl,MedShk).
259 """
260 xLvl = self.xFunc(mLvl, pLvl, MedShk)
261 dxdp = self.xFunc.derivativeY(mLvl, pLvl, MedShk)
262 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
263 dcdp = dxdp * dcdx
264 dMeddp = (dxdp - dcdp) / self.MedPrice
265 return dcdp, dMeddp
267 def derivativeZ(self, mLvl, pLvl, MedShk):
268 """
269 Evaluate the derivative of consumption and medical care with respect to
270 medical need shock at given levels of market resources, permanent income,
271 and medical need shocks.
273 Parameters
274 ----------
275 mLvl : np.array
276 Market resource levels.
277 pLvl : np.array
278 Permanent income levels; should be same size as mLvl.
279 MedShk : np.array
280 Medical need shocks; should be same size as mLvl.
282 Returns
283 -------
284 dcdShk : np.array
285 Derivative of consumption with respect to medical need for each
286 point in (xLvl,MedShk).
287 dMeddShk : np.array
288 Derivative of medical care with respect to medical need for each
289 point in (xLvl,MedShk).
290 """
291 xLvl = self.xFunc(mLvl, pLvl, MedShk)
292 dxdShk = self.xFunc.derivativeZ(mLvl, pLvl, MedShk)
293 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
294 dcdShk = dxdShk * dcdx + self.cFunc.derivativeY(xLvl, MedShk)
295 dMeddShk = (dxdShk - dcdShk) / self.MedPrice
296 return dcdShk, dMeddShk
299class cThruXfunc(MetricObject):
300 """
301 Class for representing consumption function derived from total expenditure
302 and consumption.
304 Parameters
305 ----------
306 xFunc : function
307 Optimal total spending as a function of market resources, permanent
308 income, and the medical need shock.
309 cFunc : function
310 Optimal consumption as a function of total spending and the medical
311 need shock.
312 """
314 distance_criteria = ["xFunc", "cFunc"]
316 def __init__(self, xFunc, cFunc):
317 self.xFunc = xFunc
318 self.cFunc = cFunc
320 def __call__(self, mLvl, pLvl, MedShk):
321 """
322 Evaluate optimal consumption at given levels of market resources, perma-
323 nent income, and medical need shocks.
325 Parameters
326 ----------
327 mLvl : np.array
328 Market resource levels.
329 pLvl : np.array
330 Permanent income levels; should be same size as mLvl.
331 MedShk : np.array
332 Medical need shocks; should be same size as mLvl.
334 Returns
335 -------
336 cLvl : np.array
337 Optimal consumption for each point in (xLvl,MedShk).
338 """
339 xLvl = self.xFunc(mLvl, pLvl, MedShk)
340 cLvl = self.cFunc(xLvl, MedShk)
341 return cLvl
343 def derivativeX(self, mLvl, pLvl, MedShk):
344 """
345 Evaluate the derivative of consumption with respect to market resources
346 at given levels of market resources, permanent income, and medical need
347 shocks.
349 Parameters
350 ----------
351 mLvl : np.array
352 Market resource levels.
353 pLvl : np.array
354 Permanent income levels; should be same size as mLvl.
355 MedShk : np.array
356 Medical need shocks; should be same size as mLvl.
358 Returns
359 -------
360 dcdm : np.array
361 Derivative of consumption with respect to market resources for each
362 point in (xLvl,MedShk).
363 """
364 xLvl = self.xFunc(mLvl, pLvl, MedShk)
365 dxdm = self.xFunc.derivativeX(mLvl, pLvl, MedShk)
366 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
367 dcdm = dxdm * dcdx
368 return dcdm
370 def derivativeY(self, mLvl, pLvl, MedShk):
371 """
372 Evaluate the derivative of consumption and medical care with respect to
373 permanent income at given levels of market resources, permanent income,
374 and medical need shocks.
376 Parameters
377 ----------
378 mLvl : np.array
379 Market resource levels.
380 pLvl : np.array
381 Permanent income levels; should be same size as mLvl.
382 MedShk : np.array
383 Medical need shocks; should be same size as mLvl.
385 Returns
386 -------
387 dcdp : np.array
388 Derivative of consumption with respect to permanent income for each
389 point in (xLvl,MedShk).
390 """
391 xLvl = self.xFunc(mLvl, pLvl, MedShk)
392 dxdp = self.xFunc.derivativeY(mLvl, pLvl, MedShk)
393 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
394 dcdp = dxdp * dcdx
395 return dcdp
397 def derivativeZ(self, mLvl, pLvl, MedShk):
398 """
399 Evaluate the derivative of consumption and medical care with respect to
400 medical need shock at given levels of market resources, permanent income,
401 and medical need shocks.
403 Parameters
404 ----------
405 mLvl : np.array
406 Market resource levels.
407 pLvl : np.array
408 Permanent income levels; should be same size as mLvl.
409 MedShk : np.array
410 Medical need shocks; should be same size as mLvl.
412 Returns
413 -------
414 dcdShk : np.array
415 Derivative of consumption with respect to medical need for each
416 point in (xLvl,MedShk).
417 """
418 xLvl = self.xFunc(mLvl, pLvl, MedShk)
419 dxdShk = self.xFunc.derivativeZ(mLvl, pLvl, MedShk)
420 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
421 dcdShk = dxdShk * dcdx + self.cFunc.derivativeY(xLvl, MedShk)
422 return dcdShk
425class MedThruXfunc(MetricObject):
426 """
427 Class for representing medical care function derived from total expenditure
428 and consumption.
430 Parameters
431 ----------
432 xFunc : function
433 Optimal total spending as a function of market resources, permanent
434 income, and the medical need shock.
435 cFunc : function
436 Optimal consumption as a function of total spending and the medical
437 need shock.
438 MedPrice : float
439 Relative price of a unit of medical care.
440 """
442 distance_criteria = ["xFunc", "cFunc", "MedPrice"]
444 def __init__(self, xFunc, cFunc, MedPrice):
445 self.xFunc = xFunc
446 self.cFunc = cFunc
447 self.MedPrice = MedPrice
449 def __call__(self, mLvl, pLvl, MedShk):
450 """
451 Evaluate optimal medical care at given levels of market resources,
452 permanent income, and medical need shock.
454 Parameters
455 ----------
456 mLvl : np.array
457 Market resource levels.
458 pLvl : np.array
459 Permanent income levels; should be same size as mLvl.
460 MedShk : np.array
461 Medical need shocks; should be same size as mLvl.
463 Returns
464 -------
465 Med : np.array
466 Optimal medical care for each point in (xLvl,MedShk).
467 """
468 xLvl = self.xFunc(mLvl, pLvl, MedShk)
469 Med = (xLvl - self.cFunc(xLvl, MedShk)) / self.MedPrice
470 return Med
472 def derivativeX(self, mLvl, pLvl, MedShk):
473 """
474 Evaluate the derivative of consumption and medical care with respect to
475 market resources at given levels of market resources, permanent income,
476 and medical need shocks.
478 Parameters
479 ----------
480 mLvl : np.array
481 Market resource levels.
482 pLvl : np.array
483 Permanent income levels; should be same size as mLvl.
484 MedShk : np.array
485 Medical need shocks; should be same size as mLvl.
487 Returns
488 -------
489 dcdm : np.array
490 Derivative of consumption with respect to market resources for each
491 point in (xLvl,MedShk).
492 dMeddm : np.array
493 Derivative of medical care with respect to market resources for each
494 point in (xLvl,MedShk).
495 """
496 xLvl = self.xFunc(mLvl, pLvl, MedShk)
497 dxdm = self.xFunc.derivativeX(mLvl, pLvl, MedShk)
498 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
499 dcdm = dxdm * dcdx
500 dMeddm = (dxdm - dcdm) / self.MedPrice
501 return dcdm, dMeddm
503 def derivativeY(self, mLvl, pLvl, MedShk):
504 """
505 Evaluate the derivative of medical care with respect to permanent income
506 at given levels of market resources, permanent income, and medical need
507 shocks.
509 Parameters
510 ----------
511 mLvl : np.array
512 Market resource levels.
513 pLvl : np.array
514 Permanent income levels; should be same size as mLvl.
515 MedShk : np.array
516 Medical need shocks; should be same size as mLvl.
518 Returns
519 -------
520 dMeddp : np.array
521 Derivative of medical care with respect to permanent income for each
522 point in (xLvl,MedShk).
523 """
524 xLvl = self.xFunc(mLvl, pLvl, MedShk)
525 dxdp = self.xFunc.derivativeY(mLvl, pLvl, MedShk)
526 dMeddp = (dxdp - dxdp * self.cFunc.derivativeX(xLvl, MedShk)) / self.MedPrice
527 return dMeddp
529 def derivativeZ(self, mLvl, pLvl, MedShk):
530 """
531 Evaluate the derivative of medical care with respect to medical need
532 shock at given levels of market resources, permanent income, and medical
533 need shocks.
535 Parameters
536 ----------
537 mLvl : np.array
538 Market resource levels.
539 pLvl : np.array
540 Permanent income levels; should be same size as mLvl.
541 MedShk : np.array
542 Medical need shocks; should be same size as mLvl.
544 Returns
545 -------
546 dMeddShk : np.array
547 Derivative of medical care with respect to medical need for each
548 point in (xLvl,MedShk).
549 """
550 xLvl = self.xFunc(mLvl, pLvl, MedShk)
551 dxdShk = self.xFunc.derivativeZ(mLvl, pLvl, MedShk)
552 dcdx = self.cFunc.derivativeX(xLvl, MedShk)
553 dcdShk = dxdShk * dcdx + self.cFunc.derivativeY(xLvl, MedShk)
554 dMeddShk = (dxdShk - dcdShk) / self.MedPrice
555 return dMeddShk
558def make_market_resources_grid(mNrmMin, mNrmMax, mNrmNestFac, mNrmCount, mNrmExtra):
559 """
560 Constructor for mNrmGrid that aliases make_assets_grid.
561 """
562 return make_assets_grid(mNrmMin, mNrmMax, mNrmCount, mNrmExtra, mNrmNestFac)
565def make_capital_grid(kLvlMin, kLvlMax, kLvlCount, kLvlOrder):
566 """
567 Constructor for kLvlGrid, using a simple "invertible" format.
568 """
569 base_grid = np.linspace(0.0, 1.0, kLvlCount) ** kLvlOrder
570 kLvlGrid = (kLvlMax - kLvlMin) * base_grid + kLvlMin
571 return kLvlGrid
574def reformat_bequest_motive(BeqMPC, BeqInt, CRRA):
575 """
576 Reformats interpretable bequest motive parameters (terminal intercept and MPC)
577 into parameters that are easily useable in math (shifter and scaler).
578 """
579 BeqParamDict = {
580 "BeqFac": BeqMPC ** (-CRRA),
581 "BeqShift": BeqInt / BeqMPC,
582 }
583 return BeqParamDict
586def make_lognormal_MedShkDstn(
587 T_cycle,
588 MedShkAvg,
589 MedShkStd,
590 MedShkCount,
591 MedShkCountTail,
592 RNG,
593 MedShkTailBound=[0.0, 0.9],
594):
595 r"""
596 Constructs discretized lognormal distributions of medical preference shocks
597 for each period in the cycle.
599 .. math::
600 \text{ medShk}_t \sim \exp(\mathcal{N}(\textbf{MedShkStd}^2)) \\
601 \mathbb{E}[\text{medShk}_t]=\textbf{MedShkAvg}
604 Parameters
605 ----------
606 T_cycle : int
607 Number of non-terminal periods in the agent's cycle.
608 MedShkAvg : [float]
609 Mean of medical needs shock in each period of the problem.
610 MedShkStd : [float]
611 Standard deviation of log medical needs shock in each period of the problem.
612 MedShkCount : int
613 Number of equiprobable nodes in the "body" of the discretization.
614 MedShkCountTail : int
615 Number of nodes in each "tail" of the discretization.
616 RNG : RandomState
617 The AgentType's internal random number generator.
618 MedShkTailBound : [float,float]
619 CDF bounds for the tail of the discretization.
621 Returns
622 -------
623 MedShkDstn : [DiscreteDistribuion]
624 """
625 MedShkDstn = [] # empty list for medical shock distribution each period
626 for t in range(T_cycle):
627 # get shock distribution parameters
628 MedShkAvg_t = MedShkAvg[t]
629 MedShkStd_t = MedShkStd[t]
630 MedShkDstn_t = Lognormal(
631 mu=np.log(MedShkAvg_t) - 0.5 * MedShkStd_t**2,
632 sigma=MedShkStd_t,
633 seed=RNG.integers(0, 2**31 - 1),
634 ).discretize(
635 N=MedShkCount,
636 method="equiprobable",
637 tail_N=MedShkCountTail,
638 tail_bound=MedShkTailBound,
639 )
640 MedShkDstn_t = add_discrete_outcome_constant_mean(
641 MedShkDstn_t, 0.0, 0.0, sort=True
642 ) # add point at zero with no probability
643 MedShkDstn.append(MedShkDstn_t)
644 return MedShkDstn
647def make_continuous_MedShockDstn(
648 MedShkLogMean, MedShkLogStd, MedCostLogMean, MedCostLogStd, MedCorr, T_cycle, RNG
649):
650 """
651 Construct a time-varying list of bivariate lognormals for the medical shocks
652 distribution. This representation uses fully continuous distributions, with
653 no discretization in either dimension.
655 Parameters
656 ----------
657 MedShkLogMean : [float]
658 Age-varying list of means of log medical needs (utility) shocks.
659 MedShkLogStd : [float]
660 Age-varying list of standard deviations of log medical needs (utility) shocks.
661 MedCostLogMean : [float]
662 Age-varying list of means of log medical expense shocks.
663 MedCostLogStd : [float]
664 Age-varying list of standard deviations of log medical expense shocks..
665 MedCorr : [float]
666 Age-varying correlation coefficient between log medical expenses and utility shocks.
667 T_cycle : int
668 Number of periods in the agent's sequence.
669 RNG : RandomState
670 Random number generator for this type.
672 Returns
673 -------
674 MedShockDstn : [MultivariateLognormal]
675 Age-varying list of bivariate lognormal distributions, ordered as (MedCost,MedShk).
676 """
677 MedShockDstn = []
678 for t in range(T_cycle):
679 s1 = MedCostLogStd[t]
680 s2 = MedShkLogStd[t]
681 diag = MedCorr[t] * s1 * s2
682 S = np.array([[s1, diag], [diag, s2]])
683 M = np.array([MedCostLogMean[t], MedShkLogMean[t]])
684 seed_t = RNG.integers(0, 2**31 - 1)
685 dstn_t = MultivariateLognormal(mu=M, sigma=S, seed=seed_t)
686 MedShockDstn.append(dstn_t)
687 return MedShockDstn
690def make_MedShock_solution_terminal(
691 CRRA, CRRAmed, MedShkDstn, MedPrice, aXtraGrid, pLvlGrid, CubicBool
692):
693 """
694 Construct the terminal period solution for this type. Similar to other models,
695 optimal behavior involves spending all available market resources; however,
696 the agent must split his resources between consumption and medical care.
698 Parameters
699 ----------
700 None
702 Returns:
703 --------
704 None
705 """
706 # Take last period data, whichever way time is flowing
707 MedPrice = MedPrice[-1]
708 MedShkVals = MedShkDstn[-1].atoms.flatten()
709 MedShkPrbs = MedShkDstn[-1].pmv
711 # Initialize grids of medical need shocks, market resources, and optimal consumption
712 MedShkGrid = MedShkVals
713 xLvlMin = np.min(aXtraGrid) * np.min(pLvlGrid)
714 xLvlMax = np.max(aXtraGrid) * np.max(pLvlGrid)
715 xLvlGrid = make_grid_exp_mult(xLvlMin, xLvlMax, 3 * aXtraGrid.size, 8)
716 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
718 # Make the policy functions for the terminal period
719 xFunc_terminal = TrilinearInterp(
720 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
721 trivial_grid,
722 trivial_grid,
723 trivial_grid,
724 )
725 policyFunc_terminal = MedShockPolicyFunc(
726 xFunc_terminal,
727 xLvlGrid,
728 MedShkGrid,
729 MedPrice,
730 CRRA,
731 CRRAmed,
732 xLvlCubicBool=CubicBool,
733 )
734 cFunc_terminal = cThruXfunc(xFunc_terminal, policyFunc_terminal.cFunc)
735 MedFunc_terminal = MedThruXfunc(xFunc_terminal, policyFunc_terminal.cFunc, MedPrice)
737 # Calculate optimal consumption on a grid of market resources and medical shocks
738 mLvlGrid = xLvlGrid
739 mLvlGrid_tiled = np.tile(
740 np.reshape(mLvlGrid, (mLvlGrid.size, 1)), (1, MedShkGrid.size)
741 )
742 pLvlGrid_tiled = np.ones_like(
743 mLvlGrid_tiled
744 ) # permanent income irrelevant in terminal period
745 MedShkGrid_tiled = np.tile(
746 np.reshape(MedShkVals, (1, MedShkGrid.size)), (mLvlGrid.size, 1)
747 )
748 cLvlGrid, MedGrid = policyFunc_terminal(
749 mLvlGrid_tiled, pLvlGrid_tiled, MedShkGrid_tiled
750 )
752 # Integrate marginal value across shocks to get expected marginal value
753 vPgrid = cLvlGrid ** (-CRRA)
754 vPgrid[np.isinf(vPgrid)] = 0.0 # correct for issue at bottom edges
755 PrbGrid = np.tile(np.reshape(MedShkPrbs, (1, MedShkGrid.size)), (mLvlGrid.size, 1))
756 vP_expected = np.sum(vPgrid * PrbGrid, axis=1)
758 # Construct the marginal (marginal) value function for the terminal period
759 vPnvrs = vP_expected ** (-1.0 / CRRA)
760 vPnvrs[0] = 0.0
761 vPnvrsFunc = BilinearInterp(
762 np.tile(np.reshape(vPnvrs, (vPnvrs.size, 1)), (1, trivial_grid.size)),
763 mLvlGrid,
764 trivial_grid,
765 )
766 vPfunc_terminal = MargValueFuncCRRA(vPnvrsFunc, CRRA)
767 vPPfunc_terminal = MargMargValueFuncCRRA(vPnvrsFunc, CRRA)
769 # Integrate value across shocks to get expected value
770 vGrid = utility(cLvlGrid, rho=CRRA) + MedShkGrid_tiled * utility(
771 MedGrid, rho=CRRAmed
772 )
773 # correct for issue when MedShk=0
774 vGrid[:, 0] = utility(cLvlGrid[:, 0], rho=CRRA)
775 vGrid[np.isinf(vGrid)] = 0.0 # correct for issue at bottom edges
776 v_expected = np.sum(vGrid * PrbGrid, axis=1)
778 # Construct the value function for the terminal period
779 vNvrs = utility_inv(v_expected, rho=CRRA)
780 vNvrs[0] = 0.0
781 vNvrsP = vP_expected * utility_invP(v_expected, rho=CRRA)
782 # TODO: Figure out MPCmax in this model
783 vNvrsP[0] = 0.0
784 tempFunc = CubicInterp(mLvlGrid, vNvrs, vNvrsP)
785 vNvrsFunc = LinearInterpOnInterp1D([tempFunc, tempFunc], trivial_grid)
786 vFunc_terminal = ValueFuncCRRA(vNvrsFunc, CRRA)
788 # Make and return the terminal period solution
789 solution_terminal = ConsumerSolution(
790 cFunc=cFunc_terminal,
791 vFunc=vFunc_terminal,
792 vPfunc=vPfunc_terminal,
793 vPPfunc=vPPfunc_terminal,
794 hNrm=0.0,
795 mNrmMin=0.0,
796 )
797 solution_terminal.MedFunc = MedFunc_terminal
798 solution_terminal.policyFunc = policyFunc_terminal
799 # Track absolute human wealth and minimum market wealth by permanent income
800 solution_terminal.hLvl = ConstantFunction(0.0)
801 solution_terminal.mLvlMin = ConstantFunction(0.0)
802 return solution_terminal
805###############################################################################
808def solve_one_period_ConsMedShock(
809 solution_next,
810 IncShkDstn,
811 MedShkDstn,
812 LivPrb,
813 DiscFac,
814 CRRA,
815 CRRAmed,
816 Rfree,
817 MedPrice,
818 pLvlNextFunc,
819 BoroCnstArt,
820 aXtraGrid,
821 pLvlGrid,
822 vFuncBool,
823 CubicBool,
824):
825 """
826 Class for solving the one period problem for the "medical shocks" model, in
827 which consumers receive shocks to permanent and transitory income as well as
828 shocks to "medical need"-- multiplicative utility shocks for a second good.
830 Parameters
831 ----------
832 solution_next : ConsumerSolution
833 The solution to next period's one period problem.
834 IncShkDstn : distribution.Distribution
835 A discrete approximation to the income process between the period being
836 solved and the one immediately following (in solution_next).
837 MedShkDstn : distribution.Distribution
838 Discrete distribution of the multiplicative utility shifter for medical care.
839 LivPrb : float
840 Survival probability; likelihood of being alive at the beginning of
841 the succeeding period.
842 DiscFac : float
843 Intertemporal discount factor for future utility.
844 CRRA : float
845 Coefficient of relative risk aversion for composite consumption.
846 CRRAmed : float
847 Coefficient of relative risk aversion for medical care.
848 Rfree : float
849 Risk free interest factor on end-of-period assets.
850 MedPrice : float
851 Price of unit of medical care relative to unit of consumption.
852 pLvlNextFunc : float
853 Expected permanent income next period as a function of current pLvl.
854 BoroCnstArt: float or None
855 Borrowing constraint for the minimum allowable assets to end the
856 period with.
857 aXtraGrid: np.array
858 Array of "extra" end-of-period (normalized) asset values-- assets
859 above the absolute minimum acceptable level.
860 pLvlGrid: np.array
861 Array of permanent income levels at which to solve the problem.
862 vFuncBool: boolean
863 An indicator for whether the value function should be computed and
864 included in the reported solution.
865 CubicBool: boolean
866 An indicator for whether the solver should use cubic or linear inter-
867 polation.
869 Returns
870 -------
871 solution_now : ConsumerSolution
872 Solution to this period's consumption-saving problem.
873 """
874 # Define the utility functions for this period
875 uFunc = UtilityFuncCRRA(CRRA)
876 uMed = UtilityFuncCRRA(CRRAmed) # Utility function for medical care
877 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
879 # Unpack next period's income shock distribution
880 ShkPrbsNext = IncShkDstn.pmv
881 PermShkValsNext = IncShkDstn.atoms[0]
882 TranShkValsNext = IncShkDstn.atoms[1]
883 PermShkMinNext = np.min(PermShkValsNext)
884 TranShkMinNext = np.min(TranShkValsNext)
885 MedShkPrbs = MedShkDstn.pmv
886 MedShkVals = MedShkDstn.atoms.flatten()
888 # Calculate the probability that we get the worst possible income draw
889 IncNext = PermShkValsNext * TranShkValsNext
890 WorstIncNext = PermShkMinNext * TranShkMinNext
891 WorstIncPrb = np.sum(ShkPrbsNext[IncNext == WorstIncNext])
892 # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing
894 # Unpack next period's (marginal) value function
895 vFuncNext = solution_next.vFunc # This is None when vFuncBool is False
896 vPfuncNext = solution_next.vPfunc
897 vPPfuncNext = solution_next.vPPfunc # This is None when CubicBool is False
899 # Update the bounding MPCs and PDV of human wealth:
900 PatFac = ((Rfree * DiscFacEff) ** (1.0 / CRRA)) / Rfree
901 try:
902 MPCminNow = 1.0 / (1.0 + PatFac / solution_next.MPCmin)
903 except:
904 MPCminNow = 0.0
905 mLvlMinNext = solution_next.mLvlMin
907 # TODO: Deal with this unused code for the upper bound of MPC (should be a function now)
908 # Ex_IncNext = np.dot(ShkPrbsNext, TranShkValsNext * PermShkValsNext)
909 # hNrmNow = 0.0
910 # temp_fac = (WorstIncPrb ** (1.0 / CRRA)) * PatFac
911 # MPCmaxNow = 1.0 / (1.0 + temp_fac / solution_next.MPCmax)
913 # Define some functions for calculating future expectations
914 def calc_pLvl_next(S, p):
915 return pLvlNextFunc(p) * S["PermShk"]
917 def calc_mLvl_next(S, a, p_next):
918 return Rfree * a + p_next * S["TranShk"]
920 def calc_hLvl(S, p):
921 pLvl_next = calc_pLvl_next(S, p)
922 hLvl = S["TranShk"] * pLvl_next + solution_next.hLvl(pLvl_next)
923 return hLvl
925 def calc_v_next(S, a, p):
926 pLvl_next = calc_pLvl_next(S, p)
927 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
928 v_next = vFuncNext(mLvl_next, pLvl_next)
929 return v_next
931 def calc_vP_next(S, a, p):
932 pLvl_next = calc_pLvl_next(S, p)
933 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
934 vP_next = vPfuncNext(mLvl_next, pLvl_next)
935 return vP_next
937 def calc_vPP_next(S, a, p):
938 pLvl_next = calc_pLvl_next(S, p)
939 mLvl_next = calc_mLvl_next(S, a, pLvl_next)
940 vPP_next = vPPfuncNext(mLvl_next, pLvl_next)
941 return vPP_next
943 # Construct human wealth level as a function of productivity pLvl
944 hLvlGrid = 1.0 / Rfree * expected(calc_hLvl, IncShkDstn, args=(pLvlGrid))
945 hLvlNow = LinearInterp(np.insert(pLvlGrid, 0, 0.0), np.insert(hLvlGrid, 0, 0.0))
947 # Make temporary grids of income shocks and next period income values
948 ShkCount = TranShkValsNext.size
949 pLvlCount = pLvlGrid.size
950 PermShkVals_temp = np.tile(
951 np.reshape(PermShkValsNext, (1, ShkCount)), (pLvlCount, 1)
952 )
953 TranShkVals_temp = np.tile(
954 np.reshape(TranShkValsNext, (1, ShkCount)), (pLvlCount, 1)
955 )
956 pLvlNext_temp = (
957 np.tile(
958 np.reshape(pLvlNextFunc(pLvlGrid), (pLvlCount, 1)),
959 (1, ShkCount),
960 )
961 * PermShkVals_temp
962 )
964 # Find the natural borrowing constraint for each persistent income level
965 aLvlMin_candidates = (
966 mLvlMinNext(pLvlNext_temp) - TranShkVals_temp * pLvlNext_temp
967 ) / Rfree
968 aLvlMinNow = np.max(aLvlMin_candidates, axis=1)
969 BoroCnstNat = LinearInterp(
970 np.insert(pLvlGrid, 0, 0.0), np.insert(aLvlMinNow, 0, 0.0)
971 )
973 # Define the minimum allowable mLvl by pLvl as the greater of the natural and artificial borrowing constraints
974 if BoroCnstArt is not None:
975 BoroCnstArt = LinearInterp(np.array([0.0, 1.0]), np.array([0.0, BoroCnstArt]))
976 mLvlMinNow = UpperEnvelope(BoroCnstArt, BoroCnstNat)
977 else:
978 mLvlMinNow = BoroCnstNat
980 # Make the constrained total spending function: spend all market resources
981 trivial_grid = np.array([0.0, 1.0]) # Trivial grid
982 spendAllFunc = TrilinearInterp(
983 np.array([[[0.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [1.0, 1.0]]]),
984 trivial_grid,
985 trivial_grid,
986 trivial_grid,
987 )
988 xFuncNowCnst = VariableLowerBoundFunc3D(spendAllFunc, mLvlMinNow)
990 # Define grids of pLvl and aLvl on which to compute future expectations
991 pLvlCount = pLvlGrid.size
992 aNrmCount = aXtraGrid.size
993 MedCount = MedShkVals.size
994 pLvlNow = np.tile(pLvlGrid, (aNrmCount, 1)).transpose()
995 aLvlNow = np.tile(aXtraGrid, (pLvlCount, 1)) * pLvlNow + BoroCnstNat(pLvlNow)
996 # shape = (pLvlCount,aNrmCount)
997 if pLvlGrid[0] == 0.0: # aLvl turns out badly if pLvl is 0 at bottom
998 aLvlNow[0, :] = aXtraGrid
1000 # Calculate end-of-period marginal value of assets
1001 EndOfPrd_vP = (
1002 DiscFacEff * Rfree * expected(calc_vP_next, IncShkDstn, args=(aLvlNow, pLvlNow))
1003 )
1005 # If the value function has been requested, construct the end-of-period vFunc
1006 if vFuncBool:
1007 # Compute expected value from end-of-period states
1008 EndOfPrd_v = expected(calc_v_next, IncShkDstn, args=(aLvlNow, pLvlNow))
1009 EndOfPrd_v *= DiscFacEff
1011 # Transformed value through inverse utility function to "decurve" it
1012 EndOfPrd_vNvrs = uFunc.inv(EndOfPrd_v)
1013 EndOfPrd_vNvrsP = EndOfPrd_vP * uFunc.derinv(EndOfPrd_v, order=(0, 1))
1015 # Add points at mLvl=zero
1016 EndOfPrd_vNvrs = np.concatenate(
1017 (np.zeros((pLvlCount, 1)), EndOfPrd_vNvrs), axis=1
1018 )
1019 EndOfPrd_vNvrsP = np.concatenate(
1020 (
1021 np.reshape(EndOfPrd_vNvrsP[:, 0], (pLvlCount, 1)),
1022 EndOfPrd_vNvrsP,
1023 ),
1024 axis=1,
1025 )
1026 # This is a very good approximation, vNvrsPP = 0 at the asset minimum
1028 # Make a temporary aLvl grid for interpolating the end-of-period value function
1029 aLvl_temp = np.concatenate(
1030 (
1031 np.reshape(BoroCnstNat(pLvlGrid), (pLvlGrid.size, 1)),
1032 aLvlNow,
1033 ),
1034 axis=1,
1035 )
1037 # Make an end-of-period value function for each persistent income level in the grid
1038 EndOfPrd_vNvrsFunc_list = []
1039 for p in range(pLvlCount):
1040 EndOfPrd_vNvrsFunc_list.append(
1041 CubicInterp(
1042 aLvl_temp[p, :] - BoroCnstNat(pLvlGrid[p]),
1043 EndOfPrd_vNvrs[p, :],
1044 EndOfPrd_vNvrsP[p, :],
1045 )
1046 )
1047 EndOfPrd_vNvrsFuncBase = LinearInterpOnInterp1D(
1048 EndOfPrd_vNvrsFunc_list, pLvlGrid
1049 )
1051 # Re-adjust the combined end-of-period value function to account for the
1052 # natural borrowing constraint shifter and "re-curve" it
1053 EndOfPrd_vNvrsFunc = VariableLowerBoundFunc2D(
1054 EndOfPrd_vNvrsFuncBase, BoroCnstNat
1055 )
1056 EndOfPrd_vFunc = ValueFuncCRRA(EndOfPrd_vNvrsFunc, CRRA)
1058 # Solve the first order condition to get optimal consumption and medical
1059 # spending, then find the endogenous mLvl gridpoints
1060 # Calculate endogenous gridpoints and controls
1061 cLvlNow = np.tile(
1062 np.reshape(uFunc.derinv(EndOfPrd_vP, order=(1, 0)), (1, pLvlCount, aNrmCount)),
1063 (MedCount, 1, 1),
1064 )
1065 MedBaseNow = np.tile(
1066 np.reshape(
1067 uMed.derinv(MedPrice * EndOfPrd_vP, order=(1, 0)),
1068 (1, pLvlCount, aNrmCount),
1069 ),
1070 (MedCount, 1, 1),
1071 )
1072 MedShkVals_tiled = np.tile( # This includes CRRA adjustment
1073 np.reshape(MedShkVals ** (1.0 / CRRAmed), (MedCount, 1, 1)),
1074 (1, pLvlCount, aNrmCount),
1075 )
1076 MedLvlNow = MedShkVals_tiled * MedBaseNow
1077 aLvlNow_tiled = np.tile(
1078 np.reshape(aLvlNow, (1, pLvlCount, aNrmCount)), (MedCount, 1, 1)
1079 )
1080 xLvlNow = cLvlNow + MedPrice * MedLvlNow
1081 mLvlNow = xLvlNow + aLvlNow_tiled
1083 # Limiting consumption is zero as m approaches the natural borrowing constraint
1084 x_for_interpolation = np.concatenate(
1085 (np.zeros((MedCount, pLvlCount, 1)), xLvlNow), axis=-1
1086 )
1087 temp = np.tile(
1088 BoroCnstNat(np.reshape(pLvlGrid, (1, pLvlCount, 1))),
1089 (MedCount, 1, 1),
1090 )
1091 m_for_interpolation = np.concatenate((temp, mLvlNow), axis=-1)
1093 # Make a 3D array of permanent income for interpolation
1094 p_for_interpolation = np.tile(
1095 np.reshape(pLvlGrid, (1, pLvlCount, 1)), (MedCount, 1, aNrmCount + 1)
1096 )
1098 MedShkVals_tiled = np.tile( # This does *not* have the CRRA adjustment
1099 np.reshape(MedShkVals, (MedCount, 1, 1)), (1, pLvlCount, aNrmCount)
1100 )
1102 # Build the set of cFuncs by pLvl, gathered in a list
1103 xFunc_by_pLvl_and_MedShk = [] # Initialize the empty list of lists of 1D xFuncs
1104 if CubicBool:
1105 # Calculate end-of-period marginal marginal value of assets
1106 vPP_fac = DiscFacEff * Rfree * Rfree
1107 EndOfPrd_vPP = expected(calc_vPP_next, IncShkDstn, args=(aLvlNow, pLvlNow))
1108 EndOfPrd_vPP *= vPP_fac
1109 EndOfPrd_vPP = np.tile(
1110 np.reshape(EndOfPrd_vPP, (1, pLvlCount, aNrmCount)), (MedCount, 1, 1)
1111 )
1113 # Calculate the MPC and MPM at each gridpoint
1114 dcda = EndOfPrd_vPP / uFunc.der(np.array(cLvlNow), order=2)
1115 dMedda = EndOfPrd_vPP / (MedShkVals_tiled * uMed.der(MedLvlNow, order=2))
1116 dMedda[0, :, :] = 0.0 # dMedda goes crazy when MedShk=0
1117 MPC = dcda / (1.0 + dcda + MedPrice * dMedda)
1118 MPM = dMedda / (1.0 + dcda + MedPrice * dMedda)
1120 # Convert to marginal propensity to spend
1121 MPX = MPC + MedPrice * MPM
1122 MPX = np.concatenate(
1123 (np.reshape(MPX[:, :, 0], (MedCount, pLvlCount, 1)), MPX), axis=2
1124 ) # NEED TO CALCULATE MPM AT NATURAL BORROWING CONSTRAINT
1125 MPX[0, :, 0] = 1.0
1127 # Loop over each permanent income level and medical shock and make a cubic xFunc
1128 xFunc_by_pLvl_and_MedShk = [] # Initialize the empty list of lists of 1D xFuncs
1129 for i in range(pLvlCount):
1130 temp_list = []
1131 pLvl_i = p_for_interpolation[0, i, 0]
1132 mLvlMin_i = BoroCnstNat(pLvl_i)
1133 for j in range(MedCount):
1134 m_temp = m_for_interpolation[j, i, :] - mLvlMin_i
1135 x_temp = x_for_interpolation[j, i, :]
1136 MPX_temp = MPX[j, i, :]
1137 temp_list.append(CubicInterp(m_temp, x_temp, MPX_temp))
1138 xFunc_by_pLvl_and_MedShk.append(deepcopy(temp_list))
1140 # Basic version: use linear interpolation within a pLvl and MedShk
1141 else:
1142 # Loop over pLvl and then MedShk within that
1143 for i in range(pLvlCount):
1144 temp_list = []
1145 pLvl_i = p_for_interpolation[0, i, 0]
1146 mLvlMin_i = BoroCnstNat(pLvl_i)
1147 for j in range(MedCount):
1148 m_temp = m_for_interpolation[j, i, :] - mLvlMin_i
1149 x_temp = x_for_interpolation[j, i, :]
1150 temp_list.append(LinearInterp(m_temp, x_temp))
1151 xFunc_by_pLvl_and_MedShk.append(deepcopy(temp_list))
1153 # Combine the nested list of linear xFuncs into a single function
1154 pLvl_temp = p_for_interpolation[0, :, 0]
1155 MedShk_temp = MedShkVals_tiled[:, 0, 0]
1156 xFuncUncBase = BilinearInterpOnInterp1D(
1157 xFunc_by_pLvl_and_MedShk, pLvl_temp, MedShk_temp
1158 )
1159 xFuncNowUnc = VariableLowerBoundFunc3D(xFuncUncBase, BoroCnstNat)
1160 # Re-adjust for lower bound of natural borrowing constraint
1162 # Combine the constrained and unconstrained functions into the true consumption function
1163 xFuncNow = LowerEnvelope3D(xFuncNowUnc, xFuncNowCnst)
1165 # Transform the expenditure function into policy functions for consumption and medical care
1166 aug_factor = 2
1167 xLvlGrid = make_grid_exp_mult(
1168 np.min(x_for_interpolation),
1169 np.max(x_for_interpolation),
1170 aug_factor * aNrmCount,
1171 8,
1172 )
1173 policyFuncNow = MedShockPolicyFunc(
1174 xFuncNow,
1175 xLvlGrid,
1176 MedShkVals,
1177 MedPrice,
1178 CRRA,
1179 CRRAmed,
1180 xLvlCubicBool=CubicBool,
1181 )
1182 cFuncNow = cThruXfunc(xFuncNow, policyFuncNow.cFunc)
1183 MedFuncNow = MedThruXfunc(xFuncNow, policyFuncNow.cFunc, MedPrice)
1185 # Make the marginal value function by integrating over medical shocks
1186 # Make temporary grids to evaluate the consumption function
1187 temp_grid = np.tile(
1188 np.reshape(aXtraGrid, (aNrmCount, 1, 1)), (1, pLvlCount, MedCount)
1189 )
1190 aMinGrid = np.tile(
1191 np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount, 1)),
1192 (aNrmCount, 1, MedCount),
1193 )
1194 pGrid = np.tile(np.reshape(pLvlGrid, (1, pLvlCount, 1)), (aNrmCount, 1, MedCount))
1195 mGrid = temp_grid * pGrid + aMinGrid
1196 if pLvlGrid[0] == 0:
1197 mGrid[:, 0, :] = np.tile(np.reshape(aXtraGrid, (aNrmCount, 1)), (1, MedCount))
1198 MedShkGrid = np.tile(
1199 np.reshape(MedShkVals, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
1200 )
1201 probsGrid = np.tile(
1202 np.reshape(MedShkPrbs, (1, 1, MedCount)), (aNrmCount, pLvlCount, 1)
1203 )
1205 # Get optimal consumption (and medical care) for each state
1206 cGrid, MedGrid = policyFuncNow(mGrid, pGrid, MedShkGrid)
1208 # Calculate expected marginal value by "integrating" across medical shocks
1209 vPgrid = uFunc.der(cGrid)
1210 vPnow = np.sum(vPgrid * probsGrid, axis=2)
1212 # Add vPnvrs=0 at m=mLvlMin to close it off at the bottom (and vNvrs=0)
1213 mGrid_small = np.concatenate(
1214 (np.reshape(mLvlMinNow(pLvlGrid), (1, pLvlCount)), mGrid[:, :, 0])
1215 )
1216 vPnvrsNow = np.concatenate(
1217 (np.zeros((1, pLvlCount)), uFunc.derinv(vPnow, order=(1, 0)))
1218 )
1220 # Calculate expected value by "integrating" across medical shocks
1221 if vFuncBool:
1222 # interpolation error sometimes makes Med < 0 (barely), so fix that
1223 MedGrid = np.maximum(MedGrid, 1e-100)
1224 # interpolation error sometimes makes tiny violations, so fix that
1225 aGrid = np.maximum(mGrid - cGrid - MedPrice * MedGrid, aMinGrid)
1226 vGrid = uFunc(cGrid) + MedShkGrid * uMed(MedGrid) + EndOfPrd_vFunc(aGrid, pGrid)
1227 vNow = np.sum(vGrid * probsGrid, axis=2)
1229 # Switch to pseudo-inverse value and add a point at bottom
1230 vNvrsNow = np.concatenate((np.zeros((1, pLvlCount)), uFunc.inv(vNow)), axis=0)
1231 vNvrsPnow = vPnow * uFunc.derinv(vNow, order=(0, 1))
1232 vNvrsPnow = np.concatenate((np.zeros((1, pLvlCount)), vNvrsPnow), axis=0)
1234 # Construct the pseudo-inverse value and marginal value functions over mLvl,pLvl
1235 vPnvrsFunc_by_pLvl = []
1236 vNvrsFunc_by_pLvl = []
1237 # Make a pseudo inverse marginal value function for each pLvl
1238 for j in range(pLvlCount):
1239 pLvl = pLvlGrid[j]
1240 m_temp = mGrid_small[:, j] - mLvlMinNow(pLvl)
1241 vPnvrs_temp = vPnvrsNow[:, j]
1242 vPnvrsFunc_by_pLvl.append(LinearInterp(m_temp, vPnvrs_temp))
1243 if vFuncBool:
1244 vNvrs_temp = vNvrsNow[:, j]
1245 vNvrsP_temp = vNvrsPnow[:, j]
1246 vNvrsFunc_by_pLvl.append(CubicInterp(m_temp, vNvrs_temp, vNvrsP_temp))
1248 # Combine those functions across pLvls, and adjust for the lower bound of mLvl
1249 vPnvrsFuncBase = LinearInterpOnInterp1D(vPnvrsFunc_by_pLvl, pLvlGrid)
1250 vPnvrsFunc = VariableLowerBoundFunc2D(vPnvrsFuncBase, mLvlMinNow)
1251 if vFuncBool:
1252 vNvrsFuncBase = LinearInterpOnInterp1D(vNvrsFunc_by_pLvl, pLvlGrid)
1253 vNvrsFunc = VariableLowerBoundFunc2D(vNvrsFuncBase, mLvlMinNow)
1255 # "Re-curve" the (marginal) value function
1256 vPfuncNow = MargValueFuncCRRA(vPnvrsFunc, CRRA)
1257 if vFuncBool:
1258 vFuncNow = ValueFuncCRRA(vNvrsFunc, CRRA)
1259 else:
1260 vFuncNow = NullFunc()
1262 # If using cubic spline interpolation, construct the marginal marginal value function
1263 if CubicBool:
1264 vPPfuncNow = MargMargValueFuncCRRA(vPfuncNow.cFunc, CRRA)
1265 else:
1266 vPPfuncNow = NullFunc()
1268 # Package and return the solution object
1269 solution_now = ConsumerSolution(
1270 cFunc=cFuncNow,
1271 vFunc=vFuncNow,
1272 vPfunc=vPfuncNow,
1273 vPPfunc=vPPfuncNow,
1274 mNrmMin=0.0, # Not a normalized model, mLvlMin will be added below
1275 hNrm=0.0, # Not a normalized model, hLvl will be added below
1276 MPCmin=MPCminNow,
1277 MPCmax=0.0, # This should be a function, need to make it
1278 )
1279 solution_now.hLvl = hLvlNow
1280 solution_now.mLvlMin = mLvlMinNow
1281 solution_now.MedFunc = MedFuncNow
1282 solution_now.policyFunc = policyFuncNow
1283 return solution_now
1286###############################################################################
1288# Make a constructor dictionary for the general income process consumer type
1289medshock_constructor_dict = {
1290 "IncShkDstn": construct_lognormal_income_process_unemployment,
1291 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
1292 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
1293 "aXtraGrid": make_assets_grid,
1294 "pLvlPctiles": make_basic_pLvlPctiles,
1295 "pLvlGrid": make_pLvlGrid_by_simulation,
1296 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
1297 "MedShkDstn": make_lognormal_MedShkDstn,
1298 "solution_terminal": make_MedShock_solution_terminal,
1299 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
1300 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
1301}
1303# Make a dictionary with parameters for the default constructor for kNrmInitDstn
1304default_kNrmInitDstn_params = {
1305 "kLogInitMean": -6.0, # Mean of log initial capital
1306 "kLogInitStd": 1.0, # Stdev of log initial capital
1307 "kNrmInitCount": 15, # Number of points in initial capital discretization
1308}
1310# Make a dictionary with parameters for the default constructor for pLvlInitDstn
1311default_pLvlInitDstn_params = {
1312 "pLogInitMean": 0.0, # Mean of log permanent income
1313 "pLogInitStd": 0.4, # Stdev of log permanent income
1314 "pLvlInitCount": 15, # Number of points in initial capital discretization
1315}
1317# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
1318default_IncShkDstn_params = {
1319 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
1320 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
1321 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
1322 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
1323 "UnempPrb": 0.05, # Probability of unemployment while working
1324 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
1325 "T_retire": 0, # Period of retirement (0 --> no retirement)
1326 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
1327 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
1328}
1330# Default parameters to make aXtraGrid using make_assets_grid
1331default_aXtraGrid_params = {
1332 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
1333 "aXtraMax": 30, # Maximum end-of-period "assets above minimum" value
1334 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
1335 "aXtraCount": 32, # Number of points in the grid of "assets above minimum"
1336 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
1337}
1339# Default parameters to make pLvlGrid using make_basic_pLvlPctiles
1340default_pLvlPctiles_params = {
1341 "pLvlPctiles_count": 19, # Number of points in the "body" of the grid
1342 "pLvlPctiles_bound": [0.05, 0.95], # Percentile bounds of the "body"
1343 "pLvlPctiles_tail_count": 4, # Number of points in each tail of the grid
1344 "pLvlPctiles_tail_order": np.e, # Scaling factor for points in each tail
1345}
1347# Default parameters to make pLvlGrid using make_trivial_pLvlNextFunc
1348default_pLvlGrid_params = {
1349 "pLvlInitMean": 0.0, # Mean of log initial permanent income
1350 "pLvlInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
1351 # "pLvlPctiles": pLvlPctiles, # Percentiles of permanent income to use for the grid
1352 "pLvlExtra": [
1353 0.0001
1354 ], # Additional permanent income points to automatically add to the grid, optional
1355}
1357# Default parameters to make pLvlNextFunc using make_AR1_style_pLvlNextFunc
1358default_pLvlNextFunc_params = {
1359 "PermGroFac": [1.0], # Permanent income growth factor
1360 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
1361}
1363# Default parameters to make MedShkDstn using make_lognormal_MedShkDstn
1364default_MedShkDstn_params = {
1365 "MedShkAvg": [0.1], # Average of medical need shocks
1366 "MedShkStd": [4.0], # Standard deviation of (log) medical need shocks
1367 "MedShkCount": 5, # Number of medical shock points in "body"
1368 "MedShkCountTail": 15, # Number of medical shock points in "tail" (upper only)
1369 "MedPrice": [1.5], # Relative price of a unit of medical care
1370}
1372# Make a dictionary to specify a medical shocks consumer type
1373init_medical_shocks = {
1374 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
1375 "cycles": 1, # Finite, non-cyclic model
1376 "T_cycle": 1, # Number of periods in the cycle for this agent type
1377 "pseudo_terminal": False, # Terminal period really does exist
1378 "constructors": medshock_constructor_dict, # See dictionary above
1379 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
1380 "CRRA": 2.0, # Coefficient of relative risk aversion on consumption
1381 "CRRAmed": 3.0, # Coefficient of relative risk aversion on medical care
1382 "Rfree": [1.03], # Interest factor on retained assets
1383 "DiscFac": 0.96, # Intertemporal discount factor
1384 "LivPrb": [0.99], # Survival probability after each period
1385 "BoroCnstArt": 0.0, # Artificial borrowing constraint
1386 "vFuncBool": False, # Whether to calculate the value function during solution
1387 "CubicBool": False, # Whether to use cubic spline interpolation when True
1388 # (Uses linear spline interpolation for cFunc when False)
1389 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
1390 "AgentCount": 10000, # Number of agents of this type
1391 "T_age": None, # Age after which simulated agents are automatically killed
1392 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
1393 # (The portion of PermGroFac attributable to aggregate productivity growth)
1394 "NewbornTransShk": False, # Whether Newborns have transitory shock
1395 # ADDITIONAL OPTIONAL PARAMETERS
1396 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
1397 # (Forces Newborns to follow solution path of the agent they replaced if True)
1398 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
1399}
1400init_medical_shocks.update(default_IncShkDstn_params)
1401init_medical_shocks.update(default_aXtraGrid_params)
1402init_medical_shocks.update(default_pLvlPctiles_params)
1403init_medical_shocks.update(default_pLvlGrid_params)
1404init_medical_shocks.update(default_MedShkDstn_params)
1405init_medical_shocks.update(default_pLvlNextFunc_params)
1406init_medical_shocks.update(default_pLvlInitDstn_params)
1407init_medical_shocks.update(default_kNrmInitDstn_params)
1410class MedShockConsumerType(PersistentShockConsumerType):
1411 r"""
1412 A consumer type based on GenIncShockConsumerType, with two types of consumption goods (medical and nonmedical) and random shocks to medical utility.
1414 .. math::
1415 \begin{eqnarray*}
1416 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})], \\
1417 A_t &=& M_t - X_t, \\
1418 X_t &=& C_t +med_t \textbf{ medPrice}_t,\\
1419 A_t/P_t &\geq& \underline{a}, \\
1420 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
1421 Y_{t+1} &=& P_{t+1} \theta_{t+1}
1422 M_{t+1} &=& R A_t + Y_{t+1}, \\
1423 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
1424 \eta_t &~\sim& G_t,\\
1425 U_t(C, med; \eta) &=& \frac{C^{1-\rho}}{1-\rho} +\eta \frac{med^{1-\nu}}{1-\nu}.
1426 \end{eqnarray*}
1429 Constructors
1430 ------------
1431 IncShkDstn: Constructor, :math:`\psi`, :math:`\theta`
1432 The agent's income shock distributions.
1433 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.construct_lognormal_income_process_unemployment`
1434 aXtraGrid: Constructor
1435 The agent's asset grid.
1436 Its default constructor is :func:`HARK.utilities.make_assets_grid`
1437 pLvlNextFunc: Constructor
1438 An arbitrary function used to evolve the GenIncShockConsumerType's permanent income
1439 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_trivial_pLvlNextFunc`
1440 pLvlGrid: Constructor
1441 The agent's pLvl grid
1442 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_pLvlGrid_by_simulation`
1443 pLvlPctiles: Constructor
1444 The agents income level percentile grid
1445 Its default constructor is :func:`HARK.Calibration.Income.IncomeProcesses.make_basic_pLvlPctiles`
1446 MedShkDstn: Constructor, :math:`\text{medShk}`
1447 The agent's Medical utility shock distribution.
1448 Its default constructor is :func:`HARK.ConsumptionSaving.ConsMedModel.make_lognormal_MedShkDstn`
1450 Solving Parameters
1451 ------------------
1452 cycles: int
1453 0 specifies an infinite horizon model, 1 specifies a finite model.
1454 T_cycle: int
1455 Number of periods in the cycle for this agent type.
1456 CRRA: float, :math:`\rho`
1457 Coefficient of Relative Risk Aversion for consumption.
1458 CRRAmed: float, :math:`\nu`
1459 Coefficient of Relative Risk Aversion for medical care.
1460 Rfree: float or list[float], time varying, :math:`\mathsf{R}`
1461 Risk Free interest rate. Pass a list of floats to make Rfree time varying.
1462 DiscFac: float, :math:`\beta`
1463 Intertemporal discount factor.
1464 LivPrb: list[float], time varying, :math:`1-\mathsf{D}`
1465 Survival probability after each period.
1466 PermGroFac: list[float], time varying, :math:`\Gamma`
1467 Permanent income growth factor.
1468 BoroCnstArt: float, :math:`\underline{a}`
1469 The minimum Asset/Perminant Income ratio, None to ignore.
1470 vFuncBool: bool
1471 Whether to calculate the value function during solution.
1472 CubicBool: bool
1473 Whether to use cubic spline interpoliation.
1475 Simulation Parameters
1476 ---------------------
1477 AgentCount: int
1478 Number of agents of this kind that are created during simulations.
1479 T_age: int
1480 Age after which to automatically kill agents, None to ignore.
1481 T_sim: int, required for simulation
1482 Number of periods to simulate.
1483 track_vars: list[strings]
1484 List of variables that should be tracked when running the simulation.
1485 For this agent, the options are 'Med', 'MedShk', 'PermShk', 'TranShk', 'aLvl', 'cLvl', 'mLvl', 'pLvl', and 'who_dies'.
1487 PermShk is the agent's permanent income shock
1489 MedShk is the agent's medical utility shock
1491 TranShk is the agent's transitory income shock
1493 aLvl is the nominal asset level
1495 cLvl is the nominal consumption level
1497 Med is the nominal medical spending level
1499 mLvl is the nominal market resources
1501 pLvl is the permanent income level
1503 who_dies is the array of which agents died
1504 aNrmInitMean: float
1505 Mean of Log initial Normalized Assets.
1506 aNrmInitStd: float
1507 Std of Log initial Normalized Assets.
1508 pLvlInitMean: float
1509 Mean of Log initial permanent income.
1510 pLvlInitStd: float
1511 Std of Log initial permanent income.
1512 PermGroFacAgg: float
1513 Aggregate permanent income growth factor (The portion of PermGroFac attributable to aggregate productivity growth).
1514 PerfMITShk: boolean
1515 Do Perfect Foresight MIT Shock (Forces Newborns to follow solution path of the agent they replaced if True).
1516 NewbornTransShk: boolean
1517 Whether Newborns have transitory shock.
1519 Attributes
1520 ----------
1521 solution: list[Consumer solution object]
1522 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
1523 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
1525 Unlike other models with this solution type, this model's variables are NOT normalized.
1526 The solution functions additionally depend on the permanent income level and the medical shock.
1527 For example, :math:`C=\text{cFunc}(M,P,MedShk)`.
1528 hNrm has been replaced by hLvl which is a function of permanent income.
1529 MPC max has not yet been implemented for this class. It will be a function of permanent income.
1531 This solution has two additional functions
1532 :math:`\text{Med}=\text{MedFunc}(M,P,\text{MedShk})`: returns the agent's spending on Medical care
1534 :math:`[C,Med]=\text{policyFunc}(M,P,\text{MedShk})`: returns the agent's spending on consumption and Medical care as numpy arrays
1536 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
1537 history: Dict[Array]
1538 Created by running the :func:`.simulate()` method.
1539 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
1540 Visit :class:`HARK.core.AgentType.simulate` for more information.
1541 """
1543 default_ = {
1544 "params": init_medical_shocks,
1545 "solver": solve_one_period_ConsMedShock,
1546 "model": "ConsMedShock.yaml",
1547 }
1549 time_vary_ = PersistentShockConsumerType.time_vary_ + ["MedPrice", "MedShkDstn"]
1550 time_inv_ = PersistentShockConsumerType.time_inv_ + ["CRRAmed"]
1551 shock_vars_ = PersistentShockConsumerType.shock_vars_ + ["MedShk"]
1552 state_vars = PersistentShockConsumerType.state_vars + ["mLvl"]
1553 distributions = [
1554 "IncShkDstn",
1555 "PermShkDstn",
1556 "TranShkDstn",
1557 "kNrmInitDstn",
1558 "pLvlInitDstn",
1559 "MedShkDstn",
1560 ]
1562 def pre_solve(self):
1563 self.construct("solution_terminal")
1565 def get_shocks(self):
1566 """
1567 Gets permanent and transitory income shocks for this period as well as medical need shocks
1568 and the price of medical care.
1570 Parameters
1571 ----------
1572 None
1574 Returns
1575 -------
1576 None
1577 """
1578 # Get permanent and transitory income shocks
1579 PersistentShockConsumerType.get_shocks(self)
1581 # Initialize medical shock array and relative price array
1582 MedShkNow = np.zeros(self.AgentCount)
1583 MedPriceNow = np.zeros(self.AgentCount)
1585 # Get shocks for each period of the cycle
1586 for t in range(self.T_cycle):
1587 these = t == self.t_cycle
1588 N = np.sum(these)
1589 if N > 0:
1590 MedShkNow[these] = self.MedShkDstn[t].draw(N)
1591 MedPriceNow[these] = self.MedPrice[t]
1592 self.shocks["MedShk"] = MedShkNow
1593 self.shocks["MedPrice"] = MedPriceNow
1595 def get_controls(self):
1596 """
1597 Calculates consumption and medical care for each consumer of this type
1598 using the consumption and medical care functions.
1600 Parameters
1601 ----------
1602 None
1604 Returns
1605 -------
1606 None
1607 """
1608 cLvlNow = np.zeros(self.AgentCount) + np.nan
1609 MedNow = np.zeros(self.AgentCount) + np.nan
1610 for t in range(self.T_cycle):
1611 these = t == self.t_cycle
1612 cLvlNow[these], MedNow[these] = self.solution[t].policyFunc(
1613 self.state_now["mLvl"][these],
1614 self.state_now["pLvl"][these],
1615 self.shocks["MedShk"][these],
1616 )
1617 self.controls["cLvl"] = cLvlNow
1618 self.controls["Med"] = MedNow
1619 return None
1621 def get_poststates(self):
1622 """
1623 Calculates end-of-period assets for each consumer of this type.
1625 Parameters
1626 ----------
1627 None
1629 Returns
1630 -------
1631 None
1632 """
1633 self.state_now["aLvl"] = (
1634 self.state_now["mLvl"]
1635 - self.controls["cLvl"]
1636 - self.shocks["MedPrice"] * self.controls["Med"]
1637 )
1639 # moves now to prev
1640 AgentType.get_poststates(self)
1643###############################################################################
1646class ConsMedExtMargSolution(MetricObject):
1647 """
1648 Representation of the solution to one period's problem in the extensive margin
1649 medical expense model. If no inputs are passed, a trivial object is constructed,
1650 which can be used as the pseudo-terminal solution.
1652 Parameters
1653 ----------
1654 vFunc_by_pLvl : [function]
1655 List of beginning-of-period value functions over kLvl, by pLvl.
1656 vPfunc_by_pLvl : [function]
1657 List of beginning-of-period marginal functions over kLvl, by pLvl.
1658 cFunc_by_pLvl : [function]
1659 List of consumption functions over bLvl, by pLvl.
1660 vNvrsFuncMid_by_pLvl : [function]
1661 List of pseudo-inverse value function for consumption phase over bLvl, by pLvl.
1662 ExpMedFunc : function
1663 Expected medical care as a function of mLvl and pLvl, just before medical
1664 shock is realized.
1665 CareProbFunc : function
1666 Probability of getting medical treatment as a function of mLvl and pLvl,
1667 just before medical shock is realized.
1668 pLvl : np.array
1669 Grid of permanent income levels during the period (after shocks).
1670 CRRA : float
1671 Coefficient of relative risk aversion
1672 """
1674 distance_criteria = ["cFunc"]
1676 def __init__(
1677 self,
1678 vFunc_by_pLvl=None,
1679 vPfunc_by_pLvl=None,
1680 cFunc_by_pLvl=None,
1681 vNvrsFuncMid_by_pLvl=None,
1682 ExpMedFunc=None,
1683 CareProbFunc=None,
1684 pLvl=None,
1685 CRRA=None,
1686 ):
1687 self.pLvl = pLvl
1688 self.CRRA = CRRA
1689 if vFunc_by_pLvl is None:
1690 self.vFunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1691 else:
1692 self.vFunc_by_pLvl = vFunc_by_pLvl
1693 if vPfunc_by_pLvl is None:
1694 self.vPfunc_by_pLvl = pLvl.size * [ConstantFunction(0.0)]
1695 else:
1696 self.vPfunc_by_pLvl = vPfunc_by_pLvl
1697 if cFunc_by_pLvl is not None:
1698 self.cFunc = LinearInterpOnInterp1D(cFunc_by_pLvl, pLvl)
1699 else:
1700 self.cFunc = None
1701 if vNvrsFuncMid_by_pLvl is not None:
1702 vNvrsFuncMid = LinearInterpOnInterp1D(vNvrsFuncMid_by_pLvl, pLvl)
1703 self.vFuncMid = ValueFuncCRRA(vNvrsFuncMid, CRRA, illegal_value=-np.inf)
1704 if ExpMedFunc is not None:
1705 self.ExpMedFunc = ExpMedFunc
1706 if CareProbFunc is not None:
1707 self.CareProbFunc = CareProbFunc
1710def make_MedExtMarg_solution_terminal(pLogCount):
1711 """
1712 Construct a trivial pseudo-terminal solution for the extensive margin medical
1713 spending model: a list of constant zero functions for (marginal) value. The
1714 only piece of information needed for this is how many such functions to include.
1715 """
1716 pLvl_terminal = np.arange(pLogCount)
1717 solution_terminal = ConsMedExtMargSolution(pLvl=pLvl_terminal)
1718 return solution_terminal
1721###############################################################################
1724def solve_one_period_ConsMedExtMarg(
1725 solution_next,
1726 DiscFac,
1727 CRRA,
1728 BeqFac,
1729 BeqShift,
1730 Rfree,
1731 LivPrb,
1732 MedShkLogMean,
1733 MedShkLogStd,
1734 MedCostLogMean,
1735 MedCostLogStd,
1736 MedCorr,
1737 MedCostBot,
1738 MedCostTop,
1739 MedCostCount,
1740 aNrmGrid,
1741 pLogGrid,
1742 pLvlMean,
1743 TranShkDstn,
1744 pLogMrkvArray,
1745 mNrmGrid,
1746 kLvlGrid,
1747):
1748 """
1749 Solve one period of the "extensive margin medical care" model. Each period, the
1750 agent receives a persistent and transitory shock to income, and then a medical
1751 shock with two components: utility and cost. He makes a binary choice between
1752 paying the cost in medical expenses or suffering the utility loss, then makes
1753 his ordinary consumption-saving decision (technically made simultaneously, but
1754 solved as if sequential). This version has one health state and no insurance choice
1755 and hardcodes a liquidity constraint.
1757 Parameters
1758 ----------
1759 solution_next : ConsMedExtMargSolution
1760 Solution to the succeeding period's problem.
1761 DiscFac : float
1762 Intertemporal discount factor.
1763 CRRA : float
1764 Coefficient of relative risk aversion.
1765 BeqFac : float
1766 Scaling factor for bequest motive.
1767 BeqShift : float
1768 Shifter for bequest motive.
1769 Rfree : float
1770 Risk free return factor on saving.
1771 LivPrb : float
1772 Survival probability from this period to the next one.
1773 MedShkLogMean : float
1774 Mean of log utility shocks, assumed to be lognormally distributed.
1775 MedShkLogStd : float
1776 Stdev of log utility shocks, assumed to be lognormally distributed.
1777 MedCostLogMean : float
1778 Mean of log medical expense shocks, assumed to be lognormally distributed.
1779 MedCostLogStd : float
1780 Stdev of log medical expense shocks, assumed to be lognormally distributed.
1781 MedCorr : float
1782 Correlation coefficient betwen log utility shocks and log medical expense
1783 shocks, assumed to be joint normal (in logs).
1784 MedCostBot : float
1785 Lower bound of medical costs to consider, as standard deviations of log
1786 expenses away from the mean.
1787 MedCostTop : float
1788 Upper bound of medical costs to consider, as standard deviations of log
1789 expenses away from the mean.
1790 MedCostCount : int
1791 Number of points to use when discretizing MedCost
1792 aNrmGrid : np.array
1793 Exogenous grid of end-of-period assets, normalized by income level.
1794 pLogGrid : np.array
1795 Exogenous grid of *deviations from mean* log income level.
1796 pLvlMean : float
1797 Mean income level at this age, for generating actual income levels from
1798 pLogGrid as pLvl = pLvlMean * np.exp(pLogGrid).
1799 TranShkDstn : DiscreteDistribution
1800 Discretized transitory income shock distribution.
1801 pLogMrkvArray : np.array
1802 Markov transition array from beginning-of-period (prior) income levels
1803 to this period's levels. Pre-computed by (e.g.) Tauchen's method.
1804 mNrmGrid : np.array
1805 Exogenous grid of decision-time normalized market resources,
1806 kLvlGrid : np.array
1807 Beginning-of-period capital grid (spanning zero to very high wealth).
1809 Returns
1810 -------
1811 solution_now : ConsMedExtMargSolution
1812 Representation of the solution to this period's problem, including the
1813 beginning-of-period (marginal) value function by pLvl, the consumption
1814 function by pLvl, and the (pseudo-inverse) value function for the consumption
1815 phase (as a list by pLvl).
1816 """
1817 # Define (marginal) utility and bequest motive functions
1818 u = lambda x: CRRAutility(x, rho=CRRA)
1819 uP = lambda x: CRRAutilityP(x, rho=CRRA)
1820 W = lambda x: BeqFac * u(x + BeqShift)
1821 Wp = lambda x: BeqFac * uP(x + BeqShift)
1822 n = lambda x: CRRAutility_inv(x, rho=CRRA)
1824 # Make grids of pLvl x aLvl
1825 pLvl = np.exp(pLogGrid) * pLvlMean
1826 aLvl = np.dot(
1827 np.reshape(aNrmGrid, (aNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvl.size))
1828 )
1829 aLvl = np.concatenate([np.zeros((1, pLvl.size)), aLvl]) # add zero entries
1831 # Evaluate end-of-period marginal value at each combination of pLvl x aLvl
1832 pLvlCount = pLvl.size
1833 EndOfPrd_vP = np.empty_like(aLvl)
1834 EndOfPrd_v = np.empty_like(aLvl)
1835 for j in range(pLvlCount):
1836 EndOfPrdvFunc_this_pLvl = solution_next.vFunc_by_pLvl[j]
1837 EndOfPrdvPfunc_this_pLvl = solution_next.vPfunc_by_pLvl[j]
1838 EndOfPrd_v[:, j] = DiscFac * LivPrb * EndOfPrdvFunc_this_pLvl(aLvl[:, j])
1839 EndOfPrd_vP[:, j] = DiscFac * LivPrb * EndOfPrdvPfunc_this_pLvl(aLvl[:, j])
1840 EndOfPrd_v += (1.0 - LivPrb) * W(aLvl)
1841 EndOfPrd_vP += (1.0 - LivPrb) * Wp(aLvl)
1843 # Calculate optimal consumption for each (aLvl,pLvl) gridpoint, roll back to bLvl
1844 cLvl = CRRAutilityP_inv(EndOfPrd_vP, CRRA)
1845 bLvl = aLvl + cLvl
1847 # Construct consumption functions over bLvl for each pLvl
1848 cFunc_by_pLvl = []
1849 for j in range(pLvlCount):
1850 cFunc_j = LinearInterp(
1851 np.insert(bLvl[:, j], 0, 0.0), np.insert(cLvl[:, j], 0, 0.0)
1852 )
1853 cFunc_by_pLvl.append(cFunc_j)
1855 # Construct pseudo-inverse value functions over bLvl for each pLvl
1856 v_mid = u(cLvl) + EndOfPrd_v # value of reaching consumption phase
1857 vNvrsFuncMid_by_pLvl = []
1858 for j in range(pLvlCount):
1859 b_cnst = np.linspace(0.001, 0.95, 10) * bLvl[0, j] # constrained wealth levels
1860 c_cnst = b_cnst
1861 v_cnst = u(c_cnst) + EndOfPrd_v[0, j]
1862 b_temp = np.concatenate([b_cnst, bLvl[:, j]])
1863 v_temp = np.concatenate([v_cnst, v_mid[:, j]])
1864 vNvrs_temp = n(v_temp)
1865 vNvrsFunc_j = LinearInterp(
1866 np.insert(b_temp, 0, 0.0), np.insert(vNvrs_temp, 0, 0.0)
1867 )
1868 vNvrsFuncMid_by_pLvl.append(vNvrsFunc_j)
1870 # Make a grid of (log) medical expenses (and probs), cross it with (mLvl,pLvl)
1871 bot = MedCostLogMean + MedCostBot * MedCostLogStd
1872 top = MedCostLogMean + MedCostTop * MedCostLogStd
1873 MedCostLogGrid = np.linspace(bot, top, MedCostCount)
1874 MedCostGrid = np.exp(MedCostLogGrid)
1875 mLvl_base = np.dot(
1876 np.reshape(mNrmGrid, (mNrmGrid.size, 1)), np.reshape(pLvl, (1, pLvlCount))
1877 )
1878 mLvl = np.reshape(mLvl_base, (mNrmGrid.size, pLvlCount, 1))
1879 bLvl_if_care = mLvl - np.reshape(MedCostGrid, (1, 1, MedCostCount))
1880 bLvl_if_not = mLvl_base
1882 # Calculate mean (log) utility shock for each MedCost gridpoint, and conditional stdev
1883 MedShkLog_cond_mean = MedShkLogMean + MedCorr * MedShkLogStd * MedCostLogGrid
1884 MedShkLog_cond_mean = np.reshape(MedShkLog_cond_mean, (1, MedCostCount))
1885 MedShkLog_cond_std = np.sqrt(MedShkLogStd**2 * (1.0 - MedCorr**2))
1886 MedShk_cond_mean = np.exp(MedShkLog_cond_mean + 0.5 * MedShkLog_cond_std**2)
1888 # Initialize (marginal) value function arrays over (mLvl,pLvl,MedCost)
1889 v_at_Dcsn = np.empty_like(bLvl_if_care)
1890 vP_at_Dcsn = np.empty_like(bLvl_if_care)
1891 care_prob_array = np.empty_like(bLvl_if_care)
1892 for j in range(pLvlCount):
1893 # Evaluate value function for (bLvl,pLvl_j), including MedCost=0
1894 v_if_care = u(vNvrsFuncMid_by_pLvl[j](bLvl_if_care[:, j, :]))
1895 v_if_not = np.reshape(
1896 u(vNvrsFuncMid_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1897 )
1898 cant_pay = bLvl_if_care[:, j, :] <= 0.0
1899 v_if_care[cant_pay] = -np.inf
1901 # Find value difference at each gridpoint, convert to MedShk stdev; find prob of care
1902 v_diff = v_if_not - v_if_care
1903 log_v_diff = np.log(v_diff)
1904 crit_stdev = (log_v_diff - MedShkLog_cond_mean) / MedShkLog_cond_std
1905 prob_no_care = norm.cdf(crit_stdev)
1906 prob_get_care = 1.0 - prob_no_care
1907 care_prob_array[:, j, :] = prob_get_care
1909 # Calculate expected MedShk conditional on not getting medical care
1910 crit_z = crit_stdev - MedShkLog_cond_std
1911 MedShk_no_care_cond_mean = 0.5 * MedShk_cond_mean * erfc(crit_z) / prob_no_care
1913 # Compute expected (marginal) value over MedShk for each (mLvl,pLvl_j,MedCost)
1914 v_if_care[cant_pay] = 0.0
1915 v_at_Dcsn[:, j, :] = (
1916 prob_no_care * (v_if_not - MedShk_no_care_cond_mean)
1917 + prob_get_care * v_if_care
1918 )
1919 vP_if_care = uP(cFunc_by_pLvl[j](bLvl_if_care[:, j, :]))
1920 vP_if_not = np.reshape(
1921 uP(cFunc_by_pLvl[j](bLvl_if_not[:, j])), (mNrmGrid.size, 1)
1922 )
1923 vP_if_care[cant_pay] = 0.0
1924 MedShk_rate_of_change = (
1925 norm.pdf(crit_stdev) * (vP_if_care - vP_if_not) * MedShk_no_care_cond_mean
1926 )
1927 vP_at_Dcsn[:, j, :] = (
1928 prob_no_care * vP_if_not
1929 + prob_get_care * vP_if_care
1930 + MedShk_rate_of_change
1931 )
1933 # Compute expected (marginal) value over MedCost for each (mLvl,pLvl)
1934 temp_grid = np.linspace(MedCostBot, MedCostTop, MedCostCount)
1935 MedCost_pmv = norm.pdf(temp_grid)
1936 MedCost_pmv /= np.sum(MedCost_pmv)
1937 MedCost_probs = np.reshape(MedCost_pmv, (1, 1, MedCostCount))
1938 v_before_shk = np.sum(v_at_Dcsn * MedCost_probs, axis=2)
1939 vP_before_shk = np.sum(vP_at_Dcsn * MedCost_probs, axis=2)
1940 vNvrs_before_shk = n(v_before_shk)
1941 vPnvrs_before_shk = CRRAutilityP_inv(vP_before_shk, CRRA)
1943 # Compute expected medical expenses at each state space point
1944 ExpCare_all = care_prob_array * np.reshape(MedCostGrid, (1, 1, MedCostCount))
1945 ExpCare = np.sum(ExpCare_all * MedCost_probs, axis=2)
1946 ProbCare = np.sum(care_prob_array * MedCost_probs, axis=2)
1947 ExpCareFunc_by_pLvl = []
1948 CareProbFunc_by_pLvl = []
1949 for j in range(pLvlCount):
1950 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1951 EC_temp = np.insert(ExpCare[:, j], 0, 0.0)
1952 prob_temp = np.insert(ProbCare[:, j], 0, 0.0)
1953 ExpCareFunc_by_pLvl.append(LinearInterp(m_temp, EC_temp))
1954 CareProbFunc_by_pLvl.append(LinearInterp(m_temp, prob_temp))
1955 ExpCareFunc = LinearInterpOnInterp1D(ExpCareFunc_by_pLvl, pLvl)
1956 ProbCareFunc = LinearInterpOnInterp1D(CareProbFunc_by_pLvl, pLvl)
1958 # Fixing kLvlGrid, compute expected (marginal) value over TranShk for each (kLvl,pLvl)
1959 v_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1960 vP_by_kLvl_and_pLvl = np.empty((kLvlGrid.size, pLvlCount))
1961 for j in range(pLvlCount):
1962 p = pLvl[j]
1964 # Make (marginal) value functions over mLvl for this pLvl
1965 m_temp = np.insert(mLvl_base[:, j], 0, 0.0)
1966 vNvrs_temp = np.insert(vNvrs_before_shk[:, j], 0, 0.0)
1967 vPnvrs_temp = np.insert(vPnvrs_before_shk[:, j], 0, 0.0)
1968 vNvrsFunc_temp = LinearInterp(m_temp, vNvrs_temp)
1969 vPnvrsFunc_temp = LinearInterp(m_temp, vPnvrs_temp)
1970 vFunc_temp = lambda x: u(vNvrsFunc_temp(x))
1971 vPfunc_temp = lambda x: uP(vPnvrsFunc_temp(x))
1973 # Compute expectation over TranShkDstn
1974 v = lambda TranShk, kLvl: vFunc_temp(kLvl + TranShk * p)
1975 vP = lambda TranShk, kLvl: vPfunc_temp(kLvl + TranShk * p)
1976 v_by_kLvl_and_pLvl[:, j] = expected(v, TranShkDstn, args=(kLvlGrid,))
1977 vP_by_kLvl_and_pLvl[:, j] = expected(vP, TranShkDstn, args=(kLvlGrid,))
1979 # Compute expectation over persistent shocks by using pLvlMrkvArray
1980 v_arvl = np.dot(v_by_kLvl_and_pLvl, pLogMrkvArray.T)
1981 vP_arvl = np.dot(vP_by_kLvl_and_pLvl, pLogMrkvArray.T)
1982 vNvrs_arvl = n(v_arvl)
1983 vPnvrs_arvl = CRRAutilityP_inv(vP_arvl, CRRA)
1985 # Construct "arrival" (marginal) value function by pLvl
1986 vFuncArvl_by_pLvl = []
1987 vPfuncArvl_by_pLvl = []
1988 for j in range(pLvlCount):
1989 vNvrsFunc_temp = LinearInterp(kLvlGrid, vNvrs_arvl[:, j])
1990 vPnvrsFunc_temp = LinearInterp(kLvlGrid, vPnvrs_arvl[:, j])
1991 vFuncArvl_by_pLvl.append(ValueFuncCRRA(vNvrsFunc_temp, CRRA))
1992 vPfuncArvl_by_pLvl.append(MargValueFuncCRRA(vPnvrsFunc_temp, CRRA))
1994 # Gather elements and return the solution object
1995 solution_now = ConsMedExtMargSolution(
1996 vFunc_by_pLvl=vFuncArvl_by_pLvl,
1997 vPfunc_by_pLvl=vPfuncArvl_by_pLvl,
1998 cFunc_by_pLvl=cFunc_by_pLvl,
1999 vNvrsFuncMid_by_pLvl=vNvrsFuncMid_by_pLvl,
2000 pLvl=pLvl,
2001 CRRA=CRRA,
2002 ExpMedFunc=ExpCareFunc,
2003 CareProbFunc=ProbCareFunc,
2004 )
2005 return solution_now
2008###############################################################################
2011# Define a dictionary of constructors for the extensive margin medical spending model
2012med_ext_marg_constructors = {
2013 "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
2014 "IncomeProcessDict": make_persistent_income_process_dict,
2015 "pLogGrid": get_it_from("IncomeProcessDict"),
2016 "pLvlMean": get_it_from("IncomeProcessDict"),
2017 "pLogMrkvArray": get_it_from("IncomeProcessDict"),
2018 "IncShkDstn": construct_lognormal_income_process_unemployment,
2019 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
2020 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
2021 "BeqParamDict": reformat_bequest_motive,
2022 "BeqFac": get_it_from("BeqParamDict"),
2023 "BeqShift": get_it_from("BeqParamDict"),
2024 "aNrmGrid": make_assets_grid,
2025 "mNrmGrid": make_market_resources_grid,
2026 "kLvlGrid": make_capital_grid,
2027 "solution_terminal": make_MedExtMarg_solution_terminal,
2028 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
2029 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
2030 "MedShockDstn": make_continuous_MedShockDstn,
2031}
2033# Make a dictionary with default bequest motive parameters
2034default_BeqParam_dict = {
2035 "BeqMPC": 0.1, # Hypothetical "MPC at death"
2036 "BeqInt": 1.0, # Intercept term for hypothetical "consumption function at death"
2037}
2039# Make a dictionary with parameters for the default constructor for kNrmInitDstn
2040default_kNrmInitDstn_params_ExtMarg = {
2041 "kLogInitMean": -6.0, # Mean of log initial capital
2042 "kLogInitStd": 1.0, # Stdev of log initial capital
2043 "kNrmInitCount": 15, # Number of points in initial capital discretization
2044}
2046# Default parameters to make IncomeProcessDict using make_persistent_income_process_dict;
2047# some of these are used by construct_lognormal_income_process_unemployment as well
2048default_IncomeProcess_params = {
2049 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
2050 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
2051 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
2052 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
2053 "UnempPrb": 0.05, # Probability of unemployment while working
2054 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
2055 "T_retire": 0, # Period of retirement (0 --> no retirement)
2056 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
2057 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
2058 "pLogInitMean": 0.0, # Mean of log initial permanent income
2059 "pLogInitStd": 0.4, # Standard deviation of log initial permanent income *MUST BE POSITIVE*
2060 "pLvlInitCount": 25, # Number of discrete nodes in initial permanent income level dstn
2061 "PermGroFac": [1.0], # Permanent income growth factor
2062 "PrstIncCorr": 0.98, # Correlation coefficient on (log) persistent income
2063 "pLogCount": 45, # Number of points in persistent income grid each period
2064 "pLogRange": 3.5, # Upper/lower bound of persistent income, in unconditional standard deviations
2065}
2067# Default parameters to make aNrmGrid using make_assets_grid
2068default_aNrmGrid_params = {
2069 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
2070 "aXtraMax": 40.0, # Maximum end-of-period "assets above minimum" value
2071 "aXtraNestFac": 2, # Exponential nesting factor for aXtraGrid
2072 "aXtraCount": 96, # Number of points in the grid of "assets above minimum"
2073 "aXtraExtra": [0.005, 0.01], # Additional other values to add in grid (optional)
2074}
2076# Default parameters to make mLvlGrid using make_market_resources_grid
2077default_mNrmGrid_params = {
2078 "mNrmMin": 0.001,
2079 "mNrmMax": 40.0,
2080 "mNrmNestFac": 2,
2081 "mNrmCount": 72,
2082 "mNrmExtra": None,
2083}
2085# Default parameters to make kLvlGrid using make_capital_grid
2086default_kLvlGrid_params = {
2087 "kLvlMin": 0.0,
2088 "kLvlMax": 200,
2089 "kLvlCount": 250,
2090 "kLvlOrder": 2.0,
2091}
2093# Default "basic" parameters
2094med_ext_marg_basic_params = {
2095 "constructors": med_ext_marg_constructors,
2096 "cycles": 1, # Lifecycle by default
2097 "T_cycle": 1, # Number of periods in the default sequence
2098 "T_age": None,
2099 "AgentCount": 10000, # Number of agents to simulate
2100 "DiscFac": 0.96, # intertemporal discount factor
2101 "CRRA": 1.5, # coefficient of relative risk aversion
2102 "Rfree": [1.02], # risk free interest factor
2103 "LivPrb": [0.99], # survival probability
2104 "MedCostBot": -3.1, # lower bound of medical cost distribution, in stdevs
2105 "MedCostTop": 5.2, # upper bound of medical cost distribution, in stdevs
2106 "MedCostCount": 84, # number of nodes in medical cost discretization
2107 "MedShkLogMean": [-2.0], # mean of log utility shocks
2108 "MedShkLogStd": [1.5], # standard deviation of log utility shocks
2109 "MedCostLogMean": [-1.0], # mean of log medical expenses
2110 "MedCostLogStd": [1.0], # standard deviation of log medical expenses
2111 "MedCorr": [0.3], # correlation coefficient between utility shock and expenses
2112 "PermGroFacAgg": 1.0, # Aggregate productivity growth rate (leave at 1)
2113 "NewbornTransShk": True, # Whether "newborns" have a transitory income shock
2114}
2116# Combine the dictionaries into a single default dictionary
2117init_med_ext_marg = med_ext_marg_basic_params.copy()
2118init_med_ext_marg.update(default_IncomeProcess_params)
2119init_med_ext_marg.update(default_aNrmGrid_params)
2120init_med_ext_marg.update(default_mNrmGrid_params)
2121init_med_ext_marg.update(default_kLvlGrid_params)
2122init_med_ext_marg.update(default_kNrmInitDstn_params_ExtMarg)
2123init_med_ext_marg.update(default_BeqParam_dict)
2126class MedExtMargConsumerType(PersistentShockConsumerType):
2127 r"""
2128 Class for representing agents in the extensive margin medical expense model.
2129 Such agents have labor income dynamics identical to the "general income process"
2130 model (permanent income is not normalized out), and also experience a medical
2131 shock with two components: medical cost and utility loss. They face a binary
2132 choice of whether to pay the cost or suffer the loss, then make a consumption-
2133 saving decision as normal. To simplify the computation, the joint distribution
2134 of medical shocks is specified as bivariate lognormal. This can be loosened to
2135 accommodate insurance contracts as mappings from total to out-of-pocket expenses.
2136 Can also be extended to include a health process.
2138 .. math::
2139 \begin{eqnarray*}
2140 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}], \\
2141 A_t &=& M_t - C_t - D_t med_t, \\
2142 A_t/ &\geq& 0, \\
2143 D_t &\in& \{0,1\}, \\
2144 P_{t+1} &=& \Gamma_{t+1}(P_t)\psi_{t+1}, \\
2145 Y_{t+1} &=& P_{t+1} \theta_{t+1}
2146 M_{t+1} &=& R A_t + Y_{t+1}, \\
2147 (\psi_{t+1},\theta_{t+1}) &\sim& F_{t+1},\\
2148 (med_t,\eta_t) &~\sim& G_t,\\
2149 U_t(C) &=& \frac{C^{1-\rho}}{1-\rho}.
2150 \end{eqnarray*}
2151 """
2153 default_ = {
2154 "params": init_med_ext_marg,
2155 "solver": solve_one_period_ConsMedExtMarg,
2156 "model": "ConsExtMargMed.yaml",
2157 }
2159 time_vary_ = [
2160 "Rfree",
2161 "LivPrb",
2162 "MedShkLogMean",
2163 "MedShkLogStd",
2164 "MedCostLogMean",
2165 "MedCostLogStd",
2166 "MedCorr",
2167 "pLogGrid",
2168 "pLvlMean",
2169 "TranShkDstn",
2170 "pLogMrkvArray",
2171 "pLvlNextFunc",
2172 "IncShkDstn",
2173 "MedShockDstn",
2174 ]
2175 time_inv_ = [
2176 "DiscFac",
2177 "CRRA",
2178 "BeqFac",
2179 "BeqShift",
2180 "MedCostBot",
2181 "MedCostTop",
2182 "MedCostCount",
2183 "aNrmGrid",
2184 "mNrmGrid",
2185 "kLvlGrid",
2186 ]
2187 shock_vars = ["PermShk", "TranShk", "MedShk", "MedCost"]
2189 def get_shocks(self):
2190 """
2191 Gets permanent and transitory income shocks for this period as well as
2192 medical need and cost shocks.
2193 """
2194 # Get permanent and transitory income shocks
2195 PersistentShockConsumerType.get_shocks(self)
2197 # Initialize medical shock array and cost of care array
2198 MedShkNow = np.zeros(self.AgentCount)
2199 MedCostNow = np.zeros(self.AgentCount)
2201 # Get shocks for each period of the cycle
2202 for t in range(self.T_cycle):
2203 these = t == self.t_cycle
2204 if np.any(these):
2205 N = np.sum(these)
2206 dstn_t = self.MedShockDstn[t]
2207 draws_t = dstn_t.draw(N)
2208 MedCostNow[these] = draws_t[0, :]
2209 MedShkNow[these] = draws_t[1, :]
2210 self.shocks["MedShk"] = MedShkNow
2211 self.shocks["MedCost"] = MedCostNow
2213 def get_controls(self):
2214 """
2215 Finds consumption for each agent, along with whether or not they get care.
2216 """
2217 # Initialize output
2218 cLvlNow = np.empty(self.AgentCount)
2219 CareNow = np.zeros(self.AgentCount, dtype=bool)
2221 # Get states and shocks
2222 mLvl = self.state_now["mLvl"]
2223 pLvl = self.state_now["pLvl"]
2224 MedCost = self.shocks["MedCost"]
2225 MedShk = self.shocks["MedShk"]
2227 # Find remaining resources with and without care
2228 bLvl_no_care = mLvl
2229 bLvl_with_care = mLvl - MedCost
2231 # Get controls for each period of the cycle
2232 for t in range(self.T_cycle):
2233 these = t == self.t_cycle
2234 if np.any(these):
2235 vFunc_t = self.solution[t].vFuncMid
2236 cFunc_t = self.solution[t].cFunc
2238 v_no_care = vFunc_t(bLvl_no_care[these], pLvl[these]) - MedShk[these]
2239 v_if_care = vFunc_t(bLvl_with_care[these], pLvl[these])
2240 get_care = v_if_care > v_no_care
2242 b_temp = bLvl_no_care[these]
2243 b_temp[get_care] = bLvl_with_care[get_care]
2244 cLvlNow[these] = cFunc_t(b_temp, pLvl[these])
2245 CareNow[these] = get_care
2247 # Store the results
2248 self.controls["cLvl"] = cLvlNow
2249 self.controls["Care"] = CareNow
2251 def get_poststates(self):
2252 """
2253 Calculates end-of-period assets for each consumer of this type.
2254 """
2255 self.state_now["MedLvl"] = self.shocks["MedCost"] * self.controls["Care"]
2256 self.state_now["aLvl"] = (
2257 self.state_now["mLvl"] - self.controls["cLvl"] - self.state_now["MedLvl"]
2258 )
2259 # Move now to prev
2260 AgentType.get_poststates(self)