Large-sample approximations for variance-covariance matrices of high-dimensional time series

TL;DR

This paper develops large-sample approximations for the variance-covariance matrices of high-dimensional time series, enabling better inference and analysis in complex, high-dimensional data scenarios.

Contribution

It introduces strong Brownian motion approximations for bilinear forms of sample covariance matrices in high-dimensional time series, covering dependent processes and various statistical applications.

Findings

Valid for high-dimensional, possibly dependent, time series.

Applicable to sparse principal components and portfolio selection.

Supports inference in big data contexts.

Abstract

Distributional approximations of (bi--) linear functions of sample variance-covariance matrices play a critical role to analyze vector time series, as they are needed for various purposes, especially to draw inference on the dependence structure in terms of second moments and to analyze projections onto lower dimensional spaces as those generated by principal components. This particularly applies to the high-dimensional case, where the dimension $d$ is allowed to grow with the sample size $n$ and may even be larger than $n$. We establish large-sample approximations for such bilinear forms related to the sample variance-covariance matrix of a high-dimensional vector time series in terms of strong approximations by Brownian motions. The results cover weakly dependent as well as many long-range dependent linear processes and are valid for uniformly $ \ell_1 $-bounded projection vectors, which arise, either naturally or by construction, in many statistical problems extensively studied for high-dimensional series. Among those problems are sparse financial portfolio selection, sparse principal components, the LASSO, shrinkage estimation and change-point analysis for high--dimensional time series, which matter for the analysis of big data and are discussed in greater detail.

Equations

$$Y_1^{(\nu )},\dots ,Y_n^{(\nu )},\nu =1,\dots ,d_n,$$

$${\mathbf{Y}}_{ni}=(Y_i^{(1)},\dots ,Y_i^{(d_n)}{)}^{\prime },1\le i\le n,$$

$${\mathbf{Y}}_n=(Y^{(1)},\dots ,Y^{(d_n)}{)}^{\prime }$$

$${\mathbf{\Sigma }}_n=E({\mathbf{Y}}_n{\mathbf{Y}}_n^{\prime })={\left(E(Y^{(\nu )}{Y}^{(\mu )})\right)}_{1\le \nu ,\mu \le d_n}$$

$${\hat{\mathbf{\Sigma }}}_n=\frac{1}{n}\sum _{i=1}^n{\mathbf{Y}}_{ni}{\mathbf{Y}}_{ni}^{\prime }$$

$${\hat{\sigma }}_{\nu \mu }=\frac{1}{n}\sum _{i=1}^n{Y}_i^{(\nu )}{Y}_i^{(\mu )},\nu ,\mu =1,\dots ,d_n.$$

$${\mathbf{Y}}_{ni}=(Z_i,Z_{i+1},\dots ,Z_{i+d_n-1}{)}^{\prime },i=1,\dots ,n.$$

$${\hat{\sigma }}_{\nu \mu }=\frac{1}{n}\sum _{i=1}^n{Z}_{i+\nu -1}{Z}_{i+\mu -1},1\le \nu ,\mu \le d_n.$$

$${\tilde{\gamma }}_Z(h)=\frac{T}{T-h}{\hat{\gamma }}_Z(h),{\hat{\gamma }}_Z(h)=\frac{1}{T}\sum _{i=1}^{T-h}{Z}_i{Z}_{i+h},$$

$${\mathbf{w}}_n=(w_1,\dots ,w_{d_n}{)}^{\prime },n\ge 1,$$

$$\underset{n\in \mathbb{N}}{\sup }\Vert {\mathbf{w}}_n{\Vert }_{{\ell }_1}=\underset{n\in \mathbb{N}}{\sup }\sum _{\nu =1}^{d_n}|w_{\nu }|<\infty$$

$$Q_n({\mathbf{v}}_n,{\mathbf{w}}_n)={\mathbf{v}}_n^{\prime }{\hat{\mathbf{\Sigma }}}_n{\mathbf{w}}_n$$

$$|Q_n({\mathbf{v}}_n,{\mathbf{w}}_n)|=\sigma |{\mathbf{v}}_n^{\prime }11^{\prime }{\mathbf{w}}_n|=\sigma |\sum _i{v}_{ni}\sum _i{w}_{ni}|\le \sigma \Vert {\mathbf{v}}_n{\Vert }_{{\ell }_1}\Vert {\mathbf{w}}_n{\Vert }_{{\ell }_1}.$$

$${\mathcal{U}}_{\delta }=\{\{{\mathbf{v}}_n\}:\underset{n\in \mathbb{N}}{\sup }\Vert {\mathbf{v}}_n{\Vert }_{{\ell }_1}\le \delta \},$$

$${\sigma }_k^2=E({ϵ}_k^2)$$

$$\underset{k}{\sup }E|{ϵ}_k{|}^{4+\delta }<\infty ,$$

$$Y_k^{(\nu )}=Y_{nk}^{(\nu )}=\sum _{j=0}^{\infty }{c}_{nj}^{(\nu )}{ϵ}_{k-j},k=1,\dots ,n,$$

$$\text{mbox}Cov(Y_{nt}^{(\nu )},Y_{nt}^{(\mu )})$$

$$\text{mbox}Cov(Y_{nt}^{(\nu )},Y_{n,t+h}^{(\mu )})$$

$${\mathbf{\Sigma }}_n[t]={\mathbf{C}}_n\Lambda {\mathbf{C}}_n^{\prime }=\sum _{j=0}^{\infty }{\sigma }_{t-j}^2{\mathbf{c}}_{nj}{\mathbf{c}}_{nj}^{\prime },\Lambda =\mathrm{diag}({\sigma }_0^2,{\sigma }_1^2,\dots ),$$

$${\mathbf{\Sigma }}_n(h)=E({\mathbf{Y}}_{nt}{\mathbf{Y}}_{n,t+h}^{\prime })={\mathbf{C}}_n\Lambda (L^{-h}{\mathbf{C}}_n{)}^{\prime }$$

$$\underset{n\in \mathbb{N}}{\sup }\underset{1\le \nu \le d_n}{\max }|c_{nj}^{(\nu )}{|}^2<<(j\vee 1{)}^{-3/2-\theta }$$

$$(1-L{)}^d{X}_t={ϵ}_t,$$

$${\mathbf{Y}}_n$$

$${\mathbf{\Sigma }}_n$$

$$c_j^w=\sum _{\nu =1}^{d_n}{w}_{\nu }{c}_j^{(\nu )},j\ge 0,$$

$$f_{0,j}^{(n)}=f_{0,j}^{(n)}({\mathbf{v}}_n,{\mathbf{w}}_n)=\sum _{\nu ,\mu =1}^{d_n}{v}_{\nu }{w}_{\mu }{c}_j^{(\nu )}{c}_j^{(\mu )},j=0,1,\dots ,$$

$$f_{l,j}^{(n)}=f_{l,j}^{(n)}({\mathbf{v}}_n,{\mathbf{w}}_n)=\sum _{\nu ,\mu =1}^{d_n}{v}_{\nu }{w}_{\mu }[c_j^{(\nu )}{c}_{j+l}^{(\mu )}+c_j^{(\mu )}{c}_{j+l}^{(\nu )}],l=1,2,\dots ;j=0,1,\dots ,$$

$${\tilde{f}}_{l,i}^{(n)}={\tilde{f}}_{l,i}^{(n)}({\mathbf{v}}_n,{\mathbf{w}}_n)=\sum _{j=i}^{\infty }{f}_{l,j}^{(n)}=\sum _{j=i}^{\infty }\sum _{\nu ,\mu =1}^{d_n}{v}_{\nu }{w}_{\mu }[c_j^{(\nu )}{c}_{j+l}^{(\mu )}+c_j^{(\mu )}{c}_{j+l}^{(\nu )}],l,i=0,1,\dots .$$

$${\tilde{f}}_{l,i}^{(n)}(s_1{\mathbf{v}}_n,s_2{\mathbf{w}}_n)=s_1{s}_2{\tilde{f}}_{l,i}^{(n)}({\mathbf{v}}_n,{\mathbf{w}}_n)$$

$$|f_{0,j}^{(n)}({\mathbf{v}}_n,{\mathbf{w}}_n)|$$

Full text

Keywords: Brownian motion Central limit theorem Change-points Long memory Multivariate analysis Principal component analysis Portfolio analysis Strong approximation Time Series

1. Introduction

The estimation of high-dimensional variance-covariance matrices based on a vector time series arises in diverse areas such as financial portfolio optimization, image analysis and multivariate time series analysis in general. Of particular interest is the case that the dimension $d=d_{n}$ of the time series grows even faster than the sample size $n$. Due to the lack of consistency of the sample variance-covariance estimator with respect to commonly used norms such as the Frobenius norm, various regularized modifications have been proposed and extensively studied within a high-dimensional context. For example, banding and tapering estimators, recently studied by Bickel and Levina (2008) for Gaussian samples, may achieve consistency if $\log d_{n}/n=o(1)$. Chen et al. (2013) establish bounds for thresholded sample covariance estimators for a high-dimensional vector time series in terms of scaled Frobenius and spectral norms, allowing for non-stationarity and dependence. The performance of shrinkage estimation, a widely used technique dating back to the seminal work of Stein (1956), has been investigated by Ledoit and Wolf (2003) for i.i.d. samples of growing dimension and further studied for the weak dependent case in Sancetta (2008). For results on shrinkage estimation in the frequency domain we refer to Böhm and von Sachs (2009).

However, often the estimation of the $d_{n}^{2}$-dimensional variance-covariance matrix is an intermediate step and one is mainly interested in the behavior of functions of the sample variance-covariance matrix, especially quadratic and bilinear forms which naturally arise when studying projection type statistics. In addition, one often needs distributional approximations of such functions, in order to construct statistical decision procedures. Whereas consistency and performance properties have been already investigated to some extent (see above), there are only a few results about the asymptotic distribution theory in the sense of distributional convergence (weak convergence) and strong approximations by Brownian motions, respectively, going beyond the classical results. Our framework allows to embed autocovariances and, in particular, cross-autocovariances as well. To the best of our knowledge, there are no weak convergence results addressing those issues within a high-dimensional framework, i.e. for a large number of correlated time series. In some sense close to the present paper are the following results for fixed dimension and autocovariance matrices.. Wu and Min (2005) derived a CLT for a finite number of sample auto-covariances assuming a linear process. For a one-dimensional context, Xiao and Wu (2014) study a central limit theorem (CLT), Portmanteau tests and simultaneous inference for a growing number of lags based on a Gumbel type extreme value theory, see also the review Wu and Xiao (2011) and Jirak (2011). In Wu et al. (2010) the estimation of autocovariances for long memory linear processes has been discussed and studied in depth including the case of a finite number of lags starting at a large lag $k_{n}$ with $k_{n}/n=o(1)$. Kouritzin (1995), also working within a linear process framework, established large-sample distributional asymptotics, based on strong approximations, of the sample cross-covariance matrix for two time series. His assumptions are weak enough to cover the case of long-range dependence as well.

The present paper builds upon the latter result by establishing strong approximations of bilinear forms associated to the centered sample covariance matrix of a high-dimensional vector time series. The result implies the validity of a central limit theorem (CLT) for scaled bilinear forms $\sqrt{n}{\mathbf{v}}_{n}^{\prime}[\widehat{\boldsymbol{\Sigma}}_{n}-E(\widehat{\boldsymbol{\Sigma}})]{\mathbf{w}}_{n}$, where $\widehat{\boldsymbol{\Sigma}}_{n}$ is the usual sample variance-covariance matrix and ${\mathbf{v}}_{n}$, ${\mathbf{w}}_{n}$ are weighting vectors. It turns out that $d_{n}$ may even grow faster than $n$.

Concerning the weighting (or projection) vectors, our results assume that they are uniformly bounded in the $\ell_{1}$-norm. Such projections naturally arise in many problems studied in the area of high-dimensional statistics and probability: Sparse optimal portfolio selection, as recently studied by Brodie et al. (2009), deals with explicit construction of $\ell_{1}$-bounded portfolios from historical data sets. The same applies to several approaches of sparse principal component analysis, especially those of Jolliffe et al. (2003), Shen and Huang (2008) and Witten and Tibshirani (2008), where $\ell_{1}$-bounded principal components are constructed, in order to represent high-dimensional data by only a few sparse projections. We discuss those applications in greater detail in Section 5. We also illustrate how the results can be applied to obtain distributional approximations of shrinkage estimators of a high-dimensional covariance matrix. Lastly, we discuss the application to detect the presence of a change-point. Such procedures analyze the data to identify changes in the distribution and have been thoroughly studied for various second order problems for time series. Of course, a change in a covariance $\gamma_{X}(i,j)=E(X^{(i)}X^{(j)})$, $t\geq 1$, of a vector time series can be analyzed by applying any method which is sensitive to location changes to the sequence $X_{t}^{(i)}X_{t}^{(j)}$, $t\geq 1$. Such methods are discussed in Hušková and Hlávka (2012), Aue and Horváth (2013). Detectors based on local linear estimators have been proposed by Steland (2010) and kernel detectors were studied by Steland (2004) and Steland (2005). For methods based on characteristic functions we refer to Steland and Rafajłowicz (2014) and the references therein.

The organization of the paper is as follows. Section 2 explains the general setting, discusses its basic relationship to projection-based analyses and introduces the bilinear form of interest. Notation, the specific model for the vector time series and its interpretation in terms of an infinite–dimensional latent factor model as well as assumptions are introduced and discussed in Section 3. Section 4 provides several results on strong approximations, which imply CLTs and functional central limit theorems (FCLTs) in the sense of Donsker’s theorem. We also propose estimators for the asymptotic variance parameters and show their consistency, uniformly in the dimension. Lastly, Section 5 elaborates on several statistical problems to which our results are directly applicable. Proofs of the main results are provided in an appendix.

2. Projection-based analysis of high-dimensional time series

Let us assume that we observe $d$ possibly dependent time series such that at time $n$ we are given the observations

[TABLE]

where the dimension $d=d_{n}$ may grow with the sample size $n$, such that, as time proceeds, there may be more and more time series available. Equivalently, we are given a time series of length $n$ of possibly dependent random vectors

[TABLE]

of dimension $d_{n}$, constituting the $(n\times d_{n})$–dimensional data matrix $\mathcal{Y}_{n}=\left(Y_{i}^{(j)}\right)_{1\leq i\leq n,1\leq j\leq d_{n}}.$ We are interested in the second moment structure and thus assume $E(Y_{i}^{(j)})=0$ for all $j=1,\dots,d_{n}$, $i=1,\dots,n$ and $n\geq 1$. Our assumptions on the coordinate processes, basically that they are linear processes with sufficiently fast decreasing coefficients, are weak enough to cover the common framework of correlated ARMA($p,q$)-processes and also allow for a wide class of long-range dependent series.

Let us assume for a moment that ${\mathbf{Y}}_{n1},\dots,{\mathbf{Y}}_{nn}$ is stationary, a condition that we shall relax later, and let

[TABLE]

be a generic copy. For the analysis of such high-dimensional time series, the unknown variance-covariance matrix

[TABLE]

is of substantial interest, but difficult to estimate from past data, in particular if $d_{n}>>n$. It comprises the second-order information on the dependence structure of the $d_{n}$ variables. Any conclusions on the correlation structure have to rely on estimators calculated from the time series, and inferential procedures require appropriate large-sample asymptotics. Let

[TABLE]

be the $(d_{n}\times d_{n})$-dimensional sample variance–covariance matrix with elements

[TABLE]

Before proceeding, let us observe that the above framework also covers the case of a univariate time series $\{Z_{k}:k\geq 0\}$ as an interesting special case. The embedding is given by

[TABLE]

Then

[TABLE]

It follows that the $(h+1)$th element of the first row of $\widehat{\boldsymbol{\Sigma}}_{n}$ estimates $\gamma_{Z}(h)=E(Z_{0}Z_{h})$ using the observations $Z_{1},\dots,Z_{T}$, where $T=n+h$, and can be written as

[TABLE]

for $h=0,\dots,T-1$. In a similar way, one may consider autocovariances and cross-covariances of, say, $r_{n}$ time series $\{Z_{k}^{(l)}:k\geq 0\}$, $l=1,\dots,r_{n}$.

Estimators of $\boldsymbol{\Sigma}_{n}$ are also needed and have to be evaluated in terms of their asymptotic laws when interest focuses on the analysis of (a set of) linear combinations of ${\mathbf{Y}}_{n}$. Here ${\mathbf{Y}}_{n}$ may be a generic copy when the vector time series is strictly stationary or equal (in distribution) to ${\mathbf{Y}}_{nn}$ in the general case. Typical examples are convex combinations, contrasts and, more generally, projections. Thus let

[TABLE]

be a sequence of weights $w_{j}=w_{nj}$, not necessarily non-negative, with uniformly bounded $\ell_{1}$-norm, i.e.

[TABLE]

The class of weighting vectors of the form ${\mathbf{w}}_{n}=(w_{1},w_{2},\dots,w_{d_{n}})^{\prime}$ for some sequence $\{w_{i}:i\in\mathbb{N}\}$ with $\sum_{i=1}^{\infty}|w_{i}|<\infty$ certainly satisfies the assumptions and covers, for example, the case of averaging a finite number of coordinates. However, if $w_{i}>0$ holds for infinitely many $i\geq 1$, then $w_{i}=o(1)$, as $i\to\infty$. Hence infinitely many coordinates are not really taken into account. Allowing for weights that depend on the dimension $d_{n}$ substantially widens the scope. Now, for example, one may average all coordinates by using the weights $w_{ni}=1/d_{n}$, for $i=1,\dots,d_{n}$. Although the dependence on the sample size $n$ through the dimension $d_{n}$ may be of primary importance for high-dimensional problems, several of our results even allow the weights to depend on $n$.

The variance of the projection ${\mathbf{w}}_{n}^{\prime}{\mathbf{Y}}_{n}$ is given by ${\mathbf{w}}_{n}^{\prime}\boldsymbol{\Sigma}_{n}{\mathbf{w}}_{n}$ and can be estimated nonparametrically by the quadratic form $Q_{n}({\mathbf{w}}_{n})={\mathbf{w}}_{n}^{\prime}\widehat{\boldsymbol{\Sigma}}_{n}{\mathbf{w}}_{n},$ whose random fluctuations may severely affect any inferential procedure related to ${\mathbf{w}}_{n}^{\prime}{\mathbf{Y}}_{n}$. More generally, we shall study the bilinear form

[TABLE]

for weighting vectors ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ with uniformly bounded $\ell_{1}$-norms in the sense of (2.2), which corresponds to the estimator of the covariance of two projections ${\mathbf{v}}_{n}^{\prime}{\mathbf{Y}}_{n}$ and ${\mathbf{w}}_{n}^{\prime}{\mathbf{Y}}_{n}$. Observe that (2.3) also allows us to handle the case of weighted sums of subsets $\widehat{\boldsymbol{\Sigma}}_{\mathcal{I},\mathcal{J}}=\{\widehat{\sigma}_{ij}:i\in\mathcal{I},j\in\mathcal{J}\}$, where $\mathcal{I},\mathcal{J}\subset\{1,\dots,d_{n}\}$.

It is worth mentioning that $Q_{n}$ remains bounded even for degenerate covariance matrices, if ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ have uniformly bounded $\ell_{1}$-norm, as can be seen if we put $\boldsymbol{\Sigma}_{n}=\sigma\mathbf{1}\mathbf{1}^{\prime}$, where here and throughout the article $\mathbf{1}=(1,\dots,1)^{\prime}\in\mathbb{R}^{d_{n}}$, leading to

[TABLE]

Thus, the $\ell_{1}$-norm condition is a natural one: It ensures that $Q_{n}$ maps products $\mathcal{U}_{\delta}\times\mathcal{U}_{\delta}$ of $\delta$-balls

[TABLE]

for $\delta>0$, onto bounded sets, for all covariance matrices with uniformly bounded entries, thus including cases that correspond to perfectly correlated coordinates.

The behavior of projections for high-dimensional observations has also been studied by Diaconis and Freedman (1984), but from a different perspective. There it is shown that for large dimension $d$ projections ${\mathbf{w}}^{\prime}{\mathbf{X}}$ are asymptotically normal under weak assumptions by showing that, given a (non-random) sample $X_{1},\dots,X_{n}$ and a unit vector ${\mathbf{w}}$ uniformly distributed on the $d$-dimensional unit sphere, the empirical measure of the sample ${\mathbf{w}}^{\prime}{\mathbf{X}}_{1},\dots,{\mathbf{w}}^{\prime}{\mathbf{X}}_{n}$ converges weakly to a normal law, in probability. Further, the joint distribution of two linear combinations, say, ${\mathbf{w}}_{d}^{\prime}{\mathbf{Z}}$ and ${\mathbf{v}}_{d}^{\prime}{\mathbf{Z}}$, of a $d$-dimensional random vector ${\mathbf{Z}}$ possessing a Lebesgue density and being standardized, i.e. $E({\mathbf{Z}})=0$ and $E({\mathbf{Z}}{\mathbf{Z}}^{\prime})=\operatorname{id}_{d}$, is bivariate normal in that sense with unit variances and covariance ${\mathbf{w}}^{\prime}{\mathbf{v}}$. However, in that work the projection vectors are random and the data are assumed fixed (e.g. by conditioning) and constrained to satisfy conditions that are satisfied by, e.g., i.i.d. random vectors of dimension $d_{n}$ with i.i.d. entries. Contrary, in this paper the projection vectors are fixed and the observations random. Assuming a linear process framework, we provide Gaussian approximations for the sample estimates of the variances of the projections.

3. A framework for high-dimensional time series

Our results dealing with strong approximations of the bilinear forms introduced above rely on some linear process framework. This section provides a careful introduction, discusses some interesting properties and interpretations as well as introduces required notation and some preparatory approximations used later.

3.1. Model and assumptions

Let $\{\epsilon_{k}:k\in\mathbb{Z}\}$ be a sequence of independent random variables with mean zero, variances

[TABLE]

and uniformly bounded moments of the order $(4+\delta)$,

[TABLE]

for some $\delta>0$, such that $\gamma_{k}=E\epsilon_{k}^{4}$ and $\sigma_{k}^{2}$ are finite, for all $k\in\mathbb{Z}$.

We assume that the $\nu$th coordinate of ${\mathbf{Y}}_{n}$ is given by

[TABLE]

for coefficients $\{c_{nj}^{(\nu)}:j\in\mathbb{N}_{0}\}$, $\nu=1,\dots,d_{n}$. We mainly have in mind the case that we observe, at time $n$, the first $n$ observations of $d_{n}$ sequences $\{Y_{k}^{(\nu)}:k\geq 0\}$, $\nu=1,\dots,d_{n}$, but our results also allow for arrays $\{Y_{nk}^{(\nu)}:k\geq 0,n\geq 1\}$, since the coefficients may depend on $n$. Also notice that in model (3.2) $Y_{nk}^{(\nu)}$ is well defined for $k>n$.

Model (3.2) implies that the cross-sectional as well as serial correlations have a specific structure, since

[TABLE]

for $h>0$, $1\leq\nu,\mu\leq d_{n}$ and all $t$. Consequently, the cross-sectional variance-covariance matrix ${\mbox{Var\,}}({\mathbf{Y}}_{nt})$ is given by

[TABLE]

where ${\mathbf{C}}_{n}=(c_{nj}^{(\nu)})_{1\leq\nu\leq d_{n},1\leq j}$ is the $(d_{n}\times\infty)$-dimensional matrix with column vectors ${\mathbf{c}}_{nj}=(c_{nj}^{(1)},\dots,c_{nj}^{(d_{n})})^{\prime}$, $j\geq 0.$ The lag $h$ serial covariance matrix attains the representation

[TABLE]

where $L$ denotes the lag operator that acts on all columns, i.e. $L^{-1}{\mathbf{C}}_{n}=(c_{n,i+1}^{(\nu)})_{i\geq 0,1\leq\nu\leq d_{n}}$.

We shall impose the following condition on the decay of the coefficients $c_{nj}^{(\nu)}$, which is similar to the assumption imposed in Johnstone and Lu (2009), where it controls the principal component eigenvectors within a factor model (also see our discussion below), and to condition (2.4) in Chan et al. (2013), where it controls the error terms of a panel time series model.

Assumption (A) The sequences $\{c_{nj}^{(\nu)}:j\in\mathbb{N}_{0}\}$ satisfy

[TABLE]

for some $0<\theta<1/2$.

Here and in the sequel $a_{n}<<b_{n}$ stands for $a_{n}=O(b_{n})$. Further, we shall write $a_{nm}\stackrel{{\scriptstyle n,m}}{{<<}}b_{nm}$ if there exists a constant $C$ such that $a_{nm}\leq Cb_{nm}$ for all $n,m$

Indeed, (3.6) covers not only short memory processes for which the covariances, say, $r_{k}=E(X_{0}X_{k})$, are summable, i.e. $\sum_{k}|r_{k}|<\infty$, but also many long-range dependent series. An example for the latter is fractionally integrated noise of order $d\in(-1/2,1/4-\theta/2)$, i.e. a stationary solution of the equation

[TABLE]

where $\{\epsilon_{t}\}$ is a white-noise series, that is given by $X_{t}=\sum_{k=0}^{\infty}\theta_{k}\epsilon_{t-k}$ with coefficients $\theta_{k}=\Gamma(k+d)/(\Gamma(k+1)\Gamma(d))\sim k^{d-1}/\Gamma(d)$, see e.g. Steland (2012). The growth condition $O(j^{-3/4-\theta})$ for some $\theta>0$ on the coefficients of the linear processes arises also in other works, especially in Wu et al. (2010) where the asymptotics of sample autocovariances of a linear process is studied. The case that the coefficients are $O(j^{-3/4}L(j))$ for some slowly-varying function $L$ represents a boundary case. Here one can obtain Gaussian limits for sample autocovariances, see (Wu et al., 2010, Theorem 3) and Hosking (1996) for i.i.d. Gaussian innovations, but then the convergence rate changes from $n^{-1/2}$ to $(n/\bar{L}_{n})^{-1/2}$ where $\bar{L}_{n}=\sum_{i=1}^{n}L^{4}(i)/i$.

It is worth noting that, under certain circumstances, one may interpret model (3.2) as an infinite-dimensional latent one-factor model. For that purpose, let us assume for a moment that ${\mathbf{Y}}_{1},\dots,{\mathbf{Y}}_{n}$ is a stationary series, such that $\sigma_{k}^{2}=\sigma^{2}$ for all $k$. Then, by stationarity, ${\mathbf{Y}}_{n}=(Y^{(1)},\dots,Y^{(d_{n})})^{\prime}$ satisfies

[TABLE]

where $\{\epsilon_{n}\}$ is the unobservable common factor.

3.2. Some preparatory approximations

In the main proof we shall study in detail the linear process $\sum_{j=0}^{\infty}c_{j}^{w}\epsilon_{k-j}$, $k\geq 1$, with coefficients

[TABLE]

associated to a weighting vector ${\mathbf{w}}_{n}$. Behind our main results are martingale approximations related to that linear process, whose definitions require the following quantities, which are controlled under Assumption (A) by virtue of Lemma 3.1 given below. Let

[TABLE]

[TABLE]

and

[TABLE]

Clearly, the quantities $\widetilde{f}_{l,i}^{(n)}$ have the scaling property

[TABLE]

for arbitrary $s_{1},s_{2}\in\mathbb{R}$.

In addition, under Assumption (A) we have for $\eta>0$

[TABLE]

for constants $C_{1},C_{2},C_{3}$ not depending on $n$ and the weighting vectors.

Lemma 3.1*.*

Suppose that ${\mathbf{v}}_{n},{\mathbf{w}}_{n}$ have uniformly bounded $\ell_{1}$-norm in the sense of equation (2.2). Then Assumption (A) implies

[TABLE]

[TABLE]

[TABLE]

where the constant may differ from line to line and depends on the weighting vectors only through their $\ell_{1}$–norms. There exist

[TABLE]

such that

[TABLE]

for all $n^{\prime},m^{\prime}=0,1,\cdots$.

Further, if ${\mathbf{v}}_{n},{\mathbf{w}}_{n},\widetilde{{\mathbf{v}}}_{n},\widetilde{{\mathbf{w}}}_{n}$, $n\geq 1$, have uniformly bounded $\ell_{1}$-norms, then there exist

[TABLE]

with

[TABLE]

Remark 3.1*.*

If the moments up to the order $4$ are stationary such that $\gamma_{k}=\gamma$ and $\sigma_{k}^{2}=\sigma^{2}$, say, (3.22) is a consequence of (3.19) and

[TABLE]

as well as

[TABLE]

Note also that, in general, (3.18)–(3.20) and (3.22) ensure that the quantities $\alpha_{n}^{2}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ and $\beta_{n}({\mathbf{v}}_{n},{\mathbf{w}}_{n},\widetilde{{\mathbf{v}}}_{n},\widetilde{{\mathbf{w}}}_{n})$ are bounded.

Suppose that the coefficients of the time series as well as the weighting vectors do not depend on $n$, i.e. $c_{nj}^{(\nu)}=c_{j}^{(\nu)}$, $\nu=1,2,\dots,d_{n}$, $j\geq 0$, and $w_{n\nu}=w_{\nu}$, $v_{n\nu}=v_{\nu}$, $\nu=1,2,\dots$. Since then $\widetilde{f}_{0,0}^{(n)}=\sum_{j=0}^{\infty}\left(\sum_{\nu=1}^{d_{n}}[c_{j}^{(\nu)}]^{2}\right)^{2}$ and

[TABLE]

we have $\alpha_{n}({\mathbf{v}}_{n},{\mathbf{w}}_{n})\to\alpha^{*}$, as $n\to\infty$. The limit $\alpha^{*}$ is obtained by replacing formally $d_{n}$ by $\infty$. Similarly, one obtains the convergence $\beta_{n}({\mathbf{v}}_{n},{\mathbf{w}}_{n},\widetilde{{\mathbf{v}}}_{n},\widetilde{{\mathbf{w}}}_{n})\to\beta^{*}$, as $n\to\infty$, and an explicit formula for $\beta^{*}$.

4. Asymptotics

For many sequences, $X_{n}$, $n\geq 1$, of random variables, such as sums of i.i.d. random variables with finite second moment, a strong approximation with rate holds true, i.e. after redefining the process on a new probability space, there exists a Brownian motion $B(t)$, $t\geq 0$, such that, on that new space,

[TABLE]

a.s., as $n\to\infty$, for some positive constant $\sigma$ and a constant $\lambda>0$. Results of this type date back to the seminal work of Komlós et al. (1975), Komlós et al. (1976) and have been extended and refined since then, in particular to dependent sequences attaining values in a Hilbert space and martingales, see e.g. Philipp (1986). Such a strong approximation result also yields an approximation of the rescaled càdlàg process $n^{-1/2}X_{{\lfloor tn\rfloor}}$, $t\in[0,1]$, by the Brownian motion $\sigma B(t)$, $t\in[0,1]$, and implies the FCLT, i.e. the weak convergence

[TABLE]

as $n\to\infty$, where $\Rightarrow$ signifies weak convergence in the Skorohod space $D[0,1]$. This, in turn, implies the weak convergence of continuous mappings of $n^{-1/2}X_{n}({\lfloor n\bullet\rfloor})$. Observe that when $X_{n}=\sum_{i=1}^{n}\xi_{i}$ is a sum of i.i.d. random variables with $E(\xi_{1}^{2})<\infty$, then we obtain the classical Donsker theorem and, for $t=1$, the CLT.

In order to obtain strong approximations for the bilinear form $Q({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, we shall derive a martingale approximation for a partial sum associated to $Q({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, which eases the verification of general conditions for sequences taking values in Hilbert spaces due to Philipp (1986) to obtain a strong approximation.

We need to introduce further notation. Define

[TABLE]

for $n,k\geq 1$. To be precise, our results shall deal with

[TABLE]

and the associated càdlàg processes

[TABLE]

and

[TABLE]

Observe that $\mathcal{D}_{n}^{0}(t)$ does not depend on the true variance-covariance matrices $\{\boldsymbol{\Sigma}_{nk}:1\leq k\leq n\}$. Processes of this form are therefore frequently used in change-point analysis, see, e.g., Csörgő and Horváth (1997) and Steland (2015). If the dependence of the above quantities on ${\mathbf{v}}_{n},{\mathbf{w}}_{n}$ matters, we shall indicate this in our notation and then write

[TABLE]

Recalling that $\widehat{\boldsymbol{\Sigma}}_{n}=n^{-1}\widehat{\boldsymbol{\Sigma}}_{n,n}$, cf. (2.1) and (4.1),

[TABLE]

is the centered and scaled version of the bilinear form $Q({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, where

[TABLE]

If $\{{\mathbf{Y}}_{ni}:1\leq i\leq n\}$ is stationary, then $\boldsymbol{\Sigma}_{n}$ simplifies to $\boldsymbol{\Sigma}_{n}=E({\mathbf{Y}}_{n1}{\mathbf{Y}}_{n1}^{\prime})$.

For two weighting vectors ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ we may associate the martingales

[TABLE]

and the corresponding càdlàg processes

[TABLE]

It turns out that, under the assumptions of the article, those martingales are close to $D_{nk}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ and $\mathcal{D}_{n}(t;{\mathbf{v}}_{n},{\mathbf{w}}_{n})$, respectively, which is the key to obtain large-sample asymptotics in terms of strong approximations to infer second order information, i.e. variances and covariances, of projections.

The following theorem shows that bilinear forms of uniformly bounded $\ell_{1}$-projections satisfy a strong approximation result that implies the functional central limit theorem and therefore also the central limit theorem. Recall that to a standard Brownian motion $B_{n}(t)$, $t\in[0,\infty)$, we may define the rescaled version $\overline{B}_{n}(s)=n^{-1/2}B_{n}(sn)$, $s\in[0,1]$, which satisfies $E(\overline{B}_{n}(s)\overline{B}_{n}(t))=\min(s,t)$ for $s,t\in[0,1]$. For convenience, we shall denote $\overline{B}_{n}$ also by $B_{n}$ and call it the $[0,1]$–version of $B_{n}$.

Theorem 4.1*.*

Suppose ${\mathbf{Y}}_{ni}$, $1\leq i\leq n$, $n\geq 1$, is a vector time series according to model (3.2) that satisfies Assumption (A). Let ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ be weighting vectors with uniformly bounded $\ell_{1}$–norm in the sense of (2.2). Then, for each $n\in\mathbb{N}$, there exists an equivalent version of $D_{nk}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ and thus of $D_{n}(t;{\mathbf{v}}_{n},{\mathbf{w}}_{n})$, $t\geq 0$, again denoted by $D_{nk}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ and $D_{n}(t;{\mathbf{v}}_{n},{\mathbf{w}}_{n})$, and a standard Brownian motion $\{B_{n}(t):t\geq 0\}$, which depends on $({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, i.e. $B_{n}(t)=B_{n}(t;{\mathbf{v}}_{n},{\mathbf{w}}_{n})$, both defined on some probability space $(\Omega_{n},\mathcal{F}_{n},P_{n})$, such that for some $\lambda>0$ and a constant $C_{n}$

[TABLE]

and

[TABLE]

for $k\leq n$. If

[TABLE]

as $n\to\infty$, this implies the strong approximation

[TABLE]

as $n\to\infty$, for the $[0,1]$–version of $B_{n}$. Further,

[TABLE]

and

[TABLE]

a.s., as $n\to\infty$, where $B_{n}^{0}(t)=B_{n}(t)-tB_{n}(1)$, $t\in[0,1]$, denotes the Brownian bridge associated to the $[0,1]$–version of $B_{n}$.

Clearly, (4.10) implies

[TABLE]

as $n\to\infty$, i.e. $\mathcal{D}_{n}(1;{\mathbf{v}}_{n},{\mathbf{w}}_{n})$ is asymptotically $\mathcal{N}(0,\alpha_{n}^{2}({\mathbf{v}}_{n},{\mathbf{w}}_{n}))$, suggesting the standardized statistic

[TABLE]

which is asymptotically standard normal by (4.13), to draw statistical inference on the covariance of two projections ${\mathbf{v}}_{n}^{\prime}{\mathbf{Y}}_{n}$ and ${\mathbf{w}}_{n}^{\prime}{\mathbf{Y}}_{n}$. By virtue of the scaling property (3.13), which carries over to $\alpha_{n}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, $\mathcal{D}_{n}^{*}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ depends on the weighting vectors ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ only through the associated unit vectors ${\mathbf{v}}_{n}^{*}={\mathbf{v}}_{n}/\|{\mathbf{v}}_{n}\|$ and ${\mathbf{w}}_{n}^{*}={\mathbf{w}}_{n}/\|{\mathbf{w}}_{n}\|$.

For ${\mathbf{v}}_{n}={\mathbf{w}}_{n}$ we obtain the asymptotic normality of the statistic

[TABLE]

which allows to draw inference on the variance of the projection $\pi_{n}={\mathbf{v}}_{n}^{\prime}{\mathbf{Y}}_{n}$ of ${\mathbf{Y}}_{n}$ onto $\operatorname{span}\{{\mathbf{v}}_{n}\}$. Estimation of the asymptotic variance parameter is discussed below.

Clearly, it is of interest to project onto, say $K$, vectors ${\mathbf{v}}_{n1},\dots,{\mathbf{v}}_{nK}$, as discussed in greater detail in Subsection 5.2. The following generalization to a finite number of bilinear forms provides the required multivariate extension.

Recall that ${\mathbf{B}}(t)$, $t\geq 0$, is a Brownian motion in $\mathbb{R}^{K}$ with covariance matrix ${\mathbf{V}}>0$, if ${\mathbf{B}}(t)={\mathbf{V}}^{1/2}{\mathbf{B}}_{0}(t)$ for a standard Brownian motion, i.e. ${\mathbf{B}}_{0}(s)\sim\mathcal{N}({0},sI_{K})$ and $E[({\mathbf{B}}_{0}(s)-{\mathbf{B}}_{0}(t))({\mathbf{B}}_{0}(s)-{\mathbf{B}}_{0}(t))^{\prime}]=(s-t)I_{K}$ for $0\leq s\leq t$. This means, ${\mathbf{B}}(t)$ is a Gaussian process with independent increments, $E({\mathbf{B}}(s){\mathbf{B}}(t)^{\prime})=\min(s,t){\mathbf{V}}$, $s,t\geq 0$, and thus ${\mathbf{B}}(t)$, $t\geq 0$, is characterized by ${\mathbf{V}}=E({\mathbf{B}}(1){\mathbf{B}}(1)^{\prime})$. The latter still makes sense if ${\mathbf{V}}$ is singular.

Theorem 4.2*.*

Let $\{{\mathbf{v}}_{nj},{\mathbf{w}}_{nj}:1\leq j\leq K\}$ be weighting vectors of dimension $d_{n}$ satisfying condition (2.2). Then, under the assumptions of Theorem 4.1, there exists a $K$–dimensional Brownian motion $\{{\mathbf{B}}^{(n)}(t):t\in[0,1]\}$ with coordinates $B_{ni}=B_{n}(t;{\mathbf{v}}_{ni},{\mathbf{w}}_{ni})$, $t\in[0,1]$, $i=1,\dots,K$, characterized by

[TABLE]

for $1\leq i,j\leq K$ with $i\not=j$ and $E(B_{n}^{2}(1;{\mathbf{v}}_{ni},{\mathbf{w}}_{ni}))=\alpha_{n}({\mathbf{v}}_{ni},{\mathbf{w}}_{ni}),$ for $i=1,\dots,K$, such that

[TABLE]

a.s., as $n\to\infty$, where $\|\bullet\|$ denotes an arbitrary vector norm on $\mathbb{R}^{K}$.

Having in mind the application of the above result to change-point detection, see Subsection 5.4, the following corollary dealing with the frequently used maximally selected CUSUM statistic is of interest.

Corollary 4.1*.*

Suppose that ${\mathbf{Y}}_{n1},\dots,{\mathbf{Y}}_{nn}$ is a $d_{n}$–dimensional vector time series following model (3.2) which satisfies Assumption (A). Let ${\mathbf{v}}_{n}$ and ${\mathbf{w}}_{n}$ be weighting vectors with uniformly bounded $\ell_{1}$–norm in the sense of (2.2). Then, after redefining the series on a new probability space, there exist standard Brownian motions $B_{n}$ on $[0,1]$, such that

[TABLE]

and

[TABLE]

for the Brownian bridge $B_{n}^{0}$ associated to $B_{n}$, a.s., as $n\to\infty$.

We conjecture that the strong approximation (4.7) holds with the rate $O(n^{-\lambda})$, i.e. one can select a universal constant $C$ in (4.7), such that in turn the strong approximations (4.10)–(4.16) would hold with the rate $O(n^{-\lambda})$ as well. This is also plausible from the following more general result which provides us with a strong approximation with that rate when weighting over the sample sizes with $\ell_{1}$-bounded weights: Attach to each sample size $n$ a weight $\lambda_{n}$ such that $\sum_{n=1}^{\infty}|\lambda_{n}|<\infty$. Define

[TABLE]

where

[TABLE]

for $n,m\geq 1$, based on the full set of observations $\{Y_{i}^{(\nu)}:1\leq\nu<d_{n},1\leq i\leq n,n\geq 1\}$. At this point it is worth recalling that in model (3.2) the $d_{n}$ coordinate processes of $Y_{ni}$ are well defined for all $i\geq 1$. For fixed $k$, $D_{k}(\{{\mathbf{v}}_{n},{\mathbf{w}}_{n}\})$ depends on the array $\{Y_{ni}^{(\nu)},i\leq k,1\leq\nu\leq d_{n},n\in\mathbb{N}\}$. Notice that $D_{nk}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, defined in (4.3), appears as a special case of (4.17). We may define the following càdlàg process associated to (4.17)

[TABLE]

for some sample size $N$ (which can be equal to $n$).

We have the following general large sample approximation result.

Theorem 4.3*.*

Let $\{{\mathbf{v}}_{n}\}$ and $\{{\mathbf{w}}_{n}\}$ be two sequences of weighting vectors with uniformly bounded $\ell_{1}$–norm in the sense of (2.2). Assume that $\{Y_{i}^{(\nu)}:1\leq\nu<\infty,1\leq i\leq n,n\geq 1\}$ follows model (3.2) and satisfies Assumption (A). Let $\{\lambda_{n}\}$ be weights with $\sum_{n=1}^{\infty}|\lambda_{n}|<\infty$. Then there exist constants $\alpha(\{{\mathbf{v}}_{n},{\mathbf{w}}_{n}\})$, $\lambda>0$ and $C$ (not depending on the sample size), such that for equivalent versions and a standard Brownian motion $B$ on $[0,\infty)$, defined on a new probability space,

[TABLE]

a.s., for all $t>0$, leading to the strong invariance principle with rate,

[TABLE]

a.s., for the $[0,1]$–version of $B$, where $C$ is the same constant as in (4.20) and may depend on the sequence $\{\lambda_{n}\}$.

It remains to discuss how one may estimate the asymptotic variance parameters $\alpha^{2}=\alpha^{2}({\mathbf{v}},{\mathbf{w}})$ and $\beta^{2}=\beta^{2}({\mathbf{v}}_{r},{\mathbf{w}}_{r},{\mathbf{v}}_{s},{\mathbf{w}}_{s})$ arising in the strong approximations. For a sequence of lag truncation constants $m=m_{n}$, $n\geq 1$, and weights $\{w_{mh}:h\in\mathbb{Z},m\in\mathbb{N}\}$ define the estimator

[TABLE]

where

[TABLE]

with $\widehat{\mu}_{n}^{(v,w)}=n^{-1}\sum_{j=1}^{n}Y_{j}^{(v)}Y_{j}^{(w)}$ and $Y_{i}^{(v)}={\mathbf{v}}^{\prime}{\mathbf{Y}}_{i}=\sum_{j=0}^{\infty}c_{j}^{(v)}\epsilon_{i-j}$, $i\geq 1$, with weights $c_{j}^{(v)}=\left(\sum_{\nu=1}^{d}c_{j}^{(\nu)}v_{j}\right)$, $j\geq 0$. At this point we consider the dimension as a further parameter, such that $Y_{i}^{(v)}$, $\widehat{\Gamma}_{h}=\widehat{\Gamma}_{h}(d)$ and $\alpha^{2}=\alpha^{2}(d)$ depend on the dimension $d$.

More generally, the covariance parameters $\beta_{rs}=\beta({\mathbf{v}}_{r},{\mathbf{w}}_{r};{\mathbf{v}}_{s},{\mathbf{w}}_{s})$, $1\leq r,s\leq K$, arising in Theorem 4.2 can be estimated by

[TABLE]

where $m=m_{n}$ is a lag truncation sequence and

[TABLE]

with $\widehat{\mu}_{n}(r)=n^{-1}\sum_{j=1}^{n}Y_{j}({\mathbf{v}}_{r})Y_{j}({\mathbf{w}}_{r})$, for $1\leq r,s\leq K$.

For the weights we impose the following standard assumptions.

• (W1)

$w_{mh}\to 1$, as $m\to\infty$, for all $h\in\mathbb{Z}$.

• (W2)

$0\leq w_{mh}\leq W<\infty$ for some constant $W$, for all $m\geq 1$, $h\in\mathbb{Z}$.

Typically, the weights are defined by a kernel function $w$ via $w_{mh}=w(h/b_{m})$ for some bandwidth parameter $b=b_{m}$. Examples are the truncated kernel, $k_{tr}(x)=\mathbf{1}(|x|\leq 1)$ with $b=m+1$, or Bartlett’s estimator given by the triangular weight function $w(x)=(1-x)\mathbf{1}(|x|\leq 1)$, see Newey and West (1987) and Andrews (1991) amongst others.

Jirak (2012) has studied estimation of long run variance parameters associated to $d\to\infty$ nonlinear time series, based on results of Wu (2009). Those results are not applicable to our general setting, since the coefficients of the linear processes $Y_{i}^{(v)}$ depend on $n$ if $d=d_{n}$. The following theorem provides the $L_{1}$–consistency, uniformly over $d$ and a possibly infinite number of weighting vectors, as long as they are taken from an $\ell_{1}$–bounded set. This yields the consistency of $\widehat{\alpha}_{n}^{2}(d_{n})$ for a growing dimension $d_{n}\to\infty$, since

[TABLE]

without a constraint on the growth of $d_{n}$, thus going beyond the known results.

Theorem 4.4*.*

Assume (W1) and (W2) and suppose that $c_{nj}^{(\nu)}=c_{j}^{(\nu)}$ for $\nu=1,2,\dots$, $j\geq 0,n\geq 1,$ satisfy the decay condition

[TABLE]

for some $\delta>0$. Suppose that $\epsilon_{k}$ are i.i.d. with $\max_{k}E|\epsilon_{k}|^{8}<\infty$ and assume that $m=m_{n}\to\infty$ with $m^{2}/n=o(1)$, as $n\to\infty$. If ${\mathbf{v}},{\mathbf{w}}$ are weighting vectors with $\|{\mathbf{v}}\|_{\ell_{1}},\|{\mathbf{w}}\|_{\ell_{1}}<\infty$, then

[TABLE]

Further, if ${\mathbf{v}}_{\ell},{\mathbf{w}}_{\ell}$, $\ell\geq 1$, are weighting vectors with uniformly bounded $\ell_{1}$–norm, i.e.

[TABLE]

for some constant $C_{v,w}$, then

[TABLE]

as $n\to\infty$.

Suppose that the dimension is a random variable, $D$, and drawn according to some (prior) probability measure $Q$ on $\mathbb{N}$, such that $D$ and $\{{\mathbf{Y}}_{ni}\}$ are independent, and statistical inference is conducted given $D=d$. A natural measure to evaluate the estimator $\widehat{\alpha}_{n}^{2}(D)$ in this setting is the expected mean deviation,

[TABLE]

Theorem 4.4 readily implies $EMD(\widehat{\alpha}_{n}^{2}(D))\to 0,$ as $n\to\infty$.

5. Applications

The results of the present paper have direct applications to several problems and procedures, respectively, which are extensively studied for high–dimensional time series, especially for big data. They contribute novel large-sample approximations for making inference based on the corresponding statistics, usually projections.

5.1. Optimal portfolio selection

The problem of optimal portfolio selection, dating back to Markowitz’ seminal work, see Markowitz (1952), is an intrinsically high-dimensional problem. We are given a usually large number $d_{n}$ of assets and associated returns ${\mathbf{R}}_{n}=(R_{n}^{(1)},\dots,R_{n}^{(d_{n})})^{\prime}$ corresponding to the time period $[n-1,n]$ with mean vector $\mu$ and covariance matrix $\boldsymbol{\Sigma}=(\sigma_{ij})_{ij}$. Let us assume that the asset return vector time series satisfies the standing assumptions of the paper. Since $\sigma_{ij}$ is the covariance between the return of asset $i$ and asset $j$, $1\leq i,j,\leq d_{n}$, it is not restrictive, by the very nature of the problem, to assume that the entries of $\boldsymbol{\Sigma}$ neither depend on $n$ nor $d_{n}$. Suppose that an investor holds at time $n-1$ the position $w_{nj}$ in asset $j$, where $w_{nj}>0$ represents a long position and $w_{nj}<0$ a short position. W.l.o.g. we may assume that the initial value (capital) at time $n-1$ equals $V=\sum_{j=1}^{d_{n}}w_{nj}=1$, such that the value at time instant $n$ is ${\mathbf{w}}_{n}^{\prime}{\mathbf{R}}_{n}$. A classical formulation of the portfolio optimization problem is to minimize the risk, defined as the variance, associated to the portfolio return ${\mathbf{w}}_{n}^{\prime}{\mathbf{R}}_{n}$ at time $n$, i.e. to consider the problem

[TABLE]

whose solution is known to be ${\mathbf{w}}_{n}=(\mathbf{1}^{\prime}\Sigma^{-1}\mathbf{1})^{-1}\mathbf{1}^{\prime}\Sigma^{-1}.$ Here $\mathbf{1}$ is the $n$-vector with unit entries. Obviously, if that optimal solution satisfies the no-short-sales condition ${\mathbf{w}}_{n}\geq 0$, then $\|{\mathbf{w}}_{n}\|_{\ell_{1}}={\mathbf{w}}_{n}^{\prime}\mathbf{1}=(\mathbf{1}^{\prime}\Sigma^{-1}\mathbf{1})^{-1}\mathbf{1}^{\prime}\Sigma^{-1}\mathbf{1}=1,$ such that the optimal portfolio has uniformly bounded $\ell_{1}$-norm. The mean-variance formulation adds a constraint on the target mean portfolio return and thus considers the problem

[TABLE]

for some $\mu_{0}$. The solution is

[TABLE]

with $a=\mathbf{1}^{\prime}\boldsymbol{\Sigma}^{-1}\mathbf{1}$, $b=\mathbf{1}^{\prime}\boldsymbol{\Sigma}^{-1}\mu$ and $c=\mu^{\prime}\boldsymbol{\Sigma}^{-1}\mu$. Based on estimates of $\boldsymbol{\Sigma}^{-1}$ and $\mu$ from past data, one calculates the optimal portfolio which is then held until the next rebalancing. If the dimension is large, inverting the sample covariance matrix may result in substantial numerical instability and becomes impossible if the dimension is larger than the sample size. Shrinking is a commonly applied approach for regularization, in order to obtain a stable and invertible estimator, see Ledoit and Wolf (2003), amongst others. (We also refer to our subsection 5.3 for more on this.) It is worth mentioning that adding a non-negativity constraint (i.e. no-short-sales) has been observed to have a similar regularizing effect, see Jagannathan and Ma (2003).

To obtain sparse portfolios, Brodie et al. (2009) proposed to add an $\ell_{1}$-constraint. To discuss their approach, let us assume that we are given an additional independent learning sample $(r_{t}^{(j)})_{t,j}$, $1\leq j\leq d_{N}$, $1\leq t\leq N$, of size $N$ of returns for the same assets, such that in particular $d_{N}=d_{n}$, which is used to estimate the optimal weighting vector. Markowitz’ problem is equivalent to

[TABLE]

which suggests the following empirical version

[TABLE]

where $\mathcal{X}_{N}=(r_{t}^{(j)})_{t,j}$ is the $N\times d_{n}$ dimensional data matrix of returns with rows ${\mathbf{r}}_{t}^{\prime}=(r_{t}^{(1)},\dots,r_{t}^{(d_{n})})$, $t=1,\dots,N$, and $\widehat{\mu}_{N}=N^{-1}\sum_{t=1}^{N}{\mathbf{r}}_{t}$, see (Brodie et al., 2009, formula 1). These authors examine the $\ell_{1}$-regularized version of (5.1),

[TABLE]

for some regularization parameter $\rho>0$. Whereas for large values of $\rho$ the classical solution is recovered, smaller values lead to effective $\ell_{1}$-penalization and sparse portfolio vectors with only relatively few active positions.

5.2. Projections onto lower-dimensional subspaces, sparse principal components and the LASSO

A primary goal of multivariate statistical analysis and an indispensable tool to investigate big data is to project high–dimensional data onto lower-dimensional subspaces. Therefore the results of this paper directly address various approaches that are based on $\ell_{1}$ projection vectors to ensure sparse representations of the high–dimensional data. In particular, the results apply to the recently proposed methods of Jolliffe et al. (2003), Shen and Huang (2008), Witten et al. (2009) for sparse principal component analysis and the LASSO, see Tibshirani (1996).

Recall that for $L$ vectors ${\mathbf{w}}_{n}^{(k)}$, $k=1,\dots,L$, of dimension $d_{n}$ we may define the linear mapping $\pi_{n}:\mathbb{R}^{d_{n}}\to\mathbb{R}^{d_{n}}$,

[TABLE]

onto the associated linear subspace $\text{span}\{{\mathbf{w}}_{n}^{(k)}:k=1,\dots L\}$, which is the orthogonal projection if the ${\mathbf{w}}_{n}^{(k)}$ are orthonormal. Otherwise, the latter is given by ${\mathbf{P}}_{n}({\mathbf{P}}_{n}^{\prime}{\mathbf{P}}_{n})^{-}{\mathbf{P}}_{n}^{\prime}$, where ${\mathbf{A}}^{-}$ denotes a generalized inverse of a square matrix ${\mathbf{A}}$. In both cases, inferential procedures for $\pi_{n}$ can be based on the asymptotics of the dimension-reducing statistic ${\mathbf{P}}_{n}^{\prime}{\mathbf{Y}}_{n}$, such that our results apply, provided the projection vectors in use and the vector time series satisfy our assumptions. However, in general, the columns of ${\mathbf{P}}_{n}^{\prime}$ and ${\mathbf{P}}_{n}({\mathbf{P}}_{n}^{\prime}{\mathbf{P}}_{n})^{-}{\mathbf{P}}_{n}^{\prime}$, respectively, do not necessarily satisfy the uniform $\ell_{1}$–condition.

One may apply the following ad hoc approach to transform each vector ${\mathbf{w}}_{n}^{(j)}$ such that it satisfies an $\ell_{1}$-constraint $\|{\mathbf{w}}_{n}^{(j)}\|_{\ell_{1}}\leq c$, $j=1,\dots,L,$ for some preassigned constant $c$, by calculating the $\ell_{1}$-constrained optimal projection of ${\mathbf{w}}_{n}^{(j)}$ onto $\text{span}\{{\mathbf{w}}_{n}^{(j)}\}$,

[TABLE]

whose solution is known to be

[TABLE]

Here $S(a,\delta)=\text{sgn}(a)(|a|-\delta)_{+}$ is the soft-thresholding function, $x_{+}=x$ if $x>0$ and $=0$ otherwise, and $\delta\geq 0$ is chosen such that $\|S({\mathbf{w}}_{n}^{(j)},\delta)\|_{\ell_{1}}=c$, see (Witten et al., 2009, Lemma 2.2) and Shen and Huang (2008).

Imposing $\ell_{1}$-constraints is also the basic idea behind most approaches to define a sparse principal component analysis for high-dimensional data, which aims at determining lower-dimensional subspaces generated by $\ell_{1}$-vectors in such a way that they explain a large part of the variation in the observed data and also provide low-rank approximation of the data matrix. Let $\mathcal{X}_{n}$ be a $n\times d_{n}$-dimensional data matrix with centered columns corresponding to an independent sample of size $n$ of the variables $Y^{(1)},\dots,Y^{(d_{n})}$, whose columns are assumed to be centered. The simplified component technique-lasso (SCoTLASS) approach of Jolliffe et al. (2003) defines the first sparse principal component as a solution of the optimization problem

[TABLE]

Further sparse components are obtained by maximizing the same objective function under above constraints and the additional constraints that the further component is orthogonal to the previous components. In this way, after $L$ steps we obtain $L$ orthogonal $\ell_{1}$-vectors.

In a similar way, Witten et al. (2009) propose a sparse principal component analysis by solving, for the first component, the penalized matrix decomposition problem (PMD) with $\ell_{1}$-constraints,

[TABLE]

Observing that for fixed ${\mathbf{v}}$ the solution is given by ${\mathbf{u}}=\mathcal{X}_{n}{\mathbf{v}}/\|\mathcal{X}_{n}{\mathbf{v}}\|_{\ell_{2}}$, the first sparse principal component of PMD with $\ell_{1}$-constraints also solves SCoTLASS, see (Witten et al., 2009, p. 525). However, the PMD approach does not constrain the further components to be orthogonal, so that it differs from SCoTLASS. Closely related is the sparse PCA (SPCA) method of Shen and Huang (2008). They consider the problem to determine a regularized low-rank matrix approximation,

[TABLE]

for several penalty terms $p_{\rho}({\mathbf{v}})$ including the case $\rho\|{\mathbf{v}}\|_{\ell_{1}}$, for some $\rho>0$.

The LASSO, see Tibshirani (1996) and Tibshirani (2011), is a well established approach to determine $\ell_{1}$–sparse coefficient vectors in a high–dimensional linear regression model

[TABLE]

where the conditional expectation $E(Y_{t}|{\mathbf{X}}_{t})={\mathbf{X}}_{t}^{\prime}\beta_{0}$ is the $L_{2}$–optimal predictor for $Y_{t}$ given ${\mathbf{X}}_{t}$. Given some estimator $\widehat{\beta}_{n}$ of the unknown coefficient vector $\beta_{0}\in R^{d}$, the linear projection $\pi_{n}({\mathbf{X}})={\mathbf{X}}^{\prime}\widehat{\beta}_{n}$ is used to predict the outcome of the response for some (future) observed ${\mathbf{X}}$. The LASSO minimizes the $\ell_{1}$–constrained least squares criterion,

[TABLE]

for some bound $c>0$ for the $\ell_{1}$–norm, such that the resulting estimator $\widehat{\beta}_{n}$ is $\ell_{1}$–sparse. Consequently, our results can be applied to draw inference on the variance of the LASSO-based prediction $\pi_{n}({\mathbf{X}})={\mathbf{X}}^{\prime}\widehat{\beta}_{n}$, provided the regressor vector time series satisfies the assumptions of this article. Among the diverse applications where the prediction of the response in the presence of a large number of correlated explanatory variables is of interest, is the analysis of genetic association studies, see Allison et al. (2006) for its basic analysis and Biedermann et al. (2006) for nonparametric tests and further discussion. In such studies the regressors are gene expression data and the response is a phenotype. A sparse coefficient vector with only a few nonvanishing entries may allow to identify (groups of) genes which are associated to the phenotype.

5.3. Shrinkage estimation

In many applications, from a statistical point of view, when estimating the common variance-covariance matrix $\boldsymbol{\Sigma}_{n}$ of a stationary vector time series ${\mathbf{Y}}_{n1},\dots,{\mathbf{Y}}_{nn}$ of dimension $d_{n}$, one has an interest to regularise $\widehat{\boldsymbol{\Sigma}}_{n}$ to improve its (finite-sample) properties such as its mean-squared error $E[\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}_{n}\|^{2}_{F}]$ or its condition number, defined to be the ratio of its largest to its smallest eigenvalue. This is of particular interest if one needs an invertible estimator of $\boldsymbol{\Sigma}_{n}$. One well-established possibility to regularise $\widehat{\boldsymbol{\Sigma}}_{n}$ (Ledoit and Wolf (2004), Sancetta (2008)) is to consider a shrinkage estimator defined by a linear (in fact a convex) combination of $\widehat{\boldsymbol{\Sigma}}_{n}$ with a well-conditioned ”target”. Already in the population, for situations where the dimensionality $d_{n}$ is in the order magnitude of the sample size $n$, shrinkage of the high-dimensional variance-covariance matrix $\boldsymbol{\Sigma}_{n}$ towards a target, similarly to ridge regression, can reduce a potentially large condition number. This is achieved by reducing the dispersion of the eigenvalues of $\boldsymbol{\Sigma}_{n}$ around its ”grand mean” $\mu_{n}:=d_{n}^{-1}\operatorname{tr}{\boldsymbol{\Sigma}}_{n}$: large eigenvalues are pulled down towards $\mu_{n}$, small eigenvalues are lifted up to $\mu_{n}$. Improvement of the mean-squared error $E[\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}_{n}\|^{2}_{F}]$ is achieved via a potentially tremendous variance reduction (due to stabilisation via regularisation), even if obviously a bias is introduced by adding a deliberately misspecified shrinkage target of low complexity (but high regularity) which usually underfits the true underlying variance-covariance matrix.

A comparatively straightforward, but in practice often well working, choice of the target is a multiple of the $d_{n}-$dimensional identity matrix $I_{n}$ (Ledoit and Wolf (2004), Fiecas et al. (2014)). Other choices consist in specifying, e.g. in a context of an economic time series panel, a given or a latent factor which describes the ”mean-behaviour” of the panel well in terms of a low-dimensional and hence very stable approximation to the high-dimensional panel structure (Ledoit and Wolf (2003), Böhm and von Sachs (2008)). Similarly, adding a parametric estimator of small complexity to the fully nonparametric sample estimator as in Fiecas and von Sachs (2014), follows the same aforementioned paradigm of reducing variance by adding a (model) bias.

The success of this approach, quite naturally, lies in the correct specification of the shrinkage weight $W_{n}$, the proportion with which the shrinkage target enters into the convex combination: obviously it has to be the higher the less regular the given variance-covariance matrix. More specifically, in the above mentioned literature, a theory of optimal choice of $W_{n}$ has been delivered, for various scenario, by minimising the mean squared error between the shrunken estimator $\boldsymbol{\Sigma}_{n}^{s}$ and the true variance-covariance matrix $\boldsymbol{\Sigma}_{n}$. Hence, let

[TABLE]

which shrinks the sample covariance matrix towards the shrinkage target $\mu_{n}I_{n}$. In the population, the optimal shrinkage weight is derived as

[TABLE]

leading to the MSE-optimally shrunken matrix $\boldsymbol{\Sigma}_{n}^{*}=\boldsymbol{\Sigma}_{n}^{s}(W_{n}^{*})$. A closed form solution can be derived as

[TABLE]

where one observes the trade-off between the distance of $\widehat{\boldsymbol{\Sigma}}_{n}$ towards $\boldsymbol{\Sigma}_{n}$ (being large in case of a badly conditioned sample covariance matrix) and the distance of the sample estimator to the shrinkage target. This choice leads to a true improvement on the level of the mean-squared error:

[TABLE]

It is obvious that our distributional results of Section 4 can be directly applied to the shrunken matrix $\boldsymbol{\Sigma}_{n}^{s}$, provided the stationary vector time series ${\mathbf{Y}}_{n1},\dots,{\mathbf{Y}}_{nn}$ satisfying (3.2) and Assumption A, this being the first step in the direction of some inference theory for this kind of shrinkage estimators (which is still lacking in the literature). In practice, the population quantities $\mu_{n}$ and $W^{*}_{n}$ need to be replaced by some estimators. In the situation of disposing of independent copies of the sampled data, a possibility is to use those, quite analogously to Section 5.2 and many other ”statistical learning situations”, in order to construct these estimators $\widehat{\mu}_{n}$ and $\widehat{W}_{n}$. Then, the discussed results on inference on $\boldsymbol{\Sigma}_{n}^{*}$ continue to hold (conditionally on the ”learning sample”).

5.4. Change-point analysis

Change-point analysis is concerned with the detection and analysis of possible structural changes in the distribution of observations and the determination of the time points of their occurence called change-points. For general methodological overviews we refer to Csörgő and Horváth (1997), Steland (2012), the recent review Hušková and Hlávka (2012) and the references given in these papers, amongst others. In view of Theorem 4.1 and Corollary 4.1, we are in a position to study an a posteriori (off-line) change-in-variance problem for the variance of a projection ${\mathbf{w}}_{n}^{\prime}{\mathbf{Y}}_{n}$ due to a change in the covariance structure. Off-line procedures are conducted after having observed the full sample and aim at testing for the presence of a change-point within the sample. In case that such a change-point test rejects the null hypothesis of no change, one is interested in estimating the location of the change-point as well. A test for the change in the covariance matrix has been also proposed by Jarušková (2013), but only for independent Gaussian random vectors. For a fixed number of time series see Aue et al. (2009).

Suppose that under the null hypothesis of no change ${\mathbf{Y}}_{n1},\dots,{\mathbf{Y}}_{nn}$ forms a $d_{n}$–dimensional mean zero stationary vector time series with variance-covariance matrix $\boldsymbol{\Sigma}_{n}^{(0)}$ and satisfying the assumptions of Corollary 4.1. Our change-point model is formulated in terms of the sequence of variances of the projections which are determined by the variance-covariance matrices

[TABLE]

Under the change-point alternative hypothesis we assume that this sequence is equal to $\boldsymbol{\Sigma}_{n}^{(0)}$ up to the change-point $q\in\{1,\dots,n-1\}$ and changes for $i>q$ in such a way that for appropriately selected ${\mathbf{w}}_{n}$ satisfying the uniform $\ell_{1}$–condition (2.2),

[TABLE]

changes from $\sigma_{n0}^{2}={\mathbf{w}}_{n}^{\prime}\boldsymbol{\Sigma}_{n}^{(0)}{\mathbf{w}}_{n}$ to some different value $\sigma_{n1}^{2}\not=\sigma_{n0}^{2}$ and then remains constant again. The applications discussed above in Subsections 5.1 and 5.2 provide examples for the selection of the projection vector. However, it can also be chosen in order to analyze certain elements of the variance-covariance matrix.

An appropriate change-point test statistic directly suggested by our results for the case of known $\boldsymbol{\Sigma}_{n}^{(0)}$ is given by

[TABLE]

cf. (4.1) and (4.2). Corollary 4.1 now provides us with the asymptotic null distribution,

[TABLE]

that is needed to determine critical values in order to devise such a test. Critical values $c_{1-\alpha}$, $\alpha\in(0,1)$, can be easily calculated using the well known explicit formula for the d.f. of $\sup_{t\in[0,1]}|B(t)|$, see e.g. Shorack (2000). If $\alpha_{n}$ is unknown, it can be estimated by $\widehat{\alpha}_{n_{0}}$ from a learning sample of size $n_{0}\geq n$, which satisfies the no-change null hypothesis and the assumptions of Theorem 4.4, using the first $n$ time series. Assuming that such a learning sample is given is, however, standard in the change-point literature, see e.g. Chu et al. (1996) where it has been named non-contamination assumption. We now reject the no–change null hypothesis in favor that a change has occured, if $V_{n}>c_{1-\alpha}$. In this case, the unknown (first) change-point, i.e. the onset, is estimated canonically by

[TABLE]

If $\boldsymbol{\Sigma}_{n}^{(0)}$ is unknown, one may rely on $\mathcal{D}_{n}^{0}$ defined in (4.5) and use the test statistic

[TABLE]

where $B^{0}$ is a Brownian bridge on $[0,1]$. Again, $\alpha_{n}$ can be replaced by our estimator if it is unknown, using a learning sample satisfying the assumptions of Theorem 4.4.

Acknowledgements

Parts of this paper have been written during a research stay of the first author as visiting professor at Université catholique de Louvain, Louvain-la-Neuve, a visit at the Department of Statistics of Seoul National University, Seoul, South Korea, and a sabbatical. He thanks Rainer von Sachs and Sangyeol Lee for their hospitality. Rainer von Sachs gratefully acknowledges funding by contract “Projet d’Actions de Recherche Concertées” No. 12/17-045 of the ,,Communauté française de Belgique” and by IAP research network Grant P7/06 of the Belgian government (Belgian Science Policy). He wants to thank the Institute of Stochastics of the RWTH Aachen (Prof. A. Steland and his group) for the warm hospitality during his sabbatical stay in autumn 2013, and part of this sabbatical has generously been funded by the Université catholique de Louvain. We thank two anonymous referees for helpful comments which improved the presentation of the results.

Proofs

Proof of Lemma 3.1.

The assertions follow from (3.14)–(3.16),

[TABLE]

and (Kouritzin, 1995, Remark 3.2). ∎

Proof of Theorem 4.1.

Notice that we have

[TABLE]

leading to the representation

[TABLE]

with linear processes

[TABLE]

w.r.t. $\{\epsilon_{t}\}$ given by the coefficients

[TABLE]

for $j\geq 0$ and $n\geq 1$. We may now follow the method of proof of Kouritzin (1995), however, we have to take into account that the above processes depend on $n$.

Let $\mathcal{F}_{m}=\sigma(\epsilon_{i}:i\leq m)$, $m\geq 1$. It is easy to check that, for any fixed $n\in\mathbb{N}$, the r.v.s.

[TABLE]

satisfy $E(M_{m}^{(n)}({\mathbf{v}}_{n},{\mathbf{w}}_{n})\,|\,\mathcal{F}_{m-1})=M_{m-1}^{(n)}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$, for $m\geq 0$, thus forming a martingale array $\{M_{m}^{(n)}({\mathbf{v}}_{n},{\mathbf{w}}_{n}):m\in\mathbb{N},n\in\mathbb{N}\}$ with associated martingale differences

[TABLE]

for $n^{\prime},m^{\prime}\geq 0$. Put

[TABLE]

and consider the decomposition

[TABLE]

In order to justify the approximation of $D_{n^{\prime},m^{\prime}}^{(n)}({\mathbf{v}}_{n},{\mathbf{w}}_{n})$ by the martingale differences defined above, it suffices to show that $\sup_{n}E[R_{n^{\prime},m^{\prime}}^{(n)}({\mathbf{v}}_{n},{\mathbf{w}}_{n})]^{2}$ tends to [math] sufficiently fast, as $n^{\prime},m^{\prime}\to\infty$. Using the representation (Kouritzin, 1995, (4.3)), repeating the arguments in Kouritzin (1995) leading to the bounds in (4.8), (4.9) and (4.10) therein and noting that those bounds are uniform in $n\geq 1$, we obtain

[TABLE]

and, for each $n\in\mathbb{N}$,

[TABLE]

This approximation in $L_{2}$ with a rate allows for very general conditions for the validity of a strong approximation. Again, we may follow the arguments given by Kouritzin (1995), by verifying the following sufficient conditions due to Philipp (1986). In terms of an array $\xi_{k}^{(n)}$, $k=1,\dots,n$, of r.v.s., Philipp’s result is as follows. Let $\mathcal{G}_{m}^{(n)}=\sigma(\xi_{ni}:i\leq m)$. If

[TABLE]

satisfies

• (I)

$\|E(S_{n^{\prime},m^{\prime}}^{(n)}|\mathcal{G}_{m^{\prime}}^{(n)})\|_{1}\stackrel{{\scriptstyle m^{\prime},n^{\prime}}}{{<<}}(n^{\prime})^{1/2-\varepsilon}$, a.s., for some $\varepsilon>0$,

• (II)

there exists an $\alpha_{n}^{2}\geq 0$ such that $\|E[(S_{n^{\prime},m^{\prime}}^{(n)})^{2}|\mathcal{G}_{m^{\prime}}^{(n)}]-n^{\prime}\alpha_{n}^{2}\|_{1}\stackrel{{\scriptstyle n^{\prime},m^{\prime}}}{{<<}}(n^{\prime})^{1-\varepsilon}$, a.s., for some $\varepsilon>0$.

• (III)

$\sup_{k\geq 0}E|\xi_{k}^{(n)}|^{4+\delta}<\infty$ for some $\delta>0$,

then there exists a process $\{\widetilde{S}_{n^{\prime}}^{(n)}:n^{\prime}\geq 0\}$ and a standard Brownian motion $\{\widetilde{B}_{t}^{(n)}:t\geq 0\}$ on some probability space $(\widetilde{\Omega},\widetilde{\mathcal{F}},\widetilde{P})$, such that $\{\widetilde{S}_{n^{\prime}}^{(n)}:n^{\prime}\geq 0\}\stackrel{{\scriptstyle d}}{{=}}\{S_{n^{\prime},0}^{(n)}:n^{\prime}\geq 0\}$ and for some $\lambda>0$

[TABLE]

for all $t>0$ $\widetilde{P}$-a.s. Putting, for fixed $n\geq 1$,

[TABLE]

$\mathcal{G}_{m}^{(n)}=\mathcal{F}_{m}$ (since the $\xi_{k}^{(n)}$ are $\mathcal{F}_{k}$–measurable) and repeating the arguments of Kouritzin (1995), we see that, by virtue of Assumption (A), (I)-(III) hold true, which establishes the existence of a standard Brownian motion, $B_{n}(t)$, $t\in[0,\infty)$, such that for some constant $C_{n}$ and some universal $\lambda>0$

[TABLE]

for all $t>0$, a.s. Denoting the standard Brownian motion on $[0,1]$ associated to $B_{n}$, $t\mapsto n^{-1/2}B_{n}(tn)$, $t\in[0,1]$, again by $B_{n}$, we obtain

[TABLE]

which establishes (4.10) and (4.13), provided $C_{n}n^{-\lambda}=o(1)$. It also follows that, for each fixed $n$, the conditional variance of $M_{m^{\prime}+n^{\prime}}^{(n)}-M_{m^{\prime}}^{(n)}$ satisfies

[TABLE]

and (cf. (Kouritzin, 1995, (4.22)))

[TABLE]

as well as

[TABLE]

The constants appearing in (.7), (.8) and (.9) depend on the weighting vectors only through their $\ell_{1}$–norms, as a consequence of eqns (4.22)–(4.26) in Kouritzin (1995) and Lemma 3.1. ∎

Recall that for a $K$–dimensional mean zero random vector ${\mathbf{Z}}$ with variance-covariance matrix ${\mathbf{W}}$, say, the covariance operator $C({\mathbf{u}})=E({\mathbf{u}}^{\prime}{\mathbf{Z}}{\mathbf{Z}})$, ${\mathbf{u}}\in\mathbb{R}^{K}$, can be identified with the linear mapping $\mathbb{R}^{K}\mapsto L(\mathbb{R}^{K};\mathbb{R})$, ${\mathbf{u}}\mapsto{\mathbf{u}}^{\prime}{\mathbf{W}}$, ${\mathbf{u}}\in\mathbb{R}^{K}$, induced by the variance-covariance matrix ${\mathbf{W}}$, where $L(A,B)$ denotes the set of linear mappings $A\to B$.

Proof of Theorem 4.2.

Put

[TABLE]

and notice that

[TABLE]

Also put ${\mathbf{S}}_{n^{\prime},m^{\prime}}^{(n)}=\sum_{k=m^{\prime}+1}^{m^{\prime}+n^{\prime}}\boldsymbol{\xi}_{k}^{(n)}$, $m^{\prime},n^{\prime}\geq 0$. For the Euclidean space $\mathbb{R}^{K}$ equipped with the usual inner product and the induced vector $\ell_{2}$–norm, the conditions (I) and (III) are easily checked. For instance, Jensen’s inequality yields

[TABLE]

for some constant $C$ which does not depend on $n^{\prime},m^{\prime},n$ and $\{({\mathbf{v}}_{nj},{\mathbf{w}}_{nj}):j=1,\dots,K\}$. Introduce the conditional covariance operator

[TABLE]

and the covariance operator

[TABLE]

Noting that $T^{(n)}({\mathbf{u}})=\sum_{j=1}^{K}u_{j}({\mbox{Cov\,}}(B_{n1},B_{nj}),\dots,{\mbox{Cov\,}}(B_{nK},B_{nj}))^{\prime}$, we have to check the remaining condition

• (II)

$E\|(n^{\prime})^{-1}C_{n^{\prime},m^{\prime}}^{(n)}-C^{(n)}\|\stackrel{{\scriptstyle n^{\prime},m^{\prime}}}{{<<}}(n^{\prime})^{-\theta}$, for some covariance operator $C^{(n)}$.

Here the norm is the operator norm defined as $\|L\|=\sup_{{\mathbf{u}}\in\mathbb{R}^{K},\|{\mathbf{u}}\|=1}|{\mathbf{u}}^{\prime}L({\mathbf{u}})|$ for a linear operator $L:\mathbb{R}^{K}\to\mathbb{R}^{K}$. It follows that (II) holds true with $C^{(n)}=T^{(n)}$, if for $i=1,\dots,K$

[TABLE]

and for $i,j=1,\dots,K$ with $i\not=j$

[TABLE]

The last fact follows by a lengthy but straightforward calculation using (3.1).∎

Proof of Theorem 4.3.

The proof is similar to the proof of Theorem 4.1 and Kouritzin (1995) having observed the following crucial facts. We have

[TABLE]

for the linear processes $Y_{i}(\{{\mathbf{v}}_{n}\})=\sum_{j=0}^{\infty}c_{j}^{(v)}\epsilon_{i-j}$, $i\geq 1$, and $Y_{i}(\{{\mathbf{v}}_{n}\})=\sum_{j=0}^{\infty}c_{j}^{(w)}\epsilon_{i-j}$ with coefficients $c_{j}^{(v)}=\sum_{n=1}^{\infty}\lambda_{n}\sum_{\nu=1}^{d_{n}}v_{\nu}c_{nj}^{(\nu)}$ and $c_{j}^{(w)}=\sum_{n=1}^{\infty}\lambda_{n}\sum_{\nu=1}^{d_{n}}w_{\nu}c_{nj}^{(\nu)}$ for $j\geq 0$, which do not depend on the sample size. Since $\sum_{n}|\lambda_{n}|<\infty$,

[TABLE]

i.e. Assumption (A) is satisfied. Hence the result ∎

Proof of Theorem 4.4.

First observe that $Y_{k}^{(v)}Y_{k}^{(w)}$, $k\geq 1$, as well as $Y_{k}^{(v_{r})}Y_{k}^{(w_{r})}Y_{k+h}^{(v_{s})}Y_{k+h}^{(w_{s})}$, $k\geq 1$, are strictly stationary for any fixed $r,s$ and $h$. Their dependence on $d$ will be suppressed in notation. It follows from the proof of Theorem 4.1 and (Kouritzin, 1995, p. 351) that $\alpha^{2}$ can be represented as

[TABLE]

and therefore $\alpha^{2}$ is the long-run variance parameter associated to the time series

[TABLE]

and $\widehat{\alpha}^{2}_{n}$ is the Bartlett type estimator calculated from the first $n$ observations. Analogously, by virtue of (.10),

[TABLE]

where

[TABLE]

for $\ell=1,2,\dots$. Again, we omit the dependence on $d$, but will indicate it for relevant derived quantities when appropriate. Put

[TABLE]

for $h\in\mathbb{Z}$. Observe that $E(Y_{1}^{(v_{r})}Y_{1}^{(w_{r})})E(Y_{1+h}^{(v_{s})}Y_{1+h}^{(w_{s})})=\sigma^{4}\prod_{z\in\{r,s\}}\left(\sum_{j=0}^{\infty}c_{j}^{(v_{z})}c_{j}^{(w_{z})}\right)$. Further, $\Gamma_{h}^{(2)}(r,s)$ is the sum over all terms $c_{j}^{(v_{r})}c_{k}^{(w_{r})}c_{l}^{(v_{s})}c_{m}^{(w_{s})}E(\epsilon_{i-j}\epsilon_{i-k}\epsilon_{i+h-l}\epsilon_{i+h-m})$ such that either all four indices $a=i-j,b=i-k,c=i+h-l,d=i+h-m$ are equal, denoted by $\Gamma_{h,1}^{(2)}(r,s)$, plus the sum of all terms such that $\{a,b,c,d\}=\{A,B\}$ with $A\not=B$, denoted by $\Gamma_{h,2}^{(2)}(r,s)$. Using

[TABLE]

and

[TABLE]

for $j\geq 0$ and $h\geq 1$, the former case leads to $\Gamma_{h,1}^{(2)}(r,s)=E(\epsilon_{1}^{4})\sum_{j=0}^{\infty}c_{j}^{(v_{r})}c_{j}^{(w_{r})}c_{j+h}^{(v_{s})}c_{j+h}^{(w_{s})}=O(h^{-2(1+\delta)})$, uniformly in $d\in\mathbb{N}$ and $1\leq r,s$. To discuss $\Gamma_{h,2}^{(2)}(r,s;d)$ consider the subcase $j=k,l=m,l\not=h+j$ corresponding to the sum

[TABLE]

where $\sup_{1\leq r,s}\sup_{d\in\mathbb{N}}|\overline{\Gamma}_{h}^{(2)}(r,s;d)|<<\sum_{j=0}^{\infty}(j\vee 1)^{-2(1+\delta)}h^{-2(1+\delta)}$, leading to

[TABLE]

The other subcases are straightforward, such that we arrive at

[TABLE]

Hence

[TABLE]

Next introduce the coupling dependence measure

[TABLE]

$p\geq 1$, for a time series $Z_{n}=Z(\epsilon_{n},\epsilon_{n-1},\dots)$, where $Z_{n}^{\prime}=Z(\epsilon_{n},\dots,\epsilon_{1},\epsilon_{0}^{\prime},\epsilon_{-1},\dots)$ with $\epsilon_{0}\stackrel{{\scriptstyle d}}{{=}}\epsilon_{0}^{\prime}$ such that $\epsilon_{0}^{\prime}$ is independent from $\{\epsilon_{k}\}$. Since $Y_{i}^{(v)}$ is a causal linear process with coefficients $c_{j}^{(v)}=\sum_{\nu=1}^{d}v_{\nu}c_{j}^{(\nu)}$, we have

[TABLE]

and by (4.24)

[TABLE]

such that $\sum_{k=0}^{\infty}\delta_{8}(\{Y_{i}^{(v)}\},k)\leq C<\infty$ for constants $c,C<\infty$ not depending on $d$ and uniformly over $\|{\mathbf{v}}\|_{\ell_{1}}\leq C_{v,w}$. Further, (using $[E(X^{4}Y^{4})]^{1/4}\leq[(EX^{8})^{1/2}(EY^{8})^{1/2}]^{1/4}=[EX^{8}]^{1/8}[EY^{8}]^{1/8}$),

[TABLE]

leading to $\sum_{k=0}^{\infty}\delta_{4}(\{\xi_{i}\},k)\leq C_{1}<\infty$ for some constant $C_{1}$, uniformly over $d\in\mathbb{N}$ and $\|{\mathbf{v}}\|_{\ell_{1}},\|{\mathbf{w}}\|_{\ell_{1}}\leq C_{v,w}$. Lastly,

[TABLE]

leading to $\sum_{k=0}^{\infty}\delta_{2}(\{\xi_{i}(r)\xi_{i+h}(s)\},k)<C_{2}$ for a constant $C_{2}$, uniformly over $d\in\mathbb{N}$ and $1\leq r,s$. Define

[TABLE]

By virtue of (Wu, 2007, Th. 1), we obtain

[TABLE]

for some constant $C_{3}<\infty$ not depending on $h$ or $m$, uniformly over $\|{\mathbf{v}}\|_{\ell_{1}}\leq C_{v,w}$, such that

[TABLE]

for some constant $C_{4}<\infty$. Observe that $\widehat{\Gamma}_{h}(r,s)=\frac{1}{n}\sum_{i=1}^{n-h}(\xi_{i}(r)-\overline{\xi}_{n}(r))(\xi_{i+h}(s)-\overline{\xi}_{n}(s))$, where $\overline{\xi}_{n}(\ell)=n^{-1}\sum_{j=1}^{n}\xi_{j}(\ell)$, $\ell=1,2,\dots$, and

[TABLE]

Using the Cauchy-Schwarz inequality, we thus obtain

[TABLE]

uniformly over $h\leq m$, $d\in\mathbb{N}$ and $1\leq r,s$. It follows that

[TABLE]

which implies

[TABLE]

as $n\to\infty$. Hence it suffices to show the result for

[TABLE]

Using the representation $\beta^{2}(r,s;d)=\sum_{h\in\mathbb{Z}}\Gamma_{h}(r,s;d)$, we obtain the decomposition

[TABLE]

where

[TABLE]

First observe that $\sup_{1/ler,s}\sup_{d\in\mathbb{N}}|D_{n}(r,s;d)|=o(1)$ by (.11), and, by Fubini and (.12),

[TABLE]

as $n\to\infty$, since $m^{2}/n=o(1)$ by assumption. Here $d\pi$ denotes the counting measure on $\mathbb{Z}$. Further,

[TABLE]

where the integrand is $o(1)$ point-wise in $h$ and bounded by the (uniformly over $d\in\mathbb{N}$ and $1\leq r,s$) $\pi$–integrable function $2W|\Gamma_{h}(r,s;d)|$, such that $\sup_{1\leq r,s}\sup_{d\in\mathbb{N}}E|B_{n}(r,s;d)|=o(1)$, as $n\to\infty$, follows. Similiarly, $\sup_{1\leq r,s}E|C_{n}(r,s;d)|=o(1)$, as $n\to\infty$, uniformly in $d\in\mathbb{N}$, by (.11) and (W1). Hence the assertion follows. ∎