My opinion:
Install from Julia REPL with Pkg.add("Interpolations")
In [1]:
using QuantEcon: tauchen
using CompEcon
using Interpolations
using PlotlyJS
Interpolation object is defined on the set of knot indexes - not the actual knot values
Construct interpolation object itp using
itp = interpolate(Y, options ...)
where Y is an Array of values of the function being interpolated, one for each data point (if the domain is multi-dimensional, a data point is a list of knots - one for each dimension)
Options:
BSpline(Constant()) degree 0BSpline(Linear()) degree 1BSpline(Quadratic(x())) degree 2, must specify the boundary condition: x can be Flat, Natural, Free, Periodic, ReflectBSpline(Cubic(x())) degree 3, must specify the boundary condition: x can be Flat, Natural, Free, PeriodicNoInterp() just lookup of interpolated function value at the knot index - useful if state space is finite in one dimension (e.g. Markov chain)
(BSpline(Linear()), BSpline(Quadratic(Flat())))OnGrid() if data points in Y lie on the boundaries of the interpolation intervalOnCell() if lie on half-intervals between boundariesMore natural to define the interpolant on the domain of the interpolated function: use the scale function
In 1 dimension,
itp_scaled = scale(itp, x)
where x is the vector of knots
In n dimensions,
itp_scaled = scale(itp, x1, ..., xn)
where xi is the vector of knots in the i-th dimension
To evaluate the intepolant at position z use
v = itp[z]
where z is, in general, a tuple of numbers - one for each domain dimension - and itp can be a standard or scaled interpolation object. If itp is standard interpolant, z is grid position; if itp is scaled interpolant, z is point in the domain
Example: 1 dimension
In [2]:
xmin, xmax = -2, 2.
x = linspace(xmin, xmax, 5)
Out[2]:
In [3]:
y = collect(x).^2
Out[3]:
In [4]:
itp1 = interpolate(y, BSpline(Linear()), OnGrid())
Out[4]:
In [5]:
itp1_scaled = scale(itp1, x)
Out[5]:
In [6]:
function plot_itp_1d(itp, itp_label)
pos = collect(linspace(xmin, xmax, 100))
tr_itp = scatter(; x = pos, y = [itp[p] for p in pos], name = itp_label)
tr_actual = scatter(; x = pos, y = pos.^2, name = "Actual")
tstring = string("Inspect ", itp_label, " Approximation")
plot([tr_itp, tr_actual], Layout(title = tstring, xaxis_title = "x value", yaxis_title = "f(x)"))
end
Out[6]:
In [7]:
plot_itp_1d(itp1_scaled, "B-Spline deg. 1")
Out[7]:
In [8]:
itp2 = interpolate(y, BSpline(Quadratic(Free())), OnGrid())
Out[8]:
In [9]:
itp2_scaled = scale(itp2, x)
plot_itp_1d(itp2_scaled, "B-Spline deg. 2")
Out[9]:
In [10]:
itp3 = interpolate(y, BSpline(Cubic(Natural())), OnGrid())
itp3_scaled = scale(itp3, x)
plot_itp_1d(itp3_scaled, "B-Spline deg. 3")
Out[10]:
In [11]:
itp0 = interpolate(y, BSpline(Constant()), OnGrid())
itp0_scaled = scale(itp0, x)
plot_itp_1d(itp0_scaled, "B-Spline deg. 0")
Out[11]:
Example: 2 dimensions
In [12]:
zmin, zmax = 0., 4.
z = linspace(zmin, zmax, 5)
f(x, z) = x.^2 + log(1 + z)
Y = Array(Float64, length(x), length(z))
for j in 1:length(z), i in 1:length(x)
Y[i, j] = f(x[i], z[j])
end
Y
Out[12]:
In [13]:
itp1 = interpolate(Y, BSpline(Linear()), OnGrid())
Out[13]:
In [14]:
itp1_scaled = scale(itp1, x, z)
Out[14]:
In [15]:
function plot_itp_2d(itp, xslice, zslice, itp_labelx, itp_labelz)
posx = collect(linspace(xmin, xmax, 100))
posz = collect(linspace(zmin, zmax, 100))
trx_itp = scatter(; x = posx, y = [itp[px, zslice] for px in posx], name = itp_labelx)
trz_itp = scatter(; x = posz, y = [itp[xslice, pz] for pz in posz], name = itp_labelz)
trx_actual = scatter(; x = posx, y = f(posx, zslice), name = "Actual")
trz_actual = scatter(; x = posz, y = f(xslice, posz), name = "Actual")
tstringx = string("Inspect ", itp_labelx, " Approximation Along x-Dimension")
tstringz = string("Inspect ", itp_labelz, " Approximation Along z-Dimension")
plotx = plot([trx_itp, trx_actual], Layout(title = tstringx, xaxis_title = "x value",
yaxis_title = "f(x, $zslice)"))
plotz = plot([trz_itp, trz_actual], Layout(title = tstringz, xaxis_title = "z value",
yaxis_title = "f($xslice, z)"))
[plotx; plotz]
end
Out[15]:
In [16]:
plot_itp_2d(itp1_scaled, 0., 0., "B-Spline deg. 1", "B-Spline deg. 1")
Out[16]:
In [17]:
itp2_1 = interpolate(Y, (BSpline(Quadratic(Free())), BSpline(Linear())), OnGrid())
itp2_1scaled = scale(itp2_1, x, z)
plot_itp_2d(itp2_1scaled, 0., 0., "B-Spline deg. 2", "B-Spline deg. 1")
Out[17]:
Same constructor name, different method:
itp = interpolate(knots, Y, options ...)
where knots is a tuple of vectors (knots1, knots2, ...) specifying for each dimension the knots where the function is interpolated (if only 1 dimension, (knots1, ))
Options:
Gridded(Constant()) degree 0Gridded(Linear()) degree 1NoInterp()OnGridExample
In [18]:
xmin, xmax = 1., 3.
knotsx = linspace(0, (xmax-xmin)^.4, 5).^(1/.4) + xmin
Out[18]:
In [19]:
y = log(knotsx)
Out[19]:
In [20]:
itpgd = interpolate((knotsx,), y, Gridded(Linear()))
Out[20]:
In [21]:
pos = linspace(xmin, xmax, 100)
tr_itp = scatter(; x = pos, y = [itpgd[p] for p in pos], name = "B-Spline deg. 1")
tr_actual = scatter(; x = pos, y = log(pos), name = "Actual")
tstring = string("Inspect Gridded B-Spline deg. 1 Approximation")
plot([tr_itp, tr_actual], Layout(title = tstring, xaxis_title = "x value", yaxis_title = "f(x)"))
Out[21]:
g = gradient(itp, x, y ...) or gradient!(g, itp, x, y ...) (g pre-allocated vector) to compute the gradient at positions x, y, ...extrapolate(itp, x()) where x can be Throw, Flat, Linear, Periodic, Reflect to create a new interpolation object based on itp, with user-defined behavior outside of the interpolation intervalExample 1: high-dimensional interpolation
In [33]:
# function that evaluates itp at all data points (= all possible combinations of knots)
function ongrid!(dest, itp)
for I in CartesianRange(size(itp))
dest[I] = itp[I]
end
end
# function that creates a collection of 20 knots with random endpoints
new_k(n=20) = linspace(-rand(), rand(), 20)
# create knots in 4 dimensions
knots1, knots2, knots3, knots4 = [new_k() for i=1:4]
function compare_to_ce(Y, k1, k2, k3, k4)
# create itp object given a 4d Array Y
itp = interpolate(Y, BSpline(Linear()), OnGrid())
# get scaled itp object using knots in the 4 dimensions
itp_scaled = scale(itp, k1, k2, k3, k4)
# allocate array with same shape and content as Y, which will be replaced by evaluation of the interpolant
# at all data ponts
vals = similar(Y)
# evaluate
ongrid!(vals, itp)
end
A = rand(20, 20, 20, 20)
@time compare_to_ce(A, knots1, knots2, knots3, knots4)
In [32]:
# function that creates a linear spline with 20 knots between random endpoints
new_p(n=20) = LinParams(n, -rand(), rand())
# create basis object with a 20-knot linear spline for each of the 4 dimensions
basis = Basis([new_p() for i=1:4]...)
# get all data points in a matrix (with 20 knots per dimension and 4d, it's a 20^4 by 4 matrix)
knots = nodes(basis)[1]
function compare_to_itp(basis, knots, Y)
# get interpolation coefficients and basis matrix structure (field vals is a tuple of basis matrixes,
# one for each dimension)
coeffs, bmat_struct = funfitxy(basis, knots, Y)
# evaluate
funeval(coeffs, bmat_struct)
end
A = randn(size(knots, 1))
@time compare_to_itp(basis, knots, A);
Example 2: income fluctuation problem
An agent is subject to a stochastic, persistent income process $y_t$ and every period can decide how much to save and consume. The current stock of savings is denoted $a_t$ and consumption is denoted $c_t$. Savings earn a net interest rate $r$, constant and exogenous.
The agent enjoys utility $u(c)$ from consumption.
The agent is not allowed to borrow, thus $a'\geq 0$, and consumption cannot be negative ($c\geq 0)$.
The states are the current stock of assets $a$ and the current realization of the income process $y$. The value function for the agent's problem is therefore
\begin{equation} v(a, y) = \max_{c\in [0, (1+r)a + y]}\left(u(c) + \beta \mathbb{E}_{y'}\left[v(a', y')\right]\right) \end{equation}where
\begin{equation} a' = y + (1+r)a - c, a' \geq 0 \text{ which is imposed through }c\leq (1+r)a + y, \beta \text{ is the discount factor} \end{equation}I assume CRRA utility function:
\begin{equation} u(c) = \frac{c^{(1-\gamma)}}{1-\gamma}\text{, where }\gamma \text{ is the risk-aversion parameter} \end{equation}The income process follows a Markov chain constructed from an AR(1) with persistence $\rho$ and standard deviation $\sigma$ using the Tauchen method.
The Euler equation for this problem is
\begin{equation} c^{-\gamma} \geq \beta (1+r) \sum_{y'\in S}p(y, y') c(a', y')^{-\gamma} \end{equation}
In [34]:
# type that holds all primitives and grid elements shared by the approximation methods
type IncomeFluct
# Model parameters
beta::Float64 # discount factor
risk_av::Float64 # risk-aversion parameter
r::Float64 # net interest rate
rho::Float64 # income process persistence
sigma::Float64 # income process standard deviation
# Grid parameters assets
na::Int64 # number of knots for the asset space
amin::Float64 # lower bound of the asset space (borrowing limit)
amax::Float64 # upper bound of the asset space
knotsa::Vector{Float64} # asset grid
# Grid parameters income
ny::Int64 # number of states for the income process
knotsy::Vector{Float64} # income values corresponding to the states
Py::Matrix{Float64} # stochastic matrix for the Markov chain
end
# type constructor which takes parameter values as keyword arguments
function IncomeFluct(; beta = .95, risk_av = 2., r = .02, rho = .9, sigma = sqrt(.06), na = 100, amin = 0.,
amax = 150., ny = 5)
# create an unequally-spaced grid for assets
knotsa = linspace(0, (amax-amin)^.4, na).^(1/.4) + amin
# Tauchen discretization of the income process (using QuantEcon tauchen function)
d_proc = tauchen(ny, rho, sigma)
knotsy = exp(d_proc.state_values)
Py = d_proc.p
# return the type which holds all elements
IncomeFluct(beta, risk_av, r, rho, sigma, na, amin, amax, knotsa, ny, knotsy, Py)
end
Out[34]:
Endogenous grid method using Interpolations.jl: policy function iteration exploiting the Euler equation to avoid using non-linear solver
In [35]:
# function that computes the approximate solution to the agent's problem
function solve_agent(; maxiter = 10000, tol = 1e-8)
# store the model and grid objects
holder = IncomeFluct()
# unpack all elements needed to solve the problem
beta, risk_av, r, na, ny, amin, knotsa, knotsy, Py = holder.beta, holder.risk_av, holder.r, holder.na, holder.ny,
holder.amin, holder.knotsa, holder.knotsy, holder.Py
# separate a-values and y-values from the Cartesian product of knotsa and knotsy
mesha = repmat(knotsa'', 1, ny)
meshy = repmat(knotsy', na, 1)
# initial guess for the consumption policy function
c = r * mesha + meshy
# create array to store the consumption policy function in the next algorithm iteration
cnext = similar(c)
@time begin
for n in 1:maxiter # count the algorithm iteration number
# given the guess, invert the Euler equation to obtain c(a',y)
c_Euler = (beta * (1 + r) * c.^(-risk_av) * Py').^(-1/risk_av)
# use the budget constraint to obtain a(a',y) [this is the endogenous grid]
grida_end = 1/(1+r) * (c_Euler + mesha - meshy)
# fill in the policy function for the next iteration (i.e. interpolate c_Euler on knotsa)
for j in 1:ny, i in 1:na
# take the endogenous grid for each income state
knotsa_end = grida_end[:, j]
# pull out a(a'=amin,y) i.e. the current value of assets which induces the constraint to bind
min_knotsa_end = knotsa_end[1]
# create interpolant
itpa_end = interpolate((knotsa_end,), c_Euler[:, j], Gridded(Linear()))
# pull out knots where we want to evaluate the consumption policy function
a, y = knotsa[i], knotsy[j]
if a < min_knotsa_end # if the knot where we evaluate the consumption policy function is
# below the minimum level which makes the constraint bind,
# then the constraint necessarily binds and consumption comes
# from the budget constraint
cnext[i, j] = (1+r) * a + y - amin
else
cnext[i, j] = itpa_end[a]
end
end
err = maxabs(cnext - c)
copy!(c, cnext)
if n % 50 == 0; @printf("Iter = %d, Dist = %1.2e\n", n, err); end
if err < tol; break; end
end
end
# get the asset policy function once the consumption policy function has been solved
aprime_star = meshy + (1+r)*mesha - c
# solve for the value function by iterating on the Bellman equation, given the policy we solved for
U = c.^(1-risk_av)/(1-risk_av)
V_star = copy(U)
cV = similar(V_star)
for n in 1:maxiter
# since the asset policy function maps from the knots to asset values which are not on the knots,
# we have to interpolate the value function we're solving for (which is also interpolated
# on the original knots)
for j in 1:ny
itpcV = interpolate((knotsa, ), V_star[:,j], Gridded(Linear()))
cV[:,j] = itpcV[aprime_star[:,j]]
end
V_next = U + beta * (cV * Py')
errV = maxabs(V_next - V_star)
copy!(V_star, V_next)
if errV < tol; break; end
end
aprime_star, V_star, c
end
Out[35]:
In [38]:
@time aprime_star, V_star, c_star = solve_agent();
Newton Value-function iteration using CompEcon.jl
In [39]:
# type that holds all the CompEcon.jl objects needed in the algorithm
type CEObjects
ba::Basis{1} # basis for the asset space
bmaty::SparseMatrixCSC # basis matrix evaluated on the knots for income
bmat_ay::SparseMatrixCSC # basis matrix evaluated on the Cartesian product of knotsa and knotsy
nS::Int64 # total number of data points (na*ny)
cash::Vector{Float64} # cash-on-hand value at all data points (y + (1+r)a)
coeffs_V::Vector{Float64} # interpolation coefficients for the value function
coeffs_cV::Vector{Float64} # interpolation coefficients for the continuation value
end
In [40]:
# utility function defined based on cash-on-hand
function u(aprime, holder, holder_ce)
risk_av = holder.risk_av
cash = holder_ce.cash
(cash - aprime).^(1-risk_av)/(1-risk_av)
end
Out[40]:
In [41]:
# function that solves for the optimal savings (a') at each data point using golden search
function solve_aprime_ce(holder, holder_ce)
# unpack the objects that we need
amin, amax, r = holder.amin, holder.amax, holder.r
nS, cash = holder_ce.nS, holder_ce.cash
# set the lower bound on savings (in this case, the borrowing constraint)
lb = fill(amin, nS)
function bmat_aprimey(arg, holder_ce)
# unpack
ba, bmaty = holder_ce.ba, holder_ce.bmaty
# compute the basis matrix on the vector of a' values (one for each data point) that we are considering
bmatarg = BasisStructure(ba, Direct(), arg).vals[1]
# get the basis matrix on the Cartesian product of a' values and knotsy
row_kron(bmaty, bmatarg)
end
function obj(arg, holder, holder_ce)
# unpack
beta = holder.beta
coeffs_cV = holder_ce.coeffs_cV # these are the interpolation coefficients for the continuation value
# compute the rhs of the Bellman equation given the vector of a' values
vec(u(arg, holder, holder_ce) + beta*bmat_aprimey(arg, holder_ce)*coeffs_cV)
end
# use vectorized golden method to solve for the optional savings
obj_gm(x) = obj(x, holder, holder_ce)
aprime, V = golden_method(obj_gm, lb, min(cash, amax))
# return the objects we want (optimal savings, rhs of the Bellman equation, basis matrix on
# optimal savings*knotsy)
aprime, V, bmat_aprimey(aprime, holder_ce)
end
Out[41]:
In [42]:
# function that computes the approximate solution to the agent's problem
function solve_agent_ce(; ordera = 1, maxiter = 20, tol = 1e-8) # I also allow to set the B-spline order
# on the asset space
# store the model and grid objects
holder = IncomeFluct()
# unpack
beta, r, knotsa, knotsy, Py = holder.beta, holder.r, holder.knotsa, holder.knotsy, holder.Py
# create the 2d basis
b = Basis(SplineParams(knotsa, 0, ordera), LinParams(knotsy, 0))
# pull out the basis for the asset space (to be stored later)
ba = b[1]
# get all data points (saved in S)
S, kn = nodes(b)
# compute basis matrixes on all data points, separately for a and y
bmata, bmaty = BasisStructure(b, Direct(), S).vals
# compute the basis matrix on all data points
bmat_ay = row_kron(bmaty, bmata)
# get number of data points
nS = size(S, 1)
# get actual knots number for the asset space (if B-spline order > 1, CompEcon.jl adds knots)
na = ba.n[1]
# separate a-values and y-values from the data points
Sa = S[:, 1]
Sy = S[:, 2]
# compute cash-on-hand on all data points
cash = Sy + (1+r)*Sa
# manipulate the stochastic matrix to fit the matrix representation of data points
Pexp = kron(Py, eye(na))
# store the CompEcon object that we need inside other functions
holder_ce = CEObjects(ba, bmaty, bmat_ay, nS, cash, zeros(nS), zeros(nS))
@time begin
for i in 1:maxiter
# pull out interpolation coefficients and stack them in a single vector to use Newton method
coeffs_V, coeffs_cV = holder_ce.coeffs_V, holder_ce.coeffs_cV
coeffs = [coeffs_V; coeffs_cV]
# given the current guess of the interpolation coefficients (i.e. of the value function and continuation
# value on the data points), solve for the optimal savings on all data points
aprime, V, bmat_aprimey = solve_aprime_ce(holder, holder_ce)
# create Jacobian for the linear system in the coefficients
jacobian = [bmat_ay -beta*bmat_aprimey;-Pexp*bmat_ay bmat_ay]
# compute distance between current guess of the value function at the data points and
# the rhs of the Bellman equation
updV = bmat_ay * coeffs_V - V
# compute distance between current guess of the continuation value at the data points and
# the continuation value implied by the current guess of the value function
updcV = bmat_ay * coeffs_cV - Pexp * bmat_ay * coeffs_V
# stack distances
upd = [updV; updcV]
# update coefficients using the Newton method
coeffs_next = coeffs - jacobian \ upd
err = maxabs(coeffs_next - coeffs)
copy!(holder_ce.coeffs_V, coeffs_next[1:nS])
copy!(holder_ce.coeffs_cV, coeffs_next[nS+1:end])
@printf("Iter = %d, Dist = %1.2e\n", i, err)
if err < tol; break; end
end
end
# once the value function has been solved for, get the policy functions
aprime_star, V_star, bmat_aprimey = solve_aprime_ce(holder, holder_ce)
c_star = cash - aprime_star
aprime_star, V_star, c_star, holder_ce
end
Out[42]:
In [46]:
@time aprime_star_ce, V_star_ce, c_star_ce, holder_ce = solve_agent_ce();
In [47]:
function plots(obj, obj_ce, holder_ce, tit)
# store the model and grid objects
holder = IncomeFluct()
# unpack
ny, knotsa, knotsy = holder.ny, holder.knotsa, holder.knotsy
# unpack the CompEcon objects we need
ba = holder_ce.ba
knotsa_ce = nodes(ba)[1] # get actual knots on the asset space (if B-spline order > 1, CompEcon.jl adds knots)
na = length(knotsa_ce) # knots number
# reshape the object we're plotting to have asset knots on the rows and income knots on the column
obj_ce = reshape(obj_ce, na, ny)
# create traces for Plotly.JS to store all lines
tr = GenericTrace[]
tr_ce = GenericTrace[]
# set initial value to compute distance between Interpolations.jl and CompEcon.jl solutions
maxdiff = 0.
colors = ["blue", "orange", "green", "red", "violet"]
for i in 1:ny
itp = interpolate((knotsa,), obj[:,i], Gridded(Linear())) # if B-spline order > 1, CompEcon.jl has more knots
# so have to interpolate the solution from Interpolations.jl on those points
maxdiff = max(maxabs(itp[knotsa_ce] - obj_ce[:,i]), maxdiff) # iterative compute distance
push!(tr, scatter(; x = knotsa_ce, y = itp[knotsa_ce], name = "itp, y = $(round(knotsy[i], 2))",
line_color = colors[i])) # store line for Interpolations.jl
push!(tr_ce, scatter(; x = knotsa_ce, y = obj_ce[:,i], name = "ce, y = $(round(knotsy[i], 2))",
line_color = colors[i], line_dash="dash")) # store line for CompEcon.jl
end
println(tit)
@printf("Maximum difference = %1.4e\n", maxdiff)
@printf("\n")
plot([tr; tr_ce], Layout(title = tit, xaxis_title = "Assets level"))
end
# create plot objects and compute distances between solutions
pV = plots(V_star, V_star_ce, holder_ce, "Value Function")
paprime = plots(aprime_star, aprime_star_ce, holder_ce, "Assets Policy Function")
pcons = plots(c_star, c_star_ce, holder_ce, "Consumption Policy Function")
;
In [48]:
pV
Out[48]:
In [49]:
paprime
Out[49]:
In [50]:
pcons
Out[50]:
In [ ]: