Perturbation of Neoclassical Model

Pablo Winant

Our goal here is to compute a linear approximation of solution to the neoclassical model, c# execute: # freeze: true
# execute: # enabled: falselose ot the steady-state.

1. 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.

# function f(x::Vector{T}) where T <: Number
function f(x)
    a = x[1]
    b = x[2]
    x1 = a+b
    x2 = a*exp(b)
    return [x1,x2]
end

f (generic function with 1 method)

# your code here
using ForwardDiff
ForwardDiff.jacobian(f, [1.0, 2.0])
# library uses dual numbers so f should not be too specific about the type of x

2×2 Matrix{Float64}:
 1.0      1.0
 7.38906  7.38906

# # function f(x::Vector{T}) where T <: Number
# function g(x::Vector{Float64})
#     a = x[1]
#     b = x[2]
#     x1 = a+b
#     x2 = a*exp(b)
#     return [x1,x2]
# end
# g([1.0, 2.0])
# ForwardDiff.jacobian(g, [1.0, 2.0]) # doesn't work

MethodError: MethodError: no method matching g(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}})
The function `g` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  g(!Matched::Vector{Float64})
   @ Main ~/Teaching/econobits/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X36sZmlsZQ==.jl:2

MethodError: no method matching g(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}})

The function `g` exists, but no method is defined for this combination of argument types.



Closest candidates are:

  g(!Matched::Vector{Float64})

   @ Main ~/Teaching/econobits/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X36sZmlsZQ==.jl:2





Stacktrace:

 [1] vector_mode_dual_eval!(f::typeof(g), cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}}}, x::Vector{Float64})

   @ ForwardDiff ~/.julia/packages/ForwardDiff/egQMG/src/apiutils.jl:24

 [2] vector_mode_jacobian(f::typeof(g), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}}})

   @ ForwardDiff ~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:129

 [3] jacobian(f::typeof(g), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}}}, ::Val{true})

   @ ForwardDiff ~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:22

 [4] jacobian(f::typeof(g), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}}})

   @ ForwardDiff ~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:19

 [5] jacobian(f::typeof(g), x::Vector{Float64})

   @ ForwardDiff ~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:19

 [6] top-level scope

   @ ~/Teaching/econobits/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X36sZmlsZQ==.jl:10

1. Create a NamedTuple to hold the model parameters.

# your code here
model  = (;
    α = 0.3,
    β = 0.96,
    γ = 2.0,
    δ = 0.1,
    ρ = 0.9
)

(α = 0.3, β = 0.96, γ = 2.0, δ = 0.1, ρ = 0.9)

1. 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)

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

"""
p : parameters of the model
"""
function transition(z,k,i,p)
    (;ρ, δ) = p
    #same as
    # ρ = p.ρ
    # δ = p.δ

    Z = ρ*z
    K = (1-δ)*k + i
    return (Z, K)
end
transition( 0.0, 1.0, 0.1, model)

(0.0, 1.0)

function arbitrage(z,k,i,Z,K,I,p)
    (;ρ, δ, α, β, γ) = p

    y = exp(z)*k^α # lower case for today
    c = y-i
    Y = exp(Z)*K^α # upper case for tomorrow
    C = Y-I

    r = β*(C/c)^(-γ)*( (1-δ) + α*Y/K ) - 1

    return (r, ) # tuple with one value

end
arbitrage(0.0, 1.0, 0.1, 0.0, 1.0, 0.1, model)

(0.1519999999999999,)

# overload the definitions for vectors

transition(s, x, p) = [transition(s[1], s[2], x[1], p)...]

# [t...]  #unpacks the tuple t as a vector
transition(s, x, p) = [transition(s..., x..., p)...]

transition (generic function with 2 methods)

transition([0.0, 1.0], [0.1], model)

2-element Vector{Float64}:
 0.0
 1.0

arbitrage(s, x, S, X,  p) = [arbitrage(s..., x..., S..., X..., p)...]
arbitrage([0.0, 1.0], [0.1],[0.0, 1.0], [0.1], model)

1-element Vector{Float64}:
 0.1519999999999999

1. 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.

function steady_state(p)
    (;ρ, δ, α, β, γ) = p
    k = ((1/β-(1-δ))/α)^(1/(α-1))
    i = δ*k
    z = 0.0
    return [z,k], [i]
    
end

s_, x_ = steady_state(model)

([0.0, 2.920822149964071], [0.29208221499640713])

arbitrage(s_,x_,s_,x_,model)

1-element Vector{Float64}:
 0.0

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}\]

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

# your code here

1. 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.

using ForwardDiff
A = ForwardDiff.jacobian( u-> arbitrage(u, x_, s_, x_, model ), s_)
B = ForwardDiff.jacobian( u-> arbitrage(s_, u, s_, x_, model ), x_)
C = ForwardDiff.jacobian( u-> arbitrage(s_, x_, u, x_, model ), s_)
D = ForwardDiff.jacobian( u-> arbitrage(s_, x_, s_, u, model ), x_)
E = ForwardDiff.jacobian( u-> transition(u, x_, model ), s_)
F = ForwardDiff.jacobian( u-> transition(s_, u, model ), x_)

2×1 Matrix{Float64}:
 0.0
 1.0

1. 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.

# your code here

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

# your code here

1. 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.

# your code here

1. Check blanchard Kahn Conditions

# your code here

1. 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.

# your code here

1. Make some nice plots.