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

## Description

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.

## Usage

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

## Arguments

 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 -bigm*k and bigm*k 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. 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.

## Details

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:

1. 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.

2. 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).

## Value

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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 # Example 1: # Define the drift function f(x) = -x and diffusion function g(x) = 1 # using C++ code: require(Rcpp) sourceCpp(code = '#include using namespace Rcpp; double drift(double& x) { return(-x); } double diff(double& x) { return(1.0); } typedef double (*funcPtr)(double& x); // [[Rcpp::export]] XPtr driftXPtr() { return(XPtr(new funcPtr(&drift))); } // [[Rcpp::export]] XPtr diffXPtr() { return(XPtr(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')