psMle {sdetorus}

Maximum pseudo-likelihood estimation by wrapped pseudo-likelihoods

Description

Maximum pseudo-likelihood using the Euler and Shoji–Ozaki pseudo-likelihoods.

Usage

psMle(data, delta, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2,
  start, lower, upper, circular = TRUE, maxK = 2, vmApprox = FALSE, ...)

Arguments

data

a matrix of dimension c(n, p).

delta

discretization step.

method

a string for choosing "E" (Euler), "SO" (Shoji–Ozaki) or "SO2" (Shoji–Ozaki with Ito's expansion in the drift) method.

b

drift function. Must return a matrix of the same size as x.

jac.b

jacobian of the drift function.

sigma2

diagonal of the diffusion matrix (if univariate, this is the square of the diffusion coefficient). Must return an object of the same size as x.

b1

first derivative of the drift function (univariate). Must return a vector of the same length as x.

b2

second derivative of the drift function (univariate). Must return a vector of the same length as x.

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.

circular

flag to indicate circular data.

maxK

maximum absolute winding number used if circular = TRUE.

vmApprox

flag to indicate von Mises approximation to wrapped normal. See
momentMatchWnVm and scoreMatchWnBvm.

...

further parameters passed to mleOptimWrapper.

Details

See Section 3.2 in García-Portugués et al. (2019) for details. "SO2" implements Shoji and Ozai (1998)'s expansion with for p = 1. "SO" is the same expansion, for arbitrary p, but considering null second derivatives.

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

Shoji, I. and Ozaki, T. (1998) A statistical method of estimation and simulation for systems of stochastic differential equations. Biometrika, 85(1):240–243. doi:10.1093/biomet/85.1.240

Examples

# Example in 1D

delta <- 0.5
pars <- c(0.25, 0, 2)
set.seed(12345678)
samp <- rTrajWn1D(x0 = 0, alpha = pars[1], mu = pars[2], sigma = pars[3],
                  N = 100, delta = delta)
b <- function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2],
                                 sigma = pars[3], maxK = 2, expTrc = 30)
b1 <- function(x, pars, h = 1e-4) {
  l <- length(x)
  res <- b(x = c(x + h, x - h), pars = pars)
  drop(res[1:l] - res[(l + 1):(2 * l)])/(2 * h)
}
b2 <- function(x, pars, h = 1e-4) {
  l <- length(x)
  res <- b(x = c(x + h, x, x - h), pars = pars)
  drop(res[1:l] - 2 * res[(l + 1):(2 * l)] + res[(2 * l + 1):(3 * l)])/(h^2)
}
sigma2 <- function(x, pars) rep(pars[3]^2, length(x))
lower <- c(0.1, -pi, 0.1)
upper <- c(10, pi, 10)
psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2,
      start = pars, lower = lower, upper = upper)
psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2,
      start = pars, lower = lower, upper = upper, vmApprox = TRUE)
psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1,
      b2 = b2, sigma2 = sigma2, start = pars, lower = lower, upper = upper)
psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1,
      b2 = b2, sigma2 = sigma2, start = pars, lower = lower,
      upper = upper, vmApprox = TRUE)
psMle(data = samp, delta = delta, method = "SO", b = b, b1 = b1,
      lower = lower, upper = upper, sigma2 = sigma2, start = pars)
approxMleWn1D(data = samp, delta = delta, start = pars)
mlePde1D(data = samp, delta = delta, b = b, sigma2 = sigma2,
         start = pars, lower = lower, upper = upper)

# Example in 2D

delta <- 0.5
pars <- c(1, 0.5, 0, 0, 0, 1, 2)
set.seed(12345678)
samp <- rTrajWn2D(x0 = c(0, 0), alpha = pars[1:3], mu = pars[4:5],
                  sigma = pars[6:7], N = 100, delta = delta)
b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3],
                                                     sigma = pars[6:7]),
                                 mu = pars[4:5], sigma = pars[6:7], maxK = 2,
                                 expTrc = 30)
jac.b <- function(x, pars, h = 1e-4) {
  l <- nrow(x)
  res <- b(x = rbind(cbind(x[, 1] + h, x[, 2]),
                     cbind(x[, 1] - h, x[, 2]),
                     cbind(x[, 1], x[, 2] + h),
                     cbind(x[, 1], x[, 2] - h)), pars = pars)
  cbind(res[1:l, ] - res[(l + 1):(2 * l), ],
        res[2 * l + 1:l, ] - res[2 * l + (l + 1):(2 * l), ]) / (2 * h)
}
sigma2 <- function(x, pars) matrix(pars[6:7]^2, nrow = length(x) / 2L,
                                   ncol = 2)
lower <- c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01)
upper <- c(25, 25, 25, pi, pi, 25, 25)
psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2,
      start = pars, lower = lower, upper = upper)
psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2,
      start = pars, lower = lower, upper = upper, vmApprox = TRUE)
psMle(data = samp, delta = delta, method = "SO", b = b, jac.b = jac.b,
      sigma2 = sigma2, start = pars, lower = lower, upper = upper)
approxMleWn2D(data = samp, delta = delta, start = c(pars, 0))
# Set maxit = 5 to test and avoid a very long evaluation
mlePde2D(data = samp, delta = delta, b = b, sigma2 = sigma2, start = pars,
         lower = lower, upper = upper, maxit = 5)