Summary

This report is part of a series of reports replicating and extending Krebs’ (2007 AER) model of welfare cost of business cycles with uninsurable job displacement risk. Previous reports reproduce the results by calibrating the model using closed-formula solutions (1), and reproduce the results by simulation (2). Here, I depart from the report 2, and add the following features: new mortality options (realistic and cyclical), utility function including VSL, different starting ages, and Great Recession costs estimates.

If the reader wants more basic details of the model, refer to the other reports, as I have chosen not to repeat identical descriptions from one report to the next.

1 Model and simulation steps

The definitions of \(y\), \(\theta\), \(\eta\), and \(\Delta\) are the same from previous versions. With regards to the simulation steps, the only difference is that I do not use a grid search anymore to find \(\Delta\); we opted for the optimization procedure because the results were the same between the two approaches, and optimization is more efficient.

2 Mortality rates

One of the features we add to the model is the use of age-specific mortality rates. To do so, I use data from the SSA, downloaded from this website. I use a population weighted average for men and women in 2007 to calculate unisex mortality rates, aligned with task #106.

To go from mortality rates to survival probabilities, we proceed as follows. The mortality probability \(m_x\) gives us the probability that the individual will die between age \(x\) and \(x+1\). Call the sequences of mortality probablities for a given individual \(\{m_t\}_{t=0}^T\). The survival probability in period \(t\) is given by \(q_t=1-m_t\), and the cumulative survival probability in period \(T\) is calculated recursively, \(Q_T=\prod_{t=0}^{T-1}{q_t}\); \(Q_0=1\) because everyone starts alive.

We can extend this framework to have state-dependent mortality rates: \(m'_t=m_t(1+dm^S)\). For cyclical mortality rates, for example, we could have \(dm^{S=H}=0\Leftrightarrow m'_t=m_t\) if \(S=H\) and \(dm^{S=L}=-0.023 \Leftrightarrow m'_t<m_t\) if \(S=L\), which in words means that the mortality probability is 2.3%. lower on a recession; this effect is independent of age.

As for the values of \(dm^S\), we use 2.3% for Great Recession estimates, and 1.5% for normal recessions, considered in the cycles scenario. We find 1.5% assuming that an average recession increase in unemployment is of 3pp instead of the 4.6 pp of the GR. See reference.

Below, I import mortality probabilities and calculate survival probabilities.

# Import population data processed at scripts/import_seer_population.R
seer_pop <- read_csv("../../data/seer_population.csv") 

# Import life tables by gender
mortality_male <- read_csv("../../data/PerLifeTables_M_Hist_TR2022.csv", skip = 4) %>%
  filter(Year==2007) %>%
  rename(age = "x", mortality = "q(x)") 
mortality_female <- read_csv("../../data/PerLifeTables_F_Hist_TR2022.csv", skip = 4) %>%
  filter(Year==2007) %>%
  rename(age = "x", mortality = "q(x)") 
# Unisex lifetable
life_table <- 
  data.frame(age = 0:119,
             mortality = seer_pop$s_male*mortality_male$mortality +
               seer_pop$s_female*mortality_female$mortality) %>% # unisex mortality
  mutate(q = 1-mortality, # instant survival probability
         Q = NA) # cumulative survival probability, to be filled

# Cumulative survival probability
for (i in 1:nrow(life_table)) {
  life_table$Q[i] <- if(i==1) 1 else cumprod(life_table$q[1:(i-1)])[i-1]
}
# Print table
kable(life_table) %>%
  kable_styling("striped", full_width = F) %>%
  scroll_box(height = "150px")
age mortality q Q
0 0.0067520 0.9932480 1.0000000
1 0.0004646 0.9995354 0.9932480
2 0.0002872 0.9997128 0.9927865
3 0.0002170 0.9997830 0.9925014
4 0.0001746 0.9998254 0.9922860
5 0.0001599 0.9998401 0.9921128
6 0.0001504 0.9998496 0.9919541
7 0.0001408 0.9998592 0.9918049
8 0.0001267 0.9998733 0.9916652
9 0.0001090 0.9998910 0.9915395
10 0.0000949 0.9999051 0.9914314
11 0.0000974 0.9999026 0.9913373
12 0.0001293 0.9998707 0.9912407
13 0.0001991 0.9998009 0.9911125
14 0.0002979 0.9997021 0.9909152
15 0.0004050 0.9995950 0.9906200
16 0.0005092 0.9994908 0.9902188
17 0.0006110 0.9993890 0.9897146
18 0.0007030 0.9992970 0.9891099
19 0.0007858 0.9992142 0.9884145
20 0.0008723 0.9991277 0.9876378
21 0.0009546 0.9990454 0.9867763
22 0.0010119 0.9989881 0.9858343
23 0.0010344 0.9989656 0.9848367
24 0.0010321 0.9989679 0.9838179
25 0.0010188 0.9989812 0.9828026
26 0.0010090 0.9989910 0.9818013
27 0.0010082 0.9989918 0.9808107
28 0.0010196 0.9989804 0.9798218
29 0.0010389 0.9989611 0.9788228
30 0.0010613 0.9989387 0.9778059
31 0.0010850 0.9989150 0.9767682
32 0.0011185 0.9988815 0.9757084
33 0.0011610 0.9988390 0.9746171
34 0.0012153 0.9987847 0.9734855
35 0.0012817 0.9987183 0.9723025
36 0.0013578 0.9986422 0.9710563
37 0.0014551 0.9985449 0.9697377
38 0.0015689 0.9984311 0.9683267
39 0.0017033 0.9982967 0.9668075
40 0.0018520 0.9981480 0.9651608
41 0.0020190 0.9979810 0.9633733
42 0.0022028 0.9977972 0.9614283
43 0.0024008 0.9975992 0.9593104
44 0.0026185 0.9973815 0.9570073
45 0.0028553 0.9971447 0.9545013
46 0.0031093 0.9968907 0.9517760
47 0.0033885 0.9966115 0.9488166
48 0.0036791 0.9963209 0.9456015
49 0.0039955 0.9960045 0.9421226
50 0.0043366 0.9956634 0.9383583
51 0.0047030 0.9952970 0.9342890
52 0.0050714 0.9949286 0.9298951
53 0.0054225 0.9945775 0.9251792
54 0.0057819 0.9942181 0.9201624
55 0.0061811 0.9938189 0.9148422
56 0.0066353 0.9933647 0.9091874
57 0.0071318 0.9928682 0.9031547
58 0.0076604 0.9923396 0.8967135
59 0.0082487 0.9917513 0.8898443
60 0.0089041 0.9910959 0.8825043
61 0.0096390 0.9903610 0.8746464
62 0.0104771 0.9895229 0.8662157
63 0.0114248 0.9885752 0.8571402
64 0.0124912 0.9875088 0.8473476
65 0.0136871 0.9863129 0.8367633
66 0.0150103 0.9849897 0.8253104
67 0.0164108 0.9835892 0.8129223
68 0.0178668 0.9821332 0.7995816
69 0.0194018 0.9805982 0.7852956
70 0.0211568 0.9788432 0.7700594
71 0.0231575 0.9768425 0.7537674
72 0.0253093 0.9746907 0.7363121
73 0.0276386 0.9723614 0.7176765
74 0.0301922 0.9698078 0.6978410
75 0.0331376 0.9668624 0.6767716
76 0.0365071 0.9634929 0.6543450
77 0.0401367 0.9598633 0.6304567
78 0.0440803 0.9559197 0.6051523
79 0.0483903 0.9516097 0.5784770
80 0.0533429 0.9466571 0.5504843
81 0.0589681 0.9410319 0.5211199
82 0.0651976 0.9348024 0.4903905
83 0.0721521 0.9278479 0.4584182
84 0.0799302 0.9200698 0.4253424
85 0.0872866 0.9127134 0.3913447
86 0.0971795 0.9028205 0.3571855
87 0.1082492 0.8917508 0.3224744
88 0.1205498 0.8794502 0.2875668
89 0.1341004 0.8658996 0.2529007
90 0.1489010 0.8510990 0.2189866
91 0.1649388 0.8350612 0.1863793
92 0.1821896 0.8178104 0.1556381
93 0.2006206 0.7993794 0.1272825
94 0.2201907 0.7798093 0.1017470
95 0.2398608 0.7601392 0.0793432
96 0.2593182 0.7406818 0.0603119
97 0.2782237 0.7217763 0.0446719
98 0.2962218 0.7037782 0.0322431
99 0.3129502 0.6870498 0.0226920
100 0.3306297 0.6693703 0.0155905
101 0.3493154 0.6506846 0.0104358
102 0.3690645 0.6309355 0.0067904
103 0.3899379 0.6100621 0.0042843
104 0.4120002 0.5879998 0.0026137
105 0.4353200 0.5646800 0.0015369
106 0.4599690 0.5400310 0.0008678
107 0.4860231 0.5139769 0.0004687
108 0.5135628 0.4864372 0.0002409
109 0.5426745 0.4573255 0.0001172
110 0.5734475 0.4265525 0.0000536
111 0.6059774 0.3940226 0.0000229
112 0.6403657 0.3596343 0.0000090
113 0.6767186 0.3232814 0.0000032
114 0.7151489 0.2848511 0.0000010
115 0.7557773 0.2442227 0.0000003
116 0.7955180 0.2044820 0.0000001
117 0.8352940 0.1647060 0.0000000
118 0.8770580 0.1229420 0.0000000
119 0.9209110 0.0790890 0.0000000 
# plot
ggplot(data=life_table) +
  geom_line(aes(x=age,y=Q,color="Q")) +
  geom_line(aes(x=age,y=q,color="q")) + 
  geom_line(aes(x=age,y=mortality,color="m")) +
  theme_classic() +
  theme(legend.title = element_blank(),
        axis.title.y = element_blank())

3 Shocks and income path

The explanations and code are the same as before for drawing the shocks. There are differences in the income path, explained below.

# Noncyclical component of labor income risk, theta
draw_theta <- function(periods=periods_val, s2=s2_val) {
  lambda <- rnorm(periods,mean=-s2/2,sd=sqrt(s2))
  theta <- exp(lambda) - 1
  return(theta)
}
# Cyclical component of labor income risk, eta
draw_eta <- function(periods=periods_val, pi_l=pi_l_val, d_l=d_l_val, 
                     d_h=d_h_val, p_l=p_l_val, p_h=p_h_val) {
  # Displacement
  displaced <- NULL
  for (i in 1:periods) {
    displaced <- c(displaced,
                   ifelse(low[i] == TRUE, rbernoulli(1, p_l), rbernoulli(1, p_h)))
  }
  
  # Eta
  eta <- NULL
  for (i in 1:periods) {
    eta <- c(eta,
             ifelse(low[i]==TRUE & displaced[i]==TRUE, -d_l,
                    ifelse(low[i]==TRUE & displaced[i]==FALSE, p_l*d_l/(1-p_l),
                           ifelse(low[i]==FALSE & displaced[i]==TRUE, -d_h,
                                  p_h*d_h/(1-p_h)))))
  }
  return(eta)
}
# Eta for case without cycles
draw_eta_bar <- function(periods=periods_val, pi_l=pi_l_val, d_l=d_l_val, 
                         d_h=d_h_val, p_l=p_l_val, p_h=p_h_val, 
                         method=method_val, cost=cost_val) {
  if (!(cost %in% c("cycles","recessions"))) warning("Cost must be \'cycles\' or \'recessions\'")
  
  # Parameters w/o cycles
  pi_h <- 1-pi_l
  p_bar <- ifelse(method == 1, 
                  pi_l*p_l + pi_h*p_h, # Equation 13 from the paper
                  pi_l*p_l + pi_h*p_h) # Equation 12 from the paper -- p_bar is orthogonal to method
  
  if (cost=="cycles") {
    d_bar <- ifelse(method == 1, 
                    pi_l*d_l + pi_h*d_h, # Equation 13 from the paper
                    (pi_l*p_l/p_bar)*d_l + (pi_h*p_h/p_bar)*d_h) # Equation 12 from the paper
  }
  
  if (cost=="recessions") {
    d_bar <- d_h # counterfactual is no recessions
  }
  
  # Displacement
  displaced <- rbernoulli(periods,p_bar)
  
  # Eta
  eta_bar <- NULL
  for (i in 1:periods) {
    eta_bar <- c(eta_bar,
                 ifelse(displaced[i]==1, 
                        -d_bar, p_bar*d_bar/(1-p_bar)))
  }
  return(eta_bar)
}
# Income path
income_path <- function(g=g_val, y0=y0_val, nocycles) {
  y <- c()
  for (i in 1:nrow(shocks_dataset)) {
    if (shocks_dataset$t[i]==0) {
      y[i] <- y0
    }
    else {
      y[i] <- 
        ifelse(nocycles==TRUE,
               (1+g)*(1+shocks_dataset$theta[i])*(1+shocks_dataset$eta_bar[i])*y[i-1],
               (1+g)*(1+shocks_dataset$theta[i])*(1+shocks_dataset$eta[i])*y[i-1])
    }
  }
  return(y)
}

Here, now I include an age variable, which is relevant to define different starting ages for different simulations More specifically, now we can estimate the costs for economies with agents starting at different ages, which can impact the results once we include realistic mortality rates.

For starting ages 0 and 35, the initial income equals one, \(y_0=1\), and the income path follows the process described in the model. For simulations with starting ages greater than 35, however, the initial income equals the average income at that age in the income path defined for starting age 35.

More formally, call the income path for agent \(i\) with age \(x\) in the simulation with starting age 35 \(\{y^{x_0=35}_{i,x}\}_{x=x_0}^{x_T}\) and \(i\in [1,N]\). At age 45, which corresponds to period 10 for the simulation with starting age 35, this agent will have an income level of \(y^{x_0=35}_{i,45}\). For starting ages \(x_0'>35\), we have \(y^{x_0'}_{i,x_0'}=\sum_{i=1}^{N}{y^{x_0=35}_{i,x_0'}}/N\) for all \(i\).

These averages are from the scenario with cycles (i.e. calculated using \(y\) and not \(\overline{y}\)), and realistic mortality rates (these should not affect the non-discounted income path anyway). In terms of coding, what I do is to save the average income at ages 45 and 55 when I calculate the income path for the simulation with starting age 35, realistic mortality, and without great recession. The averages are reported in the end of the first set of tables for cycles’ costs.

Another option would be to have starting incomes for initial ages greater than 35 equal to the income levels each agent had on the simulation with starting age 35. With this approach, we would have different starting incomes across agents, because each agent has its own income path, and they have different income levels when they are at age 35. We do not see a clear reason to choose one approach over the other, so we opted for the first one, since it was easier to code it.

# Create dataset with shocks and income
shocks_income_dataset <- 
  function(N=N_val, periods=periods_val, age_0=age_0_val, age_T=age_T_val) {
    
    # Shocks dataset
    id <- sort(rep(1:N,periods)) # agent_id
    age <- rep(age_0:age_T,N) # agent age
    t <- rep(c(0:(periods-1)),N) # time
    theta <- NULL # empty vector to be filled
    eta <- NULL
    eta_bar <- NULL
    for (i in 1:N) {
      theta <- c(theta,draw_theta())
      eta <- c(eta,draw_eta())
      eta_bar <- c(eta_bar,draw_eta_bar())
    }
    
    # Create dataset with shocks
    shocks_dataset <<- 
      data.frame(id = id,
                 age = age,
                 t = t,
                 theta = theta,
                 eta = eta,
                 eta_bar = eta_bar)
    
    # Income paths
    y <- income_path(nocycles=FALSE)
    y_bar <- income_path(nocycles=TRUE)
    
    # Add income paths to shocks dataset
    shocks_income <- shocks_dataset %>%
      mutate(y = y,
             y_bar = y_bar)
    rm(shocks_dataset, envir=.GlobalEnv) # delete temporary shocks_dataset
    
    return(shocks_income)
  }

4 Utility and squared difference

4.1 Utility function

Below, I define an individual-level utility calculator. For that, I use a CRRA-type utility function, similar to the one given after Equation 6 in the paper, but adding the parameter \(b\). For the case without business cycles, \(c_{it} = \overline{y}_{it} = \overline{c}_{it}\), for the case with business cycles, I calculate the utility with the compensation value, \(c'_{it} = (1+\Delta)y_{it} = (1+\Delta)c_{it}\).

This function returns a vector of utilities by period. These will be discounted and summed in the aggregate utility function calculator defined in the next step.

\[ u(c) = \begin{cases} \frac{c^{1-\gamma}}{1-\gamma}+b& \gamma \ne 1 \\ \log{(c)}+b& \gamma = 1 \end{cases} \]

# Utility function with VSLY
u <- function(gamma,nocycles,Delta) {
  # b parameter
  b <<- b_val_vector[which(gammas_val==gamma)]
  # income variable
  y <- if(nocycles==TRUE){shocks_income$y_bar} else{(1+Delta)*shocks_income$y}
  # utility function
  u_vector <- if(gamma==1){log(y)+b} else{y^(1-gamma)/(1-gamma)+b}
  # checking that u is always positive
  if (!all(u_vector>=0)) {
    warning("There is at least one negative utility.")
  }
  
  return(u_vector)
}

4.2 Aggregate utility

Here, I define a function that calculates the aggregate utility in the economy. In short, I use the individual utility calculator to create a vector of utilities for all individuals in all periods in the economy, then I discount each value by \(\beta^t Q_t\), and sum it across periods and individuals.

The function below allow for four cases of mortality and, consequently, survival rates. Only the first two cases are discussed in Krebs (2007 AER). Note that to calculate \(\overline{U}\), mortality will never be cyclical.

  1. None. In this case, \(m_t=0\) and \(q_t=Q_t=1\).
  2. Constant. In this case, \(m_t=c\), \(q_t=1-c\), and \(Q_t=(1-c)^{t-1}\).
  3. Exogenous. In this case, we follow the values of life_table, defined in the mortality rates section. This is coded as mortality by_age because we have one mortality probability for each age.
  4. Cyclical. In this case, in each period, the exogenous mortality rate is multiplied by \((1+dm^S)\), with \(dm^H=0\), and \(dm^L<0\). Note that \(dm^S\) also depends on the recession type, \(|dm^L_{GR}|>|dm^L_{NR}|\), where NR stands for normal recession. The sequences of \(\{q_t\}_{t=0}^T\) and \(\{Q_t\}_{t=0}^T\) are recalculated using the state-dependent mortality.

Formally, \(U=\sum_{i=1}^{N}{\left[\sum_{t=0}^{T}{(\beta Q_t)^t u(c_{it})}\right]}\), where \(c_{it}\) follows the same definition from above, for the cases with and without business cycles, and \(Q_t\) is defined in the mortality rates section.

# Aggregate utility function
U <- function(gamma, nocycles, Delta, initial_recession=initial_recession_val, # initial recession argument is new here
              mortality=mortality_val, dm_l=dm_l_val,
              pi_l=pi_l_val, beta=beta_val, periods=periods_val, N=N_val,
              age_0=age_0_val, age_T=age_T_val) {
  
  # Mortality options
  if (mortality == "none") {
    Q <- 1
  }
  if (mortality == "constant") {
    q <- parameters$q_val
    Q <- NULL
    for (i in 1:periods) {
      Q[i] <- if(i==1) 1 else q^(i-1)
    }
  }
  # for by_age mortality, it will always be equal to mortality vector
  # for cyclical mortality, it will be mortality vector if no cycles and GR=0
  # if nocycles and GR>0, it will be cyclical post GR 
  if ((mortality %in% c("by_age","cyclical")) & nocycles==T & initial_recession==0 |
      mortality=="by_age" & nocycles==T & initial_recession>0 |
      mortality=="by_age" & nocycles==F) { 
    q <- life_table$q[life_table$age %in% age_0:age_T]
    Q <- NULL
    for (i in 1:periods) {
      Q[i] <- if(i==1) 1 else cumprod(q[1:(i-1)])[i-1]
    }
  }
  if (mortality=="cyclical" & nocycles==T & initial_recession>0) { 
    # State dependent mortality after GR
    m <- life_table$mortality[life_table$age %in% age_0:age_T]
    for (i in (initial_recession+1):periods) {
      m[i] <- if(low[i]==TRUE) (1+dm_l)*m[i] else m[i]
    }
    q <- 1-m # cyclical survival probability 
    Q <- NULL
    for (i in 1:periods) { # cumulative survival probability
      Q[i] <- if(i==1) 1 else cumprod(q[1:(i-1)])[i-1]
    }
  }
  if (mortality=="cyclical" & nocycles == F) {
    # State dependent mortality
    m <- life_table$mortality[life_table$age %in% age_0:age_T]
    for (i in 1:periods) {
      m[i] <- if(low[i]==TRUE) (1+dm_l)*m[i] else m[i]
    }
    q <- 1-m # cyclical survival probability 
    Q <- NULL
    for (i in 1:periods) { # cumulative survival probability
      Q[i] <- if(i==1) 1 else cumprod(q[1:(i-1)])[i-1]
    }
  }
  if (!(mortality %in% c("none","constant","by_age","cyclical"))) {
    warning("Mortality must be \'none\', \'constant\',
            \'by_age\' or \'cyclical\'")
  }
  
  # Vector of utility per period
  per_period_utility <- u(gamma=gamma, nocycles=nocycles, Delta=Delta) 
  
  # Vector of discount rates
  discounting <- rep((beta)^(0:(periods-1))*Q, N) 
  
  # Aggregate utility
  agg_util <- sum(discounting*per_period_utility) # multiply vectors and sum
  
  return(agg_util) 
}

4.3 Squared error

The squared difference function \(g(\Delta)=\left[U(\Delta)-\overline{U}\right]^2\) is defined below.

# Squared difference calculator
diff_squared <- function(gamma, Delta, mortality=mortality_val, dm_l=dm_l_val,
                         beta=beta_val, periods=periods_val, N=N_val) {
  g <- gamma
  U_cycles <- U(gamma=g, nocycles=FALSE, Delta=Delta)
  diff <- (U_cycles - subset(U_bar, gamma==g)$U_bar)^2
  return(diff)
}

5 Optimization

Here, I define the optimization method to find \(\Delta\) such that \(\Delta^*=\underset{\Delta}{\arg\min}\, g(\Delta)\). Note that I add a restriction such that \(\Delta\in[-0.1,0.1]\).

# Optimize Delta
optimize_Deltas <- function(initial_recession=initial_recession_val) {
  optim_Deltas <- NULL
  for (g in gammas_val) { # optimize diff_squared over Delta, given gamma
    optim_Deltas <- c(optim_Deltas,
                      100*optimize(f = function(Delta) diff_squared(Delta, gamma=g), 
                                   interval=c(-0.1,0.1), tol = 1e-14)$minimum)
  }
  output <- data.frame(gamma=gammas_val, optimization=optim_Deltas) %>%
    mutate_if(is.numeric, ~round(.,3))
  
  return(output)
}

Calibration results

Same as before.

# Deltas from calibration
Deltas_calibration <- function(cost=cost_val,mortality=mortality_val,
                               initial_recession=0) {
  # Calibration results for beta=.96
  if (cost=="cycles" & mortality=="none" & parameters$beta_val==.96) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, .527,
                        ifelse(gammas_val[i]==1.5, .735,
                        ifelse(gammas_val[i]==2, .982,
                        ifelse(gammas_val[i]==2.5, 1.309,
                        ifelse(gammas_val[i]==3, 1.787,
                        ifelse(gammas_val[i]==3.5, 2.570,
                        ifelse(gammas_val[i]==4, 4.092, NA)))))))
    }}
  if (cost=="cycles" & mortality=="constant" & parameters$beta_val==.96) {
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 0.325, 
                        ifelse(gammas_val[i]==1.5, 0.476,
                        ifelse(gammas_val[i]==2, 0.648,
                        ifelse(gammas_val[i]==2.5, 0.863,
                        ifelse(gammas_val[i]==3, 1.153,
                        ifelse(gammas_val[i]==3.5, 1.573,
                        ifelse(gammas_val[i]==4, 2.244, NA)))))))
    }}
  if (cost=="recessions" & mortality=="none" & parameters$beta_val==.96) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 1.342, 
                        ifelse(gammas_val[i]==1.5, 1.852,
                        ifelse(gammas_val[i]==2, 2.425,
                        ifelse(gammas_val[i]==2.5, 3.165,
                        ifelse(gammas_val[i]==3, 4.217,
                        ifelse(gammas_val[i]==3.5, 5.888,
                        ifelse(gammas_val[i]==4, 9.000, NA)))))))
    }}
  if (cost=="recessions" & mortality=="constant" & parameters$beta_val==.96) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 0.831, 
                        ifelse(gammas_val[i]==1.5, 1.199,
                        ifelse(gammas_val[i]==2, 1.602,
                        ifelse(gammas_val[i]==2.5, 2.091,
                        ifelse(gammas_val[i]==3, 2.731,
                        ifelse(gammas_val[i]==3.5, 3.636,
                        ifelse(gammas_val[i]==4, 5.038, NA)))))))
    }}
  
  # Calibration results for beta=.99
  if (cost=="cycles" & mortality=="none" & parameters$beta_val==.99) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 2.168, 
                        ifelse(gammas_val[i]==1.5, 2.225,
                        ifelse(gammas_val[i]==2, 2.677,
                        ifelse(gammas_val[i]==2.5, 3.593,
                        ifelse(gammas_val[i]==3, 5.559,
                        ifelse(gammas_val[i]==3.5, 11.809, NA))))))
    }}
  if (cost=="cycles" & mortality=="constant" & parameters$beta_val==.99) {
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 0.627, 
                        ifelse(gammas_val[i]==1.5, 0.858,
                        ifelse(gammas_val[i]==2, 1.135,
                        ifelse(gammas_val[i]==2.5, 1.514,
                        ifelse(gammas_val[i]==3, 2.089,
                        ifelse(gammas_val[i]==3.5, 3.084,
                        ifelse(gammas_val[i]==4, 5.222, NA)))))))
    }}
  if (cost=="recessions" & mortality=="none" & parameters$beta_val==.99) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 5.536, 
                        ifelse(gammas_val[i]==1.5, 5.636,
                        ifelse(gammas_val[i]==2, 6.612,
                        ifelse(gammas_val[i]==2.5, 8.621,
                        ifelse(gammas_val[i]==3, 12.823,
                        ifelse(gammas_val[i]==3.5, 25.251, NA))))))
    }}
  if (cost=="recessions" & mortality=="constant" & parameters$beta_val==.99) {
    # Vector of Deltas from calibration
    Deltas_cal <- rep(NA,length(gammas_val))
    for (i in 1:length(gammas_val)) {
      Deltas_cal[i] <-  ifelse(gammas_val[i]==1.01, 1.601, 
                        ifelse(gammas_val[i]==1.5, 2.161,
                        ifelse(gammas_val[i]==2, 2.803,
                        ifelse(gammas_val[i]==2.5, 3.658,
                        ifelse(gammas_val[i]==3, 4.921,
                        ifelse(gammas_val[i]==3.5, 7.033,
                        ifelse(gammas_val[i]==4, 11.353, NA)))))))
    }}
  # non-replication cases
  if (initial_recession!=0 | (mortality %in% c("by_age","cyclical)"))) {
    Deltas_cal <- NA
  }
  return(Deltas_cal)
}

6 Simulation

6.1 Parameters

If ones want to try different parameters, there is no need to touch the rest of the code before here. In this version we are using a higher displacement probability for the GR scenario, as in this report.

The number of periods on the simulation is defined by the difference between a starting age for the agents, and a final age. We vary the starting age, and use 100 for the final age in most cases. For none and constant mortality, however, we extend the periods to 150, regardless of the starting age. We do so, to be able to make a closer comparison to Krebs’ results, with infinitely lived agents.

parameters <- data.frame(
  age_0_val = 0, # initial age
  N_val = 1000, # agents in the economy
  s2_val = 0.01, # sigma squared
  pi_l_val = 0.5, # low economy state probability
  d_l_val = 0.21, # income loss if displaced in low
  d_h_val = 0.09, # income loss if displaced in high
  p_l_rec_val = 0.05, # displacement probability in normal recession
  p_h_val = 0.03, # displacement probability in high state
  method_val = 1, # method to calculate p and d w/o cycles
  beta_val = 0.99, # beta
  g_val = 0.02, # growth rate
  y0_val = 1, # initial income
  q_val = 0.976, # constant survival probability
  dm_l_gr_val = -.023, # drop in mortality in great recession
  rounds_val = 200) %>% # number of draws 
  mutate(dm_l_rec_val = dm_l_gr_val*3.1/4.6,
         p_l_gr_val = p_l_rec_val*4.6/3.1) # smaller dm for normal recessions
# 3.1 and 4.6 are average unemployment rate increases for regular recessions
# and the GR
gammas_val <- c(1.5,2,2.5) # constant degree of relative risk aversion
csv_id <- "mort_"

To facilitate visualization and coding, I define a table with different values of \(b\) for each \(\gamma\).

# b values
b_values <- data.frame(
  gamma = gammas_val,
  b0 = rep(0,length(gammas_val)),
  b1 = c(3.563483225,2.380952381,1.88544086),
  b2 = c(6.236095645,4.761904762,4.006561828),
  b3 = c(8.908708064,7.142857143,6.127682795)
)

# Print table
kable(b_values) %>%
  kable_styling("striped", full_width = F)
gamma b0 b1 b2 b3
1.5 0 3.563483 6.236096 8.908708
2.0 0 2.380952 4.761905 7.142857
2.5 0 1.885441 4.006562 6.127683

6.2 One draw

This function uses all the other ones defined above, and the chosen parameters to run the simulation with one click for a given randomization seed. Note that there is a new parameter here: initial_recession. This allows us to estimate Great Recession costs separately, where we start with a recession for a given number of periods, and stay in the high state for the rest of the time.

simulate <- 
  function(seed, parameters, cost, mortality, initial_recession, sanity_check=F) {
    
    # Set parameters
    age_0_val <<- parameters$age_0_val
    age_T_val <<- ifelse(mortality%in%c("none","constant"),age_0_val+150,100)#longer lived agents if no mortality, to replicate krebs
    T_val <<- age_T_val - age_0_val
    periods_val <<- T_val+1
    N_val <<- parameters$N_val
    s2_val <<- parameters$s2_val
    pi_l_val <<- parameters$pi_l_val
    d_l_val <<- parameters$d_l_val
    d_h_val <<- parameters$d_h_val
    p_l_val <<- ifelse(initial_recession==0,parameters$p_l_rec_val,parameters$p_l_gr_val)
    p_h_val <<- parameters$p_h_val
    method_val <<- parameters$method_val
    beta_val <<- parameters$beta_val
    g_val <<- parameters$g_val
    y0_val <<- ifelse(age_0_val<=35,parameters$y0_val,get(paste0("y0_age",age_0_val))) # initial income depends on age
    mortality_val <<- mortality
    cost_val <<- cost
    initial_recession_val <<- initial_recession
    dm_l_val <<- ifelse(initial_recession==0,parameters$dm_l_rec_val,parameters$dm_l_gr_val)
    values <<- ls(pattern="*_val$",envir=.GlobalEnv)
    
    # Stop if age_0 \nin {0,initial_ages}
    if (!(age_0_val%in%c(0,ages))) {
      stop("Age is not in initial ages vector")
    }
    
    # Calibration Deltas for replication
    Deltas_cal <<- Deltas_calibration(initial_recession=initial_recession_val)
    
    # Set seed
    seed_val <<- seed
    set.seed(seed_val)
    
    # Economy state S
    low <<- if(initial_recession==0) rbernoulli(periods_val,pi_l_val) else
      c(rep(TRUE,initial_recession),rbernoulli(periods_val-initial_recession,pi_l_val)) 
    # this is used to draw_eta and cyclical mortality
    
    # Dataset with shocks and income path
    shocks_income <<- shocks_income_dataset()
    
    # Save income averages for ages>35
    if (age_0_val==35 & mortality_val=="by_age" & initial_recession_val==0 & parameters$beta_val==.99) {
      avg_income_by_age <- shocks_income %>% group_by(age) %>% 
        summarise(y_avg = mean(y), y_bar_avg = mean(y_bar))
      for (i in ages[ages>35]) {
        assign(paste0("y0_age",i),
               avg_income_by_age %>% filter(age==i) %>% 
                 select(y_avg) %>% as.numeric(),
               envir = .GlobalEnv)
      }
    }
    
    # Loop over possible values of VSLY
    i0 <- ifelse(mortality=="cyclical",1,0) # starts in b1 for endogenous mort
    imax <- ifelse(mortality=="cyclical",3,0) # stops in b0 if non-endogenous
    if (sanity_check==T) { # allows different VSLYs regardless of mortality, for sanity check
      i0 <- 1
      imax <- 3
    }
    output <- data.frame(gamma = gammas_val) # empty output table
    for (i in i0:imax) {
      # Extract vector of b_values
      b_val_vector <<- b_values[,paste0("b",i)]
      # Utility w/o business cycles
      U_bar <<-
        foreach(g=gammas_val, .combine=rbind, 
                .packages=c("tidyverse","foreach","doParallel"),
                .export=c("U","u","shocks_income","parameters","low","life_table",
                          "b_val_vector",values)) %dopar% {
                            data.frame(gamma=g, U_bar=U(gamma=g, nocycles=TRUE))
                          }
      
      # Run optimization and assign to output
      output[,paste0("b",i)] <- optimize_Deltas()$optimization
    }
    
    # Remove temporary objects
    rm(list = c("shocks_income","U_bar","low",values[values!="gammas_val"]), envir=.GlobalEnv)
    rm(values, envir=.GlobalEnv)
    
    return(output)
  }

6.3 Multiple draws

Same as before, but with initial recession option.

# Simulate multiple draws
simulate_multi <- function(cost, mortality, initial_recession, 
                           rounds=parameters$rounds_val, sanity_check=F) {
  
  # Simulate multiple times and append
  appended_results <- NULL
  for (i in 1:rounds) {
    appended_results <- 
      bind_rows(appended_results,
                simulate(seed=i, parameters=parameters, cost=cost,
                         mortality=mortality, initial_recession=initial_recession,
                         sanity_check=sanity_check))
  }
  
  # Calculate averages
  if (length(names(appended_results))>2) {
    output <- appended_results %>% group_by(gamma) %>%
      summarise(b1=mean(b1), b2=mean(b2), b3=mean(b3)) 
  }
  else {
    output <- appended_results %>% group_by(gamma) %>%
      summarise(b0=mean(b0))
  }
  
  return(output)
}

6.4 Table results

The following function generates results for a given starting age, with none and constant mortality, looking at cycles and recessions. If we run this function with \(\beta=0.96\), we should replicate Krebs tables 1, 2, and mortality results.

results_aggregate <- function(initial_recession,age) {
  # Specify initial age
  parameters$age_0_val <<- age
  # Simulations list
  simulations_list <- list(
    data.frame(gamma=gammas_val, cal_cycles_none=
                 Deltas_calibration(cost='cycles',mortality='none',initial_recession=0)),
    simulate_multi(cost='cycles',mortality='none',initial_recession=0) %>% rename(cycles_none=b0),
    data.frame(gamma=gammas_val, cal_cycles_constant=
                 Deltas_calibration(cost='cycles',mortality='constant',initial_recession=0)),
    simulate_multi(cost='cycles',mortality='constant',initial_recession=0) %>% rename(cycles_constant=b0),
    data.frame(gamma=gammas_val, cal_recessions_none=
                 Deltas_calibration(cost='recessions',mortality='none',initial_recession=initial_recession)),
    simulate_multi(cost='recessions',mortality='none',initial_recession=initial_recession) %>% rename(recessions_none=b0),
    data.frame(gamma=gammas_val, cal_recessions_constant=
                 Deltas_calibration(cost='recessions',mortality='constant',initial_recession=initial_recession)),
    simulate_multi(cost='recessions',mortality='constant',initial_recession=initial_recession) %>% rename(recessions_constant=b0))
  # Merge results
  merged_results <- simulations_list %>% 
    reduce(left_join, by="gamma") %>% # merge
    round(3) # round to 3 decimals
  # Save csv
  beta96 <- if(parameters$beta_val==.96) "_beta96" else NULL
  write.csv(merged_results, row.names = FALSE, file = 
              paste0("../../tables/",csv_id,"GR",initial_recession,beta96,".csv"))
  # Format table to print in report
  table <- kable(rbind(merged_results,c("Mortality",rep(c("None","None","Constant","Constant"),2))), 
                 col.names=c("gamma",rep(c("calibration","simulation"),4)),
                 caption = paste0("Results for starting age ", age)) %>%
    kable_styling("striped", full_width = F) %>%
    add_header_above(c(" ","Costs of cycles"=4, setNames(4,paste0("Costs of GR",initial_recession)))) 
  
  return(table)
}

The function below outputs tables for a given set of starting ages. Each panel will have four columns, one with realistic exogenous mortality, and three with endogenous mortality and different VSLYs. The parameters allow for costs of recessions or cycles, with or without great recession.

results_by_age <- function(cost, initial_recession, ages) {
  ages <<- ages
  results_list <- NULL
  for (age in ages) {
    # Specify initial age
    parameters$age_0_val <<- age
    # Age-specific results
    results_by_age <- left_join(
      simulate_multi(cost=cost,mortality='by_age',initial_recession=initial_recession) %>%
        rename(by_age=b0),
      simulate_multi(cost=cost,mortality='cyclical',initial_recession=initial_recession),
      by="gamma") %>%
      round(3) %>% # round to 3 decimals
      mutate(age0 = age)
    # Save age specific result
    results_list[[which(ages==age)]] <- results_by_age
  }
  # Append age-specific results
  appended_ages <- reduce(results_list,bind_rows) %>% relocate(age0, .after=gamma)
  # Save csv
  beta96 <- if(parameters$beta_val==.96) "_beta96" else NULL
  write.csv(appended_ages, row.names = FALSE, file = 
              paste0("../../tables/",csv_id,cost,"_by_age_","GR",
                     initial_recession,beta96,".csv"))
  
  # Format table to print in report
  title_cost <- ifelse(cost=='cycles', 'cycles', 
                       ifelse(initial_recession==0, 'recessions',
                              ifelse(initial_recession==5, 'GR5',
                                     ifelse(initial_recession==10, 'GR10'))))
  table <- kable(rbind(appended_ages,c("VSLY","-","-","100K","250K","450K")), 
                 col.names=c("gamma","initial age","exogenous",rep("endogenous",3)),
                 caption = paste0("Welfare costs of ",title_cost, 
                                  " by age: full dynamic model")) %>%
    kable_styling("striped", full_width = F) %>%
    add_header_above(c(" "=2,"Realistic mortality"=4)) %>%
    scroll_box(height = "500px")
  
  return(table)
}

7 Main results

7.1 Recessions by age

This should be run before the GR by age, because we use this simulation to save the average income levels for ages greater than 35.

start <- Sys.time()
results_by_age(cost='recessions',initial_recession=0,ages=c(35,45,55,65))
Welfare costs of recessions by age: full dynamic model
Realistic mortality
gamma initial age exogenous endogenous endogenous endogenous
1.5 35 1.742 1.404 0.956 0.512
2 35 2.364 2.057 1.628 1.203
2.5 35 3.091 2.826 2.435 2.049
1.5 45 1.457 0.993 0.401 -0.185
2 45 1.995 1.531 0.914 0.305
2.5 45 2.62 2.172 1.556 0.95
1.5 55 1.2 0.559 -0.229 -1.008
2 55 1.63 0.93 0.042 -0.831
2.5 55 2.119 1.377 0.409 -0.536
1.5 65 0.919 0.002 -1.079 -2.143
2 65 1.256 0.168 -1.153 -2.439
2.5 65 1.639 0.371 -1.198 -2.704
VSLY
100K 250K 450K 
Sys.time()-start # time difference

Time difference of 47.08593 mins

Initial income for ages 45, 55, and 65 are, respectively, 1.21, 1.47, and 1.79.

7.2 GR10 by age

start <- Sys.time()
results_by_age(cost='recessions',initial_recession=10,ages=c(35,45,55,65))
Welfare costs of GR10 by age: full dynamic model
Realistic mortality
gamma initial age exogenous endogenous endogenous endogenous
1.5 35 3.919 3.848 3.737 3.625
2 35 5.353 5.29 5.182 5.075
2.5 35 7.03 6.978 6.881 6.784
1.5 45 3.374 3.186 2.922 2.659
2 45 4.631 4.44 4.163 3.887
2.5 45 6.087 5.904 5.625 5.347
1.5 55 2.94 2.512 1.956 1.404
2 55 4.022 3.549 2.917 2.293
2.5 55 5.261 4.753 4.059 3.377
1.5 65 2.352 1.318 0.068 -1.159
2 65 3.239 1.999 0.466 -1.022
2.5 65 4.251 2.793 0.967 -0.78
VSLY
100K 250K 450K 
Sys.time()-start # time difference

Time difference of 46.04236 mins

8 Sanity check

8.1 Replication with different age

The replication should be independent of the starting age being 0 or 35, because the periods are extended when we have constant or no mortality. Therefore, if everything is correct, this table should be identical to the previous one.

start <- Sys.time()
results_aggregate(initial_recession=0,age=35)
Results for starting age 35
Costs of cycles
Costs of GR0
gamma calibration simulation calibration simulation calibration simulation calibration simulation
1.5 2.225 1.682 0.858 0.837 5.636 4.229 2.161 2.117
2 2.677 2.15 1.135 1.1 6.612 5.341 2.803 2.744
2.5 3.593 2.798 1.514 1.455 8.621 6.837 3.658 3.557
Mortality None None Constant Constant None None Constant Constant 
Sys.time()-start # time difference

Time difference of 56.40248 mins

8.2 Exogenous for different b’s

The estimates with exogenous mortality should be independent of the VSLY choice, because when we optimize the squared difference function, the parameter \(b\) cancel out. Below, I code one example to check that this is the case.

kable(rbind(simulate_multi(cost='recessions',mortality='by_age',
                     initial_recession=0,sanity_check=T) %>% round(3),
            c("VSLY","100K","250K","450K")),
      caption = "Costs of recessions with exogenous mortality and different VSLYs") %>%
  kable_styling("striped", full_width=F)
Costs of recessions with exogenous mortality and different VSLYs
gamma b1 b2 b3
1.5 1.742 1.742 1.742
2 2.364 2.364 2.364
2.5 3.091 3.091 3.091
VSLY 100K 250K 450K