BoundsError in Bisection Algorithm using a while loop

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

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

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)

results = NLopt.optimize(opt, [1e-5, 1e-5])
opt_inp = results[2]
land[i] = opt_inp[1]

# 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

1. write a general bisection function (a few lines, or just use one from Roots.jl),
2. implement your problem as `struct` and a residual/discrepancy calculator function,
3. just pass a closure to the bisection algorithm.
1 Like

it is sorted. but thank you!

Hi
You are missing the point. You need to have functions not a plain script.

no, it works with the plain script now! I just got right of the splicing of the “beq” array and it is fine

1 Like