rdtq: Density Tracking by Quadrature
In Rdtq: Density Tracking by Quadrature

Description Usage Arguments Details Value See Also Examples

rdtq implements density tracking by quadrature (DTQ) algorithms to compute the probability density function of a stochastic differential equation with user-specified drift and diffusion functions.

1 2	rdtq(h, k = NULL, bigm, a = NULL, b = NULL, init, fT, drift = NULL, diffusion = NULL, thresh = 0, method = "sparse")

`h`	Time step size, a positive numeric scalar.
`k`	Spatial grid spacing, a positive numeric scalar. (Must be specified if `a` and `b` are not specified.)
`bigm`	If k is specified, then bigm is a positive integer such that -`bigmk` and `bigmk` are, respectively, the minimum and maximum grid points. If `a` and `b` are specified, then `bigm` is the total number of grid points. Note that the fractional part of `bigm` is ignored, and that `floor(bigm)` must be at least 2.
`a`	Left boundary, a numeric scalar. (Must be specified if `k` is not specified.)
`b`	Right boundary, a numeric scalar. (Must be specified if `k` is not specified.)
`init`	A numeric scalar indicating either a fixed initial condition of the form X(0)=`init`, or a numeric vector giving the PDF at time t=0. In the latter case, the vector must have the same size as the spatial grid.
`fT`	The final time, a positive numeric scalar. The computation assumes an initial time of t=0 and then computes the PDF at time t=`fT`.
`drift`	When the user chooses the `method="cpp"` algorithm, this should be a pointer to a drift function that is implemented in C++ using `Rcpp`. In our C++ code, we define the type `funcPtr` using the following code: `typedef double (*funcPtr)(const double& x);` We expect the drift function to be a C++ function, implemented using Rcpp, of type `XPtr<funcPtr>`. See the first example below. When the user chooses the `method="sparse"` algorithm, this should be an R function that takes as input a numeric vector of values. The function should return as output a numeric vector containing the drift function evaluated at each element of the input vector. See the second example below.
`diffusion`	When the user chooses the `method="cpp"` algorithm, this should be a pointer to a diffusion function that is implemented in C++ using `Rcpp`. All of the details are analogous to that of the `drift` function described above. When the user chooses the `method="sparse"` algorithm, this should be an R function that takes as input a numeric vector of values. The function should return as output a numeric vector containing the diffusion function evaluated at each element of the input vector. See the second example below.
`thresh`	This is an optional numeric scalar parameter that is only used for the `method="cpp"` algorithm. When the DTQ summand drops below codethresh, the algorithm stops summing, even if it has not summed over all grid points. The default value of this parameter is zero, indicating that the full DTQ sum is evaluated. Setting this parameter to a small positive value such as 2.2 \times 10^{-16} can result in a substantial speed up for computations on large spatial grids, especially when h is also small.
`method`	A string that indicates which DTQ algorithm to use. There are two choices: `"cpp"` This DTQ method is implemented in C++. No matrices are formed; the method is highly conservative in its usage of memory. For sufficiently small `h` and `k`, it is necessary to use this method. This method also allows for approximate evaluation of the DTQ algorithm by setting a positive threshold parameter. `"sparse"` This DTQ method is implemented in R using sparse matrices from the Matrix package. The method uses more memory than the `"cpp"` method, but may be faster for larger values of `h` and `k`. This is the default method.

Consider the stochastic differential equation (SDE)

dX(t) = f(X(t)) dt + g(X(t)) dW(t)

where W(t) is standard Brownian motion, f is the drift function, and g is the diffusion function. Let p(x,t) denote the time-dependent probability density function (PDF) of X(t); then rdtq computes p(x,T) for a fixed time T.

Note that the PDF is computed on a spatial grid that can be specified in one of two ways:

Specify a real, positive value k and a positive integer M = bigm. In this case, the PDF will be computed on the grid x_j = j k where j = -M, -M+1, ..., M-1, M. In total, there will be 2M+1 grid points.
Specify a real, positive integer M and a computational domain [a,b]. In this case, there will be exactly M equispaced grid points. The grid spacing will be k = (b-a)/(M-1).

The output consists of a list with two elements:

xvec: a numeric vector that contains the spatial grid
pdf: a numeric vector that contains the PDF evaluated at the grid points

H. S. Bhat and R. W. M. A. Madushani, "Density Tracking by Quadrature for Stochastic Differential Equations," arXiv:1610.09572 [stat.CO], http://bit.ly/2fbNsp5

# Example 1:
# Define the drift function f(x) = -x and diffusion function g(x) = 1
# using C++ code:
require(Rcpp)
sourceCpp(code = '#include <Rcpp.h>
using namespace Rcpp;
double drift(double& x)
{
  return(-x);
}
double diff(double& x)
{
  return(1.0);
}
typedef double (*funcPtr)(double& x);
// [[Rcpp::export]]
XPtr<funcPtr> driftXPtr()
{
  return(XPtr<funcPtr>(new funcPtr(&drift)));
}
// [[Rcpp::export]]
XPtr<funcPtr> diffXPtr()
{
  return(XPtr<funcPtr>(new funcPtr(&diff)));
}')
# Solve for the PDF (at final time fT=1) of the SDE with drift f,
# diffusion g, and deterministic initial condition X(0) = 0.
# First we solve using the grid specified by k and bigm.
# Then we solve using the grid specified by a, b, and bigm.
# We then check that we get the same PDF either way.
k = 0.01
M = 250
test1 = rdtq(h=0.1,k,bigm=M,init=0,fT=1,
             drift=driftXPtr(),diffusion=diffXPtr(),method="cpp")
test2 = rdtq(h=0.1,a=-2.5,b=2.5,bigm=501,init=0,fT=1,
             drift=driftXPtr(),diffusion=diffXPtr(),method="cpp")
print(k*sum(abs(test1$pdf-test2$pdf)))

# Example 2:
# We again use the drift function f(x) = -x and diffusion function g(x) = 1.
# This time, we use the method="sparse" version of DTQ.
# This requires us to define the drift and diffusion functions in R:
mydrift = function(x) { -x }
mydiff = function(x) { rep(1,length(x)) }
test = rdtq(h=0.1,k=0.01,bigm=250,init=0,fT=1,
            drift=mydrift,diffusion=mydiff,method="sparse")
plot(test$xvec,test$pdf,type='l')