mlePde2D {sdetorus}

MLE for toroidal process via PDE solving in 2D

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

Usage

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

Arguments

data

a matrix of dimension c(n, p).

delta

discretization step.

b

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

sigma2

function giving the diagonal of the diffusion matrix. Must return a vector of the same size as its argument.

Mx, My

sizes of the equispaced spatial grids in [-\pi,\pi) for each component.

Mt

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

sdInitial

standard deviations 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 mleOptimWrapper.

Details

See Sections 3.4.2 and 3.4.4 in García-Portugués et al. (2019) for details. The function currently includes the region function for imposing a feasibility region on the parameters of the bivariate WN diffusion.

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 process
alpha <- c(1, 2, -0.5)
mu <- c(0, 0)
sigma <- c(0.5, 1)
set.seed(2334567)
data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = alpha, sigma = sigma),
                 mu = mu, Sigma = diag(sigma^2), N = 500, delta = 0.5)
b <- function(x, pars) sweep(-x, 2, pars[4:5], "+") %*%
                       t(alphaToA(alpha = pars[1:3], sigma = sigma))
sigma2 <- function(x, pars) repRow(sigma^2, nrow(x))

exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma,
                  start = c(1, 1, 0, 2, 2),
                  lower = c(0.1, 0.1, -25, -10, -10),
                  upper = c(25, 25, 25, 10, 10))
head(exactOu, 2)
pdeOu <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2,
                  Mx = 10, My = 10, Mt = 10,
                  start = rbind(c(1, 1, 0, 2, 2)),
                  lower = c(0.1, 0.1, -25, -10, -10),
                  upper = c(25, 25, 25, 10, 10), verbose = 2)
head(pdeOu, 2)
pdeOuLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2,
                     Mx = 10, My = 10, Mt = 10,
                     start = rbind(c(1, 1, 0, 2, 2)),
                     lower = c(0.1, 0.1, -25, -10, -10),
                     upper = c(25, 25, 25, 10, 10), verbose = 2,
                     linearBinning = TRUE)
head(pdeOuLin, 2)

# Test in WN diffusion
alpha <- c(1, 0.5, 0.25)
mu <- c(0, 0)
sigma <- c(2, 1)
set.seed(234567)
data <- rTrajWn2D(x0 = c(0, 0), alpha = alpha, mu = mu, sigma = sigma,
                    N = 200, delta = 0.5)
b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3],
                                                     sigma = sigma),
                                 mu = pars[4:5], sigma = sigma)
sigma2 <- function(x, pars) repRow(sigma^2, nrow(x))

exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma,
                  start = c(1, 1, 0, 1, 1),
                  lower = c(0.1, 0.1, -25, -25, -25),
                  upper = c(25, 25, 25, 25, 25), optMethod = "nlm")
pdeWn <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2,
                  Mx = 20, My = 20, Mt = 10, start = rbind(c(1, 1, 0, 1, 1)),
                  lower = c(0.1, 0.1, -25, -25, -25),
                  upper = c(25, 25, 25, 25, 25), verbose = 2,
                  optMethod = "nlm")
pdeWnLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2,
                     Mx = 20, My = 20, Mt = 10,
                     start = rbind(c(1, 1, 0, 1, 1)),
                     lower = c(0.1, 0.1, -25, -25, -25),
                     upper = c(25, 25, 25, 25, 25), verbose = 2,
                     linearBinning = TRUE)

head(exactOu)
head(pdeOu)
head(pdeOuLin)