# Simulation of Increments of Bivariate Brownian Motions with Multi-scale Lead-lag Relationships 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). ```r simBmllag(n, J, rho, theta, delta = 1/2^(J + 1), imaginary = FALSE) simBmllag.coef(n, J, rho, theta, delta = 1/2^(J + 1)) ``` ## Arguments - `n`: the number of increments to be simulated. - `J`: a positive integer to determine the finest time resolution: `2^(-J-1)` is regarded as the finest time resolution. - `rho`: a vector of scale-by-scale correlation coefficients. If `length(rho) < J`, zeros are appended to make the length equal to `J`. - `theta`: a vector of scale-by-scale lead-lag parameters. If `length(theta) < J`, zeros are appended to make the length equal to `J`. - `delta`: the step size of time increments. This must be smaller than or equal to `2^(-J-1)`. - `imaginary`: logical. See `Details'. ## Details Let $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 `simBmllag` simulates the increments $B(i\delta)-B((i-1)\delta)$, $i=1,\dots,n$. The parameters $R_j$ and $theta_j$ in Eq.(14) of Hayashi and Koike (2018a) are specified by `rho` and `theta`, while $\delta$ and $n$ are specified by `delta` and `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 $B(k\delta)-B((k-1)\delta)$, $k=1,\dots,n$; see Step 3 of Chan and Wood (1999, Section 3.2). If `imaginary = TRUE`, the function `simBmllag` directly returns this bivariate complex-valued sequence, so we obtain two sets of simulated increments of $B(t)$ by taking its real and complex parts. If `imaginary = FALSE` (default), the function returns only the real part of this sequence, so we directly obtain simulated increments of $B(t)$. The function `simBmllag.coef` is internally used to compute the sequence of coefficient matrices $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. 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 $(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 $(B(k\delta)-B((k-1)\delta))_{k=1}^n$. ## Returns `simBmllag` returns a `n` x 2 matrix if `imaginary = FALSE` (default). Otherwise, `simBmllag` returns a complex-valued `n` x 2 matrix. `simBmllag.coef` returns a complex-valued $m$ x 2 x 2 array, where $m$ is an integer determined by the rule described at the end of Chan and Wood (1999, Section 2.3). ## References Chan, G. and Wood, A. T. A. (1999). Simulation of stationary Gaussian vector fields, **Statistics and Computing**, 9 , 265--268. Hayashi, T. and Koike, Y. (2018a). Wavelet-based methods for high-frequency lead-lag analysis, **SIAM Journal of Financial Mathematics**, 9 , 1208--1248. Hayashi, T. and Koike, Y. (2018b). Multi-scale analysis of lead-lag relationships in high-frequency financial markets. tools:::Rd_expr_doi("10.48550/arXiv.1708.03992") . ## Author(s) Yuta Koike with YUIMA project Team ## Note There are typos in the first and second displayed equations in page 1221 of Hayashi and Koike (2018a): The $j$-th summands on their right hand sides should be multiplied by $2^j$. ## See Also `wllag` ## Examples ```r ## 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) ```

simBmllag function