I ran the following code for 7 hours without any results, how can I modify it? Thanks
using Optim, LinearAlgebra, Statistics, Plots, QuantEcon, Distributions, Interpolations, NLsolve, Random, JuMP, GLPK
import Ipopt
β = 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
grid_z = range(1, 10, length = 50)
grid_b = range(0, 5, length = 1000)
grid_b = collect(grid_b)
grid_z = collect(grid_z)
dist_z = Pareto(θ)
dist_b = Pareto(θ)
wage = 1
rate = 1
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_choice1 = zeros(1000,1)
l_choice1 = zeros(1000,1)
c_choice2 = zeros(1000,1)
k_choice2 = zeros(1000,1)
l_choice2 = zeros(1000,1)
c_choice3 = zeros(1000,1)
k_choice3 = zeros(1000,1)
l_choice3 = zeros(1000,1)
#################################### bellman begin #################################
for i in 1:50
################################### choice1 ##########################
function T1(w)
# 插值函数t1
function t1(x)
w_func = LinearInterpolation(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-δ))/(1-δ) + β * (w_func((b*(1+rate)+wage-x))))
for j in 1:Nb
set_value(b, grid_b[j])
optimize!(choice1)
Tw1[j,1] = objective_value(choice1)
c_choice1[j,1] = value(x)
end
return Tw1
end
################################### choice2 ##########################
function T2(w)
# 插值函数t2
function t2(x)
w_func2 = LinearInterpolation(grid_b, w, extrapolation_bc=Line())
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-δ))/(1-δ) + β * (w_func2(-x + (1+rate)*(b-y) + (1-δ)*y + grid_z[i]*((y^α)*(z^(1-α))^(1-v))-wage*z)))
@NLconstraint(choice2, y <= b)
for j in 1:Nb
set_value(b, 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
################################### choice3 ##########################
function T3(w)
function t3(x)
w_func3 = LinearInterpolation(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-δ))/(1-δ) + β * (w_func3(-x - wage*z + (1-δ)*y - (1+rate+χ)*(y-b+ψ) +grid_z[i]*((y^α)*(z^(1-α))^(1-v)))))
@NLconstraint(choice3, ξ*y <= b - ψ)
@NLconstraint(choice3, y >= b)
for j in 1:Nb
set_value(b, 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
################################### solve ##########################
ww1 = grid_b
ww2 = grid_b
ww3 = grid_b
for mm in 1:500
Tw1 = T1(ww1)
ww1 = [Tw1[j,1] for j in 1:1000]
Tw2 = T2(ww1)
ww2 = [Tw2[j,1] for j in 1:1000]
Tw3 = T3(ww1)
ww4 = [Tw3[j,1] for j in 1:1000]
end
func1 = LinearInterpolation(grid_b, ww1, extrapolation_bc=Line())
func2 = LinearInterpolation(grid_b, ww2, extrapolation_bc=Line())
func3 = LinearInterpolation(grid_b, ww3, extrapolation_bc=Line())
for j in 1:1000
b = 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
#################################### Bellman end #################################