In [1]:
using QuantEcon
using Optim
using BasisMatrices
using Plots
plotlyjs()
Out[1]:
In [2]:
type ProductionAdjustmentModel
alpha::Float64
beta ::Float64
kappa::Float64
sigma::Float64
delta::Float64
max_iter::Int
eps_nodes::Vector{Float64}
eps_weights::Vector{Float64}
basis::BasisMatrices.Basis
basis_values::Matrix{Float64}
snodes::Matrix{Float64}
d_vec::Vector{Float64}
l_vec::Vector{Float64}
function ProductionAdjustmentModel(
alpha::Float64 = 0.5,
beta ::Float64 = 0.5,
kappa::Float64 = 0.5,
sigma::Float64 = 0.4,
delta::Float64 = 0.9,
l_num::Int = 10,
l_min::Float64 = 0.0,
l_max::Float64 = 10.0,
eps_num::Int = 3,
d_num::Int = 15,
max_iter::Int = 100
)
eps_nodes, eps_weights = qnwlogn(eps_num, -(sigma^2/2), sigma^2)
basis = Basis(
ChebParams(d_num, minimum(eps_nodes), maximum(eps_nodes)),
ChebParams(l_num, l_min, l_max)
)
snodes, (d_vec, l_vec) = nodes(basis)
basis_values = BasisMatrix(basis, Expanded(), snodes, 0).vals[1]
return new(
alpha, beta, kappa, sigma, delta, max_iter,
eps_nodes, eps_weights, basis, basis_values, snodes, d_vec, l_vec
)
end
end
In [3]:
function cost(PAM::ProductionAdjustmentModel, quantity::Number)
return PAM.kappa * quantity
end
function adjustmentCost(PAM::ProductionAdjustmentModel, quantityDiff::Number)
return 0.5 * PAM.alpha * (quantityDiff^2)
end
function price(PAM::ProductionAdjustmentModel, demandShock::Number, quantity::Number)
return demandShock * quantity^(-PAM.beta)
end
function reward(
PAM::ProductionAdjustmentModel,
demandShock::Number,
previousProduction::Number,
currentProduction::Number
)
return price(PAM, demandShock, currentProduction) *
currentProduction - cost(PAM, currentProduction) -
adjustmentCost(PAM, currentProduction - previousProduction)
end
Out[3]:
In [11]:
function update_bellman(PAM::ProductionAdjustmentModel, coefs::Vector{Float64})
opt_values = similar(coefs)
opt_actions = similar(coefs)
for i in 1:size(PAM.snodes, 1)
d, l = PAM.snodes[i, :]
function objective(q::Float64)
expected_v = 0.0
for (error_v, error_weight) in zip(PAM.eps_nodes, PAM.eps_weights)
expected_v += error_weight * funeval(coefs, PAM.basis, [error_v, q])
end
return -1 * (reward(PAM, d, l, q) + PAM.delta * expected_v)
end
opt = optimize(objective, 0, 100)
opt_values[i] = -opt.minimum
opt_actions[i] = opt.minimizer
end
new_coefs = PAM.basis_values \ opt_values
return new_coefs, opt_actions
end
Out[11]:
In [17]:
PAM = ProductionAdjustmentModel(
0.5, # alpha
0.5, # beta
0.5, # kappa
0.4, # sigma
0.9, # delta
10, # l_num
0.0, # l_min
4.0, # l_max
3, # eps_num
15, # d_num
200 # max_iter
)
Out[17]:
In [18]:
eps = 1e-6
total_states = size(PAM.basis_values, 1)
initial_values = ones(total_states)
initial_coefs = PAM.basis_values \ initial_values
count = 0
coefs = initial_coefs
# just for scope problem
new_coefs = similar(coefs); opt_actions = similar(coefs);
while true
new_coefs, opt_actions = update_bellman(PAM, coefs)
max_error = maximum(abs, coefs - new_coefs)
if max_error < eps
println("Loop finished in $(count+1) / $(PAM.max_iter) loops")
println("Final maximum error: $(max_error)")
break
end
coefs = new_coefs
count += 1
if count >= PAM.max_iter
println("Loop not finished in $(PAM.max_iter) loops")
println("Final maximum error: $(max_error)")
break
end
end
opt_coefs = new_coefs;
In [19]:
for i in 1:20 #length(opt_actions)
d, l = PAM.snodes[i, :]
a = opt_actions[i]
@printf "d: %0.4f, l: %0.4f, opt_action: %0.4f\n" d l a
end
In [21]:
opt_action_matrix = reshape(opt_actions, length(PAM.d_vec), length(PAM.l_vec))
Plots.surface(
PAM.d_vec, PAM.l_vec, opt_action_matrix',
xlims=(0.5, 2),
ylims=(0, 4),
zlims=(0, 4),
xlabel = "demand shock",
ylabel = "lagged production",
size=(800, 450),
guidefont = font(8, "Courier")
)
Out[21]:
In [23]:
value_matrix = Array{Float64, 2}(length(PAM.d_vec), length(PAM.l_vec))
for (i, d) in enumerate(PAM.d_vec)
for (j, l) in enumerate(PAM.l_vec)
value_matrix[i, j] = funeval(opt_coefs, PAM.basis, [d, l])
end
end
Plots.surface(
PAM.d_vec, PAM.l_vec, value_matrix',
xlims=(0.5, 2),
ylims=(0, 4),
zlims=(3, 8),
xlabel = "demand shock",
ylabel = "lagged production",
guidefont = font(8, "Courier")
)
Out[23]:
In [ ]: