Coverage for HARK / Calibration / Income / IncomeProcesses.py: 90%
375 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-08 05:31 +0000
1"""
2This file has various classes and functions for constructing income processes.
3"""
5import numpy as np
6from scipy.stats import norm
7from HARK.metric import MetricObject
8from HARK.distributions import (
9 add_discrete_outcome,
10 add_discrete_outcome_constant_mean,
11 combine_indep_dstns,
12 DiscreteDistribution,
13 DiscreteDistributionLabeled,
14 IndexDistribution,
15 MeanOneLogNormal,
16 Lognormal,
17 Uniform,
18 make_tauchen_ar1,
19)
20from HARK.interpolation import IdentityFunction, LinearInterp
21from HARK.utilities import get_percentiles, make_polynomial_params
24class BinaryIncShkDstn(DiscreteDistribution):
25 """
26 A one period income shock distribution (transitory, permanent, or other)
27 with only two outcomes. One probability and value are specified, and the
28 other is implied to make it a mean one distribution.
30 Parameters
31 ----------
32 shk_prob : float
33 Probability of one of the income shock outcomes.
34 shk_val : float
35 Value of the specified income shock outcome.
36 seed : int, optional
37 Random seed. The default is 0.
39 Returns
40 -------
41 ShkDstn : DiscreteDistribution
42 Binary income shock distribuion.
43 """
45 def __init__(self, shk_prob, shk_val, seed=0):
46 if shk_prob > 1.0 or shk_prob < 0.0:
47 raise ValueError("Shock probability must be between 0 and 1!")
49 other_prob = 1.0 - shk_prob
50 other_val = (1.0 - shk_prob * shk_val) / other_prob
51 probs = [shk_prob, other_prob]
52 vals = [shk_val, other_val]
53 super().__init__(pmv=probs, atoms=vals, seed=seed)
56class LognormPermIncShk(DiscreteDistribution):
57 """
58 A one-period distribution of a multiplicative lognormal permanent income shock.
59 Parameters
60 ----------
61 sigma : float
62 Standard deviation of the log-shock.
63 n_approx : int
64 Number of points to use in the discrete approximation.
65 neutral_measure : Bool, optional
66 Whether to use Harmenberg's permanent-income-neutral measure. The default is False.
67 seed : int, optional
68 Random seed. The default is 0.
70 Returns
71 -------
72 PermShkDstn : DiscreteDistribution
73 Permanent income shock distribution.
75 """
77 def __init__(self, sigma, n_approx, neutral_measure=False, seed=0):
78 # Construct an auxiliary discretized normal
79 lognormal_dstn = MeanOneLogNormal(sigma)
80 logn_approx = lognormal_dstn.discretize(
81 n_approx if sigma > 0.0 else 1, method="equiprobable", tail_N=0
82 )
84 limit = {
85 "dist": lognormal_dstn,
86 "method": "equiprobable",
87 "N": n_approx,
88 "endpoints": False,
89 "infimum": logn_approx.limit["infimum"],
90 "supremum": logn_approx.limit["supremum"],
91 }
93 # Change the pmv if necessary
94 if neutral_measure:
95 logn_approx.pmv = (logn_approx.atoms * logn_approx.pmv).flatten()
97 super().__init__(
98 pmv=logn_approx.pmv, atoms=logn_approx.atoms, limit=limit, seed=seed
99 )
102class MixtureTranIncShk(DiscreteDistribution):
103 """
104 A one-period distribution for transitory income shocks that are a mixture
105 between a log-normal and a single-value unemployment shock.
107 Parameters
108 ----------
109 sigma : float
110 Standard deviation of the log-shock.
111 UnempPrb : float
112 Probability of the "unemployment" shock.
113 IncUnemp : float
114 Income shock in the "unemployment" state.
115 n_approx : int
116 Number of points to use in the discrete approximation.
117 seed : int, optional
118 Random seed. The default is 0.
120 Returns
121 -------
122 TranShkDstn : DiscreteDistribution
123 Transitory income shock distribution.
125 """
127 def __init__(self, sigma, UnempPrb, IncUnemp, n_approx, seed=0):
128 dstn_orig = MeanOneLogNormal(sigma)
129 dstn_approx = dstn_orig.discretize(
130 n_approx if sigma > 0.0 else 1, method="equiprobable", tail_N=0
131 )
133 if UnempPrb > 0.0:
134 dstn_approx = add_discrete_outcome_constant_mean(
135 dstn_approx, p=UnempPrb, x=IncUnemp
136 )
138 super().__init__(
139 pmv=dstn_approx.pmv,
140 atoms=dstn_approx.atoms,
141 limit=dstn_approx.limit,
142 seed=seed,
143 )
146class MixtureTranIncShk_HANK(DiscreteDistribution):
147 """
148 A one-period distribution for transitory income shocks that are a mixture
149 between a log-normal and a single-value unemployment shock. This version
150 has additional parameters that makes it useful for HANK models.
152 Parameters
153 ----------
154 sigma : float
155 Standard deviation of the log-shock.
156 UnempPrb : float
157 Probability of the "unemployment" shock.
158 IncUnemp : float
159 Income shock in the "unemployment" state.
160 n_approx : int
161 Number of points to use in the discrete approximation.
162 tax_rate : float
163 Flat tax rate on labor income.
164 labor : float
165 Intensive margin labor supply.
166 wage : float
167 Wage rate scaling factor.
168 seed : int, optional
169 Random seed. The default is 0.
170 Returns
171 -------
172 TranShkDstn : DiscreteDistribution
173 Transitory income shock distribution.
174 """
176 def __init__(
177 self,
178 sigma,
179 UnempPrb,
180 IncUnemp,
181 n_approx,
182 wage,
183 labor,
184 tax_rate,
185 seed=0,
186 ):
187 dstn_approx = MeanOneLogNormal(sigma).discretize(
188 n_approx if sigma > 0.0 else 1, method="equiprobable", tail_N=0
189 )
191 if UnempPrb > 0.0:
192 dstn_approx = add_discrete_outcome_constant_mean(
193 dstn_approx, p=UnempPrb, x=IncUnemp
194 )
195 # Rescale the transitory shock values to account for new features
196 TranShkMean_temp = (1.0 - tax_rate) * labor * wage
197 dstn_approx.atoms *= TranShkMean_temp
198 super().__init__(
199 pmv=dstn_approx.pmv,
200 atoms=dstn_approx.atoms,
201 limit=dstn_approx.limit,
202 seed=seed,
203 )
206class BufferStockIncShkDstn(DiscreteDistributionLabeled):
207 """
208 A one-period distribution object for the joint distribution of income
209 shocks (permanent and transitory), as modeled in the Buffer Stock Theory
210 paper:
211 - Lognormal, discretized permanent income shocks.
212 - Transitory shocks that are a mixture of:
213 - A lognormal distribution in normal times.
214 - An "unemployment" shock.
216 Parameters
217 ----------
218 sigma_Perm : float
219 Standard deviation of the log- permanent shock.
220 sigma_Tran : float
221 Standard deviation of the log- transitory shock.
222 n_approx_Perm : int
223 Number of points to use in the discrete approximation of the permanent shock.
224 n_approx_Tran : int
225 Number of points to use in the discrete approximation of the transitory shock.
226 UnempPrb : float
227 Probability of the "unemployment" shock.
228 IncUnemp : float
229 Income shock in the "unemployment" state.
230 neutral_measure : Bool, optional
231 Whether to use Hamenberg's permanent-income-neutral measure. The default is False.
232 seed : int, optional
233 Random seed. The default is 0.
235 Returns
236 -------
237 IncShkDstn : DiscreteDistribution
238 Income shock distribution.
240 """
242 def __init__(
243 self,
244 sigma_Perm,
245 sigma_Tran,
246 n_approx_Perm,
247 n_approx_Tran,
248 UnempPrb,
249 IncUnemp,
250 neutral_measure=False,
251 seed=0,
252 ):
253 perm_dstn = LognormPermIncShk(
254 sigma=sigma_Perm, n_approx=n_approx_Perm, neutral_measure=neutral_measure
255 )
256 tran_dstn = MixtureTranIncShk(
257 sigma=sigma_Tran,
258 UnempPrb=UnempPrb,
259 IncUnemp=IncUnemp,
260 n_approx=n_approx_Tran,
261 )
262 joint_dstn = combine_indep_dstns(perm_dstn, tran_dstn)
264 super().__init__(
265 name="Joint distribution of permanent and transitory shocks to income",
266 var_names=["PermShk", "TranShk"],
267 pmv=joint_dstn.pmv,
268 atoms=joint_dstn.atoms,
269 limit=joint_dstn.limit,
270 seed=seed,
271 )
274class IncShkDstn_HANK(DiscreteDistributionLabeled):
275 """
276 A one-period distribution object for the joint distribution of income
277 shocks (permanent and transitory), as modeled in the Buffer Stock Theory
278 paper:
279 - Lognormal, discretized permanent income shocks.
280 - Transitory shocks that are a mixture of:
281 - A lognormal distribution in normal times.
282 - An "unemployment" shock.
284 This version has additional features that make it particularly useful for HANK models.
286 Parameters
287 ----------
288 sigma_Perm : float
289 Standard deviation of the log- permanent shock.
290 sigma_Tran : float
291 Standard deviation of the log- transitory shock.
292 n_approx_Perm : int
293 Number of points to use in the discrete approximation of the permanent shock.
294 n_approx_Tran : int
295 Number of points to use in the discrete approximation of the transitory shock.
296 UnempPrb : float
297 Probability of the "unemployment" shock.
298 IncUnemp : float
299 Income shock in the "unemployment" state.
300 tax_rate : float
301 Flat tax rate on labor income.
302 labor : float
303 Intensive margin labor supply.
304 wage : float
305 Wage rate scaling factor.
306 neutral_measure : Bool, optional
307 Whether to use Harmenberg's permanent-income-neutral measure. The default is False.
308 seed : int, optional
309 Random seed. The default is 0.
310 Returns
311 -------
312 IncShkDstn : DiscreteDistribution
313 Income shock distribution.
314 """
316 def __init__(
317 self,
318 sigma_Perm,
319 sigma_Tran,
320 n_approx_Perm,
321 n_approx_Tran,
322 UnempPrb,
323 IncUnemp,
324 tax_rate,
325 labor,
326 wage,
327 neutral_measure=False,
328 seed=0,
329 ):
330 perm_dstn = LognormPermIncShk(
331 sigma=sigma_Perm, n_approx=n_approx_Perm, neutral_measure=neutral_measure
332 )
333 tran_dstn = MixtureTranIncShk_HANK(
334 sigma=sigma_Tran,
335 UnempPrb=UnempPrb,
336 IncUnemp=IncUnemp,
337 n_approx=n_approx_Tran,
338 wage=wage,
339 labor=labor,
340 tax_rate=tax_rate,
341 )
342 joint_dstn = combine_indep_dstns(perm_dstn, tran_dstn)
344 super().__init__(
345 name="HANK",
346 var_names=["PermShk", "TranShk"],
347 pmv=joint_dstn.pmv,
348 atoms=joint_dstn.atoms,
349 limit=joint_dstn.limit,
350 seed=seed,
351 )
354###############################################################################
357def construct_lognormal_wage_dstn(
358 T_cycle, WageRteMean, WageRteStd, WageRteCount, IncUnemp, UnempPrb, RNG
359):
360 """
361 Constructor for an age-dependent wage rate distribution. The distribution
362 at each age is (equiprobably discretized) lognormal with a point mass to
363 represent unemployment. This is effectively a "transitory only" income process.
365 Parameters
366 ----------
367 T_cycle : int
368 Number of periods in the agent's cycle or sequence.
369 WageRteMean : [float]
370 Age-varying list (or array) of mean wage rates.
371 WageRteStd : [float]
372 Age-varying standard deviations of (log) wage rates.
373 WageRteCount : int
374 Number of equiprobable nodes in the lognormal approximation.
375 UnempPrb : [float] or float
376 Age-varying probability of unemployment; can be specified to be constant.
377 IncUnemp : [float] or float
378 Age-varying "wage" rate when unemployed, maybe representing benefits.
379 Can be specified to be constant.
380 RNG : np.random.RandomState
381 Agent's internal random number generator.
383 Returns
384 -------
385 WageRteDstn : [DiscreteDistribution]
386 Age-varying list of discrete approximations to the lognormal wage distribution.
387 """
388 if len(WageRteMean) != T_cycle:
389 raise ValueError("WageRteMean must be a list of length T_cycle!")
390 if len(WageRteStd) != T_cycle:
391 raise ValueError("WageRteStd must be a list of length T_cycle!")
392 if not (isinstance(UnempPrb, float) or len(UnempPrb) == T_cycle):
393 raise ValueError("UnempPrb must be a single value or list of length T_cycle!")
394 if not (isinstance(IncUnemp, float) or len(IncUnemp) == T_cycle):
395 raise ValueError("IncUnemp must be a single value or list of length T_cycle!")
397 WageRteDstn = []
398 N = WageRteCount # lazy typing
399 for t in range(T_cycle):
400 # Get current period values
401 W_sig = WageRteStd[t]
402 W_mu = np.log(WageRteMean[t]) - 0.5 * W_sig**2
403 B = IncUnemp if isinstance(IncUnemp, float) else IncUnemp[t]
404 U = UnempPrb if isinstance(UnempPrb, float) else UnempPrb[t]
405 temp_dstn = Lognormal(mu=W_mu, sigma=W_sig, seed=RNG.integers(0, 2**31 - 1))
406 temp_dstn_alt = add_discrete_outcome(temp_dstn.discretize(N), B, U)
407 WageRteDstn.append(temp_dstn_alt)
408 return WageRteDstn
411def construct_lognormal_income_process_unemployment(
412 T_cycle,
413 PermShkStd,
414 PermShkCount,
415 TranShkStd,
416 TranShkCount,
417 T_retire,
418 UnempPrb,
419 IncUnemp,
420 UnempPrbRet,
421 IncUnempRet,
422 RNG,
423 neutral_measure=False,
424):
425 r"""
426 Generates a list of discrete approximations to the income process for each
427 life period, from end of life to beginning of life. Permanent shocks (:math:`\psi`) are mean
428 one lognormally distributed with standard deviation PermShkStd[t] during the
429 working life, and degenerate at 1 in the retirement period. Transitory shocks (:math:`\theta`)
430 are mean one lognormally distributed with a point mass at IncUnemp with
431 probability UnempPrb while working; they are mean one with a point mass at
432 IncUnempRet with probability UnempPrbRet. Retirement occurs
433 after t=T_retire periods of working.
435 .. math::
436 \begin{align*}
437 \psi_t &\sim \begin{cases}
438 \exp(\mathcal{N}(-\textbf{PermShkStd}_{t}^{2}/2,\textbf{PermShkStd}_{t}^{2})) & \text{if } t \leq t_{\text{retire}}\\
439 1 & \text{if } t > t_{\text{retire}}
440 \end{cases}\\
441 p_{\text{unemp}} & = \begin{cases}
442 \textbf{UnempPrb} & \text{if } t \leq t_{\text{retire}} \\
443 \textbf{UnempPrbRet} & \text{if } t > t_{\text{retire}} \\
444 \end{cases}\\
445 &\text{if } p > p_{\text{unemp}} \\
446 \theta_t &\sim\begin{cases}
447 \exp(\mathcal{N}(-\textbf{PermShkStd}_{t}^{2}/2-\ln(\frac{1-\textbf{IncUnemp }\textbf{UnempPrb}}{1-\textbf{UnempPrb}}),\textbf{PermShkStd}_{t}^{2})) & \text{if } t\leq t_{\text{retire}}\\
448 \frac{1-\textbf{UnempPrbRet }\textbf{IncUnempRet}}{1-\textbf{UnempPrbRet}} & \text{if } t > t_{\text{retire}} \\
449 \end{cases}\\
450 &\text{otherwise}\\
451 \theta_t &\sim\begin{cases}
452 \textbf{IncUnemp} & \text{if } t\leq t_{\text{retire}}\\
453 \textbf{IncUnempRet} & \text{if } t\leq t_{\text{retire}}\\
454 \end{cases}\\
455 \mathbb{E}[\psi]&=\mathbb{E}[\theta] = 1.\\
456 \end{align*}
458 All time in this function runs forward, from t=0 to t=T
460 Parameters
461 ----------
462 PermShkStd : [float]
463 List of standard deviations in log permanent income uncertainty during
464 the agent's life.
465 PermShkCount : int
466 The number of approximation points to be used in the discrete approximation
467 to the permanent income shock distribution.
468 TranShkStd : [float]
469 List of standard deviations in log transitory income uncertainty during
470 the agent's life.
471 TranShkCount : int
472 The number of approximation points to be used in the discrete approximation
473 to the permanent income shock distribution.
474 UnempPrb : float or [float]
475 The probability of becoming unemployed during the working period.
476 UnempPrbRet : float or None
477 The probability of not receiving typical retirement income when retired.
478 T_retire : int
479 The index value for the final working period in the agent's life.
480 If T_retire <= 0 then there is no retirement.
481 IncUnemp : float or [float]
482 Transitory income received when unemployed.
483 IncUnempRet : float or None
484 Transitory income received while "unemployed" when retired.
485 T_cycle : int
486 Total number of non-terminal periods in the consumer's sequence of periods.
487 RNG : np.random.RandomState
488 Random number generator for this type.
489 neutral_measure : bool
490 Indicator for whether the permanent-income-neutral measure should be used.
492 Returns
493 -------
494 IncShkDstn : [distribution.Distribution]
495 A list with T_cycle elements, each of which is a
496 discrete approximation to the income process in a period.
497 """
498 if T_retire > 0:
499 normal_length = T_retire
500 retire_length = T_cycle - T_retire
501 else:
502 normal_length = T_cycle
503 retire_length = 0
505 if all(
506 [
507 isinstance(x, (float, int)) or (x is None)
508 for x in [UnempPrb, IncUnemp, UnempPrbRet, IncUnempRet]
509 ]
510 ):
511 UnempPrb_list = [UnempPrb] * normal_length + [UnempPrbRet] * retire_length
512 IncUnemp_list = [IncUnemp] * normal_length + [IncUnempRet] * retire_length
514 elif all([isinstance(x, list) for x in [UnempPrb, IncUnemp]]):
515 UnempPrb_list = UnempPrb
516 IncUnemp_list = IncUnemp
518 else:
519 raise Exception(
520 "Unemployment must be specified either using floats for UnempPrb,"
521 + "IncUnemp, UnempPrbRet, and IncUnempRet, in which case the "
522 + "unemployment probability and income change only with retirement, or "
523 + "using lists of length T_cycle for UnempPrb and IncUnemp, specifying "
524 + "each feature at every age."
525 )
527 PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
528 TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
529 neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
531 IncShkDstn = IndexDistribution(
532 engine=BufferStockIncShkDstn,
533 conditional={
534 "sigma_Perm": PermShkStd,
535 "sigma_Tran": TranShkStd,
536 "n_approx_Perm": PermShkCount_list,
537 "n_approx_Tran": TranShkCount_list,
538 "neutral_measure": neutral_measure_list,
539 "UnempPrb": UnempPrb_list,
540 "IncUnemp": IncUnemp_list,
541 },
542 RNG=RNG,
543 seed=RNG.integers(0, 2**31 - 1),
544 )
545 return IncShkDstn
548def construct_markov_lognormal_income_process_unemployment(
549 T_cycle,
550 PermShkStd,
551 PermShkCount,
552 TranShkStd,
553 TranShkCount,
554 T_retire,
555 UnempPrb,
556 IncUnemp,
557 UnempPrbRet,
558 IncUnempRet,
559 RNG,
560 neutral_measure=False,
561):
562 """
563 Generates a nested list of discrete approximations to the income process for each
564 life period, from end of life to beginning of life, for each discrete Markov
565 state in the problem. This function calls construct_lognormal_income_process_unemployment
566 for each Markov state, then rearranges the output. Permanent shocks are mean
567 one lognormally distributed with standard deviation PermShkStd[t] during the
568 working life, and degenerate at 1 in the retirement period. Transitory shocks
569 are mean one lognormally distributed with a point mass at IncUnemp with
570 probability UnempPrb while working; they are mean one with a point mass at
571 IncUnempRet with probability UnempPrbRet. Retirement occurs after t=T_retire
572 periods of working. The problem is specified as having T=T_cycle periods and
573 K discrete Markov states.
575 Parameters
576 ----------
577 PermShkStd : np.array
578 2D array of shape (T,K) of standard deviations of log permanent income
579 uncertainty during the agent's life.
580 PermShkCount : int
581 The number of approximation points to be used in the discrete approxima-
582 tion to the permanent income shock distribution.
583 TranShkStd : np.array
584 2D array of shape (T,K) standard deviations in log transitory income
585 uncertainty during the agent's life.
586 TranShkCount : int
587 The number of approximation points to be used in the discrete approxima-
588 tion to the permanent income shock distribution.
589 UnempPrb : np.array
590 1D or 2D array of the probability of becoming unemployed during the working
591 portion of the agent's life. If 1D, it should be size K, providing one
592 unemployment probability per discrete Markov state. If 2D, it should be
593 shape (T,K).
594 UnempPrbRet : np.array or None
595 The probability of not receiving typical retirement income when retired.
596 If not None, should be size K. Can be set to None when T_retire is 0.
597 T_retire : int
598 The index value for the final working period in the agent's life.
599 If T_retire <= 0 then there is no retirement.
600 IncUnemp : np.array
601 1D or 2D array of transitory income received when unemployed during the
602 working portion of the agent's life. If 1D, it should be size K, providing
603 one unemployment probability per discrete Markov state. If 2D, it should be
604 shape (T,K).
605 IncUnempRet : np.array or None
606 Transitory income received while "unemployed" when retired. If provided,
607 should be size K. Can be None when T_retire is 0.
608 T_cycle : int
609 Total number of non-terminal periods in the consumer's sequence of periods.
610 RNG : np.random.RandomState
611 Random number generator for this type.
612 neutral_measure : bool
613 Indicator for whether the permanent-income-neutral measure should be used.
615 Returns
616 -------
617 IncShkDstn : [[distribution.Distribution]]
618 A list with T_cycle elements, each of which is a discrete approximation
619 to the income process in a period.
620 """
621 if T_retire > 0:
622 normal_length = T_retire
623 retire_length = T_cycle - T_retire
624 else:
625 normal_length = T_cycle
626 retire_length = 0
628 # Check dimensions of inputs
629 try:
630 PermShkStd_K = PermShkStd.shape[1]
631 TranShkStd_K = TranShkStd.shape[1]
632 if UnempPrb.ndim == 2:
633 UnempPrb_K = UnempPrb.shape[1]
634 else:
635 UnempPrb_K = UnempPrb.shape[0]
636 if IncUnemp.ndim == 2:
637 IncUnemp_K = IncUnemp.shape[1]
638 else:
639 IncUnemp_K = IncUnemp.shape[0]
640 K = PermShkStd_K
641 assert K == TranShkStd_K
642 assert K == UnempPrb_K
643 assert K == IncUnemp_K
644 except:
645 raise Exception(
646 "The last dimension of PermShkStd, TranShkStd, IncUnemp,"
647 + " and UnempPrb must all be K, the number of discrete states."
648 )
649 try:
650 if T_retire > 0:
651 assert K == UnempPrbRet.size
652 assert K == IncUnempRet.size
653 except:
654 raise Exception(
655 "When T_retire is not zero, UnempPrbRet and IncUnempRet"
656 + " must be specified as arrays of size K, the number of "
657 + "discrete Markov states."
658 )
659 try:
660 D = UnempPrb.ndim
661 assert D == IncUnemp.ndim
662 if T_retire > 0:
663 assert D == UnempPrbRet.ndim
664 assert D == IncUnempRet.ndim
665 except:
666 raise Exception(
667 "If any of UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet "
668 + "are 2D arrays, then they must *all* be 2D arrays."
669 )
670 try:
671 assert D == 1 or D == 2
672 except:
673 raise Exception(
674 "UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet must "
675 + "all be 1D or 2D arrays."
676 )
678 # Prepare lists that don't vary by Markov state
679 PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
680 TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
681 neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
683 # Loop through the Markov states, constructing the lifecycle income process for each one
684 IncShkDstn_by_Mrkv = []
685 for k in range(K):
686 if D == 1: # Unemployment parameters don't vary by age other than retirement
687 if T_retire > 0:
688 UnempPrb_list = [UnempPrb[k]] * normal_length + [
689 UnempPrbRet[k]
690 ] * retire_length
691 IncUnemp_list = [IncUnemp[k]] * normal_length + [
692 IncUnempRet[k]
693 ] * retire_length
694 else:
695 UnempPrb_list = [UnempPrb[k]] * normal_length
696 IncUnemp_list = [IncUnemp[k]] * normal_length
697 else: # Unemployment parameters vary by age
698 UnempPrb_list = UnempPrb[:, k].tolist()
699 IncUnemp_list = IncUnemp[:, k].tolist()
701 PermShkStd_k = PermShkStd[:, k].tolist()
702 TranShkStd_k = TranShkStd[:, k].tolist()
704 IncShkDstn_k = IndexDistribution(
705 engine=BufferStockIncShkDstn,
706 conditional={
707 "sigma_Perm": PermShkStd_k,
708 "sigma_Tran": TranShkStd_k,
709 "n_approx_Perm": PermShkCount_list,
710 "n_approx_Tran": TranShkCount_list,
711 "neutral_measure": neutral_measure_list,
712 "UnempPrb": UnempPrb_list,
713 "IncUnemp": IncUnemp_list,
714 },
715 RNG=RNG,
716 seed=RNG.integers(0, 2**31 - 1),
717 )
718 IncShkDstn_by_Mrkv.append(IncShkDstn_k)
720 # Rearrange the list that was just constructed, so that element [t][k] represents
721 # the income shock distribution in period t, Markov state k.
722 IncShkDstn = []
723 for t in range(T_cycle):
724 IncShkDstn_t = [IncShkDstn_by_Mrkv[k][t] for k in range(K)]
725 IncShkDstn.append(IncShkDstn_t)
727 return IncShkDstn
730def construct_HANK_lognormal_income_process_unemployment(
731 T_cycle,
732 PermShkStd,
733 PermShkCount,
734 TranShkStd,
735 TranShkCount,
736 T_retire,
737 UnempPrb,
738 IncUnemp,
739 UnempPrbRet,
740 IncUnempRet,
741 tax_rate,
742 labor,
743 wage,
744 RNG,
745 neutral_measure=False,
746):
747 """
748 Generates a list of discrete approximations to the income process for each
749 life period, from end of life to beginning of life. Permanent shocks are mean
750 one lognormally distributed with standard deviation PermShkStd[t] during the
751 working life, and degenerate at 1 in the retirement period. Transitory shocks
752 are mean one lognormally distributed with a point mass at IncUnemp with
753 probability UnempPrb while working; they are mean one with a point mass at
754 IncUnempRet with probability UnempPrbRet. Retirement occurs after t=T_retire
755 periods of working.
757 This version of the function incorporates additional flexibility with respect
758 to transitory income (continuous labor supply, wage rate, tax rate) and thus
759 is useful in HANK models (hence the name!).
761 Note 1: All time in this function runs forward, from t=0 to t=T
763 Note 2: All parameters are passed as attributes of the input parameters.
765 Parameters (passed as attributes of the input parameters)
766 ---------------------------------------------------------
767 PermShkStd : [float]
768 List of standard deviations in log permanent income uncertainty during
769 the agent's life.
770 PermShkCount : int
771 The number of approximation points to be used in the discrete approxima-
772 tion to the permanent income shock distribution.
773 TranShkStd : [float]
774 List of standard deviations in log transitory income uncertainty during
775 the agent's life.
776 TranShkCount : int
777 The number of approximation points to be used in the discrete approxima-
778 tion to the permanent income shock distribution.
779 UnempPrb : float or [float]
780 The probability of becoming unemployed during the working period.
781 UnempPrbRet : float or None
782 The probability of not receiving typical retirement income when retired.
783 T_retire : int
784 The index value for the final working period in the agent's life.
785 If T_retire <= 0 then there is no retirement.
786 IncUnemp : float or [float]
787 Transitory income received when unemployed.
788 IncUnempRet : float or None
789 Transitory income received while "unemployed" when retired.
790 tax_rate : [float]
791 List of flat tax rates on labor income, age-varying.
792 labor : [float]
793 List of intensive margin labor supply, age-varying.
794 wage : [float]
795 List of wage rate scaling factors, age-varying.
796 T_cycle : int
797 Total number of non-terminal periods in the consumer's sequence of periods.
799 Returns
800 -------
801 IncShkDstn : [distribution.Distribution]
802 A list with T_cycle elements, each of which is a
803 discrete approximation to the income process in a period.
804 PermShkDstn : [[distribution.Distributiony]]
805 A list with T_cycle elements, each of which
806 a discrete approximation to the permanent income shocks.
807 TranShkDstn : [[distribution.Distribution]]
808 A list with T_cycle elements, each of which
809 a discrete approximation to the transitory income shocks.
810 """
811 if T_retire > 0:
812 normal_length = T_retire
813 retire_length = T_cycle - T_retire
814 else:
815 normal_length = T_cycle
816 retire_length = 0
818 if all(
819 [
820 isinstance(x, (float, int)) or (x is None)
821 for x in [UnempPrb, IncUnemp, UnempPrbRet, IncUnempRet]
822 ]
823 ):
824 UnempPrb_list = [UnempPrb] * normal_length + [UnempPrbRet] * retire_length
825 IncUnemp_list = [IncUnemp] * normal_length + [IncUnempRet] * retire_length
827 elif all([isinstance(x, list) for x in [UnempPrb, IncUnemp]]):
828 UnempPrb_list = UnempPrb
829 IncUnemp_list = IncUnemp
831 else:
832 raise Exception(
833 "Unemployment must be specified either using floats for UnempPrb,"
834 + "IncUnemp, UnempPrbRet, and IncUnempRet, in which case the "
835 + "unemployment probability and income change only with retirement, or "
836 + "using lists of length T_cycle for UnempPrb and IncUnemp, specifying "
837 + "each feature at every age."
838 )
840 PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
841 TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
842 neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
844 IncShkDstn = IndexDistribution(
845 engine=IncShkDstn_HANK,
846 conditional={
847 "sigma_Perm": PermShkStd,
848 "sigma_Tran": TranShkStd,
849 "n_approx_Perm": PermShkCount_list,
850 "n_approx_Tran": TranShkCount_list,
851 "neutral_measure": neutral_measure_list,
852 "UnempPrb": UnempPrb_list,
853 "IncUnemp": IncUnemp_list,
854 "wage": wage,
855 "tax_rate": tax_rate,
856 "labor": labor,
857 },
858 RNG=RNG,
859 )
861 return IncShkDstn
864def construct_lognormal_income_process_with_mvg_medical_expenses(
865 T_cycle,
866 PermShkStd,
867 PermShkCount,
868 TranShkStd,
869 TranShkCount,
870 T_retire,
871 UnempPrb,
872 IncUnemp,
873 UnempPrbRet,
874 IncUnempRet,
875 age_min,
876 RNG,
877 CollegeBool=True,
878):
879 """
880 Extension of the standard permanent-transitory income process that adds medical
881 expense shocks, modeled as *negative* transitory income shocks. That is,
882 TranShk = 1.0 - ExpShk. The expense shock calibration is "hardwired" based on
883 table 12 of Mateo Velasquez-Giraldo's "Life-Cycle Portfolio Choices and Heterogeneous
884 Stock Market Expectations" (https://www.federalreserve.gov/econres/feds/files/2024097pap.pdf).
886 Expense shocks are only active after age 50, and are specified as seven equi-
887 probable point masses for each five-year age range. Because the expense shocks
888 are simply "negative transitory income shocks", they are homothetic with respect
889 to permanent income. This is not a particularly justified assumption, as income
890 elasticity of demand for medical care is around 0.15, but it's workable for a
891 representation of medical expenses that can be used with HARK's workhorse models.
892 This function assumes an annual frequency. If net income would ever be negative,
893 it is truncated to zero instead.
895 For microeconomic models that incorporate medical expenses that are not homothetic
896 in income, see ConsMedModel.
898 Relative to the standard construct_lognormal_income_process_unemployment, this
899 function takes two additional parameters: age_min and CollegeBool (default True).
900 This function assumes that the calibration does not include ages above 120 years.
902 Parameters
903 ----------
904 T_cycle : int
905 Total number of non-terminal periods in the consumer's sequence of periods.
906 PermShkStd : [float]
907 List of standard deviations in log permanent income uncertainty during
908 the agent's life.
909 PermShkCount : int
910 The number of approximation points to be used in the discrete approximation
911 to the permanent income shock distribution.
912 TranShkStd : [float]
913 List of standard deviations in log transitory income uncertainty during
914 the agent's life.
915 TranShkCount : int
916 The number of approximation points to be used in the discrete approximation
917 to the permanent income shock distribution.
918 UnempPrb : float or [float]
919 The probability of becoming unemployed during the working period.
920 UnempPrbRet : float or None
921 The probability of not receiving typical retirement income when retired.
922 T_retire : int
923 The index value for the final working period in the agent's life.
924 If T_retire <= 0 then there is no retirement.
925 IncUnemp : float or [float]
926 Transitory income received when unemployed.
927 IncUnempRet : float or None
928 Transitory income received while "unemployed" when retired.
929 age_min : int
930 Age (in years) of agents at model birth.
931 RNG : np.random.RandomState
932 Random number generator for this type.
933 CollegeBool : bool
934 Indicator for whether to use the college (True, default) or high school
935 (False) calibration.
937 Returns
938 -------
939 IncShkDstn : [distribution.Distribution]
940 A list with T_cycle elements, each of which is a
941 discrete approximation to the income process in a period.
942 """
943 # First, make the standard income process
944 IncShkDstnBase = construct_lognormal_income_process_unemployment(
945 T_cycle,
946 PermShkStd,
947 PermShkCount,
948 TranShkStd,
949 TranShkCount,
950 T_retire,
951 UnempPrb,
952 IncUnemp,
953 UnempPrbRet,
954 IncUnempRet,
955 RNG,
956 )
957 PermShkDstnBase = get_PermShkDstn_from_IncShkDstn(IncShkDstnBase, RNG)
958 TranShkDstnBase = get_TranShkDstn_from_IncShkDstn(IncShkDstnBase, RNG)
960 equiprobable_one_seventh = np.ones(7) / 7.0
961 if CollegeBool:
962 # Copy Mateo's college-educated expense shock distribution
963 exp_shks_51_to_55 = np.array([0.000, 0.003, 0.007, 0.012, 0.021, 0.039, 0.121])
964 exp_shks_56_to_60 = np.array([0.001, 0.005, 0.010, 0.016, 0.027, 0.049, 0.163])
965 exp_shks_61_to_65 = np.array([0.002, 0.007, 0.014, 0.023, 0.040, 0.078, 0.227])
966 exp_shks_66_to_70 = np.array([0.003, 0.010, 0.019, 0.031, 0.050, 0.089, 0.227])
967 exp_shks_71_to_75 = np.array([0.004, 0.013, 0.024, 0.039, 0.060, 0.103, 0.262])
968 exp_shks_76_to_80 = np.array([0.004, 0.015, 0.028, 0.047, 0.074, 0.123, 0.294])
969 exp_shks_81_to_85 = np.array([0.004, 0.017, 0.033, 0.054, 0.089, 0.155, 0.410])
970 exp_shks_86_to_90 = np.array([0.002, 0.015, 0.033, 0.057, 0.100, 0.191, 0.719])
971 exp_shks_91_plus = np.array([0.000, 0.015, 0.039, 0.075, 0.160, 0.389, 1.485])
972 exp_shks_all = (
973 50 * [None]
974 + 5 * [exp_shks_51_to_55]
975 + 5 * [exp_shks_56_to_60]
976 + 5 * [exp_shks_61_to_65]
977 + 5 * [exp_shks_66_to_70]
978 + 5 * [exp_shks_71_to_75]
979 + 5 * [exp_shks_76_to_80]
980 + 5 * [exp_shks_81_to_85]
981 + 5 * [exp_shks_86_to_90]
982 + 31 * [exp_shks_91_plus]
983 )
985 else:
986 # Copy Mateo's high school-educated expense shock distribution
987 exp_shks_51_to_55 = np.array([0.000, 0.004, 0.010, 0.019, 0.033, 0.064, 0.205])
988 exp_shks_56_to_60 = np.array([0.000, 0.005, 0.013, 0.023, 0.040, 0.077, 0.245])
989 exp_shks_61_to_65 = np.array([0.000, 0.008, 0.018, 0.032, 0.055, 0.104, 0.290])
990 exp_shks_66_to_70 = np.array([0.001, 0.011, 0.023, 0.038, 0.064, 0.111, 0.264])
991 exp_shks_71_to_75 = np.array([0.002, 0.014, 0.028, 0.046, 0.074, 0.126, 0.293])
992 exp_shks_76_to_80 = np.array([0.001, 0.015, 0.031, 0.053, 0.084, 0.143, 0.346])
993 exp_shks_81_to_85 = np.array([0.001, 0.016, 0.033, 0.059, 0.096, 0.168, 0.433])
994 exp_shks_86_to_90 = np.array([0.000, 0.016, 0.036, 0.066, 0.110, 0.229, 0.849])
995 exp_shks_91_plus = np.array([0.000, 0.011, 0.034, 0.069, 0.131, 0.301, 1.479])
996 exp_shks_all = (
997 50 * [None]
998 + 5 * [exp_shks_51_to_55]
999 + 5 * [exp_shks_56_to_60]
1000 + 5 * [exp_shks_61_to_65]
1001 + 5 * [exp_shks_66_to_70]
1002 + 5 * [exp_shks_71_to_75]
1003 + 5 * [exp_shks_76_to_80]
1004 + 5 * [exp_shks_81_to_85]
1005 + 5 * [exp_shks_86_to_90]
1006 + 31 * [exp_shks_91_plus]
1007 )
1009 # Incorporate the expense shock distribution into the transitory income shock distribution
1010 IncShkDstn = []
1011 for t in range(T_cycle):
1012 age = age_min + t
1013 if age >= len(exp_shks_all):
1014 raise ValueError(
1015 f"Age {age} is outside the supported calibration range for expense shocks "
1016 f"(maximum supported index is {len(exp_shks_all) - 1})."
1017 )
1018 exp_shks_t = exp_shks_all[age]
1020 # If age <= 50, just use the baseline income shock distribution
1021 if exp_shks_t is None:
1022 IncShkDstn.append(IncShkDstnBase[t])
1023 continue
1025 # Otherwise, make a distribution of net transitory income as the difference
1026 # between the transitory income shock and the medical expense shock
1027 seed_t = RNG.integers(0, 2**31 - 1)
1028 PermShkDstn_t = PermShkDstnBase[t]
1029 ExpShkDstn_t = DiscreteDistribution(
1030 pmv=equiprobable_one_seventh,
1031 atoms=exp_shks_t,
1032 )
1034 # Prepare these for multiplication with broadcasting
1035 TranAtomsBase = TranShkDstnBase[t].atoms[0][:, np.newaxis]
1036 TranProbsBase = TranShkDstnBase[t].pmv[:, np.newaxis]
1037 ExpAtoms = ExpShkDstn_t.atoms[0][np.newaxis, :]
1038 ExpProbs = ExpShkDstn_t.pmv[np.newaxis, :]
1040 # Calculate net income and joint probability
1041 TranAtoms = np.maximum(TranAtomsBase - ExpAtoms, 0.0).flatten()
1042 TranProbs = (TranProbsBase * ExpProbs).flatten()
1044 # Make this period's joint income distribution
1045 TranShkDstn_t = DiscreteDistribution(pmv=TranProbs, atoms=TranAtoms)
1046 IncShkDstn_t = DiscreteDistributionLabeled.from_unlabeled(
1047 combine_indep_dstns(PermShkDstn_t, TranShkDstn_t, seed=seed_t),
1048 name="Income shocks with medical expenses",
1049 var_names=["PermShk", "TranShk"],
1050 )
1051 IncShkDstn.append(IncShkDstn_t)
1053 return IncShkDstn
1056###############################################################################
1059def get_PermShkDstn_from_IncShkDstn(IncShkDstn, RNG):
1060 PermShkDstn = [
1061 this.make_univariate(0, seed=RNG.integers(0, 2**31 - 1)) for this in IncShkDstn
1062 ]
1063 return IndexDistribution(distributions=PermShkDstn, seed=RNG.integers(0, 2**31 - 1))
1066def get_TranShkDstn_from_IncShkDstn(IncShkDstn, RNG):
1067 TranShkDstn = [
1068 this.make_univariate(1, seed=RNG.integers(0, 2**31 - 1)) for this in IncShkDstn
1069 ]
1070 return IndexDistribution(distributions=TranShkDstn, seed=RNG.integers(0, 2**31 - 1))
1073def get_PermShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
1074 PermShkDstn = [
1075 [
1076 this.make_univariate(0, seed=RNG.integers(0, 2**31 - 1))
1077 for this in IncShkDstn_t
1078 ]
1079 for IncShkDstn_t in IncShkDstn
1080 ]
1081 return PermShkDstn
1084def get_TranShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
1085 TranShkDstn = [
1086 [
1087 this.make_univariate(1, seed=RNG.integers(0, 2**31 - 1))
1088 for this in IncShkDstn_t
1089 ]
1090 for IncShkDstn_t in IncShkDstn
1091 ]
1092 return TranShkDstn
1095def get_TranShkGrid_from_TranShkDstn(T_cycle, TranShkDstn):
1096 TranShkGrid = [TranShkDstn[t].atoms.flatten() for t in range(T_cycle)]
1097 return TranShkGrid
1100def make_polynomial_PermGroFac(T_cycle, PermGroFac_coeffs, age_0=0.0, age_step=1.0):
1101 """
1102 Construct the profile of permanent growth factors by age using polynomial coefficients.
1104 Parameters
1105 ----------
1106 T_cycle : int
1107 Number of non-terminal period's in this agent's cycle.
1108 PermGroFac_coeffs : [float]
1109 Arbitrary length list or 1D vector of polynomial coefficients of age on
1110 permanent income growth factor.
1111 age_0 : float, optional
1112 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
1113 The default is 0.
1114 age_step : float, optional
1115 Age increment in the model, with respect to the polynomial coefficients.
1116 The default is 1.
1118 Returns
1119 -------
1120 PermGroFac : [float]
1121 List of permanent income growth factors, one per period.
1122 """
1123 PermGroFac = make_polynomial_params(
1124 PermGroFac_coeffs, T_cycle, offset=0.0, step=1.0
1125 )
1126 return PermGroFac.tolist()
1129def make_polynomial_PermShkStd(T_cycle, PermShkStd_coeffs, age_0=0.0, age_step=1.0):
1130 """
1131 Construct the profile of (log) permanent income shock standard deviations by
1132 age using polynomial coefficients.
1134 Parameters
1135 ----------
1136 T_cycle : int
1137 Number of non-terminal period's in this agent's cycle.
1138 PermGroFac_coeffs : [float]
1139 Arbitrary length list or 1D vector of polynomial coefficients of age on
1140 (log) permanent income shock standard deviation.
1141 age_0 : float, optional
1142 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
1143 The default is 0.
1144 age_step : float, optional
1145 Age increment in the model, with respect to the polynomial coefficients.
1146 The default is 1.
1148 Returns
1149 -------
1150 PermShkStd : [float]
1151 List of (log) permanent income shock standard deviations, one per period.
1152 """
1153 PermShkStd = make_polynomial_params(
1154 PermShkStd_coeffs, T_cycle, offset=0.0, step=1.0
1155 )
1156 return PermShkStd.tolist()
1159def make_polynomial_TranShkStd(T_cycle, TranShkStd_coeffs, age_0=0.0, age_step=1.0):
1160 """
1161 Construct the profile of (log) transitory income shock standard deviations by
1162 age using polynomial coefficients.
1164 Parameters
1165 ----------
1166 T_cycle : int
1167 Number of non-terminal period's in this agent's cycle.
1168 PermGroFac_coeffs : [float]
1169 Arbitrary length list or 1D vector of polynomial coefficients of age on
1170 (log) transitory income shock standard deviation.
1171 age_0 : float, optional
1172 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
1173 The default is 0.
1174 age_step : float, optional
1175 Age increment in the model, with respect to the polynomial coefficients.
1176 The default is 1.
1178 Returns
1179 -------
1180 TranShkStd : [float]
1181 List of (log) permanent income shock standard deviations, one per period.
1182 """
1183 TranShkStd = make_polynomial_params(
1184 TranShkStd_coeffs, T_cycle, offset=0.0, step=1.0
1185 )
1186 return TranShkStd.tolist()
1189class pLvlFuncAR1(MetricObject):
1190 """
1191 A class for representing AR1-style persistent income growth functions.
1193 Parameters
1194 ----------
1195 pLogMean : float
1196 Log persistent income level toward which we are drawn.
1197 PermGroFac : float
1198 Autonomous (e.g. life cycle) pLvl growth (does not AR1 decay).
1199 Corr : float
1200 Correlation coefficient on log income.
1201 """
1203 def __init__(self, pLogMean, PermGroFac, Corr):
1204 self.pLogMean = pLogMean
1205 self.LogGroFac = np.log(PermGroFac)
1206 self.Corr = Corr
1208 def __call__(self, pLvlNow):
1209 """
1210 Returns expected persistent income level next period as a function of
1211 this period's persistent income level.
1213 Parameters
1214 ----------
1215 pLvlNow : np.array
1216 Array of current persistent income levels.
1218 Returns
1219 -------
1220 pLvlNext : np.array
1221 Identically shaped array of next period persistent income levels.
1222 """
1223 pLvlNext = np.exp(
1224 self.Corr * np.log(pLvlNow)
1225 + (1.0 - self.Corr) * self.pLogMean
1226 + self.LogGroFac
1227 )
1228 return pLvlNext
1231def make_PermGroFac_from_ind_and_agg(PermGroFacInd, PermGroFacAgg):
1232 """
1233 A very simple function that constructs *overall* permanent income growth over
1234 the lifecycle as the sum of idiosyncratic productivity growth PermGroFacInd and
1235 aggregate productivity growth PermGroFacAgg. In most HARK models, PermGroFac
1236 is treated as the overall permanent income growth factor, regardless of source.
1237 In applications in which the user has estimated permanent income growth from
1238 *cross sectional* data, or wants to investigate how a change in aggregate
1239 productivity growth affects behavior or outcomes, this function can be used
1240 as the constructor for PermGroFac.
1242 To use this function, specify idiosyncratic productivity growth in the attribute
1243 PermGroFacInd (or construct it), and put this function as the entry for PermGroFac
1244 in the constructors dictionary of your AgentType subclass instances.
1246 Parameters
1247 ----------
1248 PermGroFacInd : [float] or np.array
1249 Lifecycle sequence of idiosyncratic permanent productivity growth factors.
1250 These represent individual-based productivity factors, like experience.
1251 PermGroFacAgg : float
1252 Constant aggregate permanent growth factor, representing (e.g.) TFP.
1254 Returns
1255 -------
1256 PermGroFac : [float] or np.array
1257 Lifecycle sequence of overall permanent productivity growth factors.
1258 Returns same type as PermGroFacInd.
1259 """
1260 PermGroFac = [PermGroFacAgg * G for G in PermGroFacInd]
1261 if type(PermGroFacInd) is np.array:
1262 PermGroFac = np.array(PermGroFac)
1263 return PermGroFac
1266###############################################################################
1268# Define income processes that can be used in the ConsGenIncProcess model
1271def make_trivial_pLvlNextFunc(T_cycle):
1272 r"""
1273 A dummy function that creates default trivial permanent income dynamics:
1274 none at all! Simply returns a list of IdentityFunctions, one for each period.
1276 .. math::
1277 G_{t}(x) = x
1279 Parameters
1280 ----------
1281 T_cycle : int
1282 Number of non-terminal periods in the agent's problem.
1284 Returns
1285 -------
1286 pLvlNextFunc : [IdentityFunction]
1287 List of trivial permanent income dynamic functions.
1288 """
1289 pLvlNextFunc_basic = IdentityFunction()
1290 pLvlNextFunc = T_cycle * [pLvlNextFunc_basic]
1291 return pLvlNextFunc
1294def make_explicit_perminc_pLvlNextFunc(T_cycle, PermGroFac):
1295 r"""
1296 A function that creates permanent income dynamics as a sequence of linear
1297 functions, indicating constant expected permanent income growth across
1298 permanent income levels.
1300 .. math::
1301 G_{t+1} (x) = \Gamma_{t+1} x
1303 Parameters
1304 ----------
1305 T_cycle : int
1306 Number of non-terminal periods in the agent's problem.
1307 PermGroFac : [float]
1308 List of permanent income growth factors over the agent's problem.
1310 Returns
1311 -------
1312 pLvlNextFunc : [LinearInterp]
1313 List of linear functions representing constant permanent income growth
1314 rate, regardless of current permanent income level.
1315 """
1316 pLvlNextFunc = []
1317 for t in range(T_cycle):
1318 pLvlNextFunc.append(
1319 LinearInterp(np.array([0.0, 1.0]), np.array([0.0, PermGroFac[t]]))
1320 )
1321 return pLvlNextFunc
1324def make_AR1_style_pLvlNextFunc(T_cycle, pLogInitMean, PermGroFac, PrstIncCorr):
1325 r"""
1326 A function that creates permanent income dynamics as a sequence of AR1-style
1327 functions. If cycles=0, the product of PermGroFac across all periods must be
1328 1.0, otherwise this method is invalid.
1330 .. math::
1331 \begin{align}
1332 log(G_{t+1} (x)) &=\varphi log(x) + (1-\varphi) log(\overline{P}_{t})+log(\Gamma_{t+1}) + log(\psi_{t+1}), \\
1333 \overline{P}_{t+1} &= \overline{P}_{t} \Gamma_{t+1} \\
1334 \end{align}
1336 Parameters
1337 ----------
1338 T_cycle : int
1339 Number of non-terminal periods in the agent's problem.
1340 pLogInitMean : float
1341 Mean of log permanent income at initialization.
1342 PermGroFac : [float]
1343 List of permanent income growth factors over the agent's problem.
1344 PrstIncCorr : float
1345 Correlation coefficient on log permanent income today on log permanent
1346 income in the succeeding period.
1348 Returns
1349 -------
1350 pLvlNextFunc : [pLvlFuncAR1]
1351 List of AR1-style persistent income dynamics functions
1352 """
1353 pLvlNextFunc = []
1354 pLogMean = pLogInitMean # Initial mean (log) persistent income
1355 for t in range(T_cycle):
1356 pLvlNextFunc.append(pLvlFuncAR1(pLogMean, PermGroFac[t], PrstIncCorr))
1357 pLogMean += np.log(PermGroFac[t])
1358 return pLvlNextFunc
1361###############################################################################
1364def make_basic_pLvlPctiles(
1365 pLvlPctiles_count,
1366 pLvlPctiles_bound=[0.001, 0.999],
1367 pLvlPctiles_tail_count=0,
1368 pLvlPctiles_tail_order=np.e,
1369):
1370 """
1371 Make a relatively basic specification for pLvlPctiles by choosing the number
1372 of uniformly spaced nodes in the "body", the percentile boundaries for the
1373 body, the number of nodes in each tail, and the order/factor by which the
1374 tail percentiles approach 0 and 1 respectively.
1376 Parameters
1377 ----------
1378 pLvlPctile_count : int
1379 Number of nodes in the "body" of the percentile set.
1380 pLvlPctile_bound : [float,float], optional
1381 Percentile bounds for the "body" of the set. The default is [0.0, 1.0].
1382 pLvlPctile_tail_count : int, optional
1383 Number of nodes in each extant tail of the set. The default is 0.
1384 pLvlPctile_tail_order : float, optional
1385 Factor by which tail percentiles shrink toward 0 and 1. The default is np.e.
1387 Returns
1388 -------
1389 pLvlPctiles : np.array
1390 Array of percentiles of pLvl, usually used to construct pLvlGrid using
1391 the function below.
1392 """
1393 bound = pLvlPctiles_bound
1394 fac = 1.0 / pLvlPctiles_tail_order
1395 body = np.linspace(bound[0], bound[1], num=pLvlPctiles_count)
1397 if bound[0] > 0.0:
1398 lower = []
1399 val = bound[0]
1400 for i in range(pLvlPctiles_tail_count):
1401 val *= fac
1402 lower.append(val)
1403 lower.reverse()
1404 lower = np.array(lower)
1405 else:
1406 lower = np.array([])
1408 if bound[1] < 1.0:
1409 upper = []
1410 val = 1.0 - bound[1]
1411 for i in range(pLvlPctiles_tail_count):
1412 val *= fac
1413 upper.append(val)
1414 upper = 1.0 - np.array(upper)
1415 else:
1416 upper = np.array([])
1418 pLvlPctiles = np.concatenate((lower, body, upper))
1419 return pLvlPctiles
1422def make_pLvlGrid_by_simulation(
1423 cycles,
1424 T_cycle,
1425 PermShkDstn,
1426 pLvlNextFunc,
1427 LivPrb,
1428 pLogInitMean,
1429 pLogInitStd,
1430 pLvlPctiles,
1431 pLvlExtra=None,
1432):
1433 """
1434 Construct the permanent income grid for each period of the problem by simulation.
1435 If the model is infinite horizon (cycles=0), an approximation of the long run
1436 steady state distribution of permanent income is used (by simulating many periods).
1437 If the model is lifecycle (cycles=1), explicit simulation is used. In either
1438 case, the input pLvlPctiles is used to choose percentiles from the distribution.
1440 If the problem is neither infinite horizon nor life-cycle, this method will fail.
1441 If the problem is infinite horizon, cumprod(PermGroFac) must equal one.
1443 Parameters
1444 ----------
1445 cycles : int
1446 Number of times the sequence of periods happens for the agent; must be 0 or 1.
1447 T_cycle : int
1448 Number of non-terminal periods in the agent's problem.
1449 PermShkDstn : [distribution]
1450 List of permanent shock distributions in each period of the problem.
1451 pLvlNextFunc : [function]
1452 List of permanent income dynamic functions.
1453 LivPrb : [float]
1454 List of survival probabilities by period of the cycle. Only used in infinite
1455 horizon specifications.
1456 pLogInitMean : float
1457 Mean of log permanent income at initialization.
1458 pLogInitStd : float
1459 Standard deviaition of log permanent income at initialization.
1460 pLvlPctiles : [float]
1461 List or array of percentiles (between 0 and 1) of permanent income to
1462 use for the pLvlGrid.
1463 pLvlExtra : None or [float], optional
1464 Additional pLvl values to automatically include in pLvlGrid.
1466 Returns
1467 -------
1468 pLvlGrid : [np.array]
1469 List of permanent income grids for each period, constructed by simulating
1470 the permanent income process and extracting specified percentiles.
1471 """
1472 LivPrbAll = np.array(LivPrb)
1473 Agent_N = 100000
1475 # Simulate the distribution of persistent income levels by t_cycle in a lifecycle model
1476 if cycles == 1:
1477 pLvlNow = Lognormal(pLogInitMean, sigma=pLogInitStd, seed=31382).draw(Agent_N)
1478 pLvlGrid = [] # empty list of time-varying persistent income grids
1479 # Calculate distribution of persistent income in each period of lifecycle
1480 for t in range(T_cycle):
1481 if t > 0:
1482 PermShkNow = PermShkDstn[t - 1].draw(N=Agent_N)
1483 pLvlNow = pLvlNextFunc[t - 1](pLvlNow) * PermShkNow
1484 pLvlGrid.append(get_percentiles(pLvlNow, percentiles=pLvlPctiles))
1486 # Calculate "stationary" distribution in infinite horizon (might vary across periods of cycle)
1487 elif cycles == 0:
1488 T_long = (
1489 1000 # Number of periods to simulate to get to "stationary" distribution
1490 )
1491 pLvlNow = Lognormal(mu=pLogInitMean, sigma=pLogInitStd, seed=31382).draw(
1492 Agent_N
1493 )
1494 t_cycle = np.zeros(Agent_N, dtype=int)
1495 for t in range(T_long):
1496 # Determine who dies and replace them with newborns
1497 LivPrb = LivPrbAll[t_cycle]
1498 draws = Uniform(seed=t).draw(Agent_N)
1499 who_dies = draws > LivPrb
1500 pLvlNow[who_dies] = Lognormal(
1501 pLogInitMean, pLogInitStd, seed=t + 92615
1502 ).draw(np.sum(who_dies))
1503 t_cycle[who_dies] = 0
1505 for j in range(T_cycle): # Update persistent income
1506 these = t_cycle == j
1507 PermShkTemp = PermShkDstn[j].draw(N=np.sum(these))
1508 pLvlNow[these] = pLvlNextFunc[j](pLvlNow[these]) * PermShkTemp
1509 t_cycle = t_cycle + 1
1510 t_cycle[t_cycle == T_cycle] = 0
1512 # We now have a "long run stationary distribution", extract percentiles
1513 pLvlGrid = [] # empty list of time-varying persistent income grids
1514 for t in range(T_cycle):
1515 these = t_cycle == t
1516 pLvlGrid.append(get_percentiles(pLvlNow[these], percentiles=pLvlPctiles))
1518 # Throw an error if cycles>1
1519 else:
1520 assert False, "Can only handle cycles=0 or cycles=1!"
1522 # Insert any additional requested points into the pLvlGrid
1523 if pLvlExtra is not None:
1524 pLvlExtra_alt = np.array(pLvlExtra)
1525 for t in range(T_cycle):
1526 pLvlGrid_t = pLvlGrid[t]
1527 pLvlGrid[t] = np.unique(np.concatenate((pLvlGrid_t, pLvlExtra_alt)))
1529 return pLvlGrid
1532###############################################################################
1535def make_persistent_income_process_dict(
1536 cycles,
1537 T_cycle,
1538 PermShkStd,
1539 PermShkCount,
1540 pLogInitMean,
1541 pLogInitStd,
1542 PermGroFac,
1543 PrstIncCorr,
1544 pLogCount,
1545 pLogRange,
1546):
1547 """
1548 Constructs a dictionary with several elements that characterize the income
1549 process for an agent with AR(1) persistent income process and lognormal transitory
1550 shocks (with unemployment). The produced dictionary includes permanent income
1551 grids and transition matrices and a mean permanent income lifecycle sequence.
1553 This function only works with cycles>0 or T_cycle=1.
1555 Parameters
1556 ----------
1557 cycles : int
1558 Number of times the agent's sequence of periods repeats.
1559 T_cycle : int
1560 Number of periods in the sequence.
1561 PermShkStd : [float]
1562 Standard deviation of mean one permanent income shocks in each period,
1563 assumed to be lognormally distributed.
1564 PermShkCount : int
1565 Number of discrete nodes in the permanent income shock distribution (can
1566 be used during simulation).
1567 pLogInitMean : float
1568 Mean of log permanent income at model entry.
1569 pLogInitStd : float
1570 Standard deviation of log permanent income at model entry.
1571 PermGroFac : [float]
1572 Lifecycle sequence of permanent income growth factors, *not* offset by
1573 one period as in most other HARK models.
1574 PrstIncCorr : float
1575 Correlation coefficient of the persistent component of income.
1576 pLogCount : int
1577 Number of gridpoints in the grid of (log) persistent income deviations.
1578 pLogRange : float
1579 Upper bound of log persistent income grid, in standard deviations from
1580 the mean; grid has symmetric lower bound.
1582 Returns
1583 -------
1584 IncomeProcessDict : dict
1585 Dictionary with the following entries.
1587 pLogGrid : [np.array]
1588 Age-dependent grids of log persistent income, in deviations from mean.
1589 pLvlMean : [float]
1590 Mean persistent income level by age.
1591 pLogMrkvArray : [np.array]
1592 Age-dependent Markov transition arrays among pLog levels at the start of
1593 each period in the sequence.
1594 """
1595 if cycles == 0:
1596 if T_cycle > 1:
1597 raise ValueError(
1598 "Can't handle infinite horizon models with more than one period!"
1599 )
1600 if PermGroFac[0] != 1.0:
1601 raise ValueError(
1602 "Can't handle permanent income growth in infinite horizon!"
1603 )
1605 # The single pLogGrid and transition matrix can be generated by the basic
1606 # Tauchen AR(1) method from HARK.distributions.
1607 pLogGrid, pLogMrkvArray = make_tauchen_ar1(
1608 pLogCount,
1609 sigma=PermShkStd[0],
1610 ar_1=PrstIncCorr,
1611 bound=pLogRange,
1612 )
1613 pLogGrid = [pLogGrid]
1614 pLogMrkvArray = [pLogMrkvArray]
1615 pLvlMean = [np.exp(pLogInitMean + 0.5 * pLogInitStd**2)]
1617 else:
1618 # Start with the pLog distribution at model entry
1619 pLvlMeanNow = np.exp(pLogInitMean + 0.5 * pLogInitStd**2)
1620 pLogStdNow = pLogInitStd
1621 pLogGridPrev = np.linspace(
1622 -pLogRange * pLogStdNow, pLogRange * pLogStdNow, pLogCount
1623 )
1625 # Initialize empty lists to hold output
1626 pLogGrid = []
1627 pLogMrkvArray = []
1628 pLvlMean = []
1630 for c in range(cycles):
1631 for t in range(T_cycle):
1632 # Update the distribution of persistent income deviations from mean
1633 pLvlMeanNow *= PermGroFac[t]
1634 pLogStdNow = np.sqrt(
1635 (PrstIncCorr * pLogStdNow) ** 2 + PermShkStd[t] ** 2
1636 )
1637 pLogGridNow = np.linspace(
1638 -pLogRange * pLogStdNow, pLogRange * pLogStdNow, pLogCount
1639 )
1641 # Compute transition distances from prior grid to this one
1642 pLogCuts = (pLogGridNow[1:] + pLogGridNow[:-1]) / 2.0
1643 pLogCuts = np.concatenate(([-np.inf], pLogCuts, [np.inf]))
1644 distances = np.reshape(pLogCuts, (1, pLogCount + 1)) - np.reshape(
1645 PrstIncCorr * pLogGridPrev, (pLogCount, 1)
1646 )
1647 distances /= PermShkStd
1649 # Compute transition probabilities, ensuring that very small
1650 # probabilities are treated identically in both directions
1651 cdf_array = norm.cdf(distances)
1652 sf_array = norm.sf(distances)
1653 pLogMrkvNow = cdf_array[:, 1:] - cdf_array[:, :-1]
1654 pLogMrkvNowAlt = sf_array[:, :-1] - sf_array[:, 1:]
1655 pLogMrkvNow = np.maximum(pLogMrkvNow, pLogMrkvNowAlt)
1656 pLogMrkvNow /= np.sum(pLogMrkvNow, axis=1, keepdims=True)
1658 # Add this period's output to the lists
1659 pLogGrid.append(pLogGridNow)
1660 pLogMrkvArray.append(pLogMrkvNow)
1661 pLvlMean.append(pLvlMeanNow)
1662 pLogGridPrev = pLogGridNow
1664 # Gather and return the output
1665 IncomeProcessDict = {
1666 "pLogGrid": pLogGrid,
1667 "pLogMrkvArray": pLogMrkvArray,
1668 "pLvlMean": pLvlMean,
1669 }
1670 return IncomeProcessDict
1673###############################################################################
1676def combine_ind_and_agg_income_shocks(IncShkDstnInd, AggShkDstn, RNG, T_cycle):
1677 """
1678 Updates attribute IncShkDstn by combining idiosyncratic shocks with aggregate shocks.
1680 Parameters
1681 ----------
1682 IncShkDstnInd : [DiscreteDistribution]
1683 Age-varying list of idiosyncratic income shock distributions.
1684 AggShkDstn : DiscreteDistribution
1685 Aggregate productivity shock distribution.
1686 RNG : RandomState
1687 Internal random number generator.
1688 T_cycle : int
1689 Number of periods in the agent's sequence of problems.
1691 Returns
1692 -------
1693 IncShkDstn : [DiscreteDistribution]
1694 Combined idiosyncratic and aggregate income shocks, discretized, by age.
1695 """
1696 IncShkDstn = [
1697 combine_indep_dstns(
1698 IncShkDstnInd[t], AggShkDstn, seed=RNG.integers(0, 2**31 - 1)
1699 )
1700 for t in range(T_cycle)
1701 ]
1702 return IncShkDstn
1705def combine_markov_ind_and_agg_income_shocks(
1706 IncShkDstnInd, AggShkDstn, MrkvArray, RNG, T_cycle
1707):
1708 """
1709 Updates attribute IncShkDstn by combining idiosyncratic shocks with aggregate shocks,
1710 for a model with an aggregate Markov discrete state.
1712 Parameters
1713 ----------
1714 IncShkDstnInd : [DiscreteDistribution]
1715 Age-varying list of idiosyncratic income shock distributions.
1716 AggShkDstn : DiscreteDistribution
1717 Aggregate productivity shock distribution.
1718 MrkvArray : np.array
1719 Markov transition matrix among discrete macroeconomic states.
1720 RNG : RandomState
1721 Internal random number generator.
1722 T_cycle : int
1723 Number of periods in the agent's sequence of problems.
1725 Returns
1726 -------
1727 IncShkDstn : [[DiscreteDistribution]]
1728 Combined idiosyncratic and aggregate income shocks, discretized, by age
1729 and Markov state.
1730 """
1731 IncShkDstn = []
1732 N = MrkvArray.shape[0]
1733 for t in range(T_cycle):
1734 IncShkDstn.append(
1735 [
1736 combine_indep_dstns(
1737 IncShkDstnInd[t][n],
1738 AggShkDstn[n],
1739 seed=RNG.integers(0, 2**31 - 1),
1740 )
1741 for n in range(N)
1742 ]
1743 )
1744 return IncShkDstn