Hello, I got the following error when I ran the following code：setindex! not defined for StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
I found out he had run T1. But why is this wrong? thanks
using Optim, LinearAlgebra, Statistics, Plots, QuantEcon, Distributions, Interpolations, NLsolve, Random, JuMP, GLPK
import Ipopt
function CareerWorkerProblem(;β = 0.88,
θ = 3.4,
Nb = 1000,
Nz = 50,
v = 0.19,
α = 0.33,
δ = 1.5,
χ = 0.05,
ϵ = 0.335,
γ = 0.89,
ψ = 0.95,
ξ = 0.8 )
wage = 1
rate = 1
grid_z = range(1, 10, length = 50)
grid_b = range(0, 5, length = 1000)
dist_z = Pareto(θ)
z_probs = pdf.(dist_z, grid_z)
return (β=β, θ=θ, Nb=Nb, Nz=Nz,v=v,α=α, δ=δ, χ=χ, ϵ=ϵ, γ=γ, ψ=ψ, wage=wage, rate=rate, grid_z=grid_z, grid_b=grid_b, dist_z=dist_z, z_probs=z_probs)
end
out = zeros(50,1000)
action = zeros(50,1000)
Tw1 = zeros(1000,1)
Tw2 = zeros(1000,1)
Tw3 = zeros(1000,1)
c_choice1 = zeros(1000,1)
k_choice2 = zeros(1000,1)
l_choice3 = zeros(1000,1)
c_choice1 = zeros(1000,1)
k_choice2 = zeros(1000,1)
l_choice3 = zeros(1000,1)
c_choice1 = zeros(1000,1)
k_choice2 = zeros(1000,1)
l_choice3 = zeros(1000,1)
#################################### bellman begin #################################
function bellman(cp)
for i in 1:50
################################### choice1 ##########################
function T1(w)
# 插值函数t1
function t1(x)
w_func = LinearInterpolation(cp.grid_b, w,extrapolation_bc=Line())
return w_func(x)
end
# 优化模型choice1
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
register(choice1,:w_func,1,t1; autodiff = true)
@NLparameter(choice1, b == 0)
@NLobjective(choice1, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (w_func((b*(1+cp.rate)+cp.wage-x))))
for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice1)
Tw1[j,1] = objective_value(choice1)
c_choice1[j,1] = value(x)
end
return Tw1
end
function solve_optgrowth1(initial_w; tol = 1e-6, max_iter = 500)
fixedpoint(w -> T1(w), initial_w).zero
end
################################### choice2 ##########################
function T2(w)
# 插值函数t2
function t2(x)
w_func2 = LinearInterpolation(cp.grid_b, w)
return w_func2(x)
end
#优化模型choice2
choice2 = Model(Ipopt.Optimizer)
@variable(choice2, x >= 0)
@variable(choice2, y >= 0)
@variable(choice2, z >= 0)
register(choice2, :w_func2, 1, t2; autodiff = true)
@NLparameter(choice2, b == 1)
@NLobjective(choice2, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (w_func2(-x + (1+cp.rate)*(b-y) + (1-cp.δ)*y + cp.grid_z[i]*((y^α)*(z^(1-α))^(1-cp.v))-cp.wage*z)))
@NLconstraint(choice2, x[2] <= b)
for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice2)
Tw2[j,1] = objective_value(choice2)
c_choice2[j,1] = value(x)
k_choice2[j,1] = value(y)
l_choice2[j,1] = value(z)
end
return Tw2
end
function solve_optgrowth2(initial_w; tol = 1e-6, max_iter = 500)
fixedpoint(w -> T2(w), initial_w).zero # gets returned
end
################################### choice3 ##########################
function T3(w)
function t3(x)
w_func3 = LinearInterpolation(cp.grid_b, w,extrapolation_bc=Line())
return w_func3(x)
end
choice3 = Model(Ipopt.Optimizer)
@variable(choice3, x >= 0)
@variable(choice3, y >= 0)
@variable(choice3, z >= 0)
register(choice3, :w_func3, 1, t3; autodiff = true)
@NLparameter(choice3, b == 1)
@NLobjective(choice3, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (w_func3(-x - cp.wage*z + (1-cp.δ)*y - (1+cp.rate+cp.χ)*(y-b+cp.ψ) +cp.grid_z[i]*((y^α)*(z^(1-α))^(1-cp.v)))))
@NLconstraint(choice3, cp.ξ*y <= b - cp.ψ)
@NLconstraint(choice3, x[2] >= b)
for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice3)
Tw3[j,1] = objective_value(choice3) # 返回生产率为zi，初始财富为bj下的最优值函数
c_choice3[j,1] = value(x)
k_choice3[j,1] = value(y)
l_choice3[j,1] = value(z)
end
return Tw3
end
function solve_optgrowth3(initial_w; tol = 1e-6, max_iter = 500)
fixedpoint(w -> T3(w), initial_w).zero # gets returned
end
################################### solve ##########################
initial_w= cp.grid_b
v_star_approx1 = solve_optgrowth1(initial_w)
v_star_approx2 = solve_optgrowth2(initial_w)
v_star_approx3 = solve_optgrowth3(initial_w.(cp.grid_b))
func1 = LinearInterpolation(cp.grid_b, v_star_approx1,extrapolation_bc=Line())
func2 = LinearInterpolation(cp.grid_b, v_star_approx2,extrapolation_bc=Line())
func3 = LinearInterpolation(cp.grid_b, v_star_approx3,extrapolation_bc=Line())
for j in 1:cp.Nb
b = cp.grid_b[j]
v1 = func1(b)
v2 = func2(b)
v3 = func3(b)
if v1 > max(v2, v3)
action = 1
elseif v2 > max(v1, v3)
action = 2
else
action = 3
end
action[i,j] = action
out[i,j] = max(v1, v2, v3)
end
end
return action, out
end
#################################### Bellman end #################################
wp = CareerWorkerProblem()
action, out = bellman(wp)