Hi everyone!
I would like to ask for help/suggestions/advice with my code. I am trying to solve a simple macroeconomic model (plain vanilla deterministic optimal growth) using collocation with Chebyshev polynomials. To solve collocation equations, I tried to use NLsolve.jl (trust-region and newton algorithms). As an initial guess, I used coefficients obtained by least-squares fitting polynomials to closed-form solutions of the model.
After checking that these coefficients are correct, I used them as a guess for NLsolve. However, solver diverged even for the smallest “homotopy” from known solution (it “converged” only when it started at correct solution), regardless of used algorithm or number of iterations. Since I am pretty sure that my residual function/collocation system is correct (I solved model by using BlackBoxOptimize.jl to find zero of squared residuals, it works robustly, but pretty slowly), I would like to ask whether I am using NLsolve.jl correctly, or there is some error in my understanding, how NLsolve works.
I would also be glad for any alternative suggestions, how to solve nonlinear equations in Julia (“my” solution through BlackBoxOptimize is robust but very slow).
Thanks in advance!
Honza
#*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
#-----Solve simple deterministic growth model using ChebyshevProjection-----
#*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
#(1) Install and initialize packages
using Pkg
Pkg.add("NLsolve")
Pkg.add("Distributed")
Pkg.add("Plots")
Pkg.add("LinearAlgebra")
Pkg.add("Random")
Pkg.add("Distributions")
Pkg.add("Parameters")
Pkg.add("Optim")
using NLsolve
using Distributed
using Parameters
using LinearAlgebra
using Random
using Distributions
using Plots
using Optim
#(2) Include basic functions
function Transform(y,a,b)
let
x = 2 .*((y.-a)./(b.-a)) .-1
return x
end
end
function iTransform(x,a,b)
let
y = a .+ (b.-a).*(x.+1)./2
return y
end
end
function ChebyZeroBase(n)
let
Zrs = Array{Float64}(undef,n,1)
for i = 1:n
Zrs[i] = -cos(pi*(2*i-1)/(2*n))
end
return Zrs
end
end
function Cheby(n,x)
let
m = n+1
if(n<0)
println("Error, please set order to n >=0")
return
elseif(n==0)
y = 1.0
return y
elseif(n==1)
y = x
return y
elseif(n>=2)
T = Array{Float64}(undef,m,1)
T[1,1] = 1.0
T[2,1] = x
for i = 3:m
T[i,1] = 2*x*T[i-1,1] - T[i-2,1]
end
y = T[m,1]
return y
end
end
end
function ChebyDir(n,x)
let
y = cos.(n.*acos(x))
return y
end
end
function BuildX(Param,n,ζ)
let
@unpack a,b = Param
X = Array{Float64}(undef,length(ζ),n+1)
𝞯 = Transform(ζ,a,b)
for i in 1:length(𝞯),j in 1:(n+1)
X[i,j] = Cheby(j-1,𝞯[i])
end
return X
end
end
function NGsteady(Econ)
@unpack α,β,δ = Econ
Kss = (α/(1/β - (1-δ)))^(1/(1-α))
Css = Kss^α - δ*Kss
return Kss,Css
end
function BMPolicy(Econ,K,A=nothing)
let
@unpack α,β = Econ
if(typeof(K)!=Array{Tuple{Float64,Float64},4})
if(A==nothing)
A = 1.0
C = (1-α*β).*A.*K.^α
else
C = (1-α*β).*A.*K.^α
end
return C
elseif(typeof(K)==Array{Tuple{Float64,Float64},4})
C = Array{Float64}(undef,length(K),1)
for i in 1:length(K)
C[i] = (1-α*β)*(K[i][2])*((K[i][1])^α)
end
return C
end
end
end
function ChebyApIn(Param,θ,y)
let
@unpack a,b = Param
ϕ = Array{Float64}(undef,length(θ),1)
for i in 1:length(θ)
ϕ[i] = θ[i]*Cheby(i-1,Transform(y,a,b))
end
π = sum(ϕ)
return π
end
end
function ChebyApOt(Param,θ,ζ)
let
φ = Array{Float64}(undef,length(ζ),1)
for i in 1:length(ζ)
φ[i] = ChebyApIn(Param,θ,ζ[i])
end
return φ
end
end
#(3) Calibrate parameter values
@with_kw struct Param
a::Float64 = 0.5-1.5^(1/3)
b::Float64 = 0.5+2^(1/3)
end
@with_kw struct Econ
α::Float64 = 1/3
β::Float64 = 1/(1+0.05)
γ::Float64 = 1
δ::Float64 = 0.85
end
Econ1 = Econ(δ=1.0)
Kss = NGsteady(Econ1)[2]
Param2 = Param(a=Kss*0.25,b=Kss*5.5)
@unpack a,b = Param2
τ = 15
#(4) Build grid
ζ = iTransform(ChebyZeroBase(500),a,b)
kGrid = iTransform(ChebyZeroBase(τ+1),a,b)
#(5) Prepare initial guess
Y = log.(BMPolicy(Econ1,ζ))
X = BuildX(Param2,τ,ζ)
θ = inv(transpose(X)*X)*transpose(X)*Y
ψ = exp.(ChebyApOt(Param2,θ,ζ))
display(Plots.plot(ζ,exp.(Y),title="Policy function computed by ChebyshevProjection",
xlabel="K(t)",ylabel="C(t)",legend=:bottomright,label="BrockMirman"))
display(Plots.plot!(ζ,ψ,label="LeastSquares"))
#(6) Prepare residual function
function EulerCollocation(θ,Econ,Param,Grid)
let
@unpack α,β,γ,δ = Econ
@unpack a,b = Param
Eu = Array{Float64}(undef,length(Grid),1)
for i in 1:length(Grid)
Ct = exp.(ChebyApIn(Param,θ,Grid[i]))
Kt1 = (1-δ)*Grid[i] + Grid[i]^α - Ct
Ct1 = exp.(ChebyApIn(Param,θ,Kt1))
Eu[i] = (Ct^(-γ) - β*(Ct1^(-γ))*(α*(Kt1^(α-1))+(1-δ)))
end
return Eu
end
end
#(7) Solve it using NLsolve.jl
Econ2 = Econ(δ=0.99)
EulerCollResid(θ) = EulerCollocation(θ,Econ1,Param2,kGrid)
Sol = nlsolve(EulerCollResid,θ,ftol=10^(-14))
converged(Sol)
θ2 = Sol.zero
ψ2 = exp.(ChebyApOt(Param2,θ2,ζ))
display(Plots.plot!(ζ,ψ2,label="Collocation"))
EulerCollResid(θ) = EulerCollocation(θ2,Econ2,Param2,kGrid)
Sol2 = nlsolve(EulerCollResid,θ2,ftol=10^(-12),iterations=500000)
converged(Sol2)
Econ3 = Econ(γ=1.2)
EulerCollResid(θ) = EulerCollocation(θ2,Econ3,Param2,kGrid)
Sol3 = nlsolve(EulerCollResid,θ2,ftol=10^(-12),iterations=500000)
converged(Sol3)
