Extend `simulate` to multidimensional states

[oyamad]

So far 1dim linear interpolation QuantEcon.LinInterp is used to interpolate the optimal policy function.

QuantEcon/QuantEcon.jl#190

[sglyon]

I think it would be great to extend QuantEcon.LinInterp -- I've received multiple requests to do so.

In the short run (before someone finds the time to make that happen), we could use LinParams from BasisMatrices

[oyamad]

I didn't think about LinParams.

Suppose I have a Matrix eval_nodes generated by gridmake(vec_of_1st_dim, vec_of_2nd_dim). Does LinParams accept eval_nodes to generate a 2dim interpolation?

[sglyon]

You use LinParams the same way you use ChebParams or SplineParams

Here's a 2d example:

lin_basis = Basis(LinParams(nodes_dim1, 0), LinParams(nodes_dim2, 0))

# get coefs from somewhere...

# option one
funeval(coefs, lin_basis, eval_nodes)

# option two (more efficient)
funeval(coefs, basis, (vec_of_1st_dim, vec_of_2nd_dim))

Both of the funeval calls should produce equivalent results. The benefit of the second is never forming the Direct or Expanded forms of the basis matrix when performing the tensor-product with coefs.

[sglyon]

One note about the comment to # get coefs from somewhere...is that if you are using only LinParams, then the basis matrix at the nodes (Phi throughout this code) will be the identity matrix, so coefs will exactly equal the function values at the nodes.

With that in mind, you can just set coefs = values and proceed to one of the two forms of evaluation shown after that comment

[oyamad]

Maybe I should change the implementation, but the current code only stores eval_nodes and does not keep vec_of_1st_dim and vec_of_2nd_dim, and the optimal policy X is defined on eval_nodes (i.e., the policy at state eval_nodes[i, :] is given by X[i]). Isn't it the case that internally, in the Expanded form, gridmake(vec_of_1st_dim, vec_of_2nd_dim) is used?

[sglyon]

So here's an example of how you can use LinParams:

function simulate_lin_params!(rng::AbstractRNG, s_path::TS,
                   res::ContinuousDPs.CDPSolveResult{Algo,N,TR,TS},
                   s_init) where {Algo,N,TR,TS<:VecOrMat}
    ts_length = size(s_path)[end]
    cdf = cumsum(res.cdp.weights)
    r = rand(rng, ts_length-1)
    e_ind = Array{Int}(ts_length-1)
    for t in 1:ts_length-1
        e_ind[t] = searchsortedlast(cdf, r[t]) + 1
    end

    basis = Basis(map(LinParams, res.cdp.interp.Scoord, ntuple(i -> 0, ndims(res.cdp))))
    X_interp = Interpoland(basis, res.X)

    s_ind_front = Base.front(indices(s_path))
    e_ind_tail = Base.tail(indices(res.cdp.shocks))
    s_path[(s_ind_front..., 1)...] = s_init
    for t in 1:ts_length-1
        s = s_path[(s_ind_front..., t)...]
        e = res.cdp.shocks[(e_ind[t], e_ind_tail...)...]
        s_path[(s_ind_front..., t+1)...] = res.cdp.g(s, X_interp(s), e)
    end

    return s_path
end

I've tested this using the monetary policy example from the MF notebook and the opt_Growth example from the other example notebook and it seems to work.

The only issue is that it assumes that res.X corresponds to res.cdp.interp.Scoord, which should be true once #10 is resolved.

[sglyon]

Note that when I say

res.X corresponds to res.cdp.interp.Scoord

I really should be more precise and say what you said:

(i.e., the policy at state eval_nodes[i, :] is given by X[i])

[sglyon]

Also to answer your question

Isn't it the case that internally, in the Expanded form, gridmake(vec_of_1st_dim, vec_of_2nd_dim) is used?

When you use Expanded form at matrix like eval_nodes you evaluate the dimension i basis functions at eval_nodes[:, i], the you use a row-wise kronecker product (row_kron) in reverse order to produce Phi.

If I had two dimensions in my basis, then using the Expanded form to evaluate at eval_nodes is equivalent to

Phi1 = evalbase(basis.params[1], eval_nodes[:, 1])
Phi2 =evalbase(basis.params[2], eval_nodes[:, 2]) 
Phi = row_kron(Phi2, Phi1)

[oyamad]

• 15 simulate now works also for multidimensional cases thanks to @sglyon https://github.com/QuantEcon/ContinuousDPs.jl/issues/4#issuecomment-345730733

• where res.cdp.interp.Scoord is replaced with res.eval_nodes_coord #14

[oyamad]

Done.