Sophie Achard

Irène Gannaz

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

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

. The objective of this paper is to provide an overview of the R package multiwave with its practical application perspectives.

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(•) be the autocovariance function of a time series X. X is said to be long-memory if there exists α, 0 < α < 1, such that g(t) is asymptotically equivalent to |t| -α when t → +∞ (see [START_REF] Beran | Statistics for long-memory processes[END_REF] and references therein). This definition implies that the covariance function is not summable. Equivalently, the spectral density f (•), if it exists, is such that f (λ) is equivalent up to a constant to |λ| 1-α when λ → 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 longmemory for defining good estimators [START_REF] Beran | Statistics for long-memory processes[END_REF]. In the case of univariate time series, several very efficient approaches have been developed and validated. A web page entitled "Time Series Analysis" 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 [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF], arfima and FGN [START_REF] Veenstra | Persistence and Anti-persistence: Theory and Software[END_REF], longmemo [START_REF] Beran | Statistics for long-memory processes[END_REF] and 1 https://cran.r-project.org/web/views/TimeSeries.html 1 forecast [START_REF] Rob | Automatic time series forecasting: the forecast package for R[END_REF]. 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 multivariate 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 [START_REF] Peter M Robinson | Robust covariance matrix estimation: HAC estimates with long memory/antipersistence correction[END_REF]. One R package, waveslim, is dedicated to the wavelet correlation analysis for pairs of random variables [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF] but long-range dependence properties are not considered. [START_REF] Hurvich | The averaged periodogram estimator for a power law in coherency[END_REF] provide R code2 for bivariate longrange 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 [START_REF] Achard | Fractal connectivity of long-memory networks[END_REF] 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 [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF], 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 [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] 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 [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF]. 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.

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 [START_REF] Moulines | On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter[END_REF]. We consider N observations of a long-memory p-vector process X = {X (k), k ∈ Z, = 1, . . . , p}, namely X(1), . . . X(N ). X is said to be a multivariate M (d) process when for each = 1, . . . , p there exists D ∈ N such that the D -th order difference ∆ D X is covariance stationary. In addition, let us assume that for any , m = 1, . . . , p the generalized cross-spectral density of X and X m is

f ,m (λ) = 1 2π Ω ,m (1 -e -iλ ) -d (1 -e iλ ) -dm f S ,m (λ), λ ∈ [-π, π], (1)
with Ω = (Ω ,m ) ,m=1,...,p an Hermitian matrix. The functions f S ,m (•) correspond to the short-memory dynamics of the process. The parameters d satisfies -1/2 < d -D < 1/2. More generally, the wavelet-based procedure is available for cross-spectral density satisfying an approximation

f (λ) ∼ Λ(d)ΩΛ(d) * , when λ → 0, with Λ(d) = diag(|λ| -d e -i sign(λ)πd/2 ), (2)
where the exponent * denotes the conjugate transpose operator. Here and subsequently ∼ means that the ratio of the left-and right-hand sides converges to one. Note that, the process X 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 Ω 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 Ω ,m / Ω , Ω m,m ,m=1,...,p , which are also encountered in literature as power-law coherencies between two time series [START_REF] Hurvich | The averaged periodogram estimator for a power law in coherency[END_REF] or as fractal connectivities [START_REF] Achard | Fractal connectivity of long-memory networks[END_REF].

A parametric example: FIVARMA

Fractionally Integrated Vector Auto Regressive Moving Average (FIVARMA) processes are parametric models with a spectral density satisfying approximation [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF]. They correspond to Model A of [START_REF] Lobato | Consistency of the averaged cross-periodogram in long memory series[END_REF]. We refer to this paper for a detailed mathematical description. Let u be a p-dimensional white noise with

E[u(t) | F t-1 ] = 0 and E[u(t)u(t) T | F t-1 ] = Σ, where F t-1
is the σ-field generated by {u(s), s < t}, and Σ is a positive definite matrix. Let (A k ) k∈N be a sequence of R p×p -valued matrices with A 0 the identity matrix and

∞ k=0 A k 2 < ∞. The discrete Fourier transform of the sequence is denoted A(•), that is A(λ) = ∞ k=0 A k e ikλ .
We assume that all the roots of |A(L)| are outside the closed unit circle, where L denotes the lag operator. Let also (B k ) k∈N be a sequence in R p×p with B 0 the identity matrix and

∞ k=0 B k 2 < ∞. As defined for A, B(•) denotes the discrete Fourier transform of the sequence, B(λ) = ∞ k=0 B k e ikλ .
Let X be defined by

A(L) diag(1 1 -L) d X(t) = B(L)u(t). (3)
The spectral density satisfies

f ,m (λ) ∼ λ→0 + 1 2π Ω ,m e -iπ/2(d -dm) λ -(d +dm)
with

Ω = A(1) -1 B(1)ΣB(1) T A(1) T -1 . ( 4
)
Then X is called a F IV ARM A(d, q) process and satisfies approximation [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF].

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 [START_REF] Lobato | Consistency of the averaged cross-periodogram in long memory series[END_REF] and V ARF I models of [START_REF] Hurvich | The averaged periodogram estimator for a power law in coherency[END_REF] where

X is obtained with equation diag(1 1 -L) d A(L) X(t) = B(L)u(t)
, with similar notations than above. That is, X is obtained by fractional integration after autoregression, which is also called cointegration. The spectral density still satisfies the approximation (2) however the matrix Ω 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 [START_REF] Peter M Robinson | Multiple local Whittle estimation in stationary systems[END_REF][START_REF] Hurvich | The averaged periodogram estimator for a power law in coherency[END_REF][START_REF] Shimotsu | Exact local whittle estimation of fractionally cointegrated systems[END_REF]. Taking into account cointegration is a difficult problem that exceeds the scope of this paper. Future work is needed to handle this particular case.

Fourier-based estimation (MFW)

The discrete Fourier transform and the periodogram of X evaluated at frequency λ are defined as in [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF]'s procedure

W F (λ) = 1 √ 2πN N t=1 X(t)e itλ , I F (λ) = W F (λ)W F (λ) * .
Let λ j = 2πj/N , j = 1, . . . , m, be the Fourier frequencies used in estimation, m ∈ N.

Define Λ F j (d) = diag λ d j e i(π-λj )d/2 . The estimators ( d M F W , Ω M F W
) are minimizers of the criterion

L M F W (d, Ω) = 1 m m j=1 log det Λ F j (d)Ω(d)Λ F j (d) * + W F (λ j ) * Λ F j (d)Ω(d)Λ F j (d) * -1 W F (λ j ) .
The solution satisfies

d M F W = argmin d log det( Ω M F W (d)) -2 log(2)   1 m m j=1 λ j   , (5)
Ω M F W = Ω M F W ( d M F W ), (6)
with

Ω M F W (d) = 1 m m j=1 Re Λ F j (d) -1 I F (j)Λ F j (d) -1 .
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 [START_REF] Karim M Abadir | Nonstationarity-extended local Whittle estimation[END_REF], rather than the approximation [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF]. For an estimation based on the first-order approximation (2), one should replace Λ F j (d) by Λ

F (1) j (d) = diag λ d j e iπd/2
. [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] 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 d is decreased for the multivariate procedure with respect to a univariate one. It is worth mentioning that [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF] developed a similar estimation procedure, based on a rougher approximation of the cross-spectral density,

Λ F j (d) = diag λ d j .
Interestingly, the quality of estimation for the vector d is similar. Nevertheless, the estimation of the long-run covariance matrix Ω is biased since it does not take into account the phase-shift appearing in Λ F j (d). We refer to [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF] and to [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] for a more detailed study of these estimators and their consistency.

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 [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF]. Let (φ(•), ψ(•)) be respectively a father and a mother wavelets, satisfying regularity conditions, as stated in [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. At a given resolution j 0, for k ∈ Z, we define the dilated and translated functions

φ j,k (•) = 2 -j/2 φ(2 -j • -k) and ψ j,k (•) = 2 -j/2 ψ(2 -j • -k).
The wavelet coefficients of the process X are defined by

W j,k = R X(t)ψ j,k (t)dt j 0, k ∈ Z,
where X(t) = k∈Z X(k)φ(t -k). For given j 0 and k ∈ Z, W j,k is a p-dimensional vector W jk = W j,k (1) W j,k (2) . . . W j,k (p) where W j,k ( ) = R X (t)ψ j,k (t)dt. For any j 0, the process (W j,k ) k∈Z is covariance stationary [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. Let θ ,m (j) denote the wavelet covariance at scale j between processes X and X m , i.e. θ ,m (j) = Cov(W j,k ( ), W j,k (m)) for any position k. Let us introduce the wavelet scalogram

I W (j) = k∈Z W j,k W T j,k . (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(•), defined as

K(δ) = ∞ -∞ |λ| -δ | ψ(λ)| 2 dλ, δ ∈ (-α, M ). (8)
The wavelet Whittle procedure is described in [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. Let Λ j (d) = diag 2 jd , j 0 j j 1 . Let G(d) denote a p × p-matrix with ( , m)-th element equal to

G ,m (d) = f S (0)Ω ,m K(d + d m )cos(π(d -d m )/2). (9)
The estimators ( d

M W W , G M W W
) are defined by minimization of L M W M (d, G), with

L M W M (d, G) = 1 n j1 j=j0 n j log det Λ W j (d)G(d)Λ W j (d) + k W T j,k Λ W j (d)G(d)Λ W j (d) -1 W j,k .
The estimation is here based on a first-order approximation of the spectral density matrix around 0. The solutions of this problem satisfy

d M W W = argmin d log det( G M W W (d)) + 2 log(2)   1 n j1 j=j0 jn j   p =1 d , (10)
G M W W (d) = 1 n j1 j=j0 Λ W j (d) -1 I W (j)Λ W j (d) -1 , (11)
where I W (j) is the wavelet scalogram at scale j defined in [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF]. The long-run covariance matrix can then be estimated by

Ω M W W ,m = G M W W ,m ( d M W W )/(cos(π( d M W W -d M W W m )/2)K( d M W W + d M W W m )). (12)
This second step in the estimation of the long-run covariance matrix Ω 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 [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF]. 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. [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] established that the MWW estimators [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF] and [START_REF] Moulines | On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter[END_REF] are consistent under non-restrictive conditions. The rate of convergence for the estimation of the long-range parameters d is similar to the MFW estimator and is minimax. We refer to [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] for the detailed study of the asymptotic behaviour of MWW estimation.

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.

Simulations of FIVARMA

multiwave package proposes simulation functions for time series with a spectral density satisfying approximation [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF]. The main function is fivarma which computes a parametric FIVARMA process defined in Section 2.1.

The input parameters of a F IV ARM A(q, d, r) process are the covariance matrix Σ of the innovation process u, the vector AR (AutoRegressive) (A k ) k=0,...,q , A k ∈ R p×p , the vector MA (MovingAverage) (B k ) k=0,...,r , B k ∈ R p×p , and the vector of long-range parameters d ∈ R p . FIVAR model of [START_REF] Sela | Computationaly efficient methods for two multivariate fractionnaly integrated models[END_REF] is a subcase, corresponding to MA coefficients equal to zero. The parameters of fivarma are thus, in order, (N, d, cov_matrix, VAR, VMA) where cov_matrix= Σ and VAR and VMA denote respectively the sequences of matrices (A k ) k=0,...,q and (B k ) k=0,...,r .

The output of the function fivarma is a list with first the values X(1), X(2), . . . , X(N ) obtained by equation ( 3) with u(t) white noise with centered Gaussian distribution and covariance Σ. The second element of the list is the value of the matrix Ω defined in [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF].

fivarma is based on two other functions:

• fracdiff applies a vectorial fractional differencing procedure and corresponds to a FIVARMA(0,d,0).

• varma computes a realisation of a multivariate ARMA process and corresponds to the case d = 0.

Similar functions can be found in other packages (e.g. fracdiff and MST [START_REF] Ruey | Multivariate Time Series Analysis: with R and financial applications[END_REF]) but were re-implemented in multiwave package.

Example.

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> long_run_cov <-resp$long_run_cov R> long_run_cov [,1] [,2] [1,] 0.6049383 0.5854938 [2,] 0.5854938 0.9730806 R> par(mfrow=c(2,1),mai=c(0.5,1,0.5,0.5)) R> plot(x[,1],type='l',lty=1) R> plot(x[,2],type='l',lty=1) 0 50 100 150 200 250 -3 -2 -1 0 1 2 3 Index x[, 1] 0 50 100 150 200 250 -3 -2 -1 0 1 2 3 x[, 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 [START_REF] Maxim | Fractional Gaussian noise, functional MRI and Alzheimer's disease[END_REF]. Other data sets presenting similar features are coming from finance e.g. [START_REF] Songsiri | Topology selection in graphical models of autoregressive processes[END_REF], 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 [START_REF] Termenon | Reliability of graph analysis of resting state fmri using test-retest dataset from the human connectome project[END_REF]. The data set called brainHCP contains the time series of 1200 points in time and 89 regions of the brain. Figure 1

R> data(brainHCP) R> dim(brainHCP)

[1] 1200 89 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)

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 × N -matrix of observations, with general term x ,i = X (i), = 1, . . . , p and i = 1, . . . , 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 λ j = 2π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 [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF], the optimization is done in 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 d

M F W
and second the p × p-matrix

Ω M F W
. 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. [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF], [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] or [START_REF] Frank S Nielsen | Local Whittle estimation of multi-variate fractionally integrated processes[END_REF] propose a default value m = N 0.65 . This choice is discussed in the simulation study in Section 6.

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> # Estimation R> m <-N^(0.65) ## default value of Shimotsu R> res_mfw <-mfw(x,m) R> res_mfw $d [1] 0.1239810 0.3497609 $cov [,1] [,2] [1,] 0.5004685 0.5840763 [2,] 0.5840763 1.2547056

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 Ω 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.

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 [START_REF] Rob | Automatic time series forecasting: the forecast package for R[END_REF].

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,
R> xwav[seq(index[2]+1,index[3]),1]
Finally it is useful to compute the quantity K(δ) defined in equation [START_REF] Rob | Automatic time series forecasting: the forecast package for R[END_REF]. 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 i )) i=1,...,q * 2 J where q is the length of the filter and u i are equally spaced points on the interval [-π2 J-3 (q -1)/2; π2 J-3 (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(d) for a given value of d:

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

Estimation

mww is now described in detail. Let X be the p × N -matrix of observations, with general entries x ,i = X (i), = 1, . . . , p and i = 1, . . . , 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 [START_REF] Moulines | On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter[END_REF]. 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

M W W
and second the p × p-matrix

Ω M W W
. The quality of estimation is depending on the parameters j 0 and j 1 appearing in [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF] and [START_REF] Maxim | Fractional Gaussian noise, functional MRI and Alzheimer's disease[END_REF]. Default value for j 0 is set to 2 and for j 1 to the highest integer lower than log 2 (N ). The critical value to choose for estimation is j 0 , as it can be seen in the theoretical conditions for consistency [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] 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 ( 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 lon-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.

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.

Choice of the wavelet bases

Actually multiwave only proposes Daubechies' wavelets, which satisfy theoretical properties of [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. 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 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 [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF], 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. [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF] and references therein). For multivariate estimation [START_REF] Frank S Nielsen | Local Whittle estimation of multi-variate fractionally integrated processes[END_REF] proposes an extension of [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] based on the transform defined in [START_REF] Karim M Abadir | Nonstationarity-extended local Whittle estimation[END_REF]. However, this approach gives satisfactory results only for d < 1.5 and we decided not to implement it in the package for simplicity.

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 j 0 and j 1 defined in equations ( 10) and [START_REF] Maxim | Fractional Gaussian noise, functional MRI and Alzheimer's disease[END_REF] 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 [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF][START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF] 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 [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF], 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 [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF]. This is also true for the covariances as shown in Proposition 2 of [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. For all k ∈ Z, Cov(W j,k ( ), W j,k (m)) is equivalent to 2 j(d +dm) G ,m (d) when j goes to infinity, with G ,m (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 [START_REF] Abry | Wavelet analysis of long-range-dependent traffic. Information Theory[END_REF]. 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 , m = 1, . . . , p, for all k ∈ Z,

Cor(W j,k ( ), W j,k (m)) is equivalent to G ,m (d)/ G , (d)G m,m (
d) when j goes to infinity [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. 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.

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), q q q q q q q q q q -4 -2 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 √ bias 2 + std 2 . For clarity, all tables are displayed at the end of the paper.

A reference example

We first consider a simple example, with neither short-range dependence, nor nonstationarity. Time series were simulated using a bivariate F IV ARM A(0, d, 0) with a long-run correlation matrix Ω = 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 Σ = 1 0.8 0.8 1 , the AR coefficient is set equal to A = 0.8 0 0.2 0.6 and the MA coefficient is set equal to B = 0.4 0 0.2 0.7 .

Calculation was done on N = 512 observations for 100 replications.

q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0

(a)

n j wavelet correlation 254 124 59 26 10 2 q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0

(b)

n j wavelet correlation 254 124 59 26 10 2 q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0

(c)

n j wavelet correlation 254 124 59 26 10 2 q q q q q q q q q 0.0 0.2 0.4 0.6 0. 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 d. This illustrates that multivariate estimation improves the quality of estimation for 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.

Short-range dependence

Consider a F IV ARM A(1, d, 0) obtained with the model described in section 2.1 and given by the function fivarma. This case corresponds to a F IV AR model of [START_REF] Hurvich | The averaged periodogram estimator for a power law in coherency[END_REF]. The AR coefficient is taken equal to A = 0.8 0 0.2 0.6 and the correlation between the innovation processes equal to ρ = 0.8. More precisely let ε be a bivariate white noise process with covariance matrix Σ = 1 ρ ρ 1 and let u be the AR(1) process defined by u(t) + Au(t -1) = ε(t). The time series observation X(t) at time t satisfies (1 -L) d X(t) = u(t). The matrix Ω in equation ( 4) is equal to

Ω = (I + A) -1 Σ(I + A) -1 T 0.3086 0.2392 0.2392 0.3260 .
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 d and Ω are still satisfactory. The RMSE is very similar to the previous case with no short-range behaviour.

Nonstationarity

Nonstationary examples are simulated using values of d higher than 0.5. We consider order 1 or 2 of nonstationarity, that is d ∈ [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 Ω, 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.

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 and d m satisfying d -d m = 1 the long-run covariance matrix Ω is no longer identifiable with the wavelet-based procedure. Indeed, Proposition 2 in [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] states that in this particular case the covariance Cov(W j,k ( ), W j,k (m)) tends to 0 when the scale j tends to infinity. Figure 5 illustrates this approximation for a bivariate F IV ARM A(0, (0.2, 1.2), 0), a correlation matrix Ω = 1 ρ ρ 1 and ρ = 0.8. . The index of the horizontal axis displays the number of coefficients available. Calculation was done on N = 512 observations for 1000 replications.

q q q q q -0.5 0.0 0.5 1.0

FIVARMA(0,(0.2,1.2),0) n j
wavelet correlation 253 123 58 26 10 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0.5 0.0 0.5 When d -d m = 1, the estimator ( 12) is no longer defined. In practice, the quantity d -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 , d m ) will lead to an important bias in the estimation of Ω ,m . As it can be seen in Figure 6, the resulting bias increases in the neighbourhood of the non-identifiable lines d -d m = ±1.

When this situation occurs, say when the difference between d -d m is between 0.75 and 1.25, the estimation of d is not affected. But the user must be careful for the estimation of Ω. One solution is to differentiate or integrate one of the two processes. For example, Table 1 illustrates the non-identifiability of Ω in a bivariate F IV ARM A(0, 0.2 1.2 , 0). When differentiating the second component (with d 2 = 1.2) the estimator has again good performances.

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 F IV ARM A(0, d, 0) with a long-run correlation matrix Ω = 1 ρ ρ 1 and ρ = 0.8. The bivariate vector 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 [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] is m = N 0.65 . We also make comparisons with an optimal value computed by minimizing the RMSE. 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 [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF], 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 [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF] 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 1: Multivariate Wavelet Whittle estimation of Ω for a bivariate F IV ARM A(0, (0.2, 1.2), 0) with ρ = 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.

Without differentiation

With differentiation [START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF]. 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.

Ω

Estimation of the long-run covariance

Finally, Table 9 and Table 10 display results for the estimation of Ω 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 Ω 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.

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. q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0. q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0

(c)

n j wavelet correlation 254 124 59 26 10 2 q q q q q q 0.0 0.2 0.4 0.6 0.

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 [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] and [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF]. 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.

  displays 6 arbitrary signals from one subject in this data set.

Figure 1 : 4

 14 Figure 1: Plot of 6 arbitrary signals from a subject of fMRI data set. channel 4

Figure 2 :

 2 Figure 2: Input and Output of multiwave package in a two-dimensional case. INPUT OUTPUT

Figure 3 : 7 .

 37 Figure 3: Boxplots of the log2-variance of the wavelet coefficients at different scales for simulated bivariate F IV ARM A(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 + d m ) + log 2 (G ,m (d)), with = m = 1 for subfigure (a), = m = 2 for subfigure (b) and = 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 Σ = 1 0.8 0.8 1 , the AR coefficient is set equal to A = 0.8 0 0.2 0.6 and the MA coefficient is set equal to B = 0.4 0 0.2 0.7 . Calculation was done on N = 512 observations for 100 replications.

Figure 4 :

 4 Figure 4: Boxplots of the correlation of the wavelet coefficients at different scales for four different simulated data: (a) bivariate F IV ARM A(0, (0.4, 0.4), 0) (b) bivariate F IV ARM A(1, (0.4, 0.4), 0); (c) bivariate F IV ARM A(1, (0.4, 0.4), 1); and finally, (d) an example with non stationary time series, F IV ARM A(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

8 .

 8 The bivariate vector d is chosen in [0, 0.5) 2 , such that the time series are stationary.

Figure 5 :

 5 Figure 5: Boxplots of the covariance of the wavelet coefficients at different scales for a bivariate ARF IM A(0, (0.2, 1.2), 0) with Ω = 1 0.8 0.8 1 . The index of the horizontal axis displays the number

Figure 6 : 2 6. 1

 621 Figure 6: RMSE in the estimation of the cross-covariance term Ω 12 with respect to (d 1 , d 2 ). Estimation was done using multivariate Wavelet Whittle estimator in a bivariate F IV ARM A(0, (d 1 , d 2 ), 0) with Ω = 1 0.8 0.8 1 . 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 = ±1. Calculation was done on N = 512 observations for 1000 replications. (a) (b)

Figure 7 :

 7 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).

Figures 8

 8 Figures 8 and 9 display an example of long-memory parameter and long-run correlation estimated for one subject.

Figure 8 :Figure 9 :

 89 Figure 8: Histogram of d from a subject of fMRI data set.

  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:

	[1] 6
	R> index
	[1] 127 189 219 233 239 241
	R> length(x)
	[1] 256
	index gives the
	filter)
	R> xwav <-resw$dwt # returns the vector of the wavelet coefficients
	R> index <-resw$indmaxband
	R> Jmax <-resw$Jmax
	R> Jmax

  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

		bias	std	RMSE	bias	std	RMSE
	Ω 1,1	0.0935	0.0762	0.1206	0.0349	0.0678	0.0762
	Ω 1,2	4.1863	7.2103	8.3375	0.0261	0.06	0.0654
	Ω 2,2	0.2215	0.0819	0.2362	0.0292	0.07	0.0758
	correlation	3.5255	6.2935	7.2137	0.0003	0.0155	0.0155
	procedure than with Fourier-based procedure.

Table 2 :

 2 Multivariate wavelet Whittle estimation of d for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. For the estimation, j 0 = 1.

	d	bias	std	RMSE	ratio M/U
	0.2 -0.017	0.039	0.0425	0.7785
	0.0	0.013	0.0391	0.0412	0.9014
	0.2 -0.0313	0.0376	0.049	0.8960
	0.2 -0.0316	0.0378	0.0493	0.8805
	0.2 -0.017	0.0383	0.0419	0.7673
	0.4 -0.0442	0.0395	0.0592	0.7902

Table 3 :

 3 Wavelet Whittle estimation of Ω for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. For the estimation, j 0 = 1.

			d = (0.2, 0)	d = (0.2, 0.2)	d = (0.2, 0.4)
	Ω	bias	std RMSE	bias	std RMSE	bias	std RMSE
	Ω 1,1	0.0417 0.0724 0.0836	0.0343 0.0711 0.0789	0.0417 0.0717 0.083
	Ω 1,2	0.0382 0.0657 0.0759	0.0279 0.0626 0.0686	0.0673 0.0684 0.0959
	Ω 2,2	0.0048 0.0709 0.071	0.0323 0.0714 0.0784	0.0748 0.0748 0.1057
	correlation 0.0191 0.0227 0.0296	0.001	0.0164 0.0164	0.0194 0.0234 0.0304

Table 4 :

 4 Multivariate wavelet Whittle estimation of d for a bivariate ARF IM A(1, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. For the estimation, j 0 = 3.

	d	bias	std	RMSE	ratio M/U
	0.2	-0.0473	0.1213	0.1302	0.8472
	0.0	-0.0371	0.1266	0.132	0.8511
	0.2	-0.0623	0.1209	0.136	0.8848
	0.2	-0.0526	0.1258	0.1364	0.8714
	0.2	-0.066	0.1244	0.1408	0.9161
	0.4	-0.0584	0.1293	0.1418	0.8935

Table 5 :

 5 Wavelet Whittle estimation of Ω for a bivariate ARF IM A(1, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. For the estimation, j 0 = 3.

			d = (0.2, 0)		d = (0.2, 0.2)	d = (0.2, 0.4)
	Ω	bias	std RMSE	bias	std RMSE	bias	std RMSE
	Ω 1,1	-0.0067 0.0828 0.0831	0.0024 0.0859 0.0860	0.005	0.0889 0.089
	Ω 1,2	0.0551 0.0828 0.0995	0.0495 0.0794 0.0936	0.0507 0.0875 0.1012
	Ω 2,2	0.1363 0.131	0.1891	0.1412 0.1384 0.1977	0.1391 0.1426 0.1992
	correlation -0.0088 0.0805 0.081	-0.0386 0.0527 0.0653	-0.0358 0.0985 0.1047

Table 6 :

 6 Multivariate wavelet Whittle estimation of d for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. Nonstationary cases. For the estimation, j 0 = 2.

	d	bias	std	RMSE	ratio M/U	d	bias	std	RMSE	ratio M/U
	1.2 -0.0338	0.0762	0.0834	0.8510	2.2 -0.0421	0.0884	0.0979	0.8718
	1	-0.0276	0.0725	0.0776	0.8316	2	-0.0403	0.0862	0.0951	0.8516
	1.2 -0.043	0.0732	0.0849	0.8672	2.2 -0.0503	0.086	0.0996	0.8874
	1.2 -0.0411	0.0743	0.0849	0.8591	2.2 -0.049	0.0823	0.0958	0.8566
	1.2 -0.0338	0.0741	0.0814	0.8310	2.2 -0.0436	0.0868	0.0971	0.8651
	1.4 -0.0356	0.0797	0.0873	0.8344	2.4 -0.0429	0.0831	0.0935	0.8400
					29

Table 7 :

 7 Wavelet Whittle estimation of Ω for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. Nonstationary cases. For the estimation, j 0 = 2.

			d = (1.2, 1)	d = (1.2, 1.2)	d = (1.2, 1.4)
	Ω	bias	std RMSE	bias	std RMSE	bias	std RMSE
	Ω 1,1	-0.0049 0.1362 0.1363	0.0059 0.1369 0.137	-0.0047 0.1361 0.1361
	Ω 1,2	0.018	0.1168 0.1182	0.0093 0.1155 0.1158	0.0159 0.1266 0.1276
	Ω 2,2	-0.0027 0.1277 0.1277	0.0042 0.1386 0.1386	-0.0113 0.1487 0.1491
	correlation 0.0214 0.0475 0.0521	0.0051 0.0286 0.0291	0.0227 0.0507 0.0555
			d = (2.2, 2)	d = (2.2, 2.2)	d = (2.2, 2.4)
	Ω	bias	std RMSE	bias	std RMSE	bias	std RMSE
	Ω 1,1	-0.0383 0.1795 0.1835	-0.0253 0.1789 0.1807	-0.0361 0.1776 0.1812
	Ω 1,2	-0.0043 0.1565 0.1565	-0.0129 0.1493 0.1498	-0.0097 0.1602 0.1605
	Ω 2,2	-0.0318 0.1776 0.1804	-0.0276 0.1809 0.183	-0.0481 0.1813 0.1876
	correlation 0.0251 0.0604 0.0654	0.0087 0.0374 0.0384	0.0249 0.0626 0.0674

Table 8 :

 8 Multivariate Fourier Whittle estimation of d for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions depending on the number of frequencies m. x denotes the closest integer smaller than x.

			m = N 0.65 = 57				m = N η
	d	bias	std	RMSE	W/F	η	bias	std	RMSE	W/F
	0.2	-0.002	0.0576	0.047	0.9049 0.90 -0.0197	0.0271	0.0335	1.2689
	0	-9 × 10 -4	0.0593	0.0238	1.7295		-0.0033	0.0264	0.0267	1.5475
	0.2	-0.0033	0.0574	0.0531	0.9214 0.85 -0.013	0.0306	0.0332	1.4750
	0.2	-0.0031	0.0591	0.0522	0.9434		-0.0123	0.0293	0.0318	1.5503
	0.2	8 × 10 -4	0.0576	0.0579	0.7239 0.85 -0.0136	0.0308	0.0337	1.2455
	0.4	9 × 10 -4	0.0595	0.0889	0.6666		-0.0192	0.0299	0.0355	1.6690

Table 9 :

 9 Fourier Whittle estimation of Ω for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. The number of frequencies is m = n η with η = 0.65 as chosen in[START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF].

		d = (0.2, 0)		d = (0.2, 0.2)		d = (0.2, 0.4)
	Ω	bias RMSE	W/F	bias	RMSE	W/F	bias RMSE	W/F
	Ω 1,1	0.0186 0.2027 0.4123	0.022	0.2036 0.3876	0.0113 0.201	0.4127
	Ω 1,2	0.0124 0.0759 0.4381	0.0197	0.0686 0.3923	0.0149 0.0959 0.5545
	Ω 2,2	0.0146 0.2138 0.3322	0.0263	0.2163 0.3624	0.0291 0.2162 0.4891
	correlation -0.0014 0.0357 0.8298	-2 × 10 -4	0.0355 0.4621 -0.0018 0.0359 0.8483

Table 10 :

 10 Fourier Whittle estimation of Ω for a bivariate ARF IM A(0, d, 0) with ρ = 0.8, N = 512 with 1000 repetitions. The number of frequencies is m = n η with η such that RMSE of d is minimized.

		d = (0.2, 0)		d = (0.2, 0.2)		d = (0.2, 0.4)
		η = 0.9		η = 0.85			η = 0.85
	Ω	bias RMSE	W/F	bias RMSE	W/F	bias RMSE	W/F
	Ω 1,1	0.0622 0.094	0.8896	0.0387	0.0832 0.9484	0.0389 0.0833 0.9962
	Ω 1,2	0.0222 0.0638 1.1907	0.0304	0.0731 0.9379	0.047	0.0826 1.1615
	Ω 2,2	-0.0031 0.0637 1.1149	0.0373	0.0839 0.9341	0.0812 0.1132 0.9338
	correlation -0.0013 0.0163 1.8222	-4 × 10 -4	0.0177 0.9237 -0.0012 0.0179 1.6969

freely available http://pages.stern.nyu.edu/~rsela/VARFI/code.html

code available in matlab http://shimotsu.web.fc2.com/Site/Matlab_Codes.html

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Ãťne-Alpes and a grant from AGIR-PEPS, UniversitÃľ Grenoble AlpesâĂŞCNRS.


Assoul Abdelkarim
email: assak_maths@yahoo.fr

The application of the hypercomplex number in a three dimensional space in quantum computing

Keywords: hypercomplex number, superposition, qubit

In a publication about the hypercomplex number in a three-dimensional space, we defined a hypercomplex number in a three-dimensional space as follows: s = x + yi + zj with x, y, z real numbers and i, j purely imaginary numbers such as: i 2 = j 2 = -1 and ij = ji = 0.

The goal of the article is the use of the hypercomplex number in a three-dimensional space to facilitate the computation of the amplitude and the probability of the superposition of all qubit states that enters into the quantum computing improvement as well as quantum physics.

Introduction

The conventional computer calculates by manipulating the bits, each bit carries either a 1 or a 0 whereas the quantum computer calculates by handling the qubits [START_REF]Calculateur quantique[END_REF]2], each qubit carries either a 1 or a 0 or a superposition [START_REF]Principe de superposition quantique[END_REF] of a 1 and a 0, this way the quantum computer will be more powerful. Many articles were written around the qubit and the quantum computer [START_REF]Calculateur quantique[END_REF] where the complex numbers were used in the calculation of the amplitude noted z = x + iy and the probability noted |z| 2 = a 2 + b 2 .

Our study exactly is the replacement of complex numbers by hypercomplex numbers in a three dimensional Hilbertian space [START_REF] Abdelkarim | [END_REF] for calculating the amplitude [START_REF]Amplitude de probabilité[END_REF] and probability [START_REF]Amplitude de probabilité[END_REF][START_REF] Hanckowiak | [END_REF] for finding new results.

The states of superposition

The states of superposition of a qubit

In the general case, a qubit has infinity of state superposition (*) denoted by Dirac [START_REF]Notation bra-ket[END_REF] as follows in fact: = α + β whither and are the basic states and α, β two hypercomplex numbers in a three-dimensional Hilbert space such as: |α| 2 +|β| 2 = 1

We take: |α| 2 = p and |β| 2 = 1-p, p is a real number such as: 0 < p < 1

If , are the image vectors respectively of the numbers α and β in an orthonormal coordinate system (O, I, J, K) then M α , M β appertain respectively to the spheres S 1 (0, ), S 2 (0, ).

(2)

The ideal case is to take p = , then M α , M β S (0, ).

That is to say M α , M β appertain to the sphere of equation: x 2 +y 2 +z 2 = If we apply (1), with a simple calculation, we find an example among an infinity of examples, so we search α and β two hypercomplex numbers in a hilbertian space of dimension three such that: = α + β and |α| 2 +|β| 2 = 1, as in the following table [START_REF]Calculateur quantique[END_REF] :

Etats Amplitude : s = x+yi+zj Probabilité :| s| 2 = x 2 +y 2 +z 2 α = Β = |α| 2 +|β| 2 = + = 1 (*)
The states of superposition is a purely mathematical consequence of quantum theory while in quantum physics, there are several interpretations namely: The Copenhagen interpretation, Everett's theory David Deutsch, the interpretation of De Broglie-Bohm and Schrödinger cat.

Superposition states of two qubits

Two qubits are brought together in a superposition of states:

= α + β + λ + μ whither α, β, λ, μ are hypercomplex numbers in a hilbertian space of dimension three such that: |α| 2 +|β| 2 +|λ| 2 +|μ| 2 = 1

We take in this case: |α| 2 = p, |β| 2 = -p, |λ| 2 = q and |μ| 2 = -q,

p and q two real numbers with: 0 < p < , 0 < q < Then: |α| 2 +|β| 2 +|λ| 2 +|μ| 2 = p + -p + q + -q = 1 If , , , are vector images respectively of the numbers α, β, λ, μ then:

M α є S 1 (0, ), M β є S 2 (0, ), M λ є S 3 (0, ), M μ є S 4 (0, ) (5)
The ideal case is to take p = q = then M α , M β , M λ , M μ S (0, )

That is to say M α , M β , M λ , M μ appertain to the sphere of equation: x 2 +y 2 +z 2 = If we apply (4), with a simple calculation, we find an example among an infinity of examples, We search α, β, λ, μ four hypercomplex numbers in a Hilbert space of dimension three such that: = α + β + λ + μ and |α| 2 +|β| 2 +|λ| 2 +|μ| 2 = 1, as in the following table [START_REF]Calculateur quantique[END_REF]: = q, = -q, 0 < q < = , = -, <

Etats

Then:

If , i = are vector images respectively of the numbers α i , i = then the images Mα i , i = appertain respectively to the following spheres:

The ideal case is to take: p = q = = = whither Mα i , i = appertain to the sphere of equation:

S (0, ) (9)
that is to say the sphere of equation: x 2 +y 2 +z 2 = If we apply [START_REF]Notation bra-ket[END_REF], with a simple calculation, we find an example among infinity of examples, we search α i , i = eight hypercomplex numbers in a three-dimensional Hilbert space such as:

= α 1 + α 2 + α 3 + α 4 + α 5 + α 6 + α 7 +α 8
and As in the following table [START_REF]Calculateur quantique[END_REF]:

Etats Amplitude : s = x+yi+zj Probabilité :| s| 2 = x 2 +y 2 +z 2 α ₁ = | ₁ | 2 = α ₂ = | ₂ | 2 = α ₃ = | ₃ | 2 = α ₄ = | ₄ | 2 = α ₅ = | ₅ | 2 = α ₆ = | ₆ | 2 = α ₇ = | ₇ | 2 = α ₈ = | ₈ | 2 =
3. Generalization: the states of superposition of n qubits, nєN* n qubits are in a superposition of states: = such as = 1

Whither α i , i = are hypercomplex numbers in a hilbertian space of dimension three and , i = are all states.

In this general case, we have: Finally we find M ₂ⁿ -1 є S ₂ⁿ -1 (0, ), and M ₂ⁿ є S ₂ⁿ (0, )

| ₁ | 2 = p ₁ and | ₂ | 2 = -p ₁ such as: 0 < p ₁ < (10)
If , i = are vector images respectively of the numbers α i , i = So the images Mα i , i = appertain in the ideal case to the sphere S (0, )

That is to say the sphere of equation: x 2 +y 2 +z 2 = Note. We can take simple examples of ideal cases using spherical coordinates:

z = [8] With = , = arc cos ( ), = and =2π -

Remarque

We know that the order of magnitude of the radius of the atom r a is 10 -10 m [9] and the states of superposition of n qubits can be represented by 2 n points of the sphere of center O and radius r s = so if: r a = r s that is to say = 10 -10 , we find 2 n = 10 20 then n = ≈ 66.43.

This allows us to say that the necessary number of qubit to have a really efficient quantum computer with atoms is : n = 67

Conclusion

Our study shows that the states of superposition of a single qubit can be represented by any two points of the sphere S (0, ), the superposition states of two qubits can be represented by any four points of the sphere S (0, ), the superposition states of three qubits can be represented by any eight points of the sphere S (0, ) and the states of superposition of n qubits can be represented by any 2 n points of the sphere S (O, ).

According to my calculation in this article if we take the number of qubits n 67, we will have an ideal quantum computer!

3 .

 3 Superposition states of three qubitsThree qubits are brought together in a superposition of states: whither α i , i = are hypercomplex numbers in a hilbertian space of dimension three andIn fact if we take:

1 If

 1 1 | 2 = p ₂ⁿ⁻¹ and | ₂ⁿ | 2 = -p ₂ⁿ⁻¹ such as: 0 < p ₂ⁿ⁻¹ < Then: | ₁ | 2 + | ₃ | 2 + | ₃ | 2 + ………. + | ₂ⁿ -1 | 2 + | ₂ⁿ | 2 = . = , i = are vector images respectively of the numbers α i , i =Then the images Mα i , i = appertain respectively to the following spheres: