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_reps = 64
model_vec_baseline = Bit.ensemblerun(model, T, N_reps);

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)
    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)
    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(model, T, N_reps; 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_reps)
sem_gdp_shocked = std(data_vector_shocked.real_gdp, dims = 2) / sqrt(N_reps)

17×1 Matrix{Float64}:
   0.0
   1.7865068761941965
  21.110315985816616
  41.35949422351184
  55.40933271207089
  69.58970582798813
  86.87652447949039
 101.73678371578859
 111.77651010487111
 121.04474103015707
 129.66599186680705
 137.36395181230398
 146.50821156098633
 150.85098128866176
 162.43347921103492
 173.8016335454425
 178.78416339724694

Compute the ratio of shocked to baseline GDP

gdp_ratio = mean_gdp_shocked ./ mean_gdp_baseline

17×1 Matrix{Float64}:
 1.0
 1.0005003541765403
 1.0059634825383448
 1.008246011507844
 1.0094385770019694
 1.0061080519320071
 1.004501257752781
 1.0031727208012866
 1.00208183324786
 1.0008142987503215
 1.0002291265272192
 0.9999680473725676
 1.000396348679974
 1.0009213519792797
 1.0011619388125652
 1.0012566040096564
 1.0011844442201825

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}:
 0.0
 2.7685467795758086e-5
 0.0004019072192263169
 0.0007586222394135282
 0.0010279555991530576
 0.0013259485352288282
 0.0016244141522447002
 0.0018861392656532442
 0.0020862127251288286
 0.002269706777348419
 0.002408084065851926
 0.0025976670729668176
 0.0028071312896316585
 0.002896757987779719
 0.003026575355825149
 0.0031874822535389634
 0.00328348087495681

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