Nejib Zemzemi
email: nejib.zemzemi@inria.fr

A Domain Decomposition Approch in the Electrocardiography Inverse Problem

Keywords: Cardiac electrophysiology, Inverse problem, Bidomain equations, Heart-torso coupling, Finite element method, domain decomposition

The mostly used mathematical formulation of the inverse problem in electrocardiography is based on a least method using a transfer matrix that maps the electrical potential on the heart to the body surface potential (BSP). This mathematical model is ill based and a lot of works have been concentrating on the regularization term without thinking of reformulating the problem itself. We propose in this study to solve the inverse problem based on a domain decomposition technique on a fictitious mirror-like boundary conditions. We conduct BSP simulations to produce synthetic data and use it to evaluate the accuracy of the inverse problem solution.

Introduction

The inverse problem in cardiac electrophysiology also known as electrocardiography imaging (ECGI) is a new and a powerful diagnosis technique. It allows the reconstruction of the electrical potential on the heart surface from electrical potentials measured on the body surface. This non-invasive technology and other similar techniques like the electroencephalography imaging interest more and more medical industries. The success of these technology would be considered as a breakthrough in the cardiac and brain diagnosis. However, in many cases the quality of reconstructed electrical potential is not sufficiently accurate. The difficulty comes from the fact that the inverse problem in cardiac electrophysiology is well known as a mathematically ill-posed problem. Different methods based on Thikhnov regularization [START_REF] Ghosh | Application of l1-norm regularization to epicardial potential solution of the inverse electrocardiography problem[END_REF] have been used in order to regularize the problem, but still the reconstructed electrical potential is not sufficiently satisfactory. In this study we present a domain decomposition approach to solve the inverse problem.

Methods

The domain decomposition method that will be presented in this paper is tested on synthetical data. This data is generated by solving the forward problem of ECGs. In the following paragraphs we will present first, the forward problem then the domain decomposition method for solving the inverse problem.

Forward problem

The bidomain equations were used to simulate the electrical activity of the heart and extracellular potentials in the whole body (see e.g. [START_REF] Pullan | Mathematically modelling the electrical activity of the heart. From cell to body surface and back again[END_REF], [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF], [START_REF] Tung | A bi-domain model for describing ischemic myocardial D-C potentials[END_REF]). These equations in the heart domain Ω H are given by:

               A m C m Vm + I ion (V m ,w) -div σ i ∇V m = div σ i ∇u e + I stim , in Ω H , -div (σ i + σ e )∇u e = div(σ i ∇V m ), in Ω H , ẇ + g(V m , w) = 0, in Ω H , σ i ∇V m • n = -σ i ∇u e • n, on Σ. (1)
The state variables V m and u e stand for the transmembrane and the extracellular potentials. Constants A m and C m represent the rate of membrane surface per unit of volume and the membrane capacitance, respectively. I stim and I ion are the stimulation and the transmembrane ionic currents. The hearttorso interface is denoted by Σ. The intra-and extracellular (anisotropic) conductivity tensors, σ i and σ e , are given by σ i,e = σ t i,e I + (σ l i,e -σ t i,e )a ⊗ a, where a is a unit vector parallel to the local fibre direction and σ l i,e and σ t i,e are, respectively, the longitudinal and transverse conductivities of the intraand extra-cellular media. The field of variables w is a vector containing different chemical concentrations and various gate variables. Its time derivative is given by the vector of functions g.

The precise definition of g and I ion depend on the electrophysiological transmembrane ionic model. In the present work we make use of one of the biophysically detailed human ventricular myocyte model Ten Tusscher and [START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF]. The ion channels and transporters have been modeled on the basis of the most recent experimental data from human ventricular myocytes.

Figure 1 provides a geometrical representation of the domains considered to compute extracellular potentials in the human body. In the torso domain Ω T , the electrical potential u T is described by the Laplace equation.

div(σ T ∇u T ) = 0, in Ω T , σ T ∇u T • n T = 0, on Γ ext .
(

) 2
where σ T stands for the torso conductivity tensor and n T is the outward unit normal to the torso external boundary Γ ext .

Ω T Ω H Γ ext Σ Fig. 1
Two-dimensional geometrical description: heart domain Ω H , torso domain Ω T (extramyocardial regions), heart-torso interface Σ and torso external boundary Γext.

The heart-torso interface Σ is supposed to be a perfect conductor. Then we have a continuity of current and potential between the extra-cellular myocardial region and the torso region.

u e = u T , on Σ, (σ e + σ i )∇u e • n T = σ T ∇u T • n T , on Γ ext .
(3)

Other works [START_REF] Clements | Acivation dynamics in anisotropic cardiac tissue via decoupling[END_REF], [START_REF] Potse | ECG simulations with realistic human membrane, heart, and torso models[END_REF], [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] consider that the electrical current does not flow from the heart to the torso by assuming that the heart is isolated from torso. This approximation is appealing in terms of computational cost because it uncouples the Laplace equation (2) in the torso from the bidomain equations in the heart (1), which allows to reduce the size of the linear system to solve. It is even more appealing when the interest is only on the ECG computation, in that case the ECG solution could be an "off line" matrix vector multiplication after solving the bidomain equation, details about computing the transfer matrix could be found in [START_REF] Zemzemi | Étude théorique et numérique de l'activité électrique du coeur: Applications aux électrocardiogrammes[END_REF]. Although this approach is very appealing in terms of computational cost, numerical evidence has shown that it can compromise the accuracy of the ECG signals (see e.g. [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF], [START_REF] Pullan | Mathematically modelling the electrical activity of the heart. From cell to body surface and back again[END_REF], [START_REF] Boulakia | Mathematical modeling of electrocardiograms: a numerical study[END_REF]). Thus, in order to accurately compute ECGs we consider the state-of-the-art heart-torso full coupled electrophysiological problem (1)-( 3) representing the cardiac electrical activity from the cell to the human body surface.

Inverse problem

The inverse problem in electrocadiography imaging (ECGI) is a technique that allows to construct the electrical potential on the heart surface Σ from data measured on the body surface Γ ext . We assume that the electrical potential is governed by the diffusion equation in the torso as shown in the previous paragraph. For a given potential data d measured on the body surface Γ ext , the goal is to find the extracellular heart potential u e on the heart boundary Σ such that the electrical potential in the torso domain u T satisfies as follows,

     div(σ T ∇u T ) = 0, in Ω T , u T = d and σ T ∇u T .n = 0, on Γ ext , u T = u e , on Σ.
     div(σ T ∇u = 0, in Ω T , u = d, on Γ ext , σ T ∇u.n = -σ T ∇v.n, on Σ.
(5)

The function v is defined in the fictitious domain that we denoted -Ω T as shown in Figure 2.

     div(σ T ∇v) = 0, in -Ω T , σ T ∇v.n = 0, on -Γ ext , v = u, on Σ. (6)
In order to solve the coupled problem ( 5)-( 6), we use the domain decomposition technique. For a given initial guess u 0 = 0, we compute the incomplete boundary condition following the algorithm 1.

Algorithm 1

for i = f irst time step to last time step do u 0 = 0, on Γext. load the known boundary data u T (t i ) /Γ ext load the known boundary flux σ T ∇u T (t i ) /Γ ext .n error = 2 x tolerence. while error > tolerence) do Set v k+1 = u k on Σ compute v k+1 in -Ω T Set σ T ∇u k+1 .n = -σ T ∇v k+1 .n on Σ compute u k+1 in Ω T compute error and relative error: norm( u k+1 -u k ) end while save u T (t i ) /Σ = u k+1 /Σ end for
The domain decomposition algorithm described in Algorithm 1 is accelerated using Aitken algorithm [START_REF] Irons | A version of the aitken accelerator for computer implementation[END_REF]. At the end of the while loop we have that v k+1 = u k+1 up to the defined tolerance and σ T ∇u k+1 .n = -σ T ∇v k+1 .n on Σ. If we define vk+1 as the symmetrical image of v k+1 through the interface Σ, we obtain σ T ∇u k+1 .n = σ T ∇v k+1 .n and u k+1 = vk+1 on Σ. Following [START_REF] Faker | On cauchy's problem: I. a variational steklov-poincaré theory[END_REF], [START_REF] Zemzemi | A steklov-poincaré approach to solve the inverse problem in electrocardiography[END_REF], this allows us to conclude that u k+1 = vk+1 in the whole torso domain Ω T . Consequently, u k+1 satisfies both Dirichlet and Neumann boundary conditions on the external boundary Γ ext and the Laplace equation in the torso domain, which solves the inverse problem.

In order to test the domain decomposition approach for solving the inverse problem in electrocardiography, we generate synthetical data. We conducted numerical simulations solving the forward problem. We use the finite element LifeV 1 library for the numerical implementation of the method. For the sake of simplicity, we perform simulations on the volumes between three concentric spheres. The volume between the small and the medium sphere represents the heart domain and the volume between the medium and the large spheres represents the torso. We show the torso potential in Figure 3.The epicardium Σ is given by the internal sphere in Figure 3 while the body surface Γ ext is given by the external sphere in the same figure. In Figure 3 (top) we show the forward problem solution at times 5, 10 and 25 ms (from left to right). We extract the electrical potential at the boundary Γ ext and we use it as an input for the inverse problem after adding 10% of white noise. We show the solution of the inverse problem in Figure 3 (bottom). The three panels represent the inverse problem solution at times 5, 10 and 25 ms (from left to right). We observe that the space distribution of the electrical wave is well captured mainly in the internal sphere. This allows to construct the activation times on the heart surface with a good accuracy.

Conclusion

In this paper we presented a domain decomposition approach to solve the inverse problem in electrocardiography. The problem was formulated using a mirror-like boundary conditions at the heart torso interface, where we have continuity of potential and opposite fluxes. The preliminary results presented in this work show the capability of this method to capture the spatial distribution of the electrical wave. In particular the wave front is well captured even with a relatively high level of noise. In future work we will test this method on clinical data including CT-scans of real geometry of the torso and measurements of the electrical potential on the body surface.

Fig. 2

 2 Fig. 2 Two-dimensional geometrical description of the physical torso domain Ω T and its image (-Ω T ) through the boundary Σ. The arrows on the interface Σ show the opposite fluxes σ T ∇u.n = -σ T ∇v.n.

Fig. 3

 3 Fig. 3 Top (respectively, bottom): Snapshots of numerical solution of the forward (re-spectively, inverse) problem at times 5, 10 and 25 ms (from left to right). The unit in the color bar is mV.

www.lifev.org


Gilles De Truchis
email: gilles.de-truchis@univ-amu.fr

Florent Dubois
email: florent.dubois@univ-amu.fr

Guillaume Chevillon

Mathieu Faure

Sébastien Laurent

Keywords: Unbalanced cointegration, Fractional cointegration, No-arbitrage condition, Local Whittle likelihood, Commodity markets JEL: C22, G10

A necessary condition for two time series to be nontrivially cointegrated is the equality of their respective integration orders. Nonetheless, in some cases, the apparent unbalance of integration orders of the observables can be misleading and the cointegration theory applies all the same. This situation refers to unbalanced cointegration in the sense that balanced long run relationship can be recovered by an appropriate filtering of one of the time series. In this paper, we suggest a local Whittle estimator of bivariate unbalanced fractional cointegration systems. Focusing on a degenerating band around the origin, it estimates jointly the unbalance parameter, the long run coefficient and the integration orders of the regressor and the cointegrating errors. Its consistency is demonstrated for the stationary regions of the parameter space and a finite sample analysis is conducted by means of Monte Carlo experiments. An application to the no-arbitrage condition between crude oil spot and futures prices is proposed to illustrate the empirical relevance of the developed estimator.

The no-arbitrage condition between spot and future prices implies an analogous condition on their underlying volatilities. Interestingly, the long memory behavior of the volatility series also involves a long-run relationship that allows to test for the no-arbitrage condition by means of cointegration techniques. Unfortunately, the persistent nature of the volatility can vary with the future maturity, thereby leading to unbalanced integration orders between spot and future volatility series. Nonetheless, if a balanced long-run relationship can be recovered by an appropriate filtering of one of the time series, the cointegration theory applies all the same and unbalanced cointegration operates between the raw series. In this paper, we introduce a new estimator of unbalanced fractional cointegration systems that allows to test for the no-arbitrage condition between the crude oil spot and futures volatilities.

Introduction

This paper deals with the estimation of a general class of models, terms unbalanced cointegration systems, that encompass the well known triangular cointegration system. The former systems originate from very recent developments (see [START_REF] Hualde | Unbalanced Cointegration[END_REF], while the latter relies on the seminal definition of the cointegration (see [START_REF] Granger | Some properties of time series data and their use in econometric model specification[END_REF] which establishes that two time series y t and x t share a common stochastic trend, if (i) y t and x t are both integrated of order δ 2 , hereafter I(δ 2 ); (ii) there exists a nonnull scalar β so that e t = y tβx t ∼ I(δ 1 ) and δ 2δ 1 > 0.2 This definition is quite general and does not constrain integration orders to be integers. Nonetheless, following the paper of [START_REF] Engle | Co-integration and error correction: representation, estimation, and testing[END_REF], the literature has primarily investigated the particular case where observables are unit root processes, i.e. I(1), and a linear combination between them has short memory, i.e. I(0).

In a pioneer work, [START_REF] Cheung | A fractional cointegration analysis of purchasing power parity[END_REF] extended the model of [START_REF] Engle | Co-integration and error correction: representation, estimation, and testing[END_REF], allowing the integration order of the cointegration errors (i.e. δ 1 ) to be a real number. Their methodology is simple and operates in two steps if y t and x t are I(1). When y t and x t are fractionally integrated (i.e. δ 2 ∈ R), some new complications arise. While a part of the literature has focused on testing the homogeneity of integration orders between y t and x t (see e.g [START_REF] Robinson | Determination of cointegrating rank in fractional systems[END_REF][START_REF] Nielsen | Determining the cointegrating rank in nonstationary fractional systems by the exact local Whittle approach[END_REF][START_REF] Hualde | A simple test for the equality of integration orders[END_REF], another part has investigated the estimation of such systems with unknown integration orders (see e.g [START_REF] Robinson | Cointegration in fractional systems with unknown integration orders[END_REF]. In the latter case, many difficulties appear in the uniform treatment of the objective functions on the entire parameter space. Thereby, three cases are generally distinguished in the literature. Considering that a fractionally integrated process, is stationary when the integration order is less than 1/2 and nonstationary otherwise, and interpreting δ 2δ 1 as the cointegration strength, (i) the strong cointegration occurs when δ 2δ 1 > 1/2; (ii) the weak cointegration occurs when δ 2δ 1 < 1/2; (iii) the stationary cointegration occurs when δ 1 < δ 2 < 1/2. As [START_REF] Velasco | Gaussian semi-parametric estimation of fractional cointegration[END_REF] suggested, it is more efficient to estimate simultaneously all the parameters of interest (see also [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF]. Among recent contributions in this direction we can mention [START_REF] Hualde | Root-n-consistent estimation of weak fractional cointegration[END_REF] and [START_REF] Shimotsu | Exact local Whittle estimation of fractionally cointegrated systems[END_REF] concerning the weak and strong cointegration cases and [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] and [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] for the stationary case. 3 In this paper we are particularly interested in the stationary cointegration case which is mainly attractive in empirical finance where time series have long range dependence but are likely to be stationary (e.g. volatility, volume, closing prices of commodities). Investigating stationary cointegration is also of interest because spurious regression can occurs, whether y t and x t are stationary or not, as long as their integration orders sum up to a value greater than 1/2 (see [START_REF] Tsay | The spurious regression of fractionally integrated processes[END_REF]. Nonetheless, it is not so easy to identify whether or not the observed time series are stationary leading practitioners to possibly make a wrong decision. Accordingly, we are also interested in non-stationary region of the parameter space, although the proposed estimator is not theoretically designed to handle this case. In the present paper, this issue is investigating by means of Monte Carlo study.

Meanwhile, [START_REF] Hualde | Unbalanced Cointegration[END_REF] investigated a promising alternative avenue of research termed unbalanced cointegration. Let y t and x t be two observable time series integrated of orders δ 2 and δ 2 + ξ respectively.

Following [START_REF] Hualde | Unbalanced Cointegration[END_REF], unbalanced cointegration is likely to occur between y t and x t and cointegration theory, in the usual sense, is likely to apply between y t and x t (ξ), if there exists a linear combination between them which has less memory δ 1 . Accordingly, this type of cointegration does not differ from the original cointegration theory but is useful from an empirical point of view. [START_REF] Hualde | Unbalanced Cointegration[END_REF] suggests to estimate ξ by the difference between integration orders of y t and x t and discusses the consistency of the OLS estimator of β. With respect to the joint estimation of β and ξ, the standard OLS are unfeasible and [START_REF] Hualde | Estimation of long-run parameters in unbalanced cointegration[END_REF] investigated the consistency of the non-linear least squares estimates of ξ and β when δ 2 > 1/2 and δ 1 ≥ 0.4 Conversely, we focus here on the stationary case assuming that δ 2 ∈ (0, 1/2 -|ξ|) and δ 2 > δ 1 ≥ 0. We propose a local Whittle estimator of bivariate unbalanced fractional cointegration systems. Focusing on a degenerating band around the origin of the spectral density matrix, it estimates jointly the unbalance parameter, the long run coefficient and the integration orders of the regressor and the cointegrating errors without specifying the short run dynamics avoiding thereby the misspecification issue. The consistency of the estimator is discussed and highlight that β is (m/n) δ 1 -δ 2 -consistent when x t is appropriately filtered in the long run equation, with m the bandwidth number and n the sample size.

The finite sample properties are also investigated by means of Monte Carlo simulations for a wide range of specifications.

In the empirical part of the paper, we aim to test whether the no-arbitrage condition holds on the crude oil market. A naive approach would be to investigate the presence of cointegration between spot and future prices. But as argued by [START_REF] Brenner | Arbitrage, cointegration, and testing the unbiasedness hypothesis in financial markets[END_REF], they should not be cointegrated because commodity markets are subject to the so-called convenience yield, that is probably not a short memory process. Consequently the long-run equation is contaminated by an additive persistent component. However, as demonstrated by [START_REF] Liu | No-arbitrage conditions for storable commodities and the modeling of futures term structures[END_REF], the no-arbitrage condition can exist on commodity market if the convenience yield is non-negative and the no-arbitrage issue remains unsolved. Recently, [START_REF] Rossi | A no-arbitrage fractional cointegration model for futures and spot daily ranges[END_REF], have strengthened that the no-arbitrage condition between spot and future asset prices implies an analogous condition on their underlying volatilities. This original approach is convenient because it can be adapted to commodity markets under mild conditions on the volatility of the convenience yield. Accordingly, we aim to test for the no-arbitrage condition between the crude oil spot and futures volatilities rather than spot and future prices. The unbalanced cointegration framework is particularly appropriated here because we find some evidence that the persistent nature of the volatility can vary with the future maturity (i.e. unbalanced integration orders between spot and future volatility series).

The rest of the paper is laid out as follows. In Section 2, we introduce a bivariate stationary model of unbalanced cointegration. In Section 3, we develop the local Whittle estimator of unbalanced fractional cointegration (LWE-UFC). In Section 4 we demonstrate the consistency of the proposed estimator. Finite sample properties are investigated in Section 5. The application to the no-arbitrage condition between the volatilities of the spot and futures prices is proposed in Section 6. The Section 7 concludes the paper.

Proofs are given in 8. Additional results and simulations are reported in 10 and 11.

A stationary model of unbalanced cointegration

In the following, we say that a stochastic process ζ t , has long memory α ∈ (0, 1/2), if its spectral

density f ζ (λ) satisfies f ζ (λ) ∼ gλ -2α
as the frequency λ tends to 0, where the notation ∼ means that the ratio of the left and right sides tends to 1 in the limit. Then, ζ t has short memory when α = 0 and intermediate memory when α ∈ (-1/2, 0). 5 Now, we consider an unbalanced bivariate form of the triangular system introduced in [START_REF] Phillips | Optimal inference in cointegrated systems[END_REF] and extended to the fractional framework by [START_REF] Nielsen | Optimal Residual-Based Tests for Fractional Cointegration and Exchange Rate Dynamics[END_REF]. Let y t and x t be two unbalanced observable variables with unknown real integration orders, δ 2 and δ 2 + ξ n respectively. [START_REF] Hualde | Unbalanced Cointegration[END_REF] states that y t and x t are weakly unbalanced when δ 2 and δ 2 + ξ n does not diverge at infinity (i.e. ξ n → 0 as n → ∞) and strongly unbalanced when

ξ n = |ξ| > 0 as n → ∞.
To simplify notation, we shall use ξ to denote ξ n . Then, the triangular unbalanced cointegration system is defined by

y t = βx t (ξ) + u 1t (-δ 1 ), x t = u 2t (-δ 2 -ξ), t = 1, 2, ..., n, (1)
where generically, (-α) denotes the fractional filter (1

-L) -α = ∑ ∞ k=0 a k (α)L k with a k (α) := Γ(k + α)(Γ(α)k!) -1
, L, the lag operator and Γ(.), the gamma function. Although standard cointegration theory does not apply to y t and x t , it does to y t and x t (ξ). Thereby, the System (1) is a cointegration system in the sense that both series have a dominant common component with memory δ 2 that can be suitably recovered by filtering x t ∼ I(δ 2 + ξ) to obtain x t (ξ) ∼ I(δ 2 ).

Assumption 1. y t , x t and y tβx t (ξ) are covariance stationary processes integrated of orders δ 2 , δ 2 + ξ and δ 1 respectively, and satisfying

0 ≤ δ 1 < δ 2 < δ 2 + |ξ| < 1/2 (2)
where |ξ| < k, with k an arbitrary real number small compared to δ 2 .

Under Assumption 1, the System in (1) provides a valid data-generating process for stationary cointegration model so that z t = (y tβx t (ξ), x t ) possesses a spectral density, f z (λ j ), where λ j denotes the Fourier frequencies, λ j = 2πj/n, with j = 1, . . . , m and m = o(n), the bandwidth parameter. Assumption 1 leaves out anti-persistent processes because they clearly have limited economic relevance. Now, assume that u t = (u 1t , u 2t ) has short memory with spectral density f u (λ j ) satisfying, f u (λ j ) ∼ G in the neighborhood of the origin, with G (the long run covariance matrix) a real, symmetric, finite and positive definite matrix. Notice that when cointegration arises, rank(G) < 2, so that G has reduced-rank whether or not ξ 0 (see [START_REF] Hualde | Unbalanced Cointegration[END_REF].

5 When α ∈ [1/2, 1], the spectral density of ζ t is no longer defined although Velasco (1999b) demonstrated it has a pseudospectral density and standard local Whittle-based estimators are biased. In such a case, Velasco (1999b) proposed to use a tapered periodogram to reduce the bias, but at the cost of a higher variance. Exploiting this result, Shimotsu (2012) combined a multivariate extension of the exact local Whittle estimator of [START_REF] Shimotsu | Exact local Whittle estimation of fractional integration with unknown mean and time trend[END_REF] and a tapered version of the [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] to propose a fractional cointegration estimator which is consistent over -1/2 < δ 1 < δ 2 < ∞ and includes δ 2 = 1 and δ 1 = 0 as a special case.

Remark 1. Assumptions on u t entails only mild conditions so that u t is possibly a vector ARMA process or any other short memory processes with a Wold representation, u t = C(L)ε t , where ε t are further defined as martingale difference innovations and C(L) is a square-summable causal matrix filter satisfying G = C(1)C(1) (2π) -1 .

In the present study, we define processes only on the vicinity of the origin, in view of allowing a semiparametric treatment of the short-run dynamics. Indeed, we support that such approach is particularly of interest in cointegration analysis because empirical interest is more likely to lie in long run dynamics.

Accordingly, by Assumption 1 and considering the Remark 1, z t has a spectral density, so that

E z t -E(z t ) z t+k -E(z t ) = π -π e ikλ f z (λ)dλ, (3)
with f z (λ) = Λ(λ) -1 f u (λ) (Λ(λ) * ) -1 and Λ(λ) = diag (1 -e iλ ) δ 1 , (1 -e iλ ) δ 2 +ξ . Since (1 -e iλ ) α = (|2 sin(λ/2)|) α ∼ λ α as λ → 0,
we can avoid a parametric treatment of f z (λ) in favor of the following local power law representation around zero frequency,

f z (λ) ∼ Λ(λ; ϑ 1 ) -1 G Λ(λ; ϑ 1 ) * -1 , Λ(λ; ϑ 1 ) = diag λ δ 1 , λ δ 2 +ξ , as λ → 0 (4)
where ϑ 1 = (δ 1 , δ 2 + ξ) and the superscript * denotes the conjugate transpose. As [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] we assume that G is diagonal so that u 1t and u 2t are incoherent in the vicinity of the origin and the Equation ( 4) is correctly specified. 6 Thereby, the phase parameter modeled as ϕ = (πλ)(δ 2δ 1 )/2 in [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] and [START_REF] Shimotsu | Exact local Whittle estimation of fractionally cointegrated systems[END_REF] is null in our framework (see also [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF]. 7 Notice also that, if δ 2 ≤ δ 1 , β cannot be identified. This issue also occurs in standard regression analysis of balanced cointegration.

A matrix representation of the System (1) gives

1 -β(1 -L) ξ 0 1 y t x t = (1 -L) -δ 1 0 0 (1 -L) -δ 2 -ξ u 1t u 2t . ( 5
)
The local Whittle analysis of restricted versions of Equation ( 5) is not new. For instance, when β = 0, y t and x t do not share any long-run component, and the first matrix on the left side of (5) reduces to the identity matrix. Thereby, one can recover the so-called stationary ARFIMA model whose the local Whittle estimate has been studied in univariate framework by Robinson (1995a) and later extended to a multivariate setting by [START_REF] Lobato | A semiparametric two-step estimator in a multivariate long memory model[END_REF] and [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF]. When β 0 and ξ = 0, cointegration can arise in the usual sense. Considering this latter case, [START_REF] Robinson | Semiparametric frequency domain analysis of fractional cointegration[END_REF] and [START_REF] Nielsen | Semiparametric Estimation in Time-Series Regression with Long-Range Dependence[END_REF] discuss the estimation of β by means of frequency domain least squares (see also [START_REF] Robinson | Semiparametric analysis of long-memory time series[END_REF][START_REF] Christensen | Asymptotic normality of narrow-band least squares in the stationary fractional cointegration model and volatility forecasting[END_REF][START_REF] Nielsen | Fully modified narrow-band least squares estimation of weak fractional cointegration[END_REF]. As mention previously, approximate local

Whittle estimation of δ 1 , δ 2 and β is considered in [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF], [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] and [START_REF] Shimotsu | Exact local Whittle estimation of fractionally cointegrated systems[END_REF].

When ξ 0 and δ 2 > 1/2, Hualde (2014) derives the asymptotic properties of the non-linear least squares estimator of β and ξ.

Local Whittle estimation

In the following, we introduce a local Whittle estimator of θ = (δ 1 , δ 2 , ξ, β) . Let, I z be the periodogram matrix of z t , defined as I z (λ j ; ϑ 2 ) = w z (λ j ; ϑ 2 )w z (λ j ; ϑ 2 ) * with w z (λ j ; ϑ 2 ) = (2πn) -1/2 ∑ n t=1 z t e itλ j , the Fourier transform of z t and ϑ 2 = (β, ξ) . Updated to bandwidth m = o(n), i.e. j = 1, . . . , m, in view of the local treatment we obtain,

I z (λ j ; ϑ 2 ) = w y (λ j ) -βλ ξ j w x (λ j ) w x (λ j ) w y (λ j ) -βλ ξ j w x (λ j ) w x (λ j ) * . ( 6
)
Thereby, the presence of λ ξ j corrects for the fact that long memory parameters of y t and x t are unbalanced. Then, the discrete local Whittle approximation to the likelihood is

Q m (θ, G) = m -1 m ∑ j=1 log det Λ(λ j ; ϑ 1 ) -1 G Λ(λ j ; ϑ 1 ) * -1 + tr G -1 Λ(λ j ; ϑ 1 )I z (λ j ; ϑ 2 )Λ(λ j ; ϑ 1 ) * , (7)
where G ∈ Θ G , the set of real positive definite 2 × 2 matrices. The objective function Q is minimized over

Θ G by Ĝ(θ) = Re m -1 m ∑ j=1 Λ(λ j ; ϑ 1 )I z (λ j ; ϑ 2 )Λ(λ j ; ϑ 1 ) * , (8)
leading to the following concentrated likelihood function

R m (θ) = log det Ĝ(θ) - 2(δ 1 + δ 2 + ξ) m m ∑ j=1 log λ j . (9)
Accordingly, the local Whittle estimator of θ is defined as θ

= arg min θ∈Θ R m (θ), for m ∈ [1, n/2] and Θ a compact subset of R 4 with Θ = Θ δ × Θ ξ × Θ β and δ = (δ 1 , δ 2 )
. Vectors θ and θ are respectively the vector of unknown and estimated values, (δ 1 , δ 2 , ξ, β) and ( δ1 , δ2 , ξ, β) . Observe that Equation (8)

yields Ĝ11 (θ) = Re m -1 ∑ m j=1 λ 2δ 1 j I 11 z (λ j ; ϑ 2 ) with I 11 z (λ j ; ϑ 2 ) = I yy (λ j ) -2βλ ξ j I xy (λ j ) + β 2 λ 2ξ j I xx (λ j
), which has some similarities with the weighted least squares of [START_REF] Nielsen | Semiparametric Estimation in Time-Series Regression with Long-Range Dependence[END_REF]. In 10, we show that the Proposition 1 of Nielsen (2005, p. 297) and the Theorem 1 of [START_REF] Robinson | Semiparametric frequency domain analysis of fractional cointegration[END_REF] remain valid when ξ 0 if x t is appropriately differenced. Accordingly, we anticipate that β is also λ δ 01 -δ 02 m -consistent when ξ 0.

Consistency

To prove the consistency of this local Whittle estimator, we introduce several assumptions, fairly similar to those of [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF] and [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF]. In the following, θ 0 and G 0 will denote the true parameter values of θ and G. Then, let f ab z (λ) and G 0 ab denote the (a, b)th element of f z (λ) and G 0 respectively. Define also ϑ 01 = (δ 01 , δ 02 + ξ 0 ) and δ 0a the ath element of ϑ 01 .

Assumption 2. As λ → 0 + , elements of the spectral density f z (λ) satisfies

f ab z (λ) = G 0 ab λ -δ 0a -δ 0b + o(λ -δ 0a -δ 0b ), a, b = {1, 2}, (10)
where matrix G is finite, real, symmetric and positive definite. Also assume G 12 = G 21 = 0.

Assumption 3. The sequence z t is a linear process defined as

z t -E(z t ) = A(L)ε t = ∞ ∑ j=0 A j ε t-j , ∞ ∑ j=0 ||A j || 2 < ∞, (11)
with ||.|| the Euclidean norm, so that A j is a causal square summable matrix filter. Moreover, ε t satisfies, almost surely, E(ε t |F t-1 ) = 0, E(ε t ε t |F t-1 ) = I 2 , with F t a σ-field generated by {ε s , s ≤ t} and there exist a random variable ε such that E(ε 2 < ∞) and for all η > 0 and some K > 0,

Pr(||ε t || 2 > η) ≤ KPr(ε 2 > η).
Assumption 2 implies a zero-coherence condition that applies only in the vicinity of the origin. As argued in [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF], it is a less restrictive assumption than the traditional orthogonality condition encountered in the least squares theory. Notably, it allows for errors to be correlated away from the origin and share, for instance, a common short-and/or medium-term dynamics. As mentioned previously, [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] and [START_REF] Shimotsu | Exact local Whittle estimation of fractionally cointegrated systems[END_REF] relax this hypothesis. The present estimator could be modified to model f z correctly when G 12 is non-null by specifying the phase parameter and presumably, its consistency could be demonstrated in this case but we do not pursue that possibility further. Assumption 3 imposes uniformly square integrable martingale-difference innovations with constant conditional variance in view of the application of the standard CLT for martingale-difference arrays. 8 The latter assumption on conditional variance could also probably be relaxed assuming boundedness of higher moments as in [START_REF] Robinson | Long and short memory conditional heteroskedasticity in estimating the memory parameter of levels[END_REF] but we do not investigate this issue further.

Assumption 4. In a neighborhood of the origin, A(λ) = ∑ ∞ j=0 A j e ijλ is differentiable and

∂ ∂λ A a (λ) = O(λ -1 ||A a (λ)||) as λ → 0 + (12)
where A a (λ) is the a-th row of A(λ).

Assumption 4 implies ∂A a (λ)/∂λ = O(λ -δ a -1 ) because by the Cauchy inequality

||A a (λ)|| ≤ (A a (λ)A * a (λ)) 1/2 = (2π f aa (λ)) 1/2 .
Thereby, under Assumption 3 and 4 we have f z (λ) = (2π) -1 A(λ)A(λ) * .

Assumption 5. As n → ∞, the bandwidth parameter satisfies

1 m + m n → 0. ( 13
)
where m = n k , k ∈ (0, 4/5].

The bandwidth requirement defined by the Assumption 5 ensures that m tends to ∞ as n → ∞ but at a slower rate to remain in a neighborhood of the origin. The bandwidth parameter m is theoretically bounded by n 4/5 but in practice a too small bandwidth increases the variance of the estimator while a too large m generally increases the bias. Accordingly, Assumptions 2-5 are analogous to Assumptions 1-4 of [START_REF] Shimotsu | Gaussian semiparametric estimation of multivariate fractionally integrated processes[END_REF], [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] and natural multivariate extensions of Assumptions A1-A4 of Robinson (1995a). Under Assumptions 2-5 we may now state the following theorem which establishes the convergence rate of θ.

Theorem 1. Let Assumption 1-5 hold. Define ν 0 = δ 02 -δ 01 . Then, for θ 0 ∈ Θ and δ 01 < δ 02 ≤ δ 02 + |ξ 0 |, as n → ∞, δ1 -δ 01 δ2 -δ 02 = O p (m -1/2 ) (14) ( ξ -ξ 0 ) = O p (m -1/2 ) (15) ( β -β 0 ) = O p (m -1/2 (n/m) -ν 0 ) (16)
The convergence rate of β confirms that the system is rebalanced when x t is appropriately filtered. Notice that when δ 1 → 0 and

δ 2 → 1/2, β is almost √ n-consistent because m 1/2 m δ 1 -δ 2 →-1/2 n δ 2 -δ 1 →1/2 ( β - β 0 ) p -→ 0.
More importantly, when the presence of ξ 0 is neglected and thus x t (ξ) is trivially replaced by x t in Equation ( 1), β -

β 0 = O p ((n/m) δ 1 -(δ 2 +ξ)
) and β is likely to be inconsistent if ξ < 0.

Monte Carlo experiment

Simulation design

This section discusses the finite sample performance of the proposed estimator by means of Monte Carlo simulations. As argued in [START_REF] Hurvich | Estimation of the memory parameter for nonstationary or noninvertible fractionally integrated processes[END_REF], it is not so easy for practitioners to identify whether or not the observed time series are stationary. We deal with this issue by considering a data generating process (DGP) that accommodates the mean-reverting non-stationary regions of the parameter space, although the developed estimator is not theoretically designed to handle this case. Thereby, we generate a fractionally cointegrated system according to the following model,

y t = βx t (ξ) + u # 1t (-δ 1 ), x t = u # 2t (-δ 2 -ξ), t = 1, 2, ..., n, (17)
where for a generic process

ζ t , ζ # t = ζ t l(t ≥ 1), with l(.) the indicator function. ζ # t (-α) denotes the fractional truncated filter ζ # t (-α) := ∑ t-1 k=0 Γ(k + α)(Γ(α)k!) -1 ζ t-k .
The System in (17) provides a valid DGP for both, stationary and non-stationary regions of the parameter space, given that y t and x t are type II processes, thereby contrasting with our system in (1) which is of type I. [START_REF] Marinucci | Alternative forms of fractional Brownian motion[END_REF] shown that type I and type II processes are asymptotically equivalent in the stationary regions of the parameter space if both are generated from the same short memory sequence (see also [START_REF] Robinson | The distance between rival nonstationary fractional processes[END_REF]. Accordingly, the type II representation is often retained in simulation study for its simplicity.

Nonetheless, [START_REF] Davidson | Type I and type II fractional Brownian motions: A reconsideration[END_REF] demonstrated that both representations can substantially differ in finite sample. We account for their recommendations by performing some additional simulations with a type I representation. [START_REF] Davidson | Type I and type II fractional Brownian motions: A reconsideration[END_REF] discuss the benefits and limitations of several techniques devoted to type I simulation. Given that we are interested in generating simple linear fractional processes, we use the circulant embedding method extended to multivariate fractional Gaussian noise by [START_REF] Helgason | Fast and exact synthesis of stationary multivariate Gaussian time series using circulant embedding[END_REF].

In our experiment we consider the following settings. The stationary cointegration case is explored for δ 2 = 0.4 and δ 1 = {0.0, 0.2, 0.3}. Similarly, the strong and weak cointegration cases, are investigated for δ 1 = {0.0, 0.2, 0.4} and δ 1 = {0.4, 0.6, 0.8} respectively, with δ 2 = {0.6, 0.8, 1.0}.9 Because in practice the weakly unbalanced cointegration case is generally indistinguishable from that of balanced cointegration, we do not consider this case in the simulation. Conversely, the strongly unbalanced cointegration case has greater applicability and is investigated for ξ = {-0.1, 0.1}. The long-run coefficient is fixed as β = 1. The vector u t is generated from a bivariate normal distribution N 2 (µ, Σ) distribution. For each simulation, we report the bias, the variance and the Root Mean Squared Error (RMSE), defined by 1 I ∑ I i=1 E ( θiθ) 2 := Var( θ) + Bias( θ|θ) 2 with I the number of replications set to 10000. Following [START_REF] Shimotsu | Exact local Whittle estimation of fractionally cointegrated systems[END_REF], we use a penalty term defined as Π(β, β) = min(0, β -β + C) 4 + max(0, β -β -C) 4 to govern the objective function when the space for dimensionality reduction is weak, with β = βLSE , so that it preserves the asymptotic results obtained in Equation (2). Alternatively, the Narrow-Band Least Squares (NBLS) estimate has also been used for initialization and does not significantly modified the results. The initial estimates for δx and δy are obtained from the local Whittle estimator of Robinson (1995a) applied to x t and y t respectively. Therefore, the initial estimate for ξ, namely ξ, is based on the difference between δy and δx (see Hualde 2006, p. 777). The initial estimate, βLSE , results from the regression of y t and x t ( ξ). All computations are performed using MATLAB 2013a.

Simulation results

Table 1 reports the results for 10000 replications of the stationary cointegration model when ξ = 0.1 and ρ = 0. Accordingly, this specification satisfies all the assumptions of the model defined in Equation (1). In absence of short-run dynamics, taking frequencies away from the origin essentially impact the variance (rather than the bias) compare with m = n 0.5 . Thereby, in both cases, the approximation in Equation ( 4) is close to f z (λ) = (2π) -1 A(λ)A(λ) * . Observe that the strength of the cointegration and the variance of β are negatively linked which is in accordance with the limit theory (the convergence rate of β depends on the strength of the cointegration). In all cases, the variances (the RMSE so on) decrease as the sample size increases. Moreover, variances are lower when m = n 0.8 . Table 2 reports the results for 10000 replications of the stationary cointegration model when ξ = 0.1 and ρ = 0.4. Thereby, the off-diagonal elements of G are non-null and a correlation between y tβx t (ξ) and the regressor is introduced at all frequencies. Consequently, Assumption 2 is now violated and the estimator faces an additional complication. Nonetheless, simulations do not explicitly reproduce these theoretical results. For instance, the long memory parameters are estimated fairly precisely, especially when m = n 0.8 .

Concerning ξ < 0, we only report the stationary cointegration case with ρ = 0 because the results are very similar as we can see in the lower part of the Table 3. The simulations based on the type I representation are reported in the upper part of the Table 3. We limit our investigation to this case because the computations are highly time consuming but we support that these results are sufficiently informative. We can observe that the differences with the type II representation are sometimes substantial.

For instance, whatever the sample, the bias slightly differs when δ 1 = 0. In some cases, the variance increases but the divergences disappear when δ 1 = 3.

Tables 7, 8, 9 and 10 report the results for strong and weak fractional cointegration cases with either ρ = 0 or ρ 0. In all cases, there is weak evidence of consistency given that RMSEs slowly decrease when the sample size increases. This result suggests that practitioners should devote a particular attention to the stationarity or non-stationarity of the data. In the latter case, they should use the estimator of Hualde (2014).

Empirical illustration

In this section we aim to test whether the no-arbitrage condition holds on the crude oil market. A naive approach would be to investigate the presence of cointegration between spot and future prices. But as argued by [START_REF] Brenner | Arbitrage, cointegration, and testing the unbiasedness hypothesis in financial markets[END_REF], they should not be cointegrated because commodity markets are subject to the so-called convenience yield, that is probably not a short memory process. Consequently the long-run equation is contaminated by an additive persistent component. However, as demonstrated by [START_REF] Liu | No-arbitrage conditions for storable commodities and the modeling of futures term structures[END_REF], the no-arbitrage condition can exist on commodity market if the convenience yield is non-negative and the no-arbitrage issue remains unsolved. Recently, Rossi and Santucci de Magistris ( 2013), have strengthened that the no-arbitrage condition between spot and future asset prices implies an analogous condition on their underlying volatilities. This original approach is convenient because it can be adapted to commodity markets under mild conditions on the volatility of the convenience yield. Accordingly, we aim to test for the no-arbitrage condition between the crude oil spot and futures volatilities rather than spot and future prices.

The no-arbitrage condition states that on financial markets, risk-free arbitrage opportunities cannot arise. Paradoxically, [START_REF] Grossman | On the impossibility of informationally efficient markets[END_REF] emphasized that arbitrage opportunities, even infrequent, are necessary to make the market sufficiently incentive for market participants. Investigating the no-arbitrage condition is also interesting because it is closely related to the efficient market hypothesis issue. Denoting F t+k|t and S t the futures and the spot prices respectively, the no-arbitrage condition implies F t+k|t = S t • e k•r t+k|t with r t+k|t the return of a risk-free asset that expires at time t + k and e k•r t+k|t the cost of carry premium. [START_REF] Rossi | A no-arbitrage fractional cointegration model for futures and spot daily ranges[END_REF] show that using the rescaled daily range, σ t,X = (log 2) -1/2 (max τ log X τmin τ log X τ ), for X any price, the no-arbitrage condition is directly related with the second moments of the spot and futures prices by the following equation,

σ t,F = σ t,S + b t + u t , b t = (log 2) -1/2 (r τ max -r τ min ), t -1 < τ ≤ t (18)
where r τ max and r τ min stand for the risk-free rate in correspondence, respectively, of the highest and lowest log-price in a given day. The additional term u t is justified by the presence of market frictions whereas b t could be omitted because the risk-free asset intraday variations are negligible.

Unfortunately, Equation ( 18) is not correct with respect to our framework because we focus on commodity markets whereas Rossi and Santucci de Magistris (2013) focus on stock index futures and thus neglect the so-called convenience yield. The convenience yield is defined as "the flow of benefit of immediate ownership of a physical commodity" (see [START_REF] Liu | The stochastic behavior of commodity prices with heteroskedasticity in the convenience yield[END_REF]. Reformulating the no-arbitrage condition we obtain

F t+k|t = S t • e k•(r t+k|t -c t+k|t ) (19)
where c t+k|t represents the net convenience yield that is the difference between the convenience yield and the cost of storage applicable from time t to time t + k. Following the same reasoning as Rossi and Santucci de Magistris (2013), our modified version of Equation ( 18) should hold if the intraday volatility of the convenience yield is negligible. However, the convenience yield is unobservable and this assumption requires some explanations. Implications of this hypothesis are twofold: (i) it ensures that Equation ( 18) is not contaminated by the presence of an additional volatility component; (ii) it reduces the probability for the convenience yield to be negative. Indeed, in the latter case, the futures prices would be too high relative to the spot price, thereby offering an arbitrage opportunity. In financial literature, the net convenience yield is often modeled as an Ornstrein-Uhlenbeck process. In such a representation the convenience yield volatility is generally high and negative values may frequently occur. Nonetheless, in a recent paper, Liu and Tang (2010) consider Cox-Ingersoll-Ross representation of the convenience yield and demonstrate that its volatility is no necessarily high. Accordingly, in line with their results, we will assume that the intraday volatility of the convenience yield is sufficiently low to be neglected in Equation ( 18) and to avoid frequent arbitrage-free opportunities.

In the following, we conjecture that when substituting the daily range by the daily squared returns, the relation in Equation ( 18) remains valid.10 Accordingly, we can test whether the no-arbitrage condition holds, by estimating the Equation ( 18). Nonetheless, the persistent nature of the volatility reveals that long run components drive the underlying processes. In such a case, testing for the presence of cointegration is useful to guard against the risk of spurious regression but also because the dynamics of σ t,F and σ t,S are likely to slightly diverge in short run. Furthermore, one would expect that the more the maturity of the futures is far away, the more the volatilities are likely to drift far apart. Some evidences of this mechanism are provided by [START_REF] Caporale | Time-varying spot and futures oil price dynamics[END_REF] on the spot and futures prices.

We focus on WTI crude oil spot and futures prices traded in NYMEX. Our data set runs from January 2, 1996 to December 16, 2013 for a total of n = 4499 observations.11 We consider four different maturity contracts. The contract Futures 1 specifies the earliest delivery date. It expires on the third business day prior to the 25th calendar day of the month preceding the delivery month. If the 25th calendar day of the month is a non-business day, trading ceases on the third business day prior to the business day preceding the 25th calendar day. The contracts Futures 2-4 represent the successive delivery months following the contract Futures 1. Graphics 2 and 3 represents the daily log squared returns of the spot prices and the four futures.

To test for the presence of stationary (unbalanced) cointegration we apply a rigorous methodology.

First of all, we estimate the long memory parameters of each volatility series. Because their is no particular reasons for the volatility to be stationary, we use the two-step exact local Whittle (2S-ELW) estimator of [START_REF] Shimotsu | Exact local Whittle estimation of fractional integration with unknown mean and time trend[END_REF] which is robust to the presence of unknown mean and polynomial time trend. The results are reported in Table 4 for several bandwidths. Clearly, they are fairly homogeneous for a given Importantly, we also apply the procedure of [START_REF] Qu | A test against spurious long memory[END_REF] to test the null hypothesis that the volatility estimate is a stationary long memory process against the alternative of a process contaminated by level shifts or a smoothly varying trend. The so-called W-test statistic of [START_REF] Qu | A test against spurious long memory[END_REF] depends on a trimming parameter ε that we set to either 0.02 or 0.05 and for which there are two specific asymptotic critical value at a threshold of 10%: 1.118 and 1.022. We compute the test using the Gaussian semi-parametric estimator of Robinson (1995a). The results are reported in Table 5 and highlight that in all cases we do not reject the null of a stationary long memory process except for the volatility of the Futures 1 for which it is difficult to conclude. We turn now to the core of our analysis that is the (non) equality of the integration orders. To investigate whether or not the integration orders of the pairwise volatilities are homogeneous, we apply the procedure of [START_REF] Nielsen | Determining the cointegrating rank in nonstationary fractional systems by the exact local Whittle approach[END_REF]. Because the cointegration is not observed, the authors propose to test for H 0 : δ F i = δ S with i = 1, 2, 3, 4 that is informative in both cases. However, computing the 6 and lead to ambiguous conclusions. Indeed, in several cases we accept the alternative hypothesis at a threshold of 10% (i.e. 2.71) although the rank estimates are positives in most cases.

In such inconclusive situation we support that unbalanced stationary cointegration is likely to occurs.

The estimates β, δ2 , δ1 and ξ are also reported in Table 6. First of all, we observe that all ξ are negative and increases (in absolute value) with the maturity of the futures. According to the 2S-ELW estimates of δ S and δ F i , this is not surprising. Also notice that the cointegration strengths, v, are clearly non-null and decrease as the maturity of the futures increases. This result is very interesting with respect to hedging strategies and is consistent with the findings of [START_REF] Caporale | Time-varying spot and futures oil price dynamics[END_REF]. Moreover, it reveals that the no-arbitrage hypothesis seems valid for futures contracts maturing in one month whereas the results are more ambiguous for contracts with longer maturities. Indeed, the coefficient β also dramatically decreases as the maturity of the futures increases leading to lower long run hedging ratio for Futures 2-4. A possible explanation comes from the fact that the probability of occurrence of negative convenience yield -and hence the probability arbitrage opportunities -increases with the maturity.

Conclusion

In this paper we investigate the estimation of bivariate unbalanced stationary fractional cointegration models. The estimator we propose relies on the local behavior of the spectral density of the system in the vicinity of the origin. It allows for estimating jointly all parameters of interest of the model and notably the unbalanced parameter. We have demontrated the consistency of the estimator as well as its good finite sample properties estimator by means of Monte Carlo study. Our asymptotic results also suggest that neglecting for the presence of non-null unbalanced parameter might lead to inconsistency. In a short application we investigated the no-arbitrage hypothesis between the volatilities of spot and futures prices.

Our results reveal that the apparent unbalance of the integration orders between the daily squared returns of the observable is misleading. An unbalanced stationary cointegration is recovered and the results are consistent with the theory as well as some empirical features found in the literature.

Appendix: Proof of theorem 1

Proof 1. Let θ be the vector of admissible parameter value, θ 0 the vector of true parameter value and S(θ

) = R m (θ) -R m (θ 0 ). Then, define the neighborhoods Θ n δ = {δ : ||δ -δ 0 || < d}, Θ n ξ = {ξ : ||ξ -ξ 0 || < e}, Θ n β = {β : ||λ δ 01 -δ 02 m (β -β 0 )|| < b} and their complements Θ c δ = Θ δ \Θ n δ , Θ c ξ = Θ ξ \Θ n ξ and Θ c β = Θ β \Θ n β such that Θ = Θ δ × Θ ξ × Θ β \Ξ.
Without loss of generality with respect to Assumption 1 we set From Robinson (1995a, p. 1634) and by the fact that

max min i ||δ i -δ 0i ||, ||ξ -ξ 0 || ≥ d, δ ∈ Θ c δ , ξ ∈ Θ c ξ , ( 20
)
so that 1/2 > d ≥ e > 0.
θ 0 ∈ Θ n δ × Θ n ξ × Θ n β it follows Pr { δ ∈ Θ c δ } ∪ { ξ ∈ Θ c ξ } ∪ { β ∈ Θ c β } = Pr inf { δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ c β } R m (θ) ≤ inf { δ∈Θ n δ }∪{ ξ∈Θ n ξ }∪{ β∈Θ n β } R m (θ) ≤ Pr inf { δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ c β } S(θ) ≤ 0 .
Accordingly, to prove the Theorem 1, it suffices to show that, as n → 0, S(θ) is positive and bounded away from 0

uniformly on { δ ∈ Θ c δ } ∪ { ξ ∈ Θ c ξ } ∪ { β ∈ Θ c β } so that Pr inf { δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ c β } S(θ) ≤ 0 → 0. ( 21
)
Now, introduce ψ 1 = δ 1δ 01 and ψ 2 = (δ 2 + ξ) -(δ 02 + ξ 0 ) and develop S(θ) as

S(θ) = log det Ĝ(θ) -2(δ 1 + δ 2 + ξ)m -1 m ∑ j=1 log λ j -log det Ĝ(θ 0 ) + 2(δ 01 + δ 02 + ξ 0 )m -1 m ∑ j=1 log λ j = log det Ĝ(θ) -log det Ĝ(θ 0 ) -m -1 m ∑ j=1 2 log λ j (δ 1 + δ 2 + ξ -δ 01 -δ 02 -ξ 0 ) = log det Ĝ(θ) -log det Ĝ(θ 0 ) -m -1 m ∑ j=1 2 log λ j p ∑ i=1 ψ i , for p = 2 so that ∑ p i=1 ψ i = (ψ 1 + ψ 2 ). Then, by the fact that m ∑ j=1 log λ j = m ∑ j=1 log(2πjn -1 ) = m ∑ j=1 log j + log(2πn -1 ) = m log(2πn -1 ) + m ∑ j=1 log j m -1 m ∑ j=1 log λ j = log(2πn -1 ) + m -1 m ∑ j=1 log j = log λ m + m -1 m ∑ j=1 log j -log m,
and rearranging S(θ), we obtain

S(θ) = log det Ĝ(θ) -log det Ĝ(θ 0 ) + log det G 0 -log det G 0 -2 log λ m p ∑ i=1 ψ i + 2 p ∑ i=1 ψ i log m -m -1 m ∑ j=1 log j + p ∑ i=1 log(2ψ i + 1) - p ∑ i=1 log(2ψ i + 1),
and finally, S(θ) = S 1 (θ) + S 2 (θ) + S 3 (θ), where

S 1 (θ) = log det Ĝ(θ) -log det G 0 -2 log λ m p ∑ i=1 ψ i + p ∑ i=1 log(2ψ i + 1) S 2 (θ) = log det Ĝ0 -log det Ĝ(θ 0 ) S 3 (θ) = 2 p ∑ i=1 ψ i log m -m -1 m ∑ j=1 log j - p ∑ i=1 log(2ψ i + 1).
The way we split S(θ) has advantage that S 2 (θ) and S 3 (θ) do not depend on β so that we can treat them following the methodology of Robinson (1995a). To summarize, we have to demonstrate the boundedness of S(θ) under

Pr inf

{ δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ c β } S(θ) ≤ 0 ≤ Pr inf { δ∈Θ c δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ c β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (22) + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ n β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (23) + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ n β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (24) + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ n β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (25) + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ c β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (26) + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ c β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 (27) + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ c β } S 1 (θ) + S 2 (θ) + S 3 (θ) ≤ 0 . ( 28
)
We first consider S 2 (θ) because it does not depend on θ and no uniform bound is needed. Because | log(1 + x)| ≤ 2|x| (see Olver et al. 2010, p. 108), it follows that for ε ≤ 1

Pr (|S 2 (θ)| ≤ ε) = Pr | log det Ĝ(θ 0 ) -log det Ĝ0 | ≤ ε ≤ Pr det Ĝ(θ 0 ) -det Ĝ0 det Ĝ0 ≤ ε/2 . Accordingly, proving that det Ĝ(θ 0 ) -det Ĝ0 p -→ 0 suffices to show that S 2 (θ) is o p (1) (see Robinson 1995a, p.
1635). To simplify notation, let

F 0 jab = G 0 ab λ -δ 0a -δ 0b j , f ab z (λ j ) = f 0 jab , A ab (e iλ j ) = A jab and I 0 jab = I 0 ab (λ j ; ϑ 02 ) with ϑ 02 = (β 0 , ξ 0 ) , G 0
ab the (a, b)-th element of G 0 and δ 0a is the ath element of ϑ 01 = (δ 01 , δ 02 + ξ 0 ) , respectively δ 0b . Evaluating Equations ( 6) and (8) at the true value, we obtain

I 0 j = I yy (λ j ) -2β 0 λ ξ 0 j I xy (λ j ) + β 2 0 λ 2ξ 0 j I xx (λ j ) I xy (λ j ) -β 0 λ ξ 0 j I xx (λ j ) I xy (λ j ) -β 0 λ ξ 0 j I xx (λ j ) I xx (λ j )
, from which we implicitly take the real part. Then,

Ĝab (θ 0 ) -G 0 ab = 1 m m ∑ j=1 I 0 jab λ -δ 0a -δ 0b j -G 0 ab = G 0 ab m m ∑ j=1 I 0 jab G 0 ab λ -δ 0a -δ 0b j -1 = G 0 ab m m ∑ j=1 I 0 jab F 0 jab -1 ,
Recall that under Assumption 3 and 4, f z (λ j ) = A(e iλ j )A(e iλ j ) * /2π and f ε (λ j ) = I 2 /2π. Now, rewrite the latter expression as,

G 0 ab m m ∑ j=1 I 0 jab F 0 jab -1 = S 21 (θ) + S 22 (θ) + S 23 (θ), S 21 (θ) = G 0 ab m m ∑ j=1 1 - F 0 jab f 0 jab I 0 jab F 0 jab , S 22 (θ) = G 0 ab m m ∑ j=1 1 f 0 jab I 0 jab -A jab I jab ε Ājab ) , S 23 (θ) = G 0 ab m m ∑ j=1 2π I jab ε -1 .
From the analysis of (Robinson 1995a(Robinson , p. 1636) and under Assumption 2-5, as m → ∞, |1 - From Robinson (1995aRobinson ( , p. 1637)), it can be demonstrated that

F 0 jab ( f 0 jab ) -1 | ≤ η, E|I 0 jab (F 0 jab ) -1 | ≤ c
E I 0 jab -A jab I jab ε Ājab = O f 0 jab log(j + 1) 1/2 j -1/2 ,
and therefore |S 22 (θ)| = o(1) as m → ∞. Finally, as n → ∞ and under Assumption 3,

2π I ε -I 2 = S 231 (θ) + S 232 (θ), S 231 (θ) = n -1 n ∑ t=1 (ε t ε t -I 2 ) p -→ 0, S 232 (θ) = n ∑ s t n ∑ t cos (s -t)λ j ε s ε t = o(1),
and S 23 (θ) is o p (1) (see Robinson 1995aRobinson , p. 1638)). Thereby, as n → ∞, we have proved that

| Ĝab (θ 0 ) -G 0 ab | = |S 21 (θ)| + |S 22 (θ)| + |S 23 (θ)| = O p η + m -1 m ∑ j=1 log j j 1/2 + o p (1) = o p (1),
and by corollary that det Ĝ(θ 0 )det Ĝ0 p -→ 0. Straightforwardly, S 2 (θ) = o p (1) follows in Equations ( 22) to (28).

Now we turn to the analysis of S 3 (θ). First, adding and subtracting 2(ψ 1 + ψ 2 ) to S 3 (θ), we have

S 3 (θ) = -2 p ∑ i=1 ψ i m -1 m ∑ j=1 log j -log m - p ∑ i=1 log(2ψ i + 1) + 2 p ∑ i=1 ψ i -2 p ∑ i=1 ψ i
Then, rearranging these terms we obtain,

S 3 (θ) = -2 p ∑ i=1 ψ i m -1 m ∑ j=1 log j -(log m -1) + 2 p ∑ i=1 ψ i - p ∑ i=1 log(2ψ i + 1). ( 29
)
From Lemma 2 of Robinson (1995a) we know that

m -1 m ∑ j=1 log j -(log m -1) = O(m -1 log m),
so that analysis of S 3 (θ) reduces to study the greatest lower bound of (29) which is of the form f (x) = xlog(1 +

x) for ψ 1 and (x + y)log(x + y + 1) for ψ 2 . Because inf x f (x) ≥ x 2 /6 and inf x,y f (x, y) ≥ (x 2 + y 2 )/6 for 0 < |x| < 1, 0 < |y| < 1, and from the restriction stated in Equation (20), we can apply the analysis of Nielsen (2007, p. 437) ütkepohl (1996, sec. 8.5.2, p. 111) and by the triangular inequality,

uniformly over { δ ∈ Θ c δ } ∪ { ξ ∈ Θ c ξ }. From L
√ 2 max(|δ 1 -δ 01 |, |δ 2 -δ 02 |) + √ 2 max(|0|, |ξ -ξ 0 |) ≥ δ 1 -δ 01 δ 2 -δ 02 + 0 ξ -ξ 0 ≥ ||Ψ|| ≥ d + e, with Ψ = (ψ 1 , ψ 2 ) . Thereby, the infimum over { δ ∈ Θ c δ } ∪ { ξ ∈ Θ c ξ } ∪ { β ∈ Θ β } of 2 ∑ p i=1 ψ i -∑ p i=1 log(2ψ i + 1) is no less than 2(d + e) √ 2 -log 1 + 2(d + e) √ 2 ≥ 2(d 2 + e 2 ) 6 .
Then, given that f (x) has a unique minimum on (-1,

∞) at x = 0 inf { δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ β } S 3 (θ) ≥ 2(d 2 + e 2 ) 6 + O(m -1 log m), inf { δ∈Θ n δ }∪{ ξ∈Θ n ξ }∪{ β∈Θ β } S 3 (θ) = o(1).
The two remaining cases are

{ δ ∈ Θ c δ } ∪ { ξ ∈ Θ n ξ } ∪ { β ∈ Θ β } and { δ ∈ Θ n δ } ∪ { ξ ∈ Θ c ξ } ∪ { β ∈ Θ β }.
In the former case, S 3 (θ) is no less than 2d 2 /6 + O(m -1 log m) while in the latter, S 3 (θ) is no less than 2e 2 /6 + O(m -1 log m) (see graphic 4).

Figure 4: Plot of (x + y)log(x + y + 1) and (x 2 + y 2 )/6 (blue and purple curves respectively). Finally we turn to the analysis of S 1 (θ). Rewrite S 1 (θ) as

S 1 (θ) = log det Ĝ(θ) -log det G 0 -2 log λ m p ∑ i=1 ψ i + p ∑ i=1 log(2ψ i + 1) = log det Ĝ(θ) -2 log λ m × [δ 1 -δ 01 + (δ 2 + ξ -δ 02 + ξ 0 )] -log det G 0 + p ∑ i=1 log(2ψ i + 1) = log det Ĝ(θ) -2 log λ m ψ 1 -2 log λ m ψ 2 -log det G 0 (2ψ 1 + 1) -1 (2ψ 2 + 1) -1 = log det V m Ĝ(θ)V m -log det G 0 (2ψ 1 + 1) -1 (2ψ 2 + 1) -1 ,
where

ψ 2 = δ 2 -δ 02 = (δ 2 + ξ) -(δ 02 + ξ 0 ), det G 0 = G 11 G 22 -G 2 12 , V m = diag(λ -ψ 1 m , λ -ψ 2 m ) and det V m Ĝ(θ)V m = det λ -2ψ 1 m Ĝ11 (θ) λ -ψ 1 -ψ 2 m Ĝ12 (θ) λ -ψ 1 -ψ 2 m Ĝ21 (θ) λ -2ψ 2 m Ĝ22 (θ) = λ -2ψ 1 -2ψ 2 m Ĝ11 (θ) Ĝ22 (θ) -λ -2ψ 1 -2ψ 2 m Ĝ21 (θ) Ĝ12 (θ) = λ -2ψ 1 -2ψ 2 m ( Ĝ11 (θ) Ĝ22 (θ) -Ĝ2 12 (θ)).
Accordingly, analysis of S 1 (θ) reduces to study

S 1 (θ) = λ -2ψ 1 -2ψ 2 m Ĝ11 (θ) Ĝ22 (θ) -Ĝ2 12 (θ) -G 11 G 22 -G 2 12 (2ψ 1 + 1) -1 (2ψ 2 + 1) -1 . Because Ĝab (θ) = m -1 ∑ m j=1 λ δ a +δ b j I jab , S 1 (θ) = S 11 (θ) + S 12 (θ) + S 13 (θ), S 11 (θ) = λ -2ψ 1 -2ψ 2 m m -1 m ∑ j=1 λ 2δ 1 j I j11 × m -1 m ∑ j=1 λ 2(δ 2 +ξ) j I j22 , S 12 (θ) = -λ -2ψ 1 -2ψ 2 m m -1 m ∑ j=1 λ δ 1 +δ 2 +ξ j Re(I j12 ) 2 , S 13 (θ) = -G 11 G 22 -G 2 12 (1 + 2ψ 1 ) -1 (1 + 2ψ 2 ) -1 .
S 11 (θ) = λ -2ψ 1 -2ψ 2 m 1 m 2 m ∑ j=1 m ∑ k=1 λ 2δ 1 j λ 2(δ 2 +ξ) k I j11 I k22 λ -2δ 01 j λ 2δ 01 j λ -2(δ 02 +ξ 0 ) k λ 2(δ 02 +ξ 0 ) k = 1 m 2 m ∑ j=1 m ∑ k=1 λ 2δ 1 j λ -2δ 01 j λ 2(δ 2 +ξ) k λ -2(δ 02 +ξ 0 ) k λ 2ψ 1 +2ψ 2 m I j11 I k22 λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k = 1 m 2 m ∑ j=1 m ∑ k=1 λ j λ m 2ψ 1 λ k λ m 2ψ 2 I j11 I k22 λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k = 1 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 k m 2ψ 2 I j11 I k22 λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k . Now, rearranging S 12 (θ) we obtain S 12 (θ) = -λ -2ψ 1 -2ψ 2 m m -1 m ∑ j=1 λ δ 1 +δ 2 +ξ j Re(I j12 ) 2 = -λ -2ψ 1 -2ψ 2 m 1 m 2 m ∑ j=1 m ∑ k=1 Re(I j12 ) Re(I k12 )λ δ 1 +δ 2 +ξ j λ δ 1 +δ 2 +ξ k = -λ -2ψ 1 -2ψ 2 m 1 m 2 m ∑ j=1 m ∑ k=1 Re(I j12 ) Re(I k12 ) × λ δ 1 +δ 2 +ξ j λ δ 1 +δ 2 +ξ k λ δ 01 +δ 02 +ξ 0 j λ -δ 01 -δ 02 -ξ 0 j λ δ 01 +δ 02 +ξ 0 k λ -δ 01 -δ 02 -ξ 0 k = - 1 m 2 m ∑ j=1 m ∑ k=1 λ δ 1 +δ 2 +ξ j λ -δ 01 -δ 02 -ξ 0 j λ δ 1 +δ 2 +ξ k λ -δ 01 -δ 02 -ξ 0 k λ 2ψ 1 +2ψ 2 m Re(I j12 ) Re(I k12 ) λ -δ 01 -δ 02 -ξ 0 j λ -δ 01 -δ 02 -ξ 0 k = - 1 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 k m ψ 1 +ψ 2 Re(I j12 ) Re(I k12 ) λ -δ 01 -δ 02 -ξ 0 j λ -δ 01 -δ 02 -ξ 0 k
.

Then, we follow [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] and correct for the fact that I 11 (λ) and I 12 (λ) are based on estimated cointegration errors. Considering I 11 (λ) -I 0 11 (λ) and I 12 (λ) -I 0 12 (λ) we obtain

I 11 (λ) = I 0 11 (λ) -2 βλ ν 0 m λ ξ 0 Re(I 0 12 (λ)) + β2 λ 2ν 0 m λ 2ξ 0 I 22 (λ), I 12 (λ) = I 0 12 (λ) -βλ ν 0 m λ ξ 0 I 22 (λ).
with βλ ν 0 m = (ββ 0 ) and ν 0 = δ 02δ 01 . Substituting in S 11 (θ) and S 12 (θ) we have

S 11 (θ) = 1 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 k m 2ψ 2 I 0 j11 I k22 λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k - 2 β m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 k m 2ψ 2 Re(I 0 j12 )I k22 λ -ν 0 m λ -ξ 0 j λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k λ -ν 0 j λ -ν 0 j + β2 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 k m 2ψ 2 I j22 I k22 λ -2ν 0 m λ -2ξ 0 j λ -2δ 01 j λ -2(δ 02 +ξ 0 ) k λ -2ν 0 j λ -2ν 0 j S 12 (θ) = - 1 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 k m ψ 1 +ψ 2 Re(I 0 j12 ) Re(I 0 k12 ) λ -δ 01 -δ 02 -ξ 0 j λ -δ 01 -δ 02 -ξ 0 k + 2 β m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 k m ψ 1 +ψ 2 I j22 Re(I 0 k12 ) λ -ν 0 m λ -ξ 0 j λ -δ 01 -δ 02 -ξ 0 j λ -δ 01 -δ 02 -ξ 0 k λ -ν 0 j λ -ν 0 j - β2 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 k m ψ 1 +ψ 2 I j22 I k22 λ -2ν 0 m λ -ξ 0 j λ δ 01 -δ 02 -ξ 0 j λ -ξ 0 k λ δ 01 -δ 02 -ξ 0 k λ -ν 0 j λ -ν k λ -ν 0 j λ -ν k , since -I j12 I k12 = -I 0 j12 I 0 k12 + 2 βλ ν 0 m λ ξ 0 j I j22 I 0 k12 -β2 λ 2ν 0 m λ ξ 0 j λ ξ 0 k I j22 I k22 . Rearranging this leads to S 11 (θ) = G 11 G 22 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 k m 2ψ 2 I 0 j11 I k22 G 11 λ -2δ 01 j G 22 λ -2(δ 02 +ξ 0 ) k -2 β G 12 G 22 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 -ν 0 k m 2ψ 2 Re(I 0 j12 )I k22 G 12 λ -δ 01 -δ 02 -ξ 0 j G 22 λ -2(δ 02 +ξ 0 ) k + β2 G 2 22 m 2 m ∑ j=1 m ∑ k=1 j m 2ψ 1 -2ν 0 k m 2ψ 2 I j22 I k22 G 2 22 λ -2(δ 02 +ξ 0 ) j λ -2(δ 02 +ξ) k S 12 (θ) = - G 2 12 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 k m ψ 1 +ψ 2 Re(I 0 j12 ) Re(I 0 k12 ) G 2 12 λ -δ 01 -δ 02 -ξ 0 j λ -δ 01 -δ 02 -ξ 0 k + 2 β G 12 G 22 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 -ν 0 k m ψ 1 +ψ 2 I j22 Re(I 0 k12 ) G 12 λ -2(δ 02 +ξ 0 ) j G 22 λ -δ 01 -δ 02 -ξ k -β2 G 2 22 m 2 m ∑ j=1 m ∑ k=1 j m ψ 1 +ψ 2 -ν 0 k m ψ 1 +ψ 2 -ν 0 I j22 I k22 G 2 22 λ -2(δ 02 -ξ 0 ) j λ -2(δ 02 -ξ 0 ) k , Now, consider the fact that m -1 ∑ m j=1 (j/m) α = (1 + α) -1 because m -1 m ∑ j=1 j m α = m ∑ j=1 j/m (j-1)/m x α dx = 1 0 x α dx, 1 0 x α dx = 1 1 + α x α+1 1 0 = (1 + α) -1 ,
Moreover, by the analysis of (Robinson 1995a(Robinson , p. 1636(Robinson -1638) ) we have

G ab m -1 m ∑ j=1 Re(I 0 jab ) G ab λ -δ 0a -δ 0b j -1 = G ab m -1 m ∑ j=1 Re(I 0 jab )G -1 ab λ δ 0a +δ 0b j -m = G ab m -1 m ∑ j=1 Re(I 0 jab )G -1 ab λ δ 0a +δ 0b j -1 = o p (1),
and thus it follows,

S 11 (θ) = G 11 G 22 (1 + 2ψ 1 ) -1 (1 + 2ψ 2 ) -1 (1 + o p (1)) -2 βG 12 G 22 (1 + 2ψ 1 -ν 0 ) -1 (1 + 2ψ 2 ) -1 (1 + o p (1)) + β2 G 2 22 (1 + 2ψ 1 -2ν 0 ) -1 (1 + 2ψ 2 ) -1 (1 + o p (1)), S 12 (θ) = -G 2 12 (1 + ψ 1 + ψ 2 ) -2 (1 + o p (1)) + 2 βG 12 G 22 (1 + ψ 1 + ψ 2 -ν 0 ) -1 (1 + ψ 1 + ψ 2 ) -1 (1 + o p (1)) -β2 G 2 22 (1 + ψ 1 + ψ 2 -ν 0 ) -2 (1 + o p (1)),
Now we can rewrite S 1 (θ) = S 11 (θ) + S 12 (θ) + S 13 (θ) as

S 1 (θ) = G 2 12 1 0 x 2ψ 1 dx 1 0 x 2ψ 2 dx - 1 0 x ψ 1 +ψ 2 2 + β2 G 2 22 1 0 x 2ψ 1 -2ν 0 dx 1 0 x 2ψ 2 dx - 1 0 x ψ 1 +ψ 2 -ν 0 2 + 2 βG 12 G 22 1 0 x ψ 1 +ψ 2 -ν 0 dx 1 0 x ψ 1 +ψ 2 dx - 1 0 x 2ψ 1 -ν 0 dx 1 0 x 2ψ 2 dx + G 11 G 22 (1 + 2ψ 1 ) -1 (1 + 2ψ 2 ) -1 -(1 + 2ψ 1 ) -1 (1 + 2ψ 2 ) -1 + o p (1).
If G 12 = 0, by the Cauchy-Schwarz inequality and the fact that ν 0 > 0 under cointegration it is straightforward that S 1 (θ)

≥ o p (1) when { δ ∈ Θ c δ } ∩ { ξ ∈ Θ c ξ } ∩ { β ∈ Θ β } and S 1 (θ) is bounded away from zero when { δ ∈ Θ n δ } ∩ { ξ ∈ Θ n ξ } ∩ { β ∈ Θ c β }. A corollary of this finding is that S 1 (θ) is also bounded when { δ ∈ Θ c δ } ∩ { ξ ∈ Θ n ξ } ∩ { β ∈ Θ c β } and { δ ∈ Θ n δ } ∩ { ξ ∈ Θ c ξ } ∩ { β ∈ Θ c β }. Similarly, S 1 (θ) ≥ o p (1) when { δ ∈ Θ n δ } ∩ { ξ ∈ Θ c ξ } ∩ { β ∈ Θ β } and { δ ∈ Θ c δ } ∩ { ξ ∈ Θ n ξ } ∩ { β ∈ Θ β }.
To summarize our results, for c > 0 an arbitrary small positive number and B = 2(d 2 + e 2 )/6 + o(1), we can rewrite the Equations ( 22) to (28) as Pr inf

{ δ∈Θ c δ }∪{ ξ∈Θ c ξ }∪{ β∈Θ c β } S(θ) ≤ 0 ≤ Pr inf { δ∈Θ c δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ c β } [S 1 (θ) > c] + o p (1) + [S 3 (θ) ≥ B] ≤ 0 + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ n β } [S 1 (θ) ≥ o p (1)] + o p (1) + [S 3 (θ) ≥ B] ≤ 0 + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ n β } [S 1 (θ) ≥ o p (1)] + o p (1) + [S 3 (θ) ≥ 2d 2 6 + o(1)] ≤ 0 + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ n β } [S 1 (θ) ≥ o p (1)] + o p (1) + [S 3 (θ) ≥ 2e 2 6 + o(1)] ≤ 0 + Pr inf { δ∈Θ c δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ c β } [S 1 (θ) > c] + o p (1) + [S 3 (θ) ≥ 2d 2 6 + o(1)] ≤ 0 + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ c ξ }∩{ β∈Θ c β } [S 1 (θ) > c] + o p (1) + [S 3 (θ) ≥ 2e 2 6 + o(1)] ≤ 0 + Pr inf { δ∈Θ n δ }∩{ ξ∈Θ n ξ }∩{ β∈Θ c β } [S 1 (θ) > c] + o p (1) + [S 3 (θ) = o(1)] ≤ 0 ,
hence showing the Equation (21) and proving the Theorem 1.

Extended Theorem 1

Now we extend the proof of the Theorem 1 to the case where G 12 = G 21 0. Given that G 12 only impact S 1 (θ) with regard to the positivity of S(θ) we only focus on the following expression

S 1 (θ) = G 2 12 1 0 x 2ψ 1 dx 1 0 x 2ψ 2 dx - 1 0 x ψ 1 +ψ 2 2 + β2 G 2 22 1 0 x 2ψ 1 -2ν 0 dx 1 0 x 2ψ 2 dx - 1 0 x ψ 1 +ψ 2 -ν 0 2 + 2 βG 12 G 22 1 0 x ψ 1 +ψ 2 -ν 0 dx 1 0 x ψ 1 +ψ 2 dx - 1 0 x 2ψ 1 -ν 0 dx 1 0 x 2ψ 2 dx + o p (1), S 1 (θ) = S 1a (θ) + S 1b (θ) + S 1c (θ) + o p (1).
Now observe that S 1a (θ) > 0 by the Cauchy-Schwarz inequality. Then, we follow [START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] and consider the fact that S 1c (θ) = 0 when ψ 1 = ψ 2 . In our unbalanced framework this is a very specific case because ψ 2 depends on ξ. Also recall that under cointegration, ν 0 > 0. Accordingly, for η > 0 and Setting to zero the first derivative of the bound we obtain that S 1 (θ) > 0 when | β| ≥ 2ε/C and no less than -ε 2 /(4C) when | β| ≤ 2ε/C (see Figure 5). Thus, for η sufficiently small, there exists ε such that 22 to 27 tend to 0 when |ψ 2ψ 1 | < η.

|ψ 2 -ψ 1 | < η, we have S 1b (θ) ≥ β2 C with C > 0 and |S 1c (θ)| ≤ | β|ε with ε > 0 so that -| β|ε ≤ S 1c (θ) ≤ | β|ε. Therefore, S 1 (θ) is no less than S 1b (θ) -| β|ε ≥ β2 C -| β|ε.

Work in Progress

Appendix: Auxiliary result

In the following, we extend the Theorem 1 of [START_REF] Robinson | Semiparametric frequency domain analysis of fractional cointegration[END_REF] and the Proposition 1 of [START_REF] Nielsen | Semiparametric Estimation in Time-Series Regression with Long-Range Dependence[END_REF] to the case of unbalanced cointegration. From Equation (8), for δ 1 , δ 2 and ξ consistently pre-estimated, we have Ĝee

(β) = Ĝ11 (β) = m -1 m ∑ j λ 2δ 1 j I j11 = m -1 m ∑ j λ 2δ 1 j I jyy + β 2 λ 2δ 1 +2ξ j I jxx -2βλ 2δ
= m -1 m ∑ j=1 λ 2(δ 1 +ξ) j Re(I jxx ) -1 m -1 m ∑ j=1 λ 2(δ 1 +ξ/2) j
Re(I jxy ) , where m -1 ∑ m j=1 λ 2δ j Re(I jab ) is the (δ)-weighted periodogram and β the narrow-band generalized least squares (NBGLS) estimate of [START_REF] Nielsen | Semiparametric Estimation in Time-Series Regression with Long-Range Dependence[END_REF].

Proposition 1. Under Assumptions 2, 3, 5 and for δ satisfying (2(δ 2 + ξ)

+ 2δ 1 -1)/4 < δ ≤ δ 1 , β -β 0 = m -1 m ∑ j=1 λ 2(δ 1 +ξ) j Re(I jxx ) -1 m -1 m ∑ j=1 λ 2δ j Re(I jxe ) = O p (n/m) δ 1 -δ 2 , ( 30
)
where e t = y tβx t .

Proof 2. From the Proposition 1 of [START_REF] Nielsen | Semiparametric Estimation in Time-Series Regression with Long-Range Dependence[END_REF] we have

m -1 m ∑ j=1 λ 2δ j Re(I jab ) = λ 0 Re f ab (µ) dµ ∼ G ab λ 1-δ a -δ b +2δ 1 -δ a -δ b + 2δ = O p λ 1-δ a -δ b +2δ m .
Therefore, by the Cauchy-Schwarz inequality

m -1 m ∑ j=1 λ 2δ j Re(I jxe ) ≤ m -1 m ∑ j=1 λ 2(δ 1 +ξ) j Re(I jxx ) × m -1 m ∑ j=1 λ 2δ 1 j Re(I jee ) 1/2 ≤ O p m n 1-2(δ 2 +ξ)+2(δ 1 +ξ) 1/2 O p m n 1-2δ 1 +2δ 1 1/2 ≤ O p m n 1-δ 2 +δ 1 .
Finally, the proof is completed by substituting this result in Equation (30):

β -β 0 = O p m n -1+2(δ 2 +ξ)-2(δ 1 +ξ) O p m n 1-δ 2 +δ 1 = O p m n δ 2 -δ 1 .

Appendix: Simulation results

As specified in the simulation results section, there is weak evidence of consistency when the parameters lie in the non-stationary regions. Indeed, the RMSE decreases very slowly as the sample size increases. Mostly in the strong cointegration case, the bias increases as the sample size increases. Interestingly, results are less impacted by non-stationarity in the weak cointegration case. This suggests that practitioners should devote a particular attention to the stationarity or non-stationarity of the data. In the latter case, they should use the estimator of Hualde (2014).

Figure 1 :Finite

 1 Figure 1: Finite sample distribution of θ for 1000 replications of the stationary cointegration model with ρ = 0 and n = 4096.

Figure 2 :Figure 3 :

 23 Figure 2: Daily squared return volatility proxy of the crude oil spot prices from January 2, 1996 to December 16, 2013.

  test statistic, T0 , requires to estimate the cointegration rank, r and the authors suggest a model selection procedure based on a tuning parameter v k = m -κ G > 0 where m G is a specific bandwidth used to estimate Ĝ and fixed to m G = n k-0.05 . Here we consider v k = {m -0to the choice of κ.[START_REF] Nielsen | Determining the cointegrating rank in nonstationary fractional systems by the exact local Whittle approach[END_REF] show that T0 d -→ χ 2 1 if r = 0 and T0 p -→ 0 otherwise. The results are reported in the Table

  and thus |S 21 (θ)| ≤ cη, with η and c any arbitrary positive numbers. Then next term to study is S 22 (θ).

  Distinguishing the two summations by indexes j and k, multiplying by G 11 G 22 /G 11 G 22 and λ S 11 (θ) we have

Figure 5 :

 5 Figure 5: Plot of | β|(| β|Cε), | β| ≥ 2ε/C and | β| ≤ 2ε/C (blue, black and purple curves respectively).

Table 1 :

 1 Simulation results for the stationary model when ξ = 0.1 and ρ = 0

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.4	0	δ 2	0.020	0.031	0.176	0.034	0.019	0.142	0.034	0.012	0.115
			δ 1	0.063	0.039	0.207	0.074	0.025	0.175	0.076	0.017	0.150
			ξ	-0.014	0.005	0.070	-0.015	0.002	0.042	-0.012	0.001	0.027
			β	0.016	0.006	0.079	0.022	0.004	0.068	0.029	0.003	0.060
	0.4	0.2	δ 2	-0.013	0.035	0.188	-0.003	0.022	0.150	-0.002	0.015	0.122
			δ 1	0.008	0.031	0.176	0.009	0.019	0.139	0.014	0.011	0.105
			ξ	-0.003	0.012	0.109	-0.010	0.006	0.076	-0.013	0.003	0.052
			β	0.006	0.019	0.139	0.025	0.012	0.113	0.041	0.007	0.092
	0.4	0.3	δ 2	-0.029	0.044	0.212	-0.021	0.028	0.170	-0.014	0.018	0.135
			δ 1	-0.006	0.030	0.175	-0.005	0.019	0.138	-0.005	0.011	0.104
			ξ	0.009	0.020	0.143	0.003	0.011	0.106	-0.004	0.006	0.079
			β	-0.025	0.035	0.189	0.002	0.022	0.149	0.024	0.013	0.118
	m = n 0.8
	0.4	0	δ 2	-0.041	0.004	0.072	-0.044	0.002	0.061	-0.052	0.001	0.062
			δ 1	0.033	0.004	0.072	0.054	0.002	0.073	0.076	0.002	0.086
			ξ	0.016	0.002	0.043	0.014	0.001	0.029	0.015	0.000	0.021
			β	0.057	0.001	0.068	0.071	0.001	0.078	0.086	0.001	0.090
	0.4	0.2	δ 2	-0.044	0.005	0.081	-0.041	0.002	0.063	-0.041	0.001	0.056
			δ 1	-0.006	0.003	0.059	0.000	0.002	0.042	0.005	0.001	0.032
			ξ	0.022	0.003	0.057	0.021	0.001	0.041	0.020	0.001	0.032
			β	0.065	0.002	0.081	0.074	0.001	0.083	0.081	0.001	0.087
	0.4	0.3	δ 2	-0.034	0.006	0.084	-0.029	0.003	0.062	-0.025	0.002	0.048
			δ 1	-0.016	0.004	0.064	-0.011	0.002	0.045	-0.007	0.001	0.034
			ξ	0.015	0.004	0.064	0.015	0.002	0.046	0.014	0.001	0.035
			β	0.049	0.004	0.079	0.055	0.002	0.073	0.059	0.002	0.071

Table 2 :

 2 Simulation results for the stationary model when ξ = 0.1 and ρ = 0.4

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.4	0	δ 2	0.038	0.026	0.167	0.045	0.015	0.131	0.038	0.009	0.101
			δ 1	0.035	0.033	0.184	0.034	0.020	0.144	0.039	0.011	0.112
			ξ	0.029	0.004	0.067	0.030	0.001	0.046	0.031	0.001	0.039
			β	0.315	0.009	0.329	0.323	0.006	0.332	0.331	0.004	0.336
	0.4	0.2	δ 2	-0.003	0.037	0.192	0.011	0.023	0.152	0.015	0.014	0.118
			δ 1	-0.006	0.032	0.180	-0.013	0.020	0.141	-0.016	0.012	0.109
			ξ	0.021	0.010	0.100	0.012	0.004	0.062	0.008	0.002	0.041
			β	0.366	0.023	0.396	0.383	0.014	0.401	0.396	0.009	0.406
	0.4	0.3	δ 2	-0.024	0.043	0.209	-0.015	0.029	0.171	0.001	0.018	0.135
			δ 1	-0.008	0.033	0.181	-0.016	0.020	0.144	-0.018	0.012	0.110
			ξ	0.020	0.015	0.124	0.012	0.009	0.094	0.001	0.004	0.064
			β	0.357	0.046	0.416	0.387	0.026	0.420	0.403	0.016	0.423
	m = n 0.8
	0.4	0	δ 2	-0.072	0.003	0.092	-0.066	0.002	0.077	-0.067	0.001	0.074
			δ 1	0.037	0.003	0.069	0.043	0.002	0.062	0.054	0.001	0.066
			ξ	0.079	0.002	0.089	0.073	0.001	0.077	0.064	0.000	0.066
			β	0.378	0.002	0.380	0.397	0.001	0.399	0.415	0.001	0.416
	0.4	0.2	δ 2	-0.061	0.003	0.083	-0.054	0.002	0.069	-0.053	0.001	0.062
			δ 1	0.001	0.003	0.054	0.002	0.002	0.039	0.002	0.001	0.029
			ξ	0.054	0.002	0.071	0.050	0.001	0.059	0.046	0.000	0.050
			β	0.424	0.003	0.427	0.432	0.002	0.434	0.440	0.001	0.441
	0.4	0.3	δ 2	-0.042	0.004	0.076	-0.037	0.002	0.058	-0.032	0.001	0.047
			δ 1	-0.010	0.003	0.057	-0.007	0.002	0.041	-0.006	0.001	0.031
			ξ	0.030	0.002	0.057	0.029	0.001	0.045	0.027	0.001	0.037
			β	0.427	0.004	0.432	0.433	0.002	0.436	0.436	0.002	0.438

Table 3 :

 3 Comparison of type I and II processes when ξ = -0.1, ρ = 0 and m = n 0.8

	Type I			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.4	0	δ 2	-0.059	0.029	0.180	-0.011	0.017	0.130	0.021	0.005	0.075
			δ 1	0.045	0.017	0.138	0.050	0.013	0.126	0.046	0.008	0.100
			ξ	0.136	0.015	0.183	0.101	0.009	0.137	0.078	0.003	0.093
			β	-0.094	0.002	0.106	-0.035	0.002	0.055	0.046	0.001	0.056
	0.4	0.2	δ 2	-0.108	0.008	0.141	-0.088	0.004	0.109	-0.074	0.002	0.088
			δ 1	-0.041	0.004	0.076	-0.037	0.002	0.061	-0.030	0.001	0.047
			ξ	0.126	0.003	0.136	0.114	0.001	0.119	0.105	0.001	0.108
			β	0.006	0.004	0.060	0.044	0.002	0.061	0.068	0.001	0.075
	0.4	0.3	δ 2	-0.020	0.006	0.081	-0.018	0.003	0.055	-0.013	0.002	0.042
			δ 1	-0.020	0.004	0.066	-0.013	0.002	0.047	-0.010	0.001	0.034
			ξ	0.006	0.002	0.049	0.007	0.001	0.033	0.006	0.001	0.024
			β	0.043	0.004	0.076	0.052	0.002	0.072	0.056	0.002	0.068
	Type II
	0.4	0	δ 2	-0.033	0.004	0.069	-0.036	0.002	0.058	-0.039	0.001	0.053
			δ 1	0.029	0.004	0.071	0.043	0.002	0.065	0.055	0.002	0.068
			ξ	0.017	0.001	0.042	0.021	0.001	0.033	0.024	0.000	0.028
			β	0.050	0.001	0.063	0.067	0.001	0.074	0.083	0.001	0.087
	0.4	0.2	δ 2	-0.035	0.005	0.078	-0.035	0.003	0.062	-0.037	0.002	0.054
			δ 1	-0.002	0.004	0.064	0.006	0.002	0.046	0.011	0.001	0.036
			ξ	0.016	0.003	0.054	0.018	0.001	0.041	0.020	0.001	0.033
			β	0.060	0.002	0.077	0.069	0.001	0.079	0.078	0.001	0.085
	0.4	0.3	δ 2	-0.030	0.006	0.083	-0.029	0.003	0.063	-0.029	0.002	0.050
			δ 1	-0.010	0.004	0.065	-0.004	0.002	0.047	0.000	0.001	0.034
			ξ	0.010	0.004	0.063	0.013	0.002	0.046	0.014	0.001	0.035
			β	0.043	0.004	0.076	0.051	0.002	0.070	0.056	0.002	0.069

Table 4 :

 4 Two-step Exact local Whittle estimates of δ for m = n k and n = 4499 = o(n 2δ/(1+2δ) ), suggested by[START_REF] Frederiksen | Fractional integration and cointegration. An overview and an empirical application[END_REF] to prevent the presence of a neglected additive noise term in the volatility proxy. More precisely, we shall consider two bandwidths: m = n k for k = {0.4, 0.5}.

	k	Spot	Futures 1	Futures 2	Futures 3	Futures 4
	0.4	0.3828	0.4473	0.4716	0.4852	0.4837
	0.5	0.6037	0.6387	0.6643	0.6596	0.6357
	0.6	0.3892	0.4155	0.4026	0.4267	0.4001
	0.7	0.2585	0.3135	0.2962	0.3114	0.2947
	bandwidth but heterogeneous across the bandwidths. Such variability is likely to appear in presence

of level shifts or short run perturbations. Accordingly, in the following we shall apply the bandwidth requirement, m

Table 5 :

 5 W-test for m = n k and ε = {0.02, 0.05}

	k	Spot Prices	Futures 1	Futures 2	Futures 3	Futures 4
	ε = 0.02
	0.4	0.6261	0.6832	0.6664	0.6642	0.6600
	0.5	0.8726	1.2638	0.8522	0.8303	0.7964
	ε = 0.05
	0.4	0.6261	0.6832	0.6664	0.6642	0.6600
	0.5	0.8726	1.2638	0.8522	0.8303	0.7964

Table 6 :

 6 Unbalance stationary cointegration analysis

	k	0.4

  Re(I jxy ) , with e t = y tβx t . Thus, R m11 (β) = log Ĝ11 (β) and the derivative with respect to β is ∂R m11 (β)/∂β =

	2m -1 ∑ m j	λ	2δ 1 +ξ j	Re(βλ	ξ j I jxx -I jxy ) . Setting this equal to 0 we obtain
	β

1 +ξ j

Table 7 :

 7 Simulation results for 10000 replications of the strong cointegration model when ξ = 0.1 and ρ = 0

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.6	0	δ 2	0.024	0.027	0.166	0.032	0.016	0.131	0.032	0.009	0.101
			δ 1	0.169	0.045	0.271	0.197	0.027	0.257	0.234	0.017	0.267
			ξ	0.002	0.003	0.055	0.001	0.001	0.036	0.001	0.001	0.027
			β	0.039	0.013	0.122	0.043	0.010	0.110	0.050	0.008	0.103
	0.8	0.2	δ 2	0.007	0.029	0.170	0.004	0.018	0.133	0.001	0.011	0.105
			δ 1	0.185	0.042	0.277	0.217	0.028	0.273	0.259	0.017	0.290
			ξ	-0.007	0.005	0.071	-0.007	0.002	0.048	-0.006	0.001	0.034
			β	-0.002	0.037	0.193	0.005	0.029	0.169	0.010	0.022	0.148
	0.9	0.4	δ 2	-0.007	0.031	0.176	-0.007	0.021	0.144	-0.006	0.014	0.117
			δ 1	0.129	0.041	0.240	0.151	0.028	0.224	0.178	0.018	0.223
			ξ	-0.026	0.008	0.095	-0.024	0.004	0.069	-0.023	0.002	0.052
			β	-0.062	0.063	0.259	-0.070	0.051	0.236	-0.067	0.042	0.215
	m = n 0.8
	0.6	0	δ 2	-0.048	0.004	0.080	-0.052	0.002	0.071	-0.063	0.001	0.072
			δ 1	0.085	0.005	0.113	0.118	0.004	0.133	0.154	0.003	0.163
			ξ	0.034	0.001	0.050	0.041	0.001	0.049	0.047	0.000	0.051
			β	0.104	0.006	0.129	0.149	0.005	0.166	0.205	0.005	0.217
	0.8	0.2	δ 2	-0.051	0.005	0.087	-0.063	0.003	0.085	-0.083	0.002	0.096
			δ 1	0.092	0.006	0.122	0.133	0.005	0.151	0.180	0.005	0.192
			ξ	0.041	0.002	0.060	0.056	0.001	0.065	0.070	0.001	0.074
			β	0.104	0.015	0.160	0.183	0.013	0.216	0.287	0.012	0.308
	0.9	0.4	δ 2	-0.057	0.006	0.095	-0.066	0.005	0.095	-0.086	0.004	0.105
			δ 1	0.067	0.006	0.104	0.089	0.005	0.115	0.120	0.005	0.139
			ξ	0.013	0.004	0.062	0.031	0.002	0.058	0.046	0.002	0.060
			β	0.027	0.036	0.191	0.092	0.036	0.212	0.182	0.037	0.265

Table 8 :

 8 Simulation results for 10000 replications of the strong cointegration model when ξ = 0.1 and ρ = 0.4

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.6	0	δ 2	0.068	0.027	0.178	0.067	0.015	0.141	0.057	0.009	0.111
			δ 1	0.157	0.039	0.252	0.176	0.022	0.229	0.209	0.013	0.239
			ξ	0.028	0.003	0.065	0.026	0.001	0.045	0.024	0.001	0.034
			β	0.175	0.019	0.222	0.183	0.014	0.217	0.191	0.010	0.217
	0.8	0.2	δ 2	0.041	0.029	0.175	0.035	0.018	0.140	0.025	0.012	0.111
			δ 1	0.157	0.043	0.260	0.185	0.027	0.247	0.226	0.017	0.261
			ξ	0.006	0.006	0.075	0.007	0.002	0.050	0.008	0.001	0.035
			β	0.075	0.046	0.227	0.089	0.036	0.210	0.099	0.027	0.192
	0.9	0.4	δ 2	0.010	0.033	0.181	0.009	0.021	0.145	0.004	0.015	0.121
			δ 1	0.111	0.046	0.241	0.125	0.030	0.213	0.152	0.019	0.204
			ξ	-0.014	0.008	0.093	-0.012	0.004	0.068	-0.009	0.002	0.050
			β	0.032	0.080	0.284	0.031	0.067	0.260	0.039	0.055	0.238
	m = n 0.8
	0.6	0	δ 2	-0.067	0.004	0.091	-0.066	0.002	0.079	-0.068	0.001	0.074
			δ 1	0.108	0.005	0.130	0.134	0.004	0.148	0.159	0.003	0.167
			ξ	0.086	0.001	0.094	0.090	0.001	0.094	0.089	0.001	0.092
			β	0.295	0.007	0.306	0.361	0.006	0.369	0.433	0.005	0.439
	0.8	0.2	δ 2	-0.068	0.005	0.098	-0.077	0.003	0.094	-0.093	0.002	0.102
			δ 1	0.116	0.007	0.143	0.156	0.005	0.172	0.200	0.005	0.212
			ξ	0.082	0.002	0.094	0.096	0.001	0.101	0.108	0.001	0.111
			β	0.255	0.018	0.289	0.367	0.016	0.388	0.506	0.015	0.520
	0.9	0.4	δ 2	-0.077	0.006	0.108	-0.084	0.004	0.106	-0.102	0.003	0.116
			δ 1	0.082	0.007	0.118	0.104	0.006	0.128	0.135	0.005	0.152
			ξ	0.051	0.004	0.082	0.066	0.003	0.083	0.080	0.002	0.090
			β	0.191	0.050	0.294	0.279	0.050	0.358	0.399	0.051	0.459

Table 9 :

 9 Simulation results for 10000 replications of the weak cointegration model when ξ = 0.1 and ρ = 0

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.6	0.4	δ 2	-0.017	0.034	0.185	-0.002	0.022	0.147	0.002	0.013	0.116
			δ 1	0.003	0.031	0.176	0.003	0.018	0.135	0.003	0.010	0.102
			ξ	0.013	0.015	0.123	0.006	0.007	0.086	0.004	0.004	0.062
			β	0.094	0.074	0.288	0.117	0.053	0.258	0.138	0.037	0.236
	0.8	0.6	δ 2	-0.022	0.033	0.184	-0.010	0.021	0.144	-0.008	0.013	0.113
			δ 1	0.010	0.031	0.176	0.008	0.020	0.141	0.006	0.012	0.109
			ξ	0.020	0.021	0.146	0.013	0.011	0.107	0.012	0.006	0.079
			β	0.139	0.184	0.451	0.157	0.137	0.402	0.166	0.098	0.354
	0.9	0.8	δ 2	-0.032	0.038	0.197	-0.025	0.027	0.165	-0.015	0.016	0.128
			δ 1	-0.000	0.029	0.171	-0.001	0.019	0.139	0.005	0.012	0.111
			ξ	0.009	0.033	0.182	0.004	0.022	0.148	-0.003	0.012	0.109
			β	0.107	0.388	0.632	0.083	0.271	0.527	0.064	0.196	0.447
	m = n 0.8
	0.6	0.4	δ 2	-0.052	0.004	0.085	-0.048	0.002	0.067	-0.047	0.001	0.059
			δ 1	-0.013	0.004	0.061	-0.007	0.002	0.044	-0.002	0.001	0.032
			ξ	0.035	0.003	0.065	0.037	0.001	0.054	0.038	0.001	0.048
			β	0.179	0.016	0.218	0.203	0.012	0.230	0.224	0.010	0.245
	0.8	0.6	δ 2	-0.061	0.004	0.090	-0.059	0.002	0.077	-0.062	0.001	0.073
			δ 1	-0.003	0.004	0.061	0.007	0.002	0.047	0.014	0.001	0.040
			ξ	0.029	0.003	0.062	0.036	0.002	0.054	0.043	0.001	0.053
			β	0.239	0.042	0.314	0.281	0.033	0.334	0.332	0.028	0.372
	0.9	0.8	δ 2	-0.046	0.006	0.090	-0.036	0.003	0.066	-0.031	0.002	0.053
			δ 1	-0.014	0.004	0.064	-0.005	0.002	0.049	0.002	0.001	0.038
			ξ	-0.014	0.005	0.069	-0.010	0.002	0.050	-0.005	0.001	0.038
			β	0.058	0.093	0.311	0.072	0.078	0.288	0.091	0.073	0.286
	cesses. Journal of Econometrics 167, 426-447.

Table 10 :

 10 Simulation results for 10000 replications of the weak cointegration model when ξ = 0.1 and ρ = 0.4

	m = n 0.5			256			512			1024
	δ 2	δ 1	θ	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	0.6	0.4	δ 2	-0.007	0.034	0.185	0.009	0.022	0.147	0.014	0.013	0.114
			δ 1	-0.003	0.032	0.179	-0.011	0.020	0.141	-0.011	0.011	0.107
			ξ	0.030	0.011	0.111	0.022	0.005	0.076	0.019	0.002	0.053
			β	0.413	0.091	0.511	0.434	0.064	0.502	0.455	0.043	0.500
	0.8	0.6	δ 2	-0.019	0.034	0.184	-0.008	0.022	0.149	-0.001	0.013	0.113
			δ 1	0.002	0.033	0.182	-0.000	0.021	0.145	-0.005	0.012	0.109
			ξ	0.029	0.016	0.128	0.024	0.008	0.094	0.021	0.004	0.068
			β	0.407	0.221	0.622	0.421	0.158	0.579	0.435	0.115	0.552
	0.9	0.8	δ 2	-0.037	0.039	0.200	-0.025	0.024	0.158	-0.015	0.017	0.131
			δ 1	-0.001	0.032	0.180	-0.007	0.021	0.143	-0.005	0.013	0.116
			ξ	0.012	0.027	0.166	0.005	0.015	0.123	0.000	0.009	0.094
			β	0.374	0.454	0.771	0.350	0.325	0.669	0.336	0.225	0.581
	m = n 0.8
	0.6	0.4	δ 2	-0.064	0.004	0.088	-0.057	0.002	0.071	-0.053	0.001	0.061
			δ 1	-0.007	0.003	0.056	-0.005	0.002	0.041	-0.004	0.001	0.031
			ξ	0.065	0.002	0.081	0.064	0.001	0.073	0.061	0.001	0.067
			β	0.502	0.017	0.518	0.525	0.012	0.536	0.544	0.010	0.552
	0.8	0.6	δ 2	-0.071	0.004	0.094	-0.065	0.002	0.079	-0.064	0.001	0.073
			δ 1	-0.004	0.004	0.059	0.003	0.002	0.045	0.007	0.001	0.038
			ξ	0.059	0.003	0.079	0.064	0.001	0.074	0.068	0.001	0.074
			β	0.523	0.047	0.567	0.569	0.036	0.599	0.623	0.030	0.646
	0.9	0.8	δ 2	-0.053	0.005	0.089	-0.042	0.003	0.068	-0.035	0.002	0.053
			δ 1	-0.015	0.004	0.066	-0.007	0.003	0.051	-0.003	0.002	0.039
			ξ	-0.001	0.004	0.061	0.002	0.002	0.046	0.006	0.001	0.036
			β	0.339	0.104	0.468	0.342	0.088	0.452	0.365	0.081	0.462

See Granger (2010) for a very simple and short introduction to the cointegration theory.

See also[START_REF] Dueker | Maximum-likelihood estimation of fractional cointegration with an application to US and Canadian bond rates[END_REF][START_REF] Dueker | Maximum-likelihood estimation of fractional cointegration with an application to US and Canadian bond rates[END_REF][START_REF] De Truchis | Approximate Whittle Analysis of Fractional Cointegration and the Stock Market Synchronization Issue[END_REF] for a fully parametric approach in time and frequency domain respectively.[START_REF] Hualde | Semiparametric inference in multivariate fractionally cointegrated systems[END_REF] and[START_REF] Johansen | Likelihood Inference for a Fractionally Cointegrated Vector Autoregressive Model[END_REF] deal with extensions to multivariate case in frequency and time domain respectively, but this discussion goes beyond of the scope of the paper.

[START_REF] Hualde | Estimation of long-run parameters in unbalanced cointegration[END_REF] interestingly finds that the limiting distributions of estimates depend on a modified version of the Type II fractional Brownian motion. Different properties of this new type of fractional Brownian motion are discussed in[START_REF] Hualde | Weak convergence to a modified fractional Brownian motion[END_REF].

[START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF] propose a local Whittle treatment of multivariate stationary systems with unknown phase, considering the model of[START_REF] Nielsen | Local Whittle Analysis of Stationary Fractional Cointegration and the Implied-Realized Volatility Relation[END_REF] as a special case.

The presence of non-null off-diagonal elements in G should imply non-negligible imaginary part of the cross-spectrum element f ab z (λ) such as f ab z (λ) ∼ G ab λ -δa -δ b e i(π-λ)(δa -δ b )/2 as λ → 0, for a, b = 1, 2 and where G ab denotes the (a, b)th element of G.

We aim to investigate the asymptotic distribution of the estimator but for now the paper only discusses the consistency.

To generate nonstationary series (e.g. x t with δ 2 ≥ 1/2), we simulate an intermediate stationary (or asymptotically stationary) process, integrated of order I(δ ∆ = δ 2 -1) and cumulate the resulting series.

We employ this rough approximation because we do not have the data but we planned to solve this issue before submitting the paper.

All data were collected on the website of the US Energy Information Administration: http://www.eia.gov/petroleum/.