Coverage for HARK/Calibration/Income/IncomeProcesses.py: 86%
319 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-02 05:14 +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 == TranShkStd_K
643 assert K == UnempPrb_K
644 assert K == IncUnemp_K
645 except:
646 raise Exception(
647 "The last dimension of PermShkStd, TranShkStd, IncUnemp,"
648 + " and UnempPrb must all be K, the number of discrete states."
649 )
650 try:
651 if T_retire > 0:
652 assert K == UnempPrbRet.size
653 assert K == IncUnempRet.size
654 except:
655 raise Exception(
656 "When T_retire is not zero, UnempPrbRet and IncUnempRet"
657 + " must be specified as arrays of size K, the number of "
658 + "discrete Markov states."
659 )
660 try:
661 D = UnempPrb.ndim
662 assert D == IncUnemp.ndim
663 if T_retire > 0:
664 assert D == UnempPrbRet.ndim
665 assert D == IncUnempRet.ndim
666 except:
667 raise Exception(
668 "If any of UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet "
669 + "are 2D arrays, then they must *all* be 2D arrays."
670 )
671 try:
672 assert D == 1 or D == 2
673 except:
674 raise Exception(
675 "UnempPrb, IncUnemp, or UnempPrbRet, or IncUnempRet must "
676 + "all be 1D or 2D arrays."
677 )
679 # Prepare lists that don't vary by Markov state
680 PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
681 TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
682 neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
684 # Loop through the Markov states, constructing the lifecycle income process for each one
685 IncShkDstn_by_Mrkv = []
686 for k in range(K):
687 if D == 1: # Unemployment parameters don't vary by age other than retirement
688 if T_retire > 0:
689 UnempPrb_list = [UnempPrb[k]] * normal_length + [
690 UnempPrbRet[k]
691 ] * retire_length
692 IncUnemp_list = [IncUnemp[k]] * normal_length + [
693 IncUnempRet[k]
694 ] * retire_length
695 else:
696 UnempPrb_list = [UnempPrb[k]] * normal_length
697 IncUnemp_list = [IncUnemp[k]] * normal_length
698 else: # Unemployment parameters vary by age
699 UnempPrb_list = UnempPrb[:, k].tolist()
700 IncUnemp_list = IncUnemp[:, k].tolist()
702 PermShkStd_k = PermShkStd[:, k].tolist()
703 TranShkStd_k = TranShkStd[:, k].tolist()
705 IncShkDstn_k = IndexDistribution(
706 engine=BufferStockIncShkDstn,
707 conditional={
708 "sigma_Perm": PermShkStd_k,
709 "sigma_Tran": TranShkStd_k,
710 "n_approx_Perm": PermShkCount_list,
711 "n_approx_Tran": TranShkCount_list,
712 "neutral_measure": neutral_measure_list,
713 "UnempPrb": UnempPrb_list,
714 "IncUnemp": IncUnemp_list,
715 },
716 RNG=RNG,
717 seed=RNG.integers(0, 2**31 - 1),
718 )
719 IncShkDstn_by_Mrkv.append(IncShkDstn_k)
721 # Rearrange the list that was just constructed, so that element [t][k] represents
722 # the income shock distribution in period t, Markov state k.
723 IncShkDstn = []
724 for t in range(T_cycle):
725 IncShkDstn_t = [IncShkDstn_by_Mrkv[k][t] for k in range(K)]
726 IncShkDstn.append(IncShkDstn_t)
728 return IncShkDstn
731def construct_HANK_lognormal_income_process_unemployment(
732 T_cycle,
733 PermShkStd,
734 PermShkCount,
735 TranShkStd,
736 TranShkCount,
737 T_retire,
738 UnempPrb,
739 IncUnemp,
740 UnempPrbRet,
741 IncUnempRet,
742 tax_rate,
743 labor,
744 wage,
745 RNG,
746 neutral_measure=False,
747):
748 """
749 Generates a list of discrete approximations to the income process for each
750 life period, from end of life to beginning of life. Permanent shocks are mean
751 one lognormally distributed with standard deviation PermShkStd[t] during the
752 working life, and degenerate at 1 in the retirement period. Transitory shocks
753 are mean one lognormally distributed with a point mass at IncUnemp with
754 probability UnempPrb while working; they are mean one with a point mass at
755 IncUnempRet with probability UnempPrbRet. Retirement occurs after t=T_retire
756 periods of working.
758 This version of the function incorporates additional flexibility with respect
759 to transitory income (continuous labor supply, wage rate, tax rate) and thus
760 is useful in HANK models (hence the name!).
762 Note 1: All time in this function runs forward, from t=0 to t=T
764 Note 2: All parameters are passed as attributes of the input parameters.
766 Parameters (passed as attributes of the input parameters)
767 ---------------------------------------------------------
768 PermShkStd : [float]
769 List of standard deviations in log permanent income uncertainty during
770 the agent's life.
771 PermShkCount : int
772 The number of approximation points to be used in the discrete approxima-
773 tion to the permanent income shock distribution.
774 TranShkStd : [float]
775 List of standard deviations in log transitory income uncertainty during
776 the agent's life.
777 TranShkCount : int
778 The number of approximation points to be used in the discrete approxima-
779 tion to the permanent income shock distribution.
780 UnempPrb : float or [float]
781 The probability of becoming unemployed during the working period.
782 UnempPrbRet : float or None
783 The probability of not receiving typical retirement income when retired.
784 T_retire : int
785 The index value for the final working period in the agent's life.
786 If T_retire <= 0 then there is no retirement.
787 IncUnemp : float or [float]
788 Transitory income received when unemployed.
789 IncUnempRet : float or None
790 Transitory income received while "unemployed" when retired.
791 tax_rate : [float]
792 List of flat tax rates on labor income, age-varying.
793 labor : [float]
794 List of intensive margin labor supply, age-varying.
795 wage : [float]
796 List of wage rate scaling factors, age-varying.
797 T_cycle : int
798 Total number of non-terminal periods in the consumer's sequence of periods.
800 Returns
801 -------
802 IncShkDstn : [distribution.Distribution]
803 A list with T_cycle elements, each of which is a
804 discrete approximation to the income process in a period.
805 PermShkDstn : [[distribution.Distributiony]]
806 A list with T_cycle elements, each of which
807 a discrete approximation to the permanent income shocks.
808 TranShkDstn : [[distribution.Distribution]]
809 A list with T_cycle elements, each of which
810 a discrete approximation to the transitory income shocks.
811 """
812 if T_retire > 0:
813 normal_length = T_retire
814 retire_length = T_cycle - T_retire
815 else:
816 normal_length = T_cycle
817 retire_length = 0
819 if all(
820 [
821 isinstance(x, (float, int)) or (x is None)
822 for x in [UnempPrb, IncUnemp, UnempPrbRet, IncUnempRet]
823 ]
824 ):
825 UnempPrb_list = [UnempPrb] * normal_length + [UnempPrbRet] * retire_length
826 IncUnemp_list = [IncUnemp] * normal_length + [IncUnempRet] * retire_length
828 elif all([isinstance(x, list) for x in [UnempPrb, IncUnemp]]):
829 UnempPrb_list = UnempPrb
830 IncUnemp_list = IncUnemp
832 else:
833 raise Exception(
834 "Unemployment must be specified either using floats for UnempPrb,"
835 + "IncUnemp, UnempPrbRet, and IncUnempRet, in which case the "
836 + "unemployment probability and income change only with retirement, or "
837 + "using lists of length T_cycle for UnempPrb and IncUnemp, specifying "
838 + "each feature at every age."
839 )
841 PermShkCount_list = [PermShkCount] * normal_length + [1] * retire_length
842 TranShkCount_list = [TranShkCount] * normal_length + [1] * retire_length
843 neutral_measure_list = [neutral_measure] * len(PermShkCount_list)
845 IncShkDstn = IndexDistribution(
846 engine=IncShkDstn_HANK,
847 conditional={
848 "sigma_Perm": PermShkStd,
849 "sigma_Tran": TranShkStd,
850 "n_approx_Perm": PermShkCount_list,
851 "n_approx_Tran": TranShkCount_list,
852 "neutral_measure": neutral_measure_list,
853 "UnempPrb": UnempPrb_list,
854 "IncUnemp": IncUnemp_list,
855 "wage": wage,
856 "tax_rate": tax_rate,
857 "labor": labor,
858 },
859 RNG=RNG,
860 )
862 return IncShkDstn
865###############################################################################
868def get_PermShkDstn_from_IncShkDstn(IncShkDstn, RNG):
869 PermShkDstn = [
870 this.make_univariate(0, seed=RNG.integers(0, 2**31 - 1)) for this in IncShkDstn
871 ]
872 return IndexDistribution(distributions=PermShkDstn, seed=RNG.integers(0, 2**31 - 1))
875def get_TranShkDstn_from_IncShkDstn(IncShkDstn, RNG):
876 TranShkDstn = [
877 this.make_univariate(1, seed=RNG.integers(0, 2**31 - 1)) for this in IncShkDstn
878 ]
879 return IndexDistribution(distributions=TranShkDstn, seed=RNG.integers(0, 2**31 - 1))
882def get_PermShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
883 PermShkDstn = [
884 [
885 this.make_univariate(0, seed=RNG.integers(0, 2**31 - 1))
886 for this in IncShkDstn_t
887 ]
888 for IncShkDstn_t in IncShkDstn
889 ]
890 return PermShkDstn
893def get_TranShkDstn_from_IncShkDstn_markov(IncShkDstn, RNG):
894 TranShkDstn = [
895 [
896 this.make_univariate(1, seed=RNG.integers(0, 2**31 - 1))
897 for this in IncShkDstn_t
898 ]
899 for IncShkDstn_t in IncShkDstn
900 ]
901 return TranShkDstn
904def get_TranShkGrid_from_TranShkDstn(T_cycle, TranShkDstn):
905 TranShkGrid = [TranShkDstn[t].atoms.flatten() for t in range(T_cycle)]
906 return TranShkGrid
909def make_polynomial_PermGroFac(T_cycle, PermGroFac_coeffs, age_0=0.0, age_step=1.0):
910 """
911 Construct the profile of permanent growth factors by age using polynomial coefficients.
913 Parameters
914 ----------
915 T_cycle : int
916 Number of non-terminal period's in this agent's cycle.
917 PermGroFac_coeffs : [float]
918 Arbitrary length list or 1D vector of polynomial coefficients of age on
919 permanent income growth factor.
920 age_0 : float, optional
921 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
922 The default is 0.
923 age_step : float, optional
924 Age increment in the model, with respect to the polynomial coefficients.
925 The default is 1.
927 Returns
928 -------
929 PermGroFac : [float]
930 List of permanent income growth factors, one per period.
931 """
932 PermGroFac = make_polynomial_params(
933 PermGroFac_coeffs, T_cycle, offset=0.0, step=1.0
934 )
935 return PermGroFac.tolist()
938def make_polynomial_PermShkStd(T_cycle, PermShkStd_coeffs, age_0=0.0, age_step=1.0):
939 """
940 Construct the profile of (log) permanent income shock standard deviations by
941 age using polynomial coefficients.
943 Parameters
944 ----------
945 T_cycle : int
946 Number of non-terminal period's in this agent's cycle.
947 PermGroFac_coeffs : [float]
948 Arbitrary length list or 1D vector of polynomial coefficients of age on
949 (log) permanent income shock standard deviation.
950 age_0 : float, optional
951 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
952 The default is 0.
953 age_step : float, optional
954 Age increment in the model, with respect to the polynomial coefficients.
955 The default is 1.
957 Returns
958 -------
959 PermShkStd : [float]
960 List of (log) permanent income shock standard deviations, one per period.
961 """
962 PermShkStd = make_polynomial_params(
963 PermShkStd_coeffs, T_cycle, offset=0.0, step=1.0
964 )
965 return PermShkStd.tolist()
968def make_polynomial_TranShkStd(T_cycle, TranShkStd_coeffs, age_0=0.0, age_step=1.0):
969 """
970 Construct the profile of (log) transitory income shock standard deviations by
971 age using polynomial coefficients.
973 Parameters
974 ----------
975 T_cycle : int
976 Number of non-terminal period's in this agent's cycle.
977 PermGroFac_coeffs : [float]
978 Arbitrary length list or 1D vector of polynomial coefficients of age on
979 (log) transitory income shock standard deviation.
980 age_0 : float, optional
981 Initial age of agents (when t_age=0), with respect to the polynomial coefficients.
982 The default is 0.
983 age_step : float, optional
984 Age increment in the model, with respect to the polynomial coefficients.
985 The default is 1.
987 Returns
988 -------
989 TranShkStd : [float]
990 List of (log) permanent income shock standard deviations, one per period.
991 """
992 TranShkStd = make_polynomial_params(
993 TranShkStd_coeffs, T_cycle, offset=0.0, step=1.0
994 )
995 return TranShkStd.tolist()
998class pLvlFuncAR1(MetricObject):
999 """
1000 A class for representing AR1-style persistent income growth functions.
1002 Parameters
1003 ----------
1004 pLogMean : float
1005 Log persistent income level toward which we are drawn.
1006 PermGroFac : float
1007 Autonomous (e.g. life cycle) pLvl growth (does not AR1 decay).
1008 Corr : float
1009 Correlation coefficient on log income.
1010 """
1012 def __init__(self, pLogMean, PermGroFac, Corr):
1013 self.pLogMean = pLogMean
1014 self.LogGroFac = np.log(PermGroFac)
1015 self.Corr = Corr
1017 def __call__(self, pLvlNow):
1018 """
1019 Returns expected persistent income level next period as a function of
1020 this period's persistent income level.
1022 Parameters
1023 ----------
1024 pLvlNow : np.array
1025 Array of current persistent income levels.
1027 Returns
1028 -------
1029 pLvlNext : np.array
1030 Identically shaped array of next period persistent income levels.
1031 """
1032 pLvlNext = np.exp(
1033 self.Corr * np.log(pLvlNow)
1034 + (1.0 - self.Corr) * self.pLogMean
1035 + self.LogGroFac
1036 )
1037 return pLvlNext
1040def make_PermGroFac_from_ind_and_agg(PermGroFacInd, PermGroFacAgg):
1041 """
1042 A very simple function that constructs *overall* permanent income growth over
1043 the lifecycle as the sum of idiosyncratic productivity growth PermGroFacInd and
1044 aggregate productivity growth PermGroFacAgg. In most HARK models, PermGroFac
1045 is treated as the overall permanent income growth factor, regardless of source.
1046 In applications in which the user has estimated permanent income growth from
1047 *cross sectional* data, or wants to investigate how a change in aggregate
1048 productivity growth affects behavior or outcomes, this function can be used
1049 as the constructor for PermGroFac.
1051 To use this function, specify idiosyncratic productivity growth in the attribute
1052 PermGroFacInd (or construct it), and put this function as the entry for PermGroFac
1053 in the constructors dictionary of your AgentType subclass instances.
1055 Parameters
1056 ----------
1057 PermGroFacInd : [float] or np.array
1058 Lifecycle sequence of idiosyncratic permanent productivity growth factors.
1059 These represent individual-based productivity factors, like experience.
1060 PermGroFacAgg : float
1061 Constant aggregate permanent growth factor, representing (e.g.) TFP.
1063 Returns
1064 -------
1065 PermGroFac : [float] or np.array
1066 Lifecycle sequence of overall permanent productivity growth factors.
1067 Returns same type as PermGroFacInd.
1068 """
1069 PermGroFac = [PermGroFacAgg * G for G in PermGroFacInd]
1070 if type(PermGroFacInd) is np.array:
1071 PermGroFac = np.array(PermGroFac)
1072 return PermGroFac
1075###############################################################################
1077# Define income processes that can be used in the ConsGenIncProcess model
1080def make_trivial_pLvlNextFunc(T_cycle):
1081 r"""
1082 A dummy function that creates default trivial permanent income dynamics:
1083 none at all! Simply returns a list of IdentityFunctions, one for each period.
1085 .. math::
1086 G_{t}(x) = x
1088 Parameters
1089 ----------
1090 T_cycle : int
1091 Number of non-terminal periods in the agent's problem.
1093 Returns
1094 -------
1095 pLvlNextFunc : [IdentityFunction]
1096 List of trivial permanent income dynamic functions.
1097 """
1098 pLvlNextFunc_basic = IdentityFunction()
1099 pLvlNextFunc = T_cycle * [pLvlNextFunc_basic]
1100 return pLvlNextFunc
1103def make_explicit_perminc_pLvlNextFunc(T_cycle, PermGroFac):
1104 r"""
1105 A function that creates permanent income dynamics as a sequence of linear
1106 functions, indicating constant expected permanent income growth across
1107 permanent income levels.
1109 .. math::
1110 G_{t+1} (x) = \Gamma_{t+1} x
1112 Parameters
1113 ----------
1114 T_cycle : int
1115 Number of non-terminal periods in the agent's problem.
1116 PermGroFac : [float]
1117 List of permanent income growth factors over the agent's problem.
1119 Returns
1120 -------
1121 pLvlNextFunc : [LinearInterp]
1122 List of linear functions representing constant permanent income growth
1123 rate, regardless of current permanent income level.
1124 """
1125 pLvlNextFunc = []
1126 for t in range(T_cycle):
1127 pLvlNextFunc.append(
1128 LinearInterp(np.array([0.0, 1.0]), np.array([0.0, PermGroFac[t]]))
1129 )
1130 return pLvlNextFunc
1133def make_AR1_style_pLvlNextFunc(T_cycle, pLogInitMean, PermGroFac, PrstIncCorr):
1134 r"""
1135 A function that creates permanent income dynamics as a sequence of AR1-style
1136 functions. If cycles=0, the product of PermGroFac across all periods must be
1137 1.0, otherwise this method is invalid.
1139 .. math::
1140 \begin{align}
1141 log(G_{t+1} (x)) &=\varphi log(x) + (1-\varphi) log(\overline{P}_{t})+log(\Gamma_{t+1}) + log(\psi_{t+1}), \\
1142 \overline{P}_{t+1} &= \overline{P}_{t} \Gamma_{t+1} \\
1143 \end{align}
1145 Parameters
1146 ----------
1147 T_cycle : int
1148 Number of non-terminal periods in the agent's problem.
1149 pLogInitMean : float
1150 Mean of log permanent income at initialization.
1151 PermGroFac : [float]
1152 List of permanent income growth factors over the agent's problem.
1153 PrstIncCorr : float
1154 Correlation coefficient on log permanent income today on log permanent
1155 income in the succeeding period.
1157 Returns
1158 -------
1159 pLvlNextFunc : [pLvlFuncAR1]
1160 List of AR1-style persistent income dynamics functions
1161 """
1162 pLvlNextFunc = []
1163 pLogMean = pLogInitMean # Initial mean (log) persistent income
1164 for t in range(T_cycle):
1165 pLvlNextFunc.append(pLvlFuncAR1(pLogMean, PermGroFac[t], PrstIncCorr))
1166 pLogMean += np.log(PermGroFac[t])
1167 return pLvlNextFunc
1170###############################################################################
1173def make_basic_pLvlPctiles(
1174 pLvlPctiles_count,
1175 pLvlPctiles_bound=[0.001, 0.999],
1176 pLvlPctiles_tail_count=0,
1177 pLvlPctiles_tail_order=np.e,
1178):
1179 """
1180 Make a relatively basic specification for pLvlPctiles by choosing the number
1181 of uniformly spaced nodes in the "body", the percentile boundaries for the
1182 body, the number of nodes in each tail, and the order/factor by which the
1183 tail percentiles approach 0 and 1 respectively.
1185 Parameters
1186 ----------
1187 pLvlPctile_count : int
1188 Number of nodes in the "body" of the percentile set.
1189 pLvlPctile_bound : [float,float], optional
1190 Percentile bounds for the "body" of the set. The default is [0.0, 1.0].
1191 pLvlPctile_tail_count : int, optional
1192 Number of nodes in each extant tail of the set. The default is 0.
1193 pLvlPctile_tail_order : float, optional
1194 Factor by which tail percentiles shrink toward 0 and 1. The default is np.e.
1196 Returns
1197 -------
1198 pLvlPctiles : np.array
1199 Array of percentiles of pLvl, usually used to construct pLvlGrid using
1200 the function below.
1201 """
1202 bound = pLvlPctiles_bound
1203 fac = 1.0 / pLvlPctiles_tail_order
1204 body = np.linspace(bound[0], bound[1], num=pLvlPctiles_count)
1206 if bound[0] > 0.0:
1207 lower = []
1208 val = bound[0]
1209 for i in range(pLvlPctiles_tail_count):
1210 val *= fac
1211 lower.append(val)
1212 lower.reverse()
1213 lower = np.array(lower)
1214 else:
1215 lower = np.array([])
1217 if bound[1] < 1.0:
1218 upper = []
1219 val = 1.0 - bound[1]
1220 for i in range(pLvlPctiles_tail_count):
1221 val *= fac
1222 upper.append(val)
1223 upper = 1.0 - np.array(upper)
1224 else:
1225 upper = np.array([])
1227 pLvlPctiles = np.concatenate((lower, body, upper))
1228 return pLvlPctiles
1231def make_pLvlGrid_by_simulation(
1232 cycles,
1233 T_cycle,
1234 PermShkDstn,
1235 pLvlNextFunc,
1236 LivPrb,
1237 pLogInitMean,
1238 pLogInitStd,
1239 pLvlPctiles,
1240 pLvlExtra=None,
1241):
1242 """
1243 Construct the permanent income grid for each period of the problem by simulation.
1244 If the model is infinite horizon (cycles=0), an approximation of the long run
1245 steady state distribution of permanent income is used (by simulating many periods).
1246 If the model is lifecycle (cycles=1), explicit simulation is used. In either
1247 case, the input pLvlPctiles is used to choose percentiles from the distribution.
1249 If the problem is neither infinite horizon nor lifecycle, this method will fail.
1250 If the problem is infinite horizon, cumprod(PermGroFac) must equal one.
1252 Parameters
1253 ----------
1254 cycles : int
1255 Number of times the sequence of periods happens for the agent; must be 0 or 1.
1256 T_cycle : int
1257 Number of non-terminal periods in the agent's problem.
1258 PermShkDstn : [distribution]
1259 List of permanent shock distributions in each period of the problem.
1260 pLvlNextFunc : [function]
1261 List of permanent income dynamic functions.
1262 LivPrb : [float]
1263 List of survival probabilities by period of the cycle. Only used in infinite
1264 horizon specifications.
1265 pLogInitMean : float
1266 Mean of log permanent income at initialization.
1267 pLogInitStd : float
1268 Standard deviaition of log permanent income at initialization.
1269 pLvlPctiles : [float]
1270 List or array of percentiles (between 0 and 1) of permanent income to
1271 use for the pLvlGrid.
1272 pLvlExtra : None or [float], optional
1273 Additional pLvl values to automatically include in pLvlGrid.
1275 Returns
1276 -------
1277 pLvlGrid : [np.array]
1278 List of permanent income grids for each period, constructed by simulating
1279 the permanent income process and extracting specified percentiles.
1280 """
1281 LivPrbAll = np.array(LivPrb)
1282 Agent_N = 100000
1284 # Simulate the distribution of persistent income levels by t_cycle in a lifecycle model
1285 if cycles == 1:
1286 pLvlNow = Lognormal(pLogInitMean, sigma=pLogInitStd, seed=31382).draw(Agent_N)
1287 pLvlGrid = [] # empty list of time-varying persistent income grids
1288 # Calculate distribution of persistent income in each period of lifecycle
1289 for t in range(T_cycle):
1290 if t > 0:
1291 PermShkNow = PermShkDstn[t - 1].draw(N=Agent_N)
1292 pLvlNow = pLvlNextFunc[t - 1](pLvlNow) * PermShkNow
1293 pLvlGrid.append(get_percentiles(pLvlNow, percentiles=pLvlPctiles))
1295 # Calculate "stationary" distribution in infinite horizon (might vary across periods of cycle)
1296 elif cycles == 0:
1297 T_long = (
1298 1000 # Number of periods to simulate to get to "stationary" distribution
1299 )
1300 pLvlNow = Lognormal(mu=pLogInitMean, sigma=pLogInitStd, seed=31382).draw(
1301 Agent_N
1302 )
1303 t_cycle = np.zeros(Agent_N, dtype=int)
1304 for t in range(T_long):
1305 # Determine who dies and replace them with newborns
1306 LivPrb = LivPrbAll[t_cycle]
1307 draws = Uniform(seed=t).draw(Agent_N)
1308 who_dies = draws > LivPrb
1309 pLvlNow[who_dies] = Lognormal(
1310 pLogInitMean, pLogInitStd, seed=t + 92615
1311 ).draw(np.sum(who_dies))
1312 t_cycle[who_dies] = 0
1314 for j in range(T_cycle): # Update persistent income
1315 these = t_cycle == j
1316 PermShkTemp = PermShkDstn[j].draw(N=np.sum(these))
1317 pLvlNow[these] = pLvlNextFunc[j](pLvlNow[these]) * PermShkTemp
1318 t_cycle = t_cycle + 1
1319 t_cycle[t_cycle == T_cycle] = 0
1321 # We now have a "long run stationary distribution", extract percentiles
1322 pLvlGrid = [] # empty list of time-varying persistent income grids
1323 for t in range(T_cycle):
1324 these = t_cycle == t
1325 pLvlGrid.append(get_percentiles(pLvlNow[these], percentiles=pLvlPctiles))
1327 # Throw an error if cycles>1
1328 else:
1329 assert False, "Can only handle cycles=0 or cycles=1!"
1331 # Insert any additional requested points into the pLvlGrid
1332 if pLvlExtra is not None:
1333 pLvlExtra_alt = np.array(pLvlExtra)
1334 for t in range(T_cycle):
1335 pLvlGrid_t = pLvlGrid[t]
1336 pLvlGrid[t] = np.unique(np.concatenate((pLvlGrid_t, pLvlExtra_alt)))
1338 return pLvlGrid
1341###############################################################################
1344def make_persistent_income_process_dict(
1345 cycles,
1346 T_cycle,
1347 PermShkStd,
1348 PermShkCount,
1349 pLogInitMean,
1350 pLogInitStd,
1351 PermGroFac,
1352 PrstIncCorr,
1353 pLogCount,
1354 pLogRange,
1355):
1356 """
1357 Constructs a dictionary with several elements that characterize the income
1358 process for an agent with AR(1) persistent income process and lognormal transitory
1359 shocks (with unemployment). The produced dictionary includes permanent income
1360 grids and transition matrices and a mean permanent income lifecycle sequence.
1362 This function only works with cycles>0 or T_cycle=1.
1364 Parameters
1365 ----------
1366 cycles : int
1367 Number of times the agent's sequence of periods repeats.
1368 T_cycle : int
1369 Number of periods in the sequence.
1370 PermShkStd : [float]
1371 Standard deviation of mean one permanent income shocks in each period,
1372 assumed to be lognormally distributed.
1373 PermShkCount : int
1374 Number of discrete nodes in the permanent income shock distribution (can
1375 be used during simulation).
1376 pLogInitMean : float
1377 Mean of log permanent income at model entry.
1378 pLogInitStd : float
1379 Standard deviation of log permanent income at model entry.
1380 PermGroFac : [float]
1381 Lifecycle sequence of permanent income growth factors, *not* offset by
1382 one period as in most other HARK models.
1383 PrstIncCorr : float
1384 Correlation coefficient of the persistent component of income.
1385 pLogCount : int
1386 Number of gridpoints in the grid of (log) persistent income deviations.
1387 pLogRange : float
1388 Upper bound of log persistent income grid, in standard deviations from
1389 the mean; grid has symmetric lower bound.
1391 Returns
1392 -------
1393 IncomeProcessDict : dict
1394 Dictionary with the following entries.
1396 pLogGrid : [np.array]
1397 Age-dependent grids of log persistent income, in deviations from mean.
1398 pLvlMean : [float]
1399 Mean persistent income level by age.
1400 pLogMrkvArray : [np.array]
1401 Age-dependent Markov transition arrays among pLog levels at the start of
1402 each period in the sequence.
1403 """
1404 if cycles == 0:
1405 if T_cycle > 1:
1406 raise ValueError(
1407 "Can't handle infinite horizon models with more than one period!"
1408 )
1409 if PermGroFac[0] != 1.0:
1410 raise ValueError(
1411 "Can't handle permanent income growth in infinite horizon!"
1412 )
1414 # The single pLogGrid and transition matrix can be generated by the basic
1415 # Tauchen AR(1) method from HARK.distributions.
1416 pLogGrid, pLogMrkvArray = make_tauchen_ar1(
1417 pLogCount,
1418 sigma=PermShkStd[0],
1419 ar_1=PrstIncCorr,
1420 bound=pLogRange,
1421 )
1422 pLogGrid = [pLogGrid]
1423 pLogMrkvArray = [pLogMrkvArray]
1424 pLvlMean = [np.exp(pLogInitMean + 0.5 * pLogInitStd**2)]
1426 else:
1427 # Start with the pLog distribution at model entry
1428 pLvlMeanNow = np.exp(pLogInitMean + 0.5 * pLogInitStd**2)
1429 pLogStdNow = pLogInitStd
1430 pLogGridPrev = np.linspace(
1431 -pLogRange * pLogStdNow, pLogRange * pLogStdNow, pLogCount
1432 )
1434 # Initialize empty lists to hold output
1435 pLogGrid = []
1436 pLogMrkvArray = []
1437 pLvlMean = []
1439 for c in range(cycles):
1440 for t in range(T_cycle):
1441 # Update the distribution of persistent income deviations from mean
1442 pLvlMeanNow *= PermGroFac[t]
1443 pLogStdNow = np.sqrt(
1444 (PrstIncCorr * pLogStdNow) ** 2 + PermShkStd[t] ** 2
1445 )
1446 pLogGridNow = np.linspace(
1447 -pLogRange * pLogStdNow, pLogRange * pLogStdNow, pLogCount
1448 )
1450 # Compute transition distances from prior grid to this one
1451 pLogCuts = (pLogGridNow[1:] + pLogGridNow[:-1]) / 2.0
1452 pLogCuts = np.concatenate(([-np.inf], pLogCuts, [np.inf]))
1453 distances = np.reshape(pLogCuts, (1, pLogCount + 1)) - np.reshape(
1454 PrstIncCorr * pLogGridPrev, (pLogCount, 1)
1455 )
1456 distances /= PermShkStd
1458 # Compute transition probabilities, ensuring that very small
1459 # probabilities are treated identically in both directions
1460 cdf_array = norm.cdf(distances)
1461 sf_array = norm.sf(distances)
1462 pLogMrkvNow = cdf_array[:, 1:] - cdf_array[:, :-1]
1463 pLogMrkvNowAlt = sf_array[:, :-1] - sf_array[:, 1:]
1464 pLogMrkvNow = np.maximum(pLogMrkvNow, pLogMrkvNowAlt)
1465 pLogMrkvNow /= np.sum(pLogMrkvNow, axis=1, keepdims=True)
1467 # Add this period's output to the lists
1468 pLogGrid.append(pLogGridNow)
1469 pLogMrkvArray.append(pLogMrkvNow)
1470 pLvlMean.append(pLvlMeanNow)
1471 pLogGridPrev = pLogGridNow
1473 # Gather and return the output
1474 IncomeProcessDict = {
1475 "pLogGrid": pLogGrid,
1476 "pLogMrkvArray": pLogMrkvArray,
1477 "pLvlMean": pLvlMean,
1478 }
1479 return IncomeProcessDict