[Yuima-commits] r508 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 3 19:20:51 CET 2016
Author: lorenzo
Date: 2016-11-03 19:20:50 +0100 (Thu, 03 Nov 2016)
New Revision: 508
Added:
pkg/yuima/R/DiagnosticCarma.R
pkg/yuima/man/Diagnostic.Carma.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/ClassCogarch.R
pkg/yuima/R/MM.COGARCH.R
pkg/yuima/R/qmle.R
pkg/yuima/R/simulateForMapsIntegralAndOperator.R
Log:
Added Diagnostic Carma and Updated summary for Carma and Cogarch model.
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/DESCRIPTION 2016-11-03 18:20:50 UTC (rev 508)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.3.4
+Version: 1.3.4
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/NAMESPACE 2016-11-03 18:20:50 UTC (rev 508)
@@ -202,6 +202,7 @@
export(gmm) # Estimation COGARCH(P,Q) using Method Of Moments
export(cogarchNoise)
export(Diagnostic.Cogarch)
+export(Diagnostic.Carma)
# Methods
export(rand)# random number generator of a Levy process specified by user
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/R/AllClasses.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -10,7 +10,7 @@
drift="character",
jump="character",
measure="character",
-# Insert parameters for starting conditions
+# Insert parameters for starting conditions
xinit="character"
)
)
@@ -104,8 +104,8 @@
e = "numeric"
)
)
-
+
# Class 'yuima'
# this is the principal class of yuima project. It may contain up to
@@ -151,7 +151,9 @@
logLI = "ANY",
TypeI = "ANY",
NumbI = "ANY",
- StatI ="ANY"),
+ StatI ="ANY",
+ model = "yuima.carma",
+ Additional.Info = "ANY"),
contains="summary.mle"
)
@@ -172,9 +174,10 @@
setClass("summary.yuima.qmle",
representation(
model = "yuima.model",
-threshold = "ANY"),
+threshold = "ANY",
+Additional.Info = "ANY"),
contains="summary.mle"
)
# The yuima.carma.qmle extends the S4 class "mle". It contains three slots: Estimated Levy,
-# The description of the carma model and the mle.
\ No newline at end of file
+# The description of the carma model and the mle.
Modified: pkg/yuima/R/ClassCogarch.R
===================================================================
--- pkg/yuima/R/ClassCogarch.R 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/R/ClassCogarch.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -39,7 +39,8 @@
setClass("summary.cogarch.est",
representation(objFun = "ANY",
- objFunVal = "ANY"),
+ objFunVal = "ANY",
+ object = "ANY"),
contains="summary.mle"
)
Added: pkg/yuima/R/DiagnosticCarma.R
===================================================================
--- pkg/yuima/R/DiagnosticCarma.R (rev 0)
+++ pkg/yuima/R/DiagnosticCarma.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -0,0 +1,46 @@
+yuima.PhamBreton.Alg<-function(a){
+ p<-length(a)
+ gamma<-a[p:1]
+ if(p>2){
+ gamma[p]<-a[1]
+ alpha<-matrix(NA,p,p)
+ for(j in 1:p){
+ if(is.integer(as.integer(j)/2)){
+ alpha[p,j]<-0
+ alpha[p-1,j]<-0
+ }else{
+ alpha[p,j]<-a[j]
+ alpha[p-1,j]<-a[j+1]/gamma[p]
+ }
+ }
+ for(n in (p-1):1){
+ gamma[n]<-alpha[n+1,2]-alpha[n,2]
+ for(j in 1:n-1){
+ alpha[n-1,j]<-(alpha[n+1,j+2]-alpha[n,j+2])/gamma[n]
+ }
+ alpha[n-1,n-1]<-alpha[n+1,n+1]/gamma[n]
+ }
+ gamma[1]<-alpha[2,2]
+ }
+ return(gamma)
+}
+
+
+Diagnostic.Carma<-function(carma){
+
+ if(!is(carma at model,"yuima.carma"))
+ yuima.stop("model is not a carma")
+ if(!is(carma,"mle"))
+ yuima.stop("object does not belong
+ to yuima.qmle-class or yuima.carma.qmle-class ")
+ param<-coef(carma)
+ info <- carma at model@info
+ numb.ar<-info at p
+ name.ar<-paste(info at ar.par,c(numb.ar:1),sep="")
+ ar.par<-param[name.ar]
+
+ statCond<-FALSE
+ if(min(yuima.PhamBreton.Alg(ar.par[numb.ar:1]))>=0)
+ statCond<-TRUE
+}
+
Modified: pkg/yuima/R/MM.COGARCH.R
===================================================================
--- pkg/yuima/R/MM.COGARCH.R 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/R/MM.COGARCH.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -1143,7 +1143,8 @@
m2logL = 0,
#model = object at model,
objFun = labFun,
- objFunVal = obj
+ objFunVal = obj,
+ object = object
)
tmp
}
@@ -1176,6 +1177,20 @@
cat("\n",paste0(paste("Log.objFun", object at objFun),":"), object at objFunVal, "\n")
}
#cat("objFun", object at min, "\n")
+ dummy <- Diagnostic.Cogarch(object at object, display = FALSE)
+ info <- object at object@yuima at model@info
+ nameMod <- paste0("Cogarch(",info at p,
+ ",", info at q, ") model:", collapse = "")
+ if(dummy$stationary){
+ cat("\n", nameMod, "Stationary conditions are satisfied.\n")
+ }else{
+ cat("\n", nameMod, "Stationary conditions are not satisfied.\n")
+ }
+ if(dummy$positivity){
+ cat("\n", nameMod, "Variance process is positive.\n")
+ }else{
+ cat("\n", nameMod, "Variance process is not positive.\n")
+ }
}
)
@@ -1199,7 +1214,8 @@
logLI = object at logL.Incr,
TypeI = object at yuima@model at measure.type,
NumbI = length(data),
- StatI =summary(data)
+ StatI =summary(data),
+ object = object
)
tmp
}
@@ -1231,6 +1247,20 @@
cat("\nSummary statistics for increments:\n")
print(object at StatI)
cat("\n")
+ dummy <- Diagnostic.Cogarch(object at object, display = FALSE)
+ info <- object at object@yuima at model@info
+ nameMod <- paste0("Cogarch(",info at p,
+ ",", info at q, ") model:", collapse = "")
+ if(dummy$stationary){
+ cat("\n", nameMod, "Stationary conditions are satisfied.\n")
+ }else{
+ cat("\n", nameMod, "Stationary conditions are not satisfied.\n")
+ }
+ if(dummy$positivity){
+ cat("\n", nameMod, "Variance process is positive.\n")
+ }else{
+ cat("\n", nameMod, "Variance process is not positive.\n")
+ }
}
)
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/R/qmle.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -1059,9 +1059,9 @@
method = method, nobs=yuima.nobs, logL.Incr = NULL)
}else{
if(Est.Incr=="NoIncr"){
- final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ final_res<-new("yuima.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
- method = method, nobs=yuima.nobs)
+ method = method, nobs=yuima.nobs , model=yuima at model)
return(final_res)
}else{
yuima.warn("The variable Est.Incr is not correct. See qmle documentation for the allowed values ")
@@ -2252,33 +2252,8 @@
# return(Res)
# }
-yuima.PhamBreton.Alg<-function(a){
- p<-length(a)
- gamma<-a[p:1]
- if(p>2){
- gamma[p]<-a[1]
- alpha<-matrix(NA,p,p)
- for(j in 1:p){
- if(is.integer(as.integer(j)/2)){
- alpha[p,j]<-0
- alpha[p-1,j]<-0
- }else{
- alpha[p,j]<-a[j]
- alpha[p-1,j]<-a[j+1]/gamma[p]
- }
- }
- for(n in (p-1):1){
- gamma[n]<-alpha[n+1,2]-alpha[n,2]
- for(j in 1:n-1){
- alpha[n-1,j]<-(alpha[n+1,j+2]-alpha[n,j+2])/gamma[n]
- }
- alpha[n-1,n-1]<-alpha[n+1,n+1]/gamma[n]
- }
- gamma[1]<-alpha[2,2]
- }
- return(gamma)
-}
+
#yuima.PhamBreton.Inv<-function(gamma){
# p<-length(gamma)
# a<-gamma[p:1]
@@ -2412,10 +2387,14 @@
{
cmat <- cbind(Estimate = object at coef, `Std. Error` = sqrt(diag(object at vcov)))
m2logL <- 2 * object at min
-
+ Additional.Info <- list()
+ if(is(object at model,"yuima.carma")){
+ Additional.Info <-list(Stationarity = Diagnostic.Carma(object))
+ }
tmp <- new("summary.yuima.qmle", call = object at call, coef = cmat,
m2logL = m2logL,
- model = object at model
+ model = object at model,
+ Additional.Info = Additional.Info
)
tmp
}
@@ -2431,6 +2410,17 @@
cat("\nCoefficients:\n")
print(coef(object))
cat("\n-2 log L:", object at m2logL, "\n")
+ if(length(object at Additional.Info)>0){
+ if(is(object at model,"yuima.carma")){
+ Dummy<-paste0("\nCarma(",object at model@info at p,",",object at model@info at q,")",
+ collapse = "")
+ if(object at Additional.Info$Stationarity){
+ cat(Dummy,"model: Stationary conditions are satisfied.\n")
+ }else{
+ cat(Dummy,"model: Stationary conditions are not satisfied.\n")
+ }
+ }
+ }
}
)
@@ -2501,6 +2491,11 @@
data<-Re(coredata(object at Incr.Lev))
data<- data[!is.na(data)]
+ Additional.Info <- list()
+ if(is(object at model,"yuima.carma")){
+ Additional.Info <-list(Stationarity = Diagnostic.Carma(object))
+ }
+
tmp <- new("summary.yuima.carma.qmle", call = object at call, coef = cmat,
m2logL = m2logL,
MeanI = mean(data),
@@ -2508,7 +2503,9 @@
logLI = object at logL.Incr,
TypeI = object at model@measure.type,
NumbI = length(data),
- StatI =summary(data)
+ StatI = summary(data),
+ Additional.Info = Additional.Info,
+ model = object at model
)
tmp
}
@@ -2533,6 +2530,17 @@
cat("\nSummary statistics for increments:\n")
print(object at StatI)
cat("\n")
+ if(length(object at Additional.Info)>0){
+ if(is(object at model,"yuima.carma")){
+ Dummy<-paste0("\nCarma(",object at model@info at p,",",object at model@info at q,")",
+ collapse = "")
+ if(object at Additional.Info$Stationarity){
+ cat(Dummy,"model: Stationary conditions are satisfied.\n")
+ }else{
+ cat(Dummy,"model: Stationary conditions are not satisfied.\n")
+ }
+ }
+ }
}
)
Modified: pkg/yuima/R/simulateForMapsIntegralAndOperator.R
===================================================================
--- pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-11-02 15:54:03 UTC (rev 507)
+++ pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-11-03 18:20:50 UTC (rev 508)
@@ -141,36 +141,50 @@
res <- NULL
PosInMatr <- matrix(c(1:(info.int at dimIntegrand[2]*info.int at dimIntegrand[1])),
nrow = info.int at dimIntegrand[1], ncol = info.int at dimIntegrand[2])
+ my.fun <- function(my.list, my.env, i){
+ dum <- eval(my.list,envir = my.env)
+ if(length(dum)==1){
+ return(rep(dum,i))
+ }else{
+ return(dum[1:i])
+ }
+ }
+ res<-matrix(0,info.int at dimIntegrand[1],(length(time)-1))
+ element <- matrix(0, nrow =info.int at dimIntegrand[1], ncol = 1)
+
for(i in c(1:(length(time)-1))){
assign(info.var at upper.var,time[i+1],envir=my.env)
-# assign("jj",1,envir = my.env)
-# CondW <- paste0(info.var at var.time,"[jj] < ",info.var at upper.var)
-# while(eval(parse(text=CondW),envir = my.env)){
-#
-#
-# # *eval(parse(text = paste0(info.var at var.time,"<",info.var at upper.var)),
-# # envir= my.env)
-# my.env$jj <- my.env$jj+1
-# }
-# Inter <- eval(c(info.int at IntegrandList[[1]]),envir = my.env)[1:i]
- my.fun <- function(my.list, my.env, i){
- dum <- eval(my.list,envir = my.env)
- if(length(dum)==1){
- return(rep(dum,i))
- }else{
- return(dum[1:i])
- }
- }
-
Inter2 <- lapply(info.int at IntegrandList,
FUN = my.fun, my.env = my.env, i = i)
- element <- matrix(0, nrow =info.int at dimIntegrand[1], ncol = 1)
for(j in c(1:info.int at dimIntegrand[1])){
element[j,] <- sum(diag(matrix(unlist(Inter2[PosInMatr[j,]]),
ncol = info.int at dimIntegrand[2])%*%M.dX[,c(1:i)]))
}
- res <- cbind(res, element)
+ res[,i] <- element
}
+ # LengTime<-(length(time)-1)
+ # my.listenv<-as.list(c(1:LengTime))
+ # names(my.listenv)<- rep("i",LengTime)
+ # globalMyEnv <-new.env()
+ # globalMyEnv$time <- time
+ # globalMyEnv$my.env <- my.env
+ # my.listenv2<-lapply(X=my.listenv,
+ # globalMyEnv=globalMyEnv,
+ # FUN = function(X,globalMyEnv){
+ # assign(info.var at upper.var,globalMyEnv$time[X+1],
+ # envir= globalMyEnv$my.env)
+ # Inter2 <- lapply(info.int at IntegrandList,
+ # FUN = my.fun, my.env = globalMyEnv$my.env,
+ # i = X)
+ # for(j in c(1:info.int at dimIntegrand[1])){
+ # element[j,] <- sum(diag(matrix(unlist(Inter2[PosInMatr[j,]]),
+ # ncol = info.int at dimIntegrand[2])%*%M.dX[,c(1:X)]))
+ # }
+ #
+ # list(element)
+ # })
+ # res<-as.numeric(unlist(my.listenv2))
+ # res<-matrix(res,info.int at dimIntegrand[1],(length(time)-1))
res <- cbind(0,res)
rownames(res) <- object at Integral@variable.Integral at out.var
my.data <- zoo(x = t(res), order.by = time)
Added: pkg/yuima/man/Diagnostic.Carma.Rd
===================================================================
--- pkg/yuima/man/Diagnostic.Carma.Rd (rev 0)
+++ pkg/yuima/man/Diagnostic.Carma.Rd 2016-11-03 18:20:50 UTC (rev 508)
@@ -0,0 +1,43 @@
+\name{Diagnostic.Carma}
+\alias{Diagnostic.Carma}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Diagnostic Carma model}
+\description{
+This function verifies if the condition of stationarity is satisfied.}
+\usage{
+Diagnostic.Carma(carma)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{carma}{ An object of class \code{yuima.qmle-class} where the slot \code{model} is a carma process.
+}
+}
+\value{Logical variable. If \code{TRUE}, Carma is stationary.}
+\author{
+YUIMA TEAM
+}
+
+
+\examples{
+
+mod1 <- setCarma(p = 2, q = 1, scale.par = "sig",
+ Carma.var = "y")
+
+param1 <- list(a1 = 1.39631, a2 = 0.05029, b0 = 1,
+ b1 = 1, sig = 1)
+samp1 <- setSampling(Terminal = 100, n = 200)
+
+set.seed(123)
+
+sim1 <- simulate(mod1, true.parameter = param1,
+ sampling = samp1)
+
+est1 <- qmle(sim1, start = param1)
+
+Diagnostic.Carma(est1)
+
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ ~kwd1 }% use one of RShowDoc("KEYWORDS")
+\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
More information about the Yuima-commits
mailing list