# Wavelet-based and Fourier-based multivariate Whittle estimation: multiwave

Sophie Achard, Irène Gannaz

<figure><img src="image_1.png" /></figure>

To cite this version:

Sophie Achard, Irène Gannaz. Wavelet-based and Fourier-based multivariate Whittle estimation: multiwave. Journal of Statistical Software, 2019, 89 (6), 10.18637/jss.v089.i06. hal-01840099

HAL Id: hal-01840099

https://hal.science/hal-01840099

Submitted on 22 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wavelet-based and Fourier-based multivariate Whittle estimation: multiwave

Sophie Achardª and Irène Gannazᵇ

ªUniv. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France

ᵇ Univ. Lyon, INSA de Lyon, CNRS UMR 5208, Institut Camille Jordan 20 avenue Albert Einstein, F-69621 Villeurbanne Cedex, France

## Abstract

Multivariate time series with long-dependence are observed in many applications such as finance, geophysics or neuroscience. Many packages provide estimation tools for univariate settings but few are addressing the problem of long-dependence estimation for multivariate settings. The package **multiwave** is providing efficient estimation procedures for multivariate time series. Two semi-parametric estimation methods of the long-memory exponents and long-run covariance matrix of time series are implemented. The first one is the Fourier-based estimation proposed by [18] and the second one is a wavelet-based estimation described in [4]. The objective of this paper is to provide an overview of the R package **multiwave** with its practical application perspectives.

**Keywords.** wavelets, multivariate time series, Whittle estimation, long-memory properties, long-run covariance, R

## 1 Introduction

Time series with defined autocovariance functions are said to present long-memory or long-range dependency when their autocovariance function is decreasing very slowly, slower than an exponential decay. More precisely, let $g(\cdot)$ be the autocovariance function of a time series $X$. $X$ is said to be long-memory if there exists $\alpha$, $0 < \alpha < 1$, such that $g(t)$ is asymptotically equivalent to $|t|^{-\alpha}$ when $t \to +\infty$ (see [5] and references therein). This definition implies that the covariance function is not summable. Equivalently, the spectral density $f(\cdot)$, if it exists, is such that $f(\lambda)$ is equivalent up to a constant to $|\lambda|^{1-\alpha}$ when $\lambda \to 0^+$. In this case, when trying to estimate the expectation using the empirical mean of long-memory time series, the variance of the estimator is not decreasing to 0 as $N^{-1}$ (where $N$ is the sample size). Hence it is crucial to take into account the presence of long-memory for defining good estimators [5]. In the case of univariate time series, several very efficient approaches have been developed and validated. A web page entitled “Time Series Analysis”<fnref n="1" /> from the R software is providing a very exhaustive list of methods and softwares dealing with long-memory time series. Among others, we can cite the packages **fracdiff** [7], **arfima** and **FGN** [23], **longmemo** [5] and

<footnote n="1">¹https://cran.r-project.org/web/views/TimeSeries.html</footnote>
**forecast** [8]. For example, **fracdiff** is dedicated to simulation of fractional ARIMA time series and to estimation using regression of the periodogram. **longmemo** provides real data examples of time series with long-memory properties.

Approaches for multivariate long-memory time series are less developed. When dealing with multi-variate time series, an important quantity to estimate is the covariance or correlation between pairs of time series. The effect of the presence of long-memory on this estimation is obvious, as stated by [14]. One R package, **waveslim**, is dedicated to the wavelet correlation analysis for pairs of random variables [24] but long-range dependence properties are not considered. [17] provide R code<fnref n="2" /> for bivariate long-range dependent time series with parametric estimations. The objective of this paper is to provide an efficient R package, called **multiwave**, to estimate the long-memory parameters and covariance matrices for multivariate time series. The estimation procedures are based on a semi-parametric approach, which is robust to model misspecification.

The procedures are also suited for dealing with more than two dimensional data. Indeed they are based on Whittle approximation which provides a simple function to optimize. This function can be used for any dimension of the problem. In comparison, regression of the scalogram or periodogram [3] is based on a linear fit of pairs of time series, and thus there does not exist an easy way to extend to more than two dimensions.

**multiwave** package is based on [4], where we developed a wavelet-based approach using Whittle approximation for an efficient estimation of the long-memory parameters and the long-run covariance matrices. In addition, **multiwave** proposes an implementation of an alternative method using Fourier decomposition as described in [18]<fnref n="3" />.

As it is described in this paper, **multiwave** is a very versatile package and opens the way to estimation of the long-memory parameters and the long-run covariance matrices using multivariate data sets. It is in particular not necessary to assume the stationarity of the time series as it is the case when using Fourier decomposition [6]. The Whittle approximation is computed using either the coefficients of wavelet decomposition or the coefficients of Fourier decomposition when the time series are stationary.

The package **multiwave** is divided in three parts. A first group of functions is dedicated to the simulation of multivariate long-memory time series; the main function is **fivarma**. A second group of functions is implementing the wavelet decomposition, through **DWTexact** and associated functions. Finally the computation of the estimators are coded using the Fourier decomposition in **mfw** and its derivatives and using the wavelet decomposition in **mww** and its derivatives.

The mathematical background is detailed in a separate Section 2. The rest of the paper is dedicated to the description of the package **multiwave**. Simple examples of parametric models and real data are presented in Section 3 with corresponding functions of **multiwave** ready to apply. Core estimation functions using wavelets and Fourier transform are detailed in Section 4 along with pieces of code using simulated time series. Finally, practical considerations are discussed in the three last sections. Section 5 is discussing the practical choices of parameters. Comparaisons between wavelets and Fourier approaches are described in Section 6. And an application to real data in neuroscience is concluding the paper, Section 7.

<footnote n="2">²freely available http://pages.stern.nyu.edu/~rsela/VARFI/code.html</footnote>

<footnote n="3">³code available in matlab http://shimotsu.web.fc2.com/Site/Matlab_Codes.html</footnote>
## 2 Theoretical background

As in the univariate case, the definition of long-memory for a p-vector process is based on the asymptotic behaviour of the cross-spectral density in the neighbourhood of zero [12]. We consider N observations of a long-memory p-vector process $\mathbf{X} = \{X_\ell(k), k \in \mathbb{Z}, \ell = 1, \dots, p\}$, namely $\mathbf{X}(1), \dots, \mathbf{X}(N)$. $\mathbf{X}$ is said to be a multivariate $M(\mathbf{d})$ process when for each $\ell = 1, \dots, p$ there exists $D_\ell \in \mathbb{N}$ such that the $D_\ell$-th order difference $\Delta^{D_\ell} X_\ell$ is covariance stationary. In addition, let us assume that for any $\ell, m = 1, \dots, p$ the generalized cross-spectral density of $X_\ell$ and $X_m$ is

$$f_{\ell,m}(\lambda) = \frac{1}{2\pi} \Omega_{\ell,m} (1 - e^{-i\lambda})^{-d_\ell} (1 - e^{i\lambda})^{-d_m} f_{\ell,m}^{S}(\lambda), \quad \lambda \in [-\pi, \pi], \quad (1)$$

with $\Omega = (\Omega_{\ell,m})_{\ell,m=1,\dots,p}$ an Hermitian matrix. The functions $f_{\ell,m}^S(\cdot)$ correspond to the short-memory dynamics of the process. The parameters $d_\ell$ satisfies $-1/2 < d_\ell - D_\ell < 1/2$. More generally, the wavelet-based procedure is available for cross-spectral density satisfying an approximation

$$f(\lambda) \sim \Lambda(\mathbf{d})\Omega\Lambda(\mathbf{d})^*, \quad \text{when } \lambda \to 0, \quad \text{with } \Lambda(\mathbf{d}) = \text{diag}(|\lambda|^{-\mathbf{d}}e^{-i\text{sign}(\lambda)\pi\mathbf{d}/2}), \quad (2)$$

where the exponent $*$ denotes the conjugate transpose operator. Here and subsequently $\sim$ means that the ratio of the left- and right-hand sides converges to one. Note that, the process $X_\ell$ is not necessarily stationary.

The long-range dependence parameter measures the power-like rate of decay of the autocovariance function. The long-run covariance matrix $\Omega$ can be seen as the covariance at low frequencies between the time series. It gives a quantification of the link between the components of the multivariate time series. The long-run covariance parameter of the model is free from the difference in the autocorrelation behaviour of each component. It is linked with long-run correlations $(\Omega_{\ell,m}/\sqrt{\Omega_{\ell,\ell}\Omega_{m,m}})_{\ell,m=1,\dots,p}$, which are also encountered in literature as power-law coherencies between two time series [17] or as fractal connectivities [3].

### 2.1 A parametric example: FIVARMA

Fractionally Integrated Vector Auto Regressive Moving Average (FIVARMA) processes are parametric models with a spectral density satisfying approximation (2). They correspond to Model A of [9]. We refer to this paper for a detailed mathematical description.

Let $\mathbf{u}$ be a p-dimensional white noise with $\mathbb{E}[\mathbf{u}(t) | \mathcal{F}_{t-1}] = 0$ and $\mathbb{E}[\mathbf{u}(t)\mathbf{u}^T(t) | \mathcal{F}_{t-1}] = \mathbf{\Sigma}$, where $\mathcal{F}_{t-1}$ is the $\sigma$-field generated by $\{\mathbf{u}(s), s < t\}$, and $\mathbf{\Sigma}$ is a positive definite matrix.

Let $(\mathbf{A}_k)_{k \in \mathbb{N}}$ be a sequence of $\mathbb{R}^{p \times p}$-valued matrices with $\mathbf{A}_0$ the identity matrix and $\sum_{k=0}^{\infty} \|\mathbf{A}_k\|^2 < \infty$. The discrete Fourier transform of the sequence is denoted $\mathbf{A}(\cdot)$, that is $\mathbf{A}(\lambda) = \sum_{k=0}^{\infty} \mathbf{A}_k e^{ik\lambda}$. We assume that all the roots of $|\mathbf{A}(\mathbb{L})|$ are outside the closed unit circle, where $\mathbb{L}$ denotes the lag operator. Let also $(\mathbf{B}_k)_{k \in \mathbb{N}}$ be a sequence in $\mathbb{R}^{p \times p}$ with $\mathbf{B}_0$ the identity matrix and $\sum_{k=0}^{\infty} \|\mathbf{B}_k\|^2 < \infty$. As defined for $\mathbf{A}$, $\mathbf{B}(\cdot)$ denotes the discrete Fourier transform of the sequence, $\mathbf{B}(\lambda) = \sum_{k=0}^{\infty} B_k e^{ik\lambda}$. Let $\mathbf{X}$ be defined by

$$\mathbf{A}(\mathbb{L}) \text{diag}(\mathbb{1} - \mathbb{L})^{\mathbf{d}} \mathbf{X}(t) = \mathbf{B}(\mathbb{L}) \mathbf{u}(t). \quad (3)$$

The spectral density satisfies
$$f_{\ell,m}(\lambda) \sim_{\lambda \to 0^+} \frac{1}{2\pi} \Omega_{\ell,m} e^{-i\pi/(2(d_{\ell}-d_m))} \lambda^{-(d_{\ell}+d_m)}$$

with

$$\Omega = \mathbf{A}(1)^{-1}\mathbf{B}(1)\mathbf{\Sigma}\mathbf{B}(1)^T\mathbf{A}(1)^{-1}. \quad (4)$$

Then $\mathbf{X}$ is called a FIVARMA($d, q$) process and satisfies approximation (2).

## Limits of the model

Note that in definition (3) the operators are applied in a given order, where the lag operator is taken first. Changing the order of the lag operator and autoregression corresponds to model B of [9] and VARFI models of [17] where $\mathbf{X}$ is obtained with equation diag($\mathbf{1} - \mathbf{L}$)$^d \mathbf{A}(\mathbf{L}) \mathbf{X}(t) = \mathbf{B}(\mathbf{L})\mathbf{u}(t)$, with similar notations than above. That is, $\mathbf{X}$ is obtained by fractional integration after autoregression, which is also called cointegration. The spectral density still satisfies the approximation (2) however the matrix $\Omega$ may no longer be Hermitian. Clearly multiwave package is not built to deal with such cases. We refer to alternative methods in literature, among others [15, 17, 19]. Taking into account cointegration is a difficult problem that exceeds the scope of this paper. Future work is needed to handle this particular case.

## 2.2 Fourier-based estimation (MFW)

The discrete Fourier transform and the periodogram of $\mathbf{X}$ evaluated at frequency $\lambda$ are defined as in [18]'s procedure

$$\mathbf{W}^F(\lambda) = \frac{1}{\sqrt{2\pi N}} \sum_{t=1}^{N} \mathbf{X}(t)e^{it\lambda}, \\ \mathbf{I}^F(\lambda) = \mathbf{W}^F(\lambda)\mathbf{W}^F(\lambda)^*.$$

Let $\lambda_j = 2\pi j/N$, $j = 1, \dots, m$, be the Fourier frequencies used in estimation, $m \in \mathbb{N}$. Define $\mathbf{\Lambda}_j^F(\mathbf{d}) = \text{diag}(\lambda_j \mathbf{d} e^{i(\pi - \lambda_j)\mathbf{d}/2})$. The estimators $(\hat{\mathbf{d}}^{MFW}, \hat{\mathbf{\Omega}}^{MFW})$ are minimizers of the criterion

$$\mathcal{L}^{MFW}(\mathbf{d}, \mathbf{\Omega}) = \frac{1}{m} \sum_{j=1}^{m} \left[ \log \det \left( \mathbf{\Lambda}_j^F(\mathbf{d}) \mathbf{\Omega}(\mathbf{d}) \mathbf{\Lambda}_j^F(\mathbf{d})^* \right) + \mathbf{W}^F(\lambda_j)^* \left( \mathbf{\Lambda}_j^F(\mathbf{d}) \mathbf{\Omega}(\mathbf{d}) \mathbf{\Lambda}_j^F(\mathbf{d})^* \right)^{-1} \mathbf{W}^F(\lambda_j) \right].$$

The solution satisfies

$$\hat{\mathbf{d}}^{MFW} = \underset{\mathbf{d}}{\operatorname{argmin}} \log \det(\hat{\mathbf{\Omega}}^{MFW}(\mathbf{d})) - 2\log(2) \left( \frac{1}{m} \sum_{j=1}^{m} \lambda_j \right), \quad (5)$$

$$\hat{\mathbf{\Omega}}^{MFW} = \hat{\mathbf{\Omega}}^{MFW}(\hat{\mathbf{d}}^{MFW}), \quad (6)$$

$$\text{with } \hat{\mathbf{\Omega}}^{MFW}(\mathbf{d}) = \frac{1}{m} \sum_{j=1}^{m} \text{Re} \left( \mathbf{\Lambda}_j^F(\mathbf{d})^{-1} \mathbf{I}^F(j) \mathbf{\Lambda}_j^F(\mathbf{d})^{-1} \right).$$
The dynamics of the frequencies at the neighbourhood of the origin is given by the dynamics of the spectral density around the zero frequency. The form of the criterion is justified by a second-order approximation of the spectral density matrix (1), rather than the approximation (2). For an estimation based on the first-order approximation (2), one should replace $\Lambda_j^F(\mathbf{d})$ by $\Lambda_j^F(1)(\mathbf{d}) = \text{diag}(\lambda_j \mathbf{d} e^{i\pi \mathbf{d}/2})$.

[18] established the theoretical performance of this estimation procedure, for both the long-range dependence parameters and the long-run covariance matrix. It is shown that the variance for the estimation of the vector $\mathbf{d}$ is decreased for the multivariate procedure with respect to a univariate one. It is worth mentioning that [10] developed a similar estimation procedure, based on a rougher approximation of the cross-spectral density, $\Lambda_j^F(\mathbf{d}) = \text{diag}(\lambda_j \mathbf{d})$. Interestingly, the quality of estimation for the vector $\mathbf{d}$ is similar. Nevertheless, the estimation of the long-run covariance matrix $\Omega$ is biased since it does not take into account the phase-shift appearing in $\Lambda_j^F(\mathbf{d})$. We refer to [10] and to [18] for a more detailed study of these estimators and their consistency.

## 2.3 Wavelet-based estimation (MWW)

Wavelets are providing a very efficient tool because of their high flexibility to deal with nonstationary time series which is particularly useful for real data applications. Their good performances in comparison to Fourier have already been shown for example in univariate settings [6].

Let $(\phi(\cdot), \psi(\cdot))$ be respectively a father and a mother wavelets, satisfying regularity conditions, as stated in [4].

At a given resolution $j \ge 0$, for $k \in \mathbb{Z}$, we define the dilated and translated functions $\phi_{j,k}(\cdot) = 2^{-j/2}\phi(2^{-j}\cdot -k)$ and $\psi_{j,k}(\cdot) = 2^{-j/2}\psi(2^{-j}\cdot -k)$. The wavelet coefficients of the process $\mathbf{X}$ are defined by

$$ \mathbf{W}_{j,k} = \int_{\mathbb{R}} \tilde{\mathbf{X}}(t) \psi_{j,k}(t) dt \quad j \ge 0, k \in \mathbb{Z}, $$

where $\tilde{\mathbf{X}}(t) = \sum_{k \in \mathbb{Z}} \mathbf{X}(k)\phi(t-k)$. For given $j \ge 0$ and $k \in \mathbb{Z}$, $\mathbf{W}_{j,k}$ is a $p$-dimensional vector $\mathbf{W}_{jk} = (W_{j,k}(1) \ W_{j,k}(2) \ \dots \ W_{j,k}(p))$ where $W_{j,k}(\ell) = \int_{\mathbb{R}} \tilde{\mathbf{X}}_{\ell}(t) \psi_{j,k}(t) dt$.

For any $j \ge 0$, the process $(\mathbf{W}_{j,k})_{k \in \mathbb{Z}}$ is covariance stationary [4]. Let $\theta_{\ell,m}(j)$ denote the wavelet covariance at scale $j$ between processes $X_\ell$ and $X_m$, i.e. $\theta_{\ell,m}(j) = \text{Cov}(W_{j,k}(\ell), W_{j,k}(m))$ for any position $k$. Let us introduce the wavelet scalogram

$$ \mathbf{I}^W(j) = \sum_{k \in \mathbb{Z}} \mathbf{W}_{j,k} \mathbf{W}_{j,k}^T. \qquad (7) $$

The wavelet scalogram is the equivalent of the Fourier periodogram. Yet the scalogram is not normalized, on the contrary of the periodogram. We also introduce the function $K(\cdot)$, defined as

$$ K(\delta) = \int_{-\infty}^{\infty} |\lambda|^{-\delta} |\hat{\psi}(\lambda)|^2 d\lambda, \quad \delta \in (-\alpha, M). \qquad (8) $$

The wavelet Whittle procedure is described in [4]. Let $\Lambda_j(\mathbf{d}) = \text{diag}(2^{j\mathbf{d}}, j_0 \le j \le j_1)$. Let $\mathbf{G}(\mathbf{d})$ denote a $p \times p$-matrix with $(\ell, m)$-th element equal to

$$ G_{\ell,m}(\mathbf{d}) = f^S(0) \Omega_{\ell,m} K(d_\ell + d_m) \cos(\pi(d_\ell - d_m)/2). \qquad (9) $$
The estimators $(\hat{\mathbf{d}}^{MWW}, \hat{\mathbf{G}}^{MWW})$ are defined by minimization of $\mathcal{L}^{MWM}(\mathbf{d}, \mathbf{G})$, with

$$ \mathcal{L}^{MWM}(\mathbf{d}, \mathbf{G}) = \frac{1}{n} \sum_{j=j_0}^{j_1} \left[ n_j \log \det \left( \Lambda_j^W(\mathbf{d}) \mathbf{G}(\mathbf{d}) \Lambda_j^W(\mathbf{d}) \right) + \sum_k \mathbf{W}_{j,k}^T \left( \Lambda_j^W(\mathbf{d}) \mathbf{G}(\mathbf{d}) \Lambda_j^W(\mathbf{d}) \right)^{-1} \mathbf{W}_{j,k} \right]. $$

The estimation is here based on a first-order approximation of the spectral density matrix around 0. The solutions of this problem satisfy

$$ \hat{\mathbf{d}}^{MWW} = \underset{\mathbf{d}}{\operatorname{argmin}} \log \det(\hat{\mathbf{G}}^{MWW}(\mathbf{d})) + 2 \log(2) \left( \frac{1}{n} \sum_{j=j_0}^{j_1} j n_j \right) \left( \sum_{\ell=1}^{p} d_{\ell} \right), \quad (10) $$

$$ \hat{\mathbf{G}}^{MWW}(\mathbf{d}) = \frac{1}{n} \sum_{j=j_0}^{j_1} \Lambda_j^W(\mathbf{d})^{-1} \mathbf{I}^W(j) \Lambda_j^W(\mathbf{d})^{-1}, \quad (11) $$

where $\mathbf{I}^W(j)$ is the wavelet scalogram at scale $j$ defined in (7). The long-run covariance matrix can then be estimated by

$$ \hat{\Omega}_{\ell,m}^{MWW} = \hat{G}_{\ell,m}^{MWW}(\hat{\mathbf{d}}^{MWW})/(\cos(\pi(\hat{d}_{\ell}^{MWW} - \hat{d}_{m}^{MWW})/2)K(\hat{d}_{\ell}^{MWW} + \hat{d}_{m}^{MWW})). \quad (12) $$

This second step in the estimation of the long-run covariance matrix $\Omega$ is needed because the wavelets used in this paper are real and cannot correct the phase-shift (given by (12)). This is not the case for the Fourier Whittle estimation described in [18]. Fortunately, the phase-shift can be expressed as a multiplicative cosine term in the covariance of the wavelet coefficients and a correction is still possible. [4] established that the MWW estimators (10) and (12) are consistent under non-restrictive conditions. The rate of convergence for the estimation of the long-range parameters $\mathbf{d}$ is similar to the MFW estimator and is minimax. We refer to [4] for the detailed study of the asymptotic behaviour of MWW estimation.

# 3 Examples of multivariate long-memory time series

This section is describing specific functions of **multiwave** for the user to be able to simulate multivariate long-memory processes. Parametric models are defined and implemented. In addition a data sets of real data from neuroimaging is provided.

## 3.1 Simulations of FIVARMA

**multiwave** package proposes simulation functions for time series with a spectral density satisfying approximation (2). The main function is **fivarma** which computes a parametric FIVARMA process defined in Section 2.1.

The input parameters of a FIVARMA($q, d, r$) process are the covariance matrix $\Sigma$ of the innovation process $\mathbf{u}$, the vector AR (AutoRegressive) ($\mathbf{A}_k)_{k=0,...,q}$, $\mathbf{A}_k \in \mathbb{R}^{p \times p}$, the vector MA (MovingAverage)
<figure><img src="image_3.png" /></figure>

## 3.2 A real data set

In order to describe how parameters can be chosen in a practical point of view, we provide a real data example (see Section 7).

Noninvasive data recorded from the brain are an example where the proposed methodology is efficient. The data consist of time series recording signals from the brain: electroencephalography (EEG) for the electrical signals, magnetoencephalography (MEG) for the magnetic signals or functional Magnetic Renonance Imaging (fMRI) for the Blood Oxygen Level Dependent (BOLD) signals. These data are intrinscally correlated because of the known interactions of the brain areas (also called regions of interest). Furthermore, it has already been shown that these time series present long-memory features [11]. Other data sets presenting similar features are coming from finance e.g. [20], where time series are correlated because of links between companies for example, and they also present long-memory characteristics. In this section, we observed time series extracted using fMRI facilities. The whole description of this data sets is detailed in [21]. The data set called *brainHCP* contains the time series of 1200 points in time and 89 regions of the brain. Figure 1 displays 6 arbitrary signals from one subject in this data set.

R> data(brainHCP)

R> dim(brainHCP)

[1] 1200 89
<figure><img src="image_4.png" /><figcaption>Figure 1: Plot of 6 arbitrary signals from a subject of fMRI data set.</figcaption></figure>

# 4 Estimation of the long-memory parameters and covariance matrix

The objective of this package is to provide the implementation of two sets of methods based either on Fourier or wavelet decomposition. That is Multivariate Fourier Whittle (MFW) and Multivariate Wavelet Whittle (MWW) estimation procedures. The corresponding functions are respectively called mfw and mww in the package.

The output of the implemented methods consists in the estimation of two quantities, d and Ω, where d corresponds to the long-memory parameters of the time series and Ω is reflecting the coupling between the pairs of time series.

The two functions mfw and mww are implementing the semiparametric Whittle estimation using respectively Fourier decomposition and wavelet decomposition in order to estimate d and Ω.

A fast execution of these two functions is

R> ## Fourier decomposition

R> m <- 57 ## default value of Shimotsu 2017

R> res_mfw <- mfw(x,m)

R> ## Wavelet decomposition

R> res_filter <- scaling_filter('Daubechies',8) ## choice of filter

R> filter <- res_filter$h

R> LU <- c(2,11) ## choice of wavelet scales

R> res_mww <- mww(x,filter,LU)
<figure><img src="image_5.png" /><figcaption>Figure 2: Input and Output of **multiwave** package in a two-dimensional case.</figcaption></figure>

## 4.1 Multivariate Fourier Whittle estimation

As a first implementation, it is natural to use Fourier decomposition to approximate the spectral density of time series.

Package **multiwave** proposes functions to compute MFW estimators:

* **mfw** computes the multivariate Fourier Whittle estimators of both the long-range dependence parameters and the long-run covariance matrix.

* **mfw_cov_eval** computes the multivariate Fourier-based Whittle estimator for the long-run covariance matrix for a given value of the long-range dependence **d**.

* **mfw_eval** returns the value of the multivariate Fourier Whittle criterion with respect to **d** at a given value of **d**.

The functions **mfw_cov_eval** and **mfw_eval** are internal functions of **mfw**. In **mfw**, we apply first a minimum search of **mfw_eval** with respect to **d**, and **mfw_cov_eval** is returning the estimation of **Ω** for the estimated value **d**.

We only detail function **mfw** hereafter and refer to the package description for other functions.

Let **X** be the $p \times N$-matrix of observations, with general term $x_{\ell,i} = X_{\ell}(i)$, $\ell = 1, \dots, p$ and $i = 1, \dots, N$. Let **m** be the number of frequencies used in MFW procedure. Given **x** and **m**, the function **mfw** computes the MFW estimators defined by (5) and (6), with the frequencies $\lambda_j = 2\pi j/N$, $j =$
1, ..., m. The optimization in equation (5) is done using `optimize` function of R in one-dimensional settings and a Newton-type algorithm through `nlm` function of R otherwise. The initialization of the algorithm is set equal to the vector of univariate Fourier-based Whittle estimations. Even if it increases the computational time, such an initialization is important in high-dimensional settings. For example, in the MEG data set studied in [4], the optimization is done in $\mathbb{R}^{274}$. An initialization at the origin may not be able to reach the minimum even with a high number of iterations, due to the high dimension.

Function `mfw` returns a list with first the $p$-dimensional vector $\hat{\mathbf{d}}^{MFW}$ and second the $p \times p$-matrix $\hat{\Omega}^{MFW}$. The quality of estimation is depending on the parameter $m$. Theoretical results show that $m$ must be small enough so the short-range properties of the time series do not bias the estimation. On the contrary a too small value will introduce variance in estimation since it decreases the number of frequency used in the procedure. [10], [18] or [13] propose a default value $m = N^{0.65}$. This choice is discussed in the simulation study in Section 6.

### Example.

```R
R> # Simulation of the data
R> N <- 2^8
R> d0 <- c(0.2,0.4)
R> rho <- 0.8
R> cov <- matrix(c(1,rho,rho,1),2,2)
R> VMA <- diag(c(0.4,0.7))
R> VAR <- array(c(0.8,0.2,0,0.6),dim=c(2,2))
R> resp <- fivarma(N, d0, cov_matrix=cov, VAR=VAR, VMA=VMA)
R> x <- resp$x
R>
R> # Estimation
R> m <- N^(0.65) # default value of Shimotsu
R> res_mfw <- mfw(x,m)
R> res_mfw
```

```R
$d
[1] 0.1239810 0.3497609
```

```R
$cov
[,1]      [,2]
[1,] 0.5004685 0.5840763
[2,] 0.5840763 1.2547056
```

## 4.2 Multivariate Wavelet Whittle estimation

The functions applying MWW estimation in package `multiwave` are the following:

* `mww` computes the multivariate wavelet Whittle estimators of the long-range dependence parameters and the long-run covariance matrix.
* `mww_cov_eval` computes the multivariate wavelet-based Whittle estimator for the long-run covariance matrix for a given value of the long-range dependence `d`.

* `mww_eval` returns the value of the multivariate wavelet Whittle criterion with respect to `d` at a given value of `d`.

`mww_cov_eval` and `mww_eval` are internal functions of `mww`. In `mww`, we apply first a minimum search of `mww_eval` with respect to `d`, and `mww_cov_eval` is returning the estimation of $\Omega$ for the estimated value of `d`. MWW estimation is based on the wavelet transform of time series and `mww` needs the definition of a wavelet filter. The computation of a filter and of a wavelet transform are described below.

### 4.2.1 Wavelet transform and scalogram

The wavelet decomposition in package `multiwave` is implemented using an exact discrete wavelet transform.

* `scalingfilter` defines the wavelet filter (only Daubechies' wavelets are available).

* `computenj` computes the number of wavelet coefficients for each individual scale.

* `DWTexact` provides the wavelet transform of the data.

* `psi_hat_exact` gives the Fourier transform of the wavelet function.

* `K_eval` evaluates the value of the integral (8).

*Example.* To obtain the wavelet filter of a Daubechies' wavelet of order 4, that is with 2 vanishing moments, one should write:

```
R> res_filter <- scaling_filter('Daubechies',4);
```

```
R> filter <- res_filter$h
```

```
R> filter
```

```
[1] 0.4829629 0.8365163 0.2241439 -0.1294095
```

Next, given an $N$ dimensional vector $x$, the wavelet coefficients of $x$ are given by function `DWTexact`:

```
R> # Simulation of the data
```

```
R> N <- 2^8
```

```
R> d0 <- 0.2
```

```
R> resp <- fivarma(N, d0)
```

```
R> x <- resp$x
```

```
R>
```

```
R> # Wavelet decomposition
```

```
R> resw <- DWTexact(x,filter)
```

```
R> xwav <- resw$dwt # returns the vector of the wavelet coefficients
```

```
R> index <- resw$indmaxband
```

```
R> Jmax <- resw$Jmax
```

```
R> Jmax
```
[1] 6

R> index

[1] 127 189 219 233 239 241

R> length(x)

[1] 256

index gives the index of the last coefficient at each scale and Jmax gives the maximal scale. The vector of coefficients xwav is m-dimensional, with m maximum of index, equal to index[Jmax]. The coefficients of the third scale are for example given by:

R> xwav Hand [1:3],1]

Finally it is useful to compute the quantity K(δ) defined in equation (8). Thus one needs to recover
the Fourier transform of the wavelet, ψ(·). This is done using the function psi_hat_exact. Its inputs
are the filter defined previously and an index of precision. It returns (ψ̂(uᵢ))ᵢ=1,...,q*2<sup>J</sup> where q is the
length of the filter and uᵢ are equally spaced points on the interval [-π2<sup>J-3</sup>(q-1)/2; π2<sup>J-3</sup>(q-1)/2].

R> res_psi <- psi_hat_exact(filter, J=10)

R> psih <- res_psi$psih

R> gridh <- res_psi$grid

where res_psi$grid returns the values of the grid (u_i) and res_psi$psih returns the corresponding
values of ψ̂(u_i). It is recommended to take J ≤ 15 in practice and the default value is J = 10. Indeed,
a large value of J is increasing the computational time.

Given the function ψ̂(·), we are now able to evaluate K(·) for a given value of d:

R> K <- K_eval(psih, gridh, d0)

4.2.2 Estimation

mww is now described in detail. Let X be the p × N-matrix of observations, with general entries
x_{\ell,i} = X_{\ell}(i), \ell = 1, \dots, p \text{ and } i = 1, \dots, N. Let LU be the bivariate vector giving the lowest scale j_0
and the upper scale j_1 of the wavelet coefficients used in estimation. Given X and LU, the function
mww computes the MWW estimators defined by (10) and (12). As previously, the optimization in
equation (10) is done using optimize function of R in one-dimensional settings and a Newton-type
algorithm through nlm function of R otherwise. The initialization of the algorithm is set equal to the
vector of univariate wavelet-based Whittle estimations. The reasons are identical to the ones given for
the function mfw.

Function mww returns a list with first the p-dimensional vector d̂^MWW and second the p × p-matrix
Ω^MWW. The quality of estimation is depending on the parameters j₀ and j₁ appearing in (10) and
(11). Default value for j₀ is set to 2 and for j₁ to the highest integer lower than log₂(N). The critical
value to choose for estimation is $j_0$, as it can be seen in the theoretical conditions for consistency [4] and in simulations studies. Similarly to the choice of the parameter $m$ for MFW procedure, a compromise exists between choosing a small value of $j_0$, which would introduce a bias due to the short-range properties of the time series, and a high value that will reduce the number of frequencies and thus increase variance.

### Example.

```
R> # Simulation of the data
R> N <- 2^8
R> d0 <- c(0.2,0.4)
R> rho <- 0.8
R> cov <- matrix(c(1,rho,rho,1),2,2)
R> VMA <- diag(c(0.4,0.7))
R> VAR <- array(c(0.8,0.2,0,0.6),dim=c(2,2))
R> resp <- fivarma(N, d0, cov_matrix=cov, VAR=VAR, VMA=VMA)
R> x <- resp$x
R>
R> # Parameter of estimation
R> res_filter <- scaling_filter('Daubechies',8);
R> filter <- res_filter$h
R> LU <- c(2,8)
R>
R> # Estimation
R> res_mww <- mww(x, filter, LU)
R> res_mww
```

```
$d
[1] 0.1955683 0.4857156
```

```
$cov
[,1]      [,2]
[1,] 0.5826133 0.5051330
[2,] 0.5051330 0.9584876
```

If one wants to apply several times the estimation on the same data set, or modifying the parameters of estimation, it is useful to separate the wavelet transform and the estimation scheme. Wavelet-based estimation can be evaluated directly on the wavelet transform of the data using the following functions:

* `mww_wav` computes the multivariate wavelet Whittle estimators of the long-range dependence parameters and the long-run covariance matrix, given the wavelet transform of the data.

* `mww_wav_cov_eval` computes the MWW estimator for the long-run covariance matrix for a given value of the long-range dependence $d$, given the wavelet transform of the data.

* `mww_wav_eval` returns the value of the multivariate wavelet-based Whittle criterion with respect to $d$ at a given value of $d$, for a specific wavelet transform of the data.
We refer to the description of the functions in the package for more details.

### Example.

R> # Simulation of the data

R> N <- 2^8

R> d0 <- c(0.2,0.4)

R> rho <- 0.8

R> cov <- matrix(c(1,rho,rho,1),2,2)

R> VMA <- diag(c(0.4,0.7))

R> VAR <- array(c(0.8,0.2,0,0.6),dim=c(2,2))

R> resp <- fivarma(N, d0, cov_matrix=cov, VAR=VAR, VMA=VMA)

R> x <- resp$x

R>

R> N <- dim(x)[1]

R> k <- dim(x)[2]

R>

R> # Parameter of estimation

R> res_filter <- scaling_filter('Daubechies',8)

R> filter <- res_filter$h

R> LU <- c(2,8)

R>

R> ## Wavelet decomposition

R> xwav <- matrix(0,N,k)

R> for(j in 1:k){

R> xx <- x[,j]

R> resw <- DWTextext(xx,filter)

R> xwav_temp <- resw$dwt

R> index <- resw$indmaxband

R> Jmax <- resw$Jmax

R> xwav[1:index[Jmax],j] <- xwav_temp

R> }

R> ## we free some memory

R> new_xwav <- matrix(0,min(index[Jmax],N),k)

R> if(index[Jmax]<N){

R> new_xwav[(1:(index[Jmax]))],] <- xwav[(1:(index[Jmax]))],]

R> }

R> xwav <- new_xwav

R> index <- c(0,index)

R>

R> ###### Compute the wavelet functions

R> res_psi <- psi_hat_exact(filter,Jmax)

R> psih <- res_psi$psih

R> grid <- res_psi$grid

R>

R> # Estimation

R> res_mww_wav <- mww_wav(xwav,index,psih,grid,LU)

R> res_mww_wav

$d
[1] 0.1684136 0.3693829

$cov
[,1]      [,2]
[1,] 0.5177418 0.426430
[2,] 0.4264300 0.834756

# 5 Practical choices of parameters for MWW estimation

MWW procedure implemented in mww depends on mainly two parameters: the choice of the wavelet bases filter and the choice of the wavelet scales LU. These parameters are not fixed in the package as the performances of the estimation may be improved by a careful choice. The possibility to choose the wavelet scales is particularly of interest when dealing with short-range dependence. We show that a simple graphical representation of the scalogram is able to guide the user in the choice of the wavelet scales.

## 5.1 Choice of the wavelet bases

Actually **multiwave** only proposes Daubechies' wavelets, which satisfy theoretical properties of [4]. Other bases are possible but not implemented. The wavelet bases is imputed *via* the parameter *filter* in function *mww*, where *filter* is obtained by *filter <- scaling_filter('Daubechies', 2*M)$h*. The main parameter characterizing the Daubechies' bases is thus the number of vanishing moments, *M*. MWW estimation presents the advantage to be available even if the time series are nonstationary or with polynomial trends, as soon as $max d \le M$. For example, for stationary time series, a parameter *M* = 1 is sufficient (which is equivalent to considering Haar bases). For real data application, when nonstationarity or trends are suspected, a higher value of *M* is necessary.

As discussed in [6], when *M* increases, the quality of estimation (slightly) decreases. Depending on the data, a compromise is then needed between choosing a large enough number of vanishing moments *M* to handle nonstationarity in the data and the quality of estimation.

On the contrary, MFW estimators are only suited to stationary time series. Some extensions of Fourier-based estimation were proposed in univariate setting such as tapered Fourier (see e.g. [6] and references therein). For multivariate estimation [13] proposes an extension of [18] based on the transform defined in [1]. However, this approach gives satisfactory results only for $d < 1.5$ and we decided not to implement it in the package for simplicity.

## 5.2 Choice of wavelet scales

The second parameter we need to tune is *LU*, which corresponds to the range of scales used in estimation. *LU* is a two-dimensional vector, that is $LU<-c(j0, j1)$, with *j0* the lowest scale and *j1* the upper scale. Parameters *j0* and *j1* are respectively *j0* and *j1* defined in equations (10) and (11) in the estimation procedure.

One advantage of wavelets is to be able to qualitatively evaluate the choice of wavelet scales to estimate the long-memory parameters and correlation by inspection of wavelet scalogram. As mentioned in [2, 6]
for univariate settings, the first and last scales may have to be discarded from the analysis. The first scale may be affected by the presence of short-memory phenomena. In the example of FIVARMA model, this is driven by the AR and MA coefficients. For the last scales, the impact is different and it comes from the finite length of the time series. Indeed, as derived in [24], the variance of the estimator is increasing with the wavelet scales.

The usual log-scalogram diagram used in univariate settings is showing the linear behaviour of the log variance with respect to the wavelet scales [2]. This is also true for the covariances as shown in Proposition 2 of [4]. For all $k \in \mathbb{Z}$, $Cov(W_{j,k}(\ell), W_{j,k}(m))$ is equivalent to $2^{j(d_\ell+d_m)} G_{\ell,m}(\mathbf{d})$ when $j$ goes to infinity, with $G_{\ell,m}(\mathbf{d})$ defined in equation (9). This property is illustrated in Figure 3 with a bivariate FIVARMA processes, as described in Section 2.1. This figure represents the boxplots of the variance of the wavelet coefficients of each component of the time series at each scales (subfigures 3.(a) and 3.(b)). These plots correspond to the usual log-scalogram diagram [2]. Subfigure 3.(c) displays the analog representation of the covariance between the components. For both variance and covariance, the points satisfying the above approximation are aligned. Scales corresponding to non-aligned points should be removed from estimation as it can be seen on subfigure 3.(a) and 3.(c) where the highest frequencies are modified by the presence of short-range dependence. Thus, one may discard the first scale from estimation to improve its quality.

Using the results on the wavelet variance and covariance in terms of scales, we showed that the wavelet correlation is asymptotically constant with respect to the wavelet scales. Indeed, for all $\ell, m = 1, \dots, p$, for all $k \in \mathbb{Z}$, $Cor(W_{j,k}(\ell), W_{j,k}(m))$ is equivalent to $G_{\ell,m}(\mathbf{d})/\sqrt{G_{\ell,\ell}(\mathbf{d})G_{m,m}(\mathbf{d})}$ when $j$ goes to infinity [4]. As for the log-scalogram diagram, the correlation between wavelet coefficients with respect to the scales can be plotted and scales where the observed correlation is not equal to the value obtained for the majority of scales should be removed from estimation. The wavelet correlation spectrum is a complementary way to qualitatively evaluate the range of scales where the analysis should be carried out. Figure 4 illustrates this on four different data sets. Four different simulations of bivariate processes are applied using finite difference processes, FIVAR and FIVARMA processes, as described in Section 2.1. Figure 4 represents the boxplots of the correlation between wavelet coefficients of the two components of the time series at each scales. Again the presence of short-range dependence alters the highest frequencies (subfigures 4.(b) and 4.(c)). This is also observed with nonstationarity (subfigure 4.(d)). The first scales should then be removed from estimation. This free parameter of the package is particularly useful with the presence of short-range dependence or nonstationarity. Visual comparison to constant values may be easier for selection of the correct range of wavelet scale to use in the estimation.

A similar discussion is detailed in Section 7 for a real neuroscience data set. With real data sets a bootstrap procedure is necessary to obtain boxplots, as it will be explained in Section 7. For the Fourier procedure, the equivalent parameter is the number of frequencies $m$. However, wavelets are providing a graphical way to choose the upper and lower scales. To our knowledge, no equivalent qualitative evaluation for Fourier procedure is available.

## 5.3 Numerical examples

In order to quantify the quality of the choice of the parameters, numerical examples are provided. The first tables illustrate the quality of the estimators for a long-memory process with no short-range dependence. Then, the simulations are complexified by adding short-range behaviour or nonstationarity.

In each example we simulated 1000 Monte-Carlo replications of $N = 512$ observations of FIVARMA(q,d,r),
Figure 3: Boxplots of the log2-variance of the wavelet coefficients at different scales for simulated bivariate FIVARMA(1, (0.4, 0.4), 1); (a) for the first component and (b) for the second component. (c) Boxplots of the log2-absolute covariance for the same data. The red lines represent the theoretical linear prediction given by $j(d_l + d_m) + \log_2(G_{l,m}(\mathbf{d}))$, with $l = m = 1$ for subfigure (a), $l = m = 2$ for subfigure (b) and $l = 1, m = 2$ for subfigure (c). The horizontal axis corresponds to increasing scales, that is, decreasing frequencies. The indexes of the horizontal axis display the number of coefficients available. Parameters in FIVARMA model were the following: the white noise is Gaussian with a covariance equal to $\Sigma = \begin{pmatrix} 1 & 0.8 \\ 0.8 & 1 \end{pmatrix}$, the AR coefficient is set equal to $A = \begin{pmatrix} 0.8 & 0 \\ 0.2 & 0.6 \end{pmatrix}$ and the MA coefficient is set equal to $B = \begin{pmatrix} 0.4 & 0 \\ 0.2 & 0.7 \end{pmatrix}$. Calculation was done on $N = 512$ observations for 100 replications.

<figure><img src="image_6.png" /></figure>

with a dimension $p = 2$, for a set of different parameter values. Simulations are done using the function *fivarma*. MWW estimators are computed using function *mww*. MWW procedure is done using a Daubechies' wavelet bases with $M = 4$ vanishing moments. This choice is motivated by discussion of Section 5.1, because it can handle different settings, including nonstationary ones.

The quality of estimation is measured *via* the bias, the standard deviation (std) and the RMSE which is equal to $\sqrt{bias^2 + std^2}$. For clarity, all tables are displayed at the end of the paper.

### 5.3.1 A reference example

We first consider a simple example, with neither short-range dependence, nor nonstationarity. Time series were simulated using a bivariate FIVARMA(0, $\mathbf{d}$, 0) with a long-run correlation matrix $\Omega =$
<figure><img src="image_7.png" /><figcaption>Figure 4: Boxplots of the correlation of the wavelet coefficients at different scales for four different simulated data: (a) bivariate FIVARMA(0, (0.4, 0.4), 0) (b) bivariate FIVARMA(1, (0.4, 0.4), 0); (c) bivariate FIVARMA(1, (0.4, 0.4), 1); and finally, (d) an example with non stationary time series, FIVARMA(0, (0.8, 1.2), 0). The horizontal red lines represent the true long-run correlation for each simulation. The horizontal axis corresponds to increasing scales, that is, decreasing frequencies. The indexes of the horizontal axis display the number of coefficients available. Parameters in FIVARMA models were, if needed, the following: the white noise covariance is set equal to $\Sigma = \begin{pmatrix} 1 & 0.8 \\ 0.8 & 1 \end{pmatrix}$, the AR coefficient is set equal to $A = \begin{pmatrix} 0.8 & 0 \\ 0.2 & 0.6 \end{pmatrix}$ and the MA coefficient is set equal to $B = \begin{pmatrix} 0.4 & 0 \\ 0.2 & 0.7 \end{pmatrix}$. Calculation was done on $N = 512$ observations for 100 replications.</figcaption></figure>

$\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$ and $\rho = 0.8$. The bivariate vector $\mathbf{d}$ is chosen in $[0, 0.5)^2$, such that the time series are stationary.

As shown in Figure 4(a), and discussed in Section 5.2, all scales can be kept for estimation.

Table 2 displays results for the MWW estimation of $\mathbf{d}$. This illustrates that multivariate estimation improves the quality of estimation for $\mathbf{d}$. Indeed, the last column gives the ratio between the RMSEs of
the multivariate wavelet-based estimation and of the univariate wavelet-based estimation (ratio M/U). This ratio is always smaller than 1, that is, multivariate RMSE is always lower than univariate RMSE.

### 5.3.2 Short-range dependence

Consider a FIVARMA(1, d, 0) obtained with the model described in section 2.1 and given by the function fivarma. This case corresponds to a FIVAR model of [17]. The AR coefficient is taken equal to $A = \begin{pmatrix} 0.8 & 0 \\ 0.2 & 0.6 \end{pmatrix}$ and the correlation between the innovation processes equal to $\rho = 0.8$. More precisely let $\varepsilon$ be a bivariate white noise process with covariance matrix $\Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$ and let $u$ be the AR(1) process defined by $u(t) + Au(t-1) = \varepsilon(t)$. The time series observation $X(t)$ at time $t$ satisfies $(1 - L)^d X(t) = u(t)$. The matrix $\Omega$ in equation (4) is equal to

$$ \Omega = (I + A)^{-1} \Sigma (I + A)^{-1} T \approx \begin{pmatrix} 0.3086 & 0.2392 \\ 0.2392 & 0.3260 \end{pmatrix}. $$

The corresponding long-run correlation is thus equal to 0.754.

As explained above, the finest scales are influenced by the short-range dependence and they have to be discarded from the estimation. As it can be seen in Figure 4(c), the first two scales should be removed. We obtained accordingly that the lowest RMSE in estimation is obtained taking $j_0 = 3$.

Numerical results for estimation are given in Table 4 and Table 5. We can observe that estimation of $\hat{d}$ and $\hat{\Omega}$ are still satisfactory. The RMSE is very similar to the previous case with no short-range behaviour.

### 5.3.3 Nonstationarity

Nonstationary examples are simulated using values of d higher than 0.5. We consider order 1 or 2 of nonstationarity, that is $d \in [0.5, 2.5)^2$.

The behaviour of the wavelet correlations at each scale is illustrated in Figure 4(d). Contrary to stationary simulations where the optimal choice of $j_0$ was equal to $j_0 = 1$, the optimal choice of the parameter $j_0$ is $j_0 = 2$. Results are given in Table 6 and Table 7.

Comparing Table 2 and Table 6, the quality of estimation of d is still accurate in nonstationary settings, with similar values for the RMSE. As for the estimation of $\Omega$, Table 3 and Table 7 indicate that MWW still provides a good quality estimation of the long-run covariance matrix. The quality is slightly lower but still satisfactory.

## 5.4 Discussion on identifiability

In practical applications, it seems natural to assume that time series have the same order of stationarity. However, when two time series have long-memory parameters $d_\ell$ and $d_m$ satisfying $d_\ell - d_m = 1$ the long-run covariance matrix $\Omega$ is no longer identifiable with the wavelet-based procedure. Indeed, Proposition 2 in [4] states that in this particular case the covariance $Cov(W_{j,k}(\ell), W_{j,k}(m))$ tends to 0 when the scale $j$ tends to infinity. Figure 5 illustrates this approximation for a bivariate

$$ FIVARMA(0, (0.2, 1.2), 0), \text{ a correlation matrix } \Omega = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix} \text{ and } \rho = 0.8. $$
<figure><img src="image_8.png" /><figcaption>Figure 5: Boxplots of the covariance of the wavelet coefficients at different scales for a bivariate ARFIMA(0, (0.2, 1.2), 0) with $\Omega = \begin{pmatrix} 1 & 0.8 \\ 0.8 & 1 \end{pmatrix}$. The index of the horizontal axis displays the number of coefficients available. Calculation was done on $N = 512$ observations for 1000 replications.</figcaption></figure>

When $\hat{d}_{\ell} - \hat{d}_{m} = 1$, the estimator (12) is no longer defined. In practice, the quantity $\hat{d}_{\ell} - \hat{d}_{m}$ cannot be exactly equal to 1. Nevertheless, as dividing by a cosine function of this difference, a small error in the estimation of $(d_{\ell}, d_{m})$ will lead to an important bias in the estimation of $\Omega_{\ell,m}$. As it can be seen in Figure 6, the resulting bias increases in the neighbourhood of the non-identifiable lines $d_{\ell} - d_{m} = \pm 1$. When this situation occurs, say when the difference between $d_{\ell} - d_{m}$ is between 0.75 and 1.25, the estimation of $\mathbf{d}$ is not affected. But the user must be careful for the estimation of $\Omega$. One solution is to differentiate or integrate one of the two processes. For example, Table 1 illustrates the non-identifiability of $\Omega$ in a bivariate $FIVARMA(0, (0.2 \ 1.2), 0)$. When differentiating the second component (with $d_{2} = 1.2$) the estimator has again good performances.

# 6 MFW estimation and comparison with MWW

The comparison between Fourier-based and wavelet-based approach is presented now. Time series were simulated using a bivariate $FIVARMA(0, \mathbf{d}, 0)$ with a long-run correlation matrix $\Omega = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$ and $\rho = 0.8$. The bivariate vector $\mathbf{d}$ is chosen in $[0, 0.5)^2$, such that the time series are stationary. In such a setting, Fourier-based estimators are available.

For Fourier-based approach, the parameter to choose is $m$ corresponding to the number of frequencies taken into account in the estimation. The default value in $[18]$ is $m = N^{0.65}$. We also make comparisons with an optimal value computed by minimizing the RMSE.
Figure 6: RMSE in the estimation of the cross-covariance term $\Omega_{12}$ with respect to $(d_1, d_2)$. Estimation was done using multivariate Wavelet Whittle estimator in a bivariate FIVARMA(0, (d_1, d_2), 0) with $\Omega = \begin{pmatrix} 1 & 0.8 \\ 0.8 & 1 \end{pmatrix}$. Subfigure (b) represents an image plot of subfigure (a) (with a different colour scale to improve visual quality). Blue lines on subfigure (b) correspond to $d_2 - d_1 = \pm 1$. Calculation was done on $N = 512$ observations for 1000 replications.

(a)

<figure><img src="image_9.png" /></figure>

(b)

<figure><img src="image_10.png" /></figure>

## 6.1 Estimation of the long-memory parameters

Table 8 gives the results obtained for the estimation of $d$ using MFW procedure (to be compared to Table 2 for wavelets). Comparison between Fourier and wavelet-based procedures is summarized by the ratio between the RMSE given by MWW estimation and the RMSE given by MFW estimation, denoted by ratio W/F. Taking the same number of frequencies as [18], that is, $m = N^{0.65}$, Table 8 shows that the quality of MWW end MFW procedures are comparable, even if wavelet-based estimation slightly improves Fourier-based estimation with such a choice of $m$.

Next we also consider the number of frequencies leading to the minimal RMSE for MFW estimation. As it can be seen in Table 8, qualities of both procedures are very similar but MFW estimation then (slightly) surpasses MWW estimation.

Very precise comparisons of Fourier-based and wavelet-based approaches are described in [6] for a univariate setting. In particular, it is shown that since the time series are stationary, the use of Haar bases should improve MWW quality. The authors indeed obtained better results with the Haar-based
<table><caption>Table 1: Multivariate Wavelet Whittle estimation of $\Omega$ for a bivariate FIVARMA(0, (0.2, 1.2), 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. As the memory parameter of the second component is greater than one, one possibility is to differentiate the second component. $j_0$ is chosen to be equal to 1.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">Without differentiation</th><th colspan="3">With differentiation</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>0.0935</td><td>0.0762</td><td>0.1206</td><td>0.0349</td><td>0.0678</td><td>0.0762</td></tr><tr><td>Ω<sub>1,2</sub></td><td>4.1863</td><td>7.2103</td><td>8.3375</td><td>0.0261</td><td>0.06</td><td>0.0654</td></tr><tr><td>Ω<sub>2,2</sub></td><td>0.2215</td><td>0.0819</td><td>0.2362</td><td>0.0292</td><td>0.07</td><td>0.0758</td></tr><tr><td>correlation</td><td>3.5255</td><td>6.2935</td><td>7.2137</td><td>0.0003</td><td>0.0155</td><td>0.0155</td></tr></tbody></table>

procedure than with Fourier-based procedure. To highlight the versatility of wavelet-based procedures, we choose here wavelet bases with four vanishing moments. A fair comparison with a Fourier-based method should consider tapered Fourier of order 4 as it is detailed in [6]. They quantify the influence of the regularity of the wavelet bases and discuss the comparison with (tapered) Fourier bases. Similar results are expected to be obtained in the case of multivariate time series, however, this topic exceeds the scope of this paper.

## 6.2 Estimation of the long-run covariance

Finally, Table 9 and Table 10 display results for the estimation of $\Omega$ with MFW method (to be compared to Table 3 for wavelets). When MFW estimation is applied with the usual number $m = N^{0.65}$ of frequencies, one can see that the wavelet-based procedure still estimates better the long-run covariance and the long-run correlation, with a ratio W/F always lower than 0.6 for the estimation of $\Omega$ terms. When the number of frequencies in Fourier-based estimation is chosen optimally, MFW and MWW procedures behave similarly and none appears significantly better than the other.

To conclude, MFW and MWW estimation procedures give very similar results. The slight improvement of Fourier-based procedure for the estimation of $d$ can be explained by the choice of the wavelet bases, however wavelets are efficient for a large set of applications, including time series with trends and nonstationarity features.

## 7 Application on real neuroscience data

As already shown in Figure 4 for simulated data, the advantage of representing the wavelet correlation in terms of scale is to qualitatively determine the scales necessary to estimate the long-memory parameters and long-range covariance matrix. When dealing with real data, bootstrap is providing a way to assess the variability of the estimators. Using the real data described in Section 3.2, sliding overlapping window of the time series were extracted containing 512 points and we repeated the estimation until reaching the final point of the time series. This is illustrated in Figure 7, where an example of four pairs of fMRI data from one subject is presented. Boxplots are constructed using
the sliding window extractions. From these plots, and taking into account neuroscientific hypothesis stating that the signal of interest for resting state is occurring for frequency below 0.1Hz, we chose to compute the long-memory parameters between scales 3 and 6.

Figure 7: Boxplots of the correlation of the wavelet coefficients at different scales for real time series from fMRI data sets: (a) Time series 1 and 2; (b) Time series 13 and 14; (c) Time series 31 and 32; (d) Time series 47 and 48. boxplots were obtained using sliding windows with *N* = 512 points, extracted from two fMRI time series with length equal to 1200 points, from a single subject. The estimated long parameters *d* of the two time series are equal. fMRI data set is described in section 7. The index of the horizontal axis displays the number of coefficients available. The horizontal red lines represent the estimated long-run correlation. Calculation was done on *N* = 512 observations for 100 replications using sliding windows (with overlap).

<figure><img src="image_11.png" /></figure>

Figures 8 and 9 display an example of long-memory parameter and long-run correlation estimated for one subject.
<figure><img src="image_12.png" /><figcaption>Figure 9: Estimation of $\Omega$ from a subject of fMRI data set.</figcaption></figure>

<figure><img src="image_13.png" /></figure>
Conclusion

The R package **multiwave** provides a versatile wavelet-based approach, as well as a Fourier-based approach, for estimating long-memory parameters and long-run covariance matrices of multivariate time series. The two estimation procedures are based on semi-parametric approaches of [18] and [4]. The added value of the package is to provide estimations in long-range dependence multidimensional settings, which is not proposed presently by any R package to our knowledge. This paper describes the functions of the package **multiwave** and discusses some practical points for applications, including on a real data set. A simulation study shows first that multivariate estimation improves univariate estimation. The advantage of the wavelet-based procedure with respect to the Fourier-based estimation is its flexibility, allowing to take into account trends or nonstationarity.

**Acknowledgements.** The authors are grateful to Shimotsu (http://shimotsu.web.fc2.com/Site/Matlab_Codes.html) and to Faÿ, Moulines, Roueff and Taqqu for kindly providing the codes of their respective papers. This work was partly supported by the project Graphsip from Agence Nationale de la Recherche (ANR-14-CE27-0001). S.A. was partly funded by a grant from la RÃªgion RhÃtne-Alpes and a grant from AGIR-PEPS, UniversitÃl Grenoble AlpesaÃŞCNRS.

References

[1] Karim M Abadir, Walter Distaso, and Liudas Giraitis. Nonstationarity-extended local Whittle estimation. *Journal of Econometrics*, 141(2):1353–1384, 2007.

[2] Patrice Abry and Darryl Veitch. Wavelet analysis of long-range-dependent traffic. *Information Theory, IEEE Transactions on*, 44(1):2–15, 1998.

[3] Sophie Achard, Danielle S Bassett, Andreas Meyer-Lindenberg, and Ed Bullmore. Fractal connectivity of long-memory networks. *Physical Review E*, 77(3):036104, 2008.

[4] Sophie Achard and Irene Gannaz. Multivariate wavelet whittle estimation in long-range dependence. *Journal of Time Series Analysis*, 37:476–512, 2016.

[5] J. Beran. *Statistics for long-memory processes*. New York: Chapman and Hall, 1994.

[6] Gilles Faÿ, Eric Moulines, François Roueff, and Murad S Taqqu. Estimators of long-memory: Fourier versus wavelets. *Journal of Econometrics*, 151(2):159–177, 2009.

[7] John Haslett and Adrian E Raftery. Space-time modelling with long-memory dependence: Assessing Ireland's wind power resource. *Applied Statistics*, 38(1):1–50, 1989.

[8] Rob J Hyndman and Yeasmin Khandakar. Automatic time series forecasting: the forecast package for R. *Journal of Statistical Software*, 26(3):1–22, 2008.

[9] Ignacio N Lobato. Consistency of the averaged cross-periodogram in long memory series. *Journal of Time Series Analysis*, 18(2):137–155, 1997.

[10] Ignacio N Lobato. A semiparametric two-step estimator in a multivariate long memory model. *Journal of Econometrics*, 90(1):129–153, 1999.
[11] Voichiţa Maxim, Levent Şendur, Jalal Fadili, John Suckling, Rebecca Gould, Rob Howard, and Ed Bullmore. Fractional Gaussian noise, functional MRI and Alzheimer's disease. *Neuroimage*, 25(1):141-158, 2005.

[12] Eric Moulines, François Roueff, and Murad S Taqqu. On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter. *Journal of Time Series Analysis*, 28(2):155-187, 2007.

[13] Frank S Nielsen. Local Whittle estimation of multi-variate fractionally integrated processes. *Journal of Time Series Analysis*, 32(3):317-335, 2011.

[14] Peter M Robinson. Robust covariance matrix estimation: HAC estimates with long memory/an-tipersistence correction. *Econometric Theory*, 21(01):171-180, 2005.

[15] Peter M Robinson. Multiple local Whittle estimation in stationary systems. *The Annals of Statistics*, 36(05):2508-2530, 2008.

[16] R. J. Sela and C. M. Hurvich. Computationaly efficient methods for two multivariate fractionnaly integrated models. *Journal of Time Series Analysis*, 30:6, 2008.

[17] Rebecca J Sela and Clifford M Hurvich. The averaged periodogram estimator for a power law in coherency. *Journal of Time Series Analysis*, 33(2):340-363, 2012.

[18] Katsumi Shimotsu. Gaussian semiparametric estimation of multivariate fractionally integrated processes. *Journal of Econometrics*, 137(2):277-310, 2007.

[19] Katsumi Shimotsu. Exact local whittle estimation of fractionally cointegrated systems. *Journal of Econometrics*, 169:266-278, 2012.

[20] Jitkomut Songsiri and Lieven Vandenberghe. Topology selection in graphical models of autoregressive processes. *The Journal of Machine Learning Research*, 11:2671-2705, 2010.

[21] M. Termenon, A. Jaillard, C. Delon-Martin, and S. Achard. Reliability of graph analysis of resting state fmri using test-retest dataset from the human connectome project. *NeuroImage*, pages -, 2016.

[22] Ruey S Tsay. *Multivariate Time Series Analysis: with R and financial applications*. John Wiley & Sons, 2013.

[23] Justin Q. Veenstra. *Persistence and Anti-persistence: Theory and Software*. PhD thesis, Western University, 2012.

[24] B. Whitcher, P. Guttorp, and D. B. Percival. Wavelet analysis of covariance with application to atmospheric time series. *Journal of Geophysical Research*, 105(D11)(14):941-962, 2000.
<table><caption>Table 2: Multivariate wavelet Whittle estimation of d for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. For the estimation, $j_0 = 1$.</caption><thead><tr><th>d</th><th>bias</th><th>std</th><th>RMSE</th><th>ratio M/U</th></tr></thead><tbody><tr><td>0.2</td><td>-0.017</td><td>0.039</td><td>0.0425</td><td>0.7785</td></tr><tr><td>0.0</td><td>0.013</td><td>0.0391</td><td>0.0412</td><td>0.9014</td></tr><tr><td>0.2</td><td>-0.0313</td><td>0.0376</td><td>0.049</td><td>0.8960</td></tr><tr><td>0.2</td><td>-0.0316</td><td>0.0378</td><td>0.0493</td><td>0.8805</td></tr><tr><td>0.2</td><td>-0.017</td><td>0.0383</td><td>0.0419</td><td>0.7673</td></tr><tr><td>0.4</td><td>-0.0442</td><td>0.0395</td><td>0.0592</td><td>0.7902</td></tr></tbody></table>

<table><caption>Table 3: Wavelet Whittle estimation of $\Omega$ for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. For the estimation, $j_0 = 1$.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (0.2, 0)</th><th colspan="3">d = (0.2, 0.2)</th><th colspan="3">d = (0.2, 0.4)</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>0.0417</td><td>0.0724</td><td>0.0836</td><td>0.0343</td><td>0.0711</td><td>0.0789</td><td>0.0417</td><td>0.0717</td><td>0.083</td></tr><tr><td>Ω<sub>1,2</sub></td><td>0.0382</td><td>0.0657</td><td>0.0759</td><td>0.0279</td><td>0.0626</td><td>0.0686</td><td>0.0673</td><td>0.0684</td><td>0.0959</td></tr><tr><td>Ω<sub>2,2</sub></td><td>0.0048</td><td>0.0709</td><td>0.071</td><td>0.0323</td><td>0.0714</td><td>0.0784</td><td>0.0748</td><td>0.0748</td><td>0.1057</td></tr><tr><td>correlation</td><td>0.0191</td><td>0.0227</td><td>0.0296</td><td>0.001</td><td>0.0164</td><td>0.0164</td><td>0.0194</td><td>0.0234</td><td>0.0304</td></tr></tbody></table>

<table><caption>Table 4: Multivariate wavelet Whittle estimation of d for a bivariate ARFIMA(1, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. For the estimation, $j_0 = 3$.</caption><thead><tr><th>d</th><th>bias</th><th>std</th><th>RMSE</th><th>ratio M/U</th></tr></thead><tbody><tr><td>0.2</td><td>-0.0473</td><td>0.1213</td><td>0.1302</td><td>0.8472</td></tr><tr><td>0.0</td><td>-0.0371</td><td>0.1266</td><td>0.132</td><td>0.8511</td></tr><tr><td>0.2</td><td>-0.0623</td><td>0.1209</td><td>0.136</td><td>0.8848</td></tr><tr><td>0.2</td><td>-0.0526</td><td>0.1258</td><td>0.1364</td><td>0.8714</td></tr><tr><td>0.2</td><td>-0.066</td><td>0.1244</td><td>0.1408</td><td>0.9161</td></tr><tr><td>0.4</td><td>-0.0584</td><td>0.1293</td><td>0.1418</td><td>0.8935</td></tr></tbody></table>
<table><caption>Table 5: Wavelet Whittle estimation of $\Omega$ for a bivariate ARFIMA(1, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. For the estimation, $j_0 = 3$.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (0.2, 0)</th><th colspan="3">d = (0.2, 0.2)</th><th colspan="3">d = (0.2, 0.4)</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>-0.0067</td><td>0.0828</td><td>0.0831</td><td>0.0024</td><td>0.0859</td><td>0.0860</td><td>0.005</td><td>0.0889</td><td>0.089</td></tr><tr><td>Ω<sub>1,2</sub></td><td>0.0551</td><td>0.0828</td><td>0.0995</td><td>0.0495</td><td>0.0794</td><td>0.0936</td><td>0.0507</td><td>0.0875</td><td>0.1012</td></tr><tr><td>Ω<sub>2,2</sub></td><td>0.1363</td><td>0.131</td><td>0.1891</td><td>0.1412</td><td>0.1384</td><td>0.1977</td><td>0.1391</td><td>0.1426</td><td>0.1992</td></tr><tr><td>correlation</td><td>-0.0088</td><td>0.0805</td><td>0.081</td><td>-0.0386</td><td>0.0527</td><td>0.0653</td><td>-0.0358</td><td>0.0985</td><td>0.1047</td></tr></tbody></table>

<table><caption>Table 6: Multivariate wavelet Whittle estimation of d for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. Nonstationary cases. For the estimation, $j_0 = 2$.</caption><thead><tr><th>d</th><th>bias</th><th>std</th><th>RMSE</th><th>ratio M/U</th></tr></thead><tbody><tr><td>1.2</td><td>-0.0338</td><td>0.0762</td><td>0.0834</td><td>0.8510</td></tr><tr><td>1</td><td>-0.0276</td><td>0.0725</td><td>0.0776</td><td>0.8316</td></tr><tr><td>1.2</td><td>-0.043</td><td>0.0732</td><td>0.0849</td><td>0.8672</td></tr><tr><td>1.2</td><td>-0.0411</td><td>0.0743</td><td>0.0849</td><td>0.8591</td></tr><tr><td>1.2</td><td>-0.0338</td><td>0.0741</td><td>0.0814</td><td>0.8310</td></tr><tr><td>1.4</td><td>-0.0356</td><td>0.0797</td><td>0.0873</td><td>0.8344</td></tr></tbody></table>

<table><thead><tr><th>d</th><th>bias</th><th>std</th><th>RMSE</th><th>ratio M/U</th></tr></thead><tbody><tr><td>2.2</td><td>-0.0421</td><td>0.0884</td><td>0.0979</td><td>0.8718</td></tr><tr><td>2</td><td>-0.0403</td><td>0.0862</td><td>0.0951</td><td>0.8516</td></tr><tr><td>2.2</td><td>-0.0503</td><td>0.086</td><td>0.0996</td><td>0.8874</td></tr><tr><td>2.2</td><td>-0.049</td><td>0.0823</td><td>0.0958</td><td>0.8566</td></tr><tr><td>2.2</td><td>-0.0436</td><td>0.0868</td><td>0.0971</td><td>0.8651</td></tr><tr><td>2.4</td><td>-0.0429</td><td>0.0831</td><td>0.0935</td><td>0.8400</td></tr></tbody></table>
<table><caption>Table 7: Wavelet Whittle estimation of $\Omega$ for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. Nonstationary cases. For the estimation, $j_0 = 2$.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (1.2, 1)</th><th colspan="3">d = (1.2, 1.2)</th><th colspan="3">d = (1.2, 1.4)</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>-0.0049</td><td>0.1362</td><td>0.1363</td><td>0.0059</td><td>0.1369</td><td>0.137</td><td>-0.0047</td><td>0.1361</td><td>0.1361</td></tr><tr><td>Ω<sub>1,2</sub></td><td>0.018</td><td>0.1168</td><td>0.1182</td><td>0.0093</td><td>0.1155</td><td>0.1158</td><td>0.0159</td><td>0.1266</td><td>0.1276</td></tr><tr><td>Ω<sub>2,2</sub></td><td>-0.0027</td><td>0.1277</td><td>0.1277</td><td>0.0042</td><td>0.1386</td><td>0.1386</td><td>-0.0113</td><td>0.1487</td><td>0.1491</td></tr><tr><td>correlation</td><td>0.0214</td><td>0.0475</td><td>0.0521</td><td>0.0051</td><td>0.0286</td><td>0.0291</td><td>0.0227</td><td>0.0507</td><td>0.0555</td></tr></tbody></table>

<table><caption>Table 8: Multivariate Fourier Whittle estimation of d for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions depending on the number of frequencies $m$. $\lfloor x \rfloor$ denotes the closest integer smaller than $x$.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (2.2, 2)</th><th colspan="3">d = (2.2, 2.2)</th><th colspan="3">d = (2.2, 2.4)</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th><th>bias</th><th>std</th><th>RMSE</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>-0.0383</td><td>0.1795</td><td>0.1835</td><td>-0.0253</td><td>0.1789</td><td>0.1807</td><td>-0.0361</td><td>0.1776</td><td>0.1812</td></tr><tr><td>Ω<sub>1,2</sub></td><td>-0.0043</td><td>0.1565</td><td>0.1565</td><td>-0.0129</td><td>0.1493</td><td>0.1498</td><td>-0.0097</td><td>0.1602</td><td>0.1605</td></tr><tr><td>Ω<sub>2,2</sub></td><td>-0.0318</td><td>0.1776</td><td>0.1804</td><td>-0.0276</td><td>0.1809</td><td>0.183</td><td>-0.0481</td><td>0.1813</td><td>0.1876</td></tr><tr><td>correlation</td><td>0.0251</td><td>0.0604</td><td>0.0654</td><td>0.0087</td><td>0.0374</td><td>0.0384</td><td>0.0249</td><td>0.0626</td><td>0.0674</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">d</th><th colspan="4">m = ⌊N<sup>0.65</sup>⌋ = 57</th><th colspan="4">m = ⌊N<sup>η</sup>⌋</th></tr><tr><th>bias</th><th>std</th><th>RMSE</th><th>W/F</th><th>η</th><th>bias</th><th>std</th><th>RMSE</th><th>W/F</th></tr></thead><tbody><tr><td>0.2</td><td>-0.002</td><td>0.0576</td><td>0.047</td><td>0.9049</td><td>0.90</td><td>-0.0197</td><td>0.0271</td><td>0.0335</td><td>1.2689</td></tr><tr><td>0</td><td>-9 × 10<sup>-4</sup></td><td>0.0593</td><td>0.0238</td><td>1.7295</td><td></td><td>-0.0033</td><td>0.0264</td><td>0.0267</td><td>1.5475</td></tr><tr><td>0.2</td><td>-0.0033</td><td>0.0574</td><td>0.0531</td><td>0.9214</td><td>0.85</td><td>-0.013</td><td>0.0306</td><td>0.0332</td><td>1.4750</td></tr><tr><td>0.2</td><td>-0.0031</td><td>0.0591</td><td>0.0522</td><td>0.9434</td><td></td><td>-0.0123</td><td>0.0293</td><td>0.0318</td><td>1.5503</td></tr><tr><td>0.2</td><td>8 × 10<sup>-4</sup></td><td>0.0576</td><td>0.0579</td><td>0.7239</td><td>0.85</td><td>-0.0136</td><td>0.0308</td><td>0.0337</td><td>1.2455</td></tr><tr><td>0.4</td><td>9 × 10<sup>-4</sup></td><td>0.0595</td><td>0.0889</td><td>0.6666</td><td></td><td>-0.0192</td><td>0.0299</td><td>0.0355</td><td>1.6690</td></tr></tbody></table>
<table><caption>Table 9: Fourier Whittle estimation of $\Omega$ for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. The number of frequencies is $m = \lfloor n^{\eta} \rfloor$ with $\eta = 0.65$ as chosen in [18].</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (0.2, 0)</th><th colspan="3">d = (0.2, 0.2)</th><th colspan="3">d = (0.2, 0.4)</th></tr><tr><th>bias</th><th>RMSE</th><th>W/F</th><th>bias</th><th>RMSE</th><th>W/F</th><th>bias</th><th>RMSE</th><th>W/F</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>0.0186</td><td>0.2027</td><td>0.4123</td><td>0.022</td><td>0.2036</td><td>0.3876</td><td>0.0113</td><td>0.201</td><td>0.4127</td></tr><tr><td>Ω<sub>1,2</sub></td><td>0.0124</td><td>0.0759</td><td>0.4381</td><td>0.0197</td><td>0.0686</td><td>0.3923</td><td>0.0149</td><td>0.0959</td><td>0.5545</td></tr><tr><td>Ω<sub>2,2</sub></td><td>0.0146</td><td>0.2138</td><td>0.3322</td><td>0.0263</td><td>0.2163</td><td>0.3624</td><td>0.0291</td><td>0.2162</td><td>0.4891</td></tr><tr><td>correlation</td><td>-0.0014</td><td>0.0357</td><td>0.8298</td><td>-2 × 10<sup>-4</sup></td><td>0.0355</td><td>0.4621</td><td>-0.0018</td><td>0.0359</td><td>0.8483</td></tr></tbody></table>

<table><caption>Table 10: Fourier Whittle estimation of $\Omega$ for a bivariate ARFIMA(0, d, 0) with $\rho = 0.8$, $N = 512$ with 1000 repetitions. The number of frequencies is $m = \lfloor n^{\eta} \rfloor$ with $\eta$ such that RMSE of $\hat{d}$ is minimized.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">d = (0.2, 0)</th><th colspan="3">d = (0.2, 0.2)</th><th colspan="3">d = (0.2, 0.4)</th></tr><tr><th>bias</th><th>RMSE</th><th>W/F</th><th>bias</th><th>RMSE</th><th>W/F</th><th>bias</th><th>RMSE</th><th>W/F</th></tr></thead><tbody><tr><td>Ω<sub>1,1</sub></td><td>0.0622</td><td>0.094</td><td>0.8896</td><td>0.0387</td><td>0.0832</td><td>0.9484</td><td>0.0389</td><td>0.0833</td><td>0.9962</td></tr><tr><td>Ω<sub>1,2</sub></td><td>0.0222</td><td>0.0638</td><td>1.1907</td><td>0.0304</td><td>0.0731</td><td>0.9379</td><td>0.047</td><td>0.0826</td><td>1.1615</td></tr><tr><td>Ω<sub>2,2</sub></td><td>-0.0031</td><td>0.0637</td><td>1.1149</td><td>0.0373</td><td>0.0839</td><td>0.9341</td><td>0.0812</td><td>0.1132</td><td>0.9338</td></tr><tr><td>correlation</td><td>-0.0013</td><td>0.0163</td><td>1.8222</td><td>-4 × 10<sup>-4</sup></td><td>0.0177</td><td>0.9237</td><td>-0.0012</td><td>0.0179</td><td>1.6969</td></tr></tbody></table>