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-form solutions (1), reproduce the results by simulation (2), add realistic mortality rates with cyclical components and estimate Great Recession costs (3), and add retirement to the model (4). Here, I depart from report 4 and use a higher displacement probability in the Great Recession scenario.

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

The income path follows the same law of motion as in the previous report. The definitions of \(\theta\) and \(\Delta\), as well as the income law of motion are identical to the previous report. The only difference here lies on the value of \(\eta\) when in the Great Recession scenario, such that \(p_L^{Great\,Recession}>p_L^{Normal\,recession}\).

2 Mortality rates

# 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]
}

3 Shocks and income path

# 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)
}

3.1 Income path

# Income path calculator with retirement
income_path <- function(nocycles, initial_recession=initial_recession_val,
                        g=g_val, y0=y0_val) {
  y <- NULL
  for (i in 1:nrow(shocks_dataset)) {
    # First period, same for all cases
    if (shocks_dataset$t[i]==0) {
      y[i] <- y0
    }
    # Retirement, same for all cases
    if (shocks_dataset$t[i]!=0 & shocks_dataset$age[i]>65) { # retirement 
      y[i] <- y[i-1]
    }
    # No GR, same as Krebs
    if (initial_recession==0 & shocks_dataset$age[i]<=65) {
      # everything before retirement
      if (shocks_dataset$t[i]!=0) {
        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])
      }
    }
    # GR, different between during and post GR
    if(initial_recession>0 & shocks_dataset$age[i]<=65) {
      # during GR
      if (shocks_dataset$t[i] %in% 1:(initial_recession-1)) { # same as initial_rec=0
        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])
      }
      # post GR, before retirement
      if (shocks_dataset$t[i]>(initial_recession-1)) { # all random, no nocycles/eta_bar
        y[i] <- (1+g)*(1+shocks_dataset$theta[i])*(1+shocks_dataset$eta[i])*y[i-1]
      }
    }
  }
  return(y)
}

# 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

Same as previous report.

# 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)
}

# 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) 
}

# 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

Same as previous report.

# 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)
}

# 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

The parameters are the same from before, but here we allow for recession-specific \(p_L\)’s.

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 <- "gr_pl_"

# 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	b1	b2	b3
1.5	3.563483	6.236096	8.908708
2.0	2.380952	4.761905	7.142857
2.5	1.885441	4.006562	6.127683

6.2 Simulation and output

The simulation function here allows for a recession-specific displacement probability.

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)
  }

# 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)
}

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 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=seq(35,75,1))

Welfare costs of recessions by age: full dynamic model
		Realistic mortality
gamma	initial age	exogenous	endogenous	endogenous	endogenous
1.5	35	1.517	1.207	0.768	0.331
2	35	2.09	1.807	1.389	0.975
2.5	35	2.74	2.492	2.113	1.738
1.5	36	1.5	1.182	0.731	0.283
2	36	2.063	1.769	1.338	0.911
2.5	36	2.707	2.447	2.053	1.662
1.5	37	1.522	1.194	0.731	0.271
2	37	2.078	1.772	1.326	0.883
2.5	37	2.706	2.433	2.02	1.612
1.5	38	1.441	1.105	0.63	0.159
2	38	1.981	1.664	1.205	0.749
2.5	38	2.592	2.307	1.88	1.457
1.5	39	1.403	1.055	0.567	0.082
2	39	1.939	1.608	1.131	0.657
2.5	39	2.553	2.251	1.802	1.357
1.5	40	1.376	1.016	0.515	0.017
2	40	1.899	1.555	1.06	0.571
2.5	40	2.495	2.178	1.709	1.245
1.5	41	1.271	0.901	0.386	-0.125
2	41	1.771	1.413	0.902	0.396
2.5	41	2.338	2.005	1.515	1.032
1.5	42	1.245	0.862	0.333	-0.193
2	42	1.722	1.348	0.818	0.293
2.5	42	2.259	1.907	1.395	0.889
1.5	43	1.246	0.85	0.305	-0.236
2	43	1.716	1.325	0.774	0.229
2.5	43	2.24	1.867	1.329	0.798
1.5	44	1.229	0.819	0.258	-0.299
2	44	1.687	1.277	0.704	0.137
2.5	44	2.201	1.806	1.24	0.682
1.5	45	1.115	0.693	0.117	-0.454
2	45	1.556	1.131	0.539	-0.047
2.5	45	2.052	1.638	1.048	0.467
1.5	46	1.092	0.659	0.068	-0.517
2	46	1.518	1.079	0.468	-0.135
2.5	46	1.995	1.564	0.953	0.35
1.5	47	1.165	0.715	0.106	-0.497
2	47	1.589	1.129	0.493	-0.134
2.5	47	2.062	1.605	0.961	0.328
1.5	48	1.015	0.552	-0.074	-0.693
2	48	1.406	0.927	0.269	-0.381
2.5	48	1.84	1.358	0.686	0.025
1.5	49	1.013	0.535	-0.107	-0.743
2	49	1.405	0.909	0.228	-0.443
2.5	49	1.843	1.339	0.639	-0.049
1.5	50	0.97	0.473	-0.19	-0.846
2	50	1.335	0.812	0.102	-0.598
2.5	50	1.74	1.202	0.462	-0.265
1.5	51	0.894	0.384	-0.295	-0.968
2	51	1.236	0.695	-0.037	-0.758
2.5	51	1.613	1.052	0.284	-0.469
1.5	52	0.849	0.323	-0.376	-1.067
2	52	1.177	0.614	-0.144	-0.89
2.5	52	1.537	0.949	0.147	-0.64
1.5	53	0.791	0.252	-0.463	-1.171
2	53	1.092	0.513	-0.266	-1.032
2.5	53	1.422	0.813	-0.016	-0.827
1.5	54	0.73	0.173	-0.563	-1.291
2	54	1.015	0.411	-0.397	-1.192
2.5	54	1.331	0.688	-0.18	-1.029
1.5	55	0.681	0.107	-0.65	-1.398
2	55	0.941	0.314	-0.522	-1.344
2.5	55	1.226	0.553	-0.352	-1.236
1.5	56	0.593	0	-0.779	-1.548
2	56	0.831	0.179	-0.688	-1.539
2.5	56	1.093	0.386	-0.56	-1.483
1.5	57	0.536	-0.079	-0.882	-1.675
2	57	0.756	0.072	-0.83	-1.716
2.5	57	0.997	0.247	-0.749	-1.719
1.5	58	0.506	-0.13	-0.957	-1.774
2	58	0.704	-0.01	-0.946	-1.865
2.5	58	0.919	0.129	-0.914	-1.929
1.5	59	0.447	-0.209	-1.06	-1.9
2	59	0.619	-0.122	-1.091	-2.041
2.5	59	0.807	-0.021	-1.108	-2.166
1.5	60	0.38	-0.297	-1.172	-2.036
2	60	0.531	-0.238	-1.24	-2.222
2.5	60	0.698	-0.168	-1.302	-2.403
1.5	61	0.303	-0.398	-1.3	-2.189
2	61	0.422	-0.381	-1.422	-2.441
2.5	61	0.552	-0.362	-1.55	-2.703
1.5	62	0.23	-0.492	-1.419	-2.332
2	62	0.324	-0.508	-1.583	-2.634
2.5	62	0.426	-0.527	-1.762	-2.958
1.5	63	0.157	-0.587	-1.541	-2.48
2	63	0.22	-0.643	-1.754	-2.839
2.5	63	0.289	-0.708	-1.993	-3.236
1.5	64	0.081	-0.687	-1.667	-2.633
2	64	0.109	-0.786	-1.935	-3.057
2.5	64	0.14	-0.903	-2.241	-3.534
1.5	65	0	-0.793	-1.802	-2.796
2	65	0	-0.932	-2.122	-3.283
2.5	65	0	-1.096	-2.494	-3.843
1.5	66	0	-0.813	-1.848	-2.867
2	66	0	-0.956	-2.176	-3.365
2.5	66	0	-1.123	-2.556	-3.937
1.5	67	0	-0.834	-1.895	-2.938
2	67	0	-0.98	-2.23	-3.448
2.5	67	0	-1.152	-2.619	-4.032
1.5	68	0	-0.855	-1.942	-3.011
2	68	0	-1.005	-2.286	-3.532
2.5	68	0	-1.18	-2.684	-4.129
1.5	69	0	-0.876	-1.99	-3.085
2	69	0	-1.03	-2.342	-3.618
2.5	69	0	-1.21	-2.749	-4.227
1.5	70	0	-0.899	-2.04	-3.162
2	70	0	-1.056	-2.4	-3.706
2.5	70	0	-1.24	-2.816	-4.328
1.5	71	0	-0.921	-2.091	-3.239
2	71	0	-1.082	-2.459	-3.796
2.5	71	0	-1.271	-2.885	-4.431
1.5	72	0	-0.944	-2.142	-3.318
2	72	0	-1.109	-2.518	-3.886
2.5	72	0	-1.302	-2.954	-4.535
1.5	73	0	-0.967	-2.194	-3.398
2	73	0	-1.136	-2.579	-3.979
2.5	73	0	-1.334	-3.024	-4.64
1.5	74	0	-0.991	-2.248	-3.48
2	74	0	-1.164	-2.641	-4.073
2.5	74	0	-1.366	-3.095	-4.747
1.5	75	0	-1.015	-2.302	-3.563
2	75	0	-1.192	-2.704	-4.168
2.5	75	0	-1.399	-3.168	-4.856
VSLY			100K	250K	450K

Sys.time()-start # time difference

Time difference of 8.657715 hours

7.2 GR10 by age

start <- Sys.time()
results_by_age(cost='recessions',initial_recession=10,ages=seq(35,75,1))

Welfare costs of GR10 by age: full dynamic model
		Realistic mortality
gamma	initial age	exogenous	endogenous	endogenous	endogenous
1.5	35	2.216	2.145	2.038	1.931
2	35	3.06	2.996	2.895	2.794
2.5	35	4.007	3.952	3.862	3.772
1.5	36	2.165	2.087	1.971	1.854
2	36	3.007	2.936	2.826	2.716
2.5	36	3.949	3.888	3.788	3.689
1.5	37	2.213	2.128	2.001	1.874
2	37	3.054	2.976	2.854	2.733
2.5	37	3.989	3.92	3.809	3.699
1.5	38	2.199	2.105	1.967	1.829
2	38	3.043	2.956	2.823	2.69
2.5	38	3.986	3.909	3.786	3.665
1.5	39	2.191	2.088	1.936	1.785
2	39	3.022	2.925	2.777	2.63
2.5	39	3.95	3.862	3.725	3.589
1.5	40	2.07	1.957	1.792	1.627
2	40	2.874	2.766	2.604	2.442
2.5	40	3.768	3.67	3.518	3.367
1.5	41	2.183	2.058	1.877	1.697
2	41	3.005	2.885	2.707	2.529
2.5	41	3.918	3.808	3.639	3.47
1.5	42	2.194	2.057	1.86	1.663
2	42	3.025	2.892	2.695	2.499
2.5	42	3.95	3.826	3.637	3.449
1.5	43	2.096	1.945	1.73	1.516
2	43	2.899	2.751	2.534	2.318
2.5	43	3.782	3.642	3.432	3.223
1.5	44	2.169	2.003	1.769	1.536
2	44	2.986	2.822	2.583	2.345
2.5	44	3.9	3.742	3.508	3.275
1.5	45	2.111	1.93	1.677	1.424
2	45	2.924	2.742	2.482	2.223
2.5	45	3.825	3.65	3.392	3.136
1.5	46	2.183	1.988	1.714	1.441
2	46	3.007	2.81	2.527	2.245
2.5	46	3.928	3.735	3.453	3.173
1.5	47	2.117	1.904	1.608	1.313
2	47	2.937	2.719	2.41	2.103
2.5	47	3.848	3.633	3.321	3.011
1.5	48	2.144	1.913	1.593	1.275
2	48	2.957	2.718	2.382	2.048
2.5	48	3.857	3.619	3.276	2.936
1.5	49	2.073	1.823	1.48	1.138
2	49	2.873	2.612	2.248	1.887
2.5	49	3.754	3.491	3.117	2.745
1.5	50	2.005	1.733	1.363	0.994
2	50	2.801	2.515	2.117	1.723
2.5	50	3.679	3.385	2.971	2.561
1.5	51	2.035	1.742	1.344	0.949
2	51	2.843	2.532	2.102	1.676
2.5	51	3.737	3.414	2.962	2.515
1.5	52	2.051	1.734	1.307	0.883
2	52	2.853	2.514	2.049	1.589
2.5	52	3.744	3.39	2.897	2.409
1.5	53	2.073	1.733	1.276	0.821
2	53	2.879	2.513	2.012	1.516
2.5	53	3.77	3.384	2.849	2.321
1.5	54	2.064	1.696	1.202	0.713
2	54	2.858	2.457	1.912	1.373
2.5	54	3.733	3.305	2.717	2.137
1.5	55	2.063	1.663	1.132	0.604
2	55	2.848	2.41	1.819	1.235
2.5	55	3.713	3.241	2.598	1.964
1.5	56	1.998	1.565	0.99	0.421
2	56	2.778	2.299	1.656	1.02
2.5	56	3.637	3.115	2.407	1.711
1.5	57	1.895	1.422	0.8	0.183
2	57	2.609	2.08	1.376	0.682
2.5	57	3.391	2.808	2.026	1.258
1.5	58	1.586	1.071	0.399	-0.267
2	58	2.205	1.626	0.86	0.105
2.5	58	2.88	2.235	1.377	0.536
1.5	59	1.38	0.821	0.094	-0.625
2	59	1.931	1.298	0.465	-0.354
2.5	59	2.53	1.82	0.88	-0.039
1.5	60	1.228	0.623	-0.162	-0.938
2	60	1.695	1.006	0.103	-0.783
2.5	60	2.205	1.426	0.401	-0.598
1.5	61	1.001	0.342	-0.508	-1.347
2	61	1.389	0.631	-0.353	-1.319
2.5	61	1.812	0.947	-0.181	-1.278
1.5	62	0.733	0.02	-0.897	-1.802
2	62	1.024	0.2	-0.865	-1.909
2.5	62	1.339	0.393	-0.834	-2.025
1.5	63	0.53	-0.242	-1.232	-2.207
2	63	0.735	-0.161	-1.315	-2.444
2.5	63	0.958	-0.078	-1.415	-2.708
1.5	64	0.265	-0.571	-1.638	-2.689
2	64	0.37	-0.606	-1.857	-3.077
2.5	64	0.482	-0.655	-2.113	-3.518
1.5	65	0	-0.906	-2.057	-3.188
2	65	0	-1.064	-2.42	-3.739
2.5	65	0	-1.25	-2.841	-4.369
1.5	66	0	-0.975	-2.213	-3.427
2	66	0	-1.145	-2.601	-4.015
2.5	66	0	-1.345	-3.051	-4.686
1.5	67	0	-1.049	-2.379	-3.683
2	67	0	-1.232	-2.795	-4.31
2.5	67	0	-1.446	-3.276	-5.023
1.5	68	0	-1.128	-2.557	-3.955
2	68	0	-1.324	-3.002	-4.623
2.5	68	0	-1.554	-3.514	-5.38
1.5	69	0	-1.213	-2.747	-4.245
2	69	0	-1.423	-3.222	-4.956
2.5	69	0	-1.669	-3.768	-5.759
1.5	70	0	-1.304	-2.95	-4.555
2	70	0	-1.529	-3.457	-5.311
2.5	70	0	-1.793	-4.039	-6.162
1.5	71	0	-1.401	-3.167	-4.886
2	71	0	-1.643	-3.708	-5.689
2.5	71	0	-1.925	-4.327	-6.588
1.5	72	0	-1.505	-3.398	-5.237
2	72	0	-1.764	-3.975	-6.089
2.5	72	0	-2.065	-4.633	-7.039
1.5	73	0	-1.616	-3.644	-5.611
2	73	0	-1.892	-4.258	-6.513
2.5	73	0	-2.215	-4.957	-7.514
1.5	74	0	-1.734	-3.906	-6.007
2	74	0	-2.03	-4.559	-6.961
2.5	74	0	-2.374	-5.299	-8.015
1.5	75	0	-1.859	-4.184	-6.428
2	75	0	-2.176	-4.878	-7.436
2.5	75	0	-2.543	-5.662	-8.543
VSLY			100K	250K	450K

Sys.time()-start # time difference

Time difference of 7.759212 hours

8 Appendix

8.1 GR5 by age

start <- Sys.time()
results_by_age(cost='recessions',initial_recession=5,ages=c(35,45,55,65))

Welfare costs of GR5 by age: full dynamic model
		Realistic mortality
gamma	initial age	exogenous	endogenous	endogenous	endogenous
1.5	35	1.083	1.052	1.006	0.96
2	35	1.49	1.463	1.42	1.376
2.5	35	1.935	1.912	1.873	1.835
1.5	45	1.058	0.979	0.866	0.754
2	45	1.459	1.38	1.265	1.151
2.5	45	1.899	1.824	1.711	1.599
1.5	55	1.026	0.842	0.597	0.352
2	55	1.415	1.215	0.943	0.673
2.5	55	1.837	1.624	1.331	1.04
1.5	65	0	-0.437	-0.997	-1.552
2	65	0	-0.514	-1.178	-1.834
2.5	65	0	-0.606	-1.392	-2.162
VSLY			100K	250K	450K

Sys.time()-start # time difference

Time difference of 50.22257 mins

8.2 Cycles by age

start <- Sys.time()
results_by_age(cost='cycles',initial_recession=0,ages=c(35,45,55,65))

Welfare costs of cycles by age: full dynamic model
		Realistic mortality
gamma	initial age	exogenous	endogenous	endogenous	endogenous
1.5	35	0.557	0.254	-0.18	-0.61
2	35	0.788	0.514	0.107	-0.297
2.5	35	1.053	0.817	0.453	0.093
1.5	45	0.407	-0.009	-0.58	-1.145
2	45	0.589	0.173	-0.408	-0.982
2.5	45	0.801	0.402	-0.17	-0.734
1.5	55	0.275	-0.295	-1.047	-1.791
2	55	0.383	-0.236	-1.063	-1.876
2.5	55	0.503	-0.156	-1.045	-1.915
1.5	65	0	-0.793	-1.802	-2.796
2	65	0	-0.932	-2.122	-3.283
2.5	65	0	-1.096	-2.494	-3.843
VSLY			100K	250K	450K

Sys.time()-start # time difference

Time difference of 50.09129 mins

8.3 Initial income levels

The initial income levels for ages greater than 35 are reported below.

kable(data.frame(
  age = seq(35,75,1),
  y0 = c(parameters$y0_val,
         round(unlist(lapply(ls(pattern="^y0_age."),get)),2)))) %>%
  kable_styling("striped", full_width = F) %>%
  scroll_box(height = "250px")

age	y0
35	1.00
36	1.02
37	1.04
38	1.05
39	1.08
40	1.10
41	1.12
42	1.14
43	1.17
44	1.19
45	1.21
46	1.24
47	1.26
48	1.29
49	1.31
50	1.35
51	1.37
52	1.40
53	1.41
54	1.44
55	1.47
56	1.50
57	1.54
58	1.57
59	1.60
60	1.63
61	1.67
62	1.69
63	1.72
64	1.76
65	1.79
66	1.79
67	1.79
68	1.79
69	1.79
70	1.79
71	1.79
72	1.79
73	1.79
74	1.79
75	1.79

Extension of Krebs (2007)

New: higher displacement probability during Great Recession

Angelo Kisil Marino

June 02, 2023