[Yuima-commits] r317 - in pkg/yuima: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 7 04:19:12 CEST 2014
Author: kyuta
Date: 2014-07-07 04:19:11 +0200 (Mon, 07 Jul 2014)
New Revision: 317
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/llag.R
pkg/yuima/man/llag.Rd
pkg/yuima/src/cce_functions.c
Log:
modified llag.R, llag.Rd, cce_functions.c
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-05-13 12:10:47 UTC (rev 316)
+++ pkg/yuima/DESCRIPTION 2014-07-07 02:19:11 UTC (rev 317)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.21
-Date: 2014-05-02
+Version: 1.0.22
+Date: 2014-07-07
Depends: methods, zoo, stats4, utils, expm
Suggests: cubature, mvtnorm
Author: YUIMA Project Team
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2014-05-13 12:10:47 UTC (rev 316)
+++ pkg/yuima/NEWS 2014-07-07 02:19:11 UTC (rev 317)
@@ -16,6 +16,7 @@
2013/10/28: modify llag.R
2013/10/30: modify cce.R, cce_functions.c
2013/11/21: modify llag.R
-2013/11/22: modify cce.R, cce_functions.c
-2014/04/28: modified qmle, added carma, modified lasso
-2014/05/04: modified show method, setYuima sets the sampling from the data if sampling is missing
+2013/11/22: modify cce.R, cce_functions.c
+2014/04/28: modified qmle, added carma, modified lasso
+2014/05/04: modified show method, setYuima sets the sampling from the data if sampling is missing
+2014/07/07: modified llag.R, llag.Rd, cce_functions.c; estimated cross-correlation functions are converted to values in [-1,1]
Modified: pkg/yuima/R/llag.R
===================================================================
--- pkg/yuima/R/llag.R 2014-05-13 12:10:47 UTC (rev 316)
+++ pkg/yuima/R/llag.R 2014-07-07 02:19:11 UTC (rev 317)
@@ -1,368 +1,544 @@
-#lead-lag estimation
-
-#x:data
-
-#setGeneric( "llag", function(x,verbose=FALSE) standardGeneric("llag") )
-#setMethod( "llag", "yuima", function(x,verbose=FALSE) llag(x at data,verbose=verbose ))
-#setMethod( "llag", "yuima.data", function(x,verbose=FALSE) {
-
-#setGeneric( "llag",
-# function(x,from=FALSE,to=FALSE,division=FALSE,verbose=FALSE)
-# standardGeneric("llag") )
-#setMethod( "llag", "yuima",
-# function(x,from=FALSE,to=FALSE,division=FALSE,verbose=FALSE)
-# 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) {
-
-setGeneric( "llag",
- function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid)
- standardGeneric("llag") )
-setMethod( "llag", "yuima",
- function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid)
- llag(x at data,from=from,to=to,division=division,verbose=verbose,grid=grid))
-setMethod( "llag", "yuima.data", function(x,from=-Inf,to=Inf,division=FALSE,
- verbose=FALSE,grid) {
-
- if((is(x)=="yuima")||(is(x)=="yuima.data")){
- zdata <- get.zoo.data(x)
- }else{
- print("llag:invalid argument")
- return(NULL)
- }
-
- d <- length(zdata)
-
- # allocate memory
- ser.times <- vector(d, mode="list") # time index in 'x'
- ser.diffX <- vector(d, mode="list") # difference of data
- vol <- double(d)
-
- for(i in 1:d){
-
- # NA data must be skipped
- idt <- which(is.na(zdata[[i]]))
- if(length(idt>0)){
- zdata[[i]] <- zdata[[i]][-idt]
- }
- if(length(zdata[[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(zdata[[i]]))
- # set difference of the data
- ser.diffX[[i]] <- diff( as.numeric(zdata[[i]]) )
- vol[i] <- sum(ser.diffX[[i]]^2)
- }
-
- theta <- matrix(0,d,d)
- covmat <- diag(vol)
-
- d.size <- d*(d-1)/2
- crosscov <- vector(d.size,mode="list")
-
- if(missing(grid)){
-
- 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 <- ser.times[[i]]
- time2 <- ser.times[[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)
- }
- }
-
- y <- seq(-num2-tmptheta,num1-tmptheta,length=n)[2:(n-1)]
-
- tmp <- .C("HYcrosscov",
- as.integer(n-2),
- as.integer(length(time1)),
- as.integer(length(time2)),
- as.double(y),
- as.double(time1),
- as.double(time2),
- double(length(time2)),
- as.double(ser.diffX[[i]]),
- as.double(ser.diffX[[j]]),
- value=double(n-2),PACKAGE="yuima")$value
-
- idx <- which.max(abs(tmp))
- mlag <- -y[idx] # make the first timing of max or min
- cov <- tmp[idx]
-
- theta[i,j] <- mlag
- covmat[i,j] <- cov
- theta[j,i] <- -mlag
- covmat[j,i] <- covmat[i,j]
-
- crosscov[[num]] <- zoo(tmp,y)
- }
- }
-
- }else{
-
- if(!is.list(grid)){
- if(is.numeric(grid)){
- grid <- data.frame(matrix(grid,length(grid),d.size))
- }else{
- print("llag:invalid grid")
- return(NULL)
- }
- }
-
- 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)
-
- tmp <- .C("HYcrosscov",
- as.integer(length(grid[[num]])),
- as.integer(length(time1)),
- as.integer(length(time2)),
- as.double(grid[[num]]),
- as.double(time1),
- as.double(time2),
- double(length(time2)),
- as.double(ser.diffX[[i]]),
- as.double(ser.diffX[[j]]),
- value=double(length(grid[[num]])),PACKAGE="yuima")$value
-
- idx <- which.max(abs(tmp))
- mlag <- -grid[[num]][idx] # make the first timing of max or min
- cov <- tmp[idx]
-
- theta[i,j] <- mlag
- covmat[i,j] <- cov
- theta[j,i] <- -mlag
- covmat[j,i] <- covmat[i,j]
-
- crosscov[[num]] <- zoo(tmp,grid[[num]])
- }
- }
- }
-
- cormat <- diag(1/sqrt(diag(covmat)))%*%covmat%*%diag(1/sqrt(diag(covmat)))
- colnames(theta) <- names(zdata)
- rownames(theta) <- names(zdata)
-
- if(verbose==TRUE){
- colnames(covmat) <- names(zdata)
- rownames(covmat) <- names(zdata)
- colnames(cormat) <- names(zdata)
- rownames(cormat) <- names(zdata)
- return(list(lagcce=theta,covmat=covmat,cormat=cormat,crosscov=crosscov))
- }else{
- return(theta)
- }
-})
-
-
-## Old version
-if(0){
-setMethod( "llag", "yuima.data", function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE) {
-
-
- if(!is(x)=="yuima.data"){
- if(is(x)=="yuima"){
- dat <- x at data
- }else{
- print("llag:invalid argument")
- return(NULL)
- }
- }else{
- dat <- x
- }
-
- d <- length(dat at zoo.data)
-
- lagccep <- function(datp,theta){
- time(datp[[2]]) <- time(datp[[2]])+theta
- return(cce(setData(datp))$covmat[1,2])
- }
-
- lagcce <- function(datzoo,theta){
- d <- dim(theta)[1]
- lcmat <- cce(setData(datzoo))$covmat
-
- for(i in 1:(d-1)){
- for(j in (i+1):d){
- datp <- datzoo[c(i,j)]
- lcmat[i,j] <- lagccep(datp,theta[i,j])
- lcmat[j,i] <- lcmat[i,j]
- }
- }
- return(lcmat)
- }
-
- d.size <- d*(d-1)/2
-
- 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) != d.size){
- division <- c(division,rep(FALSE,d.size - length(division)))
- }
-
- find_lag <- function(i,j){
- datp <- dat at zoo.data[c(i,j)]
- time1 <- time(datp[[1]])
- time2 <- time(datp[[2]])
-
- # 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(datp[[1]]),length(datp[[2]])))+1
- }else{
- n <- division[num]
- }
-
- # maximum value of lagcce
-
- tmptheta <- as.numeric(time2[1])-as.numeric(time1[1]) # time lag (sec)
-
- num1 <- as.numeric(time1[length(time1)])-as.numeric(time1[1]) # elapsed time for series 1
- num2 <- as.numeric(time2[length(time2)])-as.numeric(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)
- }
- }
-
- y <- seq(-num2-tmptheta,num1-tmptheta,length=n)
- tmp <- double(n)
-
- for(i.tmp in 2:(n-1)){
- tmp[i.tmp] <- lagccep(datp,y[i.tmp])
- }
-
- mat <- cbind(y[2:(n-1)],tmp[2:(n-1)])
-
- #idx <- abs(mat[,2])==max(abs(mat[,2]))
- #mlag <- mat[,1][idx][1] # make the first timing of max or min
- #cov <- mat[,2][idx][1]
- idx <- which.max(abs(mat[,2]))
- mlag <- mat[,1][idx] # make the first timing of max or min
- cov <- mat[,2][idx]
- return(list(lag=-mlag,cov=cov))
- }
-
- theta <- matrix(numeric(d^2),ncol=d)
- covmat <- cce(dat)$covmat
-
- for(i in 1:(d-1)){
- for(j in (i+1):d){
- fl <- find_lag(i,j)
- theta[i,j] <- fl$lag
- covmat[i,j] <- fl$cov
- theta[j,i] <- -theta[i,j]
- covmat[j,i] <- covmat[i,j]
- }
- }
-
-
- # covmat <- lagcce(dat at zoo.data,theta)
- cormat <- diag(1/sqrt(diag(covmat)))%*%covmat%*%diag(1/sqrt(diag(covmat)))
- if(verbose==TRUE){
- return(list(lagcce=theta,covmat=covmat,cormat=cormat))
- }else{
- return(theta)
- }
-})
+#lead-lag estimation
+
+#x:data
+
+#setGeneric( "llag", function(x,verbose=FALSE) standardGeneric("llag") )
+#setMethod( "llag", "yuima", function(x,verbose=FALSE) llag(x at data,verbose=verbose ))
+#setMethod( "llag", "yuima.data", function(x,verbose=FALSE) {
+
+#setGeneric( "llag",
+# function(x,from=FALSE,to=FALSE,division=FALSE,verbose=FALSE)
+# standardGeneric("llag") )
+#setMethod( "llag", "yuima",
+# function(x,from=FALSE,to=FALSE,division=FALSE,verbose=FALSE)
+# 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) {
+
+setGeneric( "llag",
+ function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid,psd=TRUE,
+ plot=FALSE,ccor=FALSE)
+ standardGeneric("llag") )
+setMethod( "llag", "yuima",
+ function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid,psd=TRUE,
+ plot=FALSE,ccor=FALSE)
+ llag(x at data,from=from,to=to,division=division,verbose=verbose,grid=grid,
+ psd=psd,plot=plot))
+setMethod( "llag", "yuima.data", function(x,from=-Inf,to=Inf,division=FALSE,
+ verbose=FALSE,grid,psd=TRUE,
+ plot=FALSE,ccor=FALSE) {
+
+ if((is(x)=="yuima")||(is(x)=="yuima.data")){
+ zdata <- get.zoo.data(x)
+ }else{
+ print("llag:invalid argument")
+ return(NULL)
+ }
+
+ d <- length(zdata)
+
+ # allocate memory
+ ser.times <- vector(d, mode="list") # time index in 'x'
+ ser.diffX <- vector(d, mode="list") # difference of data
+ vol <- double(d)
+
+ for(i in 1:d){
+
+ # NA data must be skipped
+ idt <- which(is.na(zdata[[i]]))
+ if(length(idt>0)){
+ zdata[[i]] <- zdata[[i]][-idt]
+ }
+ if(length(zdata[[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(zdata[[i]]))
+ # set difference of the data
+ ser.diffX[[i]] <- diff( as.numeric(zdata[[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")
+
+ if(psd){
+
+ cormat <- diag(d)
+
+ if(missing(grid)){
+
+ 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 <- ser.times[[i]]
+ time2 <- ser.times[[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)
+ }
+ }
+
+ y <- seq(-num2-tmptheta,num1-tmptheta,length=n)[2:(n-1)]
+
+ tmp <- .C("HYcrosscorr",
+ as.integer(n-2),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(y),
+ 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(n-2))$value
+
+ idx <- which.max(abs(tmp))
+ mlag <- -y[idx] # make the first timing of max or min
+ corr <- tmp[idx]
+
+ theta[i,j] <- mlag
+ cormat[i,j] <- corr
+ theta[j,i] <- -mlag
+ cormat[j,i] <- cormat[i,j]
+
+ crosscor[[num]] <- zoo(tmp,-y)
+ }
+ }
+
+ }else{
+
+ if(!is.list(grid)){
+ if(is.numeric(grid)){
+ grid <- data.frame(matrix(grid,length(grid),d.size))
+ }else{
+ print("llag:invalid grid")
+ return(NULL)
+ }
+ }
+
+ 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)
+
+ tmp <- .C("HYcrosscorr",
+ as.integer(length(grid[[num]])),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(grid[[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(grid[[num]])))$value
+
+ idx <- which.max(abs(tmp))
+ mlag <- -grid[[num]][idx] # make the first timing of max or min
+ cor <- tmp[idx]
+
+ theta[i,j] <- mlag
+ cormat[i,j] <- cor
+ theta[j,i] <- -mlag
+ cormat[j,i] <- cormat[i,j]
+
+ crosscor[[num]] <- zoo(tmp,-grid[[num]])
+ }
+ }
+ }
+
+ covmat <- diag(sqrt(vol))%*%cormat%*%diag(sqrt(vol))
+
+ }else{# non-psd
+
+ covmat <- diag(vol)
+
+ if(missing(grid)){
+
+ 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 <- ser.times[[i]]
+ time2 <- ser.times[[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)
+ }
+ }
+
+ y <- seq(-num2-tmptheta,num1-tmptheta,length=n)[2:(n-1)]
+
+ tmp <- .C("HYcrosscov",
+ as.integer(n-2),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(y),
+ as.double(time1),
+ as.double(time2),
+ double(length(time2)),
+ as.double(ser.diffX[[i]]),
+ as.double(ser.diffX[[j]]),
+ value=double(n-2),PACKAGE="yuima")$value
+
+ idx <- which.max(abs(tmp))
+ mlag <- -y[idx] # make the first timing of max or min
+ cov <- tmp[idx]
+
+ theta[i,j] <- mlag
+ covmat[i,j] <- cov
+ theta[j,i] <- -mlag
+ covmat[j,i] <- covmat[i,j]
+
+ crosscor[[num]] <- zoo(tmp,-y)/sqrt(vol[i]*vol[j])
+ }
+ }
+
+ }else{
+
+ if(!is.list(grid)){
+ if(is.numeric(grid)){
+ grid <- data.frame(matrix(grid,length(grid),d.size))
+ }else{
+ print("llag:invalid grid")
+ return(NULL)
+ }
+ }
+
+ 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)
+
+ tmp <- .C("HYcrosscov",
+ as.integer(length(grid[[num]])),
+ as.integer(length(time1)),
+ as.integer(length(time2)),
+ as.double(grid[[num]]),
+ as.double(time1),
+ as.double(time2),
+ double(length(time2)),
+ as.double(ser.diffX[[i]]),
+ as.double(ser.diffX[[j]]),
+ value=double(length(grid[[num]])),PACKAGE="yuima")$value
+
+ idx <- which.max(abs(tmp))
+ mlag <- -grid[[num]][idx] # make the first timing of max or min
+ cov <- tmp[idx]
+
+ theta[i,j] <- mlag
+ covmat[i,j] <- cov
+ theta[j,i] <- -mlag
+ covmat[j,i] <- covmat[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)))
+ }
+
+ colnames(theta) <- names(zdata)
+ rownames(theta) <- names(zdata)
+
+ 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(zdata)
+ rownames(covmat) <- names(zdata)
+ colnames(cormat) <- names(zdata)
+ rownames(cormat) <- names(zdata)
+ if(ccor){
+ return(list(lagcce=theta,covmat=covmat,cormat=cormat,ccor=crosscor))
+ }else{
+ return(list(lagcce=theta,covmat=covmat,cormat=cormat))
+ }
+ }else{
+ return(theta)
+ }
+})
+
+
+## Old version
+if(0){
+setMethod( "llag", "yuima.data", function(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE) {
+
+
+ if(!is(x)=="yuima.data"){
+ if(is(x)=="yuima"){
+ dat <- x at data
+ }else{
+ print("llag:invalid argument")
+ return(NULL)
+ }
+ }else{
+ dat <- x
+ }
+
+ d <- length(dat at zoo.data)
+
+ lagccep <- function(datp,theta){
+ time(datp[[2]]) <- time(datp[[2]])+theta
+ return(cce(setData(datp))$covmat[1,2])
+ }
+
+ lagcce <- function(datzoo,theta){
+ d <- dim(theta)[1]
+ lcmat <- cce(setData(datzoo))$covmat
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+ datp <- datzoo[c(i,j)]
+ lcmat[i,j] <- lagccep(datp,theta[i,j])
+ lcmat[j,i] <- lcmat[i,j]
+ }
+ }
+ return(lcmat)
+ }
+
+ d.size <- d*(d-1)/2
+
+ 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) != d.size){
+ division <- c(division,rep(FALSE,d.size - length(division)))
+ }
+
+ find_lag <- function(i,j){
+ datp <- dat at zoo.data[c(i,j)]
+ time1 <- time(datp[[1]])
+ time2 <- time(datp[[2]])
+
+ # 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(datp[[1]]),length(datp[[2]])))+1
+ }else{
+ n <- division[num]
+ }
+
+ # maximum value of lagcce
+
+ tmptheta <- as.numeric(time2[1])-as.numeric(time1[1]) # time lag (sec)
+
+ num1 <- as.numeric(time1[length(time1)])-as.numeric(time1[1]) # elapsed time for series 1
+ num2 <- as.numeric(time2[length(time2)])-as.numeric(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)
+ }
+ }
+
+ y <- seq(-num2-tmptheta,num1-tmptheta,length=n)
+ tmp <- double(n)
+
+ for(i.tmp in 2:(n-1)){
+ tmp[i.tmp] <- lagccep(datp,y[i.tmp])
+ }
+
+ mat <- cbind(y[2:(n-1)],tmp[2:(n-1)])
+
+ #idx <- abs(mat[,2])==max(abs(mat[,2]))
+ #mlag <- mat[,1][idx][1] # make the first timing of max or min
+ #cov <- mat[,2][idx][1]
+ idx <- which.max(abs(mat[,2]))
+ mlag <- mat[,1][idx] # make the first timing of max or min
+ cov <- mat[,2][idx]
+ return(list(lag=-mlag,cov=cov))
+ }
+
+ theta <- matrix(numeric(d^2),ncol=d)
+ covmat <- cce(dat)$covmat
+
+ for(i in 1:(d-1)){
+ for(j in (i+1):d){
+ fl <- find_lag(i,j)
+ theta[i,j] <- fl$lag
+ covmat[i,j] <- fl$cov
+ theta[j,i] <- -theta[i,j]
+ covmat[j,i] <- covmat[i,j]
+ }
+ }
+
+
+ # covmat <- lagcce(dat at zoo.data,theta)
+ cormat <- diag(1/sqrt(diag(covmat)))%*%covmat%*%diag(1/sqrt(diag(covmat)))
+ if(verbose==TRUE){
+ return(list(lagcce=theta,covmat=covmat,cormat=cormat))
+ }else{
+ return(theta)
+ }
+})
}
\ No newline at end of file
Modified: pkg/yuima/man/llag.Rd
===================================================================
--- pkg/yuima/man/llag.Rd 2014-05-13 12:10:47 UTC (rev 316)
+++ pkg/yuima/man/llag.Rd 2014-07-07 02:19:11 UTC (rev 317)
@@ -5,7 +5,8 @@
\description{Estimate the lead-lag parameters of discretely observed processes by maximizing the shifted Hayashi-Yoshida covariation contrast functions, following Hoffmann et al. (2013).
}
\usage{
-llag(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid)
+llag(x,from=-Inf,to=Inf,division=FALSE,verbose=FALSE,grid,psd=TRUE,
+ plot=FALSE,ccor=FALSE)
}
\arguments{
\item{x}{an object of \code{\link{yuima-class}} or
@@ -15,25 +16,33 @@
\item{to}{a numeric vector each of whose component(s) indicates the upper end of a finite grid on which the contrast function is evaluated, if \code{grid} is missing.}
\item{division}{a numeric vector each of whose component(s) indicates the number of the points of a finite grid on which the contrast function is evaluated, if \code{grid} is missing.}
\item{grid}{a numeric vector or a list of numeric vectors. See `Details'.}
+\item{psd}{logical. If \code{TRUE}, the estimated cross-correlation functions are converted to the interval [-1,1]. See `Details'.}
+\item{plot}{logical. If \code{TRUE}, the estimated cross-correlation functions are plotted.}
+\item{ccor}{logical. If \code{TRUE}, the estimated cross-correlation functions are returned. This argument is ignored if \code{verbose} is \code{FALSE}.}
}
\details{
Let \eqn{d} be the number of the components of the \code{zoo.data} of the object \code{x}.
Let \eqn{Xi_{ti_0},Xi_{ti_1},\dots,Xi_{ti_n(i)}} be the observation data of the \eqn{i}-th component (i.e. the \eqn{i}-th component of the \code{zoo.data} of the object \code{x}).
-The shifted Hayashi-Yoshida covariation contrast function \eqn{Uij(\theta)} of the observations \eqn{Xi} and \eqn{Xj} \eqn{(i<j)} is defined by the same way as in Hoffmann et al. (2013). The lead-lag parameter \eqn{\theta_ij} is defined as a maximizer of \eqn{|Uij(\theta)|}. \eqn{Uij(\theta)} is evaluated on a finite grid \eqn{Gij} defined below. Thus \eqn{\theta_ij} belongs to this grid. If there exist more than two maximizers, the lowest one is selected.
+The shifted Hayashi-Yoshida covariation contrast function \eqn{Uij(\theta)} of the observations \eqn{Xi} and \eqn{Xj} \eqn{(i<j)} is defined by the same way as in Hoffmann et al. (2013), which corresponds to their cross-covariance function. The lead-lag parameter \eqn{\theta_ij} is defined as a maximizer of \eqn{|Uij(\theta)|}. \eqn{Uij(\theta)} is evaluated on a finite grid \eqn{Gij} defined below. Thus \eqn{\theta_ij} belongs to this grid. If there exist more than two maximizers, the lowest one is selected.
+If \code{psd} is \code{TRUE}, For any \eqn{i,j} the matrix \eqn{C:=(Ukl(\theta))_{k,l\in{i,j}}} is converted to \code{(C\%*\%C)^(1/2)} for ensuring the positive semi-definiteness, and \eqn{Uij(\theta)} is redefined as the \eqn{(1,2)}-component of the converted \eqn{C}. Here, \eqn{Ukk(\theta)} is set to the realized volatility of \eqn{Xk}. In this case \eqn{\theta_ij} is given as a maximizer of the cross-correlation functions.
+
The grid \eqn{Gij} is defined as follows. First, if \code{grid} is missing, \eqn{Gij} is given by
\deqn{a, a+(b-a)/(N-1), \dots, a+(N-2)(b-a)/(N-1), b,}
where \eqn{a,b} and \eqn{N} are the \eqn{(d(i-1)-(i-1)i/2+(j-i))}-th components of \code{from}, \code{to} and \code{division} respectively. If the corresponding component of \code{from} (resp. \code{to}) is \code{-Inf} (resp. \code{Inf}), \eqn{a=-(tj_n(j)-ti_0)} (resp. \eqn{b=ti_n(i)-tj_0}) is used, while if the corresponding component of \code{division} is \code{FALSE}, \eqn{N=round(2max(n(i),n(j)))+1} is used. Missing components are filled with \code{-Inf} (resp. \code{Inf}, \code{FALSE}). The default value \code{-Inf} (resp. \code{Inf}, \code{FALSE}) means that all components are \code{-Inf} (resp. \code{Inf}, \code{FALSE}). Next, if \code{grid} is a numeric vector, \eqn{Gij} is given by \code{grid}. If \code{grid} is a list of numeric vectors, \eqn{Gij} is given by the \eqn{(d(i-1)-(i-1)i/2+(j-i))}-th component of \code{grid}.
-The estimated lead-lag parameters are returned as the skew-symmetric matrix \eqn{(\theta_ij)_{i,j=1,\dots,d}}. If \code{verbose} is \code{TRUE}, the covariance matrix \eqn{(Uij(\theta_ij))_{i,j=1,\dots,d}} corresponding to the estimated lead-lag parameters, the corresponding correlation matrix and the computed contrast functions are also returned. The computed contrast functions are returned as a list with the length \eqn{d(d-1)/2}. For \eqn{i<j}, the \eqn{(d(i-1)-(i-1)i/2+(j-i))}-th component of the list consists of an object \eqn{Uij(\theta)} of class \link{zoo} indexed by \eqn{Gij}.}
+The estimated lead-lag parameters are returned as the skew-symmetric matrix \eqn{(\theta_ij)_{i,j=1,\dots,d}}. If \code{verbose} is \code{TRUE}, the covariance matrix \eqn{(Uij(\theta_ij))_{i,j=1,\dots,d}} corresponding to the estimated lead-lag parameters, the corresponding correlation matrix and the computed contrast functions are also returned. If further \code{ccor} is \code{TRUE},the computed cross-correlation functions are returned as a list with the length \eqn{d(d-1)/2}. For \eqn{i<j}, the \eqn{(d(i-1)-(i-1)i/2+(j-i))}-th component of the list consists of an object \eqn{Uij(\theta)/sqrt(Uii(\theta)*Ujj(\theta))} of class \link{zoo} indexed by \eqn{Gij}.
+
+If \code{plot} is \code{TRUE}, the computed cross-correlation functions are plotted sequentially.}
\value{
If \code{verbose} is \code{FALSE}, a skew-symmetric matrix corresponding to the estimated lead-lag parameters is returned. If \code{verbose} is \code{TRUE}, a list with the following components is returned:
\item{lagcce}{a skew-symmetric matrix corresponding to the estimated lead-lag parameters.}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 317
More information about the Yuima-commits
mailing list