qmle: Calculate quasi-likelihood and ML estimator of least squares estimator

Description

Calculate the quasi-likelihood and estimate of the parameters of the stochastic differential equation by the maximum likelihood method or least squares estimator of the drift parameter.

Usage




qmle(yuima, start, method = "L-BFGS-B", fixed = list(),
print = FALSE, envir = globalenv(), lower, upper, joint = FALSE, Est.Incr ="NoIncr",
aggregation = TRUE, threshold = NULL, rcpp =FALSE, ...)
quasilogl(yuima, param, print = FALSE, rcpp = FALSE)
lse(yuima, start, lower, upper, method = "BFGS", ...)

Value

QL: a real value.
opt: a list with components the same as 'optim' function.
carmaopt: if the model is an object of yuima.carma-class, qmle returns an object yuima.carma.qmle-class
cogarchopt: if the model is an object of yuima.cogarch-class, qmle returns an object of class cogarch.est-class. The estimates are obtained by maximizing the pseudo-loglikelihood function as shown in Iacus et al. (2015)

Arguments

yuima: a yuima object.

print: you can see a progress of the estimation when print is TRUE.
envir: an environment where the model coefficients are evaluated.
method: see Details.
param: list of parameters for the quasi loglikelihood.

lower: a named list for specifying lower bounds of parameters
upper: a named list for specifying upper bounds of parameters
start: initial values to be passed to the optimizer.
fixed: for conditional (quasi)maximum likelihood estimation.
joint: perform joint estimation or two stage estimation? by default joint=FALSE.
Est.Incr: If the yuima model is an object of yuima.carma-class or yuima.cogarch-class the qmle returns an object of yuima.carma.qmle-class, cogarch.est.incr-class,cogarch.est-class or object of class mle-class. By default Est.Incr="NoIncr", alternative values are IncrPar and Incr.
aggregation: If aggregation=TRUE, before the estimation of the levy parameters we aggregate the increments.
threshold: If the model has Compund Poisson type jumps, the threshold is used to perform thresholding of the increments.
...: passed to optim method. See Examples.
rcpp: use C++ code?

Author

The YUIMA Project Team

Details

qmle behaves more likely the standard mle function in stats4 and argument method is one of the methods available in optim.

lse calculates least squares estimators of the drift parameters. This is useful for initial guess of qmle estimation. quasilogl returns the value of the quasi loglikelihood for a given yuima object and list of parameters coef.

References

## Non-ergodic diffucion

Genon-Catalot, V., & Jacod, J. (1993). On the estimation of the diffusion coefficient for multi-dimensional diffusion processes. In Annales de l'IHP Probabilités et statistiques, 29(1), 119-151.

Uchida, M., & Yoshida, N. (2013). Quasi likelihood analysis of volatility and nondegeneracy of statistical random field. Stochastic Processes and their Applications, 123(7), 2851-2876.

## Ergodic diffusion

Kessler, M. (1997). Estimation of an ergodic diffusion from discrete observations. Scandinavian Journal of Statistics, 24(2), 211-229.

## Jump diffusion

Shimizu, Y., & Yoshida, N. (2006). Estimation of parameters for diffusion processes with jumps from discrete observations. Statistical Inference for Stochastic Processes, 9(3), 227-277.

Ogihara, T., & Yoshida, N. (2011). Quasi-likelihood analysis for the stochastic differential equation with jumps. Statistical Inference for Stochastic Processes, 14(3), 189-229.

## COGARCH

Iacus S. M., Mercuri L. and Rroji E.(2015) Discrete time approximation of a COGARCH (p, q) model and its estimation. tools:::Rd_expr_doi("10.48550/arXiv.1511.00253")

## CARMA

Iacus S. M., Mercuri L. (2015) Implementation of Levy CARMA model in Yuima package. Comp. Stat. (30) 1111-1141. tools:::Rd_expr_doi("10.1007/s00180-015-0569-7")

Examples

Run this code

#dXt^e = -theta2 * Xt^e * dt + theta1 * dWt
diff.matrix <- matrix(c("theta1"), 1, 1)
ymodel <- setModel(drift=c("(-1)*theta2*x"), diffusion=diff.matrix,
  time.variable="t", state.variable="x", solve.variable="x")
n <- 100

ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
yuima <- setYuima(model=ymodel, sampling=ysamp)
set.seed(123)
yuima <- simulate(yuima, xinit=1, true.parameter=list(theta1=0.3,
theta2=0.1))
QL <- quasilogl(yuima, param=list(theta2=0.8, theta1=0.7))
##QL <- ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)))
QL

## another way of parameter specification
##param <- list(theta2=0.8, theta1=0.7)
##QL <- ql(yuima, h=1/((n)^(2/3)), param=param)
##QL


## old code
##system.time(
##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
##)
##cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
##print(coef(opt))


system.time(
opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
 upper=list(theta1=1,theta2=1), method="L-BFGS-B")
)
cat(sprintf("\nTrue param. theta1 = .3, theta2 = .1\n"))
print(coef(opt2))

## initial guess for theta2 by least squares estimator
tmp <- lse(yuima, start=list(theta2=0.7), lower=list(theta2=0), upper=list(theta2=1))
tmp

system.time(
opt3 <- qmle(yuima, start=list(theta1=0.8, theta2=tmp), lower=list(theta1=0,theta2=0),
 upper=list(theta1=1,theta2=1), method="L-BFGS-B")
)
cat(sprintf("\nTrue param. theta1 = .3, theta2 = .1\n"))
print(coef(opt3))


## perform joint estimation? Non-optimal, just for didactic purposes
system.time(
opt4 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
 upper=list(theta1=1,theta2=1), method="L-BFGS-B", joint=TRUE)
)
cat(sprintf("\nTrue param. theta1 = .3, theta2 = .1\n"))
print(coef(opt4))

## fix theta1 to the true value
system.time(
opt5 <- qmle(yuima, start=list(theta2=0.7), lower=list(theta2=0),
upper=list(theta2=1),fixed=list(theta1=0.3), method="L-BFGS-B")
)
cat(sprintf("\nTrue param. theta1 = .3, theta2 = .1\n"))
print(coef(opt5))

## old code
##system.time(
##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newton")
##)
##cat(sprintf("\nTrue param. theta1 = .3, theta2 = .1\n"))
##print(coef(opt))


if (FALSE) {
###multidimension case
##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
diff.matrix <- matrix(c("theta1.1","theta1.2", "1", "1"), 2, 2)

drift.c <- c("-theta2.1*x1", "-theta2.2*x2", "-theta2.2", "-theta2.1")
drift.matrix <- matrix(drift.c, 2, 2)

ymodel <- setModel(drift=drift.matrix, diffusion=diff.matrix, time.variable="t",
                   state.variable=c("x1", "x2"), solve.variable=c("x1", "x2"))
n <- 100
ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
yuima <- setYuima(model=ymodel, sampling=ysamp)
set.seed(123)

##xinit=c(x1, x2) #true.parameter=c(theta2.1, theta2.2, theta1.1, theta1.2)
yuima <- simulate(yuima, xinit=c(1, 1),
true.parameter=list(theta2.1=0.5, theta2.2=0.3, theta1.1=0.6, theta1.2=0.2))

## theta2 <- c(0.8, 0.2) #c(theta2.1, theta2.2)
##theta1 <- c(0.7, 0.1) #c(theta1.1, theta1.2)
## QL <- ql(yuima, theta2, theta1, h=1/((n)^(2/3)))
## QL

## ## another way of parameter specification
## #param <- list(theta2=theta2, theta1=theta1)
## #QL <- ql(yuima, h=1/((n)^(2/3)), param=param)
## #QL

## theta2.1.lim <- c(0, 1)
## theta2.2.lim <- c(0, 1)
## theta1.1.lim <- c(0, 1)
## theta1.2.lim <- c(0, 1)
## theta2.lim <- t( matrix( c(theta2.1.lim, theta2.2.lim), 2, 2) )
## theta1.lim <- t( matrix( c(theta1.1.lim, theta1.2.lim), 2, 2) )

## system.time(
## opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
## )
## opt@coef

system.time(
opt2 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1),
 lower=list(theta1.1=.1,theta1.2=.1,theta2.1=.1,theta2.2=.1),
 upper=list(theta1.1=4,theta1.2=4,theta2.1=4,theta2.2=4), method="L-BFGS-B")
)
opt2@coef
summary(opt2)

## unconstrained optimization
system.time(
opt3 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
)
opt3@coef
summary(opt3)

quasilogl(yuima, param=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))

##system.time(
##opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, method="Newton")
##)
##opt@coef
##

# carma(p=2,q=0) driven by a brownian motion without location parameter

mod0<-setCarma(p=2,
               q=0,
               scale.par="sigma")

true.parm0 <-list(a1=1.39631,
                 a2=0.05029,
                 b0=1,
                 sigma=0.23)

samp0<-setSampling(Terminal=100,n=250)
set.seed(123)
sim0<-simulate(mod0,
               true.parameter=true.parm0,
               sampling=samp0)

system.time(
carmaopt0 <- qmle(sim0, start=list(a1=1.39631,a2=0.05029,
                              b0=1,
                               sigma=0.23))
)

summary(carmaopt0)

# carma(p=2,q=1) driven by a brownian motion without location parameter

mod1<-setCarma(p=2,
               q=1)

true.parm1 <-list(a1=1.39631,
                  a2=0.05029,
                  b0=1,
                  b1=2)

samp1<-setSampling(Terminal=100,n=250)
set.seed(123)
sim1<-simulate(mod1,
               true.parameter=true.parm1,
               sampling=samp1)

system.time(
  carmaopt1 <- qmle(sim1, start=list(a1=1.39631,a2=0.05029,
                                     b0=1,b1=2),joint=TRUE)
)

summary(carmaopt1)

# carma(p=2,q=1) driven by a compound poisson process where the jump size is normally distributed.

mod2<-setCarma(p=2,
               q=1,
               measure=list(intensity="1",df=list("dnorm(z, 0, 1)")),
               measure.type="CP")

true.parm2 <-list(a1=1.39631,
                  a2=0.05029,
                  b0=1,
                  b1=2)

samp2<-setSampling(Terminal=100,n=250)
set.seed(123)
sim2<-simulate(mod2,
               true.parameter=true.parm2,
               sampling=samp2)

system.time(
  carmaopt2 <- qmle(sim2, start=list(a1=1.39631,a2=0.05029,
                                     b0=1,b1=2),joint=TRUE)
)

summary(carmaopt2)

# carma(p=2,q=1) driven by a normal inverse gaussian process
mod3<-setCarma(p=2,q=1,
               measure=list(df=list("rNIG(z, alpha, beta, delta1, mu)")),
               measure.type="code")
#

# True param
true.param3<-list(a1=1.39631,
                 a2=0.05029,
                 b0=1,
                 b1=2,
                 alpha=1,
                 beta=0,
                 delta1=1,
                 mu=0)

samp3<-setSampling(Terminal=100,n=200)
set.seed(123)

sim3<-simulate(mod3,
               true.parameter=true.param3,
               sampling=samp3)


carmaopt3<-qmle(sim3,start=true.param3)

summary(carmaopt3)

# Simulation and Estimation of COGARCH(1,1) with CP driven noise

# Model parameters
eta<-0.053
b1 <- eta
beta <- 0.04
a0 <- beta/b1
phi<- 0.038
a1 <- phi


# Definition

cog11<-setCogarch(p = 1,q = 1,
  measure = list(intensity = "1",
                 df = list("dnorm(z, 0, 1)")),
  measure.type = "CP",
  XinExpr=TRUE)

# Parameter
paramCP11 <- list(a1 = a1, b1 =  b1,
                  a0 = a0, y01 = 50.31)
# Sampling scheme
samp11 <- setSampling(0, 3200, n=64000)

# Simulation
set.seed(125)

SimTime11 <- system.time(
  sim11 <- simulate(object = cog11,
    true.parameter = paramCP11,
    sampling = samp11,
    method="mixed")
)

plot(sim11)

# Estimation

timeComp11<-system.time(
  res11 <- qmle(yuima = sim11,
    start = paramCP11,
    grideq = TRUE,
    method = "Nelder-Mead")
)

timeComp11

unlist(paramCP11)

coef(res11)

# COGARCH(2,2) model driven by CP

cog22 <- setCogarch(p = 2,q = 2,
  measure = list(intensity = "1",
                 df = list("dnorm(z, 0, 1)")),
  measure.type = "CP",
  XinExpr=TRUE)

# Parameter

paramCP22 <- list(a1 = 0.04, a2 = 0.001,
  b1 =  0.705, b2 = 0.1, a0 = 0.1, y01 = (1 + 2 / 3),
  y02=0)

# Use diagnostic.cog for checking the stat and positivity

check22 <- Diagnostic.Cogarch(cog22, param = paramCP22)

# Sampling scheme

samp22 <- setSampling(0, 3600, n = 64000)

# Simulation

set.seed(125)
SimTime22 <- system.time(
  sim22 <- simulate(object = cog22,
    true.parameter = paramCP22,
    sampling = samp22,
    method = "Mixed")
)

plot(sim22)

timeComp22 <- system.time(
  res22 <- qmle(yuima = sim22,
    start = paramCP22,
    grideq=TRUE,
    method = "Nelder-Mead")
)

timeComp22

unlist(paramCP22)

coef(res22)

}

Run the code above in your browser using DataLab