[Yuima-commits] r393 - in pkg/yuima: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 10 08:36:25 CEST 2015
Author: kyuta
Date: 2015-10-10 08:36:25 +0200 (Sat, 10 Oct 2015)
New Revision: 393
Added:
pkg/yuima/R/mllag.R
pkg/yuima/R/spectralcov.R
pkg/yuima/man/mllag.Rd
pkg/yuima/man/spectralcov.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/R/llag.R
pkg/yuima/man/cce.Rd
pkg/yuima/man/llag.Rd
pkg/yuima/man/mpv.Rd
pkg/yuima/man/noisy.sampling.Rd
pkg/yuima/src/cce_functions.c
Log:
modified: llag.R, cce.Rd, llag.Rd, mpv.Rd, noisy.sampling.Rd, cce_functions.c
added: mllag.R, spectralcov.R, mllag.Rd, spectralcov.Rd
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2015-09-01 08:45:31 UTC (rev 392)
+++ pkg/yuima/DESCRIPTION 2015-10-10 06:36:25 UTC (rev 393)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.0.74
+Version: 1.0.75
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/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2015-09-01 08:45:31 UTC (rev 392)
+++ pkg/yuima/NAMESPACE 2015-10-10 06:36:25 UTC (rev 393)
@@ -1,188 +1,197 @@
-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)
-
-
-# 03/07/2015
-importFrom(stats, time)
-importFrom(stats, ts)
-importFrom(stats, rnorm)
-importFrom(stats, na.omit)
-importFrom(stats, dgamma)
-importFrom(stats, optimHess)
-importFrom(stats, filter)
-importFrom(utils, tail)
-importFrom(utils, head)
-importFrom(stats, acf)
-importFrom(stats, fft)
-importFrom(stats, rexp)
-importFrom(stats, approx)
-importFrom(stats, arima0)
-importFrom(stats, frequency)
-importFrom(stats, D)
-importFrom(stats, integrate)
-importFrom(stats, rpois)
-importFrom(stats, runif)
-importFrom(stats, optim)
-importFrom(stats, optimize)
-importFrom(stats, deltat)
-importFrom(stats, pchisq)
-importFrom(stats, symnum)
-importFrom(stats, rchisq)
-importFrom(stats, rgamma)
-importFrom(stats, diffinv)
-importFrom(stats, pnorm)
-importFrom(stats, approxfun)
-importFrom(stats, qnorm)
-importFrom(stats, rbinom)
-importFrom(stats, constrOptim)
-importFrom(stats, dnorm)
-importFrom(stats, deriv)
-importFrom(graphics, points)
-importFrom(stats, end)
-importFrom(stats, start)
-importFrom(utils, str)
-importFrom(stats, sd)
-
-
-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)
-
-
+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)
+
+
+# 03/07/2015
+importFrom(stats, time)
+importFrom(stats, ts)
+importFrom(stats, rnorm)
+importFrom(stats, na.omit)
+importFrom(stats, dgamma)
+importFrom(stats, optimHess)
+importFrom(stats, filter)
+importFrom(utils, tail)
+importFrom(utils, head)
+importFrom(stats, acf)
+importFrom(stats, fft)
+importFrom(stats, rexp)
+importFrom(stats, approx)
+importFrom(stats, arima0)
+importFrom(stats, frequency)
+importFrom(stats, D)
+importFrom(stats, integrate)
+importFrom(stats, rpois)
+importFrom(stats, runif)
+importFrom(stats, optim)
+importFrom(stats, optimize)
+importFrom(stats, deltat)
+importFrom(stats, pchisq)
+importFrom(stats, symnum)
+importFrom(stats, rchisq)
+importFrom(stats, rgamma)
+importFrom(stats, diffinv)
+importFrom(stats, pnorm)
+importFrom(stats, approxfun)
+importFrom(stats, qnorm)
+importFrom(stats, rbinom)
+importFrom(stats, constrOptim)
+importFrom(stats, dnorm)
+importFrom(stats, deriv)
+importFrom(graphics, points)
+importFrom(stats, end)
+importFrom(stats, start)
+importFrom(utils, str)
+importFrom(stats, sd)
+
+
+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(lmm) # Oct. 10, 2015: local methods of moment estimator
+export(mllag) # Oct. 10, 2015: multiple lead-lag detector
+
+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(print, yuima.llag) # Oct. 10, 2015
+S3method(print, yuima.mllag) # Oct. 10, 2015
+S3method(print, yuima.specv) # Oct. 10, 2015
+
+
+S3method(toLatex, yuima)
+S3method(toLatex, yuima.model)
+S3method(toLatex, yuima.carma)
+S3method(toLatex, yuima.cogarch)
+
+S3method(plot, yuima.llag) # Oct. 10, 2015
+S3method(plot, yuima.mllag) # Oct. 10, 2015
+
+useDynLib(yuima)
+
+
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2015-09-01 08:45:31 UTC (rev 392)
+++ pkg/yuima/NEWS 2015-10-10 06:36:25 UTC (rev 393)
@@ -34,4 +34,7 @@
modified cce.Rd, cce_functions.c
2015/05/14: fixed a bug in cce_functions.c
2015/09/01: fixed a bug in cce_functions.c
- add some contents to the examples of cce
\ No newline at end of file
+ add some contents to the examples of cce
+2015/10/10: modified llag.R, cce.Rd, llag.Rd, mpv.Rd, noisy.sampling.Rd, cce_functions.c
+ added mllag.R, mllag.Rd (multiple lead-lag detector);
+ spectralcov.R, spectralcov.Rd (spectral covariance estimator)
\ No newline at end of file
Modified: pkg/yuima/R/llag.R
===================================================================
--- pkg/yuima/R/llag.R 2015-09-01 08:45:31 UTC (rev 392)
+++ pkg/yuima/R/llag.R 2015-10-10 06:36:25 UTC (rev 393)
@@ -14,6 +14,532 @@
# llag(x at data,from=FALSE,to=FALSE,division=FALSE,verbose=verbose ))
#setMethod( "llag", "yuima.data", function(x,from=FALSE,to=FALSE,division=FALSE,verbose=FALSE) {
+
+## function to make the grid if it is missing
+make.grid <- function(d, d.size, x, from, to, division){
+
+ out <- vector(d.size, mode = "list")
+
+ if(length(from) != d.size){
+ from <- c(from,rep(-Inf,d.size - length(from)))
+ }
+
+ if(length(to) != d.size){
+ to <- c(to,rep(Inf,d.size - length(to)))
+ }
+
+ if(length(division) == 1){
+ division <- rep(division,d.size)
+ }
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ time1 <- as.numeric(time(x[[i]]))
+ time2 <- as.numeric(time(x[[j]]))
+
+ # calculate the maximum of correlation by substituting theta to lagcce
+
+ #n:=2*length(data)
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+
+ if(division[num]==FALSE){
+ n <- round(2*max(length(time1),length(time2)))+1
+ }else{
+ n <- division[num]
+ }
+
+ # maximum value of lagcce
+
+ tmptheta <- time2[1]-time1[1] # time lag (sec)
+
+ num1 <- time1[length(time1)]-time1[1] # elapsed time for series 1
+ num2 <- time2[length(time2)]-time2[1] # elapsed time for series 2
+
+ # modified
+
+ if(is.numeric(from[num])==TRUE && is.numeric(to[num])==TRUE){
+ num2 <- min(-from[num],num2+tmptheta)
+ num1 <- min(to[num],num1-tmptheta)
+ tmptheta <- 0
+
+ if(-num2 >= num1){
+ print("llag:invalid range")
+ return(NULL)
+ }
+ }else if(is.numeric(from[num])==TRUE){
+ num2 <- min(-from[num],num2+tmptheta)
+ num1 <- num1-tmptheta
+ tmptheta <- 0
+
+ if(-num2 >= num1){
+ print("llag:invalid range")
+ return(NULL)
+ }
+ }else if(is.numeric(to[num])==TRUE){
+ num2 <- num2+tmptheta
+ num1 <- min(to[num],num1-tmptheta)
+ tmptheta <- 0
+
+ if(-num2 >= num1){
+ print("llag:invalid range")
+ return(NULL)
+ }
+ }
+
+ out[[num]] <- seq(-num2-tmptheta,num1-tmptheta,length=n)[2:(n-1)]
+
+ }
+ }
+
+ return(out)
+}
+
+
+## function to compute asymptotic variances
+llag.avar <- function(x, grid, bw, alpha, fisher, ser.diffX, ser.times, vol, cormat, ccor, idx, G, d, d.size){
+
+ # treatment of the bandwidth
+ if(missing(bw)){
+
+ bw <- matrix(0, d, d)
+
+ for(i in 1:(d - 1)){
+
+ for(j in (i + 1):d){
+
+ Init <- min(ser.times[[i]][1], ser.times[[j]][1])
+ Term <- max(tail(ser.times[[i]], n = 1), tail(ser.times[[j]], n = 1))
+ bw[i, j] <- (Term - Init) * min(length(ser.times[[i]]), length(ser.times[[j]]))^(-0.45)
+ bw[j, i] <- bw[i, j]
+
+ }
+
+ }
+
+ }
+
+ bw <- matrix(bw, d, d)
+
+ p <- diag(d) # p-values
+ avar <- vector(d.size,mode="list") # asymptotic variances
+ CI <- vector(d.size,mode="list") # confidence intervals
+
+ names(avar) <- names(ccor)
+ names(CI) <- names(ccor)
+
+ vv <- (2/3) * sapply(ser.diffX, FUN = function(x) sum(x^4)) # AVAR for RV
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ time1 <- ser.times[[i]]
+ time2 <- ser.times[[j]]
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+
+ ## computing conficence intervals
+ N.max <- max(length(time1), length(time2))
+ tmp <- rev(as.numeric(ccor[[num]]))
+
+ d1 <- -0.5 * tmp * sqrt(vol[j]/vol[i])
+ d2 <- -0.5 * tmp * sqrt(vol[i]/vol[j])
+
+ avar.tmp <- .C("hycrossavar",
+ as.double(G[[num]]),
+ as.double(time1),
+ as.double(time2),
+ as.integer(length(G[[num]])),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(x[[i]]),
+ as.double(x[[j]]),
+ as.integer(N.max),
+ as.double(bw[i, j]),
+ as.double(vv[i]),
+ as.double(vv[j]),
+ as.double(d1^2),
+ as.double(d2^2),
+ as.double(2 * d1),
+ as.double(2 * d2),
+ as.double(2 * d1 * d2),
+ covar = double(length(G[[num]])),
+ corr = double(length(G[[num]])))$corr / (vol[i] * vol[j])
+
+ #avar[[num]][avar[[num]] <= 0] <- NA
+
+ if(fisher == TRUE){
+ z <- atanh(tmp) # the Fisher z transformation of the estimated correlation
+ se.z <- sqrt(avar.tmp)/(1 - tmp^2)
+ p[i,j] <- pchisq((atanh(cormat[i,j])/se.z[idx[num]])^2, df=1, lower.tail=FALSE)
+ CI[[num]] <- zoo(tanh(qnorm(1 - alpha/2) * se.z), -grid[[num]])
+ }else{
+ p[i,j] <- pchisq(cormat[i,j]^2/avar.tmp[idx[num]], df=1, lower.tail=FALSE)
+ c.alpha <- sqrt(qchisq(alpha, df=1, lower.tail = FALSE) * avar.tmp)
+ CI[[num]] <- zoo(c.alpha, -grid[[num]])
+ }
+
+ p[j,i] <- p[i,j]
+ avar[[num]] <- zoo(avar.tmp, -grid[[num]])
+
+ }
+ }
+
+ return(list(p = p, avar = avar, CI = CI))
+}
+
+
+## main body
+setGeneric( "llag", function(x, from = -Inf, to = Inf, division = FALSE,
+ verbose = (ci || ccor), grid, psd = TRUE, plot = ci,
+ ccor = ci, ci = FALSE, alpha = 0.01, fisher = TRUE, bw) standardGeneric("llag") )
+
+## yuima-method
+setMethod("llag", "yuima", function(x, from, to, division, verbose, grid, psd, plot,
+ ccor, ci, alpha, fisher, bw)
+ llag(x at data, from, to, division, verbose, grid, psd, plot, ccor, ci, alpha, fisher, bw))
+
+## yuima.data-method
+setMethod("llag", "yuima.data", function(x, from, to, division, verbose, grid, psd, plot,
+ ccor, ci, alpha, fisher, bw)
+ llag(x at zoo.data, from, to, division, verbose, grid, psd, plot, ccor, ci, alpha, fisher, bw))
+
+## list-method
+setMethod("llag", "list", function(x, from, to, division, verbose, grid, psd, plot,
+ ccor, ci, alpha, fisher, bw) {
+
+ d <- length(x)
+
+ # allocate memory
+ ser.times <- vector(d, mode="list") # time index in 'x'
+ ser.diffX <- vector(d, mode="list") # difference of data
+ vol <- double(d)
+
+ # Set the tolerance to avoid numerical erros in comparison
+ tol <- 1e-6
+
+ for(i in 1:d){
+
+ # NA data must be skipped
+ idt <- which(is.na(x[[i]]))
+ if(length(idt>0)){
+ x[[i]] <- x[[i]][-idt]
+ }
+ if(length(x[[i]])<2) {
+ stop("length of data (w/o NA) must be more than 1")
+ }
+
+ # set data and time index
+ ser.times[[i]] <- as.numeric(time(x[[i]]))/tol
+ # set difference of the data
+ ser.diffX[[i]] <- diff( as.numeric(x[[i]]) )
+ vol[i] <- sum(ser.diffX[[i]]^2)
+ }
+
+ theta <- matrix(0,d,d)
+
+ d.size <- d*(d-1)/2
+ crosscor <- vector(d.size,mode="list")
+ idx <- integer(d.size)
+
+ # treatment of the grid
+ if(missing(grid))
+ grid <- make.grid(d, d.size, x, from, to, division)
+
+ if(is.list(grid)){
+ G <- relist(unlist(grid)/tol, grid)
+ }else{
+ if(is.numeric(grid)){
+ G <- data.frame(matrix(grid/tol,length(grid),d.size))
+ grid <- data.frame(matrix(grid,length(grid),d.size))
+ }else{
+ print("llag:invalid grid")
+ return(NULL)
+ }
+ }
+
+ # core part
+ if(psd){ # positive semidefinite correction is implemented
+
+ cormat <- diag(d)
+ LLR <- diag(d)
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ time1 <- ser.times[[i]]
+ time2 <- ser.times[[j]]
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+
+ names(crosscor)[num] <- paste("(",i,",",j,")", sep = "")
+
+ tmp <- .C("HYcrosscorr",
+ as.integer(length(G[[num]])),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(G[[num]]),
+ as.double(time1),
+ as.double(time2),
+ double(length(time2)),
+ as.double(ser.diffX[[i]]),
+ as.double(ser.diffX[[j]]),
+ as.double(vol[i]),
+ as.double(vol[j]),
+ value=double(length(G[[num]])))$value
+
+ idx[num] <- which.max(abs(tmp))
+ mlag <- -grid[[num]][idx[num]] # make the first timing of max or min
+ cor <- tmp[idx[num]]
+
+ theta[i,j] <- mlag
+ cormat[i,j] <- cor
+ theta[j,i] <- -mlag
+ cormat[j,i] <- cormat[i,j]
+
+ LLR[i,j] <- sum(tmp[grid[[num]] < 0]^2)/sum(tmp[grid[[num]] > 0]^2)
+ LLR[j,i] <- 1/LLR[i,j]
+
+ crosscor[[num]] <- zoo(tmp,-grid[[num]])
+ }
+ }
+
+ covmat <- diag(sqrt(vol))%*%cormat%*%diag(sqrt(vol))
+
+ }else{# non-psd
+
+ covmat <- diag(vol)
+ LLR <- diag(d)
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ time1 <- ser.times[[i]]
+ time2 <- ser.times[[j]]
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+
+ names(crosscor)[num] <- paste("(",i,",",j,")", sep = "")
+
+ tmp <- .C("HYcrosscov",
+ as.integer(length(G[[num]])),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(G[[num]]),
+ as.double(time1),
+ as.double(time2),
+ double(length(time2)),
+ as.double(ser.diffX[[i]]),
+ as.double(ser.diffX[[j]]),
+ value=double(length(G[[num]])))$value
+
+ idx[num] <- which.max(abs(tmp))
+ mlag <- -grid[[num]][idx[num]] # make the first timing of max or min
+ cov <- tmp[idx[num]]
+
+ theta[i,j] <- mlag
+ covmat[i,j] <- cov
+ theta[j,i] <- -mlag
+ covmat[j,i] <- covmat[i,j]
+
+ LLR[i,j] <- sum(tmp[grid[[num]] < 0]^2)/sum(tmp[grid[[num]] > 0]^2)
+ LLR[j,i] <- 1/LLR[i,j]
+
+ crosscor[[num]] <- zoo(tmp,-grid[[num]])/sqrt(vol[i]*vol[j])
+ }
+ }
+
+ cormat <- diag(1/sqrt(diag(covmat)))%*%covmat%*%diag(1/sqrt(diag(covmat)))
+
+ }
+
+ if(ci){ # computing confidence intervals
+
+ out <- llag.avar(x = x, grid = grid, bw = bw, alpha = alpha, fisher = fisher,
+ ser.diffX = ser.diffX, ser.times = ser.times,
+ vol = vol, cormat = cormat, ccor = crosscor, idx = idx,
+ G = G, d = d, d.size = d.size)
+
+ p <- out$p
+ CI <- out$CI
+ avar <- out$avar
+
+ colnames(theta) <- names(x)
+ rownames(theta) <- names(x)
+ colnames(p) <- names(x)
+ rownames(p) <- names(x)
+
+ if(plot){
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ y.max <- max(abs(as.numeric(crosscor[[num]])),as.numeric(CI[[num]]))
+
+ plot(crosscor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)),
+ ylim=c(-y.max,y.max))
+
+ lines(CI[[num]],lty=2,col="blue")
+ lines(-CI[[num]],lty=2,col="blue")
+ }
+ }
+ }
+
+ if(verbose==TRUE){
+
+ colnames(covmat) <- names(x)
+ rownames(covmat) <- names(x)
+ colnames(cormat) <- names(x)
+ rownames(cormat) <- names(x)
+ colnames(LLR) <- names(x)
+ rownames(LLR) <- names(x)
+
+ if(ccor){
+ result <- list(lagcce=theta,p.values=p,covmat=covmat,cormat=cormat,LLR=LLR,ccor=crosscor,avar=avar)
+ }else{
+ result <- list(lagcce=theta,p.values=p,covmat=covmat,cormat=cormat,LLR=LLR)
+ }
+
+ }else{
+ return(theta)
+ }
+
+ }else{
+
+ colnames(theta) <- names(x)
+ rownames(theta) <- names(x)
+
+ if(plot){
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ plot(crosscor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)))
+ }
+ }
+ }
+
+ if(verbose==TRUE){
+
+ colnames(covmat) <- names(x)
+ rownames(covmat) <- names(x)
+ colnames(cormat) <- names(x)
+ rownames(cormat) <- names(x)
+ colnames(LLR) <- names(x)
+ rownames(LLR) <- names(x)
+
+ if(ccor){
+ result <- list(lagcce=theta,covmat=covmat,cormat=cormat,LLR=LLR,ccor=crosscor)
+ }else{
+ result <- list(lagcce=theta,covmat=covmat,cormat=cormat,LLR=LLR)
+ }
+ }else{
+ return(theta)
+ }
+
+ }
+
+ class(result) <- "yuima.llag"
+
+ return(result)
+})
+
+# print method for yuima.llag-class
+print.yuima.llag <- function(x, ...){
+
+ if(is.null(x$p.values)){
+
+ cat("Estimated lead-lag parameters\n")
+ print(x$lagcce, ...)
+ cat("Corresponding covariance matrix\n")
+ print(x$covmat, ...)
+ cat("Corresponding correlation matrix\n")
+ print(x$cormat, ...)
+ cat("Lead-lag ratio\n")
+ print(x$LLR, ...)
+
+ }else{
+
+ cat("Estimated lead-lag parameters\n")
+ print(x$lagcce, ...)
+ cat("Corresponding p-values\n")
+ print(x$p.values, ...)
+ cat("Corresponding covariance matrix\n")
+ print(x$covmat, ...)
+ cat("Corresponding correlation matrix\n")
+ print(x$cormat, ...)
+ cat("Lead-lag ratio\n")
+ print(x$LLR, ...)
+
+ }
+
+}
+
+# plot method for yuima.llag-class
+plot.yuima.llag <- function(x, alpha = 0.01, fisher = TRUE, ...){
+
+ if(is.null(x$ccor)){
+ warning("cross-correlation functions were not returned by llag. Set verbose = TRUE and ccor = TRUE to return them.",
+ call. = FALSE)
+ return(NULL)
+ }else{
+
+ d <- nrow(x$LLR)
+
+ if(is.null(x$avar)){
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+
+ plot(x$ccor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)))
+
+ }
+ }
+
+ }else{
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ G <- time(x$ccor[[num]])
+
+ if(fisher == TRUE){
+ z <- atanh(x$ccor[[num]]) # the Fisher z transformation of the estimated correlation
+ se.z <- sqrt(x$avar[[num]])/(1 - x$ccor[[num]]^2)
+ CI <- zoo(tanh(qnorm(1 - alpha/2) * se.z), G)
+ }else{
+ c.alpha <- sqrt(qchisq(alpha, df=1, lower.tail = FALSE) * x$avar[[num]])
+ CI <- zoo(c.alpha, G)
+ }
+
+ y.max <- max(abs(as.numeric(x$ccor[[num]])),as.numeric(CI))
+
+ plot(x$ccor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)),
+ ylim=c(-y.max,y.max))
+
+ lines(CI,lty=2,col="blue")
+ lines(-CI,lty=2,col="blue")
+
+ }
+ }
+ }
+
+ }
+
+}
+
+
+## Old version until Oct. 10, 2015
+if(0){
setGeneric( "llag",
function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid,psd=TRUE,
plot=FALSE,ccor=FALSE)
@@ -396,7 +922,7 @@
return(theta)
}
})
-
+}
## Old version
if(0){
Added: pkg/yuima/R/mllag.R
===================================================================
--- pkg/yuima/R/mllag.R (rev 0)
+++ pkg/yuima/R/mllag.R 2015-10-10 06:36:25 UTC (rev 393)
@@ -0,0 +1,135 @@
+
+# multiple lead-lag detector
+
+mllag <- function(x, from = -Inf, to = Inf, division = FALSE, grid, psd = TRUE,
+ plot = TRUE, alpha = 0.01, fisher = TRUE, bw) {
+
+ if((is(x) == "yuima")||(is(x) == "yuima.data")||(is(x) == "list")){
+ x <- llag(x, from = from, to = to, division = division, grid = grid, psd = psd, plot = FALSE,
+ ci = TRUE, alpha = alpha, fisher = fisher, bw = bw)
+ }else if((is(x) != "yuima.llag") && (is(x) != "yuima.mllag")){
+ cat("mllag: invalid x")
+ return(NULL)
+ }
+
+ d <- ncol(x$LLR)
+ d.size <- length(x$ccor)
+
+ CI <- vector(d.size,mode="list") # confidence intervals
+ result <- vector(d.size,mode="list")
+
+ names(CI) <- names(x$ccor)
+ names(result) <- names(x$ccor)
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ G <- time(x$ccor[[num]])
+ tmp <- x$ccor[[num]]
+ avar.tmp <- x$avar[[num]]
+
+ if(fisher == TRUE){
+
+ z <- atanh(tmp) # the Fisher z transformation of the estimated correlation
+ se.z <- sqrt(avar.tmp)/(1 - tmp^2)
+ c.alpha <- tanh(qnorm(1 - alpha/2) * se.z)
+
+ idx <- which(abs(tmp) > c.alpha)
+
+ result[[num]] <- data.frame(lagcce = G[idx],
+ p.values = pchisq(atanh(tmp[idx])^2/se.z[idx]^2, df=1, lower.tail=FALSE),
+ correlation = tmp[idx])
+
+ }else{
+
+ c.alpha <- sqrt(qchisq(alpha, df=1, lower.tail=FALSE) * avar.tmp)
+
+ idx <- which(abs(tmp) > c.alpha)
+
+ result[[num]] <- data.frame(lagcce = G[idx],
+ p.value = pchisq(tmp[idx]^2/avar.tmp[idx], df=1, lower.tail=FALSE),
+ correlation = tmp[idx])
+
+ }
+
+ CI[[num]] <- zoo(c.alpha, G)
+
+ }
+ }
+
+ if(plot){
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ y.max <- max(abs(as.numeric(x$ccor[[num]])),as.numeric(CI[[num]]))
+
+ plot(x$ccor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)),
+ ylim=c(-y.max,y.max))
+
+ lines(CI[[num]],lty=2,col="blue")
+ lines(-CI[[num]],lty=2,col="blue")
+
+ points(result[[num]]$lagcce, result[[num]]$correlation, col = "red", pch = 19)
+ }
+ }
+ }
+
+ out <- list(mlagcce = result, LLR = x$LLR, ccor = x$ccor, avar = x$avar, CI = CI)
+
+ class(out) <- "yuima.mllag"
+
+ return(out)
+}
+
+
+# print method for yuima.mllag-class
+print.yuima.mllag <- function(x, ...){
+
+ cat("Estimated lead-lag parameters\n")
+ print(x$mlagcce, ...)
+ cat("Lead-lag ratio\n")
+ print(x$LLR, ...)
+
+}
+
+
+# plot method for yuima.mllag-class
+plot.yuima.mllag <- function(x, alpha = 0.01, fisher = TRUE, ...){
+
+ d <- nrow(x$LLR)
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+
+ num <- d*(i-1) - (i-1)*i/2 + (j-i)
+ G <- time(x$ccor[[num]])
+
+ if(fisher == TRUE){
+ z <- atanh(x$ccor[[num]]) # the Fisher z transformation of the estimated correlation
+ se.z <- sqrt(x$avar[[num]])/(1 - x$ccor[[num]]^2)
+ CI <- zoo(tanh(qnorm(1 - alpha/2) * se.z), G)
+ }else{
+ c.alpha <- sqrt(qchisq(alpha, df=1, lower.tail = FALSE) * x$avar[[num]])
+ CI <- zoo(c.alpha, G)
+ }
+
+ y.max <- max(abs(as.numeric(x$ccor[[num]])),as.numeric(CI))
+
+ plot(x$ccor[[num]],
+ main=paste(i,"vs",j,"(positive",expression(theta),"means",i,"leads",j,")"),
+ xlab=expression(theta),ylab=expression(U(theta)),
+ ylim=c(-y.max,y.max))
+
+ lines(CI,lty=2,col="blue")
+ lines(-CI,lty=2,col="blue")
+
+ points(x$result[[num]]$lagcce, x$result[[num]]$correlation, col = "red", pch = 19)
+ }
+ }
+
+}
+
Added: pkg/yuima/R/spectralcov.R
===================================================================
--- pkg/yuima/R/spectralcov.R (rev 0)
+++ pkg/yuima/R/spectralcov.R 2015-10-10 06:36:25 UTC (rev 393)
@@ -0,0 +1,224 @@
+
+## Local Method of Moments estimator
+
+lmm <- function(x, block = 20, freq = 50, freq.p = 10, K = 4, interval = c(0, 1),
+ Sigma.p = NULL, noise.var = "AMZ", samp.adj = "direct", psd = TRUE){
+
+ Data <- get.zoo.data(x)
+
+ d.size <- length(Data)
+ N <- sapply(Data,"length") - 1
+ n <- min(N)
+
+ interval <- as.numeric(interval)
+ freq.p <- min(freq.p,freq)
+
+ if(d.size > 1){
+
+ # Constructing the stpectral statistics and noise level estimates
+ Spec <- array(0,dim=c(d.size,block,freq))
+ H <- matrix(0,d.size,block)
+ cH <- array(0,dim=c(d.size^2,block,freq))
+ zmat <- diag(d.size^2)
+
+ phi.br <- (block*pi*(1:freq))^2
+
+ for(p in 1:d.size){
+
+ dY <- diff(as.numeric(Data[[p]]))
+ Time <- (as.numeric(time(Data[[p]])) - interval[1])/(interval[2] - interval[1])
+
+ if((min(Time) < 0) || (max(Time) > 1)){
+ print("lmm:invalid interval")
+ return(NULL)
+ }
+
+ Time2 <- (Time[1:N[p]]+Time[-1])/2
+
+ k.idx <- findInterval(Time2,vec=seq(0,1,by=1/block))
+
+ # Noise variance
+ if(is.numeric(noise.var)){
+ sq.eta <- noise.var[p]
+ }else{
+ sq.eta <- switch(noise.var,
+ "BR" = mean(dY^2)/2, # Bandi and Russel (2006)
+ "O" = -sum(dY[-N[p]] * dY[-1])/N[p], # Oomen (2006)
+ "AMZ" = {tmp <- arima0(dY,order=c(0,0,1),include.mean=FALSE);
+ -tmp$coef * tmp$sigma2} # Ait-Sahalia et al. (2005)
+ )
+ }
+
+ # Effect of irregularity
+ if(is.numeric(samp.adj)){
+ EOI <- samp.adj[p, ]
+ }else{
+ EOI <- switch(samp.adj,
+ "QVT" = block * rowsum(diff(Time)^2,group = findInterval(Time[-1],vec=seq(0,1,by=1/block),
+ rightmost.closed = TRUE)),
+ "direct" = block * rowsum((diff(Time,lag=2)/2)^2, group = k.idx[-N[p]]))
+ }
+
+ Spec[p,,] <- sqrt(2*block)*
+ rowsum(sin(tcrossprod(pi*(block*Time2-(k.idx-1)),1:freq))*dY,
+ group = k.idx)
+ H[p, ] <- sq.eta * EOI
+ cH[d.size*(p-1)+p,,] <- tcrossprod(H[p,],phi.br)
+
+ for(q in 1:d.size){
+
+ tmp.mat <- diag(0, d.size)
+ tmp.mat[p, q] <- 1
+ zmat <- zmat + tmp.mat %x% t(tmp.mat)
+
+ }
+
+ }
+
+ #summand <- apply(Spec,c(2,3),FUN="tcrossprod") - cH
+ summand <- array(.C("krprod",
+ as.double(Spec),
+ as.integer(d.size),
+ as.integer(block*freq),
+ result=double(d.size^2*block*freq))$result,
+ dim=c(d.size^2,block,freq)) - cH
+
+ if(is.null(Sigma.p)){
+ # Pilot estimation of Sigma
+ # Using K adjacent blocks following Bibinger et al. (2014a)
+ #Sigma.p <- rollmean(t(rowMeans(array(summand[,,1:freq.p],
+ # dim=c(d.size^2,block,freq.p)),
+ # dims = 2)),k = K, fill = "e")
+ # Recent version of Bibinger et al. (2014b)
+ Sigma.p <- rollapply(t(rowMeans(array(summand[,,1:freq.p],
+ dim=c(d.size^2,block,freq.p)),
+ dims = 2)), width = 2*K + 1,
+ FUN="mean", partial = TRUE)
+ }
+
+ # Local method of moments
+ Sigma <- matrix(0,d.size^2,block)
+ inv.Ik <- array(0,dim=c(d.size^2,d.size^2,block))
+
+ for(k in 1:block){
+
+ tmp <- double(d.size^2)
+ H.inv <- diag(1/H[,k])
+ tmp.eigen <- eigen(matrix(Sigma.p[k,], d.size, d.size) %*% H.inv,
+ symmetric=FALSE)
+ inv.V <- solve(tmp.eigen$vectors)
+ iHV <- H.inv %*% tmp.eigen$vectors
+ iHVxiHV <- iHV %x% iHV
+ inv.VxV <- inv.V %x% inv.V
+ #Lambda <- apply(1/outer(tmp.eigen$values, phi.br, FUN="+"),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 393
More information about the Yuima-commits
mailing list