Rdtq (version 0.1)

studydtqconv: Study DTQ Convergence

Description

studydtqconv facilitates the generation of convergence plots, i.e., plots where one studies the error (in various norms) as a function of the time step $h$; the error is computed as the difference between the exact PDF and the approximate PDF computed via the DTQ method.

Usage

studydtqconv(method, drift, diffusion, exact, hseq, kseq, Mseq, init, fT, thresh = 0)

Arguments

method

This must be a string, either "cpp" or "sparse", that indicates which algorithm to use. See the parameter of the same name in the rdtq function.

drift

if method="cpp", then this should be a pointer to a drift function implemented in C++ using Rcpp. If method="sparse", then this should be the name of an R function. For further details, see the description of the drift parameter for the rdtq function.

diffusion

if method="cpp", then this should be a pointer to a diffusion function implemented in C++ using Rcpp. If method="sparse", then this should be the name of an R function. For further details, see the description of the diffusion parameter for the rdtq function.

exact

an R function that accepts two arguments, xvec and T, and returns the exact probability density function $p(x,T)$ for each $x$ in xvec.

hseq

a numeric vector of values of $h$, the time step, to use for the computation of the DTQ solution. Note that hseq and kseq must have the same lengths.

kseq

a numeric vector of values of $k$, the grid spacing, to use for the computation of the DTQ solution. Note that hseq and kseq must have the same lengths.

Mseq

a numeric vector of integer values of $M$. For each corresponding value of $k$, the spatial grid will cover the domain $[-y_M, y_M]$ where $y_M = Mk$. This corresponds to the parameter bigm in the rdtq function. Note that kseq and Mseq must have the same lengths.

init

a scalar initial condition.

fT

a positive numeric scalar giving the final time at which to compare the exact and DTQ solutions.

thresh

This optional positive scalar is only used when method="cpp". See the parameter of the same name in the rdtq function.

Value

The function returns the errors between the DTQ and exact solutions indexed by the corresponding value of hseq. The errors are returned in the $L^1$ norm, $L^\infty$ norm, and the Kolmogorov-Smirnov norm. The errors are returned in the form of a data frame.

See Also

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

Examples

# In this example, we will study the convergence of the DTQ method
# for the SDE with drift f(x) = x/2 + (1 + x^2)^(1/2) and
# diffusion g(x) = (1 + x^2)^(1/2).

library(Rdtq)
library(Rcpp)

# implement the drift and diffusion functions using C++
sourceCpp(code = '#include <Rcpp.h>
          using namespace Rcpp;
          double drift(double& x) { return(0.5*x + sqrt(x*x + 1.0)); }
          double diff(double& x) { return(sqrt(x*x + 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))); }')

# implement the drift and diffusion functions using R
mydrift = function(y)
{
  return(0.5*y + sqrt(y^2 + 1))
}
mydiff = function(y)
{
  return(sqrt(y^2 + 1))
}

# implement the exact solution at time t, i.e.,
# the analytical formula for the pdf p(x,t)
exactsol = function(xvec,t)
{
  transx = asinh(xvec) - t
  prefac = (1 + xvec^2)^(-1/2)
  z = prefac*dnorm(x=transx)
  return(z)
}

# define the sequence of parameters that will be used to study convergence
hseq = c(0.5,0.2,0.1,0.05)
kseq = hseq^(0.55)
Mseq = ceiling(5*(-log(hseq))/kseq)

# we will use the method="sparse" code for the three largest values in hseq,
# and then switch to the method="cpp" code for the three smallest values
firstpart = c(1:2)
errpart1 = studydtqconv(method="sparse",drift=mydrift,diffusion=mydiff,exact=exactsol,
                        hseq[firstpart],kseq[firstpart],Mseq[firstpart],
                        init=0,fT=1)
errpart2 = studydtqconv(method="cpp",drift=driftXPtr(),diffusion=diffXPtr(),exact=exactsol,
                        hseq[-firstpart],kseq[-firstpart],Mseq[-firstpart],
                        init=0,fT=1,thresh=1e-16)

# now we will put everything together into one data frame
mydat = rbind(errpart1,errpart2)

# we plot the convergence diagram, on a log-log scale, using ggplot2
library(ggplot2)
library(scales)
myplot = ggplot(data=mydat, aes(x=x,y=y,group=norm,color=norm))
myplot = myplot + theme_bw() + theme(plot.background = element_rect(fill='white'))
myxticks = sort(10^(round(log(hseq)/log(10)*10)/10))
rawyticks = round(log(mydat$y)/log(10)*10)/10
rawyticks = round(seq(from=min(rawyticks),to=max(rawyticks),length.out=length(myxticks))*1)/1
myyticks = unique(10^rawyticks)
myplot = myplot + scale_x_log10(breaks = hseq)
myplot = myplot + theme(axis.text.x = element_text(angle=90,hjust=1))
myplot = myplot + scale_y_log10(breaks = myyticks,
                                labels = trans_format("log10", math_format(10^.x)))
myplot = myplot + labs(x="h (temporal step size)", y="error")
myplot = myplot + geom_line() + geom_point()

# save the plot to a pdf (portable document format) file
ggsave(filename="example.pdf", plot=myplot, width=5, height=4)