psMle: Maximum pseudo-likelihood estimation by wrapped...

In sdetorus: Statistical Tools for Toroidal Diffusions

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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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)

sdetorus documentation built on Aug. 21, 2023, 1:08 a.m.