version 1.15.3

cran · Jan 31, 2022 · f29d5ed · f29d5ed
1 parent ba5cc03
commit f29d5ed
Show file tree

Hide file tree

Showing 8 changed files with 334 additions and 152 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,18 +1,18 @@
 Package: yuima
 Type: Package
 Title: The YUIMA Project Package for SDEs
-Version: 1.15.2
+Version: 1.15.3
 Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature,
         mvtnorm
-Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast, wavethresh,
-        coda, calculus (>= 0.2.0)
+Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast, coda, calculus
+        (>= 0.2.0)
 Author: YUIMA Project Team
 Maintainer: Stefano M. Iacus <stefano.iacus@unimi.it>
 Description: Simulation and Inference for SDEs and Other Stochastic Processes.
 License: GPL-2
 URL: https://yuimaproject.com
 LinkingTo: Rcpp, RcppArmadillo
 NeedsCompilation: yes
-Packaged: 2022-01-26 09:31:31 UTC; jago
+Packaged: 2022-01-31 07:58:59 UTC; jago
 Repository: CRAN
-Date/Publication: 2022-01-27 10:20:02 UTC
+Date/Publication: 2022-01-31 08:50:02 UTC
diff --git a/MD5 b/MD5
@@ -1,7 +1,7 @@
 94d55d512a9ba36caa9b7df079bae19f *COPYING
-15c2028568528c2e4390a390fc6013b4 *DESCRIPTION
-1976c100cd4dedaaf26e90e609235545 *NAMESPACE
-050fc01ddf4030dbcdafabb7f86cd760 *NEWS
+3bf57fa4c251682f9e0272e77a34dbab *DESCRIPTION
+e14bf414f7bf1023bc0fec2d9e443654 *NAMESPACE
+6c3baa4587c5855c80559757255c78ac *NEWS
 5142e55d1bab0b3d632328024ea4581b *R/AllClasses.R
 a57282007ddeb97e5a975ab60fe1d12b *R/AuxMethodforPPR.R
 4b22fa1f5201199bc7d2afc39b2f41e7 *R/CPoint.R
@@ -72,7 +72,7 @@ d0600e3fdc00d5a0c297316425f6c814 *R/simulateMultiProcess.R
 76ef13ea18f2507ece6c661cf3d1b805 *R/spectralcov.R
 f506cf2d30a5f912b5221112d7ae8965 *R/subsampling.R
 3009a7d539f3dbbef7b828129687062c *R/toLatex.R
-1db85d2f66897c8375a9539707f5b4e7 *R/wllag.r
+76530303d473d1152148193fc5a52c4e *R/wllag.R
 1ee5ca21b65c9c4bde9656529686669d *R/yuima.R
 edfabacecae0ebb9b40324f1d499b4b9 *R/yuima.characteristic.R
 6a4223ef8dae8e3b128da32d29e61b20 *R/yuima.data.R
@@ -81,7 +81,7 @@ edfabacecae0ebb9b40324f1d499b4b9 *R/yuima.characteristic.R
 28a257804b7d08a3fb47100a5c536234 *R/yuima.model.R
 e3c164669a627b3f16fc6da48e714da5 *R/yuima.sampling.R
 48a7ba9f63b106e188f6c3dc19d836eb *R/zzz.R
-0302a4260d695fbfc6b4440a8191acb3 *build/partial.rdb
+6d72f8f507ebad8777406d19040f9ec4 *build/partial.rdb
 b516711f126a8aafd4aae9b758da7824 *data/LogSPX.rda
 bce7107dc324a66a94865b56f9e3bbf3 *data/MWK151.rda
 ba391af8d40ef453c6d0a726d4170ffa *inst/AUTHORS
@@ -177,7 +177,7 @@ ab5b290fd33ce465c9c035d49995ce77 *man/simFunctional.Rd
 114989e6fbd3d202b0f3a6650a019abb *man/subsampling.Rd
 6536952261748e9437d0f7fb8b4ee476 *man/toLatex.Rd
 c8383c3a57f9c6f87e8d321e4db0e92e *man/variable.Integral.rd
-714b6309e166da646bfae28372358fff *man/wllag.Rd
+84e4060d13d434ffa5cab09ba38d444a *man/wllag.Rd
 d98511434eb2f68e90622e9d1bb6b8de *man/ybook.Rd
 fef5e35bf94cc6944496627e9f04a084 *man/yuima-class.Rd
 ab1e99f44258125e416cee06598f2ca3 *man/yuima.CP.qmle-class.Rd

diff --git a/NAMESPACE b/NAMESPACE
@@ -65,6 +65,7 @@ importFrom(stats, sd)
 importFrom("stats", cov2cor) # added by YK on Apr. 12, 2018
 importFrom("graphics", "abline", "hist", "par", "text") # added by YK on Jul. 19, 2019
 importFrom("stats", "mvfft") # added by YK on Apr. 15, 2020
+importFrom("stats", "convolve") # added by YK on Jan. 28, 2022
 importFrom(coda, as.mcmc)
 importFrom("stats", "cov", "var")
 

diff --git a/NEWS b/NEWS
@@ -94,4 +94,5 @@
             aeMean.Rd, aeMoment.Rd, aeSd.Rd, aeSkewness.Rd, 
             yuima.ae-class.Rd, updated NAMESPACE,  updated yuima_init.c,
 	    updated DESCRIPTION to add package calculus in Imports
-2021/01/25: fixed a protection issue in euler.c
+2022/01/25: fixed a protection issue in euler.c
+2022/01/28: modified wllag.R and unimported the wavethresh package
diff --git a/R/wllag.R b/R/wllag.R
@@ -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")
+  }
+
+}