In [1]:
using QuantEcon
using BasisMatrices
using ContinuousDPs
using Plots
import PyPlot
using LinearAlgebra

Optimal Economic Growth

In [2]:
n = 10
s_min, s_max = 5, 10
basis = Basis(ChebParams(n, s_min, s_max))
Out[2]:
1 dimensional Basis on the hypercube formed by (5.0,) × (10.0,).
Basis families are Cheb
In [3]:
alpha = 0.2
bet = 0.5
gamm = 0.9
sigma = 0.1
discount = 0.9;
In [4]:
x_star = ((discount * bet) / (1 - discount * gamm))^(1 / (1 - bet))
s_star = gamm * x_star + x_star^bet
s_star, x_star
Out[4]:
(7.416897506925212, 5.6094182825484795)
In [5]:
f(s, x) = (s - x)^(1 - alpha) / (1 - alpha)
g(s, x, e) = gamm * x .+ e * x^bet;
In [6]:
n_shocks = 3
shocks, weights = qnwlogn(n_shocks, -sigma^2/2, sigma^2)  # See Errata
Out[6]:
([0.8367708003019956, 0.9950124791926823, 1.1831792330610165], [0.16666666666666674, 0.6666666666666666, 0.16666666666666674])
In [7]:
x_lb(s) = 0
x_ub(s) = 0.99 * s;
In [8]:
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
Out[8]:
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)))

Value iteration

In [9]:
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
In [10]:
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;
In [11]:
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="")
Out[11]:
In [12]:
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), yformatter=:scientific,
     title=title, xlabel=xlabel, ylabel=ylabel, label="")
Out[12]:

Policy iteration

In [13]:
res = solve(cdp, PFI);
Compute iterate 6 with error 2.1316282072803006e-14
Converged in 6 steps
In [14]:
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
In [15]:
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="")
Out[15]:
In [16]:
title = "Approximation Residual"
ylabel = "Residual"
plot(s_nodes, resid, xlims=(s_min, s_max), yformatter=:scientific,
     title=title, xlabel=xlabel, ylabel=ylabel, label="")
Out[16]:

Simulation

In [17]:
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);
In [18]:
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="")
Out[18]:

Nonrenewable Resource Management

In [19]:
k = 3
m = 101
breaks = m - (k-1)
s_min, s_max = 0, 10
basis = Basis(SplineParams(breaks, s_min, s_max, k))
Out[19]:
1 dimensional Basis on the hypercube formed by (0.0,) × (10.0,).
Basis families are Spline
In [20]:
a = [10, 0.8]
b = [12, 1.0]
discount = 0.9;
In [21]:
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;
In [22]:
shocks, weights = [0.], [1.];
In [23]:
x_lb(s) = 0
x_ub(s) = s;
In [24]:
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
Out[24]:
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)))

Value iteration

In [25]:
res = solve(cdp, VFI);
Compute iterate 38 with error 1.4533888759160618e-8
Converged in 38 steps
In [26]:
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;
In [27]:
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")
Out[27]:
In [28]:
B1 = evalbase(cdp.interp.basis.params[1], s_nodes, 1)
shadow_prices = B1 * res.C;
In [29]:
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")
Out[29]:
In [30]:
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="")
Out[30]:

Policy iteration

In [31]:
res = solve(cdp, PFI);
Compute iterate 6 with error 1.046872007484656e-13
Converged in 6 steps
In [32]:
set_eval_nodes!(res, s_nodes)
V, X, resid = res.V, res.X, res.resid;
In [33]:
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")
Out[33]:
In [34]:
B1 = evalbase(cdp.interp.basis.params[1], s_nodes, 1)
shadow_prices = B1 * res.C;
In [35]:
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")
Out[35]:
In [36]:
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="")
Out[36]:

Simulation

In [37]:
n_yrs = 20
s_init = 10.
s_path = simulate(res, s_init, n_yrs+1);
In [38]:
title = "State Path"
xlabel = "Year"
ylabel = "Stock"
plot(0:n_yrs, s_path, ylims=(2, 10),
     title=title, xlabel=xlabel, ylabel=ylabel, label="")
Out[38]:

Monetary Policy

In [39]:
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;
In [40]:
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;
In [41]:
n_shocks = [3, 3]
mu = [0, 0]
shocks, weights = qnwnorm(n_shocks, mu, Sigma)
Out[41]:
([-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])
In [42]:
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)...)
Out[42]:
2 dimensional Basis on the hypercube formed by (-15.0, -10.0) × (15.0, 10.0).
Basis families are Spline × Spline
In [43]:
cdp = ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub, basis)
Out[43]:
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)))

Value iteration

In [44]:
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
In [45]:
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;
In [46]:
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)
In [47]:
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)