Coverage for HARK/ConsumptionSaving/ConsHealthModel.py: 97%
176 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +0000
1"""
2Classes to represent consumers who make decisions about health investment. The
3first model here is adapted from White (2015).
4"""
6import numpy as np
7from HARK.core import AgentType
8from HARK.distributions import (
9 expected,
10 combine_indep_dstns,
11 Uniform,
12 DiscreteDistribution,
13 DiscreteDistributionLabeled,
14)
15from HARK.Calibration.Income.IncomeProcesses import construct_lognormal_wage_dstn
16from HARK.rewards import CRRAutility, CRRAutility_inv
17from HARK.interpolation import Curvilinear2DInterp
18from HARK.utilities import make_assets_grid
19from HARK.ConsumptionSaving.ConsIndShockModel import make_lognormal_kNrm_init_dstn
21###############################################################################
24# Define a function that yields health produced from investment
25def eval_health_prod(n, alpha, gamma):
26 return (gamma / alpha) * n**alpha
29# Define a function that yields health produced from investment
30def eval_marg_health_prod(n, alpha, gamma):
31 return gamma * n ** (alpha - 1.0)
34# Define a function for computing expectations over next period's (marginal) value
35# from the perspective of end-of-period states, conditional on survival
36def calc_exp_next(shock, a, H, R, rho, alpha, gamma, funcs):
37 m_next = R * a + shock["WageRte"] * H
38 h_next = (1.0 - shock["DeprRte"]) * H
39 vNvrs_next, c_next, n_next = funcs(m_next, h_next)
40 dvdm_next = c_next**-rho
41 dvdh_next = dvdm_next / (gamma * n_next ** (alpha - 1.0))
42 v_next = CRRAutility(vNvrs_next, rho=rho)
43 dvda = R * dvdm_next
44 dvdH = (1.0 - shock["DeprRte"]) * (shock["WageRte"] * dvdm_next + dvdh_next)
45 return v_next, dvda, dvdH
48###############################################################################
51def solve_one_period_ConsBasicHealth(
52 solution_next,
53 DiscFac,
54 Rfree,
55 CRRA,
56 HealthProdExp,
57 HealthProdFac,
58 DieProbMax,
59 ShockDstn,
60 aLvlGrid,
61 HLvlGrid,
62 constrained_N,
63):
64 """
65 Solve one period of the basic health investment / consumption-saving model
66 using the endogenous grid method. Policy functions are the consumption function
67 cFunc and the health investment function nFunc.
69 Parameters
70 ----------
71 solution_next : Curvilinear2DInterp
72 Solution to the succeeding period's problem, represented as a multi-function
73 interpolant with entries vNvrsFunc, cFunc, and nFunc.
74 DiscFac : float
75 Intertemporal discount factor, representing beta.
76 Rfree : float
77 Risk-free rate of return on retained assets.
78 CRRA : float
79 Coefficient of relative risk aversion, representing rho. Assumed to be
80 constant across periods. Should be strictly between 0 and 1.
81 HealthProdExp : float
82 Exponent in health production function; should be strictly b/w 0 and 1.
83 This corresponds to alpha in White (2015).
84 HealthProdFac : float
85 Scaling factor in health production function; should be strictly positive.
86 This corresponds to gamma in White (2015).
87 DieProbMax : float
88 Maximum death probability at the end of this period, if HLvl were exactly zero.
89 ShockDstn : DiscreteDistribution
90 Joint distribution of income and depreciation values that could realize
91 at the start of the next period.
92 aLvlGrid : np.array
93 Grid of end-of-period assets (after all actions are accomplished).
94 HLvlGrid : np.array
95 Grid of end-of-period post-investment health.
96 constrained_N : int
97 Number of additional interpolation nodes to put in the mLvl dimension
98 on the liquidity-constrained portion of the consumption function.
100 Returns
101 -------
102 solution_now : dict
103 Solution to this period's problem, including policy functions cFunc and
104 nFunc, as well as (marginal) value functions vFunc, dvdmFunc, and dvdhFunc.
105 """
106 # Determine whether there is a liquidity-constrained portion of the policy functions
107 WageRte_min = np.min(ShockDstn.atoms[0])
108 constrained = WageRte_min > 0.0
110 # Adjust the assets grid if liquidity constraint will bind somewhere
111 aLvlGrid_temp = np.insert(aLvlGrid, 0, 0.0) if constrained else aLvlGrid
113 # Make meshes of end-of-period states aLvl and HLvl
114 (aLvl, HLvl) = np.meshgrid(aLvlGrid_temp, HLvlGrid, indexing="ij")
116 # Calculate expected (marginal) value conditional on survival
117 v_next_exp, dvdm_next_exp, dvdh_next_exp = expected(
118 calc_exp_next,
119 ShockDstn,
120 args=(aLvl, HLvl, Rfree, CRRA, HealthProdExp, HealthProdFac, solution_next),
121 )
123 # Calculate (marginal) survival probabilities
124 LivPrb = 1.0 - DieProbMax / (1.0 + HLvl)
125 MargLivPrb = DieProbMax / (1.0 + HLvl) ** 2.0
127 # Calculate end-of-period expectations
128 EndOfPrd_v = DiscFac * (LivPrb * v_next_exp)
129 EndOfPrd_dvda = DiscFac * (LivPrb * dvdm_next_exp)
130 EndOfPrd_dvdH = DiscFac * (LivPrb * dvdh_next_exp + MargLivPrb * v_next_exp)
131 vP_ratio = EndOfPrd_dvda / EndOfPrd_dvdH
133 # Invert the first order conditions to find optimal controls when unconstrained
134 cLvl = EndOfPrd_dvda ** (-1.0 / CRRA)
135 nLvl = (vP_ratio / HealthProdFac) ** (1.0 / (HealthProdExp - 1.0))
137 # If there is a liquidity constrained portion, find additional controls on it
138 if constrained:
139 # Make the grid of constrained health investment by scaling cusp values
140 N = constrained_N # to shorten next line
141 frac_grid = np.reshape(np.linspace(0.01, 0.99, num=N), (N, 1))
142 nLvl_at_cusp = np.reshape(nLvl[0, :], (1, HLvlGrid.size))
143 nLvl_cnst = frac_grid * nLvl_at_cusp
145 # Solve intraperiod FOC to get constrained consumption
146 marg_health_cnst = eval_marg_health_prod(
147 nLvl_cnst, HealthProdExp, HealthProdFac
148 )
149 cLvl_cnst = (EndOfPrd_dvdH[0, :] * marg_health_cnst) ** (-1.0 / CRRA)
151 # Define "constrained end of period states" and continuation value
152 aLvl_cnst = np.zeros((N, HLvlGrid.size))
153 HLvl_cnst = np.tile(np.reshape(HLvlGrid, (1, HLvlGrid.size)), (N, 1))
154 EndOfPrd_v_cnst = np.tile(
155 np.reshape(EndOfPrd_v[0, :], (1, HLvlGrid.size)), (N, 1)
156 )
158 # Combine constrained and unconstrained arrays
159 aLvl = np.concatenate([aLvl_cnst, aLvl], axis=0)
160 HLvl = np.concatenate([HLvl_cnst, HLvl], axis=0)
161 cLvl = np.concatenate([cLvl_cnst, cLvl], axis=0)
162 nLvl = np.concatenate([nLvl_cnst, nLvl], axis=0)
163 EndOfPrd_v = np.concatenate([EndOfPrd_v_cnst, EndOfPrd_v], axis=0)
165 # Invert intratemporal transitions to find endogenous gridpoints
166 mLvl = aLvl + cLvl + nLvl
167 hLvl = HLvl - eval_health_prod(nLvl, HealthProdExp, HealthProdFac)
169 # Calculate (pseudo-inverse) value as of decision-time
170 Value = CRRAutility(cLvl, rho=CRRA) + EndOfPrd_v
171 vNvrs = CRRAutility_inv(Value, rho=CRRA)
173 # Add points at the lower boundary of mLvl for each function
174 Zeros = np.zeros((1, HLvlGrid.size))
175 mLvl = np.concatenate((Zeros, mLvl), axis=0)
176 hLvl = np.concatenate((np.reshape(hLvl[0, :], (1, HLvlGrid.size)), hLvl), axis=0)
177 cLvl = np.concatenate((Zeros, cLvl), axis=0)
178 nLvl = np.concatenate((Zeros, nLvl), axis=0)
179 vNvrs = np.concatenate((Zeros, vNvrs), axis=0)
181 # Construct solution as a multi-interpolation
182 solution_now = Curvilinear2DInterp([vNvrs, cLvl, nLvl], mLvl, hLvl)
183 return solution_now
186###############################################################################
189def make_solution_terminal_ConsBasicHealth():
190 """
191 Constructor for the terminal period solution for the basic health investment
192 model. The trivial solution is to consume all market resources and invest
193 nothing in health. Takes no parameters because CRRA is irrelevant: pseudo-inverse
194 value is returned rather than value, and the former is just cLvl = mLvl.
196 The solution representation for this model is a multiple output function that
197 takes market resources and health capital level as inputs and returns pseudo-
198 inverse value, consumption level, and health investment level in that order.
199 """
200 return lambda mLvl, hLvl: (mLvl, mLvl, np.zeros_like(mLvl))
203def make_health_grid(hLvlMin, hLvlMax, hLvlCount):
204 """
205 Make a uniform grid of health capital levels.
207 Parameters
208 ----------
209 hLvlMin : float
210 Lower bound on health capital level; should almost surely be zero.
211 hLvlMax : float
212 Upper bound on health capital level.
213 hLvlCount : int
214 Number of points in uniform health capital level grid.
216 Returns
217 -------
218 hLvlGrid : np.array
219 Uniform grid of health capital levels
220 """
221 return np.linspace(hLvlMin, hLvlMax, hLvlCount)
224def make_uniform_depreciation_dstn(
225 T_cycle, DeprRteMean, DeprRteSpread, DeprRteCount, RNG
226):
227 """
228 Constructor for DeprRteDstn that makes uniform distributions that vary by age.
230 Parameters
231 ----------
232 T_cycle : int
233 Number of periods in the agent's sequence or cycle.
234 DeprRteMean : [float]
235 Age-varying list (or array) of mean depreciation rates.
236 DeprRteSpread : [float]
237 Age-varying list (or array) of half-widths of depreciate rate distribution.
238 DeprRteCount : int
239 Number of equiprobable nodes in each distribution.
240 RNG : np.random.RandomState
241 Agent's internal random number generator.
243 Returns
244 -------
245 DeprRteDstn : [DiscreteDistribution]
246 List of age-dependent discrete approximations to the depreciate rate distribution.
247 """
248 if len(DeprRteMean) != T_cycle:
249 raise ValueError("DeprRteMean must have length T_cycle!")
250 if len(DeprRteSpread) != T_cycle:
251 raise ValueError("DeprRteSpread must have length T_cycle!")
253 DeprRteDstn = []
254 probs = DeprRteCount**-1.0 * np.ones(DeprRteCount)
255 for t in range(T_cycle):
256 bot = DeprRteMean[t] - DeprRteSpread[t]
257 top = DeprRteMean[t] + DeprRteSpread[t]
258 vals = np.linspace(bot, top, DeprRteCount)
259 DeprRteDstn.append(
260 DiscreteDistribution(
261 pmv=probs,
262 atoms=vals,
263 seed=RNG.integers(0, 2**31 - 1),
264 )
265 )
266 return DeprRteDstn
269def combine_indep_wage_and_depr_dstns(T_cycle, WageRteDstn, DeprRteDstn, RNG):
270 """
271 Combine univariate distributions of wage rate realizations and depreciation
272 rate realizations at each age, treating them as independent.
274 Parameters
275 ----------
276 T_cycle : int
277 Number of periods in the agent's sequence of periods (cycle).
278 WageRteDstn : [DiscreteDistribution]
279 Age-dependent list of wage rate realizations; should have length T_cycle.
280 DeprRteDstn : [DiscreteDistribution]
281 Age-dependent list of health depreciation rate realizatiosn; should have
282 length T_cycle.
283 RNG : np.random.RandomState
284 Internal random number generator for the AgentType instance.
286 Returns
287 -------
288 ShockDstn : [DiscreteDistribution]
289 Age-dependent bivariate distribution with joint realizations of income
290 and health depreciation rates.
291 """
292 if len(WageRteDstn) != T_cycle:
293 raise ValueError(
294 "IncShkDstn must be a list of distributions of length T_cycle!"
295 )
296 if len(DeprRteDstn) != T_cycle:
297 raise ValueError(
298 "DeprRteDstn must be a list of distributions of length T_cycle!"
299 )
301 ShockDstn = []
302 for t in range(T_cycle):
303 temp_dstn = combine_indep_dstns(
304 WageRteDstn[t], DeprRteDstn[t], seed=RNG.integers(0, 2**31 - 1)
305 )
306 temp_dstn_alt = DiscreteDistributionLabeled.from_unlabeled(
307 dist=temp_dstn,
308 name="wage and depreciation shock distribution",
309 var_names=["WageRte", "DeprRte"],
310 )
311 ShockDstn.append(temp_dstn_alt)
312 return ShockDstn
315def make_logistic_polynomial_die_prob(T_cycle, DieProbMaxCoeffs):
316 """
317 Constructor for DieProbMax, the age-varying list of maximum death probabilities
318 (if health is zero). Builds the list as the logistic function evaluated on a
319 polynomial of model age, given polynomial coefficients. Logistic function is
320 applied to ensure probabilities are always between zero and one.
322 Parameters
323 ----------
324 T_cycle : int
325 Number of periods in the agent's sequence of periods (cycle).
326 DieProbMaxCoeffs : np.array
327 List or vector of polynomial coefficients for maximum death probability.
329 Returns
330 -------
331 DieProbMax : [float]
332 Age-varying list of maximum death probabilities (if health were zero).
333 """
334 age_vec = np.arange(T_cycle)
335 DieProbMax = (1.0 + np.exp(-np.polyval(DieProbMaxCoeffs, age_vec))) ** (-1.0)
336 return DieProbMax.tolist()
339def make_uniform_HLvl_init_dstn(HLvlInitMin, HLvlInitMax, HLvlInitCount, RNG):
340 """
341 Constructor for HLvlInitDstn that builds a uniform distribution for initial
342 health capital at model birth.
344 Parameters
345 ----------
346 HLvlInitMin : float
347 Lower bound of initial health capital distribution.
348 HLvlInitMax : float
349 Upper bound of initial health capital distribution
350 HLvlInitCount : int
351 Number of discretized nodes in initial health capital distribution.
352 RNG : np.random.RandomState
353 Agent's internal RNG.
355 Returns
356 -------
357 HLvlInitDstn : DiscreteDistribution
358 Discretized uniform distribution of initial health capital level.
359 """
360 dstn = Uniform(bot=HLvlInitMin, top=HLvlInitMax, seed=RNG.integers(0, 2**31 - 1))
361 HLvlInitDstn = dstn.discretize(HLvlInitCount, endpoints=True)
362 return HLvlInitDstn
365###############################################################################
367# Make a dictionary of default constructor functions
368basic_health_constructors = {
369 "WageRteDstn": construct_lognormal_wage_dstn,
370 "DeprRteDstn": make_uniform_depreciation_dstn,
371 "ShockDstn": combine_indep_wage_and_depr_dstns,
372 "aLvlGrid": make_assets_grid,
373 "HLvlGrid": make_health_grid,
374 "DieProbMax": make_logistic_polynomial_die_prob,
375 "HLvlInitDstn": make_uniform_HLvl_init_dstn,
376 "kLvlInitDstn": make_lognormal_kNrm_init_dstn,
377 "solution_terminal": make_solution_terminal_ConsBasicHealth,
378}
380# Make a dictionary of default parameters for depreciation rate distribution
381default_DeprRteDstn_params = {
382 "DeprRteMean": [0.05], # Mean of uniform depreciation distribution
383 "DeprRteSpread": [0.05], # Half-width of uniform depreciation distribution
384 "DeprRteCount": 7, # Number of nodes in discrete approximation
385}
387# Make a dictionary of default parameters for wage rate distribution
388default_WageRteDstn_params = {
389 "WageRteMean": [0.1], # Age-varying mean of wage rate
390 "WageRteStd": [0.1], # Age-varying stdev of wage rate
391 "WageRteCount": 7, # Number of nodes to use in discrete approximation
392 "UnempPrb": 0.07, # Probability of unemployment
393 "IncUnemp": 0.3, # Income when unemployed
394}
396# Make a dictionary of default parameters for assets grid
397default_aLvlGrid_params = {
398 "aXtraMin": 1e-5, # Minimum value of end-of-period assets grid
399 "aXtraMax": 100.0, # Maximum value of end-of-period assets grid
400 "aXtraCount": 44, # Number of nodes in base assets grid
401 "aXtraNestFac": 1, # Level of exponential nesting for assets grid
402 "aXtraExtra": [3e-5, 1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2], # Extra assets nodes
403}
405# Make a dictionary of default parameters for health capital grid
406default_hLvlGrid_params = {
407 "hLvlMin": 0.0, # Minimum value of health capital grid (leave at zero)
408 "hLvlMax": 50.0, # Maximum value of health capital grid
409 "hLvlCount": 50, # Number of nodes in health capital grid
410}
412# Make a dictionary of default parameters for maximum death probability
413default_DieProbMax_params = {
414 "DieProbMaxCoeffs": [0.0], # Logistic-polynomial coefficients on age
415}
417# Make a dictionary with parameters for the default constructor for kNrmInitDstn
418default_kLvlInitDstn_params = {
419 "kLogInitMean": -12.0, # Mean of log initial capital
420 "kLogInitStd": 0.0, # Stdev of log initial capital
421 "kNrmInitCount": 15, # Number of points in initial capital discretization
422}
424# Make a dictionary with parameters for the default constructor for HLvlInitDstn
425default_HLvlInitDstn_params = {
426 "HLvlInitMin": 1.0, # Lower bound of initial health capital
427 "HLvlInitMax": 2.0, # Upper bound of initial health capital
428 "HLvlInitCount": 15, # Number of points in initial health capital discretization
429}
431# Make a dictionary of default parameters for the health investment model
432basic_health_simple_params = {
433 "constructors": basic_health_constructors,
434 "DiscFac": 0.95, # Intertemporal discount factor
435 "Rfree": [1.03], # Risk-free asset return factor
436 "CRRA": 0.5, # Coefficient of relative risk aversion
437 "HealthProdExp": 0.35, # Exponent on health production function
438 "HealthProdFac": 1.0, # Factor on health production function
439 "constrained_N": 7, # Number of points on liquidity constrained portion
440 "T_cycle": 1, # Number of periods in default cycle
441 "cycles": 1, # Number of cycles
442 "T_age": None, # Maximum lifetime length override
443 "AgentCount": 10000, # Number of agents to simulate
444}
446# Assemble the default parameters dictionary
447init_basic_health = {}
448init_basic_health.update(basic_health_simple_params)
449init_basic_health.update(default_DeprRteDstn_params)
450init_basic_health.update(default_WageRteDstn_params)
451init_basic_health.update(default_aLvlGrid_params)
452init_basic_health.update(default_hLvlGrid_params)
453init_basic_health.update(default_DieProbMax_params)
454init_basic_health.update(default_kLvlInitDstn_params)
455init_basic_health.update(default_HLvlInitDstn_params)
458class BasicHealthConsumerType(AgentType):
459 r"""
460 A class to represent consumers who can save in a risk-free asset and invest
461 in health capital via a health production function. The model is a slight
462 alteration of the one from White (2015), which was in turn lifted from Ludwig
463 and Schoen. In this variation, survival probability depends on post-investment
464 health capital, rather than next period's health capital realization.
466 Each period, the agent chooses consumption $c_t$ and health investment $n_t$.
467 Consumption yields utility via CRRA function, while investment yields additional
468 health capital via production function $f(n_t)$. The agent faces a mortality
469 risk that depends on their post-investment health $H_t = h_t + g(n_t)$, as
470 well as income risk through wage rate $\omega_t$ and health capital depreciation
471 rate $\delta_t$. Health capital also serves as human capital in the sense that
472 the agent earns more income when $h_t$ is higher.
474 Unlike most other HARK models, this one is *not* normalized with respect to
475 permanent income-- indeed, there is no "permanent income" in this simple model.
476 As parametric restrictions, the solver requires that $\rho < 1$ so that utility
477 is positive everywhere. This restriction ensures that the first order conditions
478 are necessary and sufficient to characterize the solution when not liquidity-
479 constrained. The liquidity constrained portion of the policy function is handled.
481 .. math::
482 \newcommand{\CRRA}{\rho}
483 \newcommand{\DiePrb}{\mathsf{D}}
484 \newcommand{\PermGroFac}{\Gamma}
485 \newcommand{\Rfree}{\mathsf{R}}
486 \newcommand{\DiscFac}{\beta}
487 \newcommand{\DeprRte}{\delta}
488 \newcommand{\WageRte}{\omega}
489 \begin{align*}
490 v_t(m_t, h_t) &= \max_{c_t, n_t}u(c_t) + \DiscFac (1 - \DiePrb_{t}) v_{t+1}(m_{t+1}, h_{t+1}), \\
491 & \text{s.t.} \\
492 H_t &= h_t + g(n_t), \\
493 a_t &= m_t - c_t - n_t, \\
494 \DiePrb_t = \phi_t / (1 + H_t), \\
495 h_{t+1} &= (1-\DeprRte_{t+1}) H_t, \\
496 y_{t+1} &= \omega_{t+1} h_{t+1}, \\
497 m_{t+1} &= \Rfree_{t+1} a_t + y_{t+1}, \\
498 u(c) &= \frac{c^{1-\CRRA}}{1-\CRRA}, \\
499 g(n) = (\gamma / \alpha) n^{\alpha}, \\
500 (\WageRte_{t+1}, \DeprRte_{t+1}) \sim F_{t+1}.
501 \end{align*}
503 Solving Parameters
504 ------------------
505 cycles: int
506 0 specifies an infinite horizon model, 1 specifies a finite model.
507 T_cycle: int
508 Number of periods in the cycle for this agent type.
509 CRRA: float, :math:`\rho`
510 Coefficient of Relative Risk Aversion.
511 Rfree: list[float], time varying, :math:`\mathsf{R}`
512 Risk-free interest rate by age.
513 DiscFac: float, :math:`\beta`
514 Intertemporal discount factor.
515 DieProbMax: list[float], time varying, :math:`\phi`
516 Maximum death probability by age, if $H_t=0$.
517 HealthProdExp : float, :math:`\alpha`
518 Exponent in health production function; should be strictly b/w 0 and 1.
519 HealthProdFac : float, :math:`\gamma`
520 Scaling factor in health production function; should be strictly positive.
521 ShockDstn : DiscreteDistribution, time varying
522 Joint distribution of income and depreciation values that could realize
523 at the start of the next period.
524 aLvlGrid : np.array
525 Grid of end-of-period assets (after all actions are accomplished).
526 HLvlGrid : np.array
527 Grid of end-of-period post-investment health.
530 Simulation Parameters
531 ---------------------
532 AgentCount: int
533 Number of agents of this kind that are created during simulations.
534 T_age: int
535 Age after which to automatically kill agents, None to ignore.
536 T_sim: int, required for simulation
537 Number of periods to simulate.
538 track_vars: list[strings]
539 List of variables that should be tracked when running the simulation.
540 For this agent, the options are 'kLvl', 'yLvl', 'mLvl', 'hLvl', 'cLvl',
541 'nLvl', 'WageRte', 'DeprRte', 'aLvl', 'HLvl'.
543 kLvl : Beginning-of-period capital holdings, equivalent to aLvl_{t-1}
545 yLvl : Labor income, the wage rate times health capital.
547 mLvl : Market resources, the interest factor times capital holdings, plus labor income.
549 hLvl : Health or human capital level at decision-time.
551 cLvl : Consumption level.
553 nLvl : Health investment level.
555 WageRte : Wage rate this period.
557 DeprRte : Health capital depreciation rate this period.
559 aLvl : End-of-period assets: market resources less consumption and investment.
561 HLvl : End-of-period health capital: health capital plus produced health.
563 Attributes
564 ----------
565 solution: list[Consumer solution object]
566 Created by the :func:`.solve` method. Finite horizon models create a list with T_cycle+1 elements, for each period in the solution.
567 Infinite horizon solutions return a list with T_cycle elements for each period in the cycle.
569 Visit :class:`HARK.ConsumptionSaving.ConsIndShockModel.ConsumerSolution` for more information about the solution.
570 history: Dict[Array]
571 Created by running the :func:`.simulate()` method.
572 Contains the variables in track_vars. Each item in the dictionary is an array with the shape (T_sim,AgentCount).
573 Visit :class:`HARK.core.AgentType.simulate` for more information.
574 """
576 default_ = {
577 "params": init_basic_health,
578 "solver": solve_one_period_ConsBasicHealth,
579 }
580 time_vary_ = ["Rfree", "DieProbMax", "ShockDstn"]
581 time_inv_ = [
582 "DiscFac",
583 "CRRA",
584 "HealthProdExp",
585 "HealthProdFac",
586 "aLvlGrid",
587 "HLvlGrid",
588 "constrained_N",
589 ]
590 state_vars = ["kLvl", "yLvl", "mLvl", "hLvl", "aLvl", "HLvl"]
591 shock_vars_ = ["WageRte", "DeprRte"]
592 distributions = ["ShockDstn", "kLvlInitDstn", "HLvlInitDstn"]
594 def sim_death(self):
595 """
596 Draw mortality shocks for all agents, marking some for death and replacement.
598 Returns
599 -------
600 which_agents : np.array
601 Boolean array of size AgentCount, indicating who dies now.
602 """
603 # Calculate agent-specific death probability
604 phi = np.array(self.DieProbMax)[self.t_cycle]
605 DieProb = phi / (1.0 + self.state_now["hLvl"])
607 # Draw mortality shocks and mark who dies
608 N = self.AgentCount
609 DeathShks = self.RNG.random(N)
610 which_agents = DeathShks < DieProb
611 if self.T_age is not None: # Kill agents that have lived for too many periods
612 too_old = self.t_age >= self.T_age
613 which_agents = np.logical_or(which_agents, too_old)
614 return which_agents
616 def sim_birth(self, which_agents):
617 """
618 Makes new consumers for the given indices. Initialized variables include
619 kLvl and HLvl, as well as time variables t_age and t_cycle.
621 Parameters
622 ----------
623 which_agents : np.array(Bool)
624 Boolean array of size self.AgentCount indicating which agents should be "born".
626 Returns
627 -------
628 None
629 """
630 N = np.sum(which_agents)
631 kLvl_newborns = self.kLvlInitDstn.draw(N)
632 HLvl_newborns = self.HLvlInitDstn.draw(N)
633 self.state_now["aLvl"][which_agents] = kLvl_newborns
634 self.state_now["HLvl"][which_agents] = HLvl_newborns
635 self.t_age[which_agents] = 0
636 self.t_cycle[which_agents] = 0
638 def get_shocks(self):
639 """
640 Draw wage and depreciation rate shocks for all simulated agents.
641 """
642 WageRte_now = np.empty(self.AgentCount)
643 DeprRte_now = np.empty(self.AgentCount)
644 for t in range(self.T_cycle):
645 these = self.t_cycle == t
646 dstn = self.ShockDstn[t - 1]
647 N = np.sum(these)
648 Shocks = dstn.draw(N)
649 WageRte_now[these] = Shocks[0, :]
650 DeprRte_now[these] = Shocks[1, :]
651 self.shocks["WageRte"] = WageRte_now
652 self.shocks["DeprRte"] = DeprRte_now
654 def transition(self):
655 """
656 Find current market resources and health capital from prior health capital
657 and the drawn shocks.
658 """
659 kLvlNow = self.state_prev["aLvl"]
660 HLvlPrev = self.state_prev["HLvl"]
661 RfreeNow = np.array(self.Rfree)[self.t_cycle - 1]
662 hLvlNow = (1.0 - self.shocks["DeprRte"]) * HLvlPrev
663 yLvlNow = self.shocks["WageRte"] * hLvlNow
664 mLvlNow = RfreeNow * kLvlNow + yLvlNow
665 return kLvlNow, yLvlNow, mLvlNow, hLvlNow
667 def get_controls(self):
668 """
669 Evaluate consumption and health investment functions conditional on
670 current state and model age, yielding controls cLvl and nLvl.
671 """
672 # This intentionally has the same bug with cycles > 1 as all our other
673 # models. It will be fixed all in one PR.
674 mLvl = self.state_now["mLvl"]
675 hLvl = self.state_now["hLvl"]
676 cLvl = np.empty(self.AgentCount)
677 nLvl = np.empty(self.AgentCount)
678 for t in range(self.T_cycle):
679 these = self.t_cycle == t
680 func_t = self.solution[t]
681 trash, cLvl[these], nLvl[these] = func_t(mLvl[these], hLvl[these])
682 self.controls["cLvl"] = cLvl
683 self.controls["nLvl"] = nLvl
685 def get_poststates(self):
686 """
687 Calculate end-of-period retained assets and post-investment health.
688 """
689 self.state_now["aLvl"] = (
690 self.state_now["mLvl"] - self.controls["cLvl"] - self.controls["nLvl"]
691 )
692 self.state_now["HLvl"] = (
693 self.state_now["hLvl"]
694 + (self.HealthProdFac / self.HealthProdExp)
695 * self.controls["nLvl"] ** self.HealthProdExp
696 )