The CompEcon toolbox allows us to compute interpolations for:
Think on 3 main "theoretical" categories:
Basis: Family of basis functions; Domain; Interpolation Nodes [see basis.jl].BasisStructure representation: it evaluates the basis functions at the desired interpolation nodes [see basis_structure.jl].Basis to the real line [obtained by solving linear system of equations].BasisParams: abstract type whose fields are all info needed to construct univariate basis.type ChebParams <: BasisParams
n::Int
a::Float64
b::Float64
end
type SplineParams <: BasisParams
breaks::Vector{Float64}
evennum::Int
k::Int
end
type LinParams <: BasisParams
breaks::Vector{Float64}
evennum::Int
end
Basis{N}:type Basis{N, BF<:BasisFamily, BP<:BasisParams}
basistype::Vector{BF} # Basis family
n::Vector{Int} # number of points and/or basis functions
a::Vector{Float64} # lower bound of domain
b::Vector{Float64} # upper bound of domain
params::Vector{BP} # params to construct basis
end
Finally note that inside this file there are two additional useful functions:
nodes: given the chosen Basis, it computes the nodes.Basis: used to define the basis and also to create multidimensional basis.Group 2: type to represent the BasisStructure in basis_structure.jl.
AbstractBasisStructureRep [ABSR]: it groups the types of representation (Tensor, Direct, Expanded).An important function is BasisStructure.
basis, the type of representation (Tensor, etc) and the nodes.Another interesting function is Base.convert which converts from a Tensor or Direct BasisStructure to an Expanded one.
Φ_direct = BasisStructure(basis, CompEcon.Direct(),snodes,0)
Φ = convert(Expanded, Φ_direct, [0 0]).vals[1]
using CompEcon
f(x) = exp(-x)
a,b = -1.0,1.0
Option 1: Using funfitf
n = 10
basis_c = Basis(Cheb, n, a, b)
c_c = funfitf(basis_c, f)
Option 2: Using funfitxy
xgrid = collect(linspace(a,b,n))
basis_s = Basis(Spline, xgrid, 0, 1)
x_s = nodes(basis_s)[1]
y_s = f(x_s)
c_s = funfitxy(basis_s, x_s, y_s)[1]
Option 3: Using BasisStructure
x_c = nodes(basis_c)[1]
y_c = f(x_c)
phi_c = BasisStructure(basis_c).vals[1]
c_c2 = phi_c\y_c
In [1]:
using CompEcon
f(x) = exp(-x)
# Set the endpoints of approximation interval:
a = -1.0 # left endpoint
b = 1.0 # right endpoint
n = 10 # order of approximation
Out[1]:
In [2]:
# Option 1: (with Cheb)
basis_c = Basis(Cheb, n, a, b)
c_c = funfitf(basis_c, f)
# Option 2: (with Spline)
xgrid = collect(linspace(a,b,n))
basis_s = Basis(Spline, xgrid, 0, 1)
x_s = nodes(basis_s)[1]
y_s = f(x_s)
c_s = funfitxy(basis_s, x_s, y_s)[1]
Out[2]:
In [3]:
# Option 3: (with Cheb again)
x_c = nodes(basis_c)[1]
y_c = f(x_c)
phi_c = BasisStructure(basis_c, CompEcon.Expanded(),x_c).vals[1]
c_c2 = phi_c\y_c
# Compare different approaches
print("Coefficients from options 1 and 3")
[c_c c_c2]
Out[3]:
In [4]:
# Having computed the basis coefficients, one may now evaluate the
# approximant at any point x using funeval:
xvals = collect(linspace(-1.0,1.0,50))
y_c = funeval(c_c, basis_c, xvals)
y_s = funeval(c_s, basis_s, xvals)
yy = f(xvals)
yvals = [y_c y_s yy]
Method = ["Cheb", "Spline", "True Function"]
using PyPlot
fig, ax = subplots()
for i=1:3
meth = Method[i]
ax[:plot](xvals, yvals[:,i], linewidth=2, alpha=0.6, label=L"$Method$ ="" $meth")
end
ax[:legend](loc=1)
Out[4]:
Let the "state" be given by $s=(a,y)$. I take the function $f(s) = \sqrt{a}\exp{(y)}$. To generate the interpolation, I use a linear spline for the first variable and a Chebichev polynomial for the second. Then I plot the true value and the interpolation, for different values of $s$.
In this example, we get a better approximation by mixing splines and Chebychev, than by just using Chebychev in both arguments.
In [5]:
fun2(st) = sqrt(st[:,1]).*exp(st[:,2])
agrid0 = collect(linspace(0.0,10.0,10))
a_basis = Basis(Spline, agrid0, 0, 1)
y_basis = Basis(Cheb, 10, 0.0, 10.0)
basis2 = Basis(a_basis, y_basis)
c_2d = funfitf(basis2, fun2)
Out[5]:
In [6]:
avals = collect(linspace(0.0,10.0,50))
yvals = collect(linspace(0.0,10.0,50))
svals = [avals yvals]
y_val = funeval(c_2d, basis2, svals)
y_true = fun2(svals)
yvals2d = [y_val y_true]
Out[6]:
In [7]:
Method = ["Interpolation", "True Function"]
using PyPlot
fig, ax = subplots()
for i=1:2
meth = Method[i]
ax[:plot](1:1:50, yvals2d[:,i], linewidth=2, alpha=0.6, label=L"$Method$ ="" $meth")
end
ax[:legend](loc=1)
Out[7]:
To approximate the solution, I follow Simon Mongey's notes given during Gianlucca Violante's course.
Define the set of collocation nodes $s = [\mathbf{1}_{Nz}\otimes{}\mathbf{X},\mathbf{Z}\otimes{}\mathbf{1}_{Nx}]$ and let $N=N_{x}\text{x}N_{z}$
We can define the following system: \begin{equation} V(s_{i}) = \max_{x'\in B(s_{i})}F(s_{i},x')+\beta V_{e}([x',s_{i,2}]) \\ V_{e}(s_{i}) = \sum_{k=1}^{Nz}P(z,z_{k}')V([s_{i,1},z'_{k}]) \end{equation}
If we substitute for the interpolants:
\begin{equation} \sum_{j=1}^{N}\phi(s_{i})c_{j}= \max_{x'\in B(s_{i})}F(s_{i},x')+\beta\sum_{j=1}^{N}\phi([x',s_{i,2}])c_{j}^{e} \\ \sum_{j=1}^{N}\phi(s_{i})c_{j}^{e}= \sum_{k=1}^{Nz}P(z,z_{k}')\sum_{j=1}^{N}\phi([s_{i,1},z'_{k}])c_{j} \end{equation}
... and write it in a stacked form, we get:
\begin{equation} \Phi(s) = \max_{x'\in\mathbf{B}(s)}F(s,x')+\beta\Phi([x',s_{2}])c^{e} \\ \Phi(s)c^{e} = (P\otimes{}\mathbf{I}_{Nx})\Phi(s)c \end{equation}
Direct specification.
In [1]:
using CompEcon
# Step 0: Define the Economic Growth model type EconGrowth
type EconGrowth
a::Float64
b::Float64
nx::Int64
m::Int64
σ::Float64
δ::Float64
α::Float64
β::Float64
γ::Float64
basis::Basis
bs::BasisStructure
Φ::SparseMatrixCSC{Float64,Int}
snodes::Matrix{Float64}
ϵ::Vector{Float64}
P::Matrix{Float64}
tol::Float64 # Tolerance Level.
maxit::Int # Maximum iterations.
end
# Step 1: Define the model's functions
f(eg::EconGrowth, x) = x.^eg.β
g(eg::EconGrowth, sn=eg.snodes) = eg.γ*sn[:,1] + sn[:,2].*f(eg, sn[:,1])
u(eg::EconGrowth, sn, xprime) = (g(eg, sn)-xprime).^(1-eg.α)./(1-eg.α)
# Step 2: Discretize stochastic process and build interpolation basis.
# Wrap up everything in the EconGrowth function.
function EconGrowth(;a::Float64=2.0,
b::Float64=8.0,
nx::Int64=10,
m::Int64=3,
σ::Float64=0.1,
δ::Float64=0.9,
α::Float64=2.0,
β::Float64=0.5,
γ::Float64=0.9,
tol::Float64=1e-9,
maxit::Int=10_000)
# Discretize stochastic process
ϵ, ω = qnwlogn(m, 0, σ^2) # quadrature for lognormal
P = repmat(ω', m, 1) # Since the process is i.i.d but I wanted to keep the general structure.
# Build interpolation basis
# x_params = SplineParams(nx, a, b, 1)
x_params = ChebParams(nx, a, b)
z_params = LinParams(ϵ, 0)
basis = Basis(x_params, z_params)
snodes, (xnodes, znodes) = nodes(basis)
bs = BasisStructure(basis, Direct(), snodes, [0 0])
Φ = convert(Expanded, bs).vals[1]
EconGrowth(a, b, nx, m, σ, δ, α, β, γ, basis, bs, Φ, snodes, ϵ, P, tol, maxit)
end
# Step 3: Compute the Steady State
function steady_state(eg::EconGrowth)
δ, γ, β, α = eg.δ, eg.γ, eg.β, eg.α
r_ss = 1/δ-1
x_ss = ((1/β)*(1/δ-γ))^(1/(β-1))
s_ss = γ*x_ss+x_ss^β
c_ss = s_ss-x_ss
lambda_ss = (c_ss)^(-α)
vals_ss = [r_ss,x_ss,s_ss,c_ss,lambda_ss]
end
Out[1]:
In [2]:
# Step 4: Write the Optimization Problem
function obj_fun(eg::EconGrowth, ce::Array{Float64,1})
function obj(xp)
Φ_xp = BasisStructure(eg.basis[1], Expanded(), xp, [0]).vals[1]
Φ = row_kron(eg.bs.vals[2], Φ_xp)
u(eg, eg.snodes, xp) + eg.δ*Φ*ce
end
x_sol,f_sol = golden_method(obj, zeros(size(eg.snodes, 1)), g(eg, eg.snodes))
end
# Step 5: Write down the iterative procedure (one-step)
# Bellman Iteration
function one_step(eg::EconGrowth, cc::Vector{Float64}, ce::Vector{Float64}, P_Φ)
xsol,fsol = obj_fun(eg, ce)
f2 = P_Φ*cc
cc_out = eg.Φ\fsol
ce_out = eg.Φ\f2
return cc_out, ce_out
end
Out[2]:
In [3]:
# Step 6: Define the iteration function
function iter(eg::EconGrowth)
iteration = 0
di = 1
n_x = length(nodes(eg.basis.params[1]))
P_Φ = kron(eg.P, eye(n_x))*eg.Φ
cc = zeros(size(eg.snodes, 1))
ce = zeros(size(eg.snodes, 1))
while di>eg.tol
iteration += 1
if iteration > eg.maxit
break
else
cc0 = copy(cc)
ce0 = copy(ce)
cc,ce = one_step(eg, cc0, ce0, P_Φ)
resid1 = norm(cc0-cc)
resid2 = norm(ce0-ce)
di = max(resid1,resid2)
# @printf("Iteration %d with distance %.3f\n", iteration, di)
end
end
@printf("The total number of interations was %d.\n", iteration)
return cc, ce
end
Out[3]:
In [4]:
# Computation of the Steady State
eg = EconGrowth()
vals_ss = steady_state(eg)
@printf("The steady values are r = %.2f, x = %.2f, s = %.2f, c = %.2f.\n"
,vals_ss[1],vals_ss[2],vals_ss[3],vals_ss[4])
# Solution to the Model
cc, ce = iter(eg)
println("The vector of coefficients cc and ce are:\n $([cc ce])")
In [8]:
using PyPlot
xsol,fsol = obj_fun(eg,ce)
s = g(eg,eg.snodes)
yplot = (xsol)./s
plot(s[2*eg.nx+1:3*eg.nx],yplot[2*eg.nx+1:3*eg.nx])
xlabel("Wealth"), ylabel("Investment in % of Wealth")
Out[8]:
In [6]:
x_params = ChebParams(100, eg.a, eg.b)
z_params = LinParams(eg.ϵ, 0)
basis2 = Basis(x_params, z_params)
snodes2, (xnodes2, znodes2) = nodes(basis2)
bs2 = BasisStructure(basis2, Direct(), snodes2, [0 0])
function obj2(xp)
Φ_xp2 = BasisStructure(eg.basis[1], Expanded(), xp, [0]).vals[1]
Φ2 = row_kron(bs2.vals[2], Φ_xp2)
u(eg, snodes2, xp) + eg.δ*Φ2*ce
end
x_sol,f_sol = golden_method(obj2, zeros(size(snodes2, 1)), g(eg, snodes2))
Out[6]:
In [9]:
s2 = g(eg,snodes2)
yplot2 = (x_sol)./s2
plot(s2[201:300],yplot2[201:300])
xlabel("Wealth"), ylabel("Investment in % of Wealth")
Out[9]:
In [ ]: