Model solving - method error using NLsolve

axelle · August 8, 2018, 12:22pm

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)))
           
         return quadgk(val1, 0, θlow)[1] + quadgk(val2, θlow, θhigh)[1] + quadgk(val3, θhigh, 1)[1] + val4
        
        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:
 [1] evalrule(::#val1#292{#c21#289{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64},Float64,##324#328}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Base.LinAlg.#vecnorm) at C:\Users\arquie\.julia\v0.6\QuadGK\src\QuadGK.jl:66
 [2] do_quadgk(::#val1#292{#c21#289{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64},Float64,##324#328}, ::Array{Float64,1}, ::Int32, ::Type{Float64}, ::Float64, ::Float64, ::Int32, ::Base.LinAlg.#vecnorm) at C:\Users\arquie\.julia\v0.6\QuadGK\src\QuadGK.jl:131
 [3] #quadgk#15(::Array{Any,1}, ::Function, ::Function, ::Int32, ::Float64) at C:\Users\arquie\.julia\v0.6\QuadGK\src\QuadGK.jl:241
 [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.

mzaffalon · August 8, 2018, 1:51pm

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] + ...

axelle · August 8, 2018, 3:06pm

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.

mzaffalon · August 8, 2018, 3:20pm

Is there any chance you can post u_prime?

axelle · August 8, 2018, 3:23pm

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)

pkofod · August 8, 2018, 4:24pm

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

mzaffalon · August 8, 2018, 5:05pm

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

axelle · August 8, 2018, 6:21pm

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 ?

Many thanks for your help.

pkofod · August 8, 2018, 6:24pm

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

axelle · August 9, 2018, 8:04am

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…

Many thanks for your answers!

pkofod · August 9, 2018, 8:34am

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.

axelle · August 9, 2018, 8:53am

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!!

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)

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!

mzaffalon · August 9, 2018, 11:34am

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.

axelle · August 9, 2018, 1:36pm

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 * θ

        return  quadgk(c1(θ), 0, θlow)[1]
        
        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

Thanks for all your help.

mzaffalon · August 9, 2018, 2:35pm

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
using QuadGK

function bk_opti_decent()

    function foc1(Sb, c_bar)
        c1(θ) = Sb * θ + c_bar
        quadgk(c1, 0, 2)[1]
    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

axelle · August 10, 2018, 9:12am

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 )