Scenario analysis via custom shocks

In this tutorial we will illustrate how to perform a scenario analysis by running the model multiple times under a specific shock and comparing the results with the unshocked model.

import BeforeIT as Bit
import StatsBase: mean, std
using Plots

parameters = Bit.AUSTRIA2010Q1.parameters;
initial_conditions = Bit.AUSTRIA2010Q1.initial_conditions;

Initialise the model

model = Bit.Model(parameters, initial_conditions);

Simulate the baseline model for T quarters, N_reps times, and collect the data

T = 16
n_sims = 64
model_vec_baseline = Bit.ensemblerun!((deepcopy(model) for _ in 1:n_sims), T);

Now, apply a shock to the model and simulate it again. A shock is simply a function that takes the model and changes some of its parameters for a specific time period. We do this by first defining a "struct" with useful attributes. For example, we can define an productivity and a consumption shock with the following structs

struct ProductivityShock
    productivity_multiplier::Float64    # productivity multiplier
end

struct ConsumptionShock
    consumption_multiplier::Float64    # productivity multiplier
    final_time::Int
end

and then by making the structs callable functions that change the parameters of the model, this is done in Julia using the syntax below

A permanent change in the labour productivities by the factor s.productivity_multiplier

function (s::ProductivityShock)(model::Bit.Model)
    return model.firms.alpha_bar_i .= model.firms.alpha_bar_i .* s.productivity_multiplier
end

A temporary change in the propensity to consume model.prop.psi by the factor s.consumption_multiplier

function (s::ConsumptionShock)(model::Bit.Model)
    return if model.agg.t == 1
        model.prop.psi = model.prop.psi * s.consumption_multiplier
    elseif model.agg.t == s.final_time
        model.prop.psi = model.prop.psi / s.consumption_multiplier
    end
end

Define specific shocks, for example a 2% increase in productivity

productivity_shock = ProductivityShock(1.02)

Main.ProductivityShock(1.02)

or a 4 quarters long 2% increase in consumption

consumption_shock = ConsumptionShock(1.02, 4)

Main.ConsumptionShock(1.02, 4)

Simulate the model with the shock

model_vec_shocked = Bit.ensemblerun!((deepcopy(model) for _ in 1:n_sims), T; shock! = consumption_shock);

extract the data vectors from the model vectors

data_vector_baseline = Bit.DataVector(model_vec_baseline);
data_vector_shocked = Bit.DataVector(model_vec_shocked);

Compute mean and standard error of GDP for the baseline and shocked simulations

mean_gdp_baseline = mean(data_vector_baseline.real_gdp, dims = 2)
mean_gdp_shocked = mean(data_vector_shocked.real_gdp, dims = 2)
sem_gdp_baseline = std(data_vector_baseline.real_gdp, dims = 2) / sqrt(n_sims)
sem_gdp_shocked = std(data_vector_shocked.real_gdp, dims = 2) / sqrt(n_sims)

17×1 Matrix{Float64}:
   1.8333689901882087e-12
  54.80482761698779
  78.58555684765781
  88.9229421284792
  92.97893122561804
 116.50946684591669
 111.56809062010664
 127.56909078143217
 147.70617991541206
 161.5836832868594
 176.93561908212243
 188.4030894155531
 207.51353617212425
 213.8656455301067
 224.73370912410866
 228.4591691109302
 228.10520088260898

Compute the ratio of shocked to baseline GDP

gdp_ratio = mean_gdp_shocked ./ mean_gdp_baseline

17×1 Matrix{Float64}:
 1.0
 0.9999809901007637
 1.0041619974480396
 1.0071734576833524
 1.0094240585599044
 1.0072344623035758
 1.0063133016925605
 1.0043002784526611
 1.0039838798750387
 1.0029320451054775
 1.003886021632547
 1.0036552799557728
 1.0023135786037045
 1.0023629635521634
 1.0030409364300668
 1.0050577765244546
 1.005599633262487

the standard error on a ratio of two variables is computed with the error propagation formula

sem_gdp_ratio = gdp_ratio .* ((sem_gdp_baseline ./ mean_gdp_baseline) .^ 2 .+ (sem_gdp_shocked ./ mean_gdp_shocked) .^ 2) .^ 0.5

17×1 Matrix{Float64}:
 3.580093467121083e-17
 0.0010637355208865994
 0.0014575669828150056
 0.0017729068567932623
 0.0019328453000787714
 0.0022623033265474646
 0.0025024457237129423
 0.0027526134994003185
 0.0030454732722509964
 0.003349317791252117
 0.003599619886486856
 0.0038848691489997136
 0.004222853197850258
 0.004273288168786942
 0.0043981519658269665
 0.004505900493961222
 0.004602322877531018

Finally, we can plot the impulse response curve

plot(
    1:(T + 1),
    gdp_ratio,
    ribbon = sem_gdp_ratio,
    fillalpha = 0.2,
    label = "",
    xlabel = "quarters",
    ylabel = "GDP change",
)

We can save the figure using: savefig("gdp_shock.png")