# Domain issue when working with nlsolve

I am quite stumped on this problem so appreciate in advance any suggestions. I am currently working with a trade model and trying to solve a system of 4 equations using “nlsolve.” I am running into issues with the domain and have tried a few workarounds, but new errors pop up. Below is a MWE of the larger problem. Thanks again.

Function used in code:

``````function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j)  # k_e: (1000, 10, 11)  fe: (11, 11)  Fe:  (11, 11)
μ = param.μ
σ = param.σ
μj = μ[j]
σj = σ[j]
ρ = param.ρ
δ = param.δ
β =  sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
disc = ρ + δ - μj
F[1] =       Z_AB[3]*Z_AB[1]^α +        k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*Z_AB[1]^β
F[2] =       Z_AB[3]*Z_AB[2]^α +        k_e*Z_AB[2]/disc - fe/disc      - Z_AB[4]*Z_AB[2]^β
F[3] =     α*Z_AB[3]*Z_AB[1]^(α-1) +    k_e/disc    - β*Z_AB[4]*Z_AB[1]^(β-1)
F[4] =     α*Z_AB[3]*Z_AB[2]^(α-1) +    k_e/disc    - β*Z_AB[4]*Z_AB[2]^(β-1)
end
``````

And the main code that calls the function:

``````N = 10
USA = N + 1
ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
η = 5.0 # elasticity of substitution
H = η^(-η)*(η-1)^(η-1)
P = ones(1,N+1) # price indexes
n = 1000
T = 31 # number of periods
s = 100
b = 1 # lower bound of the (Pareto) productivity distribution
θ = 3 # shape parameter of the (Pareto) productivity distribution
z=10*ones(1000,1)
Y_simu=ones(1,11)
Q=10*ones(1,11)
fh=10*ones(1,10)
Fh=10*ones(1,10)
τ=ones(11,11)
w=ones(1,11)
fe=ones(11,11)
Fe=ones(11,11)
μ_USj=10*ones(11,1)
σ_USj=10*ones(11,1)
struct Parameter
μ::Matrix{Float64}
σ::Matrix{Float64}
ρ::Float64
δ::Float64
η::Float64
end
β = sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
α = -sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
disc = ρ + δ .- μ_USj
disc = disc'

# Preallocation (n firms, N host countries + the US)
k_e     = zeros(n,N,N+1)
Bo_0    = zeros(n,N,N+1)
Z_HE_0  = 10*zeros(n,N,N+1)
Z_EH_0  = 10*zeros(n,N,N+1)
init_HE           = zeros(n,N,N+1,4)
thresholds_HE_sol = zeros(n,N,N+1,4)
init_OH           = zeros(n,N,4)
thresholds_OH_sol = zeros(n,N,T,4)
status_exp = zeros(n,N,N+1,T)
status     = zeros(n,N,T)
k_h = H*(1 ./ repeat(z, 1, N)) .^(1-η) .* repeat((P[:,1:10] .^η) .* Q[:,1:10],n)
Z_OH_0   = repeat(β[1:N]' ./ (β[1:N].-1)', n, 1) .*  (k_h.^-1) .* (repeat(fh,n,1) + (ρ.+δ.-repeat(μ_USj[1:10]',n,1)) .* repeat(Fh,n,1))
Z_HO_0   = (repeat(β[1:N]' ./ (β[1:N].-1)', n, 1).* (k_h.^-1) .* repeat(fh,n,1))
Bo_0_H= 100*ones(1000,10) # k_h ./ repeat(disc[:,1:N],n,1) ./ repeat(β[1:N]',n,1) ./ Z_OH_0.^(repeat(β[1:N]',n,1).-1)
for c in 1:N
for d in 1:N+1
if d != c
k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d]
Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d]))
Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d])
Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1))
end
end
end

# SOLVING FOR THE THRESHOLDS AND THE PARAMETERS OF THE VALUE FUNCTIONS
for j in 1:2
for c in 1:N
for d in 1:N+1
if d != c
k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] # 1000 x 1 vector
Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d])) # 1000 x 1 vector
Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d]) # 1000 x 1 vector
Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1)) # 1000 x 1 vector
end
end
end
# initialization
k = N+1
thresholds_HE_sol = zeros(n,N,N+1,4)
for i in eachindex(z)
if k != j
# creating initial Hessian
init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
if i > 20 && j != 3 && isnan(sum(thresholds_HE_sol[i-1,j,k,:])) == 0
init_HE[i,j,k,:] = thresholds_HE_sol[i-1,j,k,:]
end
if init_HE[i,j,k,3] < 0
init_HE[i,j,k,3] = 0.05
end

# initial matrices and function to be minimized
initial_x = init_HE[i,j,k,:]
f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e[i,j,k], fe[j,k], Fe[j,k], param, j)
result = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true).zero
end
end
end
``````

The analogous code in Matlab (below) runs without error. Thank you in advance for helping me understand the difference.

``````f_thresholds_HE  = @(Z_thresholds_exp)sys_tHE_AB(Z_thresholds_exp, k_e(i,j,k), fe(j,k), Fe(j,k),param,k);
[thresholds_HE_sol(i,j,k,:),fval] = fsolve(f_thresholds_HE,init_HE(i,j,k,:),options);
``````

It might take a bit for someone to attempt to help, as the example you are sharing is very heavy. In the meantime, something that can help the helpers:

• Any chance you can simplify this further? Could you give just one single `nlsolve` call that fails in an obvious way. Currently you have a very large amount of code dedicated to just setup and the nlsolve is called multiple times inside of a deep loop.

• Could you elaborate on how it fails? Is there an error being raised (in which case, could you share the exception stack trace)? Is there an answer, but the answer is obviously wrong (in which case, could you explain why it is obviously wrong)?

2 Likes

It’s hard to see where exactly it is because the error message is missing from your post, but probably the simplest way to solve it could be to use NaNMath.jl on whichever of the `^` calls is failing so that it throws a NaN instead of a DomainError, which would help the system recover and re-adapt its steps.

1 Like

Thank you both @Krastanov and @ChrisRackauckas for your feedback. @ChrisRackauckas, thanks so much for the NaNMath.jl suggestion!! I implemented it, but now am getting a new error, which I think reflects the fact that nlsolve is not ignoring the “NaN”? Any ideas on how to fix this new problem? I am elaborating below on the error code and including the stack trace. Thanks a lot again.

Code:

``````# Description: this file calibrates trade costs from each country to US
using MAT, Statistics, DataFrames, LinearAlgebra
using NLsolve
using NaNMath; nm=NaNMath

function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j)  # k_e: (1000, 10, 11)  fe: (11, 11)  Fe:  (11, 11)
μ = param.μ
σ = param.σ
μj = μ[j]
σj = σ[j]
ρ = param.ρ
δ = param.δ
β =  sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
disc = ρ + δ - μj

F[1] =       Z_AB[3]*nm.Z_AB[1]^α +        k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*nm.Z_AB[1]^β
F[2] =       Z_AB[3]*nm.Z_AB[2]^α +        k_e*Z_AB[2]/disc - fe/disc      - Z_AB[4]*nm.Z_AB[2]^β
F[3] =     α*Z_AB[3]*nm.Z_AB[1]^(α-1) +    k_e/disc    - β*Z_AB[4]*nm.Z_AB[1]^(β-1)
F[4] =     α*Z_AB[3]*nm.Z_AB[2]^(α-1) +    k_e/disc    - β*Z_AB[4]*nm.Z_AB[2]^(β-1)
end

N = 10
USA = N + 1
ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
η = 5.0 # elasticity of substitution
H = η^(-η)*(η-1)^(η-1)
P = ones(1,N+1) # price indexes
n = 1000
T = 31 # number of periods
s = 100
b = 1 # lower bound of the (Pareto) productivity distribution
θ = 3 # shape parameter of the (Pareto) productivity distribution
z=10*ones(1000,1)
Y_simu=ones(1,11)
Q=10*ones(1,11)
fh=10*ones(1,10)
Fh=10*ones(1,10)
τ=ones(11,11)
w=ones(1,11)
fe=ones(11,11)
Fe=ones(11,11)
μ_USj=10*ones(11,1)
σ_USj=10*ones(11,1)

# creating a structure of the parameters
struct Parameter
μ::Matrix{Float64}
σ::Matrix{Float64}
ρ::Float64
δ::Float64
η::Float64
end
β = sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
α = -sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
disc = ρ + δ .- μ_USj
disc = disc'

# Preallocation (n firms, N host countries + the US)
k_e     = zeros(n,N,N+1)
Bo_0    = zeros(n,N,N+1)
Z_HE_0  = 10*zeros(n,N,N+1)
Z_EH_0  = 10*zeros(n,N,N+1)
init_HE           = zeros(n,N,N+1,4)
thresholds_HE_sol = zeros(n,N,N+1,4)
init_OH           = zeros(n,N,4)
thresholds_OH_sol = zeros(n,N,T,4)
status_exp = zeros(n,N,N+1,T)
status     = zeros(n,N,T)
k_h = H*(1 ./ repeat(z, 1, N)) .^(1-η) .* repeat((P[:,1:10] .^η) .* Q[:,1:10],n)
Z_OH_0   = repeat(β[1:N]' ./ (β[1:N].-1)', n, 1) .*  (k_h.^-1) .* (repeat(fh,n,1) + (ρ.+δ.-repeat(μ_USj[1:10]',n,1)) .* repeat(Fh,n,1))
Z_HO_0   = (repeat(β[1:N]' ./ (β[1:N].-1)', n, 1).* (k_h.^-1) .* repeat(fh,n,1))
Bo_0_H= 100*ones(1000,10) # k_h ./ repeat(disc[:,1:N],n,1) ./ repeat(β[1:N]',n,1) ./ Z_OH_0.^(repeat(β[1:N]',n,1).-1)
for c in 1:N
for d in 1:N+1
if d != c
k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d]
Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d]))
Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d])
Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1))
end
end
end

# SOLVING FOR THE THRESHOLDS AND THE PARAMETERS OF THE VALUE FUNCTIONS
for j in 1:2
for c in 1:N
for d in 1:N+1
if d != c
k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] # 1000 x 1 vector
Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d])) # 1000 x 1 vector
Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d]) # 1000 x 1 vector
Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1)) # 1000 x 1 vector
end
end
end
# initialization
k = N+1
thresholds_HE_sol = zeros(n,N,N+1,4)
for i in eachindex(z)
if k != j
# creating initial Hessian
init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
if i > 20 && j != 3 && isnan(sum(thresholds_HE_sol[i-1,j,k,:])) == 0
init_HE[i,j,k,:] = thresholds_HE_sol[i-1,j,k,:]
end
if init_HE[i,j,k,3] < 0
init_HE[i,j,k,3] = 0.05
end

# initial matrices and function to be minimized
initial_x = init_HE[i,j,k,:]
Z_AB=ones(4)
f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e[i,j,k], fe[j,k], Fe[j,k], param, j)
result = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true).zero
end
end
end
``````

Stack trace:

``````ERROR: UndefVarError: Z_AB not defined
Stacktrace:
[1] getproperty(x::Module, f::Symbol)
@ Base ./Base.jl:31
[2] sys_tHE_AB!(F::Vector{Float64}, Z_AB::Vector{Float64}, k_e::Float64, fe::Float64, Fe::Float64, param::Parameter, j::Int64)
@ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:17
[3] (::var"#f_thresholds_HE!#15"{Int64, Int64, Int64})(F::Vector{Float64}, Z_AB::Vector{Float64})
@ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:116
[4] (::NLSolversBase.var"#fj_finitediff!#21"{var"#f_thresholds_HE!#15"{Int64, Int64, Int64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:central}(), Float64}})(F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:138
[5] value_jacobian!!(obj::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:124
[6] value_jacobian!!
@ ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:122 [inlined]
[7] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
[8] trust_region (repeats 2 times)
@ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
[9] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
[10] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Pairs{Symbol, Real, NTuple{5, Symbol}, NamedTuple{(:xtol, :ftol, :iterations, :store_trace, :extended_trace), Tuple{Float64, Float64, Int64, Bool, Bool}}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
[11] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:117
``````

Check which line of your code causes that error to appear (and generally share that when reporting errors). On that line you are probably trying to access a variable called `Z_AB`. The error message is saying that you have not defined such a variable in that scope. Maybe due to a typo or a mixup between local and global variables?

1 Like

That doesn’t make sense. It would be `nm.^`, `Z_AB` should be unchanged.

1 Like

Thanks very much @ChrisRackauckas!! The (correct) use of `nm.` solved the domain error.

Now when I run the code, it gives a method matching error. I’ve been reading online and experimenting with changes, but it is not clear to me how to resolve it. Any suggestions on what to do?
Thank you so much again.

Updated function:

``````function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j)
μ = param.μ
σ = param.σ
μj = μ[j]
σj = σ[j]
ρ = param.ρ
δ = param.δ
β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
disc = ρ + δ - μj

F[1] =   Z_AB[3]*Z_AB[1]nm.^α +        k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*Z_AB[1]nm.^β
F[2] =   Z_AB[3]*Z_AB[2]nm.^α +        k_e*Z_AB[2]/disc - fe/disc      - Z_AB[4]*Z_AB[2]nm.^β
F[3] = α*Z_AB[3]*Z_AB[1]nm.^(α-1) +    k_e/disc    - β*Z_AB[4]*Z_AB[1]nm.^(β-1)
F[4] = α*Z_AB[3]*Z_AB[2]nm.^(α-1) +    k_e/disc    - β*Z_AB[4]*Z_AB[2]nm.^(β-1)

F = [F[1] F[2] F[3] F[4]]
end
``````

Stack trace:

``````ERROR: MethodError: no method matching length(::Module)
Closest candidates are:
length(::Union{Base.KeySet, Base.ValueIterator}) at abstractdict.jl:58
length(::Union{HDF5.Attribute, HDF5.Dataset}) at ~/.julia/packages/HDF5/HtnQZ/src/dataspaces.jl:289
length(::Union{HDF5.File, HDF5.Group}) at ~/.julia/packages/HDF5/HtnQZ/src/groups.jl:70
...
Stacktrace:
[1] _similar_shape(itr::Module, #unused#::Base.HasLength)
@ Base ./array.jl:663
[2] _collect(cont::UnitRange{Int64}, itr::Module, #unused#::Base.HasEltype, isz::Base.HasLength)
@ Base ./array.jl:718
[3] collect(itr::Module)
@ Base ./array.jl:712
[6] sys_tHE_AB!(F::Vector{Float64}, Z_AB::Vector{Float64}, k_e::Float64, fe::Float64, Fe::Float64, param::Parameter, j::Int64)
@ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:15
[7] (::var"#f_thresholds_HE!#11"{Int64, Int64, Int64})(F::Vector{Float64}, Z_AB::Vector{Float64})
@ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:114
[8] (::NLSolversBase.var"#fj_finitediff!#21"{var"#f_thresholds_HE!#11"{Int64, Int64, Int64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:central}(), Float64}})(F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:138
[9] value_jacobian!!(obj::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:124
[10] value_jacobian!!
@ ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:122 [inlined]
[11] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
[12] trust_region (repeats 2 times)
@ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
[13] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
[14] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Pairs{Symbol, Real, NTuple{5, Symbol}, NamedTuple{(:xtol, :ftol, :iterations, :store_trace, :extended_trace), Tuple{Float64, Float64, Int64, Bool, Bool}}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
[15] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:115
``````

Just do `const ^ = NanMath.^` before the function definition.

1 Like

Sorry for being dense about this @ChrisRackauckas, but when I put `const ^ = NanMath .^` 1) before the function definition, or 2) in the function, before `F[1]=...`, or 3) in the main file beginning, I get the following errors. Do I need to define it as a function?

``````julia> const ^ = NanMath .^
ERROR: syntax: incomplete: premature end of input
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5
``````
``````ERROR: syntax: invalid function name ".^" around /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:15
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:2
``````
``````julia> # Description: this file calibrates trade costs from each country to US
using MAT, Statistics, DataFrames, LinearAlgebra
using NLsolve
using NaNMath
nm = NaNMath
const ^ = NanMath .^
include("/Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl")
ERROR: cannot assign a value to variable Base.^ from module Main
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:6
``````

You have a space between `NaNMath` and `.^`. As Chris said, it should be `NanMath.^`.

1 Like

Apologies, that was a typo when I was writing about it here. However even without a space yields same errors (depending on placement):

``````julia> const ^ = NanMath.^
ERROR: syntax: incomplete: premature end of input
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5
``````

or (placed in file before function definition)

``````julia> using NaNMath
include("/Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl")
ERROR: LoadError: UndefVarError: ^ not defined
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:1
in expression starting at /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:1
``````

Any thoughts? Thanks!

You have

`using NaNMath`
`const ^ = NanMath.^`

where the second has a lower case n. The package name should be the same in both lines.

1 Like

Thanks @DanDeepPhase! The error persists with the correct spelling (and no space):
`const ^ = NaNMath.^`.

``````ERROR: UndefVarError: ^ not defined
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5
``````

Or:

``````ERROR: syntax: invalid function name ".^" around /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:3
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:2
``````

NaNMath doesn’t have a `^` function, but it does have a `pow` function.

Maybe try: `const ^ = NaNMath.pow`

debugging the debugging…Maybe that gets you to the fail state Chris was proposing.

1 Like

Thanks a lot for the followup, @DanDeepPhase!! I’m not having much success with `const ^ = NaNMath.pow` either (see error below). However now I’ve edited the code to use `max(var,0)` as a workaround (based on this post: DifferentialEquations.jl, Error: Exponentiation yielding a complex result requires a complex argument - #6 by Gregers).

``````ERROR: cannot assign a value to variable Base.^ from module Main
Stacktrace:
[1] top-level scope
@ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5
``````

Updated function:

``````function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, ρ, δ, μ, σ, μj, σj, j) # k_e: (1000, 10, 11)  fe: (11, 11)  Fe:  (11, 11)

ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
μj = μ[j]
σj = σ[j]
β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
disc = ρ + δ - μj

F[1] =   Z_AB[3]*max(Z_AB[1],0)^α +        k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*max(Z_AB[1],0)^β
F[2] =   Z_AB[3]*max(Z_AB[2],0)^α +        k_e*Z_AB[2]/disc - fe/disc      - Z_AB[4]*max(Z_AB[2],0)^β
F[3] = α*Z_AB[3]*max(Z_AB[1],0)^(α-1) +    k_e/disc    - β*Z_AB[4]*max(Z_AB[1],0)^(β-1)
F[4] = α*Z_AB[3]*max(Z_AB[2],0)^(α-1) +    k_e/disc    - β*Z_AB[4]*max(Z_AB[2],0)^(β-1)

N = [F[1] F[2] F[3] F[4]]

end
``````

I see you found another solution but for future reference:

The error `cannot assign a value to variable Base.^ from module Main` occurs when you try to define `^` after you have already used the standard `^` in Main. To make it work, restart Julia and make sure you define ` const ^ = NaNMath.pow` before using `^`.

1 Like

That is incredibly helpful @sijo – thank you!! This works, and likely more “correct” that using the max operator. Thanks again.

1 Like