[Yuima-commits] r359 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 6 19:25:59 CET 2015
Author: lorenzo
Date: 2015-02-06 19:25:59 +0100 (Fri, 06 Feb 2015)
New Revision: 359
Added:
pkg/yuima/R/MM.COGARCH.R
pkg/yuima/man/MM.COGARCH.rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/sim.euler.R
Log:
Added Method of Moments for COGARCH(P,Q)
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-11-28 10:48:18 UTC (rev 358)
+++ pkg/yuima/DESCRIPTION 2015-02-06 18:25:59 UTC (rev 359)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.54
-Date: 2014-11-27
+Version: 1.0.55
+Date: 2015-02-06
Depends: methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2014-11-28 10:48:18 UTC (rev 358)
+++ pkg/yuima/NAMESPACE 2015-02-06 18:25:59 UTC (rev 359)
@@ -123,12 +123,12 @@
export(qmleL)
export(CarmaNoise) # Estimates the Levy in carma model
+export(MM.COGARCH) # Estimation COGARCH(P,Q) using Method Of Moments
+
export(qgv)
export(mmfrac)
-
-
export(cbind.yuima)
S3method(print, phitest)
Added: pkg/yuima/R/MM.COGARCH.R
===================================================================
--- pkg/yuima/R/MM.COGARCH.R (rev 0)
+++ pkg/yuima/R/MM.COGARCH.R 2015-02-06 18:25:59 UTC (rev 359)
@@ -0,0 +1,322 @@
+# In this function we consider different kind estimation procedures for COGARCH(P,Q)
+# Model.
+is.COGARCH <- function(obj){
+ if(is(obj,"yuima"))
+ return(is(obj at model, "yuima.cogarch"))
+ if(is(obj,"yuima.cogarch"))
+ return(is(obj, "yuima.cogarch"))
+ return(FALSE)
+}
+
+# The estimation procedure for cogarch(p,q) implemented in this code are based on the
+# Chadraa phd's thesis
+MM.COGARCH<-function(yuima, data = NULL, start, method="BFGS", fixed = list(),
+ lower, upper){
+ print <- FALSE
+ call <- match.call()
+# First step we check if all inputs are inserted correctly.
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+ if( missing(start) )
+ yuima.stop("Starting values for the parameters are missing.")
+
+ if( !is.COGARCH(yuima) )
+ yuima.stop("The model is not an object of class yuima.")
+
+ if( !is(yuima,"yuima") && missing(data) )
+ yuima.stop("data are missing.")
+
+
+ if(is(yuima,"yuima")){
+ model<-yuima at model
+ if(is.null(data)){
+ observ<-yuima at data
+ }else{observ<-data}
+ }else{
+ if(is(yuima,"yuima.cogarch")){
+ model<-yuima
+ if(is.null(data)){
+ yuima.stop("Missing data")
+ }
+ observ<-data
+ }
+ }
+
+ if(!is(observ,"yuima.data")){
+ yuima.stop("Data must be an object of class yuima.data-class")
+ }
+
+ if( !missing(upper)){
+ yuima.warn("The upper constraints will be implemented as soon as possible")
+ upper <- list()
+ }
+
+ if( !missing(lower)){
+ yuima.warn("The lower constraints will be implemented as soon as possible")
+ lower <- list()
+ }
+
+ if( !missing(fixed)){
+ yuima.warn("The fixed constraints will be implemented as soon as possible")
+ fixed <- list()
+ }
+
+ # We identify the model parameters
+ info <- model at info
+
+ numb.ar <- info at q
+ ar.name <- paste(info at ar.par,c(numb.ar:1),sep="")
+
+ numb.ma <- info at p
+ ma.name <- paste(info at ma.par,c(1:numb.ma),sep="")
+
+ loc.par <- info at loc.par
+ xinit.name <- paste(info at Latent.var, c(1:numb.ar), sep = "")
+
+ meas.par <- model at parameter@measure
+
+ fixed.name <- names(fixed)
+ if(info at q==1){
+ # nm <- c(names(start), "EL1", "phi1","phi2")
+ nm <- c(names(start), "EL1")
+ }else{
+ # nm <- c(names(start), "EL1", "M2Lev","M4Lev")
+ nm <- c(names(start), "EL1")
+ }
+ # We identify the index of parameters
+ if(length(meas.par)!=0){
+ if(info at q==1){
+ # fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, measure.par, "EL1", "phi1","phi2")
+ fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, meas.par, "EL1")
+ }else{
+ # fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, measure.par, "EL1", "M2Lev","M4Lev")
+ fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, meas.par, "EL1")
+ }
+ }else{
+ if(info at q==1){
+ # fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, "EL1", "phi1","phi2")
+ fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, "EL1")
+ }else{
+ # fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, "EL1", "M2Lev","M4Lev")
+ fullcoeff <- c(ar.name, ma.name, loc.par, xinit.name, "EL1")
+ }
+ }
+ oo <- match(nm, fullcoeff)
+
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
+ if(info at q==1){
+ # start <- c(start, EL1 = 1, phi1=-1, phi2=-1)
+ start <- c(start, EL1 = 1)
+ }else{
+ # start <- c(start, EL1 = 1, M2Lev=1, M4Lev=2)
+ start <- c(start, EL1 = 1)
+ }
+ start <- start[order(oo)]
+
+ nm <- names(start)
+
+ ar.idx <- match(ar.name, fullcoeff)
+ ma.idx <- match(ma.name, fullcoeff)
+ loc.idx <- match(loc.par, fullcoeff)
+ meas.idx <- match(meas.par, fullcoeff)
+ fixed.idx <- match(fixed.name, fullcoeff)
+ EL1.idx <- match("EL1",fullcoeff) # We decide to pass EL1 as parameter !!!
+# if(info at q==1){
+# phi1.idx <- match("phi1",fullcoeff) # We decide to pass phi1 as parameter !!!
+# phi2.idx <- match("phi2",fullcoeff) # We decide to pass phi1 as parameter !!!
+# }else{
+# M2Lev.idx <- match("M2Lev",fullcoeff)
+# M4Lev.idx <- match("M4Lev",fullcoeff)
+# }
+ # We build the environment used as input in the
+ # objective function.
+
+ env <- new.env()
+ n <- length(observ)[1]
+
+ # Data
+ assign("Data", as.matrix(onezoo(observ)[,1]), envir=env)
+ assign("deltaData", frequency(onezoo(observ)[,1]), envir=env)
+ assign("time.obs",length(env$Data),envir=env)
+
+ # Order
+ assign("p", info at p, envir=env)
+ assign("q", info at q, envir=env)
+
+ # Idx
+ assign("ar.idx", ar.idx, envir=env)
+ assign("ma.idx", ma.idx, envir=env)
+ assign("loc.idx", loc.idx, envir=env)
+ assign("meas.idx", meas.idx, envir=env)
+ assign("EL1.idx", EL1.idx, envir=env)
+# if(info at q==1){
+# assign("phi1.idx", phi1.idx, envir=env)
+# assign("phi2.idx", phi1.idx, envir=env)
+# }else{
+# assign("M2Lev.idx", M2Lev.idx, envir=env)
+# assign("M4Lev.idx", M4Lev.idx, envir=env)
+# }
+
+objectiveFun <- function(p) {
+ mycoef <- as.list(p)
+ # names(mycoef) <- nm
+
+ if(length(c(fixed.idx, meas.idx))>0){ ## SMI 2/9/14
+ names(mycoef) <- nm[-c(fixed.idx,meas.idx)] ## SMI 2/9/14
+ }else{
+ names(mycoef) <- nm
+ }
+ ErrTerm(yuima=yuima, param=mycoef, print=print, env)
+}
+ out<- optim(start, objectiveFun, method=method)
+# , method = "L-BFGS-B",
+# lower=c(0.0001,0.00001,0.00001,0, 0.9,-100,-100),
+# upper=c(1,1,1,1,1.1,0,0))
+
+ return(out)
+
+}
+
+ErrTerm <- function(yuima, param, print, env){
+ # sample moments
+ typeacf <- "correlation"
+ param <- as.numeric(param)
+ G_i <- diff(env$Data)
+ r <- 1/env$deltaData
+ mu_G2 <- mean(G_i^2)
+ var_G2 <- mean(G_i^4) - mu_G2^2
+ d <- floor(sqrt(length(G_i)))
+ CovQuad <- log(abs(acf(G_i^2,plot=FALSE,lag.max=min(d,10),type=typeacf)$acf[-1]))
+ h <- seq(1, d, by = 1)*r
+ cost <- env$loc.idx
+ b <- env$ar.idx
+ a <- env$ma.idx
+ meanL1 <- param[env$EL1.idx]
+ #meanL1 <- 1
+# if(env$q == 1){
+#
+# beta <- param[cost]*param[b]
+# eta <- param[b]
+# phi <- param[a]
+#
+# # phi1 <- param[env$phi1.idx]
+# # phi2 <- param[env$phi2.idx]
+# #theo_mu_G2 <- meanL1*r*beta/abs(phi1)
+# phi1 <- meanL1*r*beta/mu_G2
+#
+# termA <- (6*mu_G2/r*beta/abs(phi1)*(2*eta/phi-meanL1)*(r-(1-exp(-r*abs(phi1)))/abs(phi1))+2*beta^2/phi^2*r)
+# phi2 <-2*termA*abs(phi1)/((var_G2-2*mu_G2^2)*abs(phi1)+termA)
+# if(typeacf == "covariance"){
+# TheoCovQuad <- meanL1*beta^2/abs(phi1)^3*(2*eta/phi-meanL1)*
+# (2/abs(phi2)-1/abs(phi1))*(1-exp(-r*abs(phi1)))*(exp(r*abs(phi1))-1)*
+# exp(-h*abs(phi1))
+# }else{
+# TheoCovQuad <- meanL1*beta^2/abs(phi1)^3*(2*eta/phi-meanL1)*
+# (2/abs(phi2)-1/abs(phi1))*(1-exp(-r*abs(phi1)))*(exp(r*abs(phi1))-1)*
+# exp(-h*abs(phi1))/(var_G2)
+# }
+#
+# }
+ if(env$q >= 1){
+# M2Lev <- param[env$M2Lev.idx]
+# M4Lev <- param[env$M4Lev.idx]
+# }else{
+# M2Lev <- param[env$phi1.idx]
+# M4Lev <- param[env$phi2.idx]
+# }
+ TheoCovQuad <- numeric(length = length(h))
+ for(i in c(1:length(h))){
+ MomentCog <- MM_Cogarch(p = env$p, q = env$q, cost=param[cost] , acoeff=param[a],
+ b=param[b], meanL1 = meanL1, r = r, h = h[i],
+ type = typeacf, m2=mu_G2, var=var_G2)
+
+
+ TheoCovQuad[i] <- MomentCog$acfG2
+ }
+ theo_mu_G2 <- MomentCog$meanG2
+ }
+ res <- sum((c(log(abs(TheoCovQuad)))-c(CovQuad))^2)
+ return(res)
+# if(env$q == 3){
+# MomentCog <- MM_Cogarch(p = env$p, q = env$q, a, b, r, h)
+# }
+# if(env$q >= 4){
+# MomentCog <- MM_Cogarch(p = env$p, q = env$q, a, b, r, h)
+# }
+}
+
+MM_Cogarch <- function(p, q, cost, acoeff, b, meanL1, r, h, type, m2, var){
+ # The code developed here is based on the acf for squared returns derived in the
+ # Chaadra phd Thesis
+ a <- e <- matrix(0,nrow=q,ncol=1)
+ e[q,1] <- 1
+ a[1:p,1] <- acoeff
+# mu <- M2Lev
+ bq <- b[1]
+ a1 <- a[1]
+ mu <-1/acoeff[1]*(bq-cost*bq*r*meanL1/m2)
+# rho <- M4Lev
+ B_tilde <- MatrixA(b[c(q:1)])+mu*e%*%t(a)
+ #beta <- as.numeric(B_tilde[q,])
+ meanG2 <- cost*bq*r/(bq-mu*a1)*meanL1
+ Inf_eps <- IdentM <- diag(q)
+ if(q==1){
+ Inf_eps[q,q] <- 1/(-2*B_tilde[1,1])
+ }
+ if(q==2){
+ Inf_eps[q,q] <- -B_tilde[2,1]
+ Inf_eps <- 1/(2*B_tilde[2,1]*B_tilde[2,2])*Inf_eps
+ }
+ if(q==3){
+ Inf_eps[1,1] <- B_tilde[3,3]/B_tilde[3,1]
+ Inf_eps[q,q] <- -B_tilde[3,2]
+ Inf_eps[1,3] <- -1
+ Inf_eps[3,1] <- -1
+ Inf_eps <- 1/(2*(B_tilde[3,3]*B_tilde[3,2]+B_tilde[3,1]))*Inf_eps
+ }
+ if(q>=4){
+ Inf_eps <- round(AsympVar(B_tilde,e,lower=0,upper=100,delta=0.01)*10^5)/10^5
+ }
+
+ term <- expm(B_tilde*h)
+ invB <- solve(B_tilde)
+ term1 <- invB%*%(IdentM-expm(-B_tilde*r))
+ term2 <- (expm(B_tilde*r)-IdentM)
+
+ P0_overRho <- 2*mu^2*(3*invB%*%(invB%*%term2-r*IdentM)-IdentM)%*%Inf_eps
+ Q0_overRho <- 6*mu*((r*IdentM-invB%*%term2)%*%Inf_eps
+ -invB%*%(invB%*%term2-r*IdentM)%*%Inf_eps%*%t(B_tilde))%*%e
+ m_overRho <- as.numeric(t(a)%*%Inf_eps%*%a)
+ Den<- (m_overRho*meanL1^2/m2^2*var*r^2+t(a)%*%Q0_overRho+t(a)%*%P0_overRho%*%a+1)
+ num <-(meanL1^2/m2^2*var-2*mu^2)*r^2
+ rh0 <- as.numeric(num/Den)
+
+
+ Inf_eps1 <- Inf_eps*rh0
+ Ph <- mu^2*term%*%term1%*%invB%*%term2%*%Inf_eps1
+ Qh <- mu*term%*%term1%*%(-term2%*%Inf_eps1-invB%*%term2%*%Inf_eps1%*%t(B_tilde))%*%e
+ m <- m_overRho*rh0
+
+
+
+ if(type=="correlation"){
+ VarTheo<-as.numeric(rh0*t(a)%*%P0_overRho%*%a+rh0*t(a)%*%Q0_overRho+2*r^2*mu^2+rh0)
+ acfG2 <- (t(a)%*%Ph%*%a+t(a)%*%Qh)/VarTheo
+ }else{
+ coeff <- cost^2*bq^2/((1-m)*(bq-mu*a1)^2)
+ acfG2 <- coeff*(t(a)%*%Ph%*%a+t(a)%*%Qh)
+ }
+ res <- list(acfG2=acfG2, meanG2=meanG2)
+ return(res)
+}
+AsympVar<-function(B_tilde,e,lower=0,upper=100,delta=0.1){
+ part<-seq(lower,upper-delta, by=delta)+delta/2
+ last <- length(part)
+ Integrand <- as.matrix(matrix(0,length(e),length(e)))
+ for(i in c(1:last)){
+ Integrand<-Integrand+expm(B_tilde*part[i])%*%e%*%t(e)%*%expm(t(B_tilde)*part[i])*delta
+ }
+ return(Integrand)
+}
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2014-11-28 10:48:18 UTC (rev 358)
+++ pkg/yuima/R/sim.euler.R 2015-02-06 18:25:59 UTC (rev 359)
@@ -187,7 +187,7 @@
# }
###
- if(sdeModel at measure.type == "CP"){ ##:: Compound-Poisson type
+ if(sdeModel at measure.type == "CP" ){ ##:: Compound-Poisson type
##:: delete 2010/09/13 for simulate func bug fix by s.h
## eta0 <- eval(sdeModel at measure$intensity)
@@ -256,17 +256,21 @@
if(JAMP==FALSE || sum(i==ij)==0){
Pi <- 0
}else{
- J <- eta0*randJ[j]/lambda
- j <- j+1
+ if(is.null(dL)){
+ J <- eta0*randJ[j]/lambda
+ j <- j+1
##cat(paste(J,"\n"))
##Pi <- zeta(dX, J)
- assign(sdeModel at jump.variable, J, env)
+ assign(sdeModel at jump.variable, J, env)
- if(sdeModel at J.flag){
- J <- 1
+ if(sdeModel at J.flag){
+ J <- 1
+ }
+
+ Pi <- p.b.j(t=i*delta,X=dX) * J
+ }else{# we add this part since we allow the user to specify the increment of CP LM 05/02/2015
+ Pi <- p.b.j(t=i*delta,X=dX) %*% dL[, i]
}
-
- Pi <- p.b.j(t=i*delta,X=dX) * J
##Pi <- p.b.j(t=i*delta, X=dX)
}
dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi
Added: pkg/yuima/man/MM.COGARCH.rd
===================================================================
--- pkg/yuima/man/MM.COGARCH.rd (rev 0)
+++ pkg/yuima/man/MM.COGARCH.rd 2015-02-06 18:25:59 UTC (rev 359)
@@ -0,0 +1,76 @@
+\name{MM.COGARCH}
+\alias{MM.COGARCH}
+\alias{Method of Moment COGARCH}
+\alias{Estimation COGARCH}
+\alias{est. COGARCH}
+\alias{est. yuima Cog}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Method of Moments for COGARCH(P,Q).
+}
+\description{
+The function returns the estimated parameters of a COGARCH(P,Q) model. The parameters are abtained by matching theoretical vs empirical autocorrelation function. The theoretical autocorrelation function is computed according the methodology developed in Chadraa (2009).
+}
+\usage{
+MM.COGARCH(yuima, data = NULL, start,
+ method="BFGS", fixed = list(), lower, upper)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuima}{a yuima object or an object of \code{\link{yuima.cogarch-class}}.}
+ \item{data}{an object of class \code{\link{yuima.data-class}} contains the observations available at uniformly spaced time. If \code{data=NULL}, the default, the function uses the data in an object of \code{\link{yuima-class}}.}
+ \item{start}{a \code{list} containing the starting values for the optimization routine.}
+ \item{method}{a string indicating one of the methods available in \code{\link{optim}}.}
+ \item{fixed}{a list of fixed parameters in optimization routine.}
+ \item{lower}{a named list for specifying lower bounds of parameters}
+ \item{upper}{a named list for specifying upper bounds of parameters}
+}
+%\details{
+%Please complete !!!
+%}
+\value{ The function returns a list with the same components of the object obtained when the function \code{\link{optim}} is used.
+}
+\references{
+Chadraa, E. (2009) Statistical Modeling with COGARCH(P,Q) Processes. Phd Thesis
+}
+\author{
+The YUIMA Project Team.
+}
+%\note{
+%% ~~further notes~~
+%}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+%\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+%}
+\examples{
+\dontrun{
+# Example COGARCH(1,1): the parameters are the same used in Haugh et al. 2005. In this case
+# we assume the underlying noise is a symmetric variance gamma.
+# As first step we define the COGARCH(1,1) in yuima:
+
+mod1 <- setCogarch(p = 1, q = 1, work = FALSE,
+ measure=list(df="rbgamma(z,1,sqrt(2),1,sqrt(2))"),
+ measure.type = "code", Cogarch.var = "y",
+ V.var = "v", Latent.var="x")
+
+param <- list(a1 = 0.038, b1 = 0.053,
+ a0 = 0.04/0.053, x1 = 0)
+
+# We generate a trajectory
+
+samp <- setSampling(Terminal=5000, n=5000)
+set.seed(100)
+sim1 <- simulate(mod1, sampling = samp, true.parameter = param)
+
+# We estimate the model
+
+res1 <- MM.COGARCH(yuima = sim1, start = param)
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ Method of Moments }
+\keyword{ Estimation COGARCH }% __ONLY ONE__ keyword per line
More information about the Yuima-commits
mailing list