using QuantEcon
using BasisMatrices
using ContinuousDPs
using Plots
import PyPlot
using LinearAlgebra
n = 10
s_min, s_max = 5, 10
basis = Basis(ChebParams(n, s_min, s_max))
1 dimensional Basis on the hypercube formed by (5.0,) × (10.0,). Basis families are Cheb
alpha = 0.2
bet = 0.5
gamm = 0.9
sigma = 0.1
discount = 0.9;
x_star = ((discount * bet) / (1 - discount * gamm))^(1 / (1 - bet))
s_star = gamm * x_star + x_star^bet
s_star, x_star
(7.416897506925212, 5.6094182825484795)
f(s, x) = (s - x)^(1 - alpha) / (1 - alpha)
g(s, x, e) = gamm * x .+ e * x^bet;
n_shocks = 3
shocks, weights = qnwlogn(n_shocks, -sigma^2/2, sigma^2) # See Errata
([0.8367708003019956, 0.9950124791926823, 1.1831792330610165], [0.16666666666666674, 0.6666666666666666, 0.16666666666666674])
x_lb(s) = 0
x_ub(s) = 0.99 * s;
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}(f, g, 0.9, [0.8367708003019956, 0.9950124791926823, 1.1831792330610165], [0.16666666666666674, 0.6666666666666666, 0.16666666666666674], x_lb, x_ub, ContinuousDPs.Interp{1,Array{Float64,1},Array{Float64,2},LU{Float64,Array{Float64,2}}}(1 dimensional Basis on the hypercube formed by (5.0,) × (10.0,). Basis families are Cheb , [5.030779148512155, 5.27248368952908, 5.732233047033631, 6.365023750651133, 7.108913837399423, 7.891086162600577, 8.634976249348867, 9.267766952966369, 9.72751631047092, 9.969220851487844], ([5.030779148512155, 5.27248368952908, 5.732233047033631, 6.365023750651133, 7.108913837399423, 7.891086162600577, 8.634976249348867, 9.267766952966369, 9.72751631047092, 9.969220851487844],), 10, (10,), (5.0,), (10.0,), [1.0 -0.9876883405951379 … 0.30901699437495467 -0.1564344650402394; 1.0 -0.8910065241883679 … -0.8090169943749478 0.4539904997395474; … ; 1.0 0.8910065241883679 … -0.8090169943749478 -0.4539904997395474; 1.0 0.9876883405951375 … 0.3090169943749398 0.15643446504022196], LU{Float64,Array{Float64,2}}([1.0 -0.9876883405951379 … 0.30901699437495467 -0.1564344650402394; 1.0 1.9753766811902755 … -1.4876988529977098e-14 0.31286893008046135; … ; 1.0 0.04894348370484648 … -5.257311121191343 1.6448493055872506; 1.0 0.5791922201622681 … -0.459649548425358 5.062325628940014], [1, 10, 5, 7, 5, 8, 7, 9, 10, 10], 0)))
res = solve(cdp, VFI);
Compute iterate 50 with error 0.00856936646624007 Compute iterate 100 with error 4.416458845568627e-5 Compute iterate 150 with error 2.2761436113682976e-7 Compute iterate 176 with error 1.470634813927063e-8 Converged in 176 steps
s_min, s_max = cdp.interp.lb[1], cdp.interp.ub[1]
grid_length = 500
s_nodes = range(s_min, stop=s_max, length=grid_length)
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Investment Policy"
xlabel = "Wealth"
ylabel = "Investment (% of Wealth)"
plot(s_nodes, X./s_nodes, xlims=(s_min, s_max), ylims=(0.65, 0.9),
title=title, xlabel=xlabel, ylabel=ylabel, label="Chebychev")
plot!([s_star], [x_star/s_star], m=(7,:star8), label="")
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), yformatter=:scientific,
title=title, xlabel=xlabel, ylabel=ylabel, label="")
res = solve(cdp, PFI);
Compute iterate 6 with error 2.1316282072803006e-14 Converged in 6 steps
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Investment Policy"
xlabel = "Wealth"
ylabel = "Investment (% of Wealth)"
plot(s_nodes, X./s_nodes, xlims=(s_min, s_max), ylims=(0.65, 0.9),
title=title, xlabel=xlabel, ylabel=ylabel, label="Chebychev")
plot!([s_star], [x_star/s_star], m=(7,:star8), label="")
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), yformatter=:scientific,
title=title, xlabel=xlabel, ylabel=ylabel, label="")
n_yrs = 20
n_paths = 2000
s_paths = Array{Float64}(undef, n_yrs+1, n_paths)
s_init = 5.
for i in 1:n_paths
s_paths[:, i] = simulate(res, s_init, n_yrs+1)
end
mean_s_path = mean(s_paths, dims=2);
title = "Expected Wealth"
xlabel = "Year"
ylabel = "Wealth"
plot(0:n_yrs, mean_s_path, ylims=(5, 7.5),
title=title, xlabel=xlabel, ylabel=ylabel, label="")
k = 3
m = 101
breaks = m - (k-1)
s_min, s_max = 0, 10
basis = Basis(SplineParams(breaks, s_min, s_max, k))
1 dimensional Basis on the hypercube formed by (0.0,) × (10.0,). Basis families are Spline
a = [10, 0.8]
b = [12, 1.0]
discount = 0.9;
p(x) = a[1] - a[2] * x / 2
c(s, x) = b[1] * x - b[2] * x * (2*s - x) / 2
f(s, x) = p(x) * x - c(s, x)
g(s, x, e) = s - x;
shocks, weights = [0.], [1.];
x_lb(s) = 0
x_ub(s) = s;
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}(f, g, 0.9, [0.0], [1.0], x_lb, x_ub, ContinuousDPs.Interp{1,Array{Float64,1},SparseArrays.SparseMatrixCSC{Float64,Int64},SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}}(1 dimensional Basis on the hypercube formed by (0.0,) × (10.0,). Basis families are Spline , [0.0, 0.034013605442176874, 0.10204081632653061, 0.20408163265306123, 0.30612244897959184, 0.40816326530612246, 0.5102040816326532, 0.6122448979591838, 0.7142857142857144, 0.8163265306122449 … 9.18367346938775, 9.285714285714297, 9.387755102040822, 9.48979591836735, 9.591836734693876, 9.693877551020401, 9.795918367346948, 9.897959183673473, 9.96598639455783, 10.0], ([0.0, 0.034013605442176874, 0.10204081632653061, 0.20408163265306123, 0.30612244897959184, 0.40816326530612246, 0.5102040816326532, 0.6122448979591838, 0.7142857142857144, 0.8163265306122449 … 9.18367346938775, 9.285714285714297, 9.387755102040822, 9.48979591836735, 9.591836734693876, 9.693877551020401, 9.795918367346948, 9.897959183673473, 9.96598639455783, 10.0],), 101, (101,), (0.0,), (10.0,), [1 , 1] = 1.0 [2 , 1] = 0.296296 [1 , 2] = 0.0 [2 , 2] = 0.564815 [3 , 2] = 0.25 [1 , 3] = 0.0 [2 , 3] = 0.132716 [3 , 3] = 0.583333 [4 , 3] = 0.166667 [1 , 4] = 0.0 [2 , 4] = 0.00617284 [3 , 4] = 0.166667 ⋮ [100, 98] = 0.00617284 [101, 98] = 0.0 [98 , 99] = 0.166667 [99 , 99] = 0.583333 [100, 99] = 0.132716 [101, 99] = 0.0 [98 , 100] = 1.64861e-40 [99 , 100] = 0.25 [100, 100] = 0.564815 [101, 100] = 0.0 [99 , 101] = 4.22045e-41 [100, 101] = 0.296296 [101, 101] = 1.0, SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}(Ptr{Nothing} @0x00007fb4d4697d90, Ptr{Nothing} @0x00007fb4ef81cf10, 101, 101, [0, 2, 5, 9, 14, 18, 22, 26, 30, 35 … 367, 371, 375, 380, 384, 388, 393, 397, 401, 404], [0, 1, 0, 1, 2, 0, 1, 2, 3, 0 … 98, 99, 100, 97, 98, 99, 100, 98, 99, 100], [0.9999999999999999, 0.2962962962962961, 0.0, 0.5648148148148147, 0.24999999999999997, 0.0, 0.13271604938271606, 0.5833333333333333, 0.16666666666666666, 0.0 … 0.5833333333333246, 0.13271604938266363, 0.0, 1.6486136302935093e-40, 0.2500000000000261, 0.5648148148147707, 0.0, 4.220450893551383e-41, 0.29629629629639687, 0.9999999999999999], 0)))
res = solve(cdp, VFI);
Compute iterate 38 with error 1.4533888759160618e-8 Converged in 38 steps
s_min, s_max = cdp.interp.lb[1], cdp.interp.ub[1]
grid_length = 500
s_nodes = range(s_min, stop=s_max, length=grid_length)
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Harvest Policy"
xlabel = "Available Stock"
ylabel = "Harvest"
plot(s_nodes, X, xlims=(s_min, s_max), ylims=(0, 10),
title=title, xlabel=xlabel, ylabel=ylabel, label="Cubic Spline")
B1 = evalbase(cdp.interp.basis.params[1], s_nodes, 1)
shadow_prices = B1 * res.C;
title = "Shadow Price Function"
ylabel = "Price"
plot(s_nodes, shadow_prices, xlims=(s_min, s_max), ylims=(-0.5, 6.5),
title=title, xlabel=xlabel, ylabel=ylabel, label="Cubic Spline")
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), ylims=(-1e-5, 1.5e-5), yformatter=:scientific,
title=title, xlabel=xlabel, ylabel=ylabel, label="")
res = solve(cdp, PFI);
Compute iterate 6 with error 1.046872007484656e-13 Converged in 6 steps
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Harvest Policy"
xlabel = "Available Stock"
ylabel = "Harvest"
plot(s_nodes, X, xlims=(s_min, s_max), ylims=(0, 10),
title=title, xlabel=xlabel, ylabel=ylabel, label="Cubic Spline")
B1 = evalbase(cdp.interp.basis.params[1], s_nodes, 1)
shadow_prices = B1 * res.C;
title = "Shadow Price Function"
ylabel = "Price"
plot(s_nodes, shadow_prices, xlims=(s_min, s_max), ylims=(-0.5, 6.5),
title=title, xlabel=xlabel, ylabel=ylabel, label="Cubic Spline")
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), ylims=(-1e-5, 1.5e-5), yformatter=:scientific,
title=title, xlabel=xlabel, ylabel=ylabel, label="")
n_yrs = 20
s_init = 10.
s_path = simulate(res, s_init, n_yrs+1);
title = "State Path"
xlabel = "Year"
ylabel = "Stock"
plot(0:n_yrs, s_path, ylims=(2, 10),
title=title, xlabel=xlabel, ylabel=ylabel, label="")
bet = [0.8 0.5; 0.2 0.6]
gamm = [-0.8, 0.0]
Omega = [0.3 0.0; 0.0 1.0]
s_target = [0., 1.]
alpha = [0.9, 0.4]
Sigma = 0.04 * Matrix(I, 2, 2)
discount = 0.9;
f(s::Vector{Float64}, x::Float64) =
-(s - s_target)' * Omega * (s - s_target) / 2
g(s::Vector{Float64}, x::Float64, e::Vector{Float64}) =
alpha + bet * s + gamm * x + e
x_lb(s) = 0.
x_ub(s) = 100;
n_shocks = [3, 3]
mu = [0, 0]
shocks, weights = qnwnorm(n_shocks, mu, Sigma)
([-0.34641016151377546 -0.34641016151377546; 0.0 -0.34641016151377546; … ; 0.0 0.34641016151377546; 0.34641016151377546 0.34641016151377546], [0.027777777777777804, 0.11111111111111116, 0.027777777777777804, 0.11111111111111116, 0.4444444444444444, 0.11111111111111116, 0.027777777777777804, 0.11111111111111116, 0.027777777777777804])
k = [3, 3]
m = [10, 10]
breaks = m - (k.-1)
s_min = [-15, -10]
s_max = [15, 10]
basis = Basis(map(SplineParams, breaks, s_min, s_max, k)...)
2 dimensional Basis on the hypercube formed by (-15.0, -10.0) × (15.0, 10.0). Basis families are Spline × Spline
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
ContinuousDP{2,Array{Float64,2},Array{Float64,2},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}(f, g, 0.9, [-0.34641016151377546 -0.34641016151377546; 0.0 -0.34641016151377546; … ; 0.0 0.34641016151377546; 0.34641016151377546 0.34641016151377546], [0.027777777777777804, 0.11111111111111116, 0.027777777777777804, 0.11111111111111116, 0.4444444444444444, 0.11111111111111116, 0.027777777777777804, 0.11111111111111116, 0.027777777777777804], x_lb, x_ub, ContinuousDPs.Interp{2,Array{Float64,2},SparseArrays.SparseMatrixCSC{Float64,Int64},SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}}(2 dimensional Basis on the hypercube formed by (-15.0, -10.0) × (15.0, 10.0). Basis families are Spline × Spline , [-15.0 -10.0; -13.571428571428575 -10.0; … ; 13.571428571428575 10.0; 15.0 10.0], ([-15.0, -13.571428571428575, -10.714285714285714, -6.428571428571431, -2.142857142857139, 2.142857142857139, 6.428571428571431, 10.714285714285714, 13.571428571428575, 15.0], [-10.0, -9.04761904761905, -7.142857142857143, -4.285714285714287, -1.4285714285714282, 1.4285714285714282, 4.285714285714287, 7.142857142857142, 9.047619047619047, 10.0]), 100, (10, 10), (-15.0, -10.0), (15.0, 10.0), [1 , 1] = 1.0 [2 , 1] = 0.296296 [11 , 1] = 0.296296 [12 , 1] = 0.0877915 [1 , 2] = 0.0 [2 , 2] = 0.564815 [3 , 2] = 0.25 [4 , 2] = 1.78017e-47 [11 , 2] = 0.0 [12 , 2] = 0.167353 [13 , 2] = 0.0740741 [14 , 2] = 5.27457e-48 ⋮ [88 , 99] = 0.0740741 [89 , 99] = 0.167353 [90 , 99] = 0.0 [97 , 99] = 1.78017e-47 [98 , 99] = 0.25 [99 , 99] = 0.564815 [100, 99] = 0.0 [88 , 100] = 0.0 [89 , 100] = 0.0877915 [90 , 100] = 0.296296 [98 , 100] = 0.0 [99 , 100] = 0.296296 [100, 100] = 1.0, SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}(Ptr{Nothing} @0x00007fb4d5bddb40, Ptr{Nothing} @0x00007fb4d5bef270, 100, 100, [0, 4, 12, 20, 32, 40, 46, 58, 66, 74 … 1524, 1532, 1540, 1552, 1560, 1566, 1578, 1586, 1594, 1600], [0, 1, 10, 11, 0, 1, 2, 3, 10, 11 … 96, 97, 98, 99, 87, 88, 89, 97, 98, 99], [1.0, 0.29629629629629745, 0.29629629629629695, 0.08779149519890314, 0.0, 0.5648148148148143, 0.2499999999999999, 1.7801680491237497e-47, 0.0, 0.1673525377229083 … 1.7801680491237497e-47, 0.2499999999999999, 0.5648148148148143, 0.0, 0.0, 0.08779149519890289, 0.2962962962962961, 0.0, 0.29629629629629745, 1.0], 0)))
res = solve(cdp, VFI);
Compute iterate 50 with error 0.000823893836127354 Compute iterate 100 with error 3.319121503864153e-6 Compute iterate 150 with error 1.7104639482568018e-8 Compute iterate 152 with error 1.3854446478944737e-8 Converged in 152 steps
s_min, s_max = cdp.interp.lb, cdp.interp.ub
grid_length = collect(cdp.interp.size) * 5
grids = [range(s_min[i], stop=s_max[i], length=grid_length[i]) for i in 1:2]
set_eval_nodes!(res, grids...)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Monetary Policy"
xlabel = "GDP Gap"
ylabel = "Inflation Rate"
zlabel = "Nominal Interest Rate"
PyPlot.surf(grids..., permutedims(reshape(X, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Approximation Residual"
zlabel = "Residual"
PyPlot.surf(grids..., permutedims(reshape(resid, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
point = ([0., 1.0], 0., [0., 0.])
res_lqa = solve(cdp, LQA, point=point);
res = solve(cdp, PFI, v_init=res_lqa.V);
Compute iterate 4 with error 1.1937117960769683e-12 Converged in 4 steps
set_eval_nodes!(res, grids...)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Monetary Policy"
xlabel = "GDP Gap"
ylabel = "Inflation Rate"
zlabel = "Nominal Interest Rate"
PyPlot.surf(grids..., permutedims(reshape(X, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Approximation Residual"
zlabel = "Residual"
PyPlot.surf(grids..., permutedims(reshape(resid, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
N = 2
n_yrs = 20
n_paths = 5000
s_paths = Array{Float64}(undef, N, n_yrs+1, n_paths)
s_init = [15., 10.]
for i in 1:n_paths
s_paths[:, :, i] = simulate(res, s_init, n_yrs+1)
end
mean_s_paths = [mean(s_paths[i, :, :], dims=2) for i in 1:N];
xlabel = "Year"
ylabel = "GDP Gap"
title = "Expected State Path: $ylabel"
plot(0:n_yrs, mean_s_paths[1], ylims=(-5, 15),
title=title, xlabel=xlabel, ylabel=ylabel, label="")
xlabel = "Year"
ylabel = "Inflation Rate"
title = "Expected State Path: $ylabel"
plot(0:n_yrs, mean_s_paths[2], ylims=(0, 10),
title=title, xlabel=xlabel, ylabel=ylabel, label="")
alpha = 0.5
bet = 0.5
kappa = 0.5
sigma = 0.4
discount = 0.9;
prod_cost(q::Float64) = kappa * q
adj_cost(dq::Float64) = 0.5 * alpha * dq^2
f(s::Vector{Float64}, x::Float64) =
s[1] * x^(1-bet) - prod_cost(x) - adj_cost(x - s[2])
g(s::Vector{Float64}, x::Float64, e::Float64) = [e, x]
x_lb(s) = 0.
x_ub(s) = 100;
n_shocks = 3
shocks, weights = qnwlogn(n_shocks, -sigma^2/2, sigma^2) # See Errata
([0.46170906161702985, 0.9231163463866358, 1.8456293363222578], [0.16666666666666674, 0.6666666666666666, 0.16666666666666674])
n = [15, 10]
x_star = ((1 - bet) / kappa)^(1/bet)
s_min = [shocks[1], x_star - 1]
s_max = [shocks[end], x_star + 3]
basis = Basis(map(ChebParams, n, s_min, s_max)...)
2 dimensional Basis on the hypercube formed by (0.46170906161702985, 0.0) × (1.8456293363222578, 4.0). Basis families are Cheb × Cheb
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
ContinuousDP{2,Array{Float64,1},Array{Float64,2},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}(f, g, 0.9, [0.46170906161702985, 0.9231163463866358, 1.8456293363222578], [0.16666666666666674, 0.6666666666666666, 0.16666666666666674], x_lb, x_ub, ContinuousDPs.Interp{2,Array{Float64,2},Array{Float64,2},LU{Float64,Array{Float64,2}}}(2 dimensional Basis on the hypercube formed by (0.46170906161702985, 0.0) × (1.8456293363222578, 4.0). Basis families are Cheb × Cheb , [0.46549969165043137 0.02462331880972446; 0.49557600132395074 0.02462331880972446; … ; 1.8117623966153367 3.9753766811902755; 1.841838706288856 3.9753766811902755], ([0.46549969165043137, 0.49557600132395074, 0.5544141416161106, 0.6394426034594216, 0.7469452350595032, 0.8722236555605718, 1.0098025968341728, 1.1536691989696437, 1.2975358011051146, 1.4351147423787156, 1.5603931628797842, 1.6678957944798656, 1.7529242563231768, 1.8117623966153367, 1.841838706288856], [0.02462331880972446, 0.2179869516232642, 0.5857864376269049, 1.0920190005209063, 1.687131069919538, 2.3128689300804615, 2.9079809994790935, 3.414213562373095, 3.7820130483767356, 3.9753766811902755]), 150, (15, 10), (0.46170906161702985, 0.0), (1.8456293363222578, 4.0), [1.0 -0.9945218953682733 … 0.03252455412868575 -0.016351854232752223; 1.0 -0.9510565162951535 … -0.09194987150091097 0.048340908203385165; … ; 1.0 0.9510565162951535 … -0.09194987150091097 -0.048340908203385165; 1.0 0.9945218953682732 … 0.032524554128683575 0.016351854232749843], LU{Float64,Array{Float64,2}}([1.0 -0.9945218953682733 … 0.03252455412868575 -0.016351854232752223; 1.0 1.9890437907365466 … -0.06504910825736931 2.3800406090401793e-15; … ; 1.0 0.9353978569087952 … -23.594422471056248 -35.068145938985964; 1.0 0.12638086253384775 … 0.747238274932304 56.74145205172091], [1, 15, 8, 5, 11, 8, 13, 9, 9, 12 … 141, 145, 147, 146, 150, 147, 147, 150, 149, 150], 0)))
res = solve(cdp, VFI);
Compute iterate 50 with error 0.0029795551247095986 Compute iterate 100 with error 1.535595732971018e-5 Compute iterate 150 with error 7.914115407459121e-8 Compute iterate 166 with error 1.4665014980153046e-8 Converged in 166 steps
s_min, s_max = cdp.interp.lb, cdp.interp.ub
grid_length = collect(cdp.interp.size) * 5
grids = [range(s_min[i], stop=s_max[i], length=grid_length[i]) for i in 1:2]
set_eval_nodes!(res, grids...)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Production Policy"
xlabel = "Demand Shock"
ylabel = "Lagged Production"
zlabel = "Production"
zlims = (0, 4)
PyPlot.surf(grids..., permutedims(reshape(X, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_zlim](zlims)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Value Function"
zlabel = "Value"
PyPlot.surf(grids..., permutedims(reshape(V, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Approximation Residual"
zlabel = "Residual"
PyPlot.surf(grids..., permutedims(reshape(resid, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
ax[:ticklabel_format](style="sci", axis="z", scilimits=(0,0))
res = solve(cdp, PFI);
Compute iterate 5 with error 6.6928129704990624e-12 Converged in 5 steps
set_eval_nodes!(res, grids...)
V, X, resid = res.V, res.X, res.resid;
title = "Optimal Production Policy"
xlabel = "Demand Shock"
ylabel = "Lagged Production"
zlabel = "Production"
zlims = (0, 4)
PyPlot.surf(grids..., permutedims(reshape(X, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_zlim](zlims)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Value Function"
zlabel = "Value"
PyPlot.surf(grids..., permutedims(reshape(V, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
title = "Approximation Residual"
zlabel = "Residual"
PyPlot.surf(grids..., permutedims(reshape(resid, grid_length...)), cmap=PyPlot.cm[:jet])
ax = PyPlot.gca()
ax[:set_xlabel](xlabel)
ax[:set_ylabel](ylabel)
ax[:set_zlabel](zlabel)
ax[:set_title](title, y=1.05)
ax[:view_init](ax[:elev], 230)
ax[:ticklabel_format](style="sci", axis="z", scilimits=(0,0))
n_yrs = 20
n_paths = 5000
s_paths = Array{Float64}(undef, ndims(cdp), n_yrs+1, n_paths)
q_init = 0.3
d_rv = DiscreteRV(cdp.weights)
for i in 1:n_paths
s_init = [cdp.shocks[rand(d_rv)], q_init]
s_paths[:, :, i] = simulate(res, s_init, n_yrs+1)
end
mean_q_path = mean(s_paths[2, :, :], dims=2);
title = "Expected Policy Path"
xlabel = "Year"
ylabel = "Production"
plot(0:n_yrs, mean_q_path, ylims=(0.2, 1.2),
title=title, xlabel=xlabel, ylabel=ylabel, label="")