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. The previous report reproduces the results by calibrating the model using closed-formula solutions (1). This report reproduces the results by adapting the model for a simulation procedure.

First, I outline the theoretical model and its version for the simulation. Second, I specify the steps followed in the simulation. Third, I define the functions to be used in the simulation. Finally, I perform the simulation and compare its results to the calibration results from Krebs (2007 AER).

1 Model

1.1 Income path and random shocks \(\theta\), \(\eta\)

The labor income law of motion is given by Equation 1 from the paper – copied below –, where random shocks are represented by \(\theta_{i,t+1}\) and \(\eta_{i,t+1}\).

\[ y_{i,t+1} = (1+g)(1+\theta_{i,t+1})(1+\eta_{i,t+1})y_{it} \]

The sequence of random variables \(\{\theta_{it}\}_t\) is i.i.d., with log-normal distribution function \(\log{(1+\theta_{i,t+1})} \sim N(-\frac{\sigma^2}{2}, \sigma^2)\) (Equation 2). On the other hand, \(\{\eta_{it}\}_t\) is independently, but not identically distributed; \(\eta_{i,t+1}\) depends on realization of state \(S_{t+1} = S\) as follows (Equation 3)

\[ \eta_{i,t+1} = \begin{cases} -d_S & prob = p_S \\ \frac{p_Sd_S}{1-p_S} & prob = 1-p_S \end{cases} . \] In the equation above, \(S\in\{L,H\}\), and \(P(S=H)=\pi_H\), \(p_S\) is a job displacement probability, and \(d_S\) is the income amount one loses when displaced. These parameters will be further discussed below.

Importantly, agents are hand-to-mouth, i.e. there are no savings and \(y_{it}=c_{it}\).

1.2 Definition of \(\Delta\)

The welfare cost of business cycles is defined as the compensation value the agents need to receive to be, in expectation, indifferent between an economy with, and without business cycles, where \(c_{it}=\overline{c}_{it}\). Formally, the cost of business cycles is given by the \(\Delta\) that satisfies Equation 14,

\[ \mathbb{E}_0\left[ \sum_{t=0}^{\infty}{\beta^t u((1+\Delta)c_{it})} \right] = \mathbb{E}_0\left[ \sum_{t=0}^{\infty}{\beta^t u(\overline{c}_{it})} \right] \]

In this simulation I use a variation Equation 14. I have an economy with \(T\) periods, \(t\in\mathbb{Z}\cap[0,T]\), and \(N\) individuals, \(i\in\mathbb{N}\cap[1,N]\). I present two approaches to find the value of \(\Delta\) that satisfies the equation below.

\[ \sum_{i=1}^{N}{\left[\sum_{t=0}^{T}{\beta^t u((1+\Delta)c_{it})}\right]} = \sum_{i=1}^{N}{\left[\sum_{t=0}^{T}{\beta^t u(\overline{c}_{it})}\right]} \]

2 Simulation steps

I organize the simulation in three main parts, described below.

A. Shocks and income path

First, I draw random shocks \(\theta\) and \(\eta\) and generate income paths for \(T\) periods and \(N\) individuals; following the steps below.

  1. Generate a sequence of random shocks \(\{\theta_{it}\}_t\), drawing from its log-normal distribution function \(\log{(1+\theta_{i,t+1})} \sim N(-\frac{\sigma^2}{2}, \sigma^2)\);
  2. Generate a sequence of random shocks \(\{\eta_{it}\}_t\), according to Equation 3;
  3. Generate income/consumption paths using these random shocks and Equation 1.

I perform steps 1-3 and save it on a dataset with \(N\times T\) rows, and the following variables: agent ID, \(t\), \(\theta\), \(\eta\), \(\overline{\eta}\), \(y\), and \(\overline{y}\).

B. Utility and squared difference

Second, I define a utility function and an aggregate utility calculator, which allows me to calculate the aggregate utility in the scenarios with, and without business cycles, given by \(\overline{U}(\overline{y})\) and \(U(y,\Delta)\), respectively. I also define a squared difference function comparing these two cases: \(g(\Delta) = \left[U(\Delta)-\overline{U}\right]^2\), given \(y\) and \(\overline{y}\). I follow the steps below.

  1. Define utility function \(u(c)\), which equals \(\log{(c)}\) if \(\gamma=1\), and \(c^{1-\gamma}/({1-\gamma})\) otherwise;
  2. Create aggregate utility calculator: \(U=\sum_{i=1}^{N}{\left[\sum_{t=0}^{T}{\beta^t u(c_{it})}\right]}\);
  3. Define squared error function \(g(\Delta) = \left[U(\Delta)-\overline{U}\right]^2\).

C. Compensation value

Lastly, I present two approaches to find \(\Delta\): a grid search, and an optimization. For the grid search, I look for a set of values of \(\Delta\) and choose the one with the smallest \(g(\cdot)\). For the optimization, I minimize \(g(\cdot)\) over \(\Delta\). The steps are outlined below.

  1. Perform grid search for \(\Delta\)
    1. Define grid: finite set of values between 0 and 10%, \(X=\left\{x|x\in[0, .1]\right\}\);
    2. Choose \(\Delta^*=x^*\) such that \(g(x^*)\le g(x)\,\forall\,x\in X\);
  2. Find \(\Delta\) by optimization, i.e. \(\Delta^*=\underset{\Delta}{\arg\min}\, g(\Delta)\).

The simulator function can perform the analysis with multiple draws and calculate the average \(\Delta\) in the end.

3 Functions

In all functions I refer to parameters as values from a parameters’ table that is defined in the fourth section of this report.

3.1 Shocks and income path

3.1.1 Noncyclical component of labor income risk, \(\theta\)

I find \(\lambda = \log{(1+\theta_{i,t+1})}\) by drawing from \(N(-\frac{\sigma^2}{2}, \sigma^2)\) with \(\sigma^2 = 0.01\), and then calculate \(\theta = e^\lambda - 1\). Note that I do not have to draw a specific \(\overline{\theta}\) for the economy without business cycles, because this is a noncyclical component of labor income risk; formally \(\theta_i = \overline{\theta}_i\).

# Noncyclical component of labor income risk, theta
draw_theta <- function(periods=periods_val, s2=parameters$s2_val) {
  lambda <- rnorm(periods,mean=-s2/2,sd=sqrt(s2))
  theta <- exp(lambda) - 1
  return(theta)
}

3.1.2 Cyclical component of labor income risk, \(\eta\)

To draw \(\eta_{it}\) I follow Equation 3, copied below. Recall that \(\pi_S\) is the probability of state \(S\in\{L,H\}\), \(p_S\) is the probability of displacement in state \(S\), and \(d_S\) is an income shock. Intuitively, \(\pi_L=1-\pi_H\), \(p_L\ge p_H\), and \(d_L\ge d_H\).

Within this framework, we have four possible values for \(\eta\), which depends on the combination of the business cycle and the displacement outcomes. Explicitly, the worker can be displaced or not, on a high or low state of the economy. To replicate this in the simulation, I am draw the state of the economy \(S \in \{L,H\}\) with probability \(\pi_L\), the displacement outcome \(displaced \in \{0,1\}\) with probabilities \(p_L,p_H\), and, finally, I use \(d_S\) to find values for \(\eta\), following Equation 3.

\[ \eta_{i,t+1} = \begin{cases} -d_S & prob = p_S \\ \frac{p_Sd_S}{1-p_S} & prob = 1-p_S \end{cases} \]

# Cyclical component of labor income risk, eta
draw_eta <- function(periods=periods_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) {
  
  # 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)
}

Before calculating the income and consumption paths, one also needs to draw the sequence of shocks in the absence of business cycles, denoted by \(\overline{\eta}_i\). Note that in this case, although the state of the economy does not vary, there is still a chance, now constant, of worker displacement in each period. Formally,

\[ \overline{\eta}_i = \begin{cases} -\overline{d} & prob = \overline{p} \\ \frac{\overline{p}\overline{d}}{1-\overline{p}} & prob = (1-\overline{p}) \end{cases} \]

Krebs (2007 AER) presents two methods to calculate \(\overline{p}\) and \(\overline{d}\). Method 1 is used in column 1 of Table 1, and is defined by Equation 13, copied below.

\[\overline{p} = \sum_S{\pi_Sp_S} \;\; and \;\; \overline{d} = \sum_S{\frac{\pi_Sp_S}{\overline{p}}d_S}\]

Method 2 is used in column 2 of Table 1 and is defined by Equation 12, copied below. The function draw_eta_bar can use both methods to draw \(\overline{\eta}\), but I will proceed using only method 1, to simplify the example.

\[\overline{p} = \sum_S{\pi_Sp_S} \;\; and \;\; \overline{d} = \sum_S{\pi_Sd_S}\]

Finally, to replicate Table 2, \(\overline{d}=d_H\). Intuitively, this estimates the costs of recessions, assuming that macroeconomic stabilization policy can remove recessions without affecting economic boom times (Krebs, 2007 AER).

# Eta for case without cycles
draw_eta_bar <- function(periods=periods_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, method=parameters$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.3 Income and consumption path, \(y\) and \(c\)

The income path is given by Equation 1, copied below. For the first period income, I use \(y_0=1\), however, this is an argument of the function and can be easily changed.

\[y_{i,t+1} = (1+g)(1+\theta_{i,t+1})(1+\eta_{i,t+1})y_{it}\]

Remember that the individuals consume everything they earn, so \(y_{it}=c_{it}\). The income path equals the consumption path, which will be used to calculate the utilities.

One should also define the income and consumption paths in the absence of business cycles, i.e. using the shocks \(\overline{\theta}_i\), and \(\overline{\eta}_i\) instead of \(\theta_i\), and \(\eta_i\). This is performed by the function income_path, defined below.

I use this function inside the later defined shocks_income_dataset function, so they will make more sense if looked together. For example, below I refer to a shocks_dataset, that is only defined in the shocks_income_dataset function before I apply the income_path function.

# Income path
income_path <- function(g=parameters$g_val, y0=parameters$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)
}

Create shocks_income dataset

Here, I use the functions defined above to create a dataset with shocks and income paths.

# Create dataset with shocks and income
shocks_income_dataset <- function(N=parameters$N_val,
                                  periods=parameters$periods_val) {
  
  # Shocks dataset
  id <- sort(rep(1:N,periods)) # agent_id
  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,
               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)
}

3.2 Utility and squared difference

3.2.1 Utility function

Below, I define an individual-level utility calculator. For that, I use the functional form of the utility function, given after Equation 6 in the paper. 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}& \gamma \ne 1 \\ \log{(c)}& \gamma = 1 \end{cases} \]

# Utility function as in Krebs
u <- function(gamma,nocycles,Delta) {
  y <- if(nocycles==TRUE){shocks_income$y_bar} else{(1+Delta)*shocks_income$y}
  u_vector <- if(gamma==1){log(y)} else{y^(1-gamma)/(1-gamma)}
  return(u_vector)
}

3.2.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\), and sum it across periods and individuals. Formally, \(U=\sum_{i=1}^{N}{\left[\sum_{t=0}^{T}{\beta^t u(c_{it})}\right]}\), where \(c_{it}\) follows the same definition from above, for the cases with and without business cycles.

# U function
U <- function(gamma, nocycles, Delta, mortality=mortality_val, pi_l=parameters$pi_l_val,
              b=parameters$b_val, periods=parameters$periods_val, N=parameters$N_val) {
  
  # Mortality options
  if (mortality == "none") {
    q <- 1
  }
  if (mortality == "constant") {
    q <- parameters$q_val
  }
  if (!(mortality %in% c("none","constant"))) {
    warning("Mortality must be \'none\' or \'constant\'")
  }
  
  # Vector of utility per period
  per_period_utility <- u(gamma=gamma, nocycles=nocycles, Delta=Delta) 
  
  # Vector of discount rates
  discounting <- rep((b*q)^(0:(periods-1)), N) 
  
  # Aggregate utility
  agg_util <- sum(discounting*per_period_utility) # multiply vectors and sum
  
  return(agg_util) 
}

3.2.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, b=parameters$b_val,
                         periods=parameters$periods_val, N=parameters$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)
}

3.3 Compensation value

3.3.2 Optimization

Lastly, I define the optimization method to find \(\Delta\) such that \(\Delta^*=\underset{\Delta}{\arg\min}\, g(\Delta)\).

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

Calibration results

The calibration results come from Krebs (2007 AER) and the replication report, for \(\gamma\notin\{1,2,3,4\}\).

For the results of business cycles costs and no mortality, we consider column 1 of Table 1 from the paper. For recessions costs and no mortality, column 1 of Table 2. For cycles costs and constant mortality we consider the results from Section V, subsection D.

Deltas_calibration <- function(cost=cost_val,mortality=mortality_val) {
  if (cost=="cycles" & mortality=="none") {
    # 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=="recessions") {
    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=="cycles" & mortality=="constant") {
    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)))))))
    }}
  return(Deltas_cal)
}

Plot results

The function below combine the results from the grid search and the optimization and plot it against the calibration results from Krebs (2007 AER) for each value of \(\gamma\).

# Function to plot results
plot_results <- function(results_table, cost=cost_val, mortality=mortality_val) {
  
  # Dataset for plot
  df_plot <- results_table %>%
    select(c("gamma","calibration","grid_search","optimization")) %>% # order
    pivot_longer(cols = `calibration`:`optimization`, # reshape
                 names_to = "method",
                 values_to = "Delta") %>%
    mutate(Delta = round(Delta,3))
  
  # Plot
  plot <- ggplot(data=df_plot,
                 aes(x=gamma, y=Delta, color=method, shape=method)) +
    geom_point(size=2) + # scatter plot
    ggtitle(paste0("Cost: ", cost, ". Mortality: ", mortality,
                   ". N = ", parameters$N_val, ", T = ", parameters$T_val)) + # title
    scale_x_continuous(breaks=gammas_val) + # x labels
    theme_classic() + 
    theme(plot.title = element_text(size=10, hjust=.5), # title format
          legend.title.align=0.5,
          legend.text = element_text(size=11)) +
    expand_limits(x=4.15)
  
  return(plot)
}

4 Simulation and results

4.1 Parameters

All the parameters used in the simulation are set below, if one wants to change anything, there is no need to touch the rest of the code.

parameters <- data.frame(
  T_val = 100, # last period
  N_val = 1000, # agents in the economy
  grid_values_val = 200, # values in the grid search
  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_val = 0.05, # displacement probability in low state
  p_h_val = 0.03, # displacement probability in high state
  method_val = 1, # method to calculate p and d w/o cycles
  b_val = 0.96, # beta
  g_val = 0.02, # growth rate
  y0_val = 1, # initial income
  neighborhood_val = 1.5, # interval around calibration Delta
  grid_ub = .1, # grid upper bound
  q_val = 0.976) %>% # constant survival probability
  mutate(periods_val = T_val + 1) # because we start in t=0
gammas_val <- c(seq(1,4,by=.5)) # constant degree of relative risk aversion
rounds_val <- 150 # number of draws for multi_draw simulation

4.2 Simulator

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.

simulate <- function(seed, parameters, cost, plot=T, mortality="none") {
  
  # Set parameters equal to chosen ones
  T_val <<- parameters$T_val
  periods_val <<- parameters$periods_val
  N_val <<- parameters$N_val
  grid_values_val <<- parameters$grid_values_val
  neighborhood_val <<- parameters$neighborhood_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 <<- parameters$p_l_val
  p_h_val <<- parameters$p_h_val
  method_val <<- parameters$method_val
  b_val <<- parameters$b_val
  g_val <<- parameters$g_val
  y0_val <<- parameters$y0_val
  mortality_val <<- mortality
  
  # Table argument to be used in other functions
  cost_val <<- cost
  
  # Calibration Deltas
  Deltas_cal <<- Deltas_calibration()
  
  # Set seed
  seed_val <<- seed
  set.seed(seed_val)
  
  # Economy state S
  low <<- rbernoulli(periods_val,pi_l_val) # used to draw_eta and cyclical mortality
  
  # Dataset with shocks and income path
  shocks_income <<- shocks_income_dataset()
  
  # Utility w/o business cycles
  U_bar <<-
    foreach(g=gammas_val, .combine=rbind,
            .export=c("U","u","shocks_income","parameters","mortality_val","low")) %dopar% {
    data.frame(gamma=g, U_bar=U(gamma=g, nocycles=TRUE))
    }
  
  # Plot results or output table
  table_output <- left_join(grid_Deltas(),optimize_Deltas(), by=c("gamma","calibration"))
    
  output <- if (plot==T) plot_results(table_output) else table_output
  
  # Remove temporary objects
  rm(list = c("shocks_income","U_bar","cost_val","Deltas_cal","seed_val","low"), envir=.GlobalEnv)
   
  return(output)
}

If we want to perform a simulation with multiple draws, and then calculate the average cost of business cycles across these simulations, we can run the simulator function multiple times, but with different seeds, and average the results. In short, this is what the function below does.

simulate_multi <- function(rounds, cost, plot=T, mortality="none") {
  
  cost_val <- cost
  mortality_val <- mortality
  
  # 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_val,
                         plot=F, mortality=mortality_val))
   }
  # Calculate averages
  results_multidraw <- appended_results %>% group_by(gamma,calibration) %>%
    summarise(grid_search = mean(grid_search), 
              optimization = mean(optimization))
  
  # Output
  output <- if (plot==T) plot_results(results_multidraw,cost=cost_val) else
    results_multidraw
  
  return(output)
}

5 Results

5.1 One draw

Cost of cycles, no mortality

start <- Sys.time()
simulate_multi(rounds=1, cost="cycles")

Sys.time()-start # time difference
## Time difference of 6.895994 secs

Cost of recessions, no mortality

start <- Sys.time()
simulate_multi(rounds=1, cost="recessions")

Sys.time()-start # time difference
## Time difference of 10.07373 secs

Cost of cycles, constant mortality

start <- Sys.time()
simulate_multi(rounds=1, cost="cycles", mortality="constant")

Sys.time()-start # time difference
## Time difference of 9.75564 secs

5.2 Multiple draws

Cost of recessions, no mortality

start <- Sys.time()
simulate_multi(rounds=rounds_val, cost="cycles")

Sys.time()-start # time difference
## Time difference of 15.3024 mins

Cost of recessions, no mortality

start <- Sys.time()
simulate_multi(rounds=rounds_val, cost="recessions")

Sys.time()-start # time difference
## Time difference of 15.18078 mins

Cost of cycles, constant mortality

start <- Sys.time()
simulate_multi(rounds=rounds_val, cost="cycles", mortality="constant")

Sys.time()-start # time difference
## Time difference of 14.89904 mins

6 Examples

6.0.1 \(\theta\)

draw_theta() %>% plot(ylab = "theta", xlab = "t") # example

6.0.2 \(\eta\) and \(\overline{\eta}\)

low <- rbernoulli(parameters$periods_val, parameters$periods_val)
draw_eta() %>% plot(ylab = "eta", xlab = "t") # example

draw_eta_bar(cost="cycles") %>% plot(ylab = "eta_bar", xlab = "t") # example

6.0.3 \(y\) and \(\overline{y}\)

cost_val="cycles"
shocks_income <- shocks_income_dataset(N=500,periods=101)

# example
ggplot(data=
         shocks_income %>% group_by(t) %>% 
         summarise(y_mean = mean(y), y_bar_mean = mean(y_bar))) +
  geom_line(aes(x=t, y=y_mean, color="y")) +
  geom_line(aes(x=t, y=y_bar_mean, color="y_bar")) +
  scale_color_manual(name = "", values = c("y" = "magenta4", "y_bar" = "deeppink")) +
  theme_classic() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom")

# Print sample dataset
kable(head(shocks_income,n=300L)) %>%
  kable_styling("striped", full_width = F) %>%
  scroll_box(height = "150px")
id t theta eta eta_bar y y_bar
1 0 0.0404833 0.0110526 0.00625 1.0000000 1.0000000
1 1 0.0043255 0.0110526 0.00625 1.0357345 1.0308146
1 2 0.0313845 0.0110526 0.00625 1.1016483 1.0912072
1 3 -0.1653728 0.0110526 0.00625 0.9482207 0.9347723
1 4 0.1502806 0.0110526 0.00625 1.1248308 1.1036102
1 5 -0.0281153 0.0110526 0.00625 1.1273943 1.1008712
1 6 -0.1212987 0.0110526 0.00625 1.0216239 0.9928504
1 7 -0.0627162 0.0110526 0.00625 0.9874977 0.9551268
1 8 -0.0436876 0.0110526 0.00625 0.9738898 0.9374904
1 9 0.0200335 0.0110526 0.00625 1.0244675 0.9814934
1 10 0.1322345 0.0110526 0.00625 1.1962129 1.1405906
1 11 -0.0263754 0.0110526 0.00625 1.2010857 1.1397967
1 12 -0.0947626 0.0110526 0.00625 1.1212705 1.0590000
1 13 -0.0742639 0.0110526 -0.15000 1.0704627 0.8499674
1 14 -0.1395658 0.0110526 0.00625 0.9498677 0.7506301
1 15 0.0148888 0.0110526 0.00625 0.9941582 0.7818987
1 16 -0.0183512 0.0110526 -0.15000 1.0064347 0.6654658
1 17 0.1135956 0.0110526 -0.15000 1.1558116 0.6424988
1 18 0.0152615 0.0110526 0.00625 1.2101491 0.6695088
1 19 0.0189748 0.0110526 0.00625 1.2716754 0.7002060
1 20 -0.0767179 0.0110526 0.00625 1.2108340 0.6635388
1 21 0.0162708 0.0110526 0.00625 1.2690186 0.6921207
1 22 -0.0744374 0.0110526 -0.15000 1.2112888 0.5554010
1 23 0.1484835 0.0110526 0.00625 1.4346515 0.6546927
1 24 -0.0105324 0.0110526 0.00625 1.4639354 0.6648829
1 25 -0.2338105 0.0110526 0.00625 1.1567300 0.5228624
1 26 -0.0482287 0.0110526 0.00625 1.1353729 0.5107708
1 27 0.0130190 0.0110526 0.00625 1.1861239 0.5310675
1 28 0.1088603 -0.2100000 0.00625 1.0598251 0.6044113
1 29 0.0698611 0.0110526 0.00625 1.1693258 0.6636912
1 30 0.0625353 0.0110526 0.00625 1.2813060 0.7237949
1 31 -0.0996791 0.0110526 0.00625 1.1896634 0.6688349
1 32 0.0011579 0.0110526 0.00625 1.2282891 0.6872703
1 33 -0.1082380 0.0110526 0.00625 1.1295970 0.6290463
1 34 -0.0318543 0.0110526 0.00625 1.1278158 0.6250710
1 35 -0.0612762 0.0110526 0.00625 1.0918172 0.6022451
1 36 -0.0400377 0.0110526 0.00625 1.0808815 0.5933809
1 37 0.0536944 0.0110526 0.00625 1.1745369 0.6417329
1 38 0.1000276 0.0110526 0.00625 1.3324294 0.7245426
1 39 0.3347880 0.0110526 0.00625 1.8341314 0.9926183
1 40 0.1418545 0.0110526 0.00625 2.1598080 1.1633198
1 41 -0.1193602 0.0110526 0.00625 1.9614957 1.0514860
1 42 -0.1131555 0.0110526 0.00625 1.7939435 0.9570993
1 43 -0.0664352 0.0110526 0.00625 1.7271385 0.9170806
1 44 0.0340403 -0.2100000 0.00625 1.4391031 0.9733098
1 45 0.2396791 0.0110526 0.00625 1.8398191 1.2384156
1 46 -0.0092929 0.0110526 0.00625 1.8797251 1.2592668
1 47 -0.0409856 0.0110526 0.00625 1.8590600 1.2395069
1 48 0.1101597 -0.2100000 0.00625 1.6630532 1.4123440
1 49 0.0020342 0.0110526 0.00625 1.7185518 1.4525434
1 50 -0.0104371 0.0110526 0.00625 1.7537996 1.4752941
1 51 0.0750461 0.0110526 0.00625 1.9443794 1.6278401
1 52 -0.1183029 0.0110526 0.00625 1.7679678 1.4731170
1 53 -0.1063641 0.0110526 0.00625 1.6293295 1.3511511
1 54 0.0164883 0.0110526 0.00625 1.7079896 1.4096535
1 55 -0.0538152 0.0110526 0.00625 1.6666145 1.3689715
1 56 -0.0557282 0.0110526 0.00625 1.6229536 1.3267757
1 57 -0.0658844 0.0110526 0.00625 1.5634379 1.2720500
1 58 0.1128348 -0.2100000 0.00625 1.4019696 1.4529174
1 59 -0.0153040 0.0110526 0.00625 1.4236876 1.4684161
1 60 0.2159037 0.0110526 0.00625 1.7852039 1.8325440
1 61 -0.0757923 0.0110526 0.00625 1.7014976 1.7383213
1 62 -0.0258967 0.0110526 0.00625 1.7092685 1.7379654
1 63 0.0028978 0.0110526 0.00625 1.7678316 1.7889733
1 64 0.0777354 0.0110526 0.00625 1.9648390 1.9788920
1 65 0.0277807 0.0110526 0.00625 2.0825784 2.0875102
1 66 0.0268709 0.0110526 0.00625 2.2054192 2.2001410
1 67 -0.1045025 0.0110526 0.00625 2.0367112 2.0221852
1 68 -0.0789566 0.0110526 0.00625 1.9345657 1.9116444
1 69 0.0065492 0.0110526 0.00625 2.0081328 1.9749140
1 70 -0.0368240 0.0110526 0.00625 1.9946745 1.9523600
1 71 0.0911876 0.0110526 0.00625 2.2446331 2.1865800
1 72 0.1070142 0.0110526 0.00625 2.5625508 2.4844177
1 73 -0.1386600 0.0110526 0.00625 2.2762555 2.1963689
1 74 0.0125407 0.0110526 0.00625 2.3768810 2.2825686
1 75 -0.0949203 0.0110526 0.00625 2.2185448 2.1203948
1 76 0.0448427 0.0110526 0.00625 2.3905238 2.2739124
1 77 0.0060112 0.0110526 0.00625 2.4801035 2.3479162
1 78 0.1828968 0.0110526 0.00625 3.0254544 2.8505951
1 79 0.1344288 0.0110526 0.00625 3.5394989 3.3190884
1 80 0.0709630 0.0110526 0.00625 3.9092206 3.6483740
1 81 -0.0665064 0.0110526 0.00625 3.7633575 3.4955602
1 82 -0.1186092 0.0110526 0.00625 3.4207231 3.1622147
1 83 0.0719702 -0.2100000 0.00625 2.9547987 3.4792060
1 84 -0.0190590 0.0110526 0.00625 2.9891296 3.5029111
1 85 0.0038324 0.0110526 0.00625 3.0944246 3.6090791
1 86 -0.1850718 0.0110526 0.00625 2.6005976 3.0187128
1 87 0.0498880 0.0110526 0.00625 2.8157239 3.2529008
1 88 -0.0967799 -0.2100000 0.00625 2.0493253 3.0155773
1 89 -0.0556159 0.0110526 0.00625 1.9958758 2.9229757
1 90 0.0849490 0.0110526 0.00625 2.2331442 3.2549220
1 91 0.0600466 0.0110526 0.00625 2.4412692 3.5413725
1 92 -0.1208812 0.0110526 0.00625 2.2132840 3.1954000
1 93 -0.0151204 0.0110526 0.00625 2.2479892 3.2300885
1 94 0.0577634 0.0110526 0.00625 2.4522046 3.5067842
1 95 -0.0307751 0.0110526 0.00625 2.4510672 3.4885077
1 96 -0.0012570 0.0110526 0.00625 2.5245437 3.5760162
1 97 0.0609225 0.0110526 0.00625 2.7621070 3.8939397
1 98 -0.1118080 0.0110526 0.00625 2.5300044 3.5497856
1 99 -0.1272243 0.0110526 0.00625 2.2771826 3.1798808
1 100 0.0406252 0.0110526 0.00625 2.4438027 3.3963407
2 0 0.0021334 0.0110526 0.00625 1.0000000 1.0000000
2 1 0.0461184 0.0110526 0.00625 1.0788343 1.0737097
2 2 -0.0205983 0.0110526 0.00625 1.0896563 1.0793289
2 3 0.0128233 -0.2100000 0.00625 0.8893044 1.1220017
2 4 0.1415758 0.0110526 0.00625 1.0469577 1.3146324
2 5 -0.0105650 0.0110526 0.00625 1.0682929 1.3350504
2 6 -0.0544506 0.0110526 0.00625 1.0417140 1.2956508
2 7 -0.1472702 0.0110526 0.00625 0.9160810 1.1339802
2 8 -0.0809858 0.0110526 0.00625 0.8682205 1.0696305
2 9 0.2708349 0.0110526 0.00625 1.1378713 1.3951759
2 10 0.1008262 0.0110526 0.00625 1.2917719 1.5763542
2 11 0.0753446 0.0110526 0.00625 1.4325422 1.7398328
2 12 0.0312338 0.0110526 0.00625 1.5234862 1.8414958
2 13 -0.0347551 0.0110526 0.00625 1.5165264 1.8243759
2 14 -0.0213664 0.0110526 0.00625 1.5305377 1.8324853
2 15 -0.1264216 0.0110526 0.00625 1.3788589 1.6430411
2 16 0.1136964 0.0110526 0.00625 1.5836550 1.8781113
2 17 -0.0234877 0.0110526 0.00625 1.5948220 1.8823704
2 18 -0.0557555 0.0110526 0.00625 1.5529971 1.8242974
2 19 -0.0504489 0.0110526 0.00625 1.5207677 1.7779520
2 20 0.0669370 0.0110526 0.00625 1.6733069 1.9469952
2 21 0.0204722 0.0110526 0.00625 1.7609649 2.0392577
2 22 -0.1097188 0.0110526 0.00625 1.6167835 1.8633970
2 23 -0.0345877 0.0110526 0.00625 1.6096766 1.8463935
2 24 0.1287566 0.0110526 0.00625 1.8737552 2.1390978
2 25 0.0433912 0.0110526 0.00625 2.0162017 2.2907826
2 26 0.0277829 0.0110526 0.00625 2.1370234 2.4165251
2 27 -0.0511171 0.0110526 0.00625 2.0912014 2.3534774
2 28 0.0098066 0.0110526 0.00625 2.1777498 2.4392386
2 29 0.1054008 -0.2100000 0.00625 1.9397914 2.7674522
2 30 -0.0733696 0.0110526 0.00625 1.8536832 2.6320415
2 31 0.1181622 0.0110526 -0.15000 2.1375400 2.5516238
2 32 -0.1664519 0.0110526 0.00625 1.8374641 2.1829982
2 33 -0.1236486 0.0110526 0.00625 1.6606232 1.9635309
2 34 0.0608439 0.0110526 0.00625 1.8167556 2.1379389
2 35 0.1368553 0.0110526 0.00625 2.1299806 2.4946324
2 36 -0.0557878 0.0110526 0.00625 2.0740499 2.4175877
2 37 0.1094223 0.0110526 0.00625 2.3729580 2.7528669
2 38 -0.1303347 0.0110526 0.00625 2.1282180 2.4572165
2 39 0.1410152 0.0110526 -0.15000 2.5042718 2.4308263
2 40 0.1140415 0.0110526 0.00625 2.8771119 2.7794659
2 41 0.0579105 0.0110526 0.00625 3.1389155 3.0179799
2 42 -0.1135647 0.0110526 0.00625 2.8694628 2.7458034
2 43 0.0386613 0.0110526 0.00625 3.0736079 2.9271801
2 44 0.0970415 0.0110526 0.00625 3.4773262 3.2959343
2 45 0.0231226 0.0110526 0.00625 3.6689946 3.4610853
2 46 -0.0628593 0.0110526 0.00625 3.5458946 3.3290720
2 47 0.0984162 0.0110526 0.00625 4.0166750 3.7531523
2 48 -0.0078844 0.0110526 0.00625 4.1096316 3.8217698
2 49 -0.0646564 0.0110526 0.00625 3.9641311 3.6689498
2 50 0.0850588 0.0110526 0.00625 4.4358333 4.0860258
2 51 -0.0284873 0.0110526 0.00625 4.4442415 4.0743249
2 52 -0.1313126 0.0110526 0.00625 3.9813935 3.6326641
2 53 0.1430762 0.0110526 0.00625 4.6933637 4.2619316
2 54 0.1138346 0.0110526 0.00625 5.3911181 4.8722913
2 55 -0.0285429 0.0110526 0.00625 5.4010277 4.8580606
2 56 0.2850270 0.0110526 0.00625 7.1575202 6.4073913
2 57 0.0838743 0.0110526 0.00625 8.0004689 7.1279761
2 58 0.0939201 0.0110526 -0.15000 9.0255775 6.7603775
2 59 0.1589067 0.0110526 0.00625 10.7869186 8.0412855
2 60 -0.1521116 0.0110526 -0.15000 9.4321359 5.9113039
2 61 -0.0652862 0.0110526 0.00625 9.0920669 5.6711090
2 62 -0.1220653 0.0110526 0.00625 8.2318753 5.1101810
2 63 -0.2364661 0.0110526 -0.15000 6.4818812 3.3828577
2 64 -0.0555337 0.0110526 0.00625 6.3133736 3.2792632
2 65 -0.1015436 0.0110526 0.00625 5.8496841 3.0239830
2 66 -0.0597564 0.0110526 0.00625 5.6721371 2.9182720
2 67 0.1627502 0.0110526 0.00625 6.8015371 3.4827175
2 68 -0.1220865 0.0110526 0.00625 6.1579015 3.1381670
2 69 0.0013296 0.0110526 0.00625 6.3589253 3.2252186
2 70 0.0327336 0.0110526 0.00625 6.7724522 3.4186412
2 71 0.0220603 0.0110526 0.00625 7.1383267 3.5862133
2 72 0.0394461 -0.2100000 0.00625 5.9789602 3.8259929
2 73 0.1814320 -0.2100000 0.00625 5.6919575 4.6393693
2 74 -0.0214822 0.0110526 0.00625 5.7438660 4.6594401
2 75 0.0075759 0.0110526 0.00625 5.9683736 4.8185632
2 76 0.1367392 0.0110526 0.00625 6.9966601 5.6219175
2 77 0.0071206 0.0110526 0.00625 7.2668503 5.8112831
2 78 -0.0137627 0.0110526 0.00625 7.3909721 5.8824672
2 79 0.0266689 0.0110526 0.00625 7.8253882 6.1986336
2 80 0.0374934 0.0110526 0.00625 8.3726930 6.6006602
2 81 0.0522403 0.0110526 0.00625 9.0856089 7.1286678
2 82 -0.0090893 0.0110526 0.00625 9.2845848 7.2501828
2 83 -0.0844936 0.0110526 0.00625 8.7659264 6.8126553
2 84 -0.1913084 0.0110526 0.00625 7.3106280 5.6546459
2 85 0.2082263 0.0110526 0.00625 9.1091303 7.0122884
2 86 0.1453940 0.0110526 -0.15000 10.7598384 6.9635992
2 87 0.0788971 -0.2100000 0.00625 9.3543380 7.7111629
2 88 0.1302577 0.0110526 0.00625 10.9034636 8.9454754
2 89 -0.1271002 0.0110526 0.00625 9.8152820 8.0144526
2 90 0.2308689 0.0110526 0.00625 12.4591529 10.1249229
2 91 -0.0957474 0.0110526 0.00625 11.6185570 9.3969633
2 92 0.0256738 0.0110526 0.00625 12.2895333 9.8924271
2 93 0.0198499 0.0110526 0.00625 12.9254471 10.3548824
2 94 0.0500073 0.0110526 0.00625 13.9962542 11.1594695
2 95 0.0260268 0.0110526 -0.15000 14.8096381 9.9270758
2 96 0.0761462 0.0110526 0.00625 16.4357553 10.9647488
2 97 0.1415869 0.0110526 0.00625 19.3496257 12.8473547
2 98 -0.0859181 0.0110526 0.00625 18.2402838 12.0532695
2 99 -0.0918266 0.0110526 0.00625 17.0833997 11.2351716
2 100 0.0102722 0.0110526 0.00625 17.7986318 11.6499525
3 0 0.0609486 0.0110526 0.00625 1.0000000 1.0000000
3 1 0.0086967 0.0110526 0.00625 1.0402424 1.0353011
3 2 -0.1438944 0.0110526 0.00625 0.9184083 0.9097039
3 3 0.0507559 0.0110526 0.00625 0.9952028 0.9810880
3 4 -0.0369223 -0.2100000 0.00625 0.7723251 0.9697847
3 5 -0.0339688 0.0110526 0.00625 0.7694231 0.9615515
3 6 0.1639584 0.0110526 0.00625 0.9235845 1.1487250
3 7 0.0019660 0.0110526 0.00625 0.9543409 1.1813405
3 8 0.0128330 -0.2100000 0.00625 0.7788765 1.2280584
3 9 0.0791489 0.0110526 -0.15000 0.8668101 1.1489986
3 10 0.0657714 0.0110526 0.00625 0.9527127 1.2568678
3 11 -0.0836091 0.0110526 0.00625 0.9003609 1.1821605
3 12 0.1530919 0.0110526 0.00625 1.0706672 1.3990925
3 13 -0.0361073 0.0110526 0.00625 1.0642829 1.3841436
3 14 0.0697117 -0.2100000 0.00625 0.9173839 1.5196864
3 15 -0.0089268 0.0110526 0.00625 0.9376284 1.5458443
3 16 0.0590275 0.0110526 0.00625 1.0240282 1.6802699
3 17 0.0877537 0.0110526 0.00625 1.1487260 1.8759260
3 18 0.1264109 0.0110526 0.00625 1.3344037 2.1687956
3 19 0.0096394 0.0110526 0.00625 1.3894006 2.2474549
3 20 0.2150147 0.0110526 0.00625 1.7409365 2.8027126
3 21 -0.0347360 0.0110526 0.00625 1.7330176 2.7767114
3 22 0.1051849 0.0110526 0.00625 1.9752035 3.1497186
3 23 0.0438546 0.0110526 0.00625 2.1263061 3.3745652
3 24 0.0042201 0.0110526 0.00625 2.2020574 3.4781859
3 25 0.0413529 0.0110526 0.00625 2.3648331 3.7175497
3 26 0.2392523 0.0110526 0.00625 3.0222763 4.7284912
3 27 -0.0301084 0.0110526 0.00625 3.0229523 4.7070829
3 28 -0.0540662 0.0110526 0.00625 2.9489403 4.5700259
3 29 -0.0710182 0.0110526 0.00625 2.8251864 4.3574450
3 30 -0.1688141 0.0110526 0.00625 2.4216938 3.7173731
3 31 0.0333455 0.0110526 0.00625 2.5807070 3.9426458
3 32 0.0645403 0.0110526 0.00625 2.8331837 4.3078039
3 33 0.0583156 0.0110526 -0.15000 3.0921737 3.9526671
3 34 0.0729871 0.0110526 0.00625 3.4216243 4.3530214
3 35 -0.1034311 0.0110526 0.00625 3.1636607 4.0057193
3 36 0.1257340 0.0110526 0.00625 3.6728197 4.6283090
3 37 0.0289541 0.0110526 0.00625 3.8973512 4.8879237
3 38 -0.0496439 0.0110526 0.00625 3.8197051 4.7677869
3 39 -0.1203187 0.0110526 -0.15000 3.4652067 3.6363134
3 40 0.0027584 0.0110526 0.00625 3.5834337 3.7425160
3 41 0.1110294 -0.2100000 0.00625 3.2081317 4.2677138
3 42 0.0183515 0.0110526 0.00625 3.3691771 4.4606595
3 43 -0.0853382 0.0110526 0.00625 3.1780323 4.1876046
3 44 0.1238400 0.0110526 0.00625 3.6832970 4.8303235
3 45 0.0354631 0.0110526 0.00625 3.9331935 5.1335398
3 46 -0.0843527 0.0110526 0.00625 3.7140478 4.8244881
3 47 -0.0340745 -0.2100000 0.00625 2.8908023 4.7830061
3 48 -0.0791669 0.0110526 0.00625 2.7451952 4.5205149
3 49 -0.2537232 0.0110526 0.00625 2.1127450 3.4625327
3 50 -0.0576552 0.0110526 0.00625 2.0531981 3.3489587
3 51 0.0835079 0.0110526 0.00625 2.2942296 3.7243282
3 52 -0.0888027 0.0110526 0.00625 2.1558732 3.4831038
3 53 -0.1149822 0.0110526 0.00625 1.9676559 3.1639125
3 54 -0.1934355 0.0110526 0.00625 1.6366741 2.6192060
3 55 -0.1070058 0.0110526 0.00625 1.5072482 2.4006252
3 56 -0.0131647 0.0110526 0.00625 1.5339225 2.4315047
3 57 0.0637892 0.0110526 0.00625 1.6828017 2.6548304
3 58 -0.0965316 0.0110526 0.00625 1.5679054 2.4618173
3 59 0.0624370 0.0110526 0.00625 1.7178964 2.6845102
3 60 0.0677139 0.0110526 0.00625 1.8915847 2.9418873
3 61 -0.1743663 0.0110526 0.00625 1.6105980 2.4929842
3 62 -0.0521291 0.0110526 0.00625 1.5743826 2.4253520
3 63 0.0127768 0.0110526 0.00625 1.6443640 2.5211262
3 64 -0.0582179 0.0110526 0.00625 1.5970640 2.4369750
3 65 -0.0281186 -0.2100000 0.00625 1.2507279 2.4309185
3 66 0.0005330 0.0110526 0.00625 1.2905302 2.4963637
3 67 -0.0157391 0.0110526 0.00625 1.3099429 2.5218786
3 68 -0.1873931 0.0110526 0.00625 1.0977585 2.1033461
3 69 0.1264038 -0.2100000 0.00625 0.9963872 2.4317051
3 70 0.0038114 0.0110526 0.00625 1.0314644 2.5053540
3 71 0.0324231 0.0110526 0.00625 1.0982112 2.6548064
3 72 0.0798162 0.0110526 0.00625 1.2229526 2.9423123
3 73 -0.1597889 0.0110526 0.00625 1.0596734 2.5373669
3 74 0.0699744 0.0110526 0.00625 1.1692822 2.7865236
3 75 0.0300194 0.0110526 0.00625 1.2420488 2.9458741
3 76 -0.1041085 0.0110526 0.00625 1.1475405 2.7087921
3 77 0.1307905 0.0110526 0.00625 1.3382095 3.1438649
3 78 0.1126757 0.0110526 0.00625 1.5355595 3.5903647
3 79 0.0298023 0.0110526 0.00625 1.6307766 3.7948839
3 80 0.0624391 0.0110526 0.00625 1.7867855 4.1381725
3 81 0.0434144 0.0110526 0.00625 1.9226631 4.4317114
3 82 -0.0723415 0.0110526 0.00625 1.8393538 4.2195456
3 83 0.0378904 0.0110526 0.00625 1.9687505 4.4949331
3 84 0.0207931 0.0110526 0.00625 2.0725372 4.7094155
3 85 0.0552823 0.0110526 0.00625 2.2555109 5.1008404
3 86 0.0361545 0.0110526 0.00625 2.4101461 5.4246574
3 87 0.0008703 0.0110526 0.00625 2.4876835 5.5725786
3 88 0.0539160 0.0110526 0.00625 2.7038032 6.0279311
3 89 -0.0022921 -0.2100000 0.00625 2.1737309 6.1727371
3 90 -0.1090240 0.0110526 0.00625 1.9973112 5.6448170
3 91 0.1549421 0.0110526 0.00625 2.3789202 6.6913868
3 92 0.0969131 0.0110526 0.00625 2.6910765 7.5334590
3 93 0.2624434 0.0110526 0.00625 3.5035789 9.7614069
3 94 -0.0591422 0.0110526 0.00625 3.3994592 9.4263264
3 95 -0.0212879 0.0110526 0.00625 3.4311423 9.4689867
3 96 0.0168427 0.0110526 0.00625 3.5980438 9.8824210
3 97 -0.1295871 0.0110526 0.00625 3.2297261 8.8286586 
rm(list=c("shocks_income","low",ls(pattern="*_val$"))) # removes example dataset