1. MOTIVATION
Stationary long-memory processes have extensive applications in economics and finance, particularly with regard to modeling financial variables such as volatility and trading volume. The autocovariances of such processes decay according to a power law, and the spectra are undefined at the origin. Correspondingly, conventional formulas and estimation procedures for long-run variance matrices that apply under weak dependence are no longer relevant under long-range dependence. Nonetheless, some modified versions of these (typically infinite-dimensional) quantities do exist and are useful in the development of asymptotics involving long-memory time series, for instance, in the estimation of fractionally cointegrated systems (Kim and Phillips, 1999; Robinson and Hualde, 2003; Velasco, 2003; Davidson, 2004; Davidson and Hashimzade, 2008). Henry and Zaffaroni (2003) provide a recent survey of the many applications of fractional integration and long-range dependence in macroeconomics and finance.
This note shows how to define long-run covariance matrices for general multivariate fractionally integrated processes, and we focus here on cases of long-range dependence, some results on overdifferenced processes being given in earlier work (Phillips, 1995, Lem. 8.1). We first develop a general form of asymptotic expansion for the autocovariance matrix of such a multivariate process. The approximation induced by this expansion is of independent interest. It has a simple form that gives the power law decay structure of the elements of the autocovariance matrix and reveals an interesting asymmetric structure for the cross autocovariances. In the scalar case, the result reduces to a formula obtained recently in Lieberman and Phillips (2006). The expansion is particularly useful in developing a limiting form of a standardized sum of the autocovariance matrices, allowing us to define the long-run variance matrix of a multivariate fractional process. In contrast to the autocovariance function, the long-run variance matrix has a simple symmetric form that depends on the long-memory parameters of the constituent processes and the long-run variance matrix of the short-memory components.
2. RESULTS
Let Xt be a real-valued covariance stationary m-vector time series generated by the system

where ut = (u1t,…,umt)′ is a covariance stationary process whose spectral density matrix fuu(λ) is assumed to be continuously differentiable on [−π,π] and bounded away from zero (in the sense of positive definite matrices) at the zero frequency λ = 0. The smoothness condition on fuu(λ) is needed to develop an asymptotic expansion of the autocovariance function defined by a Fourier integral inversion of fuu(λ).
Here Xt is a multivariate extension of a scalar fractionally integrated process (the so-called I(d) process), and each component Xat exhibits long-range dependence whenever da > 0. The time series Xt reduces to a multivariate autoregressive fractionally integrated moving average (ARFIMA) process when ut is a vector autoregressive moving average (ARMA) process, but the specification (1) does not require ut to be of this or any other parametric form.
Let fxx(λ) denote the spectral density of Xt, so that the autocovariance matrix is given by

Define

and then the spectral density of Xt satisfies (e.g., Hannan, 1970, p. 61)

where the affix * signifies complex conjugate transpose. As is well known, the memory parameters, da, govern the long-run dynamics of Xt and the behavior of its spectrum fxx(λ) around the origin. Often, when attention is focused on long-run dynamics, it is useful to specify the spectral density only locally in the vicinity of the origin and to avoid short-run dynamic specifications concerning ut altogether. In the multivariate case, we also are interested in the behavior of cross spectra at the origin and the corresponding cross autocovariances at long lags, in addition to the individual spectra and autocovariance functions.
A first-order approximation to the behavior of fxx(λ) at the origin is easily seen to be given by

and higher order approximations may be similarly obtained (e.g., Phillips and Shimotsu, 2004; Shimotsu, 2006). The factors involving the complex exponentials eiπda/2 turn out to be important in the off diagonal elements of fxx(λ), and these figure in the analysis that follows. When 0 < da < ½, the individual time series Xat have long memory, and the j-lag autocovariances decrease slowly, according to the power law j2da−1 as j → ∞. In this case, the autocovariances are not summable, and the usual formula for the long-run variance of Xat is undefined. However, as shown subsequently, upon suitable standardization, we may define the long-run variances and covariances of the elements of Xt.
We start with the following result, which gives an asymptotic approximation to the autocovariance matrix function for long lags. Lieberman and Phillips (2006) gave a complete asymptotic series expansion of the autocovariance function of a scalar long-memory time series. Under stronger smoothness conditions on the spectrum fuu(λ), similar asymptotic series may be developed here. The theorem that follows gives the leading term in the corresponding expansion for the multivariate case, which is sufficient for the present purpose of developing a formula for the long-run variance matrix. The proof is given in the Appendix.
THEOREM 1. If da ∈ [0,½) for all a = 1,…, m, and the spectral density matrix fuu(λ) of ut is continuously differentiable, then

Remark 1. Note that the asymptotic approximation (4) is asymmetric. Suppose for example that da < db and fuaub(0) > 0. Then, for large k

because sin{πda} < sin{πdb}. In particular, when da = 0 < db < ½, sin{πda} = 0, and we have the interesting phenomenon that γba(k) decays faster than the power law kdb−1, corresponding to the short-memory property of Xat, whereas γab(k) decays according to the power law kdb−1 as k → ∞, corresponding to the long-memory property of Xbt. This asymmetry is explained by the fact that γab(k) is dominated by the slow decay in the impulse responses affecting Xbt+k, whereas the impulse responses and autocovariances of Xat+k decay faster than any power rate,1 thereby determining the different behavior of γba(k) when 0 = da < db < ½.
Remark 2. When da = db, we have

corresponding to the leading term given in the asymptotic expansion of the autocovariance function for the scalar case in Lieberman and Phillips (2006).
We now define the standardization matrix Dn = diag(nd1,…,ndm) and the partial sum

and let d = mina≤m da for the following theorem, whose proof is in the Appendix.
THEOREM 2. If da ∈ (0,½) for all a = 1,…, m, and the spectral density matrix fuu(λ) of ut is continuously differentiable, then as n → ∞

Remark 3. The diagonal elements of (6) are

and, as da → 0, this formula tends to 2πfuu(0) = lrvar(ut), the limiting variance of the standardized partial sum

.
Remark 4. In the scalar case with ut ∼ iid(0,σu2), Sowell (1990) showed that

We may compare this formula with (7). By the reflection formula for the gamma function we have Γ(d)Γ(1 − d) = π/sin(πd), so that (7) may be written in the alternate form

which clearly reduces to (8) in the case of iid ut.
Remark 5. Formula (7) for the asymptotic variance corresponds to that delivered by the covariance kernel of the limiting fractional Brownian motion. In particular, it is well known (e.g., Chan and Terrin, 1995; Marinucci and Robinson, 2000) that under certain regularity conditions we have the weak convergence

It is often convenient to define the limiting fractional Brownian motions Bda(r) in terms of their harmonizable representations (see Samorodnitsky and Taqqu, 1994; Chan and Terrin, 1995; Davidson and Hashimzade, 2008) as follows:

where {Wa(λ) : a = 1,…,m} are complex-valued Gaussian random measures satisfying



for λ ∈ [−π,π], where ωab = 2πfuaua(0) and where the bar denotes complex conjugation. Observe that

using Gradshteyn and Ryzhik (2000, formula 3.823) in (11); cf. Davidson and Hashimzade (2008). Applying the reflection formula Γ(1 − z)Γ(z) = π/sin(πz) with z = 2 + 2da so that 1 − z = −1 − 2da, we have
using sin(2πda) = 2 sin πda cos πda. It now follows from (12) and (13) that
where we use the reflection formula again in the form Γ(1 − da)Γ(da) = π/sin(πda), leading to the stated correspondence with (7) and (9). Expression (12) was also obtained in Davidson and Hashimzade (2008, formula 2.6 with κ = 1). Similar arguments show that

corresponding to (6).
3. DISCUSSION AND APPLICATION
Some estimation procedures such as fully modified estimation in a fractional cointegration model (Kim and Phillips, 1999; Davidson, 2004) involve unknown long-run variance and covariance matrices for fractional processes such as those given in Theorem 2, which need to be estimated consistently for these procedures to be implemented. Consistent estimation of these long-run covariances can be accomplished by a stepwise procedure that involves separate consistent estimation of the memory parameters and the long-run variance matrix, Ωuu = 2πfuu(0) of ut. The memory parameters {da : a = 1,…,m} can be estimated by any consistent semiparametric method, such as the local Whittle (Robinson, 1995) or exact local Whittle (Shimotsu and Phillips, 2005) procedures. Using estimates

obtained in this way, estimates of the residuals uat can be constructed by the truncated filtering operation

where

. Using

, the long-run variance matrix Ωuu of ut may then be consistently estimated by any conventional heteroskedasticity and autocorrelation consistent (HAC) procedure.
2 These estimates of Ωuu and da may then be plugged into formulas such as those given in Theorem 2 to produce consistent estimates of the required long-run covariances of the fractional processes.APPENDIX: Proofs
Proof of Theorem 1. The derivation of (4) uses an asymptotic expansion of the Fourier inversion formula for Γxx(k), which can be written as

where

In what follows, we will work with the abth element of Γxx(k) denoted by

Observe that |1 − eiw| = |2 sin(w/2)| and

so that

and then
It follows that for 0 ≤ w < π
and for −π < w ≤ 0
Hence,
where Fuaub(w) = ((2 sin(w/2))/w)−da−dbfuaub(w)e−(iw(da−db))/2 ∈ C1[−π,π].
When at least one of da or db > 0, the two integrals in (A.2) have critical points (singularities in the integrand) at the origin w = 0. For Fourier integrals of this type asymptotic expansions for large k were originally developed by Erdélyi (1956), a convenient reference being Bleistein and Handelsman (1986). In particular, if F(w) ∈ C∞[a,b], and α and β are not integers, then Erdélyi's result implies that the integral

has the following complete asymptotic series representation as k → ∞:

where

and

These expansions hold to the first term n = 0 provided F(w) ∈ C1[a,b], and this degree of smoothness is all that is required for the present application.
Specializing the expansion of (A.3) to the present case, we first consider the integral

and set a = 0, b = π, β = 1, and α = 1 − da − db in formula (A.4). We deduce that
Using the same settings a = 0, b = π, β = 1, and α = 1 − da − db again, we next find that
Combining results (A.6) and (A.7) in (A.2) we obtain
which gives the stated result because Fuaub(0) = fuaub(0). █
Proof of Theorem 2. From Theorem 1, as k → ∞

Correspondingly, as k → −∞, we have, because γab(k) = γba(−k),

Then,
For some L such that (1/L) + (L/n2d) → 0 as n → ∞, we decompose the first sum in (A.10) as follows:

The abth element of this matrix is
by Euler summation. It follows that

Next consider the second sum in (A.10). The abth element is
by Euler summation again.
We now combine (A.11) and (A.12) in (A.10) giving, as n → ∞,
which corresponds to the stated result. █
