Lag-Based Moment Estimation for Toeplitz Covariance

Why Estimate Covariance by Lag?

For ordered measurements, a Toeplitz covariance model assigns the same covariance to every pair of coordinates separated by the same lag. This structure can replace many unrelated covariance parameters with a short sequence of lag-specific quantities. Estimation then becomes a pooling problem: combine all observed products at lag zero, all products at lag one, and so on.

Pooling equal-lag products is natural, but it has two separate consequences that must not be confused. A lag estimate can be unbiased for its target and still be noisy; and a collection of individually reasonable lag estimates need not assemble into a positive-definite covariance matrix. We derive the estimators and their finite-sample bias before examining that matrix constraint.

Definition: Let \(Y_r=(Y_{r1},\ldots,Y_{rp})^T\), \(r=1,\ldots,M\), be independent and identically distributed \(p\)-dimensional random vectors with finite second moments. Suppose their covariance matrix is Toeplitz: \[ \operatorname{Cov}(Y_{ri},Y_{rj})=\gamma_{|i-j|},\qquad 1\leq i,j\leq p. \] The number \(\gamma_h\) is the covariance at lag \(h\), for \(h=0,\ldots,p-1\). The mean may vary by coordinate; write it as \(\mu_i=E[Y_{ri}]\).

The target covariance matrix is therefore \(\Gamma_p=(\gamma_{|i-j|})_{i,j=1}^p\). At lag \(h\), there are \(p-h\) coordinate pairs, namely \((i,i+h)\) for \(i=1,\ldots,p-h\). A lag-based estimator averages products or sample covariances over these pairs and over the \(M\) replicate vectors.

Estimators and Their Finite-Sample Bias

If the mean vector is known, center each observation using its corresponding coordinate mean. The natural estimator averages over all \(M(p-h)\) products at lag \(h\). If the mean is unknown, first estimate the covariance between each coordinate pair across the independent replicate vectors, then average those sample covariances over positions at the same lag.

Definition: With known means, define \[ \widetilde{\gamma}_h= \frac{1}{M(p-h)} \sum_{r=1}^{M}\sum_{i=1}^{p-h} (Y_{ri}-\mu_i)(Y_{r,i+h}-\mu_{i+h}). \] With unknown means and \(M\geq 2\), let \(\bar Y_i=M^{-1}\sum_{r=1}^M Y_{ri}\) and \[ S_{ij}=\frac{1}{M-1}\sum_{r=1}^M (Y_{ri}-\bar Y_i)(Y_{rj}-\bar Y_j). \] The lag-averaged sample covariance is \[ \widehat{\gamma}_h=\frac{1}{p-h}\sum_{i=1}^{p-h}S_{i,i+h}. \]

These formulas assume the replicate vectors are independent and identically distributed, with the stated common covariance across replicates. Gaussianity is not needed. Finite second moments suffice for the expectations below. In particular, the coordinate pairs within one vector need not be independent.

Theorem: Under the conditions above, for every \(h=0,\ldots,p-1\), \[ E[\widetilde{\gamma}_h]=\gamma_h. \] For \(M\geq2\), the unknown-mean estimator also satisfies \[ E[\widehat{\gamma}_h]=\gamma_h. \] Thus both are unbiased lag-by-lag estimators.

Proof. For the known-mean estimator, each product has expectation \[ E[(Y_{ri}-\mu_i)(Y_{r,i+h}-\mu_{i+h})] =\operatorname{Cov}(Y_{ri},Y_{r,i+h})=\gamma_h. \] There are exactly \(M(p-h)\) such products. Linearity of expectation therefore gives \(E[\widetilde{\gamma}_h]=\gamma_h\).

For the estimated-mean version, fix a coordinate pair \(i,j\), and regard the \(M\) observations of these coordinates as paired replicate data. The identity \[ \sum_{r=1}^M(Y_{ri}-\bar Y_i)(Y_{rj}-\bar Y_j) = \sum_{r=1}^M(Y_{ri}-\mu_i)(Y_{rj}-\mu_j) -M(\bar Y_i-\mu_i)(\bar Y_j-\mu_j) \] follows by expanding both sides and using \(\sum_r(Y_{ri}-\bar Y_i)=0\). The expectation of the first sum on the right is \(M\operatorname{Cov}(Y_{ri},Y_{rj})\). Independence across replicates gives \[ E[(\bar Y_i-\mu_i)(\bar Y_j-\mu_j)] =\frac{1}{M}\operatorname{Cov}(Y_{ri},Y_{rj}), \] so the expected second term is one copy of that covariance. The expected left side is consequently \((M-1)\operatorname{Cov}(Y_{ri},Y_{rj})\), and \(E[S_{ij}]=\operatorname{Cov}(Y_{ri},Y_{rj})\). Averaging this identity over the \(p-h\) pairs at lag \(h\) proves \(E[\widehat{\gamma}_h]=\gamma_h\). \(\square\)

The divisor matters. If the sample means are used for centering but the covariance sum is divided by \(M\) instead of \(M-1\), its expectation is \((M-1)\operatorname{Cov}(Y_{ri},Y_{rj})/M\). The resulting lag estimate has expectation \((M-1)\gamma_h/M\) and bias \(-\gamma_h/M\). This bias applies to every lag, including lag zero; it is not removed by averaging over coordinate pairs.

Worked Example: Pooling Lag Products With a Known Mean

Let \(p=3\), \(M=2\), and suppose the known mean is zero. The two observed vectors are \(Y_1=(1,0,1)^T\) and \(Y_2=(-1,0,-1)^T\). Calculate the known-mean lag estimates and assemble their Toeplitz matrix.

At lag zero, the six squared coordinates sum to \(2+2=4\), so \[ \widetilde{\gamma}_0=\frac{4}{2\cdot3}=\frac{2}{3}. \] At lag one, the products within the two vectors are \(1\cdot0+0\cdot1=0\) and \((-1)\cdot0+0\cdot(-1)=0\). There are \(2(3-1)=4\) such products, giving \(\widetilde{\gamma}_1=0/4=0\).

At lag two, the single product per vector is \(1\cdot1=1\) and \((-1)\cdot(-1)=1\). Hence \(\widetilde{\gamma}_2=2/[2(1)]=1\). The assembled estimate is \[ \widehat{\Gamma}_{\mathrm{known}}= \begin{pmatrix} 2/3&0&1\\ 0&2/3&0\\ 1&0&2/3 \end{pmatrix}. \] The vector \((1,0,-1)^T\) is an eigenvector with eigenvalue \(2/3-1=-1/3\). Thus this estimate is not positive definite, despite each lag estimate being unbiased for its corresponding covariance.

A Different Divisor for a Single Series

A related finite-sample issue appears when the observations are one stationary series rather than independent replicate vectors. Let \(X_1,\ldots,X_n\) be consecutive observations from a weakly stationary process with mean zero and autocovariance \(\gamma_h=\operatorname{Cov}(X_t,X_{t+h})\). Here “weakly stationary” means the mean is constant and the covariance depends only on the lag. If the mean is known to be zero, the lag-\(h\) product sum contains \(n-h\) terms, not \(n\).

Theorem: For \(0\leq h<n\), define \[ \check{\gamma}_h=\frac{1}{n-h}\sum_{t=1}^{n-h}X_tX_{t+h}. \] Under weak stationarity and finite second moments, \(E[\check{\gamma}_h]=\gamma_h\). If instead the sum is divided by \(n\), then its expectation is \((n-h)\gamma_h/n\), so its bias is \(-h\gamma_h/n\).

Proof. For every \(t=1,\ldots,n-h\), weak stationarity gives \(E[X_tX_{t+h}]=\operatorname{Cov}(X_t,X_{t+h})=\gamma_h\), since the mean is zero. Summing the \(n-h\) expectations and dividing by \(n-h\) yields \(E[\check{\gamma}_h]=\gamma_h\). Dividing the same sum by \(n\) instead gives expectation \((n-h)\gamma_h/n\); subtracting \(\gamma_h\) gives the stated bias. No independence among the terms is needed for these expectation calculations. \(\square\)

Worked Example: The Single-Series Divisor Effect

Take \(n=4\), lag \(h=1\), and observed values \(X_1=1,X_2=2,X_3=0,X_4=-1\), with the mean known to be zero. The three adjacent products are \(1\cdot2=2\), \(2\cdot0=0\), and \(0\cdot(-1)=0\), whose sum is 2.

The divisor \(n-h=3\) gives \(\check{\gamma}_1=2/3\). Dividing by \(n=4\) instead gives \(2/4=1/2\). The expectation of the latter estimator is \((3/4)\gamma_1\), with bias \(-\gamma_1/4\). The observed difference between the two numerical values is a divisor effect; it does not establish the true value of \(\gamma_1\) from this one realization.

Unbiased Lag Estimates May Violate Positive Definiteness

A covariance matrix must be positive semidefinite: for every vector \(a\), the variance of \(a^TY_r\) is \(a^T\Gamma_pa\geq0\). A valid nonsingular covariance matrix must be positive definite, meaning \(a^T\Gamma_pa>0\) for every nonzero \(a\). The unbiasedness results above concern each lag separately; they impose no constraint ensuring that a particular realized vector of lag estimates satisfies all the joint matrix inequalities.

Worked Example: Unbiased Estimates Can Give an Indefinite Matrix

For the known-mean sample in the first example, the estimated lag sequence is \((2/3,0,1)\). The Toeplitz matrix is \[ \begin{pmatrix} 2/3&0&1\\ 0&2/3&0\\ 1&0&2/3 \end{pmatrix}. \] The vector \(v=(1,0,-1)^T\) gives \(v^T\widehat{\Gamma}_{\mathrm{known}}v=2(2/3)-2(1)=-2/3<0\). Equivalently, \(v\) is an eigenvector with eigenvalue \(-1/3\). This directly violates positive semidefiniteness, so the matrix cannot be a covariance matrix.

There is no contradiction with unbiasedness. The expectation of each estimated entry equals the corresponding true covariance entry, but a finite-sample realization can fall outside the set of valid covariance matrices. Even an estimator whose expectation is positive definite need not itself be positive definite for every data set.

In the replicate setting, the full sample covariance matrix \(S=(S_{ij})\) is positive semidefinite. However, the Toeplitz estimate replaces all entries on each diagonal by their average. Such averaging need not preserve positive semidefiniteness: it changes the matrix, and positive semidefiniteness is a joint restriction on all entries, not a separate restriction on each lag. The explicit eigenvalue check above is a simple diagnostic in a small problem; for larger matrices, one can check the smallest eigenvalue or attempt a Cholesky factorization.

Practical Interpretation and Common Pitfall

Lag pooling can reduce the number of distinct quantities to estimate and is useful when the ordered coordinates are plausibly governed by a common covariance pattern. Its finite-sample behavior still depends on the design. For replicated vectors, there are \(p-h\) coordinate pairs at lag \(h\), so fewer positions contribute at larger lags. For a single series of length \(n\), there are \(n-h\) products. These counts identify the number of terms in the average, but overlapping pairs are generally dependent; one should not treat the terms as independent observations when calculating uncertainty.

Pitfall: “Each lag estimate is unbiased, so the Toeplitz matrix is a valid covariance estimate.” Unbiasedness is an expectation statement, not a positive-definiteness guarantee. Check the assembled matrix itself.
Pitfall: using the wrong divisor without stating the convention. Centered covariance across \(M\) independent replicates uses \(M-1\) for unbiasedness. For a known-mean series at lag \(h\), the unbiased divisor is \(n-h\), not \(n\).
Pitfall: confusing replication with within-series observations. Independent replicate vectors and consecutive observations from one stationary series have different sampling dependence. Their formulas may look similar, but their assumptions and uncertainty calculations are not interchangeable.

Key takeaway: Lag-based moment estimation averages the covariance information at each shared lag. Choose the divisor to match the sampling setup, state whether the mean is known or estimated, and check positive definiteness after assembling the Toeplitz matrix. Unbiased lag estimates alone do not guarantee a feasible covariance matrix.

Check Your Understanding

Use the assumptions stated in each question, and distinguish an expectation statement from a claim about every realized estimate.

For \(M\) independent replicate vectors with known means, how many centered products enter the lag-\(h\) estimator, and why is its divisor \(M(p-h)\)?
Explain why replacing \(M-1\) by \(M\) after estimating the replicate means produces bias \(-\gamma_h/M\).
For a mean-zero stationary series, what is the expectation of the lag-\(h\) product sum divided by \(n\)?
Why does unbiasedness of every \(\widehat{\gamma}_h\) not imply that the assembled Toeplitz estimate is positive semidefinite?
What calculation involving a vector \(a\) can certify that a proposed covariance estimate is not positive semidefinite?