[Yuima-commits] r785 - pkg/yuima/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Jan 28 03:44:00 CET 2022


Author: kyuta
Date: 2022-01-28 03:43:59 +0100 (Fri, 28 Jan 2022)
New Revision: 785

Added:
   pkg/yuima/R/wllag.R
Removed:
   pkg/yuima/R/wllag.R
Log:
replaced wllag.R to fix a svn issue

Deleted: pkg/yuima/R/wllag.R
===================================================================
--- pkg/yuima/R/wllag.R	2022-01-28 02:06:02 UTC (rev 784)
+++ pkg/yuima/R/wllag.R	2022-01-28 02:43:59 UTC (rev 785)
@@ -1,318 +0,0 @@
-# Scale-by-scale lead-lag estimation by wavelets
-
-## function to compute Daubechies' extremal phase wavelet filter
-## The function is based on implementation of wavethresh package
-daubechies.wavelet <- function(N){
-  
-  switch (N,
-    "1" = {
-      H <- rep(0, 2)
-      H[1] <- 1/sqrt(2)
-      H[2] <- H[1]
-    },
-    "2" = {
-      H <- rep(0, 4)
-      H[1] <- 0.482962913145
-      H[2] <- 0.836516303738
-      H[3] <- 0.224143868042
-      H[4] <- -0.129409522551
-    },
-    "3" = {
-      H <- rep(0, 6)
-      H[1] <- 0.33267055295
-      H[2] <- 0.806891509311
-      H[3] <- 0.459877502118
-      H[4] <- -0.13501102001
-      H[5] <- -0.085441273882
-      H[6] <- 0.035226291882
-    },
-    "4" = {
-      H <- rep(0, 8)
-      H[1] <- 0.230377813309
-      H[2] <- 0.714846570553
-      H[3] <- 0.63088076793
-      H[4] <- -0.027983769417
-      H[5] <- -0.187034811719
-      H[6] <- 0.030841381836
-      H[7] <- 0.032883011667
-      H[8] <- -0.010597401785
-    },
-    "5" = {
-      H <- rep(0, 10)
-      H[1] <- 0.160102397974
-      H[2] <- 0.603829269797
-      H[3] <- 0.724308528438
-      H[4] <- 0.138428145901
-      H[5] <- -0.242294887066
-      H[6] <- -0.032244869585
-      H[7] <- 0.07757149384
-      H[8] <- -0.006241490213
-      H[9] <- -0.012580752
-      H[10] <- 0.003335725285
-    },
-    "6" = {
-      H <- rep(0, 12)
-      H[1] <- 0.11154074335
-      H[2] <- 0.494623890398
-      H[3] <- 0.751133908021
-      H[4] <- 0.315250351709
-      H[5] <- -0.226264693965
-      H[6] <- -0.129766867567
-      H[7] <- 0.097501605587
-      H[8] <- 0.02752286553
-      H[9] <- -0.031582039318
-      H[10] <- 0.000553842201
-      H[11] <- 0.004777257511
-      H[12] <- -0.001077301085
-    },
-    "7" = {
-      H <- rep(0, 14)
-      H[1] <- 0.077852054085
-      H[2] <- 0.396539319482
-      H[3] <- 0.729132090846
-      H[4] <- 0.469782287405
-      H[5] <- -0.143906003929
-      H[6] <- -0.224036184994
-      H[7] <- 0.071309219267
-      H[8] <- 0.080612609151
-      H[9] <- -0.038029936935
-      H[10] <- -0.016574541631
-      H[11] <- 0.012550998556
-      H[12] <- 0.000429577973
-      H[13] <- -0.001801640704
-      H[14] <- 0.0003537138
-    },
-    "8" = {
-      H <- rep(0, 16)
-      H[1] <- 0.054415842243
-      H[2] <- 0.312871590914
-      H[3] <- 0.675630736297
-      H[4] <- 0.585354683654
-      H[5] <- -0.015829105256
-      H[6] <- -0.284015542962
-      H[7] <- 0.000472484574
-      H[8] <- 0.12874742662
-      H[9] <- -0.017369301002
-      H[10] <- -0.044088253931
-      H[11] <- 0.013981027917
-      H[12] <- 0.008746094047
-      H[13] <- -0.004870352993
-      H[14] <- -0.000391740373
-      H[15] <- 0.000675449406
-      H[16] <- -0.000117476784
-    },
-    "9" = {
-      H <- rep(0, 18)
-      H[1] <- 0.038077947364
-      H[2] <- 0.243834674613
-      H[3] <- 0.60482312369
-      H[4] <- 0.657288078051
-      H[5] <- 0.133197385825
-      H[6] <- -0.293273783279
-      H[7] <- -0.096840783223
-      H[8] <- 0.148540749338
-      H[9] <- 0.030725681479
-      H[10] <- -0.067632829061
-      H[11] <- 0.000250947115
-      H[12] <- 0.022361662124
-      H[13] <- -0.004723204758
-      H[14] <- -0.004281503682
-      H[15] <- 0.001847646883
-      H[16] <- 0.000230385764
-      H[17] <- -0.000251963189
-      H[18] <- 3.934732e-05
-    },
-    "10" = {
-      H <- rep(0, 20)
-      H[1] <- 0.026670057901
-      H[2] <- 0.188176800078
-      H[3] <- 0.527201188932
-      H[4] <- 0.688459039454
-      H[5] <- 0.281172343661
-      H[6] <- -0.249846424327
-      H[7] <- -0.195946274377
-      H[8] <- 0.127369340336
-      H[9] <- 0.093057364604
-      H[10] <- -0.071394147166
-      H[11] <- -0.029457536822
-      H[12] <- 0.033212674059
-      H[13] <- 0.003606553567
-      H[14] <- -0.010733175483
-      H[15] <- 0.001395351747
-      H[16] <- 0.001992405295
-      H[17] <- -0.000685856695
-      H[18] <- -0.000116466855
-      H[19] <- 9.358867e-05
-      H[20] <- -1.3264203e-05
-    }
-  )
-  
-  return(H)
-}
-
-## function to compute autocorrelation wavelets
-autocorrelation.wavelet <- function(J, N){
-  
-  h <- daubechies.wavelet(N)
-  Phi1 <- convolve(h, h, type = "o")
-  
-  out <- vector(mode = "list", J)
-  out[[1]] <- (-1)^((-2*N+1):(2*N-1)) * Phi1
-  
-  if(J > 1){
-    
-    Phi1 <- Phi1[seq(1,4*N-1,by=2)]
-    
-    for(j in 2:J){
-      
-      Lj <- (2^j - 1) * (2 * N - 1) + 1
-      out[[j]] <- double(2*Lj - 1)
-      out[[j]][seq(2*N,2*Lj-2*N,by=2)] <- out[[j-1]]
-      out[[j]][seq(1,2*Lj-1,by=2)] <- 
-        convolve(out[[j-1]], Phi1, type = "o")
-      
-    }
-  }
-  
-  return(out)
-}
-
-
-## main function
-wllag <- function(x, y, J = 8, N = 10, #family = "DaubExPhase", 
-                  tau = 1e-3, from = -to, to = 100, 
-                  verbose = FALSE, in.tau = FALSE, tol = 1e-6){
-  
-  time1 <- as.numeric(time(x))
-  time2 <- as.numeric(time(y))
-  
-  grid <- seq(from, to, by = 1) * tau
-  
-  Lj <- (2^J - 1) * (2 * N - 1) + 1
-  #Lj <- (2^J - 1) * (length(wavethresh::filter.select(N, family)$H) - 1) + 1
-  grid2 <- seq(from - Lj + 1, to + Lj - 1, by = 1) * tau
-  
-  dx <- diff(as.numeric(x))
-  dy <- diff(as.numeric(y))
-  
-  tmp <- .C("HYcrosscov2",
-            as.integer(length(grid2)),
-            as.integer(length(time2)),
-            as.integer(length(time1)),
-            as.double(grid2/tol),
-            as.double(time2/tol),
-            as.double(time1/tol),
-            as.double(dy),
-            as.double(dx),
-            value=double(length(grid2)),
-            PACKAGE = "yuima")$value
-  
-  #if(missing(J)) J <- floor(log2(length(grid)))
-  
-  #if(J < 2) stop("J must be larger than 1")
-  
-  #acw <- wavethresh::PsiJ(-J, filter.number = N, family = "DaubExPhase")
-  #acw <- wavethresh::PsiJ(-J, filter.number = N, family = family)
-  acw <- autocorrelation.wavelet(J, N)
-  
-  theta <- double(J)
-  #covar <- double(J)
-  #LLR <- double(J)
-  corr <- double(J)
-  crosscor <- vector("list", J)
-  
-  for(j in 1:J){
-    
-    wcov <- try(stats::filter(tmp, filter = acw[[j]], method = "c", 
-                              sides = 2)[Lj:(length(grid) + Lj - 1)],
-                silent = TRUE)
-    #Mj <- (2^J - 2^j) * (2 * N - 1)
-    #wcov <- try(convolve(tmp, acw[[j]], conj = FALSE, type = "filter")[(Mj + 1):(length(grid) + Mj)],
-    #            silent = TRUE)
-    
-    if(class(wcov) == "try-error"){
-      
-      theta[j] <- NA
-      crosscor[[j]] <- NA
-      corr[j] <- NA
-      
-    }else{
-      
-      #tmp.grid <- grid[-attr(wcov, "na.action")]
-      crosscor[[j]] <- zoo(wcov, grid)
-      
-      obj <- abs(wcov)
-      idx1 <- which(obj == max(obj, na.rm = TRUE))
-      idx <- idx1[which.max(abs(grid[idx1]))]
-      # if there are multiple peaks, take the lag farthest from zero
-      theta[j] <- grid[idx]
-      corr[j] <- crosscor[[j]][idx]
-      
-    }
-    
-  }
-  
-  if(verbose == TRUE){
-    
-    #obj0 <- tmp[(Lj + 1):(length(grid) + Lj)]
-    obj0 <- tmp[Lj:(length(grid) + Lj - 1)]/sqrt(sum(dx^2)*sum(dy^2))
-    obj <- abs(obj0)
-    idx1 <- which(obj == max(obj, na.rm = TRUE))
-    idx <- idx1[which.max(abs(grid[idx1]))]
-    # if there are multiple peaks, the lag farthest from zero
-    theta.hy <- grid[idx]
-    corr.hy <- obj0[idx]
-    
-    if(in.tau == TRUE){
-      theta <- round(theta/tau)
-      theta.hy <- round(theta.hy/tau)
-    }
-    
-    result <- list(lagtheta = theta, obj.values = corr, 
-                   obj.fun = crosscor, theta.hry = theta.hy, 
-                   cor.hry = corr.hy, ccor.hry = zoo(obj0, grid))
-    
-    class(result) <- "yuima.wllag"
-    
-  }else{
-    if(in.tau == TRUE){
-      result <- round(theta/tau)
-    }else{
-      result <- theta
-    }
-  }
-  
-  return(result)
-}
-
-# print method for yuima.wllag-class
-print.yuima.wllag <- function(x, ...){
-  
-  cat("Estimated scale-by-scale lead-lag parameters\n")
-  print(x$lagtheta, ...)
-  cat("Corresponding values of objective functions\n")
-  print(x$obj.values, ...)
-  cat("Estimated lead-lag parameter in the HRY sense\n")
-  print(x$theta.hry, ...)
-  cat("Corresponding correlation coefficient\n")
-  print(x$cor.hry, ...)
-  
-}
-
-# plot method for yuima.wllag class
-plot.yuima.wllag <- function(x, selectJ = NULL, xlab = expression(theta), 
-                             ylab = "", ...){
-  
-  J <- length(x$lagtheta)
-  
-  if(is.null(selectJ)) selectJ <- 1:J
-  
-  for(j in selectJ){
-    #plot(x$ccor[[j]], main=paste("j=",j), xlab=expression(theta),
-    #     ylab=expression(U[j](theta)), type = type, pch = pch, ...)
-    plot(x$obj.fun[[j]], main = paste("j = ", j, sep =""), 
-         xlab = xlab, ylab = ylab, ...)
-    abline(0, 0, lty = "dotted")
-  }
-  
-}

Added: pkg/yuima/R/wllag.R
===================================================================
--- pkg/yuima/R/wllag.R	                        (rev 0)
+++ pkg/yuima/R/wllag.R	2022-01-28 02:43:59 UTC (rev 785)
@@ -0,0 +1,318 @@
+# Scale-by-scale lead-lag estimation by wavelets
+
+## function to compute Daubechies' extremal phase wavelet filter
+## The function is based on implementation of wavethresh package
+daubechies.wavelet <- function(N){
+  
+  switch (N,
+    "1" = {
+      H <- rep(0, 2)
+      H[1] <- 1/sqrt(2)
+      H[2] <- H[1]
+    },
+    "2" = {
+      H <- rep(0, 4)
+      H[1] <- 0.482962913145
+      H[2] <- 0.836516303738
+      H[3] <- 0.224143868042
+      H[4] <- -0.129409522551
+    },
+    "3" = {
+      H <- rep(0, 6)
+      H[1] <- 0.33267055295
+      H[2] <- 0.806891509311
+      H[3] <- 0.459877502118
+      H[4] <- -0.13501102001
+      H[5] <- -0.085441273882
+      H[6] <- 0.035226291882
+    },
+    "4" = {
+      H <- rep(0, 8)
+      H[1] <- 0.230377813309
+      H[2] <- 0.714846570553
+      H[3] <- 0.63088076793
+      H[4] <- -0.027983769417
+      H[5] <- -0.187034811719
+      H[6] <- 0.030841381836
+      H[7] <- 0.032883011667
+      H[8] <- -0.010597401785
+    },
+    "5" = {
+      H <- rep(0, 10)
+      H[1] <- 0.160102397974
+      H[2] <- 0.603829269797
+      H[3] <- 0.724308528438
+      H[4] <- 0.138428145901
+      H[5] <- -0.242294887066
+      H[6] <- -0.032244869585
+      H[7] <- 0.07757149384
+      H[8] <- -0.006241490213
+      H[9] <- -0.012580752
+      H[10] <- 0.003335725285
+    },
+    "6" = {
+      H <- rep(0, 12)
+      H[1] <- 0.11154074335
+      H[2] <- 0.494623890398
+      H[3] <- 0.751133908021
+      H[4] <- 0.315250351709
+      H[5] <- -0.226264693965
+      H[6] <- -0.129766867567
+      H[7] <- 0.097501605587
+      H[8] <- 0.02752286553
+      H[9] <- -0.031582039318
+      H[10] <- 0.000553842201
+      H[11] <- 0.004777257511
+      H[12] <- -0.001077301085
+    },
+    "7" = {
+      H <- rep(0, 14)
+      H[1] <- 0.077852054085
+      H[2] <- 0.396539319482
+      H[3] <- 0.729132090846
+      H[4] <- 0.469782287405
+      H[5] <- -0.143906003929
+      H[6] <- -0.224036184994
+      H[7] <- 0.071309219267
+      H[8] <- 0.080612609151
+      H[9] <- -0.038029936935
+      H[10] <- -0.016574541631
+      H[11] <- 0.012550998556
+      H[12] <- 0.000429577973
+      H[13] <- -0.001801640704
+      H[14] <- 0.0003537138
+    },
+    "8" = {
+      H <- rep(0, 16)
+      H[1] <- 0.054415842243
+      H[2] <- 0.312871590914
+      H[3] <- 0.675630736297
+      H[4] <- 0.585354683654
+      H[5] <- -0.015829105256
+      H[6] <- -0.284015542962
+      H[7] <- 0.000472484574
+      H[8] <- 0.12874742662
+      H[9] <- -0.017369301002
+      H[10] <- -0.044088253931
+      H[11] <- 0.013981027917
+      H[12] <- 0.008746094047
+      H[13] <- -0.004870352993
+      H[14] <- -0.000391740373
+      H[15] <- 0.000675449406
+      H[16] <- -0.000117476784
+    },
+    "9" = {
+      H <- rep(0, 18)
+      H[1] <- 0.038077947364
+      H[2] <- 0.243834674613
+      H[3] <- 0.60482312369
+      H[4] <- 0.657288078051
+      H[5] <- 0.133197385825
+      H[6] <- -0.293273783279
+      H[7] <- -0.096840783223
+      H[8] <- 0.148540749338
+      H[9] <- 0.030725681479
+      H[10] <- -0.067632829061
+      H[11] <- 0.000250947115
+      H[12] <- 0.022361662124
+      H[13] <- -0.004723204758
+      H[14] <- -0.004281503682
+      H[15] <- 0.001847646883
+      H[16] <- 0.000230385764
+      H[17] <- -0.000251963189
+      H[18] <- 3.934732e-05
+    },
+    "10" = {
+      H <- rep(0, 20)
+      H[1] <- 0.026670057901
+      H[2] <- 0.188176800078
+      H[3] <- 0.527201188932
+      H[4] <- 0.688459039454
+      H[5] <- 0.281172343661
+      H[6] <- -0.249846424327
+      H[7] <- -0.195946274377
+      H[8] <- 0.127369340336
+      H[9] <- 0.093057364604
+      H[10] <- -0.071394147166
+      H[11] <- -0.029457536822
+      H[12] <- 0.033212674059
+      H[13] <- 0.003606553567
+      H[14] <- -0.010733175483
+      H[15] <- 0.001395351747
+      H[16] <- 0.001992405295
+      H[17] <- -0.000685856695
+      H[18] <- -0.000116466855
+      H[19] <- 9.358867e-05
+      H[20] <- -1.3264203e-05
+    }
+  )
+  
+  return(H)
+}
+
+## function to compute autocorrelation wavelets
+autocorrelation.wavelet <- function(J, N){
+  
+  h <- daubechies.wavelet(N)
+  Phi1 <- convolve(h, h, type = "o")
+  
+  out <- vector(mode = "list", J)
+  out[[1]] <- (-1)^((-2*N+1):(2*N-1)) * Phi1
+  
+  if(J > 1){
+    
+    Phi1 <- Phi1[seq(1,4*N-1,by=2)]
+    
+    for(j in 2:J){
+      
+      Lj <- (2^j - 1) * (2 * N - 1) + 1
+      out[[j]] <- double(2*Lj - 1)
+      out[[j]][seq(2*N,2*Lj-2*N,by=2)] <- out[[j-1]]
+      out[[j]][seq(1,2*Lj-1,by=2)] <- 
+        convolve(out[[j-1]], Phi1, type = "o")
+      
+    }
+  }
+  
+  return(out)
+}
+
+
+## main function
+wllag <- function(x, y, J = 8, N = 10, #family = "DaubExPhase", 
+                  tau = 1e-3, from = -to, to = 100, 
+                  verbose = FALSE, in.tau = FALSE, tol = 1e-6){
+  
+  time1 <- as.numeric(time(x))
+  time2 <- as.numeric(time(y))
+  
+  grid <- seq(from, to, by = 1) * tau
+  
+  Lj <- (2^J - 1) * (2 * N - 1) + 1
+  #Lj <- (2^J - 1) * (length(wavethresh::filter.select(N, family)$H) - 1) + 1
+  grid2 <- seq(from - Lj + 1, to + Lj - 1, by = 1) * tau
+  
+  dx <- diff(as.numeric(x))
+  dy <- diff(as.numeric(y))
+  
+  tmp <- .C("HYcrosscov2",
+            as.integer(length(grid2)),
+            as.integer(length(time2)),
+            as.integer(length(time1)),
+            as.double(grid2/tol),
+            as.double(time2/tol),
+            as.double(time1/tol),
+            as.double(dy),
+            as.double(dx),
+            value=double(length(grid2)),
+            PACKAGE = "yuima")$value
+  
+  #if(missing(J)) J <- floor(log2(length(grid)))
+  
+  #if(J < 2) stop("J must be larger than 1")
+  
+  #acw <- wavethresh::PsiJ(-J, filter.number = N, family = "DaubExPhase")
+  #acw <- wavethresh::PsiJ(-J, filter.number = N, family = family)
+  acw <- autocorrelation.wavelet(J, N)
+  
+  theta <- double(J)
+  #covar <- double(J)
+  #LLR <- double(J)
+  corr <- double(J)
+  crosscor <- vector("list", J)
+  
+  for(j in 1:J){
+    
+    wcov <- try(stats::filter(tmp, filter = acw[[j]], method = "c", 
+                              sides = 2)[Lj:(length(grid) + Lj - 1)],
+                silent = TRUE)
+    #Mj <- (2^J - 2^j) * (2 * N - 1)
+    #wcov <- try(convolve(tmp, acw[[j]], conj = FALSE, type = "filter")[(Mj + 1):(length(grid) + Mj)],
+    #            silent = TRUE)
+    
+    if(class(wcov) == "try-error"){
+      
+      theta[j] <- NA
+      crosscor[[j]] <- NA
+      corr[j] <- NA
+      
+    }else{
+      
+      #tmp.grid <- grid[-attr(wcov, "na.action")]
+      crosscor[[j]] <- zoo(wcov, grid)
+      
+      obj <- abs(wcov)
+      idx1 <- which(obj == max(obj, na.rm = TRUE))
+      idx <- idx1[which.max(abs(grid[idx1]))]
+      # if there are multiple peaks, take the lag farthest from zero
+      theta[j] <- grid[idx]
+      corr[j] <- crosscor[[j]][idx]
+      
+    }
+    
+  }
+  
+  if(verbose == TRUE){
+    
+    #obj0 <- tmp[(Lj + 1):(length(grid) + Lj)]
+    obj0 <- tmp[Lj:(length(grid) + Lj - 1)]/sqrt(sum(dx^2)*sum(dy^2))
+    obj <- abs(obj0)
+    idx1 <- which(obj == max(obj, na.rm = TRUE))
+    idx <- idx1[which.max(abs(grid[idx1]))]
+    # if there are multiple peaks, the lag farthest from zero
+    theta.hy <- grid[idx]
+    corr.hy <- obj0[idx]
+    
+    if(in.tau == TRUE){
+      theta <- round(theta/tau)
+      theta.hy <- round(theta.hy/tau)
+    }
+    
+    result <- list(lagtheta = theta, obj.values = corr, 
+                   obj.fun = crosscor, theta.hry = theta.hy, 
+                   cor.hry = corr.hy, ccor.hry = zoo(obj0, grid))
+    
+    class(result) <- "yuima.wllag"
+    
+  }else{
+    if(in.tau == TRUE){
+      result <- round(theta/tau)
+    }else{
+      result <- theta
+    }
+  }
+  
+  return(result)
+}
+
+# print method for yuima.wllag-class
+print.yuima.wllag <- function(x, ...){
+  
+  cat("Estimated scale-by-scale lead-lag parameters\n")
+  print(x$lagtheta, ...)
+  cat("Corresponding values of objective functions\n")
+  print(x$obj.values, ...)
+  cat("Estimated lead-lag parameter in the HRY sense\n")
+  print(x$theta.hry, ...)
+  cat("Corresponding correlation coefficient\n")
+  print(x$cor.hry, ...)
+  
+}
+
+# plot method for yuima.wllag class
+plot.yuima.wllag <- function(x, selectJ = NULL, xlab = expression(theta), 
+                             ylab = "", ...){
+  
+  J <- length(x$lagtheta)
+  
+  if(is.null(selectJ)) selectJ <- 1:J
+  
+  for(j in selectJ){
+    #plot(x$ccor[[j]], main=paste("j=",j), xlab=expression(theta),
+    #     ylab=expression(U[j](theta)), type = type, pch = pch, ...)
+    plot(x$obj.fun[[j]], main = paste("j = ", j, sep =""), 
+         xlab = xlab, ylab = ylab, ...)
+    abline(0, 0, lty = "dotted")
+  }
+  
+}



More information about the Yuima-commits mailing list