Model solving - method error using NLsolve

Hi, I am stuck with a code (solving an economic model) and would need some help. I am sorry if the mistake is obvious, but I am new to this and would really need some help there. I have the following code:

``````function bk_opti_decent(grid::Any,
β::AbstractFloat,
Re::AbstractFloat,
Rl::AbstractFloat,
u::Function,
u_prime::Function,
E::AbstractFloat)

function foc1(Sb::Number,
c_bar::Number,
Lb::Number,
K::Number,
Lf::Number,
Sf::Number,
θhigh::Number,
P_star::Number,
θlow::Number,
D::Number)

c1_ND = D * c_bar
c1_D = Lb + Lf
c21(θ) = (Lb - θ* c_bar * D) * Rl + Sb * Re + Lf * Rl + Sf * Re
c22(θ) = (Lb - θ* c_bar * D) * Re + Sb * Re + Lf * Rl + Sf * Re
c2_D(θ) =  (Sb + Sf) * Re / (1 - θ)
E_UD = θhigh * u(Lb + Lf) + (1-θhigh) * u(Lb + Lf + Re*(Sb + Sf)/ (1-θhigh))
E_UND = θhigh * u(c1_ND) + (1-θhigh) * u( Re*(Sb + Sf)/ (1-θhigh))
Diff = E_UND - E_UD

val1(θ) = (Re - 1) * u_prime(c21(θ))
val2(θ) = (Re - Rl) * u_prime(c22(θ))
val3(θ) = (P_star-1) * θ * u_prime(c1_D) + ((P_star - 1)*(1- θ) + Re ) * u_prime(c2_D(θ))
val4 = Diff * (P_star - 1) * (Re / ( c_bar * D * (Re - P_star)))

end

Bk_policyfunc = Array{Real}(0,2) #solution matrix

for (i,D) in enumerate(grid)

function f!(F::Any ,x::Any)

F[1]= foc1(x[1], x[2],x[3],x[4],x[5],x[6],x[7],x[8],x[9], D)
F[2]= x[3] - D - x[1]
F[3]= x[4] - E - D
F[4]= x[5] - x[2] * (Rl/(Re-Rl))
F[5]= x[6] - x[4] + x[5]
F[6]= x[7] -  x[5] / x[1]
F[7]= x[9] - x[3]/ (D * x[2])
F[8]= x[8] - (Re * x[3] - Re * x[1] * x[7] + x[2] * x[7] * D) / ( x[2] * D * (Re - x[7]))

end

initial_x = [0.5 0.5 0.5 0.5 0.5 0.5 0.8 0.3 0.3 ]
results1_bk = nlsolve(f!, initial_x)

Bk_policyfunc[1, i] = results.zero[1]
Bk_policyfunc[2, i] = results.zero[2]
Bk_policyfunc[3, i] = results.zero[3]
Bk_policyfunc[4, i] = results.zero[4]
Bk_policyfunc[5, i] = results.zero[5]
Bk_policyfunc[6, i] = results.zero[6]
Bk_policyfunc[7, i] = results.zero[7]
Bk_policyfunc[8, i] = results.zero[8]

return Bk_policyfunc

end

end
``````

Then I create a structure for a model (m) with specific parameters and when i tried solving part of it with

``````bank = bk_opti_decent(m.grid,
m.β,
m.Re,
m.Rl,
m.u,
m.u_prime,
m.E)
``````

I get the following error that I am not able to solve. I changed the types of my unknow to “Number”, but it did not solve the problem. Here is the error:

``````MethodError: no method matching //(::#c21#289{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64}, ::Int32)
Closest candidates are:
//(::Integer, ::Integer) at rational.jl:40
//(::Rational, ::Integer) at rational.jl:43
//(::Complex, ::Real) at rational.jl:56
...

Stacktrace:
[4] (::#foc1#288{Float64,Float64,##323#327,##324#328})(::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64) at .\In[38]:47
[5] f! at .\In[38]:93 [inlined]
[6] (::NLsolve.#fj!#1{#f!#302{Float64,Float64,Float64,#foc1#288{Float64,Float64,##323#327,##324#328},#foc2#295{Float64,Float64,##323#327,##324#328}}})(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\objectives\autodiff.jl:5
[7] value_jacobian!!(::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLSolversBase\src\interface.jl:124
[8] trust_region_(::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,2}}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\solvers\trust_region.jl:118
[9] #nlsolve#38(::Symbol, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int32, ::Float64, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\nlsolve\nlsolve.jl:26
[10] (::NLsolve.#kw##nlsolve)(::Array{Any,1}, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}) at .\<missing>:0
[11] #nlsolve#39(::Symbol, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int32, ::Float64, ::Symbol, ::Bool, ::NLsolve.#nlsolve, ::#f!#302{Float64,Float64,Float64,#foc1#288{Float64,Float64,##323#327,##324#328},#foc2#295{Float64,Float64,##323#327,##324#328}}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\nlsolve\nlsolve.jl:59
[12] bk_opti_decent(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::##323#327, ::##324#328, ::Float64) at .\In[38]:106

``````

Many thanks.

Is `u_prime(c21)` in the line `val1(θ) = (Re - 1) * u_prime(c21)` the derivative of `c21` with respect to `θ`?

EDIT: And you may want to keep in mind that `quadgk` returns a tuple with the estimate of the integral and its error, so you should do `return quadgk(val1, 0, θlow)[1] + ...`

1 Like

Thanks a lot! I did not know for quadgk, i will correct.

`u_prime(c21)` is not the derivative with respect to `θ` (sorry I did not give details about the problem I am trying to solve), i am writing the first order condition of the model by piece, so that it is the derivative of the utility function evaluated at c21 (which itself depends on `θ`).

the error is indeed about this c21 which is a function of `θ` plus of my unkown. I do not understand what is the method error, could it come from the fact that with the nested functions i did not correctly specified the function? (i tried many ways and fail each time so far).

thanks so much for your time.

Is there any chance you can post `u_prime`?

1 Like

yes, sure, it is defined in the structure of the model below: (thanks!) (just in case you wonder, I posted only the part of the model where the error occurs)

``````
struct Model{TF <: AbstractFloat, TR <: Real, TI <: Integer}
β::TF
γ::TR
grid_min::TR
grid_max::TR
E::TR
grid_size::TI
u::Function
u_prime::Function
Re::TF
Rl::TF
grid::Vector{TR}
end

function Model(; β=0.8,
γ=3.0,
grid_min=0.0,
grid_max=4.0,
E=4.0,
grid_size=200,
u::Function= c-> c^γ,
u_prime::Function= c-> γ*((c//1)^(γ-1)),
Re=4.0,
Rl=3.0)

grid = collect(linspace(grid_min, grid_max, grid_size))

m= Model(β, γ, grid_min, grid_max, E, grid_size, u, u_prime, Re, Rl, grid)

return m
end
``````

and then i defined a specific instance where the utility function is CRRA and `u_prime` its derivative

``````m = Model(γ =3.0)
``````

Don’t know if that is the problem, but NLsolve does not support Rationals I think (the c//1 part in u_prime).

1 Like

You probably want `val1(θ) = (Re - 1) * u_prime(c21(θ))`: you forgot to pass the parameter `θ`.

Yes sorry i had corrected the error above but the method error remains. (The `val1(θ) = (Re - 1) * u_prime(c21(θ))` )

About the nlsolve not supporting rationals, do you have any suggestions ?

Remove it? I mean, what is that `c//1` for anyway?

edit: read my own response and realized that it might be read as if I have an attitude. I don’t - let me prove it with a smiley

In Julia `2//1` creates a special “rational” number type. If you just want to do “real”/floating point division, it’s `2/1` or `c/1` in your case.

``````julia> 1//2
1//2

julia> 1/2
0.5
``````
1 Like

OK so that `c//1` seems weird but i put it at some point due to a preceeding error: i had negative number in the exponent `^(γ-1)` (for other values of the parameters than the ones above) and found this solution. The problem was something along those lines : Exponentiation by a negative power throws a DomainError · Issue #3024 · JuliaLang/julia · GitHub
Is there a better solution? 'ill go through all the post again, probably i picked a stupid solution…

How sure are you that you didn’t at some point try to evaluate for integer CRRA and integer consumption? Because that shouldn’t be a problem as long as c and/or gamma-1 are floats.

2 Likes

So -checking if i understood correctly (sorry if not)- if c and gamma are integers then a negative exponent would be a problem? Whereas if they are floats there is no problem? Is that a common problem, would it happen in matlab for instance?

edit: i tried the solutions you implied (as I understood), and I am making progress. thanks!!

1. for `gamma`: i defined it as follows so it can be an integer at this point?
``````struct Model{TF <: AbstractFloat, TR <: Real, TI <: Integer}
β::TF
γ::TR
grid_min::TR
grid_max::TR
E::TR
grid_size::TI
u::Function
u_prime::Function
Re::TF
Rl::TF
grid::Vector{TR}
end
``````

(once the code works (…) i want to loop over different values of gamma to see how risk aversion affects the model)

1. for `c` : I am not sure all that all `c` (all the possible consumptions) could not be integers. The possible consumptions are defined as follows and a potential problem is that `c_bar` is a choice variable here
``````c1_ND = D * c_bar
c1_D = Lb + Lf
c21(θ) = (Lb - θ* c_bar * D) * Rl + Sb * Re + Lf * Rl + Sf * Re
c22(θ) = (Lb - θ* c_bar * D) * Re + Sb * Re + Lf * Rl + Sf * Re
c2_D(θ) =  (Sb + Sf) * Re / (1 - θ)
``````

Would a way of making sure that cannot be integers to change that

``````      function foc1(Sb::Number,
c_bar::Number,
Lb::Number,
K::Number,
Lf::Number,
Sf::Number,
θhigh::Number,
P_star::Number,
θlow::Number,
D::Number)
``````

replacing number with float. I had changed to Number to try solving the latest (of the many) method error I got.

Sorry for being so ignorant. I am really confused so far with the types problem.

and thanks a lot for your time!

You have a comma in your return statement: `return quadgk(val1, 0, θlow)[1], + quadgk(val2, θlow, θhigh)[1], + quadgk(val3, θhigh, 1)[1], + val4`. It wasn’t there the first time you posted (but also it makes no difference because the error is elsewhere).

You may also be able to proceed a bit further (and get a bit more help) if you could make a minimal working example of your code: since the error is localized within the call to `quadgk`, you could for instance leave out the definition of `val2` and the corresponding term `quadgk(val2, θlow, θhigh)[1]`, then simplify the definition of `c21(x) = x` and so on and remove all unnecessary terms until you figure it out.

1 Like

Edit: reading ForwardDiff documentation, i wonder if it is coming from the fact that my function has more than one argument. However I had alos several arguments in the first example below which is working. should i maybe use OnceDifferentiable manually to compute the jacobian and take it as an input to NLsolve?

Thanks a lot. I was able to solve the method error I had when i first posted thanks to all your comments. It was indeed linked to the exponential in the utility function and the negative number and incorrect types - which seems solved.

I now have another one (but the code goes further so it is progress). I used simplified code and try isolate the error as you suggested. The error still comes from the `quadgk` but is different from before as it is now able to evaluate it.

This one shows no method error:

``````function bk_opti_decent(grid::Any,
β::AbstractFloat,
Re::AbstractFloat,
Rl::AbstractFloat,
u::Function,
u_prime::Function,
E::AbstractFloat)

function foc1(Sb::AbstractFloat,
c_bar::AbstractFloat,
Lb::AbstractFloat,
K::AbstractFloat,
Lf::AbstractFloat,
Sf::AbstractFloat,
θhigh::AbstractFloat,
P_star::AbstractFloat,
θlow::AbstractFloat,
D::AbstractFloat)

c1_ND = D * c_bar + Lb + K + Lf + Sf + θhigh + θlow + P_star + Sb

return  c1_ND

end #end of foc1

for (i,D) in enumerate(grid)

function f!(F ,x)
F[1]= foc1(x[1], x[2],x[3],x[4],x[5],x[6],x[7],x[8],x[9], D)
F[2]= x[2] - E

end #end of function f!

initial_x = [0.5 0.5 ]
essai = nlsolve(f!, initial_x, autodiff = :forward)

end #end of loop

end

bank = bk_opti_decent(m.grid,
m.β,
m.Re,
m.Rl,
m.u,
m.u_prime,
m.E)
``````

But putting back the `quadgk` as follows does not work generating a method error in the function foc1 (see below). As from example just above which is working, it seems that the type in the two functions are correctly defined, I am not sure where the error comes from now. The “closest candidates” suggested in the error are exactly the types I have for the arguments of the function foc1. Why evaluating the integral would change the types of foc1? I am for sure completely missing something.

``````function bk_opti_decent(grid::Any,
β::AbstractFloat,
Re::AbstractFloat,
Rl::AbstractFloat,
u::Function,
u_prime::Function,
E::AbstractFloat)

function foc1(Sb::AbstractFloat,
c_bar::AbstractFloat,
Lb::AbstractFloat,
K::AbstractFloat,
Lf::AbstractFloat,
Sf::AbstractFloat,
θhigh::AbstractFloat,
P_star::AbstractFloat,
θlow::AbstractFloat,
D::AbstractFloat)

c1(θ) = D * c_bar + Lb + K + Lf + Sf + θhigh + θlow + P_star + Sb * θ

end #end of foc1

for (i,D) in enumerate(grid)

function f!(F ,x)
F[1]= foc1(x[1], x[2],x[3],x[4],x[5],x[6],x[7],x[8],x[9], D)
F[2]= x[2] - E

end #end of function f!

initial_x = [0.5 0.5 0.5 0.5 0.5 0.5 0.8 0.3 0.3 ]
essai = nlsolve(f!, initial_x, autodiff = :forward)

end #end of loop

end

bank = bk_opti_decent(m.grid,
m.β,
m.Re,
m.Rl,
m.u,
m.u_prime,
m.E)
``````

with the following error:

``````MethodError: no method matching (::#foc1#220)(::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9}, ::Float64)
Closest candidates are:
foc1(::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat, ::AbstractFloat) at In[37]:22

Stacktrace:
[1] (::#f!#222{Float64,#foc1#220})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2}) at .\In[37]:31
[2] vector_mode_dual_eval(::#f!#222{Float64,#foc1#220}, ::Array{Float64,2}, ::Array{Float64,2}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2}}}) at C:\Users\arquie\.julia\v0.6\ForwardDiff\src\apiutils.jl:42
[3] vector_mode_jacobian!(::DiffResults.MutableDiffResult{1,Array{Float64,2},Tuple{Array{Float64,2}}}, ::#f!#222{Float64,#foc1#220}, ::Array{Float64,2}, ::Array{Float64,2}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2}}}) at C:\Users\arquie\.julia\v0.6\ForwardDiff\src\jacobian.jl:163
[4] jacobian!(::DiffResults.MutableDiffResult{1,Array{Float64,2},Tuple{Array{Float64,2}}}, ::Function, ::Array{Float64,2}, ::Array{Float64,2}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{#f!#222{Float64,#foc1#220},Float64},Float64,9},2}}}, ::Val{false}) at C:\Users\arquie\.julia\v0.6\ForwardDiff\src\jacobian.jl:74
[5] (::NLsolve.#fg!#4{#f!#222{Float64,#foc1#220}})(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\objectives\autodiff.jl:24
[6] value_jacobian!!(::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLSolversBase\src\interface.jl:124
[7] trust_region_(::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,2}}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\solvers\trust_region.jl:118
[8] #nlsolve#38(::Symbol, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int32, ::Float64, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\nlsolve\nlsolve.jl:26
[9] (::NLsolve.#kw##nlsolve)(::Array{Any,1}, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,2},Array{Float64,2},Array{Float64,2},Val{false}}, ::Array{Float64,2}) at .\<missing>:0
[10] #nlsolve#39(::Symbol, ::Float64, ::Float64, ::Int32, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int32, ::Float64, ::Symbol, ::Bool, ::NLsolve.#nlsolve, ::#f!#222{Float64,#foc1#220}, ::Array{Float64,2}) at C:\Users\arquie\.julia\v0.6\NLsolve\src\nlsolve\nlsolve.jl:59
[11] (::NLsolve.#kw##nlsolve)(::Array{Any,1}, ::NLsolve.#nlsolve, ::Function, ::Array{Float64,2}) at .\<missing>:0
[12] bk_opti_decent(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Function, ::Function, ::Float64) at .\In[37]:38

``````

The function that you pass to `quadgk` should not have the argument. Moreover, I believe initial_x should be a column vector, not a row vector, i.e. you should put commas between the numbers: `initial_x=[0.5, 0.5, 1.2]`

For a minimal working example, you should post something really minimal, like the following:

``````using NLsolve
using ForwardDiff

function bk_opti_decent()

function foc1(Sb, c_bar)
c1(θ) = Sb * θ + c_bar
end

E = 4

function f!(F, x)
F[1] = foc1(x[1], x[2])
F[2] = x[2] - E
end

initial_x = [0.5, 0.2]
essai = nlsolve(f!, initial_x, autodiff = :forward)
end
``````
1 Like

Thanks a lot! I corrected the initial points error. Thanks also for your patience (first time coding and therefore posting about it but Ill make my future posts are a bit more efficient…)

I tried solving a minimalist version of my problem, using 1. Oncedifferentiable (i was getting method error using forwardiff) and 2. `mcpsolve` because i was still getting some domain error (still from this negative number / exponential problem I thought) and I realized anyway that I forgot to make sure that consumptions had to be non null anyway in my setting.

`````` initial_x = [0.5; 0.5 ]

initial_F = similar(initial_x)

df = OnceDifferentiable(f!, initial_x, initial_F)

results1_bk = mcpsolve(df, [0., -Inf ],  [Inf, Inf], initial_x )
``````