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);