@stevengj @odow
Sorry for not providing this in the first place; I was hoping the issue was basic enough to not have to adapt the code to something reproducible without loading the data and a bunch of more complicated function. I’m still getting the issue with this simpler version:
using Parameters, NLopt, Distributions
@with_kw mutable struct Primitives
αₖ::Float64 = 0.75 # capital share
αₗ::Float64 = 0.2 # land share
b::Float64 = 1.2 # capital cost shape parameter
σₙ::Float64 = 5.0 # cost scale standard deviation
end
@with_kw mutable struct Model
prim::Primitives # model primitives
L::Vector{Float64} # land sqft
K::Vector{Float64} # gross sqft
u::Vector{Int64} # unit count
K̅::Float64 = 1000.0 # normalizing level of capital
end
# save mapping between primary primitives and their respective indices
PRIMITIVE_MAPPING = Dict(
"αₖ" => 1,
"αₗ" => 2,
"b" => 3,
"σₙ" => 4
)
# function that parses the parameters into primitives
function parse_primitives(θ::Vector{Float64})
prim = Primitives()
prim.αₖ = θ[PRIMITIVE_MAPPING["αₖ"]]
prim.αₗ = θ[PRIMITIVE_MAPPING["αₗ"]]
prim.b = θ[PRIMITIVE_MAPPING["b"]]
prim.σₙ = θ[PRIMITIVE_MAPPING["σₙ"]]
return prim
end
# η likelihood
function ℓₙ(mod::Model; tol::Float64=1e-8)
@unpack K, L, u, K̅ = mod
@unpack αₖ, αₗ, b, σₙ = mod.prim
# compute price from the capital FOC
p = (K.^(b-αₖ)).*(b/(αₖ*(K̅^b))).*(u.^(αₖ+αₗ-1)).*(L.^(-αₗ))
# compute the upper and lower bounds of the integral
ζ = (b*(1-αₖ-αₗ))/(b-αₖ)
A = ((b-αₖ)/b)*((αₖ*((K̅^b)/b))^(αₖ/(b-αₖ)))*((p.*(L.^αₗ)).^(b/(b-αₖ)))
c₁ = A.*(((u .+ 1).^(ζ)) .- (u.^ζ))
c₂ = A.*((u.^ζ) .- ((u.-1).^ζ))
ϕ = Normal(0, σₙ)
out = cdf(ϕ, c₂) .- cdf(ϕ, c₁)
return log.(out)
end
# ν likelihood
function ℓᵥ(mod::Model)
@unpack K, L, u, K̅ = mod
@unpack αₖ, αₗ, b, σᵥ = mod.prim
μ = (b-αₖ)*log.(K) .- log((αₖ*(K̅^b))/b) .- (1-αₖ-αₗ)*log.(u) .- αₗ*log.(L)
out = log(b - αₖ) .- 0.5 * log(2 * π * σᵥ^2) .- ((μ .^ 2) / (2 * σᵥ^2)) .- log.(K)
return out
end
function ℓ(θ::Vector{Float64}, mod::Model)
mod.prim = parse_primitives(θ, mod; cov=cov)
return sum(ℓₙ(mod)) + sum(ℓᵥ(mod))
end
L = [1487, 84223, 11760, 688995, 15380, 12040, 3872, 18460, 16807, 28088102, 22711, 65568, 61450, 9440, 5030, 12845, 1388, 20010, 14044051, 14044051, 7220, 12295, 8410, 10030,
6180, 10695, 4198, 9600, 7690, 10170, 10010, 5510, 63752, 14645, 10550, 8490, 8385, 2695, 3910, 3910, 93230, 9060, 31460, 8400, 5173, 9950, 11000, 18480, 5970, 10000, 7050, 9600, 2203, 5020, 19950, 3296, 13610, 3790, 5020, 12300, 6131, 2750, 4480, 2400, 4950, 4000, 4150, 8640, 20125, 14955, 6720, 10163, 20000, 64040, 10400, 7125, 7203, 8100,
14640, 5013, 9500, 13680, 54835, 5940, 5000, 8240, 5770, 6720, 4810, 6810, 8960, 3600, 3690, 4120, 5100, 3330, 6575, 9150, 5500, 5130, 15380, 5190, 35186, 4400, 7655, 14044051, 5000, 4160, 4779, 5134, 4800, 4600, 1949, 5973, 625, 10800, 9657, 5000, 8000, 4400, 3603, 920, 5500, 3608, 8000, 9000, 8250, 10798, 3264, 6110, 4798, 8000, 7228, 4277, 4070, 12000, 15388, 3360, 7501, 8467, 9323, 32957, 8365, 7215, 4000, 5998, 8014, 5969, 5746, 5000, 5003, 4004, 4121, 5115, 2649, 1031, 778, 4960, 2842, 8400, 3600, 1271, 1026, 3200, 4780, 4780, 24203, 30000, 30000, 30000, 30000, 84223, 3228, 5015, 4210, 7210, 6300, 7200, 2824, 5986, 5950, 7210, 1806, 7200, 3990, 4800, 6020, 3599, 5770, 4090, 4400, 3599, 14645, 14645, 14645, 5970, 5000, 2708, 5690, 3990, 4815, 4650, 3450, 3790, 1137, 9380, 5210, 6889, 4090, 4700, 4000, 4320, 9530, 9120, 6500, 3680, 3290, 1730,
4080, 5000, 5100, 4500, 5144, 4501, 4000, 3510, 4000, 3000, 5000, 3300, 1152, 3300, 4120, 2550, 3416, 6087, 3750, 1061, 2700, 2661, 2183, 4927, 4170, 4900, 4201, 5020, 3839, 1758, 8145, 53264, 5500, 1340, 4317, 4316, 4453, 2183, 1812, 1812, 697, 1400, 5400, 4178, 5400]
K = [3659.0, 3028.0, 35653.0, 31571.0, 84389.0, 84577.0, 13640.0, 9619.0, 13770.0, 61576.0, 4042.0, 14663.0, 354555.0, 18244.0, 7079.0, 14395.0, 2934.0, 121218.0, 15938.0, 28291.0, 22367.0, 34896.0, 16371.0, 10485.0, 9138.0, 31761.0, 12989.0, 30106.0, 18291.0, 25690.0, 23567.0, 16120.0, 280830.0, 6008.0, 23453.0, 22694.0, 14666.0, 6268.0, 2515.0, 2119.0, 119878.0, 2581.0, 175844.0, 11250.0, 10010.0, 36685.0, 30266.0, 62766.0, 11717.0, 25688.0, 14727.0, 6006.0, 2406.0, 12630.0, 67190.0, 6033.0, 36956.0, 8208.0, 8198.0, 58322.0, 8823.0, 5656.0, 3423.0, 4321.0, 10400.0, 3601.0, 9620.0, 22156.0, 49460.0, 54089.0, 11634.0, 24452.0, 8768.0, 63084.0, 9374.0, 11911.0, 5184.0, 13746.0, 36656.0, 8444.0, 21583.0, 4135.0, 133072.0, 20267.0, 13832.0, 20260.0, 14296.0, 12166.0, 10082.0, 11802.0, 24528.0, 2400.0, 8079.0, 9582.0, 5360.0, 7510.0, 11730.0, 31600.0, 10010.0, 7209.0, 30035.0, 3564.0, 7294.0, 6194.0, 17484.0, 209143.0, 11500.0, 5400.0, 10104.0, 39416.0, 11870.0, 8770.0, 1581.0, 15672.0, 2540.0, 6714.0, 11650.0, 4460.0, 20178.0, 6889.0, 4615.0, 9688.0, 13182.0, 2160.0, 13413.0, 6342.0, 14515.0, 31390.0, 6045.0, 12178.0, 3785.0, 20588.0, 17384.0, 8097.0, 9562.0, 32509.0, 38819.0, 4772.0, 17385.0, 23980.0, 8390.0, 60050.0, 3949.0, 20850.0, 9209.0, 8475.0, 17904.0, 14958.0, 10500.0, 20868.0, 10638.0, 6712.0, 10750.0, 18270.0, 3931.0, 2092.0, 2330.0, 3132.0, 2624.0,
22615.0, 4068.0, 2974.0, 2937.0, 3440.0, 2224.0, 3676.0, 6325.0, 7991.0, 4817.0, 4613.0, 3057.0, 2905.0, 2102.0, 4201.0, 4420.0, 7195.0, 3691.0, 18674.0, 4656.0, 15054.0, 14029.0, 7195.0, 2125.0, 11000.0, 7869.0, 4421.0, 14228.0, 2300.0, 9457.0, 9228.0, 11894.0, 2300.0, 2916.0, 3177.0, 2944.0, 11717.0, 10105.0, 3572.0, 17843.0, 7789.0, 11000.0, 12702.0, 8902.0, 8992.0, 3403.0, 4970.0, 8331.0, 16491.0, 7538.0, 7091.0, 15069.0, 8068.0, 17243.0, 4081.0, 4534.0, 5589.0, 4706.0, 2217.0, 7016.0, 7806.0, 4370.0, 7328.0, 11229.0, 1723.0, 5215.0, 5188.0, 10453.0, 6495.0, 2151.0, 5272.0, 2975.0, 1678.0, 4062.0, 5484.0, 5233.0, 3425.0, 8529.0, 4148.0, 4857.0, 4028.0, 2170.0, 11059.0, 11765.0, 12421.0, 4713.0, 2420.0, 9470.0, 2474.0, 9811.0, 12222.0, 9089.0, 4617.0, 9956.0, 9881.0, 5533.0, 2170.0, 2170.0, 2170.0, 4117.0, 2878.0, 4903.0, 10280.0, 5064.0]
u = [3, 2, 21, 25, 75, 63, 8, 10, 3, 52, 2, 10, 234, 12, 5, 16, 2, 111, 15, 39, 22, 23, 50, 6, 30, 75, 11, 34, 23, 35, 60, 23, 193, 3, 29, 60, 8, 8, 4, 4, 112, 3, 150, 47, 16,
34, 7, 64, 6, 47, 20, 12, 2, 30, 15, 2, 49, 4, 24, 45, 30, 16, 8, 8, 8, 10, 24, 28, 3, 61, 35, 57, 20, 62, 16, 34, 18, 41, 51, 27, 34, 8, 111, 36, 27, 37, 41, 35, 3, 35, 49, 2, 23, 28, 9, 20, 29, 75, 8, 20, 67, 2, 20, 4, 41, 148, 22, 10, 27, 32, 22, 20, 2, 20, 2, 4, 6, 30, 64, 27, 12, 4, 33, 2, 42, 16, 6, 31, 3, 23, 8, 38, 42, 31, 20, 42, 66, 8, 41, 52, 4, 59, 2, 29, 31, 19, 40, 36, 22, 46, 19, 15, 18, 45, 3, 2, 2, 6, 2, 8, 3, 2, 2, 2, 2, 3, 2, 6, 4, 4, 4, 2, 2, 2, 5, 4, 2, 8, 16, 8, 7, 4, 2, 8, 8, 8, 3, 3, 4,
5, 5, 3, 3, 3, 3, 6, 6, 8, 8, 8, 5, 6, 6, 6, 2, 7, 6, 7, 12, 6, 8, 8, 60, 2, 2, 2, 8, 2, 7, 13, 2, 8, 8, 2, 8, 6, 13, 5, 2, 2, 2, 2, 6, 4, 3, 2, 8, 2, 14, 3, 2, 4, 26, 8, 9, 8, 26, 3, 8, 8, 5, 2, 30, 30, 5, 2, 2, 2, 2, 2, 2, 20, 2]
mod = Model(prim=Primitives(), L=L, u=u, K=K)
θ₀ = [0.7, 0.2, 1.2, 1.0]
θ_lb = [1e-8, 1e-8, 1 + 1e-8, 1e-8]
θ_ub = [1 - 1e-8, 1 - 1e-8, Inf, Inf]
opt = Opt(:LN_COBYLA, length(θ₀))
function c(result::Vector, x::Vector, grad::Matrix)
result[1] = -(x[PRIMITIVE_MAPPING["αₖ"]] + x[PRIMITIVE_MAPPING["αₗ"]])
result[2] = x[PRIMITIVE_MAPPING["αₖ"]] + x[PRIMITIVE_MAPPING["αₗ"]] - 1.0
result[3] = 1e-12 + x[PRIMITIVE_MAPPING["αₖ"]] - x[PRIMITIVE_MAPPING["b"]]
end
opt.lower_bounds = θ_lb
opt.upper_bounds = θ_ub
inequality_constraint!(opt, c, ones(length(θ₀)) * 1e-12)
opt.xtol_rel = 1e-3
opt.max_objective = (x, g!) -> ℓ(x)
(ml, θ₊, ret) = optimize(opt, θ₀)