I present here a replication of some results from Krebs (2007 AER). More specifically, I reproduce the results of the costs of business cycles and recessions, from Tables 1 and 2, and the results with stochastic constant mortality rate, discussed briefly in Section V-D of the paper. In the end, I present tables with the values that will be used for reference when comparing the simulation with the calibration.
More details on the model can be seen in the paper or in the simulation report, here I only copy and use the closed formula solution derived by Krebs.
I still need to understand some divergences between my and Krebs’ results, although most of them are considerably close. Noteworthy, the most important results, used in other reports, are from column 1, which are correct here.
The equations used to generate the results closed formula calibration results are presented below. This is Equation 16 from Krebs (2007 AER).
If \(\gamma \neq 1\):
\[ \Delta = -1 + \frac{\left\{1 - \beta(1+g)^{1-\gamma} e^{\frac{1}{2}\gamma(\gamma-1)\sigma^2} \times \sum_S \pi_s\left[p_S(1-d_S)^{1-\gamma} + (1-p_S) \left(1 + \frac{p_Sd_S}{1-p_S}\right)^{1-\gamma}\right]\right\}^{\frac{1}{1-\gamma}}} {\left\{1 - \beta(1+g)^ {1-\gamma} e^{\frac{1}{2}\gamma(\gamma-1)\sigma^2} \times \left[\overline{p}(1-\overline{d})^{1-\gamma} + (1-\overline{p}) \left(1+\frac{\overline{p}\overline{d}}{1-\overline{p}}\right)^{1-\gamma}\right]\right\}^{\frac{1}{1-\gamma}}} \]
if \(\gamma = 1\), log-utility case:
\[ \begin{aligned} \Delta = \;&\frac{\beta}{1-\beta} \left[\overline{p}\log{(1-\overline{d})} + (1-\overline{p})\log{\left(1+\frac{\overline{p}\overline{d}}{1-\overline{p}}\right)}\right] \\ &- \frac{\beta}{1-\beta}\sum_S\pi_Sp_S\log{(1-d_S)} \\ &- \frac{\beta}{1-\beta}\sum_S\pi_S(1-p_S)\log{\left(1+\frac{p_Sd_S}{1-p_S}\right)} \hfill . \end{aligned} \]
The author choose the following values for the calibration exercise. Here, they are all stored on the table below, which can be modified.
parameters <- data.frame(
s2_val = 0.01, # sigma squared
s2_log_val = 0.01042, # sigma squared for log normal case
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_val = 0.05, # displacement probability in low state
p_h_val = 0.03, # displacement probability in high state
beta_val = 0.96, # beta
g_val = 0.02, # growth rate
q_val = 0.976) # constant survival probability
gammas_val <- seq(1,4,.5) # constant degree of relative risk aversionThe function below calculates the costs of business cycles, \(\Delta\), given a set of parameters. Here, I use the fact that \(\overline{x} = \sum_S\pi_Sx_S\) with \(x\in\{p,d\}\).
calibrate_Deltas <-
function(method, cost="cycles", mortality="none", gammas=gammas_val, g=parameters$g_val,
beta=parameters$beta_val, pi_l=parameters$pi_l_val,
d_l=parameters$d_l_val, d_h=parameters$d_h_val,
p_l=parameters$p_l_val, p_h=parameters$p_h_val,
s2=parameters$s2_val) {
# Update beta for constant mortality
if (mortality=="constant") {
beta=parameters$beta_val*parameters$q_val
}
# State probs and displacement shocks
pi_h <- 1-pi_l # probability if high state
p_bar <- ifelse(method == 1, # displacement prob if no cycles
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
d_bar <- ifelse(method == 1, # shock if no cycles
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 # if cost of recessions, always in high
}
# Vector of Deltas
Deltas <- NULL
for (i in 1:length(gammas)) {
Deltas[i] <- 100*
ifelse(gammas[i] != 1, # Equation 16 from the paper
# Delta if gamma different than 1
-1
+ ((1-beta*(1+g)^(1-gammas[i])*exp(1/2*gammas[i]*(gammas[i]-1)*s2)
*(pi_l*(p_l*(1-d_l)^(1-gammas[i])+(1-p_l)*(1+p_l*d_l/(1-p_l))^(1-gammas[i]))
+pi_h*(p_h*(1-d_h)^(1-gammas[i])+(1-p_h)*(1+p_h*d_h/(1-p_h))^(1-gammas[i]))))
^(1/(1-gammas[i]))) # numerator ends here
/((1-beta*(1+g)^(1-gammas[i])*exp(1/2*gammas[i]*(gammas[i]-1)*s2)
*(p_bar*(1-d_bar)^(1-gammas[i])+(1-p_bar)*(1+p_bar*d_bar/(1-p_bar))^(1-gammas[i]))))
^(1/(1-gammas[i])), # denominator ends here
# Delta if gamma equals 1
beta/(1-beta)*(p_bar*log(1-d_bar)+(1-p_bar)*log(1+p_bar*d_bar/(1-p_bar))) # first line ends here
- beta/(1-beta)*(pi_l*p_l*log(1-d_l) + pi_h*p_h*log(1-d_h)) # second line ends here
- beta/(1-beta)*(pi_l*(1-p_l)*log(1+p_l*d_l/(1-p_l)) +
pi_h*(1-p_h)*log(1+p_h*d_h/(1-p_h)))) # third line ends here
}
return(Deltas)
}The following function copies Krebs’ results, for reference.
krebs_results <- function(table,gammas=gammas_val) {
# This could be more automated, but I didn't think it was worth it
if (table==1) {
krebs <- data.frame(
gammas = gammas,
kc1 = c(.527,NA,.982,NA,1.787,NA,4.092),
kc2 = c(.247,NA,.479,NA,.908,NA,2.183),
kc3 = c(.207,NA,.477,NA,.895,NA,2.068),
kc4 = c(.003,NA,.005,NA,.008,NA,.014))
}
if (table==2) {
krebs <- data.frame(
gammas = gammas,
kc1 = c(1.351,NA,2.425,NA,4.217,NA,9),
kc2 = c(1.069,NA,1.914,NA,3.304,NA,6.880))
}
if (table=="mortality") {
krebs <- data.frame(
gammas = gammas,
kc1 = c(.33,NA,.65,NA,1.15,NA,2.25))
}
return(krebs)
}The function below generate the formatted results tables.
results <- function(table) {
# Krebs results
krebs <- krebs_results(table=table)
# Our results
if (table==1) {
us <- data.frame(
gammas = gammas_val,
c1 = calibrate_Deltas(method=1),# Column 1 - Method 1
c2 = calibrate_Deltas(method=2),# Column 2 - Method 2
c3 = calibrate_Deltas(method=1,p_l=.04,p_h=.04),# Column 3 - p_l = p_h = .04
c4 = calibrate_Deltas(method=1,s2=parameters$s2_log_val)) # Column 4 - log-normal
}
if (table==2) {
us <- data.frame(
gammas = gammas_val,
c1 = calibrate_Deltas(method=1,cost="recessions"), # Column 1 - Method 1
c2 = calibrate_Deltas(method=1,p_l=.04,p_h=.04,cost="recessions")) # Column 3 - p_l = p_h = .04
}
if (table=="mortality") {
us <- data.frame(
gammas = gammas_val,
c1 = calibrate_Deltas(method=1,mortality="constant")) # changes beta to beta*q
}
# Order columns
if (table==1) {
column_order <- c("gammas","kc1","c1","kc2","c2","kc3","c3","kc4","c4")
header <- c(" "=1,"Column 1"=2,"Column 2"=2,"Column 3"=2,"Column 4"=2)
}
if (table==2) {
column_order <- c("gammas","kc1","c1","kc2","c2")
header <- c(" "=1,"Column 1"=2,"Column 2"=2)
}
if (table=="mortality") {
column_order <- c("gammas","kc1","c1")
header <- c(" "=1,"Column 1"=2)
}
# Merge results
merged_results <- left_join(krebs,us,by="gammas") %>%
select(column_order) %>% round(3)
# Format
if (parameters$beta_val==.96) {
table <-
kable(merged_results,
col.names=c("gamma",rep(c("krebs","us"),
ifelse(table==1,4,ifelse(table==2,2,1)))),
caption=paste0("Table ",table, ", beta ", parameters$beta_val)) %>%
kable_styling("striped", full_width = F) %>%
add_header_above(header)
}
else {
table <-
kable(merged_results %>% select(-contains("k")),
col.names=c("gamma",paste(1:ifelse(table==1,4,ifelse(table==2,2,1)))),
caption=paste0("Table ", table, ", beta ", parameters$beta_val)) %>%
kable_styling("striped", full_width = F) %>%
add_header_above(
c(" "=1, "Reference columns"=ncol(merged_results %>% select(-contains("k")))-1))
}
return(table)
}Table 1 is almost perfectly reproduced. The exceptions are:
I am not sure why the Column 3 estimate, for \(\gamma=1\) is different, I have to further investigate it. Regarding Column 4, I must be misunderstanding something, so I need to re-read the relevant sections of the paper. I am simply replacing \(\sigma^2=0.01\) by \(\sigma_{\log{y}}^2=0.01074\), but, intuitively, such a small change is not enough to change the results as in the paper.
results(table=1)| gamma | krebs | us | krebs | us | krebs | us | krebs | us |
|---|---|---|---|---|---|---|---|---|
| 1.0 | 0.527 | 0.526 | 0.247 | 0.247 | 0.207 | 0.247 | 0.003 | 0.526 |
| 1.5 | NA | 0.735 | NA | 0.352 | NA | 0.351 | NA | 0.738 |
| 2.0 | 0.982 | 0.982 | 0.479 | 0.479 | 0.477 | 0.477 | 0.005 | 0.990 |
| 2.5 | NA | 1.309 | NA | 0.651 | NA | 0.647 | NA | 1.331 |
| 3.0 | 1.787 | 1.787 | 0.908 | 0.908 | 0.895 | 0.895 | 0.008 | 1.839 |
| 3.5 | NA | 2.570 | NA | 1.335 | NA | 1.299 | NA | 2.698 |
| 4.0 | 4.092 | 4.092 | 2.183 | 2.183 | 2.068 | 2.068 | 0.014 | 4.469 |
Table 2 is also almost identical to the one in Krebs (2007). However, I would like to understand what is causing the small difference in the estimates from the first row, i.e. when \(\gamma=1\).
results(table=2)| gamma | krebs | us | krebs | us |
|---|---|---|---|---|
| 1.0 | 1.351 | 1.342 | 1.069 | 1.064 |
| 1.5 | NA | 1.852 | NA | 1.464 |
| 2.0 | 2.425 | 2.425 | 1.914 | 1.914 |
| 2.5 | NA | 3.165 | NA | 2.490 |
| 3.0 | 4.217 | 4.217 | 3.304 | 3.304 |
| 3.5 | NA | 5.888 | NA | 4.576 |
| 4.0 | 9.000 | 9.000 | 6.880 | 6.881 |
In part V.D, Krebs suggests a robustness check where workers have a stochastic constant mortality rate. In other words, workers face a constant probability, \(1-q\), of death in any period. Assuming that the chance of death independent of business cycles, he solves the model the same way as before, but using \(\tilde{\beta} = \beta * q\).
Krebs choose \(q\) such that the implied expected lifespan of any worker equals 40 years, which yields \(q=0.976\).
results(table="mortality")| gamma | krebs | us |
|---|---|---|
| 1.0 | 0.33 | 0.325 |
| 1.5 | NA | 0.476 |
| 2.0 | 0.65 | 0.648 |
| 2.5 | NA | 0.863 |
| 3.0 | 1.15 | 1.153 |
| 3.5 | NA | 1.573 |
| 4.0 | 2.25 | 2.244 |
Below, I report a table that contains the results for the four cases we are interested in to compare simulation and calibration, namely, cost of cycles and recessions without and with constant mortality. These results will be used as reference values for the replication exercises.
reference_table <- function(beta) {
# Calculate results
parameters$beta_val <<- beta # defines beta value
results <- data.frame(
gammas = gammas_val,
cycles_none = calibrate_Deltas(method=1),
cycles_constant = calibrate_Deltas(method=1, mortality="constant"),
rec_none = calibrate_Deltas(method=1, cost="recessions"),
rec_constant = calibrate_Deltas(method=1, cost="recessions", mortality="constant")) %>%
round(3)
# Format table
table_output <-
kable(results,
col.names = c("gammas",rep(c("No mortality","Constant\nmortality"),2)),
caption = paste0("Reference values, beta ",parameters$beta_val)) %>%
kable_styling("striped", full_width = F) %>%
add_header_above(c(" "=1,"Cost of cycles"=2, "Cost of recessions"=2))
return(table_output)
}reference_table(beta=.96)| gammas | No mortality | Constant mortality | No mortality | Constant mortality |
|---|---|---|---|---|
| 1.0 | 0.526 | 0.325 | 1.342 | 0.831 |
| 1.5 | 0.735 | 0.476 | 1.852 | 1.199 |
| 2.0 | 0.982 | 0.648 | 2.425 | 1.602 |
| 2.5 | 1.309 | 0.863 | 3.165 | 2.091 |
| 3.0 | 1.787 | 1.153 | 4.217 | 2.731 |
| 3.5 | 2.570 | 1.573 | 5.888 | 3.636 |
| 4.0 | 4.092 | 2.244 | 9.000 | 5.038 |
reference_table(beta=.99)| gammas | No mortality | Constant mortality | No mortality | Constant mortality |
|---|---|---|---|---|
| 1.0 | 2.168 | 0.627 | 5.536 | 1.601 |
| 1.5 | 2.225 | 0.858 | 5.636 | 2.161 |
| 2.0 | 2.677 | 1.135 | 6.612 | 2.803 |
| 2.5 | 3.593 | 1.514 | 8.621 | 3.658 |
| 3.0 | 5.559 | 2.089 | 12.823 | 4.921 |
| 3.5 | 11.809 | 3.084 | 25.251 | 7.033 |
| 4.0 | NaN | 5.222 | NaN | 11.353 |