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)
μ=matread("parameters.mat")["mu"]
σ=matread("parameters.mat")["sigma"]
μ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