Interpolations.jl

Alberto Polo

April 8, 2016

My opinion:

  • The package is easy to learn even with basic knowledge of interpolation theory
  • Efficient to use: only one function to create interpolation objects, with a few methods
  • Fast ...
  • ... doing a lot of things under the hood allows to achieve those results
  • CompEcon in turn allows to manipulate the matrixes which underlie interpolation, which can be useful in designing certain algorithms
  • limited types and degrees of polynomials available, mostly with the irregularly-spaced mode

Install from Julia REPL with Pkg.add("Interpolations")


In [1]:
using QuantEcon: tauchen
using CompEcon
using Interpolations
using PlotlyJS


Plotly javascript loaded.

Mode 1: Equally-spaced knots

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 0
  • BSpline(Linear()) degree 1
  • BSpline(Quadratic(x())) degree 2, must specify the boundary condition: x can be Flat, Natural, Free, Periodic, Reflect
  • BSpline(Cubic(x())) degree 3, must specify the boundary condition: x can be Flat, Natural, Free, Periodic
  • NoInterp() just lookup of interpolated function value at the knot index - useful if state space is finite in one dimension (e.g. Markov chain)

    • If the domain is multi-dimensional, can specify different polynomial degrees in different dimensions by using a tuple, e.g. (BSpline(Linear()), BSpline(Quadratic(Flat())))
  • OnGrid() if data points in Y lie on the boundaries of the interpolation interval
  • OnCell() if lie on half-intervals between boundaries

More 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]:
linspace(-2.0,2.0,5)

In [3]:
y = collect(x).^2


Out[3]:
5-element Array{Float64,1}:
 4.0
 1.0
 0.0
 1.0
 4.0

In [4]:
itp1 = interpolate(y, BSpline(Linear()), OnGrid())


Out[4]:
5-element Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0}:
 4.0
 1.0
 0.0
 1.0
 4.0

In [5]:
itp1_scaled = scale(itp1, x)


Out[5]:
5-element Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{LinSpace{Float64}}}:
  1.0
  4.0
  7.0
 10.0
 13.0

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]:
plot_itp_1d (generic function with 1 method)

In [7]:
plot_itp_1d(itp1_scaled, "B-Spline deg. 1")


Out[7]:

In [8]:
itp2 = interpolate(y, BSpline(Quadratic(Free())), OnGrid())


Out[8]:
5-element Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.Free}},Interpolations.OnGrid,1}:
  4.0        
  1.0        
 -2.77556e-17
  1.0        
  4.0        

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]:
5x5 Array{Float64,2}:
 4.0  4.69315   5.09861  5.38629  5.60944
 1.0  1.69315   2.09861  2.38629  2.60944
 0.0  0.693147  1.09861  1.38629  1.60944
 1.0  1.69315   2.09861  2.38629  2.60944
 4.0  4.69315   5.09861  5.38629  5.60944

In [13]:
itp1 = interpolate(Y, BSpline(Linear()), OnGrid())


Out[13]:
5x5 Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0}:
 4.0  4.69315   5.09861  5.38629  5.60944
 1.0  1.69315   2.09861  2.38629  2.60944
 0.0  0.693147  1.09861  1.38629  1.60944
 1.0  1.69315   2.09861  2.38629  2.60944
 4.0  4.69315   5.09861  5.38629  5.60944

In [14]:
itp1_scaled = scale(itp1, x, z)


Out[14]:
5x5 Interpolations.ScaledInterpolation{Float64,2,Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{LinSpace{Float64},LinSpace{Float64}}}:
  1.69315   2.09861   2.38629   2.60944   2.83258
  4.69315   5.09861   5.38629   5.60944   5.83258
  7.69315   8.09861   8.38629   8.60944   8.83258
 10.6931   11.0986   11.3863   11.6094   11.8326 
 13.6931   14.0986   14.3863   14.6094   14.8326 

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]:
plot_itp_2d (generic function with 1 method)

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]:

Mode 2: irregularly-spaced knots (Gridded method)

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 0
  • Gridded(Linear()) degree 1
  • NoInterp()
  • Necessarily OnGrid

Example


In [18]:
xmin, xmax = 1., 3.
knotsx = linspace(0, (xmax-xmin)^.4,  5).^(1/.4) + xmin


Out[18]:
5-element Array{Float64,1}:
 1.0    
 1.0625 
 1.35355
 1.97428
 3.0    

In [19]:
y = log(knotsx)


Out[19]:
5-element Array{Float64,1}:
 0.0      
 0.0606246
 0.302733 
 0.680203 
 1.09861  

In [20]:
itpgd = interpolate((knotsx,), y, Gridded(Linear()))


Out[20]:
5-element Interpolations.GriddedInterpolation{Float64,1,Float64,Interpolations.Gridded{Interpolations.Linear},Tuple{Array{Float64,1}},0}:
 0.0     
 0.690695
 1.09861 
 1.50653 
 1.91445 

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]:

Other

  • 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 interval

Comparing Interpolations.jl and CompEcon.jl

Example 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)


  0.028287 seconds (21.44 k allocations: 3.362 MB)

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);


  0.618889 seconds (7.67 M allocations: 718.329 MB, 21.06% gc time)

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]:
IncomeFluct

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]:
solve_agent (generic function with 1 method)

In [38]:
@time aprime_star, V_star, c_star = solve_agent();


Iter = 50, Dist = 3.05e-02
Iter = 100, Dist = 6.25e-03
Iter = 150, Dist = 3.79e-04
Iter = 200, Dist = 4.53e-06
Iter = 250, Dist = 5.16e-08
  0.188479 seconds (1.49 M allocations: 375.945 MB, 13.32% gc time)
  0.196879 seconds (1.51 M allocations: 389.655 MB, 12.75% gc time)

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]:
u (generic function with 1 method)

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]:
solve_aprime_ce (generic function with 1 method)

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]:
solve_agent_ce (generic function with 1 method)

In [46]:
@time aprime_star_ce, V_star_ce, c_star_ce, holder_ce = solve_agent_ce();


Iter = 1, Dist = 5.41e+01
Iter = 2, Dist = 2.35e+01
Iter = 3, Dist = 1.93e+01
Iter = 4, Dist = 3.07e+00
Iter = 5, Dist = 5.55e-01
Iter = 6, Dist = 1.77e-01
Iter = 7, Dist = 4.37e-02
Iter = 8, Dist = 7.21e-03
Iter = 9, Dist = 6.61e-04
Iter = 10, Dist = 3.01e-06
Iter = 11, Dist = 2.47e-11
  0.932270 seconds (2.21 M allocations: 767.083 MB, 10.89% gc time)
  0.972177 seconds (2.42 M allocations: 816.660 MB, 11.51% gc time)

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")

;


Value Function
Maximum difference = 5.7997e-01

Assets Policy Function
Maximum difference = 5.1406e-02

Consumption Policy Function
Maximum difference = 5.1406e-02


In [48]:
pV


Out[48]:

In [49]:
paprime


Out[49]:

In [50]:
pcons


Out[50]:

In [ ]: