mlePde1D {sdetorus}

MLE for toroidal process via PDE solving in 1D

Description

Maximum Likelihood Estimation (MLE) for arbitrary diffusions, based on the transition probability density (tpd) obtained as the numerical solution of the Fokker–Planck Partial Differential Equation (PDE) in 1D.

Usage

mlePde1D(data, delta, b, sigma2, Mx = 500, Mt = ceiling(100 * delta),
  sdInitial = 0.1, linearBinning = FALSE, start, lower, upper, ...)

Arguments

data

a vector of size N with the discretized trajectory of the diffusion.

delta

time discretization step.

b

drift function. Must return a vector of the same size as its argument.

sigma2

function giving the squared diffusion coefficient. Must return a vector of the same size as its argument.

Mx

size of the equispaced spatial grid in [-\pi,\pi).

Mt

size of the time grid in [0, t].

sdInitial

the standard deviation of the concentrated WN giving the initial condition.

linearBinning

flag to indicate whether linear binning should be applied for the initial values of the tpd, instead of usual simple binning (cheaper). Linear binning is always done in the evaluation of the tpd.

start

starting values, a matrix with p columns, with each entry representing a different starting value.

lower, upper

bound for box constraints as in method "L-BFGS-B" of optim.

...

Further parameters passed to crankNicolson1D.

Details

See Sections 3.4.1 and 3.4.4 in García-Portugués et al. (2019) for details.

Value

Output from mleOptimWrapper.

References

García-Portugués, E., Sørensen, M., Mardia, K. V. and Hamelryck, T. (2019) Langevin diffusions on the torus: estimation and applications. Statistics and Computing, 29(2):1–22. doi:10.1007/s11222-017-9790-2

Examples

# Test in OU
alpha <- 2
mu <- 0
sigma <- 1
set.seed(234567)
traj <- rTrajOu(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500,
                delta = 0.5)
b <- function(x, pars) pars[1] * (pars[2] - x)
sigma2 <- function(x, pars) rep(pars[3]^2, length(x))

exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2),
                 lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10))
pdeOu <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b,
                  sigma2 = sigma2, start = c(1, 1, 2),
                  lower = c(0.1, -pi, -10), upper = c(10, pi, 10),
                  verbose = 2)
pdeOuLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b,
                     sigma2 = sigma2, start = c(1, 1, 2),
                     lower = c(0.1, -pi, -10), upper = c(10, pi, 10),
                     linearBinning = TRUE, verbose = 2)
head(exactOu)
head(pdeOu)
head(pdeOuLin)

# Test in WN diffusion
alpha <- 2
mu <- 0
sigma <- 1
set.seed(234567)
traj <- rTrajWn1D(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500,
                 delta = 0.5)

exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2),
                 lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10))
pdeWn <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100,
                  b = function(x, pars)
                    driftWn1D(x = x, alpha = pars[1], mu = pars[2],
                              sigma = pars[3]),
                  sigma2 = function(x, pars) rep(pars[3]^2, length(x)),
                  start = c(1, 1, 2), lower = c(0.1, -pi, -10),
                  upper = c(10, pi, 10), verbose = 2)
pdeWnLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100,
                     b = function(x, pars)
                       driftWn1D(x = x, alpha = pars[1], mu = pars[2],
                                 sigma = pars[3]),
                     sigma2 = function(x, pars) rep(pars[3]^2, length(x)),
                     start = c(1, 1, 2), lower = c(0.1, -pi, -10),
                     upper = c(10, pi, 10), linearBinning = TRUE,
                     verbose = 2)
head(exactOu)
head(pdeWn)
head(pdeWnLin)