Hi all,
I am programming a bisection search. My algorithm is supposed to do this:
take a parameter value (a price) as input for an optimization, use outputs as input for a function whose roots I want to find (excess demand function so that the outputs of that optimization give me aggregate supply and demand), then update price using bisection method and optimize again. I get a first update, but then the loop breaks because it can go into the optimization loop within the while loop because it tries to call an element out of the array length. Do you have any idea where it goes wrong. The code is:
const pn=1.0; const r=0.05; const q=3.0; const Aa=2.9 ; An = 1
const γ=0.54; const α₁= 0.1*10.0e-1; const α₂=0.1*10.0e-2; const β₁= 3.9;
const β₂= 0.9; const μ = 3.0 ; const λ = 10000.0 ; β = 0.96 ; η = 0.7 ; const a = 0.01 ; const τ = 0.5 ; const θ = 1.01
m = 2 #number of inputs in utility function
P = ((1/m)*(pa^(1-θ) + pn^(1-θ)))^(1/(1-θ)) # price index Blanchard, Kiyotaki (1987)
iter = 100
beq = Array{Float64}(iter + 1) # bequests
occupation = Array{Float64}(iter) # indicator: 1=worker, 2=farmer
beq[1] = 0 #set initial bequests
yₐ = zeros(Array{Float64}(iter)) # production farmers
land = zeros(Array{Float64}(iter)) #land demand
adap = zeros(Array{Float64}(iter)) #demand adaptive capital ∈ non-agricultural demand
damage = zeros(Array{Float64}(iter)) #calculate climate damage
Inc = zeros(Array{Float64}(iter)) #income of worker or farmer
cₐmid = zeros(Array{Float64}(iter))
cₐold = zeros(Array{Float64}(iter))
ρ = 0.75
σ = 0.02
n = 1000
s = linspace(0.1, 0.5, n) ### state values
default_mc = tauchen(n, ρ, σ)
mc = MarkovChain(default_mc.p, s)
s_vec = simulate(mc, iter, init=Int(n/2))
top = 5.0
bottom = 10.0e-3
pa = (top+bottom)/2
m = 2 #number of inputs in utility function
P = ((1/m)*(pa^(1-θ) + pn^(1-θ)))^(1/(1-θ)) # price index Blanchard, Kiyotaki (1987)
agg_demand_cₐ = 10; #arbitray to get alg. stared
agg_supply_yₐ = 9;
excess_dem_a = agg_demand_cₐ - agg_supply_yₐ
### now comes the whole loop
while abs(excess_dem_a) >1e-5
let excess_dem_a = (agg_demand_cₐ - agg_supply_yₐ)
let pa = (top+bottom)/2
println(pa)
println(excess_dem_a)
#now comes the whole optimzation
for i = 1:iter
# println(i)
function obj(x::Vector, grad::Vector)
if length(grad) > 0
grad[1] = pa*((1 + β₁*x[2]^β₂)/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s_vec[i])^(1-γ)*γ*x[1]^(γ-1) - (1+r)*q + q
grad[2] = pa*((β₂*β₁*x[2]^(β₂-1)*(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2) - (β₂*β₁^(2)*x[2]^(2*β₂ -1) + β₂*β₁*x[2]^(β₂-1)))/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2)^2)*(Aa*s_vec[i])^(1-γ)*x[1]^γ - (1+r)*pn
end
pa*((1 + β₁*x[2]^β₂)/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s_vec[i])^(1-γ)*x[1]^γ - (1+r)*(q*x[1] + pn*x[2]) + q*x[1]
end
function myconstraint(x::Vector, grad::Vector)
if length(grad) > 0
grad[1] = λ*q + (1+r)*q
grad[2] = (1+r)*pn
end
-1* λ*q*x[1] - (1+r)*(beq[i] - q*x[1] - pn*x[2])
end
opt = Opt(:LD_MMA, 2) #algorithm, dimensionality of problem
lower_bounds!(opt, [1e-5, 1e-5])
xtol_rel!(opt,1e-4)
max_objective!(opt, obj)
inequality_constraint!(opt, (x,grad) -> myconstraint(x,grad), 1e-8)
results = NLopt.optimize(opt, [1e-5, 1e-5])
opt_inp = results[2]
land[i] = opt_inp[1]
adap[i] = opt_inp[2]
# yₐ = ((1 + β₁* opt_inp[2]^β₂)/(1 + β₁*opt_inp[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s_vec[i])^(1-γ)*opt_inp[1]^γ
π = pa*((1 + β₁* opt_inp[2]^β₂)/(1 + β₁*opt_inp[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s_vec[i])^(1-γ)*opt_inp[1]^γ - (1+r)*(q*opt_inp[1] + pn*opt_inp[2]) + q*opt_inp[1]
damage[i] = (1 + β₁* opt_inp[2]^β₂)/(1 + β₁*opt_inp[2]^β₂ + α₁*μ + α₂*μ^2)
if π < An
beq[i+1] = ((1-η)/η)*(Aa + (1+r)*beq[i])/(1 + ((1-η)/η) + β/η)
occupation[i] = 1.0 #for worker
Inc[i] = An
else
beq[i+1] = ((1-η)/η)*(π + (1+r)*beq[i])/(1 + ((1-η)/η) + β/η)
occupation[i] = 2.0 #for farmer
yₐ[i] = ((1 + β₁* opt_inp[2]^β₂)/(1 + β₁*opt_inp[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s_vec[i])^(1-γ)*opt_inp[1]^γ
Inc[i] = π
end
#aggregates vectors
cₐmid[i] = ((Inc[i] + (1+r)*beq[i])/(P*(1 + (1-η)/(η) + β/η)) -pa*a)*(((pn/pa)*(τ/(1-τ)))^θ/(pn + pn^θ*pa*(1-θ)*(τ/(1-τ))^θ)) + a
cₐold[i] = (((β*(1+r))/η)*(Inc[i] + (1+r)*beq[i])/(P*(1 + (1-η)/(η) + β/η)) -pa*a)*(((pn/pa)*(τ/(1-τ)))^θ/(pn + pn^θ*pa*(1-θ)*(τ/(1-τ))^θ)) + a
end
splice!(beq,1);
distrib = [occupation Inc s_vec beq yₐ land adap damage]'
#Now I can compute aggregate demand and uptdate the price
#LHS: aggregte demand
agg_demand_cₐ = sum(cₐmid) + sum(cₐold)
#RHS: aggregate supply
agg_supply_yₐ = sum(yₐ)
if agg_demand_cₐ - agg_supply_yₐ >0.0
top = pa
else
bottom = pa
end
end
end
end
Thanks a lot