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::Vector{<:Number}) # equivalent
    a = x[1]
    b = x[2]
    x1 = a+b
    x2 = a*exp(b)
    return [x1,x2]
end

f (generic function with 4 methods)

using ForwardDiff

ForwardDiff.jacobian(f, [0.2, 0.4])

2×2 Matrix{Float64}:
 1.0      1.0
 1.49182  0.298365

Create a NamedTuple to hold the model parameters.

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

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

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::Number, k::Number, i::Number, p)

        Z = p.ρ * z
        K = (1-p.δ) * k + i

        (;Z, K)

end

transition (generic function with 1 method)

function arbitrage(z::Number, k::Number, i::Number, Z::Number, K::Number, I::Number, p)

    # positional unpacking (error-prone)
    # α, β, γ, δ, ρ = p
    
    (;α, β, γ, δ, ρ) = p

    # define auxiliary variables today
    y = exp(z)k^α
    c = y - i

    # define auxiliary variables tomorrow
    Y = exp(Z)K^α
    C = Y - I

    residual = β*(C/c)^(-γ)*( (1-δ) + α*K^(α-1)*exp(Z)) - 1

    return residual


end

arbitrage (generic function with 1 method)

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

# this returns a number
# arbitrage(s::Vector{T}, x::Vector{T}, S::Vector{T}, X::Vector{T}, p) where T<:Number = arbitrage(s[1],s[2],x[1],S[1],S[2],X[1],p)

[2.4]   # create a vector from  a number

1-element Vector{Float64}:
 2.4

arbitrage(s::Vector{T}, x::Vector{T}, S::Vector{T}, X::Vector{T}, p) where T<:Number = [
    arbitrage(s[1],s[2],x[1],S[1],S[2],X[1],p)
]

arbitrage (generic function with 2 methods)

# this returns a tuple
# transition(s::Vector{T}, x::Vector{T}, p) where T<: Number = transition(s[1], s[2], x[1], p)

t = (1,2,3)

(1, 2, 3)

# to convert into a tuple into a vector
[t...]

3-element Vector{Int64}:
 1
 2
 3

transition(s::Vector{T}, x::Vector{T}, p) where T<: Number = [transition(s[1], s[2], x[1], p)...]

transition (generic function with 2 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(p)
    (;α, β, γ, δ, ρ) = p

    # ...
    z = 0.0

    k = ((1/β - (1-δ))/α)^ (1/(α-1))
    i = δ*k

    s = [z,k] # vector of states
    x = [i]  # vector controls

    return (;
        s,
        x
    )


end

steady_state (generic function with 1 method)

q = steady_state(parameters)

(s = [0.0, 2.920822149964071], x = [0.29208221499640713])

# check the steady-state is correct using the functions representing the model

methods(transition)

# 2 methods for generic function transition from [35mMain[39m:

transition(z::Number, k::Number, i::Number, p) in Main at /home/pablo/Teaching/ensae/mie37/tutorials/3_perturbation_neoclassical.ipynb:1
transition(s::Vector{T}, x::Vector{T}, p) where T<:Number in Main at /home/pablo/Teaching/ensae/mie37/tutorials/3_perturbation_neoclassical.ipynb:1

@assert maximum(transition(q.s, q.x, parameters) - q.s) == 0.0

@assert maximum(arbitrage(q.s, q.x, q.s, q.x, parameters)) == 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}\]

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

struct PerturbedModel
    A::Matrix
    B::Matrix
    C::Matrix
    D::Matrix
    E::Matrix
    F::Matrix

end

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 need to loosen the constraint on the arbitrage arguments:

# brutal
arbitrage(s, x, S, X, p) where T<:Number = [
    arbitrage(s[1],s[2],x[1],S[1],S[2],X[1],p)
]


# more precise
# arbitrage(s::Vector{<:Number}, x::Vector{<:Number}, S::Vector{<:Number}, X::Vector{<:Number}, p) = [
#     arbitrage(s[1],s[2],x[1],S[1],S[2],X[1],p)
# ]

WARNING: method definition for arbitrage at /home/pablo/Teaching/ensae/mie37/tutorials/3_perturbation_neoclassical.ipynb:4 declares type variable T but does not use it.

arbitrage (generic function with 3 methods)

arbitrage(s,x,s,x,parameters)

UndefVarError: UndefVarError: `s` not defined

transition(s, x, p) where T<: Number = [transition(s[1], s[2], x[1], p)...]

WARNING: method definition for transition at /home/pablo/Teaching/ensae/mie37/tutorials/3_perturbation_neoclassical.ipynb:1 declares type variable T but does not use it.

transition (generic function with 3 methods)

using ForwardDiff

function first_order_model(s, x, parameters)

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

    return PerturbedModel(A,B,C,D,E,F)

end

first_order_model (generic function with 1 method)

@time model = first_order_model(q.s, q.x, parameters)

  1.343581 seconds (7.27 M allocations: 474.892 MiB, 7.63% gc time, 99.95% compilation time)

PerturbedModel([5.074626865671642 0.5212190203776081], [-3.679193085018409;;], [-4.938626865671642 -0.5538125831185546], [3.679193085018409;;], [0.9 0.0; 0.0 0.9], [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.

function residual(X::Matrix, M::PerturbedModel)
    (;A,B,C,D,E,F) = M # keyword unpacking
    return A + B*X + (C+D*X)*(E+F*X)
end

residual (generic function with 1 method)

X0 = zeros(1, 2)
residual(X0, model)

1×2 Matrix{Float64}:
 0.629863  0.0227877

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

function T(X::Matrix, M::PerturbedModel)
    (;A,B,C,D,E,F) = M # keyword unpacking

    C_DX = (C+D*X)

    return -(B+C_DX*F) \ (A+C_DX*E)
    # return  -solve(B+C_DX*F ,   (A+C_DX*E))
    # return  -inv(B+C_DX*F)*(A+C_DX*E)) # not ok
end

T (generic function with 1 method)

T(X0, model)

1×2 Matrix{Float64}:
 0.148798  0.00538334

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.

A = rand(10000, 10000);
B = rand(10000, 10000);

distance(A,B) = sum(abs(e1-e2) for (e1, e2) in zip(A,B))

distance (generic function with 1 method)

function linear_time_iteration(X_0, M; N=5, τ_η=1e-8, verbose=true)

    η_0 = 1.0

    for n in 1:N

        X = T(X_0, M)

        # successive approximation error
        η = distance(X, X_0)


        # ratio of successive approximation errors
        λ = η/η_0

        if verbose
            println(n, " : ", η, " : ", λ)
        end

        # η_0 will be the value from last iteration
        η_0 = η

        if η<τ_η
            return X
        end

        X_0 = X
    end

    error("No convergence")

end

linear_time_iteration (generic function with 1 method)

@time sol = linear_time_iteration(X0, model; N=500)

1 : 0.15418131543048902 : 0.15418131543048902
2 : 0.12190031209691013 : 0.7906296022741327
3 : 0.0971923112531961 : 0.7973097819136732
4 : 0.07801352484927856 : 0.802671773552593
5 : 0.06295682833902616 : 0.8069988948795506
6 : 0.05102689044058986 : 0.8105060529067174
7 : 0.04150316170984929 : 0.8133586301554283
8 : 0.033853524549360574 : 0.815685435871905
9 : 0.027678228788813052 : 0.8175878038476156
10 : 0.022672513753932016 : 0.8191461211960124
11 : 0.018601088097577975 : 0.8204246030880477
12 : 0.01528032593947912 : 0.8214748438006029
13 : 0.01256560026933672 : 0.8223384971698489
14 : 0.010342108838667052 : 0.8230493264937326
15 : 0.008518120644127804 : 0.8236347902547955
16 : 0.007019930424960394 : 0.8241172810577381
17 : 0.005788038650448812 : 0.8245151020113518
18 : 0.004774224564889556 : 0.8248432419364298
19 : 0.0039392795027168095 : 0.8251139947808338
20 : 0.003251234928149705 : 0.8253374572450166
21 : 0.0026839657341018452 : 0.8255219304097795
22 : 0.0022160813815156506 : 0.8256742451509851
23 : 0.0018300400639542705 : 0.8258000266680848
24 : 0.001511437244487067 : 0.8259039101150711
25 : 0.0012484316214957068 : 0.8259897167741055
26 : 0.001031280172156878 : 0.8260605982739636
27 : 0.0008519603042240476 : 0.8261191548386015
28 : 0.0007038619422192541 : 0.8261675323715003
29 : 0.000581536017211115 : 0.826207502251863
30 : 0.00048048862532256223 : 0.8262405269872226
31 : 0.00039701228622349843 : 0.8262678142629821
32 : 0.0003280474254938187 : 0.8262903614754736
33 : 0.0002710685376324101 : 0.8263089924402341
34 : 0.00022399054337329805 : 0.8263243876611255
35 : 0.00018509169812225973 : 0.8263371093028231
36 : 0.00015295008455855352 : 0.8263476218016245
37 : 0.00012639126731726208 : 0.826356308870663
38 : 0.00010444512845603215 : 0.826363487549012
39 : 8.631026020181412e-5 : 0.8263694197872312
40 : 7.132458275916145e-5 : 0.8263743220375822
41 : 5.894109266629463e-5 : 0.8263783731524711
42 : 4.870784159012598e-5 : 0.8263817209140318
43 : 4.0251404707945804e-5 : 0.8263844874642432
44 : 3.3263228473082856e-5 : 0.8263867736898273
45 : 2.7488354905323004e-5 : 0.8263886630116805
46 : 2.2716107776505484e-5 : 0.826390224323923
47 : 1.8772398710901333e-5 : 0.8263915145849502
48 : 1.551337101909217e-5 : 0.8263925808311001
49 : 1.2820148383725638e-5 : 0.826393461997911
50 : 1.0594496141231452e-5 : 0.8263941901546545
51 : 8.755236434329627e-6 : 0.8263947919388228
52 : 7.235286145202552e-6 : 0.8263952892046078
53 : 5.979209359931875e-6 : 0.8263957001751017
54 : 4.94119493606518e-6 : 0.8263960397803294
55 : 4.083385313793236e-6 : 0.8263963204505664
56 : 3.374495545591105e-6 : 0.8263965524371217
57 : 2.788672131517206e-6 : 0.8263967439994704
58 : 2.304550011482387e-6 : 0.826396902467187
59 : 1.9044732927631092e-6 : 0.8263970333792274
60 : 1.5738512852726816e-6 : 0.8263971415368372
61 : 1.300626344003658e-6 : 0.8263972309037538
62 : 1.074834105267708e-6 : 0.8263973048239941
63 : 8.882400730839124e-7 : 0.826397365631303
64 : 7.340393016162794e-7 : 0.8263974164864485
65 : 6.066082128101946e-7 : 0.8263974578397988
66 : 5.012995060102443e-7 : 0.8263974925230678
67 : 4.1427266909629945e-7 : 0.826397521101554
68 : 3.4235391643128166e-7 : 0.8263975443470543
69 : 2.8292044249381143e-7 : 0.8263975637930232
70 : 2.338047691008427e-7 : 0.8263975803231571
71 : 1.9321569839914066e-7 : 0.8263975929242252
72 : 1.596729903574623e-7 : 0.8263976047516254
73 : 1.3195337798116435e-7 : 0.826397612305991
74 : 1.0904595768024272e-7 : 0.8263976212553532
75 : 9.011532065353323e-8 : 0.8263976269324892
76 : 7.44710876507404e-8 : 0.826397632618539
77 : 6.154273075875683e-8 : 0.8263976356486713
78 : 5.085876750646201e-8 : 0.8263976407843324
79 : 4.202956567153637e-8 : 0.8263976445397774
80 : 3.473313383434151e-8 : 0.8263976388855188
81 : 2.870338035310116e-8 : 0.8263976550460705
82 : 2.372040580275736e-8 : 0.8263976406595807
83 : 1.9602487667802482e-8 : 0.8263976523337474
84 : 1.619944979894261e-8 : 0.8263976528627315
85 : 1.33871871400848e-8 : 0.8263976435149435
86 : 1.1063140292755236e-8 : 0.8263976724153838
87 : 9.142552898772083e-9 : 0.8263976282357315
  0.188574 seconds (380.90 k allocations: 25.696 MiB, 5.04% gc time, 99.17% compilation time)

1×2 Matrix{Float64}:
 0.768674  0.0278097

residual(sol, model)

1×2 Matrix{Float64}:
 3.01193e-8  1.08968e-9

Check blanchard Kahn Conditions

# check that solution is not diverging
# let's check the transition matrix has spectral radius smaller than 1

using LinearAlgebra

function bk_check(X, M)
    (;A,B,C,D,E,F) = M # keyword unpacking

    P = E+F*X
    maximum(abs,eigvals(P))

    # can check the solution is unique by computing the first "rejected" eigenvalue


end

bk_check (generic function with 1 method)

bk_check(sol, model)

4.13068831324392

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.

P = model.E + model.F*sol

2×2 Matrix{Float64}:
 0.9       0.0
 0.768674  0.92781

Make some nice plots.

Other Formats