Simple RBC Model

You can follow the tutorial by reading this page and copying and pasting code into your Julia REPL session. In this case, you will need the model file, simple_RBC.jl.

Simple RBC Model

Part 1: The model

The simple RBC model

In this tutorial, we will use the simple RBC model presented by Villemot (2013).

In period $t=1,\ldots,\infty$, a representative household consumes $C_t$; it also provides labour $L_t$ and rents capital $K_{t-1}$ to the firm. $K_t$ is the capital stock at the end of period $t$, and is available to be rented out during $t+1$.

Consumption, labour and capital are chosen in order to maximize the sum of discounted expected utility

\[\displaystyle\sum_{t=1}^\infty\beta^{t-1}E_t\left[\log(C_t)-\frac{L_t^{1+\gamma}}{1+\gamma}\right]\]

under the sequence of budget constraints

\[\displaystyle K_t=K_{t-1}(1-\delta)+w_tL_t+r_tK_{t-1}-C_t,\]

where

$\beta=\frac{1}{1+\rho}$ is the discount rate and $\rho \in (0,\infty)$ is the rate of time preference,
$\gamma \in (0,\infty)$ is the labour supply parameter,
$\delta \in (0,1)$ is the rate of depreciation of capital,
$w_t$ is the real wage, and
$r_t$ is the real rental rate.

In period $t$, the firm uses labour $L_t$ and capital $K_{t-1}$ to produce $Y_t$ according to the production function

\[\displaystyle Y_t=A_tK_{t-1}^\alpha((1+g)^tL_t)^{1-\alpha},\]

where

$g \in (0,\infty)$ is the growth rate, and
$\alpha$ is the output elasticity of labour.

The firm chooses labour and capital in order to maximize profits

\[\displaystyle \max_{L_t,K_{t-1}} A_t K_{t-1}^\alpha ((1+g)^t L_t)^{1-\alpha} - r_t K_{t-1} - w_t L_t.\]

$A_t$ is a technological shock that follows the AR(1) process

\[\displaystyle \log(A_t)=\lambda\log(A_{t-1})+e_t,\]

where

$e_t$ is an i.i.d. zero-mean normally distributed error term with standard deviation $\gamma$, and
$\lambda \in (0,1)$ is a parameter governing the persistence of the shock.

The household problem

The Lagrangian of the constrained maximization is

\[\displaystyle \mathcal{L}(C_t,L_t,K_t) = \sum_{t=1}^\infty\beta^{t-1}E_t\left[\log(C_t)-\frac{L_t^{1+\gamma}}{1+\gamma}-\mu_t(K_t-K_{t-1}(1-\delta)-w_tL_t-r_tK_{t-1}+C_t)\right]\]

The first order conditions are

\[\begin{aligned} \frac{\partial\mathcal{L}}{\partial C_t} &= \beta^{t-1} \left(\frac{1}{C_t} - \mu_t \right) = 0 \\ \frac{\partial\mathcal{L}}{\partial L_t} &= \beta^{t-1} \left(L_t^\gamma - \mu_t w_t \right) = 0 \\ \frac{\partial\mathcal{L}}{\partial K_t} &= -\beta^{t-1} \mu_t + \beta^t E_t \left(\mu_{t+1}(1-\delta+r_{t+1}) \right) = 0. \end{aligned}\]

Once we eliminate the Lagrange multiplier $\mu_t$, we get

\[\begin{aligned} L_t^\gamma &= \frac{w_t}{C_t} \\ \frac{1}{C_t} &= \beta E_t \left(\frac{1}{C_{t+1}}(r_{t+1}+1-\delta) \right). \end{aligned}\]

The firm problem

The first order conditions are

\[\begin{aligned} r_t &= \alpha A_t K_{t-1}^{\alpha-1}((1+g)^t L_t)^{1-\alpha} \\ w_t &= (1-\alpha)A_t K_{t-1}^\alpha ((1+g)^t)^{1-\alpha} L_t^{-\alpha}. \end{aligned}\]

The goods market equilibrium

Aggregate demand must equal aggregate supply to clear the goods market.

\[\displaystyle K_t + C_t = K_{t-1}(1-\delta)+A_t K_{t-1}^\alpha ((1+g)^t L_t)^{1-\alpha}\]

The dynamic equilibrium

The dynamic equations are obtained by combining the first order conditions of the household and firm problems with the goods market equilibrium.

Based on the goods market equilibrium, consumption and capital must be growing at the same rate: $g_c=g_k=g$.

Thus, we can define stationary variables as

\[\begin{aligned} \hat{C}_t &= \frac{C_t}{(1+g)^t} \\ \hat{K}_t &= \frac{K_t}{(1+g)^t} \\ \hat{w}_t &= \frac{w_t}{(1+g)^t} \end{aligned}\]

Once stationarized (see Villemot (2013)), the dynamic equations can be written as:

\[\begin{aligned} \frac{1}{\hat{C}_t} &= \frac{1}{1+\rho} E_t \left(\frac{1}{\hat{C}_{t+1}(1+g)}(r_{t+1}+1-\delta) \right) \\ L_t^\gamma &= \frac{\hat{w}_t}{\hat{C}_t} \\ r_t &= \alpha A_t \left( \frac{\hat{K}_{t-1}}{1+g} \right)^{\alpha-1} L_t^{1-\alpha} \\ \hat{w}_t &= (1-\alpha) A_t \left( \frac{\hat{K}_{t-1}}{1+g} \right)^\alpha L_t^{-\alpha} \\ \hat{K}_t + \hat{C}_t &= \frac{\hat{K}_{t-1}}{1+g} (1-\delta) + A_t \left( \frac{\hat{K}_{t-1}}{1+g} \right)^\alpha L_t^{1-\alpha} \end{aligned}\]

The next part will discuss how to implement the simple RBC model in StateSpaceEcon.

Part 2: Implementation of the model in `StateSpaceEcon`

Installing the packages

We start by installing the packages needed for this tutorial.

julia> using StateSpaceEcon, ModelBaseEcon, TimeSeriesEcon, Test, Plots, Random, Distributions
julia> # Fix the random seed for reproducibility.
       Random.seed!(1234);

Writing the model file

In StateSpaceEcon, a model is written in its own dedicated module, which is contained in its own file, simple_RBC.jl.

A docstring can be added to the model to provide more details:

    """
    Simple RBC Model
    Model available at: https://archives.dynare.org/DynareShanghai2013/order1.pdf
    Presentation: Villemot, S., 2013. First order approximation of stochastic models. Shanghai Dynare Workshop.
    """
    module simple_RBC
        ...
    end

Then, the module is created with the same name as the model and the associated file name. The model will be constructed with macros taken from the package ModelBaseEcon. So, we need to load ModelBaseEcon within the module simple_RBC with using ModelBaseEcon. The model itself will be a global variable called model within the module simple_RBC. The command const declares global variables that will not change and the function Model() constructs a new model object.

module simple_RBC
    using ModelBaseEcon
    const model = Model()
    # Write the rest of the model below.
end # module

Flags and Options

The model object has flags and options.

Flags are (usually boolean) values which characterize the type of model we have. For example, we can specify that the model is stationary by setting the flag ssZeroSlope to true.

model.flags.ssZeroSlope = true

Once we call @initialize at the end of the model file, the flags must not be changed after that.

Options are values that adjust the operations of the algorithms. They can be assigned in the model file as well, but they can also be changed at any time after that.

We can preset model options with the function setoption!. Below, we set the desired accuracy with tol and we set maxiter for the maximum number of iterations for the iterative solvers. Auxiliary variables will not be created and substituted to help the solver (substitutions). We will opt for QR factorization, which is slower but more robust than LU factorization. Finally, we will set verbose to true to provide more information from the commands.

setoption!(model) do o
    o.tol = 1e-12
    o.maxiter = 100
    o.substitutions = false
    o.factorization = :default
    o.verbose = true
end # options

Many functions in StateSpaceEcon have optional arguments of the same name as a model option. When the argument is not explicitly given in the function call, these functions will use the value from the model option of the same name.

Model Parameters

The rest of the model file deals with the economics of the model. Different elements of the model are declared with macros, which do not have to be in any particular order.

The model object holds model parameters, which are values that can be used in the model equations. The macro @parameters declares the parameters and assigns their values. A link between parameters can be created with the macro @link. Below, the parameter $\beta$ depends on the parameter $\rho$.

@parameters model begin
    α = 0.33
    δ = 0.1
    ρ = 0.03
    λ = 0.97
    γ = 0.5
    g = 0.015
    β = @link 1/(1+ρ)
end # parameters

We must not declare new parameters after @initialize has been called. However, the values of existing parameters can be changed at any time and the new values take effect immediately (with some exceptions that we discuss below).

Variables

Similarly, model variables are specified with the macro @variables. Variables can be declared one line at a time (as with the parameters previously), or over one line by separating them with semicolons ;.

In this model all variables are strictly positive. We can take advantage of this by declaring the variables with the macro @logvariables. This indicates to the solver to work with the log of the variables. For instance, instead of working with $C_t$, the solver will work directly with the logarithm as a standalone variable ($logC_t=\log(C_t)$). The substitution $C_t = e^{logC_t}$ is done automatically in all equations. This is completely invisible to the user, but it can help the solver to avoid computing the logarithm of a negative value.

@logvariables model begin
    "Consumption" C
    "Capital Stock" K
    "Labour" L
    "Real Wage" w
    "Real Rental Rate" r
    "Technological shock" A
end # variables

Similarly, the macro @shocks declares model shocks.

@shocks model begin
    ea
end # shocks

In this case, the begin … end block is not mandatory and the technology shock can be declared in one line.

@shocks model ea

The macro @autoexogenize links a variable with a shock. This can be useful to back out historical shocks with the command autoexogenize! (see below for an example).

@autoexogenize model begin
    A = ea
end # autoexogenize

Equations

The dynamic equations of the model are defined with the macro @equations. The variables have to be indexed with t. For instance, K[t-1] refers to the capital stock at the end of t-1 and C[t+1] refers to the expectation at t for consumption at t+1, or $E_t(C_{t+1})$. When variables can be separated by a logarithm, it will help the solver to put the macro @log in front of an equation. In this way, the residual will be computed as the difference between the logarithm of the left-hand side and the logarithm of the right-hand side. For instance, the solver will prefer to work with the second equation below, because it is linear.

\[\begin{aligned} L_t^\gamma &= \frac{w_t}{C_t} \\ \gamma\log(L_t) &= \log(w_t) - \log(C_t) \end{aligned}\]

We also take the time to rearrange the equations a bit to improve their numerical stability. As a simple rule-of-thumb, we try to make sure we don't have a model variable in the denominator, or raised to a negative power.

@equations model begin
    @log C[t+1] * (1 + g) = β * C[t] * (r[t+1] + 1 - δ)
    @log (L[t])^γ * C[t] = w[t]
    @log r[t] * (K[t-1] / (1 + g))^(1 - α) = α * A[t] * L[t]^(1 - α)
    @log w[t] * L[t]^(α) = (1 - α) * A[t] * (K[t-1] / (1 + g))^α
    @lin K[t] + C[t] = A[t] * (K[t-1] / (1 + g))^α * (L[t])^(1 - α) + (1 - δ) * (K[t-1] / (1 + g))
    log(A[t]) = λ * log(A[t-1]) + ea[t]
end # equations

Initialization of the Model

Once the parameters, the variables, the shocks and the equations have been specified, the macro @initialize constructs the model within the module simple_RBC.

@initialize(model)

Loading the model

We load the module that contains the model with using simple_RBC; the model itself is a global variable called model within that module, which we assign to m in the Main module.

julia> unique!(push!(LOAD_PATH, realpath(".")));
julia> using simple_RBC
julia> m = simple_RBC.model6 variable(s): C, K, L, w, r, A
1 shock(s): ea
7 parameter(s): g = 0.015, α = 0.33, γ = 0.5, δ = 0.1, λ = 0.97, ρ = 0.03, β = @link 1 / (1 + ρ)
6 equations(s):
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 - α))
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
Maximum lag: 1
Maximum lead: 1

Important note

For the "using simple_RBC" command to work, we need the model file simple_RBC.jl to be on the search path for modules. We can do this by:

Running the file that contains the model module with include("/[path to file]/simple_RBC.jl") (in this case we don't need using simple_RBC);
Adding the file path to LOAD_PATH global variable (which we did above);
Putting the model in its own standalone package and adding it to the current Julia environment with using Pkg; Pkg.add("/[path to the package]/simple_RBC").

Examining the model

If the model has more than 20 equations the display gets truncated. We can see the entire model with fullprint.

julia> fullprint(m)6 variable(s): C, K, L, w, r, A
1 shock(s): ea
7 parameter(s): g = 0.015, α = 0.33, γ = 0.5, δ = 0.1, λ = 0.97, ρ = 0.03, β = @link 1 / (1 + ρ)
6 equations(s): 
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 - α))
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]

We can see the flags and the options of the model.

julia> m.flagsModelFlags
        linear = false
   ssZeroSlope = true
julia> m.options8 Options:
    shift = 10
    tol = 1.0e-12
    maxiter = 100
    verbose = true
    variant = :default
    substitutions = false
    warn = Options(no_t=true)
    factorization = :default

We can also examine individual components using the commands parameters, variables, shocks and equations.

julia> parameters(m)Parameters{ModelParam} with 7 entries:
  :g => 0.015
  :α => 0.33
  :γ => 0.5
  :δ => 0.1
  :λ => 0.97
  :ρ => 0.03
  :β => @link 1 / (1 + ρ)
julia> variables(m)6-element Vector{ModelVariable}:
 "Consumption" @log C
 "Capital Stock" @log K
 "Labour" @log L
 "Real Wage" @log w
 "Real Rental Rate" @log r
 "Technological shock" @log A
julia> shocks(m)1-element Vector{ModelVariable}:
 ea
julia> equations(m)OrderedCollections.OrderedDict{Symbol, Equation} with 6 entries:
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 -…
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)…
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]

Setting the model parameters

We must not change any part of the model in the active Julia session except for the values of the model parameters and steady state constraints if any (see the Smets and Wouters (2007) tutorial). If we want to add variables, shocks, or equations, we must do so in the model module file and restart a new Julia session to load the new model.

When it comes to the model parameters, we can access them by their names from the model object using the dot notation.

julia> m.β # read a parameter value0.970873786407767
julia> m.α = 0.33 # modify a parameter value0.33

Parameters can be linked to other parameters with the macro @link:

julia> parameters(m)[:β]@link 1 / (1 + ρ)

If $\rho$ changes, $\beta$ will automatically be updated.

julia> m.ρ = 0.050.05
julia> m.β0.9523809523809523
julia> m.ρ = 0.030.03
julia> m.β0.970873786407767

However, the dynamic links only work with parameter values. Otherwise, the function update_links! needs to be called to refresh all the links.

update_links!(m)

Important note

When do we need to call update_links!? Links will not be automatically updated if:

Links contain a reference outside the model parameters, such as a global variable, the steady state or another model object;
A parameter is not a number, such as if an element of a parameter vector is updated.

Part 3: The steady state solution

The steady state is a special solution of the dynamic system that remains constant over time. It is important on its own, but also it can be useful in several ways. For example, linearizing the model requires a particular solution about which to linearize, and the steady state is typically used for this purpose.

In addition to the steady state, we also consider another kind of special solution which grows linearly in time. If we know that the steady state solution is constant (i.e., its slope is zero), we can set the model flag ssZeroSlope to true. This is not required; however in a large model it might help the steady state solver converge faster to the solution.

The model object m stores information about the steady state. This includes the steady state solution itself, as well as a (possibly empty) set of additional constraints that apply only to the steady state. This information can be accessed via m.sstate.

julia> m.sstateSteady state not solved.
No additional constraints.

Solving for the steady state

The steady state solution is stored within the model object. Before solving, we have to specify an initial condition. If the model is linear, this makes no difference, but in a non-linear model a good or a bad initial guess might be the difference between success and failure of the steady state solver.

We specify the initial guess by calling clear_sstate!. This call removes any previously stored solution, sets the initial guess, and runs the pre-solve pass of the steady state solver. The initial guess can be given with the lvl and slp arguments; if not provided, an initial guess is chosen automatically.

Once that's done, we call sssolve! to find the steady state. We can see below that sssolve! cannot find a steady state solution.

julia> clear_sstate!(m)┌ Info: Presolved equation _EQ6 for #A#lvl# = 0.0
└     :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
julia> sssolve!(m);[ Info: Presolved 1 equations for 9 variables.
[ Info: Solving 5 equations for 5 variables.
[ Info: 1, NR, || Fx || = 1.1986872541507982, || Δx || = 1.1816556883242957
[ Info: 2, NR, || Fx || = 0.16988440829072382, || Δx || = 1.0183612282718562
[ Info: 3, NR, || Fx || = 0.0345712391855286, || Δx || = 0.3592981191065505
[ Info: 4, NR, || Fx || = 0.005946907573081361, || Δx || = 0.038180676933075884
[ Info: 5, NR, || Fx || = 4.49976637018068e-5, || Δx || = 0.0003268145232176359
[ Info: 6, NR, || Fx || = 4.12055989329474e-9, || Δx || = 2.475594632382498e-8
[ Info: 7, NR, || Fx || = 8.881784197001252e-16, || Δx || = 5.408231995098024e-16

Sometimes the Newton-Raphson solution algorithm, which is used by default because it is the fastest, fails to converge. If this happens, we can use method=:auto, which starts with the Levenberg-Marquardt algorithm and automatically switches to Newton-Raphson when it starts to converge.

julia> clear_sstate!(m)┌ Info: Presolved equation _EQ6 for #A#lvl# = 0.0
└     :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
julia> sssolve!(m; method = :auto);[ Info: Presolved 1 equations for 9 variables.
[ Info: Solving 5 equations for 5 variables.
[ Info: 0, || Fx || = 1.1986872541507982, || Δx || = 1.0436327779027528, lambda = 0.01
[ Info: 1, LM, || Fx || = 0.20958325671575825, || Δx || = 0.5859341014962702, lambda = 0.001
[ Info: 2, LM, || Fx || = 0.07179443905872815, || Δx || = 0.5655225835051327, lambda = 0.0001
[ Info: 3, LM, || Fx || = 0.014088859092551864, || Δx || = 0.16743352172313836, lambda = 1.0e-5
[ Info: 4, LM, || Fx || = 0.0013172231846176885, || Δx || = 0.010145059776807729, lambda = 1.0000000000000002e-6
[ Info: 5, LM, || Fx || = 4.166822614060095e-6, || Δx || = 4.225487059724547e-5, lambda = 1.0000000000000002e-7
[ Info: 6, NR, || Fx || = 6.396606579797924e-10, || Δx || = 8.183348643522827e-9
[ Info: 7, NR, || Fx || = 2.2204460492503136e-16, || Δx || = 2.35638061058079e-16

The function sssolve! returns a Vector{Float64} containing the steady state solution, and it also writes that solution into the model object. The vector is of length 2*nvariables(m) and contains the level and the slope for each variable.

If in doubt, we can use check_sstate to make sure the steady state solution stored in the model object indeed satisfies the steady state system of equations. This function returns the number of equations that are not satisfied. A value of 0 is what we want to see. In verbose mode, it also lists the problematic equations and their residuals.

julia> check_sstate(m)[ Info: All steady state equations are satisfied.
0

Examining the steady state

We can access the steady state solution via m.sstate using the dot notation.

julia> true_C = m.sstate.C.level1.0363974355242054

We can also assign new values to the steady state solution, but we should be careful to make sure it remains a valid steady state solution.

julia> m.sstate.C.level = true_C + 0.11.1363974355242055
julia> @test check_sstate(m) > 0┌ Warning: The following 2 steady state equations are not satisfied: 
│ E5  res= 0.1        :_EQ5 => K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
│ E2  res= 0.0921124  :_EQ2 => L[t] ^ γ * C[t] = w[t]
└ @ StateSpaceEcon.SteadyStateSolver ~/.julia/packages/StateSpaceEcon/zFS92/src/steadystate/diagnose.jl:49
Test Passed
  Expression: check_sstate(m) > 0
   Evaluated: 2 > 0

As the code above shows, a wrong steady state solution (based on the specified precision in the tol option) will result in one or more equation not being satisfied. Let's put back the correct value.

julia> m.sstate.C.level = true_C1.0363974355242054
julia> @test check_sstate(m) == 0[ Info: All steady state equations are satisfied.
Test Passed
  Expression: check_sstate(m) == 0
   Evaluated: 0 == 0

We can examine the entire steady state solution with printsstate.

julia> printsstate(m)Steady State Solution:
   C = 1.0364 
   K = 3.22922
   L = 0.93667
   w = 1.00304
   r = 0.14545
   A = 1.0    
  ea = 0.0

Part 4: Impulse response

Simulation plan

Before we can simulate the model, we have to decide on the length of the simulation and what data is available for each period, i.e., what values are known (exogenous). This is done with an object of type Plan.

To create a plan, all we need is the model object and a range for the simulation.

julia> sim_rng = 2000Q1:2039Q42000Q1:2039Q4
julia> p = Plan(m, sim_rng)Plan{MIT{Quarterly{3}}} with range 1999Q4:2040Q1
  1999Q4:2040Q1 → ea

The plan shows us the list of exogenous values (variables or shocks) for each period or sub-range of the simulation. By default, all shocks are exogenous and all variables are endogenous.

We also see that the range of the plan has been extended before and after the simulation range. This is necessary because we need to set initial and final conditions. The number of periods for initial conditions is equal to the largest lag in the model. Similarly, final conditions have to be imposed over as many periods as the largest lead.

julia> p.range          # the full range of the plan1999Q4:2040Q1
julia> init_rng = first(p.range):first(sim_rng)-1   # the range for initial conditions1999Q4:1999Q4
julia> final_rng = last(sim_rng)+1:last(p.range)     # the range for final conditions2040Q1:2040Q1
julia> @test length(init_rng) == m.maxlagTest Passed
  Expression: length(init_rng) == m.maxlag
   Evaluated: 1 == 1
julia> @test length(final_rng) == m.maxleadTest Passed
  Expression: length(final_rng) == m.maxlead
   Evaluated: 1 == 1

The function exportplan can be used to save a plan to a TXT or CSV file which can be opened to visualize the plan. Alternatively, the function importplan can load the plan back into Julia from the TXT or CSV file.

Exogenous data

We have to provide the data for the simulation. We start with all zeros and fill in the external data, which must include initial conditions for all variable and shocks, exogenous values (according to the plan), and possibly final conditions.

Initial conditions

In this example, we want to simulate an impulse response, so it makes sense to start from the steady state.

julia> exog = steadystatedata(m, p)162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C       K        L        w        r        A     ea
  1999Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
           ⋮       ⋮        ⋮        ⋮        ⋮        ⋮     ⋮
  2038Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2040Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0

Final conditions

For the final conditions, we can use the steady state again, because we expect that the economy will eventually return to it if the simulation is sufficiently long past the last shock. We can do this by assigning the values of the steady state to the final periods after the simulation, similarly to what we did with the initial conditions.

Alternatively, we can specify that we want to use the steady state in the call to simulate by passing fctype=fclevel. Yet another possibility is to set the final condition so that the solution slope matches the slope of the steady state by setting fctype=fcslope. In both cases, we do not need to set anything in the exogenous data array because those values would be ignored.

Pro tip

In the simple RBC model, the two ways of using the steady state for final conditions (level or slope) are equivalent, because the steady state here is stationary and unique. In models where the steady state has non-zero slope, or the steady state has zero slope but the level is not unique, we should use fctype=fcslope.

If the steady state is not solved, or if we prefer not to depend on it, we can use fctype=fcnatural. The final conditions will be constructed assuming that in the last two periods of the simulation the solution grows at the same rate, i.e., it has settled into its balanced growth. For a stationary model, the simulation needs to be long enough so that variables do not change anymore. In a model where the steady state has non-zero slope, non-stationary variables have to grow at a stable pace by the end of the simulation.

We can set the default option for the simple RBC model outside the model dedicated module.

julia> m.options.fctype = fcnaturalFCConstRate()

Otherwise, this option can be set within the model module but the StateSpaceEcon package must be loaded within the module in addition to ModelBaseEcon. For instance:

module simple_RBC
    using ModelBaseEcon
    using StateSpaceEcon
    const model = Model()
    model.flags.ssZeroSlope = true
    setoption!(model) do o
        o.tol = 1e-12
        o.maxiter = 100
        o.substitutions = false
        o.factorization = :default
        o.verbose = true
        o.fctype = fcnatural   # requires StateSpaceEcon
    end # options
    # Rest of the model...
end

A quick sanity check

If we were to run a simulation where the economy started in the steady state and there were no shocks at all, we'd expect that the economy would remain in steady state forever.

julia> ss = simulate(m, p, exog);[ Info: Simulating 2000Q1:2039Q4
[ Info: 1, || Fx || = 8.881784197001252e-16
julia> @test ss ≈ steadystatedata(m, p)Test Passed
  Expression: ss ≈ steadystatedata(m, p)
   Evaluated: 162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C       K        L        w        r        A     ea
  1999Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
           ⋮       ⋮        ⋮        ⋮        ⋮        ⋮     ⋮
  2038Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2040Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0   ≈ 162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C       K        L        w        r        A     ea
  1999Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2000Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2001Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
           ⋮       ⋮        ⋮        ⋮        ⋮        ⋮     ⋮
  2038Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2038Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q2 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q3 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2039Q4 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0
  2040Q1 : 1.0364  3.22922  0.93667  1.00304  0.14545  1.0   0.0

The simulated data, ss, should equal (up to the accuracy of the solution) the steady state data. Similar to steadystatedata, we can use zerodata to create a data set containing zeros to work in the deviation from the steady state solution.

julia> zz = simulate(m, p, zerodata(m, p); deviation = true);[ Info: Simulating 2000Q1:2039Q4
[ Info: 1, || Fx || = 8.881784197001252e-16
julia> @test zz ≈ zerodata(m, p)Test Passed
  Expression: zz ≈ zerodata(m, p)
   Evaluated: 162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C     K     L     w     r     A     ea
  1999Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
           ⋮     ⋮     ⋮     ⋮     ⋮     ⋮     ⋮
  2038Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2040Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0   ≈ 162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C     K     L     w     r     A     ea
  1999Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2000Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2001Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
           ⋮     ⋮     ⋮     ⋮     ⋮     ⋮     ⋮
  2038Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2038Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q2 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q3 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2039Q4 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0
  2040Q1 : 1.0   1.0   1.0   1.0   1.0   1.0   0.0

Exogenous data

All shocks are exogenous by default. All we have left to do is to set the value of the shock.

Let's say that we want to shock ea for the first four quarters by 0.1.

julia> exog[sim_rng[1:4], :ea] .= 0.1;
julia> exog[shocks(m)]162×1 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (ea):
           ea
  1999Q4 : 0.0
  2000Q1 : 0.1
  2000Q2 : 0.1
  2000Q3 : 0.1
  2000Q4 : 0.1
  2001Q1 : 0.0
  2001Q2 : 0.0
  2001Q3 : 0.0
  2001Q4 : 0.0
           ⋮
  2038Q1 : 0.0
  2038Q2 : 0.0
  2038Q3 : 0.0
  2038Q4 : 0.0
  2039Q1 : 0.0
  2039Q2 : 0.0
  2039Q3 : 0.0
  2039Q4 : 0.0
  2040Q1 : 0.0

Running the simulation

We call simulate, providing the model, the exogenous data, and the plan. We also specify the type of final condition we want to impose if we want to diverge from the option setting saved in the model.

julia> irf = simulate(m, p, exog; fctype=fcslope)[ Info: Simulating 2000Q1:2039Q4
[ Info: 1, || Fx || = 0.1, || Δx || = 0.4765323665515438
[ Info: 2, || Fx || = 0.1019555935796923, || Δx || = 0.02844812771090943
[ Info: 3, || Fx || = 0.0010720268544233136, || Δx || = 0.0003612138510853425
[ Info: 4, || Fx || = 1.9245162352632406e-7, || Δx || = 5.531157642768839e-8
[ Info: 5, || Fx || = 6.217248937900877e-15
162×7 MVTSeries{Quarterly} with range 1999Q4:2040Q1 and variables (C,K,L,w,r,…):
           C        K        L         w        r         A        ea
  1999Q4 : 1.0364   3.22922  0.93667   1.00304  0.14545   1.0      0.0
  2000Q1 : 1.2142   3.12762  0.873101  1.13454  0.153353  1.10517  0.1
  2000Q2 : 1.25465  3.20196  0.931434  1.21087  0.180277  1.21774  0.1
  2000Q3 : 1.32341  3.45397  0.987461  1.31509  0.202752  1.33788  0.1
  2000Q4 : 1.41709  3.90852  1.04609   1.44938  0.219453  1.46574  0.1
  2001Q1 : 1.48661  4.27475  1.02296   1.50358  0.196733  1.44902  0.0
  2001Q2 : 1.53718  4.56342  1.00461   1.54072  0.181016  1.43299  0.0
  2001Q3 : 1.57287  4.78587  0.989975  1.56497  0.169724  1.4176   0.0
  2001Q4 : 1.59686  4.95278  0.978242  1.57939  0.161391  1.40284  0.0
           ⋮        ⋮        ⋮         ⋮        ⋮         ⋮        ⋮
  2038Q1 : 1.04378  3.25138  0.936002  1.00983  0.145261  1.0041   0.0
  2038Q2 : 1.04361  3.2496   0.935867  1.00959  0.145277  1.00397  0.0
  2038Q3 : 1.04345  3.24758  0.935696  1.00935  0.145295  1.00385  0.0
  2038Q4 : 1.04332  3.24525  0.935479  1.0091   0.145316  1.00374  0.0
  2039Q1 : 1.04321  3.24253  0.935204  1.00885  0.145341  1.00363  0.0
  2039Q2 : 1.04313  3.2393   0.934855  1.00858  0.145371  1.00352  0.0
  2039Q3 : 1.04309  3.23542  0.934413  1.0083   0.145407  1.00341  0.0
  2039Q4 : 1.04309  3.23071  0.933853  1.008    0.14545   1.00331  0.0
  2040Q1 : 1.04309  3.23071  0.933853  1.008    0.14545   1.00331  0.0

We can now take a look at how some of the variables in the model have responded to this shock. We use plot from the Plots package. We specify the variables we want to plot using vars and the names of the datasets being plotted (for the legend) in the labels option.

julia> plot(ss, irf,
            vars=m.variables, # variables to plot are taken from the model
            legend= :none,
            linewidth=1.5,
           );

Part 5: Stochastic shocks simulation

Now let's run a simulation with stochastic shocks. We will have random shocks over two years and then have no shocks for several years afterwards to allow time for the economy to return to its steady state.

julia> sim_rng = 2000Q1:2049Q4      # simulate 50 years starting 20002000Q1:2049Q4
julia> shk_rng = 2004Q1 .+ (0:7)    # shock 8 quarters starting in 20042004Q1:2005Q4
julia> p = Plan(m, sim_rng)Plan{MIT{Quarterly{3}}} with range 1999Q4:2050Q1
  1999Q4:2050Q1 → ea
julia> exog = steadystatedata(m, p);

The distribution of the shock is assumed normal with mean zero. We use packages Distributions and Random to draw the necessary random values.

julia> shk_dist = (ea = Normal(0.0, 0.10),);
julia> for (shk, dist) in pairs(shk_dist)
           exog[shk_rng, shk] .= rand(dist, length(shk_rng))
       end
julia> exog[shk_rng, shocks(m)]8×1 MVTSeries{Quarterly} with range 2004Q1:2005Q4 and variables (ea):
           ea
  2004Q1 :  0.0970656
  2004Q2 : -0.0979218
  2004Q3 :  0.0901861
  2004Q4 : -0.00328031
  2005Q1 : -0.0600792
  2005Q2 : -0.144518
  2005Q3 :  0.270742
  2005Q4 :  0.152445

Now we are ready to simulate. We can set the shocks to be anticipated or unanticipated by setting the anticipate parameter in simulate.

julia> sim_a = simulate(m, p, exog; fctype=fcnatural, anticipate=true);[ Info: Simulating 2000Q1:2049Q4
[ Info: 1, || Fx || = 0.270742394171578, || Δx || = 0.47990762680967675
[ Info: 2, || Fx || = 0.10620587424829253, || Δx || = 0.035407593176925926
[ Info: 3, || Fx || = 0.001349189701278064, || Δx || = 0.0007724978977938846
[ Info: 4, || Fx || = 5.379286349693757e-7, || Δx || = 1.9070319491169272e-7
[ Info: 5, || Fx || = 5.3290705182007514e-14
julia> sim_u = simulate(m, p, exog; fctype=fcnatural, anticipate=false);[ Info: Simulating 2000Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 8.881784197001252e-16
[ Info: Simulating 2004Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.09706563288552145, || Δx || = 0.12026920951626652
[ Info: 2, || Fx || = 0.006588590235284464, || Δx || = 0.0019856757942954876
[ Info: 3, || Fx || = 4.14290387418248e-6, || Δx || = 1.6343216914677161e-6
[ Info: 4, || Fx || = 5.444533712761768e-13
[ Info: Simulating 2004Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.09792184115351996, || Δx || = 0.12130438008032718
[ Info: 2, || Fx || = 0.006495775556829564, || Δx || = 0.002014711089265727
[ Info: 3, || Fx || = 4.001809668530143e-6, || Δx || = 1.668354210102624e-6
[ Info: 4, || Fx || = 6.430411758628907e-13
[ Info: Simulating 2004Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.0901860883594094, || Δx || = 0.11200686388299859
[ Info: 2, || Fx || = 0.005735183748366346, || Δx || = 0.0017273431414442885
[ Info: 3, || Fx || = 3.142854740012524e-6, || Δx || = 1.2364076727779598e-6
[ Info: 4, || Fx || = 3.1796787425264483e-13
[ Info: Simulating 2004Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.003280312924463874, || Δx || = 0.004055609485314191
[ Info: 2, || Fx || = 7.710715913766819e-6, || Δx || = 2.2515584608206066e-6
[ Info: 3, || Fx || = 5.38857847232066e-12, || Δx || = 2.0956364093493184e-12
[ Info: 4, || Fx || = 1.7763568394002505e-15
[ Info: Simulating 2005Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.060079222335556216, || Δx || = 0.07415933021700599
[ Info: 2, || Fx || = 0.002523016511584153, || Δx || = 0.0007649398679566867
[ Info: 3, || Fx || = 5.965414908715161e-7, || Δx || = 2.4096778326116305e-7
[ Info: 4, || Fx || = 1.509903313490213e-14
[ Info: Simulating 2005Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.1445177115286231, || Δx || = 0.17947994854278057
[ Info: 2, || Fx || = 0.013372088709408203, || Δx || = 0.004548211167857031
[ Info: 3, || Fx || = 1.8820916293904588e-5, || Δx || = 8.479751541221714e-6
[ Info: 4, || Fx || = 1.7321255540991842e-11, || Δx || = 6.867922130518875e-12
[ Info: 5, || Fx || = 1.7763568394002505e-15
[ Info: Simulating 2005Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.270742394171578, || Δx || = 0.3408674094796361
[ Info: 2, || Fx || = 0.05248823260022917, || Δx || = 0.015834622495266314
[ Info: 3, || Fx || = 0.0002725127323515153, || Δx || = 0.00010430822931571414
[ Info: 4, || Fx || = 1.7437651322893544e-9, || Δx || = 4.994253044464431e-10
[ Info: 5, || Fx || = 2.6645352591003757e-15
[ Info: Simulating 2005Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.15244478634355946, || Δx || = 0.18939421102003787
[ Info: 2, || Fx || = 0.01895380518281531, || Δx || = 0.004665682536218989
[ Info: 3, || Fx || = 2.7841937441763775e-5, || Δx || = 9.092485539887768e-6
[ Info: 4, || Fx || = 1.7984724820507836e-11, || Δx || = 4.325661718091807e-12
[ Info: 5, || Fx || = 2.6645352591003757e-15
[ Info: Simulating 2005Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 2.6645352591003757e-15

As before, we can review the responses of variables to the shock using plot.

julia> observed = collect(keys(m.autoexogenize)); # the observed variable is from the autoexogenize list
julia> ss = steadystatedata(m, p);
julia> plot(ss, sim_a, sim_u,
            vars=m.variables,
            labels=("SS", "Anticipated", "Unanticipated"),
            legend=[true (false for i = 2:length(m.variables))...],
            linewidth=1.5,
           );

We see that when the shock is anticipated, the variables start to react to them right away; in the unanticipated case, there is no movement until the technology shock actually hit.

Part 6: Backing out historical shocks

Now let's pretend that the simulated values for A are historical data and that we do not know the magnitude of the shock ea. We can treat the observed (simulated) values of the variable A as known by making them exogenous. At the same time we will make the shock endogenous, so that we can solve for its values during the simulation.

We use exogenize! and endogenize! to set up a plan in which the observed variable is exogenous and the shock is endogenous throughout the stochastic range.

julia> endogenize!(p, shocks(m), shk_rng);
julia> exogenize!(p, observed, shk_rng);
julia> pPlan{MIT{Quarterly{3}}} with range 1999Q4:2050Q1
  1999Q4:2003Q4 → ea
  2004Q1:2005Q4 → A
  2006Q1:2050Q1 → ea

Another possibility is to use the autoexogenize! command, which will use the default pairing provided in the model definition under @autoexogenize.

julia> autoexogenize!(p, m, shk_rng)Plan{MIT{Quarterly{3}}} with range 1999Q4:2050Q1
  1999Q4:2003Q4 → ea
  2004Q1:2005Q4 → A
  2006Q1:2050Q1 → ea

As we can see above, the plan now reflects our intentions.

Finally, we need to set up the exogenous data. This time we do not specify the shocks; instead, we assign the known data for the observed variables for the historic range. We start with initial conditions.

julia> exog = steadystatedata(m, p);

We take the observed data from the simulation above. We show the anticipated version first.

julia> for v in observed
           exog[shk_rng, v] .= sim_a[v]
       end
julia> back_a = simulate(m, p, exog; fctype=fcnatural, anticipate=true);[ Info: Simulating 2000Q1:2049Q4
[ Info: 1, || Fx || = 0.4823964779371739, || Δx || = 0.48237002414960234
[ Info: 2, || Fx || = 0.0681435972815212, || Δx || = 0.030919482533576785
[ Info: 3, || Fx || = 0.0010838529549586084, || Δx || = 0.000517597269556976
[ Info: 4, || Fx || = 3.2100832392245593e-7, || Δx || = 7.4978492747408e-8
[ Info: 5, || Fx || = 1.1546319456101628e-14

Now we show the unanticipated case.

julia> for v in observed
           exog[shk_rng, v] .= sim_u[v]
       end
julia> back_u = simulate(m, p, exog; fctype=fcnatural, anticipate=false);[ Info: Simulating 2000Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.4823964779371739, || Δx || = 0.2956585489358212
[ Info: 2, || Fx || = 0.07025272951507233, || Δx || = 0.012858030504466655
[ Info: 3, || Fx || = 0.00022810083683744153, || Δx || = 5.116289681930143e-5
[ Info: 4, || Fx || = 3.5166634049232925e-9, || Δx || = 8.75191444432339e-10
[ Info: 5, || Fx || = 1.7763568394002505e-15
[ Info: Simulating 2004Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.14293714212341868, || Δx || = 0.11968849679571193
[ Info: 2, || Fx || = 0.0039160310197603465, || Δx || = 0.0008535631486192461
[ Info: 3, || Fx || = 1.3306363788601061e-6, || Δx || = 2.7409302394847197e-7
[ Info: 4, || Fx || = 1.2523315717771766e-13
[ Info: Simulating 2004Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.14785801040855695, || Δx || = 0.12131434143649766
[ Info: 2, || Fx || = 0.003873480525874662, || Δx || = 0.0008500424359180719
[ Info: 3, || Fx || = 1.2821444217436806e-6, || Δx || = 2.705656920017249e-7
[ Info: 4, || Fx || = 1.2168044349891716e-13
[ Info: Simulating 2004Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.13239204547879435, || Δx || = 0.1110923468341227
[ Info: 2, || Fx || = 0.0033906968559094253, || Δx || = 0.0007196690906217691
[ Info: 3, || Fx || = 9.592011638304143e-7, || Δx || = 1.8930570100862075e-7
[ Info: 4, || Fx || = 6.039613253960852e-14
[ Info: Simulating 2004Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.005143404106888738, || Δx || = 0.004055288179051868
[ Info: 2, || Fx || = 4.610862807119531e-6, || Δx || = 9.79144748827941e-7
[ Info: 3, || Fx || = 1.7905676941154525e-12, || Δx || = 3.6634272523407645e-13
[ Info: Simulating 2005Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.0912471181427108, || Δx || = 0.07450048682764063
[ Info: 2, || Fx || = 0.0015088512737859716, || Δx || = 0.0003139171072043189
[ Info: 3, || Fx || = 1.809300140820369e-7, || Δx || = 3.578320092502948e-8
[ Info: 4, || Fx || = 1.7763568394002505e-15
[ Info: Simulating 2005Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.195105802605549, || Δx || = 0.18169741864850192
[ Info: 2, || Fx || = 0.00826152894861476, || Δx || = 0.0018867571290378028
[ Info: 3, || Fx || = 4.510120714229515e-6, || Δx || = 1.2558364186981774e-6
[ Info: 4, || Fx || = 1.581845765485923e-12, || Δx || = 5.646681657926779e-13
[ Info: Simulating 2005Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.37472853602921496, || Δx || = 0.3327491641860331
[ Info: 2, || Fx || = 0.03117492241188291, || Δx || = 0.006312276459952772
[ Info: 3, || Fx || = 7.625890824236592e-5, || Δx || = 1.3806830379011638e-5
[ Info: 4, || Fx || = 3.285496319449521e-10, || Δx || = 8.090093668500239e-11
[ Info: 5, || Fx || = 2.6645352591003757e-15
[ Info: Simulating 2005Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.2787566910532986, || Δx || = 0.19011224797378276
[ Info: 2, || Fx || = 0.013066833757724439, || Δx || = 0.002399013543160369
[ Info: 3, || Fx || = 1.1703631660253677e-5, || Δx || = 2.3851362059146357e-6
[ Info: 4, || Fx || = 1.0667022820598504e-11, || Δx || = 2.5415753302203105e-12
[ Info: 5, || Fx || = 2.6645352591003757e-15
[ Info: Simulating 2005Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 2.6645352591003757e-15

If we did everything correctly, the shocks we recovered must match the shocks we used when we simulated the data.

julia> @test sim_a[:ea] ≈ back_a[:ea]Test Passed
  Expression: sim_a[:ea] ≈ back_a[:ea]
   Evaluated: 202-element TSeries{Quarterly} with range 1999Q4:2050Q1:
  1999Q4 : 0.0
  2000Q1 : 0.0
  2000Q2 : 0.0
  2000Q3 : 0.0
  2000Q4 : 0.0
  2001Q1 : 0.0
  2001Q2 : 0.0
  2001Q3 : 0.0
  2001Q4 : 0.0
    ⋮
  2047Q4 : 0.0
  2048Q1 : 0.0
  2048Q2 : 0.0
  2048Q3 : 0.0
  2048Q4 : 0.0
  2049Q1 : 0.0
  2049Q2 : 0.0
  2049Q3 : 0.0
  2049Q4 : 0.0
  2050Q1 : 0.0 ≈ 202-element TSeries{Quarterly} with range 1999Q4:2050Q1:
  1999Q4 : 0.0
  2000Q1 : 0.0
  2000Q2 : 0.0
  2000Q3 : 0.0
  2000Q4 : 0.0
  2001Q1 : 0.0
  2001Q2 : 0.0
  2001Q3 : 0.0
  2001Q4 : 0.0
    ⋮
  2047Q4 : 0.0
  2048Q1 : 0.0
  2048Q2 : 0.0
  2048Q3 : 0.0
  2048Q4 : 0.0
  2049Q1 : 0.0
  2049Q2 : 0.0
  2049Q3 : 0.0
  2049Q4 : 0.0
  2050Q1 : 0.0
julia> @test sim_u[:ea] ≈ back_u[:ea]Test Passed
  Expression: sim_u[:ea] ≈ back_u[:ea]
   Evaluated: 202-element TSeries{Quarterly} with range 1999Q4:2050Q1:
  1999Q4 : 0.0
  2000Q1 : 0.0
  2000Q2 : 0.0
  2000Q3 : 0.0
  2000Q4 : 0.0
  2001Q1 : 0.0
  2001Q2 : 0.0
  2001Q3 : 0.0
  2001Q4 : 0.0
    ⋮
  2047Q4 : 0.0
  2048Q1 : 0.0
  2048Q2 : 0.0
  2048Q3 : 0.0
  2048Q4 : 0.0
  2049Q1 : 0.0
  2049Q2 : 0.0
  2049Q3 : 0.0
  2049Q4 : 0.0
  2050Q1 : 0.0 ≈ 202-element TSeries{Quarterly} with range 1999Q4:2050Q1:
  1999Q4 : 0.0
  2000Q1 : 0.0
  2000Q2 : 0.0
  2000Q3 : 0.0
  2000Q4 : 0.0
  2001Q1 : 0.0
  2001Q2 : 0.0
  2001Q3 : 0.0
  2001Q4 : 0.0
    ⋮
  2047Q4 : 0.0
  2048Q1 : 0.0
  2048Q2 : 0.0
  2048Q3 : 0.0
  2048Q4 : 0.0
  2049Q1 : 0.0
  2049Q2 : 0.0
  2049Q3 : 0.0
  2049Q4 : 0.0
  2050Q1 : 0.0

Moreover, we must have the unobserved variables match as well. In fact, all the data must match over the entire simulation range.

julia> @test sim_a ≈ back_aTest Passed
  Expression: sim_a ≈ back_a
   Evaluated: 202×7 MVTSeries{Quarterly} with range 1999Q4:2050Q1 and variables (C,K,L,w,r,…):
           C        K        L        … r         A         ea
  1999Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q1 : 1.03751  3.22689  0.935458 … 0.145324  1.0        0.0
  2000Q2 : 1.03744  3.22438  0.935271 … 0.145375  1.0        0.0
  2000Q3 : 1.03741  3.22155  0.93501  … 0.145423  1.0        0.0
  2000Q4 : 1.03743  3.21827  0.93466  … 0.145472  1.0        0.0
  2001Q1 : 1.03751  3.21436  0.934202 … 0.145524  1.0        0.0
  2001Q2 : 1.03764  3.20961  0.93361  … 0.145581  1.0        0.0
  2001Q3 : 1.03783  3.20375  0.93285  … 0.145646  1.0        0.0
  2001Q4 : 1.0381   3.19648  0.931881 … 0.145722  1.0        0.0
           ⋮        ⋮        ⋮        ⋱ ⋮         ⋮         ⋮
  2048Q1 : 1.03941  3.24034  0.936644 … 0.145353  1.00172    0.0
  2048Q2 : 1.03932  3.24008  0.936654 … 0.145356  1.00167    0.0
  2048Q3 : 1.03923  3.23984  0.936667 … 0.145358  1.00162    0.0
  2048Q4 : 1.03914  3.23963  0.936682 … 0.145359  1.00157    0.0
  2049Q1 : 1.03905  3.23946  0.936701 … 0.145361  1.00152    0.0
  2049Q2 : 1.03896  3.23933  0.936725 … 0.145362  1.00148    0.0
  2049Q3 : 1.03888  3.23925  0.936754 … 0.145362  1.00143    0.0
  2049Q4 : 1.03879  3.23923  0.936791 … 0.145362  1.00139    0.0
  2050Q1 : 1.0387   3.23921  0.936828 … 0.145362  1.00135    0.0        ≈ 202×7 MVTSeries{Quarterly} with range 1999Q4:2050Q1 and variables (C,K,L,w,r,…):
           C        K        L        … r         A         ea
  1999Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q1 : 1.03751  3.22689  0.935458 … 0.145324  1.0        0.0
  2000Q2 : 1.03744  3.22438  0.935271 … 0.145375  1.0        0.0
  2000Q3 : 1.03741  3.22155  0.93501  … 0.145423  1.0        0.0
  2000Q4 : 1.03743  3.21827  0.93466  … 0.145472  1.0        0.0
  2001Q1 : 1.03751  3.21436  0.934202 … 0.145524  1.0        0.0
  2001Q2 : 1.03764  3.20961  0.93361  … 0.145581  1.0        0.0
  2001Q3 : 1.03783  3.20375  0.93285  … 0.145646  1.0        0.0
  2001Q4 : 1.0381   3.19648  0.931881 … 0.145722  1.0        0.0
           ⋮        ⋮        ⋮        ⋱ ⋮         ⋮         ⋮
  2048Q1 : 1.03941  3.24034  0.936644 … 0.145353  1.00172    0.0
  2048Q2 : 1.03932  3.24008  0.936654 … 0.145356  1.00167    0.0
  2048Q3 : 1.03923  3.23984  0.936667 … 0.145358  1.00162    0.0
  2048Q4 : 1.03914  3.23963  0.936682 … 0.145359  1.00157    0.0
  2049Q1 : 1.03905  3.23946  0.936701 … 0.145361  1.00152    0.0
  2049Q2 : 1.03896  3.23933  0.936725 … 0.145362  1.00148    0.0
  2049Q3 : 1.03888  3.23925  0.936754 … 0.145362  1.00143    0.0
  2049Q4 : 1.03879  3.23923  0.936791 … 0.145362  1.00139    0.0
  2050Q1 : 1.0387   3.23921  0.936828 … 0.145362  1.00135    0.0
julia> @test sim_u ≈ back_uTest Passed
  Expression: sim_u ≈ back_u
   Evaluated: 202×7 MVTSeries{Quarterly} with range 1999Q4:2050Q1 and variables (C,K,L,w,r,…):
           C        K        L        … r         A         ea
  1999Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q1 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q2 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q3 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q1 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q2 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q3 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
           ⋮        ⋮        ⋮        ⋱ ⋮         ⋮         ⋮
  2048Q1 : 1.03941  3.24034  0.936644 … 0.145353  1.00172    0.0
  2048Q2 : 1.03932  3.24008  0.936654 … 0.145356  1.00167    0.0
  2048Q3 : 1.03923  3.23984  0.936667 … 0.145358  1.00162    0.0
  2048Q4 : 1.03914  3.23963  0.936682 … 0.145359  1.00157    0.0
  2049Q1 : 1.03905  3.23946  0.936701 … 0.145361  1.00152    0.0
  2049Q2 : 1.03896  3.23933  0.936725 … 0.145362  1.00148    0.0
  2049Q3 : 1.03888  3.23925  0.936754 … 0.145362  1.00143    0.0
  2049Q4 : 1.03879  3.23923  0.936791 … 0.145362  1.00139    0.0
  2050Q1 : 1.0387   3.23921  0.936828 … 0.145362  1.00135    0.0        ≈ 202×7 MVTSeries{Quarterly} with range 1999Q4:2050Q1 and variables (C,K,L,w,r,…):
           C        K        L        … r         A         ea
  1999Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q1 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q2 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q3 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2000Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q1 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q2 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q3 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
  2001Q4 : 1.0364   3.22922  0.93667  … 0.14545   1.0        0.0
           ⋮        ⋮        ⋮        ⋱ ⋮         ⋮         ⋮
  2048Q1 : 1.03941  3.24034  0.936644 … 0.145353  1.00172    0.0
  2048Q2 : 1.03932  3.24008  0.936654 … 0.145356  1.00167    0.0
  2048Q3 : 1.03923  3.23984  0.936667 … 0.145358  1.00162    0.0
  2048Q4 : 1.03914  3.23963  0.936682 … 0.145359  1.00157    0.0
  2049Q1 : 1.03905  3.23946  0.936701 … 0.145361  1.00152    0.0
  2049Q2 : 1.03896  3.23933  0.936725 … 0.145362  1.00148    0.0
  2049Q3 : 1.03888  3.23925  0.936754 … 0.145362  1.00143    0.0
  2049Q4 : 1.03879  3.23923  0.936791 … 0.145362  1.00139    0.0
  2050Q1 : 1.0387   3.23921  0.936828 … 0.145362  1.00135    0.0

Part 7: Model variants and solvers

As of version 0.4 of StateSpaceEcon, there are two keyword arguments that control how the model is handled and which solver is used. These are variant and solver. The currently available variants are :default (the model is taken as given), :linearize and :selective_linearize. Currently there are two solvers available, namely :stackedtime (which is the default) and :firstorder.

In order to use the :linearize variant, you must first solve for the steady state, as already explained. Once the steady state solution is stored in the model instance, you all linearize!, which creates the linearization of the model about its steady state and sets the default variant to :linearize.

You can check or change the default variant via m.variant.

Once the linearized model is available, you can use either continue to use the stacked time solver or you can start using the first order solver. For this you must first call solve! with solver=:firstorder, after which you can pass solver=:firstorder to simulate (and other functions that use a solver).

julia> m.variant:default
julia> linearize!(m)6 variable(s): C, K, L, w, r, A
1 shock(s): ea
7 parameter(s): g = 0.015, α = 0.33, γ = 0.5, δ = 0.1, λ = 0.97, ρ = 0.03, β = @link 1 / (1 + ρ)
6 equations(s):
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 - α))
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
Maximum lag: 1
Maximum lead: 1
julia> m.variant:linearize
julia> solve!(m, solver = :firstorder)6 variable(s): C, K, L, w, r, A
1 shock(s): ea
7 parameter(s): g = 0.015, α = 0.33, γ = 0.5, δ = 0.1, λ = 0.97, ρ = 0.03, β = @link 1 / (1 + ρ)
6 equations(s):
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 - α))
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
Maximum lag: 1
Maximum lead: 1
julia> m.variant:linearize

Instead of linearizing all the equations, you can linearize only selected equations by creating the variant :selective_linearize with the command selective_linearize!. Once again, the steady state solution must e available for this call to succeed.

julia> selective_linearize!(m)6 variable(s): C, K, L, w, r, A
1 shock(s): ea
7 parameter(s): g = 0.015, α = 0.33, γ = 0.5, δ = 0.1, λ = 0.97, ρ = 0.03, β = @link 1 / (1 + ρ)
6 equations(s):
  :_EQ1 => @log C[t + 1] * (1 + g) = β * (C[t] * ((r[t + 1] + 1) - δ))
  :_EQ2 => @log L[t] ^ γ * C[t] = w[t]
  :_EQ3 => @log r[t] * (K[t - 1] / (1 + g)) ^ (1 - α) = α * (A[t] * L[t] ^ (1 - α))
  :_EQ4 => @log w[t] * L[t] ^ α = (1 - α) * (A[t] * (K[t - 1] / (1 + g)) ^ α)
  :_EQ5 => @lin K[t] + C[t] = A[t] * ((K[t - 1] / (1 + g)) ^ α * L[t] ^ (1 - α)) + (1 - δ) * (K[t - 1] / (1 + g))
  :_EQ6 => log(A[t]) = λ * log(A[t - 1]) + ea[t]
Maximum lag: 1
Maximum lead: 1
julia> m.options.variant:selective_linearize

The equation that will be linearized with this call must be specified in the model file by marking them with the macro @lin. For instance:

@lin K[t] + C[t] = A[t] * (K[t-1]/(1+g)) ^ α * (L[t]) ^ (1-α) + (1-δ) * (K[t-1]/(1+g))

The model variant can be reset back to the original by assigning it directly.

julia> m.variant = :default:default

In total, there are three variants:

:default: the model as given through its equations
:linearize: first-order approximation around its steady state
:selective_linearize: first-order approximation around its steady state for the equations preceded by the macro @lin..

In addition to the variant, the command simulate requires a solver. StateSpaceEcon.jl currently has two solvers:

The solver :stackedtime can be used with any variant. It is the only solver that can be used with default and :selective_linearize.
The solver :firstorder can only be used for the variant :linearize. In fact, when this solver is specified the variant is ignored.

To get the default solver, simply omit the solver= argument of the command simulate.

For demonstration purposes, we compare the three models for an unanticipated shock ea for the first four quarters by 0.1.

julia> p = Plan(m, sim_rng)Plan{MIT{Quarterly{3}}} with range 1999Q4:2050Q1
  1999Q4:2050Q1 → ea
julia> exog = zerodata(m,p);
julia> exog[sim_rng[1:4], :ea] .= 0.1;
julia> exog1 = simulate(m, p, exog; deviation = true, anticipate = false, variant = :default, solver = :stackedtime);[ Info: Simulating 2000Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.1, || Δx || = 0.1239050382107035
[ Info: 2, || Fx || = 0.007004521609296432, || Δx || = 0.002106828907350947
[ Info: 3, || Fx || = 4.674852061903323e-6, || Δx || = 1.8400259451160068e-6
[ Info: 4, || Fx || = 6.883382752675971e-13
[ Info: Simulating 2000Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.10000000000000003, || Δx || = 0.12390401277884854
[ Info: 2, || Fx || = 0.007582069866832519, || Δx || = 0.002046876729129827
[ Info: 3, || Fx || = 4.866789573121366e-6, || Δx || = 1.7423865582596924e-6
[ Info: 4, || Fx || = 6.536993168992922e-13
[ Info: Simulating 2000Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.09999999999999987, || Δx || = 0.12452576143119584
[ Info: 2, || Fx || = 0.00831434315355395, || Δx || = 0.0020002153740544255
[ Info: 3, || Fx || = 5.178028269270385e-6, || Δx || = 1.6707932693965105e-6
[ Info: 4, || Fx || = 6.465938895416912e-13
[ Info: Simulating 2000Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.10000000000000009, || Δx || = 0.12495060992984192
[ Info: 2, || Fx || = 0.009217677295228377, || Δx || = 0.0019668143569366242
[ Info: 3, || Fx || = 5.618262279405428e-6, || Δx || = 1.6213370956826984e-6
[ Info: 4, || Fx || = 6.67021993194794e-13
[ Info: Simulating 2000Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 6.67021993194794e-13
julia> exog2 = simulate(m, p, exog; deviation = true, anticipate = false, variant = :selective_linearize, solver = :stackedtime);[ Info: Simulating 2000Q1:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.1, || Δx || = 0.1239050382107035
[ Info: 2, || Fx || = 0.0005552102774712315, || Δx || = 0.000837774742334313
[ Info: 3, || Fx || = 1.2127288136193285e-8, || Δx || = 4.850684578325929e-8
[ Info: 4, || Fx || = 6.833658530540354e-16
[ Info: Simulating 2000Q2:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.10000000000000003, || Δx || = 0.12473573118137929
[ Info: 2, || Fx || = 0.0005825813604441, || Δx || = 0.0008608191404192145
[ Info: 3, || Fx || = 1.2876433779974526e-8, || Δx || = 4.985179133918102e-8
[ Info: 4, || Fx || = 6.970369648227198e-16
[ Info: Simulating 2000Q3:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.09999999999999987, || Δx || = 0.12534011781910437
[ Info: 2, || Fx || = 0.0006026344983122634, || Δx || = 0.0008773591256577829
[ Info: 3, || Fx || = 1.3394085930146538e-8, || Δx || = 5.085015039002847e-8
[ Info: 4, || Fx || = 7.979727989493313e-16
[ Info: Simulating 2000Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 0.10000000000000009, || Δx || = 0.12575293942981894
[ Info: 2, || Fx || = 0.0006166132243901846, || Δx || = 0.0008888249394585545
[ Info: 3, || Fx || = 1.3724292907499989e-8, || Δx || = 5.159662523869351e-8
[ Info: 4, || Fx || = 6.661338147750941e-16
[ Info: Simulating 2000Q4:2049Q4 with (1.0e-12, 100)
[ Info: 1, || Fx || = 6.661338147750941e-16
julia> exog3 = simulate(m, p, exog; deviation = true, anticipate = false, variant = :linearize, solver = :firstorder);
julia> plot(exog1, exog2, exog3,
            vars=m.variables,
            labels=("Stacked-Time", "Selective linearization", "Linearized"),
            linestyle=[:solid :dash :solid],
            legend=[true (false for i = 2:length(m.variables))...],
            linewidth=1.5,
           );

The first-order approximation

For anticipated shocks, the :firstorder solver is only available for empty plans. Empty plans have all the shocks as exogenous and all the variables as endogenous. For anticipated shocks with non-empty plans, use the :stackedtime solver.

Appendix

References

Villemot, S., 2013. First order approximation of stochastic models.