Perturbation of Neoclassical Model

Our goal here is to compute a linear approximation of solution to the neoclassical model, close to the steady-state.

Warm-up: install the ForwardDiff library. Use it to differentiate the function below. Check the jacobian function.

Note: the signature of function f needs to be fixed first to allow for dual numbers as arguments.

Create a NamedTuple to hold the model parameters.

Define two functions: - transition(z::Number, k::Number, i::Number, p)::Tuple{Number} which returns productivity and capital at date t+1 as a function of productivity, capital and investment at date t - arbitrage(z::Number, k::Number, i::Number, Z::Number, K::Number, I::Number, p)::Number which returns the residual of the euler equation (lower case variable for date t, upper case for date t+1)

Using multiple dispatch, define two variants of the same functions, that take vectors as input and output arguments: - arbitrage(s::Vector{T}, x::Vector{T}, S::Vector{T}, X::Vector{T}, p) where T<:Number - transition(s::Vector{T}, x::Vector{T}, p) where T<:Number

Write a function steady_state(p)::Tuple{Vector,Vector} which computes the steady-state of the model computed by hand. It returns two vectors, one for the states, one for the controls. Check that the steady-state satisfies the model equations.

The first order system satisfies: \[\begin{align}A s_t + B x_t + C s_{t+1} + D x_{t+1} & = & 0 \\\\ s_{t+1} & = & E s_t + F x_t \end{align}\]

Define a structure PerturbedModel to hold matrices A,B,C,D,E,F.

Write a function first_order_model(s::Vector, x::Vector, p)::PerturbedModel, which returns the first order model, given the steady-state and the calibration. Suggestion: use ForwardDiff.jl library.

We look for a linear solution \(x_t = X s_t\) . Write the matrix equation which X must satisfy. Write a function residual(X::Array, M::PerturbedModel) which computes the residual of this equation for a given X.

Write a function T(X, M::PerturbedModel) which implements the time iteration step.

Write function linear_time_iteration(X_0::Matrix, m::PerturbedModel)::Matrix which implements the time iteration algorithm. Apply it to X0 = rand(1,2) and check that the result satisfies the first order model.

Check blanchard Kahn Conditions

Define two linear operators L_S(U::Union{Vector, Matrix}, X_0::Matrix, m::PerturbedModel)::Matrix and L_T(U::Matrix, X_0::Matrix, m::PerturbedModel)::Matrix which implement the derivatives of the simulation and the time-iteration operator respectively.

Implement a function spectral_radius(f::Function)::Float64 which implements the power iteration method to compute the biggest eigenvalues of the two previously defined operators. Check that Blanchard-Kahn conditions are met.

Write a function simulate(s0::Vector, X::Matrix, p, T::Int64)::Tuple{Matrix, Matrix} to simulate the model over \(T\) periods (by using the formula \(\Delta s_t = (E + F X) \Delta s_{t-1}\). Return a matrix for the states (one line per date) and another matrix for the controls. Bonus: add a keyword option to compute variables levels or log-deviations. If possible, return a DataFrame object.

Make some nice plots.