Toeplitz Covariance Models for Ordered Measurements

Why Model Covariance by Lag?

For measurements recorded in a meaningful order—successive visits, equally spaced times, or positions along a sequence—pairs of measurements may be related according to how far apart they are. A covariance model that assigns the same covariance to every pair separated by the same number of positions expresses that idea compactly. It replaces a general covariance matrix, whose entries may all differ, with a small collection of lag-specific quantities.

This restriction is useful when the order is scientifically meaningful and there is reason to treat locations similarly across the sequence. It is a modeling assumption, not a consequence of ordering alone. In particular, the assumption concerns covariance; it does not by itself specify the mean vector or the full distribution.

Definition: A $p\times p$ covariance matrix is Toeplitz if its entries are constant along each diagonal: $\Sigma_{ij}$ depends only on $i-j$. For a real covariance matrix, symmetry gives the equivalent representation $$ \Sigma_{ij}=\gamma_{|i-j|},\qquad i,j=1,\ldots,p, $$ where $\gamma_k$ is the covariance at lag $k$, for $k=0,\ldots,p-1$. In particular, $\gamma_0=\operatorname{Var}(X_i)$ is the common marginal variance.

Writing out the matrix makes the restriction visible:

\Sigma= \begin{pmatrix} \gamma_0&\gamma_1&\gamma_2&\cdots&\gamma_{p-1}\\ \gamma_1&\gamma_0&\gamma_1&\cdots&\gamma_{p-2}\\ \gamma_2&\gamma_1&\gamma_0&\cdots&\gamma_{p-3}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \gamma_{p-1}&\gamma_{p-2}&\gamma_{p-3}&\cdots&\gamma_0 \end{pmatrix}.

There are $p$ lag-specific entries, rather than $p(p+1)/2$ free entries in an unrestricted symmetric covariance matrix. The reduction can be substantial for a long sequence, but the parameters are not unconstrained: a covariance matrix used in a nonsingular Gaussian model must be positive definite. The precise restrictions on the $\gamma_k$'s depend on $p$; the next topic treats positive definiteness of finite Toeplitz matrices.

Worked Example: Reading Lag Parameters From a Matrix

Suppose four ordered measurements have marginal variance $2$, covariance $0.8$ at lag one, covariance $0.3$ at lag two, and covariance $-0.1$ at lag three. The Toeplitz matrix is

\Sigma= \begin{pmatrix} 2&0.8&0.3&-0.1\\ 0.8&2&0.8&0.3\\ 0.3&0.8&2&0.8\\ -0.1&0.3&0.8&2 \end{pmatrix}.

For example, $\operatorname{Cov}(X_1,X_3)=0.3$ and $\operatorname{Cov}(X_2,X_4)=0.3$, because both pairs are separated by two positions. The matrix has four distinct lag entries, compared with ten entries needed to specify a general symmetric $4\times4$ matrix.

The listed values specify a symmetric Toeplitz matrix, but that fact alone does not verify that it is a valid covariance matrix. A covariance matrix must also be positive semidefinite, and a nonsingular Gaussian likelihood requires positive definiteness. Matrix pattern and covariance validity are separate checks.

Toeplitz Structure Is Not the Whole Stationarity Assumption

For a weakly stationary process, the mean is constant over time and the covariance between two observations depends only on their separation. Consequently, the covariance matrix of any finite consecutive block is Toeplitz. The converse needs care: saying that the covariance matrix of one finite vector is Toeplitz specifies only its second moments within that vector. It does not assert a constant mean, define observations outside the vector, or impose stationarity on all finite-dimensional distributions of a process.

For a Gaussian vector, the mean vector and covariance matrix together determine its distribution. Thus, a Gaussian vector with constant mean and Toeplitz covariance has a distributional symmetry for the shifts that keep coordinates within the observed vector. That finite-dimensional statement is still not the same as specifying a stationary process indexed over an extended range. The key distinction is that “Toeplitz” describes a matrix pattern, whereas “stationary” describes a process (or a collection of distributions) under shifts.

Theorem: For a real random vector $X=(X_1,\ldots,X_p)^\mathsf{T}$ with finite second moments, its covariance matrix is Toeplitz if and only if there are numbers $\gamma_0,\ldots,\gamma_{p-1}$ such that $\operatorname{Cov}(X_i,X_j)=\gamma_{|i-j|}$ for every $i,j$.

Proof. If the covariance matrix is Toeplitz, entries on each diagonal are equal. By symmetry of covariance, the entries on the $k$th diagonal above the main diagonal equal those on the corresponding diagonal below it. Call their common value $\gamma_k$, and call the common main-diagonal value $\gamma_0$. Then $\operatorname{Cov}(X_i,X_j)=\gamma_{|i-j|}$. Conversely, if this lag representation holds, any two matrix entries on the same diagonal have the same absolute index difference and therefore the same value. The matrix is Toeplitz. $\square$

Worked Example: Toeplitz Covariance Does Not Imply Constant Means

Let $X=(X_1,X_2,X_3)^\mathsf{T}$ have mean vector $(0,1,2)^\mathsf{T}$ and covariance matrix

\Sigma= \begin{pmatrix} 1&0.2&0\\ 0.2&1&0.2\\ 0&0.2&1 \end{pmatrix}.

The covariance matrix is Toeplitz: the diagonal entries are all $1$, the lag-one entries are all $0.2$, and the lag-two entries are both $0$. But the mean changes from $0$ to $1$ to $2$. Therefore this vector does not have a constant mean, and Toeplitz covariance alone is not the full weak-stationarity condition.

This example concerns the logical distinction between the assumptions; it does not claim that every arbitrary Toeplitz pattern is a valid covariance matrix. In this particular case, positive definiteness can be checked directly: the leading principal minors are $1$, $1-0.2^2=0.96$, and $1-2(0.2^2)=0.92$, all positive.

Gaussian Likelihood for a Finite Ordered Vector

Let $X$ be a $p$-dimensional Gaussian vector with specified mean $\mu$ and covariance $T_p(\gamma)$, the Toeplitz matrix formed from $\gamma=(\gamma_0,\ldots,\gamma_{p-1})^\mathsf{T}$. Assume $T_p(\gamma)$ is positive definite. For an observed vector $x$, write $r=x-\mu$. Its log-likelihood as a function of the lag parameters is

\ell(\gamma;x) =-\frac{p}{2}\log(2\pi) -\frac{1}{2}\log|T_p(\gamma)| -\frac{1}{2}r^\mathsf{T}T_p(\gamma)^{-1}r.

The three terms are, respectively, the normalizing constant, the covariance-volume penalty, and the quadratic penalty for the observed residual. Maximizing this expression over the admissible parameter values gives the Gaussian maximum-likelihood estimate. If the mean is unknown, it must also be included in the likelihood and estimated or otherwise handled; it is not automatically made constant by the Toeplitz covariance restriction.

For $m$ independent replicate vectors $x^{(1)},\ldots,x^{(m)}$ with the same mean and covariance, set $r^{(a)}=x^{(a)}-\mu$. Independence makes the joint likelihood a product, so the log-likelihood is

\ell(\gamma) =-\frac{mp}{2}\log(2\pi) -\frac{m}{2}\log|T_p(\gamma)| -\frac{1}{2}\sum_{a=1}^{m} \bigl(r^{(a)}\bigr)^\mathsf{T}T_p(\gamma)^{-1}r^{(a)}.

A single vector can be used to write a likelihood, but it may provide limited information for estimating many covariance parameters. Replicates contribute additional quadratic-form information while sharing the same covariance parameter. In either case, maximum likelihood is a constrained optimization over covariance matrices in the Toeplitz family, not an unconstrained optimization over arbitrary lag values.

Theorem: Suppose $x^{(1)},\ldots,x^{(m)}$ are independent $N(\mu,T_p(\gamma))$ vectors, with known $\mu$, and $T_p(\gamma)$ is positive definite and differentiable in each parameter $\gamma_k$ on the parameter region considered. Let $D_k=\partial T_p(\gamma)/\partial\gamma_k$. Then the score for $\gamma_k$ is $$ \frac{\partial\ell}{\partial\gamma_k} =-\frac{m}{2}\operatorname{tr}\!\left(T_p^{-1}D_k\right) +\frac{1}{2}\sum_{a=1}^{m} \bigl(r^{(a)}\bigr)^\mathsf{T}T_p^{-1}D_kT_p^{-1}r^{(a)}. $$

Proof. Write $T=T_p(\gamma)$, which is invertible by positive definiteness. The log-likelihood is the expression above. For a differentiable invertible matrix, differentiating $TT^{-1}=I$ gives $d(T^{-1})=-T^{-1}(dT)T^{-1}$. Also, differentiating the determinant gives $d\log|T|=\operatorname{tr}(T^{-1}dT)$. Apply these identities with $dT=D_k\,d\gamma_k$, holding the residual vectors fixed because $\mu$ is known. The derivative of the determinant term is $-\frac{m}{2}\operatorname{tr}(T^{-1}D_k)$. For each quadratic term, the inverse-matrix identity gives

-\frac{1}{2}\bigl(r^{(a)}\bigr)^\mathsf{T} \frac{\partial T^{-1}}{\partial\gamma_k}r^{(a)} = \frac{1}{2}\bigl(r^{(a)}\bigr)^\mathsf{T} T^{-1}D_kT^{-1}r^{(a)}.

Adding the determinant contribution and all $m$ quadratic contributions yields the stated score. $\square$

Here $D_k$ has ones at positions with $|i-j|=k$ and zeros elsewhere; for $k=0$, its ones are on the main diagonal. The score therefore describes how the likelihood changes when an entire lag diagonal is varied together. A maximum in the interior of the admissible parameter region must have zero score in every free parameter, though boundary maxima must also be considered.

Worked Example: Comparing Two Gaussian Likelihood Values

Suppose $X=(X_1,X_2)^\mathsf{T}$ has known mean zero and candidate covariance \[ \Sigma(\rho)= \begin{pmatrix} 2&\rho\\ \rho&2 \end{pmatrix}, \] and the observed vector is $x=(1,2)^\mathsf{T}$. Compare $\rho=0$ with $\rho=1$. Both matrices are positive definite: their eigenvalues are $2+\rho$ and $2-\rho$, positive for these two candidates.

At $\rho=0$, $|\Sigma|=4$ and $x^\mathsf{T}\Sigma^{-1}x=(1^2+2^2)/2=2.5$. Hence

\ell(0;x)=-\log(2\pi)-\frac12\log(4)-\frac12(2.5) \approx -3.7810.

At $\rho=1$, the determinant is $4-1=3$, and \[ \Sigma^{-1}=\frac{1}{3} \begin{pmatrix} 2&-1\\ -1&2 \end{pmatrix}. \] Multiplying gives $\Sigma^{-1}x=(0,1)^\mathsf{T}$, so the quadratic form is $x^\mathsf{T}\Sigma^{-1}x=2$. Thus

\ell(1;x)=-\log(2\pi)-\frac12\log(3)-\frac12(2) \approx -3.3872.

Among these two candidates, $\rho=1$ has the larger log-likelihood. This comparison does not establish that $\rho=1$ maximizes likelihood over all admissible values; it only compares the two specified matrices.

Worked Example: Maximizing a Likelihood With Two Replicates

Let two independent observations be $x^{(1)}=(1,0)^\mathsf{T}$ and $x^{(2)}=(0,1)^\mathsf{T}$, each with known mean zero and covariance \[ \Sigma(\rho)= \begin{pmatrix} 1&\rho\\ \rho&1 \end{pmatrix},\qquad -1<\rho<1. \] The determinant is $1-\rho^2$, and \[ \Sigma(\rho)^{-1}=\frac{1}{1-\rho^2} \begin{pmatrix} 1&-\rho\\ -\rho&1 \end{pmatrix}. \] Each observed vector has quadratic form $1/(1-\rho^2)$; their sum is $2/(1-\rho^2)$. The log-likelihood is therefore

\ell(\rho) =-2\log(2\pi)-\log(1-\rho^2)-\frac{1}{1-\rho^2}.

Set $t=1-\rho^2$, so $0<t\leq1$. Apart from the constant $-2\log(2\pi)$, the objective is $g(t)=-\log t-1/t$. Its derivative is \[ g'(t)=-\frac{1}{t}+\frac{1}{t^2}=\frac{1-t}{t^2}\geq0 \] for $0<t\leq1$. Thus $g$ is increasing in $t$ on this interval, and its maximum occurs at $t=1$, which is equivalent to $\rho=0$. The maximum-likelihood estimate is $\hat{\rho}=0$. As a check, the log-likelihood there is $-2\log(2\pi)-1\approx-4.6758$.

Interpretation and Common Pitfalls

A Toeplitz model is most useful when the sequence order supports a common lag interpretation. Its parameter savings can make estimation more manageable, but the assumption also forces, for example, the covariance of measurements $1$ and $3$ to equal that of measurements $2$ and $4$. If the measurement process changes across the sequence, this equality may be scientifically inappropriate.

Mistake: assuming that ordered data are automatically Toeplitz. Order alone does not impose equal covariances at equal lags. The model must be justified or assessed.
Mistake: treating Toeplitz covariance as equivalent to stationarity. Toeplitz describes the covariance matrix of a finite vector. Stationarity also concerns the mean and, depending on its definition, the distributions of a process under shifts.
Mistake: maximizing over arbitrary lag values. A symmetric Toeplitz pattern is not sufficient for a valid nonsingular Gaussian covariance. The likelihood domain must contain only positive-definite matrices.
Mistake: interpreting a larger likelihood at one candidate as a global maximum. A comparison between two parameter values is not an optimization over the full admissible region.

The practical workflow is to specify which observations form the ordered vector, decide whether a common mean is part of the model, define the lag parameters, and then maximize the Gaussian likelihood over the valid covariance family. The likelihood score can guide numerical optimization, but the matrix constraint remains part of the model throughout.

Check Your Understanding

Use the definitions and likelihood calculations above to answer the following questions.

How many distinct lag parameters specify a $p\times p$ real symmetric Toeplitz covariance matrix?
In a four-measurement sequence, which covariance entries must be equal at lag two?
Why does Toeplitz covariance alone not imply a constant mean?
What condition on $T_p(\gamma)$ is needed for the nonsingular Gaussian likelihood written here?
For independent replicate vectors with known mean, what changes in the log-likelihood compared with one vector?
Why does comparing the likelihood at two candidate parameter values not by itself identify the maximum-likelihood estimate?