BoundsError in Bisection Algorithm using a while loop


#1

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


#2

Organize your code:

  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.

#3

it is sorted. but thank you!


#4

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


#5

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


#6

https://docs.julialang.org/en/v0.6.4/manual/performance-tips/