Perturbation of Neoclassical Model

Author

Pablo Winant

Our goal here is to compute a linear approximation of solution to the neoclassical model, close ot 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.

# 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 2 methods)

using ForwardDiff
x0 = [0.21, 1//2]
f(x0)

ForwardDiff.jacobian(f,x0)

2×2 Matrix{Float64}:
 1.0      1.0
 1.64872  0.346231

Create a NamedTuple to hold the model parameters.

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

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

typeof(model)

@NamedTuple{α::Float64, β::Float64, γ::Float64, δ::Float64, ρ::Float64}

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)

function transition(z, k, i, model)
    (;α, β, γ, δ, ρ) = model
    z_new = ρ*z
    k_new = (1-δ)*k + i
    return (z_new, k_new)
end

transition (generic function with 1 method)

function arbitrage(z, k, i, Z, K, I, model)

    (;α, β, γ, δ, ρ) = model

    y = exp(z)*k^α
    c = y - i
    Y = exp(Z)*K^α
    C = Y - I
    
    r = β*(C/c)^(-γ)*( (1-δ) + α*exp(Z)*K^(α-1) ) - 1

    return r

end

arbitrage (generic function with 2 methods)

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

function transition(s::Vector{T},x::Vector{U},model) where T<:Number where U<:Number
    r = transition(s[1],s[2],x[1],model)
    # return [r[1], r[2]]
    return [r...]
end

transition (generic function with 2 methods)

transition(m.s,m.x, model)

2-element Vector{Float64}:
 0.0
 2.920822149964071

function arbitrage(s,x,S,X,model)
    r = arbitrage(s[1],s[2],x[1],S[1],S[2],X[1],model)
    return [r]
end

arbitrage (generic function with 3 methods)

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(model)
    (;α, β, γ, δ, ρ) = model
    k = (( 1/β - (1-δ)) / α )^ (1/(α-1))
    z = 0.0
    i = δ*k
    return (;
        s = [z,k],
        x = [i],
    )
end

steady_state (generic function with 1 method)

m = steady_state(model)

r1 = arbitrage(m.s, m.x, m.s, m.x, model)
r2 = transition(m.s, m.x, model) - m.s
[r1,r2]

@assert maximum(abs, [r1; r2]) < 1e-12

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.

import ForwardDiff: jacobian

A = jacobian(u->arbitrage(u, m.x, m.s, m.x, model), m.s)
B = jacobian(u->arbitrage(m.s, u, m.s, m.x, model), m.x)
C = jacobian(u->arbitrage(m.s, m.x, u, m.x, model), m.s)
D = jacobian(u->arbitrage(m.s, m.x, m.s, u, model), m.x)
E = jacobian(u->transition(u, m.x, model), m.s)
F = jacobian(u->transition(m.s, u, model), m.x)

2×1 Matrix{Float64}:
 0.0
 1.0

(A = [5.074626865671642 0.5212190203776081], B = [-3.679193085018409;;], C = [-4.938626865671642 -0.5538125831185546], D = [3.679193085018409;;], E = [0.9 0.0; 0.0 0.9], F = [0.0; 1.0;;])

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.

function first_order_model(model)
    m = steady_state(model)
    A = jacobian(u->arbitrage(u, m.x, m.s, m.x, model), m.s)
    B = jacobian(u->arbitrage(m.s, u, m.s, m.x, model), m.x)
    C = jacobian(u->arbitrage(m.s, m.x, u, m.x, model), m.s)
    D = jacobian(u->arbitrage(m.s, m.x, m.s, u, model), m.x)
    E = jacobian(u->transition(u, m.x, model), m.s)
    F = jacobian(u->transition(m.s, u, model), m.x)
    pm = (;A, B, C, D, E, F)
    pm
end

first_order_model (generic function with 1 method)

first_order_model(merge(model, (;α=0.2)))

(A = [4.6575342465753415 0.6053173918792945], B = [-4.272828648559724;;], C = [-4.521534246575342 -0.6760184632507961], D = [4.272828648559724;;], E = [0.9 0.0; 0.0 0.9], F = [0.0; 1.0;;])

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.

residual(X::Array, M)  = M.A + M.B*X + (M.C+M.D*X)*(M.E+M.F*X)

residual (generic function with 1 method)

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

T(X, M) = - (B+(C+D*X)*F) \ (A+(C+D*X)*E)

T (generic function with 1 method)

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.

function linear_time_iteration(M; TT=200)
    nx, ns = size(M.A)
    X0 = rand(nx, ns)
    for t=1:TT

        X = T(X0, M)
        r = residual(X,M)
        ϵ= maximum(abs,r)
        println(t, ϵ)
        if ϵ<1e-10 
            return X
        end
        X0 = X
    end

    error("No convergence")




end

linear_time_iteration (generic function with 1 method)

Check blanchard Kahn Conditions

M = first_order_model(model)

(A = [5.074626865671642 0.5212190203776081], B = [-3.679193085018409;;], C = [-4.938626865671642 -0.5538125831185546], D = [3.679193085018409;;], E = [0.9 0.0; 0.0 0.9], F = [0.0; 1.0;;])

linear_time_iteration(M)

132.50570793314876
210.455512505432662
31.9466036686389936
40.9287841848793645
50.5970565928183014
60.4386101972004859
70.34512711663380324
80.28217302112220377
90.23593757247082392
100.1999610141445829
110.1708632158035548
120.14671240459370383
130.12631996092750652
140.10890566120863943
150.09392679583513797
160.08098577929560458
170.06977783639770285
180.060060235812492646
190.051633602213620566
200.04433022883245341
210.038006555509233486
220.03253817699504857
230.02781641354569553
240.023745859597579688
250.02024255316404
260.0172325456218263
270.014650735639557322
280.012439883089379045
290.010549751162809429
300.008936344979571409
310.007561227293772177
320.006390899327732402
330.005396239121463253
340.0045519922426331405
350.003836311007448856
360.00323033899607017
370.0027178379089787263
380.002284853889691796
390.0019194204454136
400.001611295093300047
410.001351726877772741
420.0011332519568290067
430.0009495145444757824
440.000795110618227568
450.0006654519482678367
460.0005566481706718029
470.0004654048037391334
480.0003889352874364427
490.0003248853059747425
500.00027126782837472163
510.00022640746855007166
520.00018889292281487613
530.00015753638752835641
540.00013133899222594891
550.00010946140386947079
569.119886607455996e-5
577.596003383314454e-5
586.324905002186298e-5
595.2650385693464585e-5
604.381603262348932e-5
613.645469470381357e-5
623.0322675377281172e-5
632.521620221651588e-5
642.096496770853662e-5
651.742669806636954e-5
661.4482590021458464e-5
671.2033479718898121e-5
689.996628493436077e-6
698.303027978939781e-6
706.895142039997637e-6
715.725015829405322e-6
724.752693176346412e-6
733.944892712937076e-6
743.2739010031868077e-6
752.7166475389250877e-6
762.253932073426057e-6
771.8697794725852646e-6
781.550901268299043e-6
791.2862464315865907e-6
801.0666267296066678e-6
818.844043821731873e-7
827.332317570885039e-7
836.078344969573379e-7
845.038308832361338e-7
854.1758143876080567e-7
863.4606372745216163e-7
872.867681665463806e-7
882.3761134482214175e-7
891.9686392294104849e-7
901.6309067474296057e-7
911.3510063912036685e-7
921.1190568782737387e-7
939.268608858192806e-8
947.6761896394828e-8
956.356918635930242e-8
965.2640311665186346e-8
974.35875118220963e-8
983.6089302657416056e-8
992.9879175134084335e-8
1002.4736217962839646e-8
1012.047733405063923e-8
1021.6950784331726254e-8
1031.4030824946331677e-8
1041.1613267236754155e-8
1059.611787366026192e-9
1067.954870984150375e-9
1076.583275702354285e-9
1085.447930107038701e-9
1094.508190265539724e-9
1103.730394215750721e-9
1113.086666477969402e-9
1122.553921962800132e-9
1132.1130466265617542e-9
1141.7482140179936323e-9
1151.4463208408699302e-9
1161.1965197721508503e-9
1179.898299957455947e-10
1188.188179023704834e-10
1196.773301919338337e-10
1205.602740493770852e-10
1214.634341799203412e-10
1223.833213746418096e-10
1233.1704905367746505e-10
1242.6222757298910437e-10
1252.1687940332526523e-10
1261.793694082152797e-10
1271.4834311556910507e-10
1281.2268097648870935e-10
1291.014557327039256e-10
1308.390088623855263e-11

1×2 Matrix{Float64}:
 0.768674  0.0278097

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.

Other Formats