[Yuima-commits] r448 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat May 28 00:02:39 CEST 2016
Author: lorenzo
Date: 2016-05-28 00:02:39 +0200 (Sat, 28 May 2016)
New Revision: 448
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/MM.COGARCH.R
pkg/yuima/R/PseudoLogLikCOGARCH.R
pkg/yuima/R/qmle.R
pkg/yuima/R/setMultiModel.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/qmle.Rd
pkg/yuima/man/rng.Rd
Log:
Fixed issue # uncoding bugs in Rd files
Improved Levy estimation in CP-COGARCH # Underlying CP-Levy is now estimated using the same qmle of an object of class yuima.poisson
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/DESCRIPTION 2016-05-27 22:02:39 UTC (rev 448)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.0.88
+Version: 1.0.89
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Modified: pkg/yuima/R/MM.COGARCH.R
===================================================================
--- pkg/yuima/R/MM.COGARCH.R 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/R/MM.COGARCH.R 2016-05-27 22:02:39 UTC (rev 448)
@@ -791,7 +791,8 @@
measure,
measure.type,
aggregation,
- dt){
+ dt,
+ noCPFFT = TRUE){
fixed.carma <- unlist(fixed)
lower.carma <- unlist(lower)
@@ -800,8 +801,23 @@
Dummy <- TRUE
+ if(noCPFFT && measure.type=="CP"){
+ L<-cumsum(c(0,Increment.lev))
+ mod1 <- setPoisson(intensity=as.character(measure$intensity),
+ df=list(as.character(measure$df$expr)))
+ Yui1<-setYuima(data=setData(L,delta=dt), model=mod1)
+ Noise <- qmle(yuima=Yui1,
+ start=param0,
+ upper = upper.carma,
+ lower=lower.carma,
+ fixed = fixed.carma)
+ result.Lev <- list(estLevpar = coef(Noise), covLev = Noise at vcov, value = Noise at min)
+
+ return(result.Lev)
+ }
+
CPlist <- c("dnorm","dgamma", "dexp")
- codelist <- c("rngamma","rNIG","rIG", "rgamma")
+ codelist <- c("rvgamma","rNIG","rIG", "rgamma")
if(measure.type=="CP"){
tmp <- regexpr("\\(", measure$df$exp)[1]
measurefunc <- substring(measure$df$exp, 1, tmp-1)
Modified: pkg/yuima/R/PseudoLogLikCOGARCH.R
===================================================================
--- pkg/yuima/R/PseudoLogLikCOGARCH.R 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/R/PseudoLogLikCOGARCH.R 2016-05-27 22:02:39 UTC (rev 448)
@@ -162,6 +162,7 @@
method = character(), model, objFun, observ=yuima at data,
fixed, meas.par, lower, upper, env, yuima, start, aggregation)
+
return(res)
}
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/R/qmle.R 2016-05-27 22:02:39 UTC (rev 448)
@@ -1232,7 +1232,7 @@
}
if(yuima at model@measure.type=="code"){
- # # "rIG", "rNIG", "rgamma", "rbgamma", "rngamma"
+ # # "rIG", "rNIG", "rgamma", "rbgamma", "rvgamma"
name.func.dummy <- as.character(model at measure$df$expr[1])
name.func<- substr(name.func.dummy,1,(nchar(name.func.dummy)-1))
names.measpar<-as.vector(strsplit(name.func,', '))[[1]][-1]
@@ -1306,7 +1306,7 @@
length(coef[ names.measpar]))
)
}
- if(measurefunc=="rngamma"){
+ if(measurefunc=="rvgamma"){
# result.Lev<-yuima.Estimation.VG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
# fixed.carma=fixed.carma,
# lower.carma=lower.carma,
@@ -2557,7 +2557,7 @@
}
-minusloglik.Lev<-function(par,env){
+minusloglik.Lev <- function(par,env){
if(env$measure.type=="code"){
if(env$measure=="rNIG"){
alpha<-par[1]
@@ -2569,12 +2569,12 @@
v1<-v[!is.infinite(v)]
-sum(v1)
}else{
- if(env$measure=="rngamma"){
+ if(env$measure=="rvgamma"){
lambda<-par[1]
alpha<-par[2]
beta<-par[3]
mu<-par[4]
- f<-dngamma(env$data,lambda,alpha,beta,mu)
+ f<-dvgamma(env$data,lambda,alpha,beta,mu)
v<-log(as.numeric(na.omit(f)))
v1<-v[!is.infinite(v)]
-sum(v1)
@@ -2644,13 +2644,13 @@
v1<-v[!is.infinite(v)]
return(sum(v1))
}else{
- if(env$measure=="rngamma"){
+ if(env$measure=="rvgamma"){
lambda<-params[1]
alpha<-params[2]
beta<-params[3]
mu<-params[4]
- #return(sum(log(dngamma(env$data,lambda,alpha,beta,mu))))
- f<-dngamma(env$data,lambda,alpha,beta,mu)
+ #return(sum(log(dvgamma(env$data,lambda,alpha,beta,mu))))
+ f<-dvgamma(env$data,lambda,alpha,beta,mu)
v<-log(as.numeric(na.omit(f)))
v1<-v[!is.infinite(v)]
return(sum(v1))
@@ -2706,7 +2706,7 @@
if(env$measure=="rNIG"){
Matr.dum<-diag(c(1,1,1/env$dt,1/env$dt))
}else{
- if(env$measure=="rngamma"){
+ if(env$measure=="rvgamma"){
Matr.dum<-diag(c(1/env$dt,1,1,1/env$dt))
}else{
if(env$measure=="rIG"){
@@ -2753,7 +2753,7 @@
param0[3]<-param0[3]*dt
param0[4]<-param0[4]*dt
}else{
- if(env$measure=="rngamma"){
+ if(env$measure=="rvgamma"){
#Matr.dum<-diag(c(1/env$dt,1,1,1/env$dt))
param0[1]<-param0[1]*dt
param0[4]<-param0[4]*dt
@@ -2782,7 +2782,7 @@
ui<-rbind(c(1, -1, 0, 0),c(1, 1, 0, 0),c(1, 0, 0, 0),c(0, 0, 1, 0))
ci<-c(0,0,0,10^(-6))
}else{
- if(measure=="rngamma"){
+ if(measure=="rvgamma"){
ui<-rbind(c(1,0, 0, 0),c(0, 1, 1, 0),c(0, 1,-1, 0),c(0, 1,0, 0))
ci<-c(10^-6,10^-6,10^(-6), 0)
}else{
@@ -2899,7 +2899,7 @@
paramLev[3]<-paramLev[3]/dt
paramLev[4]<-paramLev[4]/dt
}else{
- if(env$measure=="rngamma"){
+ if(env$measure=="rvgamma"){
#Matr.dum<-diag(c(1/env$dt,1,1,1/env$dt))
paramLev[1]<-paramLev[1]/dt
paramLev[4]<-paramLev[4]/dt
@@ -3056,7 +3056,7 @@
# alpha<-par[2]
# beta<-par[3]
# mu<-par[4]
-# -sum(log(dngamma(data,lambda,alpha,beta,mu)))
+# -sum(log(dvgamma(data,lambda,alpha,beta,mu)))
# }
#
# data<-Increment.lev
@@ -3134,7 +3134,7 @@
# beta<-params[3]
# mu<-params[4]
#
-# return(sum(log(dngamma(data,lambda,alpha,beta,mu))))
+# return(sum(log(dvgamma(data,lambda,alpha,beta,mu))))
# }
# # hessian <- tsHessian(param = params, fun = logLik.VG)
# #hessian<-optimHess(par, fn, gr = NULL,data=data)
Modified: pkg/yuima/R/setMultiModel.R
===================================================================
--- pkg/yuima/R/setMultiModel.R 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/R/setMultiModel.R 2016-05-27 22:02:39 UTC (rev 448)
@@ -67,7 +67,7 @@
##::make type function list ####
CPlist <- c("dnorm", "dgamma", "dexp", "dconst")
- codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rngamma", "rstable")
+ codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rvgamma", "rstable")
##::end make type function list ####
jump.dimens <- dim(jump.coeff)[2]
numbMeasure <- length(measure.type)
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/R/yuima.model.R 2016-05-27 22:02:39 UTC (rev 448)
@@ -159,7 +159,7 @@
##::make type function list ####
CPlist <- c("dnorm", "dgamma", "dexp", "dconst")
- codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rngamma", "rstable")
+ codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rvgamma", "rstable")
##::end make type function list ####
if(!length(measure.type)){
Modified: pkg/yuima/man/qmle.Rd
===================================================================
--- pkg/yuima/man/qmle.Rd 2016-05-27 03:59:26 UTC (rev 447)
+++ pkg/yuima/man/qmle.Rd 2016-05-27 22:02:39 UTC (rev 448)
@@ -1,445 +1,450 @@
-\name{qmle}
-%\alias{ql}
-%\alias{ml.ql}
-\alias{qmle}
-\alias{quasilogl}
-%\alias{LSE}
-%\alias{ml.ql2}
-\alias{rql}
-\alias{lse}
-\alias{pseudologlikelihood}
-\alias{pseudologlikelihood.COGARCH}
-%\alias{ql,ANY-method}
-%\alias{ml.ql,ANY-method}
-%\alias{ml.ql2,ANY-method}
-%\alias{rql,ANY-method}
-\title{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{
-%ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
-%ql(yuima,theta2,theta1,h,print=FALSE,param)
-%rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
-qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper,
- joint=FALSE, Est.Incr="Carma.IncPar",aggregation=TRUE, threshold=NULL, ...)
-quasilogl(yuima, param, print=FALSE)
-lse(yuima, start, lower, upper, method = "BFGS", ...)
-}
-\arguments{
- \item{yuima}{a yuima object.}
-% \item{theta2,theta1}{parameters of the sdeModel.}
-% \item{h}{ time span of observations.}
-% \item{theta2.lim, theta1.lim}{matrixes to specify the domains of the
-% parameters. Vector can be available only if theta is a scalar.}
-% \item{ptheta2,ptheta1}{}
- \item{print}{you can see a progress of the estimation when print is TRUE.}
- \item{method}{see Details.}
- \item{param}{\code{list} of parameters for the quasi loglikelihood.}
-% \item{interval}{}
-% \item{prevparam}{}
- \item{lower}{a named list for specifying lower bounds of parameters}
- \item{upper}{a named list for specifying upper bounds of parameters}
- \item{start}{initial values to be passed to the optimizer.}
- \item{fixed}{for conditional (quasi)maximum likelihood estimation.}
- \item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
- \item{Est.Incr}{If the yuima model is an object of \code{\link{yuima.carma-class}} the \code{qmle} returns an object of \code{\link{yuima.carma.qmle-class}} or object of class \code{mle-class}. By default \code{Est.Incr="Carma.IncPar"}.}
- \item{aggregation}{If \code{aggregation=TRUE}, before the estimation of the levy parameters we aggregate the increments.}
- \item{threshold}{If the model has Compund Poisson type jumps, the threshold is
- used to perform thresholding of the increments.}
- \item{...}{passed to \code{\link{optim}} method. See Examples.}
-}
-\details{
-% A function ql calculate the quasi-likelihood of a time series data X with any
-% parameters. A function ml.pl estimates parameters of the sdeModel by
-% maximizing the quasi-likelihood.
- \code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
- argument \code{method} is one of the methods available in \code{\link{optim}}.
-
- \code{lse} calculates least squares estimators of the drift parameters. This is
- useful for initial guess of \code{qmle} estimation.
-
- \code{quasilogl} returns the value of the quasi loglikelihood for a given
- \code{yuima} object and list of parameters \code{coef}.
-
-}
-\value{
- \item{QL}{a real value.}
- \item{opt}{a list with components the same as 'optim' function.}
- \item{carmaopt}{if the model is an object of \code{\link{yuima.carma-class}}, \code{qmle} returns an object \code{\link{yuima.carma.qmle-class}}}
- \item{cogarchopt}{if the model is an object of \code{\link{yuima.cogarch-class}}, \code{qmle} returns an object of class \code{\link{cogarch.gmm-class}}. The estimates are obtained by maximizing the pseudo-loglikelihood function as shown in Iacus et al. (2015)}
-}
-\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. \href{http://arxiv.org/abs/1511.00253}{http://arxiv.org/abs/1511.00253}
-}
-
-\author{The YUIMA Project Team}
-\note{
- %The function ml.ql uses the function optim internally.
- The function qmle uses the function optim internally.
-
- The function qmle uses the function \code{\link{CarmaNoise}} internally for estimation of underlying Levy if the model is an object of \code{\link{yuima.carma-class}}.
-}
-\examples{
-#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. theta2 = .3, theta1 = .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. theta2 = .3, theta1 = .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. theta2 = .3, theta1 = .1\n"))
-print(coef(opt4))
-
-## 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. theta2 = .3, theta1 = .1\n"))
-##print(coef(opt))
-
-
-\dontrun{
-
-###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 at 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 at 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 at 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 at 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)
-
-plot(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)
-
-plot(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)
-
-plot(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)
-
-}
-}
-
-
-% Add one or more standard keywords, see file 'KEYWORDS' in the
-% R documentation directory.
-\keyword{ts}
+\encoding{UTF-8}
+\name{qmle}
+%\alias{ql}
+%\alias{ml.ql}
+\alias{qmle}
+\alias{quasilogl}
+%\alias{LSE}
+%\alias{ml.ql2}
+\alias{rql}
+\alias{lse}
+\alias{pseudologlikelihood}
+\alias{pseudologlikelihood.COGARCH}
+%\alias{ql,ANY-method}
+%\alias{ml.ql,ANY-method}
+%\alias{ml.ql2,ANY-method}
+%\alias{rql,ANY-method}
+\title{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{
+%ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
+%ql(yuima,theta2,theta1,h,print=FALSE,param)
+%rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
+qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper,
+ joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL, ...)
+quasilogl(yuima, param, print=FALSE)
+lse(yuima, start, lower, upper, method = "BFGS", ...)
+}
+\arguments{
+ \item{yuima}{a yuima object.}
+% \item{theta2,theta1}{parameters of the sdeModel.}
+% \item{h}{ time span of observations.}
+% \item{theta2.lim, theta1.lim}{matrixes to specify the domains of the
+% parameters. Vector can be available only if theta is a scalar.}
+% \item{ptheta2,ptheta1}{}
+ \item{print}{you can see a progress of the estimation when print is TRUE.}
+ \item{method}{see Details.}
+ \item{param}{\code{list} of parameters for the quasi loglikelihood.}
+% \item{interval}{}
+% \item{prevparam}{}
+ \item{lower}{a named list for specifying lower bounds of parameters}
+ \item{upper}{a named list for specifying upper bounds of parameters}
+ \item{start}{initial values to be passed to the optimizer.}
+ \item{fixed}{for conditional (quasi)maximum likelihood estimation.}
+ \item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
+ \item{Est.Incr}{If the yuima model is an object of \code{\link{yuima.carma-class}} the \code{qmle} returns an object of \code{\link{yuima.carma.qmle-class}} or object of class \code{mle-class}. By default \code{Est.Incr="Carma.IncPar"}.}
+ \item{aggregation}{If \code{aggregation=TRUE}, before the estimation of the levy parameters we aggregate the increments.}
+ \item{threshold}{If the model has Compund Poisson type jumps, the threshold is
+ used to perform thresholding of the increments.}
+ \item{...}{passed to \code{\link{optim}} method. See Examples.}
+}
+\details{
+% A function ql calculate the quasi-likelihood of a time series data X with any
+% parameters. A function ml.pl estimates parameters of the sdeModel by
+% maximizing the quasi-likelihood.
+ \code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
+ argument \code{method} is one of the methods available in \code{\link{optim}}.
+
+ \code{lse} calculates least squares estimators of the drift parameters. This is
+ useful for initial guess of \code{qmle} estimation.
+
+ \code{quasilogl} returns the value of the quasi loglikelihood for a given
+ \code{yuima} object and list of parameters \code{coef}.
+
+}
+\value{
+ \item{QL}{a real value.}
+ \item{opt}{a list with components the same as 'optim' function.}
+ \item{carmaopt}{if the model is an object of \code{\link{yuima.carma-class}}, \code{qmle} returns an object \code{\link{yuima.carma.qmle-class}}}
+ \item{cogarchopt}{if the model is an object of \code{\link{yuima.cogarch-class}}, \code{qmle} returns an object of class \code{\link{cogarch.est-class}}. The estimates are obtained by maximizing the pseudo-loglikelihood function as shown in Iacus et al. (2015)}
+}
+\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 \enc{Probabilités}{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. \href{http://arxiv.org/abs/1511.00253}{http://arxiv.org/abs/1511.00253}
+}
+
+\author{The YUIMA Project Team}
+\note{
+ %The function ml.ql uses the function optim internally.
+ The function qmle uses the function optim internally.
+
+ The function qmle uses the function \code{\link{CarmaNoise}} internally for estimation of underlying Levy if the model is an object of \code{\link{yuima.carma-class}}.
+}
+\examples{
+#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. theta2 = .3, theta1 = .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. theta2 = .3, theta1 = .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. theta2 = .3, theta1 = .1\n"))
+print(coef(opt4))
+
+## 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. theta2 = .3, theta1 = .1\n"))
+##print(coef(opt))
+
+
+\dontrun{
+
+###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 at 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 at 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 at 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 at 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)
+
+plot(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)
+
+plot(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)
+
+plot(carmaopt3)
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 448
More information about the Yuima-commits
mailing list