[Yuima-commits] r736 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 29 07:19:08 CEST 2020
Author: kyuta
Date: 2020-04-29 07:19:07 +0200 (Wed, 29 Apr 2020)
New Revision: 736
Added:
pkg/yuima/R/simBmllag.r
pkg/yuima/R/wllag.r
pkg/yuima/man/simBmllag.Rd
pkg/yuima/man/wllag.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/R/IC.R
pkg/yuima/R/llag.R
pkg/yuima/man/IC.Rd
pkg/yuima/man/cce.factor.Rd
pkg/yuima/man/llag.Rd
Log:
New functions simBmllag and wllag are added. See NEWS for other revisions.
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/DESCRIPTION 2020-04-29 05:19:07 UTC (rev 736)
@@ -1,10 +1,10 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.10.1
+Version: 1.10.2
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature,
mvtnorm
-Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast
+Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast, wavethresh
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Description: Simulation and Inference for SDEs and Other Stochastic Processes.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/NAMESPACE 2020-04-29 05:19:07 UTC (rev 736)
@@ -63,6 +63,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
exportClasses("yuima",
"yuima.data",
@@ -172,6 +173,9 @@
export(mllag) # Oct. 10, 2015: multiple lead-lag detector
export(llag.test) # Sep. 9, 2017: testing the absence of lead-lag effects
export(cce.factor) # Jul. 19, 2019: high-dimensional covariance estimator
+export(simBmllag) # added by YK on Apr. 15, 2020
+export(simBmllag.coef) # added by YK on Apr. 15, 2020
+export(wllag) # added by YK on Apr. 15, 2020
export(get.zoo.data)
@@ -244,6 +248,7 @@
S3method(print, yuima.mllag) # Oct. 10, 2015
S3method(print, yuima.specv) # Oct. 10, 2015
S3method(print, yuima.ic)
+S3method(print, yuima.wllag) # added by YK on Apr. 15, 2020
S3method(toLatex, yuima)
S3method(toLatex, yuima.model)
@@ -254,5 +259,7 @@
S3method(plot, yuima.llag) # Oct. 10, 2015
S3method(plot, yuima.mllag) # Oct. 10, 2015
+S3method(plot, yuima.wllag) # added by YK on Apr. 15, 2020
+
useDynLib(yuima)
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/NEWS 2020-04-29 05:19:07 UTC (rev 736)
@@ -66,4 +66,8 @@
2020/2/6: fixed length(class(matrix))>2 issue. Code for simPoi fastened.
2020/3/5: should have fixed qmleLevy if() condition and rhck PROTECT problems
2020/4/25: modified IC.R, IC.Rd, and adaBayes.R
-
+2020/4/28: new functions simBmllag and wllag are added
+ the default value of tol of llag is changed to 1e-6
+ the wavethresh package is imported for implementation of wllag
+ fixed length(class(matrix))>2 issue
+ cce.factor.Rd, IC.Rd and llag.Rd are revised
\ No newline at end of file
Modified: pkg/yuima/R/IC.R
===================================================================
--- pkg/yuima/R/IC.R 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/R/IC.R 2020-04-29 05:19:07 UTC (rev 736)
@@ -567,13 +567,15 @@
print(x$BIC)
cat("\nQBIC:\n")
print(x$QBIC)
- if(class(x$CIC) == "matrix"){
+ #if(class(x$CIC) == "matrix"){
+ if(is.matrix(x$CIC)){ # fixed by YK
if(!is.null(x$CIC)){
cat("\nCIC:\n")
print(x$CIC)
}
}
- if(class(x$CIC) == "list"){
+ #if(class(x$CIC) == "list"){
+ if(is.list(class(x$CIC))){ # fixed by YK
if(!is.null(x$CIC$first)){
cat("\nCIC:\n")
print(x$CIC)
Modified: pkg/yuima/R/llag.R
===================================================================
--- pkg/yuima/R/llag.R 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/R/llag.R 2020-04-29 05:19:07 UTC (rev 736)
@@ -196,7 +196,7 @@
## 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, tol = 1e-7) standardGeneric("llag") )
+ ccor = ci, ci = FALSE, alpha = 0.01, fisher = TRUE, bw, tol = 1e-6) standardGeneric("llag") )
## yuima-method
setMethod("llag", "yuima", function(x, from, to, division, verbose, grid, psd, plot,
Added: pkg/yuima/R/simBmllag.r
===================================================================
--- pkg/yuima/R/simBmllag.r (rev 0)
+++ pkg/yuima/R/simBmllag.r 2020-04-29 05:19:07 UTC (rev 736)
@@ -0,0 +1,74 @@
+
+# sinc function
+sinc <- function(x){
+
+ out <- rep(1, length(x))
+
+ idx <- abs(x) > 1e-7
+ out[idx] <- sin(x[idx])/x[idx]
+
+ return(out)
+}
+
+# Littlewood-Paley wavelet function
+psi.lp <- function(x)
+ 2 * sinc(2 * pi * x) - sinc(pi * x)
+
+# coefficient matrices for circulant embedding
+simBmllag.coef <- function(n, J, rho, theta, delta = 1/2^(J+1)){
+
+ if(length(rho) < J + 1)
+ rho <- append(rho, double(J + 1 - length(rho)))
+
+ if(length(theta) < J + 1)
+ theta <- append(theta, double(J + 1 - length(theta)))
+
+ m <- 3^ceiling(log(2 * n - 2, base = 3))
+ M <- as.integer((m - 1)/2)
+
+ tl <- ((-M):M) * delta
+
+ c12 <- double(m)
+ for(j in 1:(J + 1)){
+ c12 <- c12 + 2^(J - j + 1) * delta^2 * rho[j] * psi.lp(2^(J - j + 1) * (tl - theta[j]))
+ }
+
+ c12 <- c(c12[-(1:M)], c12[1:M])
+
+ A12 <- fft(c12)
+
+ s <- sqrt(delta^2 - Mod(A12)^2)
+ t <- sqrt(2 * (delta + s))
+
+ a <- array(0, dim = c(m, 2, 2))
+
+ a[ ,1,1] <- (delta + s)/t
+ a[ ,2,2] <- a[ ,1,1]
+ a[ ,1,2] <- A12/t
+ a[ ,2,1] <- Conj(a[ ,1,2])
+
+ return(a)
+}
+
+# main function
+simBmllag <- function(n, J, rho, theta, delta = 1/2^(J+1),
+ imaginary = FALSE){
+
+ a <- simBmllag.coef(n, J, rho, theta, delta)
+ m <- dim(a)[1]
+
+ z1 <- rnorm(m) + 1i * rnorm(m)
+ z2 <- rnorm(m) + 1i * rnorm(m)
+
+ y1 <- a[ ,1,1] * z1 + a[ ,1,2] * z2
+ y2 <- a[ ,2,1] * z1 + a[ ,2,2] * z2
+
+ dW <- mvfft(cbind(y1, y2))[1:n, ]/sqrt(m)
+
+ if(imaginary == TRUE){
+ return(dW)
+ }else{
+ return(Re(dW))
+ }
+
+}
Added: pkg/yuima/R/wllag.r
===================================================================
--- pkg/yuima/R/wllag.r (rev 0)
+++ pkg/yuima/R/wllag.r 2020-04-29 05:19:07 UTC (rev 736)
@@ -0,0 +1,138 @@
+# Scale-by-scale lead-lag estimation by wavelets
+
+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)
+
+ 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")
+ }
+
+}
Modified: pkg/yuima/man/IC.Rd
===================================================================
--- pkg/yuima/man/IC.Rd 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/man/IC.Rd 2020-04-29 05:19:07 UTC (rev 736)
@@ -10,7 +10,9 @@
}
\usage{
-IC(drif = NULL, diff = NULL, data = NULL, Terminal = 1, add.settings = list(), start, lower, upper, ergodic = TRUE, stepwise = FALSE, weight = FALSE, rcpp = FALSE, ...)
+IC(drif = NULL, diff = NULL, data = NULL, Terminal = 1,
+ add.settings = list(), start, lower, upper, ergodic = TRUE,
+ stepwise = FALSE, weight = FALSE, rcpp = FALSE, ...)
}
\arguments{
@@ -141,7 +143,8 @@
}
## Candidate coefficients
-diffusion <- c("exp((theta11*cos(x)+theta12*sin(x)+theta13)/2)", "exp((theta11*cos(x)+theta12*sin(x))/2)",
+diffusion <- c("exp((theta11*cos(x)+theta12*sin(x)+theta13)/2)",
+ "exp((theta11*cos(x)+theta12*sin(x))/2)",
"exp((theta11*cos(x)+theta13)/2)", "exp((theta12*sin(x)+theta13)/2)")
drift <- c("theta21*x + theta22", "theta21*x")
@@ -158,8 +161,9 @@
ic1
## Ex.1.2 Stepwise
-ic2 <- IC(drif=drift, diff=diffusion, data=Xt, Terminal=Ter, start=para.init, lower=para.low,
- upper=para.upp, stepwise = TRUE, weight = FALSE, rcpp = TRUE)
+ic2 <- IC(drif=drift, diff=diffusion, data=Xt, Terminal=Ter,
+ start=para.init, lower=para.low, upper=para.upp,
+ stepwise = TRUE, weight = FALSE, rcpp = TRUE)
ic2
\dontrun{
@@ -204,12 +208,14 @@
## Ex.2.1 Joint
ic3 <- IC(drif=drift, diff=diffusion, data=Xt, Terminal=Ter, add.settings=modsettings,
- start=para.init, lower=para.low, upper=para.upp, weight=FALSE, rcpp=FALSE)
+ start=para.init, lower=para.low, upper=para.upp,
+ weight=FALSE, rcpp=FALSE)
ic3
## Ex.2.2 Stepwise
ic4 <- IC(drif=drift, diff=diffusion, data=Xt, Terminal=Ter, add.settings=modsettings,
- start=para.init, lower=para.low, upper=para.upp, stepwise = TRUE, weight=FALSE, rcpp=FALSE)
+ start=para.init, lower=para.low, upper=para.upp,
+ stepwise = TRUE, weight=FALSE, rcpp=FALSE)
ic4
}
Modified: pkg/yuima/man/cce.factor.Rd
===================================================================
--- pkg/yuima/man/cce.factor.Rd 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/man/cce.factor.Rd 2020-04-29 05:19:07 UTC (rev 736)
@@ -63,7 +63,7 @@
}
\item{weight}{
%% ~~Describe \code{weight} here~~
-logical. If \code{TRUE}. A weighted version is used for \code{regularize = "glasso"} as in Koike (2019).
+logical. If \code{TRUE}, a weighted version is used for \code{regularize = "glasso"} as in Koike (2020).
}
\item{nlambda}{
%% ~~Describe \code{nlambda} here~~
@@ -145,10 +145,10 @@
\item \code{regularize = "glasso"} (the default).
-This performs the glaphical Lasso. When \code{weight=TRUE} (the default), a weighted version of the graphical Lasso is performed as in Koike (2019). Otherwise, the standard graphical Lasso is performed as in Brownlees et al. (2018).
+This performs the glaphical Lasso. When \code{weight=TRUE} (the default), a weighted version of the graphical Lasso is performed as in Koike (2020). Otherwise, the standard graphical Lasso is performed as in Brownlees et al. (2018).
If \code{lambda="aic"} (resp.~\code{lambda="bic"}), the penalty parameter for the graphical Lasso is selected by minimizing the formally defined AIC (resp.~BIC).
-The minimization is carried out by grid search, where the grid is determined as in Section 5.1 of Koike (2019).
+The minimization is carried out by grid search, where the grid is determined as in Section 5.1 of Koike (2020).
The optimization problem in the graphical Lasso is solved by the GLASSOFAST algorithm of Sustik and Calderhead (2012), which is available from the package \pkg{glassoFast}.
@@ -207,9 +207,9 @@
A blocking and regularization approach to high-dimensional realized covariance estimation,
\emph{Journal of Applied Econometrics}, \bold{27}, 625--645.
-Koike, Y. (2019).
-De-biased graphical Lasso for high-frequency data.
-\href{https://arxiv.org/abs/1905.01494}{arXiv:1905.01494}.
+Koike, Y. (2020).
+ De-biased graphical Lasso for high-frequency data,
+ \emph{Entropy}, \bold{22}, 456.
Pourahmadi, M. (2011).
Covariance estimation: The GLM and regularization perspectives.
Modified: pkg/yuima/man/llag.Rd
===================================================================
--- pkg/yuima/man/llag.Rd 2020-04-27 07:20:36 UTC (rev 735)
+++ pkg/yuima/man/llag.Rd 2020-04-29 05:19:07 UTC (rev 736)
@@ -7,7 +7,7 @@
\usage{
llag(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, tol = 1e-7)
+ fisher = TRUE, bw, tol = 1e-6)
}
\arguments{
Added: pkg/yuima/man/simBmllag.Rd
===================================================================
--- pkg/yuima/man/simBmllag.Rd (rev 0)
+++ pkg/yuima/man/simBmllag.Rd 2020-04-29 05:19:07 UTC (rev 736)
@@ -0,0 +1,174 @@
+\name{simBmllag}
+\alias{simBmllag}
+\alias{simBmllag.coef}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+%% ~~function to do ... ~~
+Simulation of increments of bivariate Brownian motions with multi-scale lead-lag relationships
+}
+\description{
+%% ~~ A concise (1-5 lines) description of what the function does. ~~
+This function simulates increments of bivariate Brownian motions with multi-scale lead-lag relationships introduced in Hayashi and Koike (2018a) by the multi-dimensional circulant embedding method of Chan and Wood (1999).
+}
+\usage{
+simBmllag(n, J, rho, theta, delta = 1/2^(J + 1), imaginary = FALSE)
+simBmllag.coef(n, J, rho, theta, delta = 1/2^(J + 1))
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{n}{
+%% ~~Describe \code{n} here~~
+the number of increments to be simulated.
+}
+ \item{J}{
+%% ~~Describe \code{J} here~~
+a positive integer to determine the finest time resolution: \code{2^(-J-1)} is regarded as the finest time resolution.
+}
+ \item{rho}{
+%% ~~Describe \code{rho} here~~
+a vector of scale-by-scale correlation coefficients. If \code{length(rho) < J}, zeros are appended to make the length equal to \code{J}.
+}
+ \item{theta}{
+%% ~~Describe \code{theta} here~~
+a vector of scale-by-scale lead-lag parameters. If \code{length(theta) < J}, zeros are appended to make the length equal to \code{J}.
+}
+ \item{delta}{
+%% ~~Describe \code{delta} here~~
+the step size of time increments. This must be smaller than or equal to \code{2^(-J-1)}.
+}
+ \item{imaginary}{
+%% ~~Describe \code{imaginary} here~~
+logical. %If \code{TRUE}, the function returns complex-valued increments whose real and imaginary parts are independent and have the same law as
+See `Details'.
+}
+}
+\details{
+%% ~~ If necessary, more details than the description above ~~
+Let \eqn{B(t)} be a bivariate Gaussian process with stationary increments such that its marginal processes are standard Brownian motions and its cross-spectral density is given by Eq.(14) of Hayashi and Koike (2018a). The function \code{simBmllag} simulates the increments \eqn{B(i\delta)-B((i-1)\delta)}, \eqn{i=1,\dots,n}. The parameters \eqn{R_j} and \eqn{theta_j} in Eq.(14) of Hayashi and Koike (2018a) are specified by \code{rho} and \code{theta}, while \eqn{\delta} and \eqn{n} are specified by \code{delta} and \code{n}, respecitively.
+
+Simulation is implemented by the multi-dimensional circulant embedding algorithm of Chan and Wood (1999). The last step of this algorithm returns a bivariate complex-valued sequence whose real and imaginary parts are independent and has the same law as \eqn{B(k\delta)-B((k-1)\delta)}, \eqn{k=1,\dots,n}; see Step 3 of Chan and Wood (1999, Section 3.2).
+If \code{imaginary = TRUE}, the function \code{simBmllag} directly returns this bivariate complex-valued sequence, so we obtain two sets of simulated increments of \eqn{B(t)} by taking its real and complex parts. If \code{imaginary = FALSE} (default), the function returns only the real part of this sequence, so we directly obtain simulated increments of \eqn{B(t)}.
+
+The function \code{simBmllag.coef} is internally used to compute the sequence of coefficient matrices \eqn{R(k)\Lambda(k)^{1/2}} in Step 2 of Chan and Wood (1999, Section 3.2). This procedure can be implemented before generating random numbers. %Thus, when we conduct a Monte Carlo simulation for \eqn{(B(k\delta)-B((k-1)\delta))_{k=1}^n} with a fixed set of parameters, it is enough to implement the function \code{simBmllag.coef} only one time.
+Since this step typically takes the most computational cost, this function is useful to reduce computational time when we conduct a Monte Carlo simulation for \eqn{(B(k\delta)-B((k-1)\delta))_{k=1}^n} with a fixed set of parameters. See `Examples' for how to use this function to simulate \eqn{(B(k\delta)-B((k-1)\delta))_{k=1}^n}.
+}
+\value{
+%% ~Describe the value returned
+%% If it is a LIST, use
+%% \item{comp1 }{Description of 'comp1'}
+%% \item{comp2 }{Description of 'comp2'}
+%% ...
+\code{simBmllag} returns a \code{n} x 2 matrix if \code{imaginary = FALSE} (default). Otherwise, \code{simBmllag} returns a complex-valued \code{n} x 2 matrix.
+
+\code{simBmllag.coef} returns a complex-valued \eqn{m} x 2 x 2 array, where \eqn{m} is an integer determined by the rule described at the end of Chan and Wood (1999, Section 2.3).
+}
+\references{
+%% ~put references to the literature/web site here ~
+Chan, G. and Wood, A. T. A. (1999).
+ Simulation of stationary Gaussian vector fields,
+ \emph{Statistics and Computing}, \bold{9}, 265--268.
+
+Hayashi, T. and Koike, Y. (2018a).
+ Wavelet-based methods for high-frequency lead-lag analysis,
+ \emph{SIAM Journal of Financial Mathematics}, \bold{9}, 1208--1248.
+
+Hayashi, T. and Koike, Y. (2018b).
+Multi-scale analysis of lead-lag relationships in high-frequency financial markets.
+\href{https://arxiv.org/abs/1708.03992}{arXiv:1708.03992}.
+}
+\author{
+%% ~~who you are~~
+Yuta Koike with YUIMA project Team
+}
+\note{
+%% ~~further notes~~
+There are typos in the first and second displayed equations in page 1221 of Hayashi and Koike (2018a): The \eqn{j}-th summands on their right hand sides should be multiplied by \eqn{2^j}.
+}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+\code{\link{wllag}}
+}
+\examples{
+## Example 1
+## Simulation setting of Hayashi and Koike (2018a, Section 4).
+
+n <- 15000
+J <- 13
+
+rho <- c(0.3,0.5,0.7,0.5,0.5,0.5,0.5,0.5)
+theta <- c(-1,-1, -2, -2, -3, -5, -7, -10)/2^(J + 1)
+
+set.seed(123)
+
+dB <- simBmllag(n, J, rho, theta)
+str(dB)
+n/2^(J + 1) # about 0.9155
+sum(dB[ ,1]^2) # should be close to n/2^(J + 1)
+sum(dB[ ,2]^2) # should be close to n/2^(J + 1)
+
+# Plot the sample path of the process
+B <- apply(dB, 2, "diffinv") # construct the sample path
+Time <- seq(0, by = 1/2^(J+1), length.out = n) # Time index
+plot(zoo(B, Time), main = "Sample path of B(t)")
+
+# Using simBmllag.coef to implement the same simulation
+a <- simBmllag.coef(n, J, rho, theta)
+m <- dim(a)[1]
+
+set.seed(123)
+
+z1 <- rnorm(m) + 1i * rnorm(m)
+z2 <- rnorm(m) + 1i * rnorm(m)
+y1 <- a[ ,1,1] * z1 + a[ ,1,2] * z2
+y2 <- a[ ,2,1] * z1 + a[ ,2,2] * z2
+dW <- mvfft(cbind(y1, y2))[1:n, ]/sqrt(m)
+dB2 <- Re(dW)
+
+plot(diff(dB - dB2)) # identically equal to zero
+
+
+## Example 2
+## Simulation Scenario 2 of Hayashi and Koike (2018b, Section 5).
+
+# Simulation of Bm driving the log-price processes
+n <- 30000
+J <- 14
+
+rho <- c(0.3,0.5,0.7,0.5,0.5,0.5,0.5,0.5)
+theta <- c(-1,-1, -2, -2, -3, -5, -7, -10)/2^(J + 1)
+
+dB <- simBmllag(n, J, rho, theta)
+
+# Simulation of Bm driving the volatility processes
+R <- -0.5 # leverage parameter
+delta <- 1/2^(J+1) # step size of time increments
+dW1 <- R * dB[ ,1] + sqrt(1 - R^2) * rnorm(n, sd = sqrt(delta))
+dW2 <- R * dB[ ,2] + sqrt(1 - R^2) * rnorm(n, sd = sqrt(delta))
+
+# Simulation of the model by the simulate function
+dW <- rbind(dB[,1], dB[,2], dW1, dW2) # increments of the driving Bm
+
+# defining the yuima object
+drift <- c(0, 0, "kappa*(eta - x3)", "kappa*(eta - x4)")
+diffusion <- diag(4)
+diag(diffusion) <- c("sqrt(max(x3,0))", "sqrt(max(x4,0))",
+ "xi*sqrt(max(x3,0))", "xi*sqrt(max(x4,0))")
+xinit <- c(0,0,"rgamma(1, 2*kappa*eta/xi^2,2*kappa/xi^2)",
+ "rgamma(1, 2*kappa*eta/xi^2,2*kappa/xi^2)")
+mod <- setModel(drift = drift, diffusion = diffusion,
+ xinit = xinit, state.variable = c("x1","x2","x3","x4"))
+samp <- setSampling(Terminal = n * delta, n = n)
+yuima <- setYuima(model = mod, sampling = samp)
+
+# simulation
+result <- simulate(yuima, increment.W = dW,
+ true.parameter = list(kappa = 5, eta = 0.04, xi = 0.5))
+
+plot(result)
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}% use one of RShowDoc("KEYWORDS")
Added: pkg/yuima/man/wllag.Rd
===================================================================
--- pkg/yuima/man/wllag.Rd (rev 0)
+++ pkg/yuima/man/wllag.Rd 2020-04-29 05:19:07 UTC (rev 736)
@@ -0,0 +1,153 @@
+\name{wllag}
+\alias{wllag}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+%% ~~function to do ... ~~
+Scale-by-scale lead-lag estimation
+}
+\description{
+%% ~~ A concise (1-5 lines) description of what the function does. ~~
+This function estimates lead-lag parameters on a scale-by-scale basis from non-synchronously observed bivariate processes, using the estimatiors proposed in Hayashi and Koike (2018b).
+}
+\usage{
+wllag(x, y, J = 8, N = 10, tau = 1e-3, from = -to, to = 100,
+ verbose = FALSE, in.tau = FALSE, tol = 1e-6)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{x}{
+%% ~~Describe \code{x} here~~
+a \code{\link{zoo}} object for observation data of the first process.
+}
+ \item{y}{
+%% ~~Describe \code{y} here~~
+a \code{\link{zoo}} object for observation data of the second process.
+}
+ \item{J}{
+%% ~~Describe \code{J} here~~
+a positive integer. Scale-by scale lead-lag parameters are estimated up to the level \code{J}.
+}
+ \item{N}{
+%% ~~Describe \code{N} here~~
+The number of vanishing moments of Daubechies' compactly supported wavelets, passed to \code{\link{PsiJ}}. This should be an integer between 1 and 10.
+}
+\item{tau}{
+%% ~~Describe \code{dgrid} here~~
+the step size of a finite grid on which objective functions are evaluated. Note that this value is identified with the finest time resolution of the underlying model.
+The default value \code{1e-3} corresponds to 1 mili-second if the unit time corresponds to 1 second.
+}
+ \item{from}{
+%% ~~Describe \code{from} here~~
+a negative integer. \code{from*tau} gives the lower end of a finite grid on which objective functions are evaluated.
+}
+ \item{to}{
+%% ~~Describe \code{to} here~~
+a positive integer. \code{to*tau} gives the upper end of a finite grid on which objective functions are evaluated.
+}
+ \item{verbose}{
+%% ~~Describe \code{verbose} here~~
+a logical. If \code{FALSE} (default), the function returns only the estimated scale-by-scale lead-lag parameters. Otherwise, the function also returns some other statistics such as values of the signed objective functions. See `Value'.
+}
+\item{in.tau}{
+%% ~~Describe \code{delta} here~~
+a logical. If \code{TRUE}, the estimated lead-lag parameters are returned in increments of \code{tau}. That is, the estimated lead-lag parameters are divided by \code{tau}.
+}
+\item{tol}{
+%% ~~Describe \code{delta} here~~
+tolelance parameter to avoid numerical errors in comparison of time stamps. All time stamps are divided by \code{tol} and rounded to integers. A reasonable choice of \code{tol} is the minimum unit of time stamps. The default value \code{1e-6} supposes that the minimum unit of time stamps is greater or equal to 1 micro-second.
+}
+}
+\details{
+%% ~~ If necessary, more details than the description above ~~
+Hayashi and Koike (2018a) introduced a bivariate continuous-time model having different lead-lag relationships at different time scales. The wavelet cross-covariance functions of this model, computed based on the Littlewood-Paley wavelets, have unique maximizers in absolute values at each time scale. These maximizer can be considered as lead-lag parameters at each time scale. To estimate these parameters from discrete observation data, Hayashi and Koike (2018b) constructed objective functions mimicking behavior of the wavelet cross-covariance functions of the underlying model. Then, estimates of the scale-by-scale lead-lag parameters can be obtained by maximizing these objective functions in absolute values.
+}
+\value{
+If \code{verbose} is \code{FALSE}, a numeric vector with length \code{J}, corresponding to the estimated scale-by-scale lead-lag parameters, is returned. Note that their positive values indicate that the first process leads the second process.
+
+Otherwise, an object of class \code{"yuima.wllag"}, which is a list with the following components, is returned:
+\item{lagtheta}{the estimated scale-by-scale lead-lag parameters. The \eqn{j} th component corresponds to the estimate at the level \eqn{j}. A positive value indicates that the first process leads the second process.}
+\item{obj.values}{the values of the objective functions evaluated at the estimated lead-lag parameters.}
+\item{obj.fun}{a list of values of the objective functions. The \eqn{j} th component of the list corresponds to a \code{\link{zoo}} object for values of the signed objective function at the level \eqn{j} indexed by the search grid.}
+\item{theta.hry}{the lead-lag parameter estimate in the sense of Hoffmann, Rosenbaum and Yoshida (2013).}
+\item{cor.hry}{the correltion coefficient in the sense of Hoffmann, Rosenbaum and Yoshida (2013), evaluated at the estimated lead-lag parameter.}
+\item{ccor.hry}{a \code{\link{zoo}} object for values of the cross correltion function in the sense of Hoffmann, Rosenbaum and Yoshida (2013) indexed by the search grid.}
+%% ~Describe the value returned
+%% If it is a LIST, use
+%% \item{comp1 }{Description of 'comp1'}
+%% \item{comp2 }{Description of 'comp2'}
+%% ...
+}
+\references{
+%% ~put references to the literature/web site here ~
+Hayashi, T. and Koike, Y. (2018a).
+ Wavelet-based methods for high-frequency lead-lag analysis,
+ \emph{SIAM Journal of Financial Mathematics}, \bold{9}, 1208--1248.
+
+Hayashi, T. and Koike, Y. (2018b).
+Multi-scale analysis of lead-lag relationships in high-frequency financial markets.
+\href{https://arxiv.org/abs/1708.03992}{arXiv:1708.03992}.
+
+Hoffmann, M., Rosenbaum, M. and Yoshida, N. (2013)
+ Estimation of the lead-lag parameter from non-synchronous data,
+ \emph{Bernoulli}, \bold{19}, no. 2, 426--461.
+}
+\author{
+%% ~~who you are~~
+Yuta Koike with YUIMA Project Team
+
+}
+\note{
+%% ~~further notes~~
+Smaller levels correspond to finer time scales. In particular, the first level corresponds to the finest time resolution, which is defined by the argument \code{tau}.
+
+If there are multiple maximizers in an objective function, \code{wllag} takes a maximizer farthest from zero (if there are two such values, the function takes the negative one). This behavior is different from \code{\link{llag}}.
+
+The objective functions themselves do NOT consitently estimate the corresponding wavelet covariance functions. This means that values in \code{obj.values} and \code{obj.fun} cannot be interpreted as covaraince estimates (their scales depend on the degree of non-synchronicity of observation times).
+}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+\code{\link{simBmllag}}, \code{\link{PsiJ}}, \code{\link{llag}}
+}
+\examples{
+## An example from a simulation setting of Hayashi and Koike (2018b)
+set.seed(123)
+
+# Simulation of Bm driving the log-price processes
+n <- 15000
+J <- 13
+tau <- 1/2^(J+1)
+
+rho <- c(0.3,0.5,0.7,0.5,0.5,0.5,0.5,0.5)
+theta <- c(-1,-1, -2, -2, -3, -5, -7, -10) * tau
+
+dB <- simBmllag(n, J, rho, theta)
+
+Time <- seq(0, by = tau, length.out = n) # Time index
+x <- zoo(diffinv(dB[ ,1]), Time) # simulated path of the first process
+y <- zoo(diffinv(dB[ ,2]), Time) # simulated path of the second process
+
+# Generate non-synchronously observed data
+x <- x[as.logical(rbinom(n + 1, size = 1, prob = 0.5))]
+y <- y[as.logical(rbinom(n + 1, size = 1, prob = 0.5))]
+
+# Estimation of scale-by-scale lead-lag parameters (compare with theta/tau)
+wllag(x, y, J = 8, tau = tau, tol = tau, in.tau = TRUE)
+# Estimation with other information
+out <- wllag(x, y, tau = tau, tol = tau, in.tau = TRUE, verbose = TRUE)
+out
+
+# Plot of the HRY cross-correlation function
+plot(out$ccor.hry, xlab = expression(theta),
+ ylab = expression(U(theta)))
+
+# Plot of the objective functions
+op <- par(mfrow = c(4,2))
+plot(out)
+par(op)
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}% use one of RShowDoc("KEYWORDS")
More information about the Yuima-commits
mailing list