[Yuima-commits] r376 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 25 18:31:22 CEST 2015
Author: lorenzo
Date: 2015-04-25 18:31:21 +0200 (Sat, 25 Apr 2015)
New Revision: 376
Added:
pkg/yuima/man/Diagnostic.Cogarch.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/DiagnosticCogarch.R
pkg/yuima/R/MM.COGARCH.R
pkg/yuima/R/cogarchNoise.R
pkg/yuima/R/simulate.R
pkg/yuima/man/cogarch.gmm.incr.rd
pkg/yuima/man/gmm.rd
Log:
Add routines for COGARCH(p,q) model
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/DESCRIPTION 2015-04-25 16:31:21 UTC (rev 376)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.64
-Date: 2015-04-21
+Version: 1.0.65
+Date: 2015-04-25
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 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/NAMESPACE 2015-04-25 16:31:21 UTC (rev 376)
@@ -1,146 +1,149 @@
-import("methods")
-
-##importFrom("stats", "end", "time", "start")
-importFrom("graphics", "plot")
-import("zoo")
-importFrom(stats, confint)
-import("stats4")
-import("expm")
-import("mvtnorm")
-import("cubature")
-
-importFrom(utils, toLatex)
-
-
-
-
-
-
-exportClasses("yuima",
- "yuima.data",
- "yuima.sampling",
- "yuima.characteristic",
- "yuima.model",
- "model.parameter",
- "yuima.carma",
- "carma.info",
- "yuima.carma.qmle",
- "yuima.poisson",
- "yuima.qmle",
- "yuima.CP.qmle",
- "cogarch.info",
- "yuima.cogarch"
- )
-
-exportMethods(
- "dim",
- "length",
- ## "start",
- "plot",
- ## "time",
- ## "end",
- "simulate",
- "cce",
- "llag",
- "poisson.random.sampling",
- "get.zoo.data",
- "initialize",
-## "ql",
-## "rql",
-## "ml.ql",
- "adaBayes",
- "limiting.gamma",
- "getF",
- "getf",
- "getxinit",
- "gete",
- "simFunctional",
- "F0",
- "Fnorm",
- "asymptotic_term",
- "cbind.yuima"
- )
-
-## function which we want to expose to the user
-## all other functions are used internally by the
-## package
-export(setYuima)
-export(setModel) ## builds sde model
-export(setData)
-export(setSampling)
-export(setCharacteristic)
-export(setCarma)
-export(setPoisson)
-export(dconst)
-export(rconst)
-
-export(setCogarch)
-
-export(dim)
-export(length)
-#export(start)
-export(plot)
-#export(time)
-#export(end)
-
-export(simulate) # simulates couple of processes
-export(subsampling)
-export(cce)
-export(llag)
-export(poisson.random.sampling)
-export(noisy.sampling)
-export(mpv)
-export(bns.test)
-export(hyavar) # asymptotic variance estimator for the Hayashi-Yoshida estimator
-
-export(get.zoo.data)
-
-##export(ql,rql,ml.ql)
-##export(rql)
-export(adaBayes)
-export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
-export(limiting.gamma)
-
-export(setFunctional)
-export(getF)
-export(getf)
-export(getxinit)
-export(gete)
-
-export(simFunctional)
-export(F0)
-export(Fnorm)
-export(asymptotic_term)
-
-##export(LSE)
-export(lse)
-
-export(qmle)
-export(quasilogl)
-export(phi.test)
-export(lasso)
-export(CPoint)
-export(qmleR)
-export(qmleL)
-
-export(CarmaNoise) # Estimates the Levy in carma model
-export(gmm) # Estimation COGARCH(P,Q) using Method Of Moments
-
-
-export(qgv)
-export(mmfrac)
-
-export(cbind.yuima)
-
-S3method(print, phitest)
-S3method(print, qgv)
-S3method(print, mmfrac)
-S3method(print, yuima.lasso)
-
-S3method(toLatex, yuima)
-S3method(toLatex, yuima.model)
-S3method(toLatex, yuima.carma)
-S3method(toLatex, yuima.cogarch)
-
-useDynLib(yuima)
-
+import("methods")
+
+##importFrom("stats", "end", "time", "start")
+importFrom("graphics", "plot")
+import("zoo")
+importFrom(stats, confint)
+import("stats4")
+import("expm")
+import("mvtnorm")
+import("cubature")
+
+importFrom(utils, toLatex)
+
+
+
+
+
+
+exportClasses("yuima",
+ "yuima.data",
+ "yuima.sampling",
+ "yuima.characteristic",
+ "yuima.model",
+ "model.parameter",
+ "yuima.carma",
+ "carma.info",
+ "yuima.carma.qmle",
+ "yuima.poisson",
+ "yuima.qmle",
+ "yuima.CP.qmle",
+ "cogarch.info",
+ "yuima.cogarch"
+ )
+
+exportMethods(
+ "dim",
+ "length",
+ ## "start",
+ "plot",
+ ## "time",
+ ## "end",
+ "simulate",
+ "cce",
+ "llag",
+ "poisson.random.sampling",
+ "get.zoo.data",
+ "initialize",
+## "ql",
+## "rql",
+## "ml.ql",
+ "adaBayes",
+ "limiting.gamma",
+ "getF",
+ "getf",
+ "getxinit",
+ "gete",
+ "simFunctional",
+ "F0",
+ "Fnorm",
+ "asymptotic_term",
+ "cbind.yuima"
+ )
+
+## function which we want to expose to the user
+## all other functions are used internally by the
+## package
+export(setYuima)
+export(setModel) ## builds sde model
+export(setData)
+export(setSampling)
+export(setCharacteristic)
+export(setCarma)
+export(setPoisson)
+export(dconst)
+export(rconst)
+
+export(setCogarch)
+
+export(dim)
+export(length)
+#export(start)
+export(plot)
+#export(time)
+#export(end)
+
+export(simulate) # simulates couple of processes
+export(subsampling)
+export(cce)
+export(llag)
+export(poisson.random.sampling)
+export(noisy.sampling)
+export(mpv)
+export(bns.test)
+export(hyavar) # asymptotic variance estimator for the Hayashi-Yoshida estimator
+
+export(get.zoo.data)
+
+##export(ql,rql,ml.ql)
+##export(rql)
+export(adaBayes)
+export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
+export(limiting.gamma)
+
+export(setFunctional)
+export(getF)
+export(getf)
+export(getxinit)
+export(gete)
+
+export(simFunctional)
+export(F0)
+export(Fnorm)
+export(asymptotic_term)
+
+##export(LSE)
+export(lse)
+
+export(qmle)
+export(quasilogl)
+export(phi.test)
+export(lasso)
+export(CPoint)
+export(qmleR)
+export(qmleL)
+
+
+export(CarmaNoise) # Estimates the Levy in carma model
+export(gmm) # Estimation COGARCH(P,Q) using Method Of Moments
+export(cogarchNoise)
+export(Diagnostic.Cogarch)
+
+
+export(qgv)
+export(mmfrac)
+
+export(cbind.yuima)
+
+S3method(print, phitest)
+S3method(print, qgv)
+S3method(print, mmfrac)
+S3method(print, yuima.lasso)
+
+S3method(toLatex, yuima)
+S3method(toLatex, yuima.model)
+S3method(toLatex, yuima.carma)
+S3method(toLatex, yuima.cogarch)
+
+useDynLib(yuima)
+
Modified: pkg/yuima/R/DiagnosticCogarch.R
===================================================================
--- pkg/yuima/R/DiagnosticCogarch.R 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/R/DiagnosticCogarch.R 2015-04-25 16:31:21 UTC (rev 376)
@@ -1,14 +1,14 @@
-StartValCog<-function(yuima.cogarch, param=list(), mu=1, rho=3){
+# We write a Diagnostic function that evaluates the following quantity
+Diagnostic.Cogarch <- function(yuima.cogarch, param = list(),
+ matrixS = NULL , mu = 1,
+ display =TRUE){
if(missing(yuima.cogarch))
yuima.stop("yuima.cogarch or yuima object is missing.")
- if(length(param)==0)
- yuima.stop("missing values parameters")
-
- if(!is.COGARCH(yuima.cogarch)){
- yuima.warn("The model does not belong to the class yuima.cogarch")
+ if((length(param)==0 && !is(yuima.cogarch, "cogarch.gmm"))){
+ yuima.stop("missing values parameters")
}
@@ -17,9 +17,20 @@
}else{
if(is(yuima.cogarch,"yuima.cogarch")){
model<-yuima.cogarch
+ }else{
+ if(is(yuima.cogarch,"cogarch.gmm")){
+ model<-yuima.cogarch at model
+ if(length(param)==0){
+ param<-coef(yuima.cogarch)
+ }
+ }
}
}
+ if(!is.COGARCH(model)){
+ yuima.warn("The model does not belong to the class yuima.cogarch")
+ }
+
info <- model at info
numb.ar <- info at q
ar.name <- paste(info at ar.par,c(numb.ar:1),sep="")
@@ -27,30 +38,117 @@
ma.name <- paste(info at ma.par,c(1:numb.ma),sep="")
loc.par <- info at loc.par
- nm <- c(names(param))
- param<-as.numeric(param)
- names(param)<-nm
+ if(is.list(param)){
+ param <- unlist(param)
+ }
-
xinit.name0 <- model at parameter@xinit
idx <- na.omit(match(c(loc.par, ma.name), xinit.name0))
xinit.name <- xinit.name0[-idx]
+
+
fullcoeff <- c(ar.name, ma.name, loc.par,xinit.name)
-
+ nm<-names(param)
oo <- match(nm, fullcoeff)
-
- if(any(is.na(oo)))
- yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima.cogarch model")
-
+ if(!is(yuima.cogarch,"cogarch.gmm")){
+ if(length(na.omit(oo))!=length(fullcoeff))
+ yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima.cogarch model")
+ }
acoeff <- param[ma.name]
b <- param[ar.name]
cost<- param[loc.par]
- res<-StationaryMoments(cost,b,acoeff,mu,rho)
+ # Check Strictly Stationary condition
+ Amatr<-MatrixA(b[c(info at q:1)])
+ if(is.null(matrixS)){
+ Dum <- yuima.carma.eigen(Amatr)
+ matrixS <- Dum$vectors
+ lambda.eig <- Dum$values
+ }else{
+ if(!is.matrix(matrixS)){
+ if(dim(matrixS)[1]!=dim(matrixS)[2]){
+ yuima.stop("matrixS must be an object of class matrix")
+ }else{
+ if(dim(matrixS)[1]!=numb.ar)
+ yuima.stop("matrixS must be an square matrix. Its row number is equal to the dimension of autoregressive coefficients")
+ }
+ }
+ lambda.eig<-diag(solve(matrixS)%*%Amatr%*%matrixS)
+ }
+ lambda1<- max(Re(lambda.eig)) # we find lambda1
+ ev.dum<-matrix(0,info at q,1)
+ ev.dum[info at q,1] <- 1
+ av.dum <-matrix(0,1,info at q)
+ av.dum[1,c(1:info at p)]<-acoeff
+ if(display==TRUE){
+ cat(paste0("\n COGARCH(",info at p,info at q,") model \n"))
+ }
+ matrforck<-solve(matrixS)%*%(av.dum%*%ev.dum)%*%matrixS
+ if(matrforck*mu < -lambda1){
+ if(display==TRUE){
+ cat("\n The process is strictly stationary\n The unconditional first moment of the Variance process exists \n")
+ }
+ res.stationarity <- TRUE
+ }else{
+ if(display==TRUE){
+ cat("\n We are not able to establish if the process is stationary \n
+ Try a new S matrix such that the Companion matrix is diagonalizable")
+ }
+ res.stationarity <- "Try a different S matrix"
+ }
+ # Check the positivity
+ massage <- "\n We are not able to establish the positivity of the Variance \n"
+ res.pos <- " "
+ if(info at p==1){
+ if(is.numeric(lambda.eig) && all(lambda.eig<0)){
+ massage <- "\n the Variance is a positive process \n"
+ res.pos <- TRUE
+ }
+ if(is.complex(lambda.eig) && all(Im(lambda.eig)==0) && all(Re(lambda.eig) < 0)){
+ massage <- "\n the Variance is a positive process \n"
+ res.pos <- TRUE
+ }
+ }
+ if((info at q==2 && info at p==2 )){
+ if(is.numeric(lambda.eig)|| (is.complex(lambda.eig) && all(Im(lambda.eig)==0))){
+ if(acoeff[1]>=-acoeff[2]*lambda1){
+ massage <- "\n the Variance is a positive process \n"
+ res.pos <- TRUE
+ }else{
+ massage <- "\n the Variance process is not strictly positive \n"
+ res.pos <- FALSE
+ }
+ }
+ }
+ if(info at p>=2 && info at q>info at p){
+ gamma.eig<-polyroot(acoeff)
+ if(all(Im(gamma.eig)==0)){
+ gamm.eig<-as.numeric(gamma.eig)
+ }
+ if(is.numeric(lambda.eig) && all(is.numeric(gamma.eig)<0) && all(is.numeric(lambda.eig)<0)){
+ if(cumsum(sort(lambda.eig[c(1:(info at p-1))]))>=cumsum(sort(gamma.eig[c(1:(info at p-1))]))){
+ massage <- "\n The Variance process is strictly positive. \n"
+ }
+ }
+ }
+ if(display==TRUE){
+ cat(massage)
+ }
+
+ res1<-StationaryMoments(cost,b,acoeff,mu)
+ res <- list(meanVarianceProc=res1$ExpVar,
+ meanStateVariable=res1$ExpStatVar,
+ stationary = res.stationarity,
+ positivity = res.pos)
return(res)
}
-StationaryMoments<-function(cost,b,acoeff,mu=1,rho=3){
+yuima.norm2<-function(x){
+ res<-sqrt(max(eigen(t(x) %*% x)$values))
+ return(res)
+}
+
+StationaryMoments<-function(cost,b,acoeff,mu=1){
# We obtain stationary mean of State process
# stationary mean of the variance
# E\left(Y\right)=-a_{0}m_{2}\left(A+m_{2}ea'\right)^{-1}e
Modified: pkg/yuima/R/MM.COGARCH.R
===================================================================
--- pkg/yuima/R/MM.COGARCH.R 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/R/MM.COGARCH.R 2015-04-25 16:31:21 UTC (rev 376)
@@ -44,10 +44,11 @@
# The estimation procedure for cogarch(p,q) implemented in this code are based on the
# Chadraa phd's thesis
gmm<-function(yuima, data = NULL, start, method="BFGS", fixed = list(),
- lower, upper, lag.max = NULL, aggr.G = TRUE, aggregation=TRUE,
+ lower, upper, lag.max = NULL, equally.spaced = TRUE, aggregation=TRUE,
Est.Incr = "NoIncr", objFun = "L2"){
print <- FALSE
+ aggr.G <- equally.spaced
call <- match.call()
if(objFun=="L1" && method!="Nelder-Mead"){
@@ -214,14 +215,24 @@
assign("objFun",objFun, envir=env)
if(aggr.G==TRUE){
+ if(floor(n/index(observ at zoo.data[[1]])[n])!=env$deltaData){
+ yuima.stop("the n/Terminal in sampling information is not an integer. equally.spaced=FALSE is recommended")
+ }
+ }
+
+ if(aggr.G==TRUE){
# Aggregate returns G
#dt<-round(deltat(onezoo(observ)[,1])*10^5)/10^5
# Time<-index(observ at zoo.data[[1]])[n]
G_i <- diff(env$Data[seq(1,length(env$Data),by=env$deltaData)])
r<-1
}else{
- G_i <- diff(env$Data)
- r <- 1/env$deltaData
+ dummydata<-index(onezoo(observ)[,1])
+ unitarytime<-floor(dummydata)
+ index<-!duplicated(unitarytime)
+ G_i <- diff(env$Data[index])
+
+ r <- 1
}
d <- min(floor(sqrt(length(G_i))),env$lag)
assign("d", d, envir=env)
@@ -601,7 +612,7 @@
# res <- log(sum((TheoCovQuad[CovQuad>0]-CovQuad[CovQuad>0])^2))
# emp <- log(CovQuad[CovQuad>0])
# theo <- log(TheoCovQuad[CovQuad>0])
- res <- sum((log(TheoCovQuad[CovQuad>0])-log(CovQuad[CovQuad>0]))^2)
+ res <- sum((log(abs(TheoCovQuad[CovQuad>0]))-log(abs(CovQuad[CovQuad>0])))^2)
# res <- sum((TheoCovQuad[CovQuad>0]-CovQuad[CovQuad>0])^2)
# res <- sum((log(abs(TheoCovQuad))-log(abs(CovQuad)))^2)
# res <- sum((log(TheoCovQuad[CovQuad>0]))-log(CovQuad[CovQuad>0]))^2)
@@ -648,6 +659,7 @@
mu<-1 # we assume variance of the underlying L\'evy is one
meanL1<-mu
+
cost<-(bq-mu*a1)*m2/(bq*r)
B<- MatrixA(b[c(q:1)])
@@ -765,10 +777,14 @@
yuima.warn("Complex increments. We plot only the real part")
Incr.L<-Re(Incr.L)
}
- plot(x=Time,y=Incr.L, type=type,...)
+# model <- x at model
+# EndT <- Time[length(Time)]
+# numb <- (length(Incr.L)+1)
+ plot(x= Time,y= Incr.L, type=type,...)
}
)
+
setMethod("summary", "cogarch.gmm",
function (object, ...)
{
@@ -836,7 +852,7 @@
function (object)
{
- cat("Two Stage GMM estimation \n\nCall:\n")
+ cat("Two Stages GMM estimation \n\nCall:\n")
print(object at call)
cat("\nCoefficients:\n")
print(coef(object))
Modified: pkg/yuima/R/cogarchNoise.R
===================================================================
--- pkg/yuima/R/cogarchNoise.R 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/R/cogarchNoise.R 2015-04-25 16:31:21 UTC (rev 376)
@@ -49,11 +49,11 @@
- oo <- match(nm, fullcoeff)
+# oo <- match(nm, fullcoeff)
+#
+# if(any(is.na(oo)))
+# yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima.cogarch model")
- if(any(is.na(oo)))
- yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima.cogarch model")
-
acoeff <- param[ma.name]
b <- param[ar.name]
cost<- param[loc.par]
@@ -73,19 +73,24 @@
p<-length(acoeff)
a[1:p,1] <- acoeff
B <- MatrixA(b[c(q:1)])
- DeltaG <- diff(Data)
+ DeltaG <- c(0,diff(Data))
squaredG <- DeltaG^2
Process_Y <- ExpY0
+# Process_Y <- as.matrix(50.33)
var_V<-cost + sum(acoeff*Process_Y)
delta <- 1/freq
for(t in c(2:(length(Data)))){
# Y_t=e^{A\Delta t}Y_{t-\Delta t}+e^{A\left(\Delta t\right)}e\left(\Delta G_{t}\right)^{2}
- Process_Y <- cbind(Process_Y, (expm(B*delta)%*%(Process_Y[,t-1]+e*squaredG[t-1])))
- var_V[t] <- cost + sum(a*Process_Y[,t])
+ Process_Y <- cbind(Process_Y, (expm(B*delta)%*%(Process_Y[,t-1]+e*squaredG[t])))
+ #Process_Y <- cbind(Process_Y, (Process_Y[,t-1]+delta*B%*%Process_Y[,t-1]+e*squaredG[t]))
+# sim[t,3:ncolsim]<-sim[t-1,3:ncolsim]+(AMatrix*Delta)%*%sim[t-1,3:ncolsim]+evect*sim[t-1,2]*incr.L[2,t-1]
+# sim[t,2]<-value.a0+tavect%*%sim[t-1,3:ncolsim]
+# sim[t,1]<-sim[t-1,1]+sqrt(sim[t,2])*incr.L[1,t]
+ var_V[t] <- cost + t(a)%*%Process_Y[,t-1]
}
#\Delta L_{t}=\frac{\Delta G_{t}}{\sqrt{V_{t}}}.
- incr.L<- DeltaG/sqrt(var_V[c(2:(length(Data)))])
+ incr.L<- DeltaG/sqrt(var_V)
return(incr.L)
}
\ No newline at end of file
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/R/simulate.R 2015-04-25 16:31:21 UTC (rev 376)
@@ -481,35 +481,42 @@
model<-yuimaCogarch at model
info<-model at info
samp <- yuimaCogarch at sampling
+ if(model at measure.type=="CP" && !is.null(increment.L)){
+ method="euler"
+ }
+
if(method=="euler"||(method=="mixed" && model at measure.type=="code")){
- aux.Noise<-setModel(drift="0",
- diffusion="0",
- jump.coeff="1",
- measure=info at measure,
- measure.type=info at measure.type)
+ if(length(increment.L)==0){
+ aux.Noise<-setModel(drift="0",
+ diffusion="0",
+ jump.coeff="1",
+ measure=info at measure,
+ measure.type=info at measure.type)
+
+
+ # aux.samp<-setSampling(Initial = samp at Initial, Terminal = samp at Terminal[1], n = samp at n[1], delta = samp at delta,
+ # grid=samp at grid, random = samp at random, sdelta=samp at sdelta,
+ # sgrid=samp at sgrid, interpolation=samp at interpolation )
- # aux.samp<-setSampling(Initial = samp at Initial, Terminal = samp at Terminal[1], n = samp at n[1], delta = samp at delta,
- # grid=samp at grid, random = samp at random, sdelta=samp at sdelta,
- # sgrid=samp at sgrid, interpolation=samp at interpolation )
-
- aux.samp<-setSampling(Initial = samp at Initial,
- Terminal = samp at Terminal[1],
- n = samp at n[1])
- auxModel<-setYuima(model=aux.Noise, sampling= aux.samp)
-
- if(length(model at parameter@measure)==0){
- aux.incr2<-aux.simulate(object=auxModel, nsim=nsim, seed=seed,
- space.discretized=space.discretized, increment.W=increment.W,
- increment.L=increment.L,
- hurst=0.5,methodfGn=methodfGn)
- }else{
- aux.incr2<-aux.simulate(object=auxModel, nsim=nsim, seed=seed,
- true.parameter = true.parameter[model at parameter@measure],
- space.discretized=space.discretized, increment.W=increment.W,
- increment.L=increment.L,
- hurst=0.5,methodfGn=methodfGn)
- }
- increment<-diff(as.numeric(get.zoo.data(aux.incr2)[[1]]))
+ aux.samp<-setSampling(Initial = samp at Initial,
+ Terminal = samp at Terminal[1],
+ n = samp at n[1])
+ auxModel<-setYuima(model=aux.Noise, sampling= aux.samp)
+
+ if(length(model at parameter@measure)==0){
+ aux.incr2<-aux.simulate(object=auxModel, nsim=nsim, seed=seed,
+ space.discretized=space.discretized, increment.W=increment.W,
+ increment.L=increment.L,
+ hurst=0.5,methodfGn=methodfGn)
+ }else{
+ aux.incr2<-aux.simulate(object=auxModel, nsim=nsim, seed=seed,
+ true.parameter = true.parameter[model at parameter@measure],
+ space.discretized=space.discretized, increment.W=increment.W,
+ increment.L=increment.L,
+ hurst=0.5,methodfGn=methodfGn)
+ }
+ increment<-diff(as.numeric(get.zoo.data(aux.incr2)[[1]]))
+ } else{increment<-increment.L}
# Using the simulated increment for generating the quadratic variation
# As first step we compute it in a crude way. A more fine approach is based on
# the mpv function.
@@ -600,8 +607,8 @@
# sim[t,2]<-value.a0+tavect%*%sim[t,3:ncolsim]
# sim[t,1]<-sim[t-1,1]+sqrt(sim[t,2])*incr.L[1,t]
# sim[t,3:ncolsim]<-expm(AMatrix*Delta)%*%sim[t-1,3:ncolsim]+expm(AMatrix)%*%evect*sim[t-1,2]*incr.L[2,t]
- sim[t,3:ncolsim]<-sim[t-1,3:ncolsim]+(AMatrix*Delta)%*%sim[t-1,3:ncolsim]+evect*sim[t-1,2]*incr.L[2,t-1]
sim[t,2]<-value.a0+tavect%*%sim[t-1,3:ncolsim]
+ sim[t,3:ncolsim]<-sim[t-1,3:ncolsim]+(AMatrix*Delta)%*%sim[t-1,3:ncolsim]+evect*sim[t-1,2]*incr.L[2,t]
sim[t,1]<-sim[t-1,1]+sqrt(sim[t,2])*incr.L[1,t]
}
X <- ts(sim[-(samp at n[1]+1),])
@@ -670,13 +677,14 @@
if(yuimaCogarch at model@measure.type=="code"){
for(t in c(2:n)){
- sim[t,3:ncolsim] <- value.a0*expm(AMatrix*Delta)%*%evect*incr.L[2,t-1]+
- expm(AMatrix*Delta)%*%(Indent+evect%*%tavect*incr.L[2,t-1])%*%sim[t-1,3:ncolsim]
+
# sim[t,2]<-value.a0+tavect%*%sim[t,3:ncolsim]
# sim[t,1]<-sim[t-1,1]+sqrt(sim[t,2])*incr.L[1,t]
# sim[t,3:ncolsim]<-expm(AMatrix*Delta)%*%sim[t-1,3:ncolsim]+expm(AMatrix)%*%evect*sim[t-1,2]*incr.L[2,t]
# sim[t,3:ncolsim]<-sim[t-1,3:ncolsim]+AMatrix*Delta%*%sim[t-1,3:ncolsim]+evect*sim[t-1,2]*incr.L[2,t-1]
sim[t,2]<-value.a0+tavect%*%sim[t-1,3:ncolsim]
+ sim[t,3:ncolsim] <- value.a0*expm(AMatrix*Delta)%*%evect*incr.L[2,t]+
+ expm(AMatrix*Delta)%*%(Indent+evect%*%tavect*incr.L[2,t])%*%sim[t-1,3:ncolsim]
sim[t,1]<-sim[t-1,1]+sqrt(sim[t,2])*incr.L[1,t]
}
@@ -795,3 +803,48 @@
return(result)
}
+# Simulate method for an object of class cogarch.gmm.incr
+
+setMethod("simulate","cogarch.gmm.incr",
+ function(object, nsim=1, seed=NULL, xinit, ...){
+
+ out <-aux.simulategmm(object=object, nsim=nsim, seed=seed, xinit=xinit, ...)
+# out <- simulate(object = model, nsim = nsim, seed=seed, xinit=xinit,
+# sampling = samp,
+# method = "euler",
+# increment.L = t(as.matrix(c(0,Incr.L))),
+# true.parameter = true.parameter,
+# )
+
+ return(out)
+ }
+ )
+
+aux.simulategmm<-function(object, nsim=1, seed=NULL, xinit, ...){
+ Time<-index(object at Incr.Lev)
+ Incr.L<-coredata(object at Incr.Lev)
+
+ model <- object at model
+ EndT <- Time[length(Time)]
+ numb <- (length(Incr.L)+1)
+ valpar<-coef(object)
+
+ idx <- na.omit(match(names(valpar),model at parameter@xinit))
+ solnam <- model at parameter@xinit[-idx]
+ solval <- as.numeric(Diagnostic.Cogarch(object, display=FALSE)$meanStateVariable)
+
+# solval <-50.33
+
+ names(solval) <- solnam
+
+ true.parameter <- as.list(c(valpar,solval))
+
+ samp <- setSampling(Initial = 0, Terminal = EndT, n = numb)
+ out <- simulate(object = model, nsim = nsim, seed=seed, xinit=xinit,
+ sampling = samp,
+ method = "euler",
+ increment.L = t(as.matrix(c(0,Incr.L))),
+ true.parameter = true.parameter,
+ )
+ return(out)
+}
\ No newline at end of file
Added: pkg/yuima/man/Diagnostic.Cogarch.Rd
===================================================================
--- pkg/yuima/man/Diagnostic.Cogarch.Rd (rev 0)
+++ pkg/yuima/man/Diagnostic.Cogarch.Rd 2015-04-25 16:31:21 UTC (rev 376)
@@ -0,0 +1,89 @@
+\name{Diagnostic.Cogarch}
+\alias{Diagnostic.Cogarch}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Function for checking the statistical properties of the COGARCH(p,q) model
+}
+\description{
+The function check the statistical properties of the COGARCH(p,q) model. We verify if the process has a strict positive stationary variance model.}
+\usage{
+Diagnostic.Cogarch(yuima.cogarch, param = list(), matrixS = NULL, mu = 1, display = TRUE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuima.cogarch}{ an object of class \code{yuima.cogarch}, \code{yuima} or a class \code{cogarch.gmm-class}
+}
+ \item{param}{ a list containing the values of the parameters}
+ \item{matrixS}{ a Square matrix.}
+ \item{mu}{ first moment of the Levy measure.}
+ \item{display}{ a logical variable, if \code{TRUE} the function displays the result in the \code{console}. }
+}
+
+\value{The functon returns a List with entries:
+%% ~Describe the value returned
+%% If it is a LIST, use
+\item{meanVarianceProc }{ Unconditional Stationary mean of the variance process. }
+\item{meanStateVariable}{ Unconditional Stationary mean of the state process.}
+\item{stationary}{ If \code{TRUE}, the COGARCH(p,q) has stationary variance.}
+\item{positivity}{ If \code{TRUE}, the variance process is strictly positive.}
+}
+\author{YUIMA Project Team}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\examples{
+\dontrun{
+# Definition of the COGARCH(1,1) process driven by a Variance Gamma nois:
+param.VG <- list(a1 = 0.038, b1 = 0.053,
+ a0 = 0.04/0.053,lambda = 1, alpha = sqrt(2), beta = 0, mu = 0,
+ x01 = 50.33)
+
+cog.VG <- setCogarch(p = 1, q = 1, work = FALSE,
+ measure=list("rngamma(z, lambda, alpha, beta, mu)"),
+ measure.type = "code",
+ Cogarch.var = "y",
+ V.var = "v", Latent.var="x",
+ XinExpr=TRUE)
+
+# Verify the stationarity and the positivity of th variance process
+
+test <- Diagnostic.Cogarch(cog.VG,param=param.VG)
+show(test)
+
+# Simulate a sample path
+
+set.seed(210)
+
+Term=800
+num=24000
+
+samp.VG <- setSampling(Terminal=Term, n=num)
+
+sim.VG <- simulate(cog.VG,
+ true.parameter=param.VG,
+ sampling=samp.VG,
+ method="euler")
+plot(sim.VG)
+
+# Estimate the model
+
+res.VG <- gmm(sim.VG, start = param.VG, Est.Incr = "IncrPar")
+
+summary(res.VG)
+
+# Check if the estimated COGARCH(1,1) has a positive and stationary variance
+
+test1<-Diagnostic.Cogarch(res.VG)
+show(test1)
+
+# Simulate a COGARCH sample path using the estimated COGARCH(1,1)
+# and the recovered increments of underlying Variance Gamma Noise
+
+esttraj<-simulate(res.VG)
+plot(esttraj)
+
+
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
Modified: pkg/yuima/man/cogarch.gmm.incr.rd
===================================================================
--- pkg/yuima/man/cogarch.gmm.incr.rd 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/man/cogarch.gmm.incr.rd 2015-04-25 16:31:21 UTC (rev 376)
@@ -2,8 +2,9 @@
\docType{class}
\alias{cogarch.gmm.incr-class}
\alias{plot,cogarch.gmm.incr,ANY-method}
-\alias{gmm.cogarch.class}
-\alias{cogarch.gmm.class}
+\alias{gmm.cogarch.incr-class}
+\alias{cogarch.gmm.incr-class}
+\alias{simulate,cogarch.gmm.incr-method}
%%\alias{setSampling,yuima.carma-method}
\title{Class for Generalized Method of Moments Estimation for COGARCH(p,q) model with underlying increments}
@@ -27,6 +28,7 @@
}
\section{Methods}{
\describe{
+ \item{simulate}{simulation method. For more information see \code{\link{simulate}}.}
\item{plot}{Plot method for estimated increment of the noise.}
\item{Methods mle}{All methods for \code{mle-class} are available.}
}
Modified: pkg/yuima/man/gmm.rd
===================================================================
--- pkg/yuima/man/gmm.rd 2015-04-21 10:51:29 UTC (rev 375)
+++ pkg/yuima/man/gmm.rd 2015-04-25 16:31:21 UTC (rev 376)
@@ -15,7 +15,7 @@
\usage{
gmm(yuima, data = NULL, start,
method="BFGS", fixed = list(), lower, upper, lag.max = NULL,
- aggr.G = TRUE, aggregation=TRUE, Est.Incr = "NoIncr", objFun = "L2")
+ equally.spaced = TRUE, aggregation=TRUE, Est.Incr = "NoIncr", objFun = "L2")
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -27,7 +27,7 @@
\item{lower}{a named list for specifying lower bounds of parameters.}
\item{upper}{a named list for specifying upper bounds of parameters.}
\item{lag.max}{maximum lag at which to calculate the theoretical and empirical acf. Default is \code{sqrt{N}} where \code{N} is the number of observation.}
- \item{aggr.G}{Logical variable. If \code{aggr.G = TRUE.}, the function use the returns of COGARCH(P,Q) evaluated at unitary length for the computation of the empirical autocorrelations. If \code{aggr.G = FALSE}, the increments are evaluated on the interval with frequency specified in an object of class \code{\link{yuima.data-class}} that contains the observed time series.}
+ \item{equally.spaced}{Logical variable. If \code{equally.spaced = TRUE.}, the function use the returns of COGARCH(P,Q) evaluated at unitary length for the computation of the empirical autocorrelations. If \code{equally.spaced = FALSE}, the increments are evaluated on the interval with frequency specified in an object of class \code{\link{yuima.data-class}} that contains the observed time series.}
\item{aggregation}{If \code{aggregation=TRUE}, before the estimation of the levy parameters we aggregate the estimated increments}
\item{Est.Incr}{ a string variable, If \code{Est.Incr = "NoIncr"}, default value, \code{gmm} returns an object of class \code{\link{cogarch.gmm-class}} that contains the COGARCH parameters.
If \code{Est.Incr = "Incr"} or \code{Est.Incr = "IncrPar"} the output is an object of class \code{\link{cogarch.gmm.incr-class}}. In the first case the object contains the increments of underlying noise while in the second case also the estimated parameter of levy measure.}
@@ -78,6 +78,9 @@
# We estimate the model
res1 <- gmm(yuima = sim1, start = param)
+
+summary(res1)
+
}
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
More information about the Yuima-commits
mailing list