NLsolve.jl collocation divergence

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)

Are the equations you’re solving polynomials? If so you could try HomotopyContinuation.jl

And if it’s a low dimensional system you can try IntervalRootFinding.jl

Could you try to add the keyword inplace=false?

Hi, thank you for your response!

Unfortunately, it isn’t polynomial, it is a “freestyle” nonlinear problem. I am using Chebyshev polynomials to approximate the unknown function © in the (economics) Euler equation by forcing this approximate Euler equation to hold exactly at some pre-specified grid points. The approximation of C is a linear combination of Chebyshev polynomials, however, Euler equation is nonlinear in C, C also influences next periods capital stock…

How fast is IntervalRootFinding.jl? The dimensionality of my problem easily exceeds 100 equations, and I expect it to grow over time. :-/ Solution for 160 equations using BlackBoxOptim takes on my Intel Core i7-8750H about 36 minutes.

Výstřižek

Thank you for the suggestion!

Unfortunately, it didn’t work. I want to ask, is my function EulerCollocation, EulerCollResid correctly specified for NLsolve.jl? I worry, that I made some mistake there…

So I have used Chebyshev projection to solve these kinds of models but what I found was that doing it strictly by colocation was generally not very stable. As you say, unless the initial guess is very good, it can/will fail. There were two hybrid methods that I started to use. You can use Cheb. colocation but iteratively increase the order, that way when you reach higher orders you start with a very good approximation already (which helps convergence). Alternatively, and this is the one prefer for what I am doing, is that you can use Chebyshev approximations within the standard value function iteration. Convergence of VFI is very robust, and the Chebyshev helps speed up each iteration relative to other interpolation techniques.

I’m not sure how this could fit into your model if it truly is as high dimensional as you say though…

2 Likes

You are right, but this quess should be good enought, people are ushually solving this model using coefficients from linear approximation as initial guess (much worse than mine). I am literally starting at correct approximation, and when I try epsilon-deviation from this known solution, NLsove diverges. Actually, codes, that i saw in Matlab (but not my own code, off course) were solved simply by fsolve, so I hoped, that NLsolve.jl would be able to do this.

Thank you all, I managed to solve the problem using this reference. Parametrize NLsolve function

1 Like

Great. Also, I noticed that you have a bunch of let in your function definitions. Those are not necessary and could conceivably confuse things. Potentially harmless, and definetely not helpful to always use that pattern. Functions introduce their own scope already.

Thank you for this suggestion!