Data analysis for the physical sciences

Methods for analysing fields and time series.

In this chapter 6

Tick a unit to add it to your PDF basket.

Fields and time series

Fields and time series

Method

Principal components and empirical orthogonal functions

Adapted

Add to PDF

Builds the sample covariance matrix of a centred data matrix, shows that its eigenvectors are the directions of successively largest sample variance and that they are the right singular vectors of the data matrix itself, reads the eigenvalues as the variance each direction explains, and identifies the truncated expansion as the best approximation of its rank — which is what an EOF analysis of a space-time field computes.

Prerequisites: Singular value decomposition, Spectral theorem for Hermitian matrices, Quadratic forms and the Rayleigh quotient, Orthogonal projection and least squares, Expectation, variance and covariance, The multivariate normal distribution

Note(What this unit is, and what is ours)Original

Principal component analysis takes a table of numbers — $n$ observations of $p$ variables — and finds the few directions in the $p$ -dimensional variable space along which the data actually vary. In the atmospheric and ocean sciences the same computation, applied to a field sampled at $p$ grid points on $n$ occasions, is called empirical orthogonal function (EOF) analysis: the eigenvectors are maps, and the coefficients that multiply them are time series. The two names describe one piece of linear algebra.

None of the books in this handbook's library treats principal components under that name. What they supply is the machinery: the singular value decomposition and the variational property of the eigenvectors of a symmetric matrix from Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §8.17.1 and §8.18.3, pp. 290-291 and 302-304), and the covariance matrix of a random vector from Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.9, pp. 132-134). The assembly of those pieces into the statement below, and every number in the example, are ours.

Definition(Data matrix, sample covariance, principal components, EOFs)Original

Let $X$ be a real $n \times p$ data matrix with $n \ge 2$ : row $t$ is one observation $x_t \in \mathbb{R}^p$ , column $j$ is one variable. For a space-time field, $t$ indexes time and $j$ indexes a grid point, so $X_{tj}$ is the value of the field at grid point $j$ at time $t$ .

The column mean is $\bar x = n^{-1}\sum_{t=1}^{n} x_t$ , and $X$ is centred when $\bar x = 0$ , that is when every column sums to zero. We assume from here on that $X$ has been centred by subtracting $\bar x$ from every row; the matrix before centring is called the raw data and is never used below.

The sample covariance matrix of a centred $X$ is the $p \times p$ matrix

S \;=\; \frac{1}{n-1}\,X^{\mathsf T}X , \qquad S_{jk} \;=\; \frac{1}{n-1}\sum_{t=1}^{n} X_{tj}X_{tk} .

It is symmetric and positive semidefinite, since $a^{\mathsf T}Sa = \lVert Xa\rVert^{2}/(n-1) \ge 0$ for every $a \in \mathbb{R}^p$ . Its diagonal entries are the sample variances of the $p$ variables, and $\operatorname{tr}S$ is the total variance of the data set.

Given a unit vector $a \in \mathbb{R}^p$ , the numbers $z_t = x_t \cdot a$ ( $t = 1,\dots,n$ ) are the data projected onto the direction $a$ . Their mean is zero, because $X$ is centred, and their sample variance is

\frac{1}{n-1}\sum_{t=1}^{n} (x_t\cdot a)^{2} \;=\; \frac{\lVert Xa\rVert^{2}}{n-1} \;=\; a^{\mathsf T}Sa .

A principal direction (in the geosciences, an empirical orthogonal function, or EOF) is a unit vector $e_k \in \mathbb{R}^p$ produced by the following rule: $e_1$ maximises $a^{\mathsf T}Sa$ over all unit $a$ , and for $k \ge 2$ , $e_k$ maximises $a^{\mathsf T}Sa$ over all unit $a$ orthogonal to $e_1, \dots, e_{k-1}$ . The corresponding principal component (in the geosciences, the PC time series, or expansion coefficient) is the column vector

z_k \;=\; X e_k \;\in\; \mathbb{R}^{n} ,

whose $t$ th entry is the amount of pattern $e_k$ present in observation $t$ . The sample variance of $z_k$ is $e_k^{\mathsf T}Se_k$ , written $\lambda_k$ , and, provided $X \ne 0$ — equivalently $\operatorname{tr}S > 0$ , since $\operatorname{tr}S = \lVert X\rVert_F^2/(n-1)$ — the explained variance fraction of mode $k$ is $\lambda_k/\operatorname{tr}S$ . For the zero data matrix every $\lambda_k$ vanishes and the fractions are not defined; that case is excluded wherever they are used.

The population counterpart of $S$ is the covariance matrix $V$ of a random vector, $V_{jk} = \operatorname{cov}(X_j, X_k)$ . Grimmett and Stirzaker (Grimmett & Stirzaker (2020), Thm (4.9.5) and Exercise (4.9.3), pp. 132, 134) state this naming of $V$ at their Theorem (4.9.5), and give the variance $\sum_i a_i^2\operatorname{var}(X_i) + 2\sum_{i<j}a_ia_j\operatorname{cov}(X_i,X_j)$ of a linear combination at their Exercise (4.9.3), which is the population form of the identity $\operatorname{var}(a\cdot X) = a^{\mathsf T}Va$ just computed for a sample. Both are stated there for a vector with the multivariate normal distribution $N(\mu, V)$ , which is the only distribution that section treats; neither uses normality, since both are identities between second moments that hold for any random vector with finite variances. Everything below is in any case a statement about the matrix $S$ actually formed from the data, not about $V$ .

Method(Principal components, the SVD, and optimal truncation)Original

Let $X$ be a centred real $n \times p$ data matrix with $n \ge 2$ and sample covariance $S = X^{\mathsf T}X/(n-1)$ , and let $r = \operatorname{rank}X$ .

1. The principal directions exist and are the eigenvectors of $S$ . Write the eigenvalues of $S$ , which are real and non-negative, in decreasing order $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0$ , with an orthonormal basis $e_1,\dots,e_p$ of eigenvectors, $Se_k = \lambda_ke_k$ . Then this $e_1,\dots,e_p$ is a valid choice of principal directions, the successive maxima are

\max_{\lVert a\rVert = 1,\; a \perp e_1,\dots,e_{k-1}} a^{\mathsf T}Sa \;=\; \lambda_k ,

and the sample variance of the $k$ th principal component is $\lambda_k$ . Exactly $r$ of the $\lambda_k$ are non-zero. The choice is unique up to sign whenever $\lambda_{k-1} > \lambda_k > \lambda_{k+1}$ , read with the end conventions $\lambda_0 = +\infty$ and $\lambda_{p+1} = -\infty$ (so at $k = 1$ the condition is $\lambda_1 > \lambda_2$ and at $k = p$ it is $\lambda_{p-1} > \lambda_p$ ), and within any block of equal eigenvalues no direction is singled out at all.

2. The principal components are uncorrelated, and the variances add up. For $j \ne k$ the sample covariance of $z_j = Xe_j$ and $z_k = Xe_k$ is $z_j^{\mathsf T}z_k/(n-1) = e_j^{\mathsf T}Se_k = \lambda_k\,e_j\cdot e_k = 0$ , and

\sum_{k=1}^{p}\lambda_k \;=\; \operatorname{tr}S ,

so, when $X \ne 0$ and hence $\operatorname{tr}S > 0$ , the explained variance fractions $\lambda_k/\operatorname{tr}S$ are non-negative and sum to one.

3. They come straight out of the SVD of $X$ . Let $X = U\Sigma V^{\mathsf T}$ be a singular value decomposition (Singular value decomposition). An $n \times p$ matrix has $q = \min(n,p)$ singular values $\sigma_1 \ge \dots \ge \sigma_q \ge 0$ ; extend the list by $\sigma_k = 0$ for $q < k \le p$ , which is consistent because $S$ then has $p - q$ further zero eigenvalues. With that convention one may take, for $k = 1,\dots,p$ ,

e_k = v_k \quad(\text{the }k\text{th column of }V), \qquad \lambda_k = \frac{\sigma_k^{2}}{n-1}, \qquad z_k = Xv_k = \begin{cases} \sigma_k u_k, & k \le q,\\ 0, & k > q . \end{cases}

The second case is empty when $p \le n$ ; when $p > n$ — the case the remark below stresses — $U$ supplies only $u_1,\dots,u_n$ , and the remaining principal components vanish. In particular the principal directions never require $S$ to be formed.

4. Truncation is optimal. For $1 \le m < r$ let

\widehat X^{(m)} \;=\; \sum_{k=1}^{m} z_k e_k^{\mathsf T} \;=\; \sum_{k=1}^{m}\sigma_k u_kv_k^{\mathsf T}

be the field reconstructed from the first $m$ modes. Then $\widehat X^{(m)}$ minimises $\lVert X - B\rVert_F$ over all $n \times p$ matrices $B$ of rank at most $m$ , the minimum value being $\bigl(\sum_{k>m}\sigma_k^{2}\bigr)^{1/2}$ , and the fraction of the total variance it retains is

\frac{\sum_{k\le m}\lambda_k}{\operatorname{tr}S} \;=\; \frac{\sum_{k\le m}\sigma_k^{2}}{\sum_{k\le p}\sigma_k^{2}} .

Equivalently, among all $m$ -dimensional subspaces of $\mathbb{R}^p$ , the span of $e_1,\dots,e_m$ is one onto which the orthogonal projection of the data loses the least squared length.

Idea of the proofOriginal

Two short observations do all the work.

The first is that $a^{\mathsf T}Sa$ depends on $a$ only through its components in an eigenbasis: writing $a = \sum_k c_ke_k$ with $\sum_k c_k^2 = 1$ gives $a^{\mathsf T}Sa = \sum_k \lambda_kc_k^2$ , a weighted average of the eigenvalues. An average of numbers cannot exceed the largest of them, and it reaches that value only by putting all the weight on entries equal to the largest — on one eigenvector when the largest eigenvalue is simple, and anywhere in its eigenspace when it is not. So maximising a quadratic form over the unit sphere is not a calculus problem at all once the spectral theorem has been applied; it is the remark that a convex combination is largest at an extreme point. Imposing orthogonality to $e_1,\dots,e_{k-1}$ deletes the first $k-1$ weights and repeats the argument on what is left.

The second is that $S$ and $X$ carry the same information: $X^{\mathsf T}X = (n-1)S$ , so the right singular vectors of $X$ are eigenvectors of $S$ and the squared singular values are $(n-1)$ times its eigenvalues. Truncating the principal component expansion is then literally truncating the dyadic sum $\sum_k\sigma_ku_kv_k^{\mathsf T}$ , and its optimality is the Eckart-Young theorem already proved in Singular value decomposition.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §8.17.1, Eqs. (8.113)-(8.116), pp. 290-291; Riley, Hobson & Bence (2006), §8.18.3, Eqs. (8.131)-(8.136), pp. 302-304

Part 1. $S$ is real symmetric, so by the spectral theorem (Spectral theorem for Hermitian matrices) it has an orthonormal eigenbasis $e_1,\dots,e_p$ with real eigenvalues, which we order decreasingly; they are non-negative because $\lambda_k = e_k^{\mathsf T}Se_k = \lVert Xe_k\rVert^2/(n-1) \ge 0$ .

Fix $k$ and let $a$ be a unit vector orthogonal to $e_1,\dots,e_{k-1}$ . Expand $a = \sum_{i=1}^{p}c_ie_i$ with $c_i = a\cdot e_i$ ; orthogonality to the first $k-1$ eigenvectors says $c_1 = \dots = c_{k-1} = 0$ , and $\lVert a \rVert = 1$ says $\sum_{i\ge k}c_i^2 = 1$ . Then, using $Se_i = \lambda_ie_i$ and $e_i\cdot e_j = \delta_{ij}$ ,

a^{\mathsf T}Sa \;=\; \sum_{i\ge k}\sum_{j \ge k} c_ic_j\,e_i^{\mathsf T}Se_j \;=\; \sum_{i\ge k}\lambda_ic_i^{2} \;\le\; \lambda_k\sum_{i\ge k}c_i^{2} \;=\; \lambda_k ,

the inequality because $\lambda_i \le \lambda_k$ for every $i \ge k$ , by the decreasing ordering, and because every $c_i^2 \ge 0$ . The value $\lambda_k$ is attained at $a = e_k$ , which is a unit vector orthogonal to $e_1,\dots,e_{k-1}$ . So the maximum is $\lambda_k$ and $e_k$ achieves it, which is what the defining rule for $e_k$ asks. The variance of $z_k = Xe_k$ is $\lVert Xe_k\rVert^2/(n-1) = e_k^{\mathsf T}Se_k = \lambda_k$ .

The number of non-zero eigenvalues is $\operatorname{rank}S$ , and $\operatorname{rank}S = \operatorname{rank}(X^{\mathsf T}X) = \operatorname{rank}X = r$ : indeed $Xa = 0$ implies $X^{\mathsf T}Xa = 0$ , and conversely $X^{\mathsf T}Xa = 0$ implies $\lVert Xa\rVert^2 = a^{\mathsf T}X^{\mathsf T}Xa = 0$ , so the two matrices have the same null space.

For uniqueness, suppose $\lambda_{k-1} > \lambda_k > \lambda_{k+1}$ (read $\lambda_0 = +\infty$ at the left end and $\lambda_{p+1} = -\infty$ at the right, so that for $k = p$ the hypothesis is only $\lambda_{p-1} > \lambda_p$ ) and let $a$ be any maximiser. In the display above, equality in $\sum_{i\ge k}\lambda_ic_i^2 \le \lambda_k\sum_{i\ge k}c_i^2$ forces $c_i = 0$ for every $i > k$ with $\lambda_i < \lambda_k$ , which under the assumed strict gap is every $i>k$ ; hence $a = c_ke_k$ with $|c_k| = 1$ . If instead $\lambda_k = \lambda_{k+1}$ , every unit vector in the corresponding eigenspace is a maximiser.

That argument used only the right-hand gap $\lambda_k > \lambda_{k+1}$ , and it determines $e_k$ up to sign given $e_1,\dots,e_{k-1}$ . The left-hand gap $\lambda_{k-1} > \lambda_k$ is what makes the conclusion a statement about the rule and not about one run of it: the constraint set $\{a : \lVert a\rVert = 1,\ a \perp e_1,\dots,e_{k-1}\}$ depends on the earlier choices, and those are themselves arbitrary inside any block of equal eigenvalues. If $\lambda_{k-1} = \lambda_k$ , then $e_{k-1}$ may be taken to be any unit vector of the shared eigenspace, and $e_k$ is then forced to be a vector of that same eigenspace orthogonal to it — so $e_k$ moves with $e_{k-1}$ and is not determined up to sign. With both gaps, $\lambda_k$ is a simple eigenvalue, and the left gap also pins down the constraint set: any admissible $e_1,\dots,e_{k-1}$ span the sum of the eigenspaces of $\lambda_1,\dots,\lambda_{k-1}$ , and since $\lambda_{k-1} > \lambda_k$ that sum contains every block of equal eigenvalues it meets, so it is the same subspace — the sum of all eigenspaces for eigenvalues exceeding $\lambda_k$ — however the ties among the earlier modes were resolved. Hence $e_k$ is the same line whatever admissible $e_1,\dots,e_{k-1}$ preceded it, which is what "unique up to sign" asserts.

Part 2. The $z_k$ have mean zero because the columns of $X$ do, so their sample covariance is $z_j^{\mathsf T}z_k/(n-1) = e_j^{\mathsf T}X^{\mathsf T}Xe_k/(n-1) = e_j^{\mathsf T}Se_k = \lambda_k\,e_j\cdot e_k$ , which vanishes for $j \ne k$ by orthonormality. The trace identity is the fact that a square matrix has trace equal to the sum of its eigenvalues, applied to $S$ ; concretely, $\operatorname{tr}S = \operatorname{tr}(Q\Lambda Q^{\mathsf T}) = \operatorname{tr}(\Lambda Q^{\mathsf T}Q) = \operatorname{tr}\Lambda$ with $Q = (e_1 \cdots e_p)$ orthogonal.

Part 3. Let $X = U\Sigma V^{\mathsf T}$ be a singular value decomposition. Then

X^{\mathsf T}X \;=\; V\Sigma^{\mathsf T}U^{\mathsf T}U\Sigma V^{\mathsf T} \;=\; V(\Sigma^{\mathsf T}\Sigma)V^{\mathsf T} ,

using $U^{\mathsf T}U = I$ . This is the computation Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (8.132), p. 302) make at their Eq. (8.132); since $\Sigma^{\mathsf T}\Sigma$ is the $p \times p$ diagonal matrix with entries $\sigma_1^2,\dots,\sigma_q^2$ followed by $p - q$ zeros — that is, with entries $\sigma_k^2$ under the convention $\sigma_k = 0$ for $k > q$ — and $V$ is orthogonal, it exhibits an orthogonal diagonalisation of $X^{\mathsf T}X = (n-1)S$ . Hence the columns $v_k$ of $V$ are orthonormal eigenvectors of $S$ with eigenvalues $\sigma_k^2/(n-1)$ , and because the SVD convention orders $\sigma_1 \ge \sigma_2 \ge \dots$ , these eigenvalues come out already in decreasing order. By part 1 they are an admissible set of principal directions and variances. Finally $Xv_k = \sigma_ku_k$ for $k \le q$ , which is their Eq. (8.134), so $z_k = \sigma_ku_k$ there; for $k > q$ the $k$ th column of $\Sigma$ is zero, so $z_k = Xv_k = U\Sigma V^{\mathsf T}v_k = 0$ , in agreement with $\lambda_k = 0$ .

Part 4. Substituting $z_k = \sigma_ku_k$ and $e_k = v_k$ into the definition of $\widehat X^{(m)}$ gives $\sum_{k\le m}\sigma_ku_kv_k^{\mathsf T}$ , which is the rank- $m$ truncation of the dyadic form of $X$ (Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), Eq. (8.136), p. 303) Eq. (8.136)). Its optimality in the Frobenius norm, with the stated minimum, is the Eckart-Young statement of Singular value decomposition. Dividing $\lVert X\rVert_F^2 = \sum_k\sigma_k^2$ by $n-1$ turns it into $\operatorname{tr}S = \sum_k\lambda_k$ , which gives the variance fraction. For the last sentence, let $P$ be the orthogonal projection of $\mathbb{R}^p$ onto an $m$ -dimensional subspace; then $XP$ has rank at most $m$ , so $\lVert X - XP\rVert_F \ge \lVert X - \widehat X^{(m)}\rVert_F$ , and $\widehat X^{(m)} = XP_m$ with $P_m$ the projection onto $\operatorname{span}(e_1,\dots,e_m)$ , because $\sum_{k\le m}Xe_ke_k^{\mathsf T} = X\sum_{k\le m}e_ke_k^{\mathsf T}$ and the latter sum is exactly $P_m$ (Orthogonal projection and least squares). $\square$

Remark(Where the hypotheses are used)Original

Centring is what makes $S$ a covariance matrix, and it is the hypothesis most often dropped by accident. Parts 3 and 4 — the identification with the SVD and the Eckart-Young optimality — are facts about an arbitrary matrix $X$ and hold whether or not it is centred; what centring buys is the reading of $a^{\mathsf T}Sa$ as a sample variance, and with it the zero mean of the principal components used in part 2. If the columns of $X$ are not centred, then $X^{\mathsf T}X/(n-1)$ is a matrix of second moments, not of covariances; its leading eigenvector is dominated by the mean $\bar x$ and typically just reproduces it. Concretely, take a single variable with the two values $10$ and $12$ , so that $n = 2$ and $n - 1 = 1$ . Centred, the formula gives the sample variance $\bigl((-1)^2 + 1^2\bigr)/1 = 2$ ; applied to the raw column it gives $(10^2 + 12^2)/1 = 244$ instead, of which the mean alone contributes $2\cdot 11^2/1 = 242$ . So more than $99\%$ of the quantity being maximised has nothing to do with variability, and the "first mode" of an uncentred analysis is just the mean.
The decreasing ordering is used once, in the inequality $\lambda_i \le \lambda_k$ for $i \ge k$ . Without it the argument still shows that stationary directions are eigenvectors — this is the Lagrange-multiplier computation of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §8.17.1, Eqs. (8.113)-(8.116), pp. 290-291), who note that stationarity is necessary but not sufficient for a maximum — but it no longer identifies which eigenvector is the maximiser. The convex-combination argument above settles that and is why we do not need the multiplier.
Symmetry of $S$ is what gives an orthonormal eigenbasis, and with it both the uncorrelatedness in part 2 and the completeness used in the expansion $a = \sum c_ie_i$ . It is automatic here, because $S$ is built as $X^{\mathsf T}X$ ; there is no hypothesis to check and no non-symmetric variant of the method.
Degenerate eigenvalues destroy the modes, not the subspace. When $\lambda_k = \lambda_{k+1}$ , the individual patterns $e_k$ and $e_{k+1}$ are arbitrary inside a two-dimensional eigenspace and will change completely under a tiny perturbation of the data; only their span, and the sum of the two explained variances, is stable. Near-degeneracy is the usual situation for a propagating signal, whose two phases have nearly equal variance, and it is the standard reason an individual EOF should not be read as a physical mode.
Rank is limited by the smaller dimension. Since $r \le \min(n,p)$ , a field on $p = 10^4$ grid points observed on $n = 50$ occasions has at most $50$ non-zero eigenvalues, and at most $49$ after centring, because centring puts the rows in the hyperplane $\sum_t x_t = 0$ . Part 3 is then the only practical route: the $p \times p$ matrix $S$ need never be built.
Scaling is a modelling choice, not mathematics. $S$ is not invariant under rescaling the individual variables, so for variables in different units — hectopascals and metres per second — PCA of $S$ is meaningless and one works with the correlation matrix instead, that is, with each column divided by its sample standard deviation. For a field in one unit this choice is usually inappropriate, but a weight proportional to the square root of the area each grid point represents is needed, or dense high-latitude grid points are over-counted.
$n-1$ versus $n$ changes every $\lambda_k$ by the same factor and therefore changes no eigenvector and no explained variance fraction. We use $n-1$ for consistency with the unbiased sample variance of Sampling distributions.
Row vectors. Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.9, pp. 132-133) write random vectors as rows, so their covariance matrix is $\mathbb{E}[(X-\mu)'(X-\mu)]$ and their linear map acts on the right, $Y = XD$ with covariance $D'VD$ . We use column vectors, in which the same statements read $V = \mathbb{E}[(X-\mu)(X-\mu)^{\mathsf T}]$ and $\operatorname{cov}(D^{\mathsf T}X) = D^{\mathsf T}VD$ . The transposes are not a disagreement, only a transcription.

Method(Computing an EOF analysis)Original

Lay out the data matrix. Rows are observations (times), columns are variables (grid points). Decide first what the "observations" are: EOFs of a field in time and EOFs of the same field transposed are different analyses, related by part 3.
Handle missing values and weights before anything else. Fill or drop; then, for a field on a sphere, multiply column $j$ by $\sqrt{w_j}$ with $w_j$ the area weight of grid point $j$ , so that $\operatorname{tr}S$ is the area-weighted variance. For variables in different units, divide each column by its sample standard deviation instead.
Centre. Subtract each column's mean. If a trend or a seasonal cycle is to be excluded rather than analysed, remove it here too; otherwise it will be mode 1.
Factor, do not square. Compute the singular value decomposition $X = U\Sigma V^{\mathsf T}$ (Singular value decomposition). Forming $S = X^{\mathsf T}X/(n-1)$ and diagonalising it gives the same answer in exact arithmetic but squares the condition number (Matrix norms and condition numbers), and is impossible anyway when $p$ is large.
Read off the modes. EOF $k$ is the column $v_k$ , plotted as a map; the PC time series is $z_k = \sigma_ku_k$ ; the variance of mode $k$ is $\lambda_k = \sigma_k^2/(n-1)$ and its explained fraction is $\sigma_k^2/\sum_i\sigma_i^2$ .
Fix the sign. Each pair $(v_k, z_k)$ may be replaced by $(-v_k, -z_k)$ ; software makes an arbitrary choice. Fix it by a stated rule, for example by requiring the value of $v_k$ at a chosen grid point, or the projection of $z_k$ on a chosen index, to be positive, and say which rule you used.
Decide how many modes to keep. Plot $\lambda_k$ against $k$ and truncate where the spectrum flattens; check whether the gap you are truncating at is larger than the sampling uncertainty in the eigenvalues, which for $n$ independent observations is of order $\lambda_k\sqrt{2/n}$ . Keep whole near-degenerate pairs or none of them.
Reconstruct and check. Form $\widehat X^{(m)} = \sum_{k\le m}z_kv_k^{\mathsf T}$ and verify that $\lVert X - \widehat X^{(m)}\rVert_F^2$ equals $\sum_{k>m}\sigma_k^2$ . This is an arithmetic check on the whole computation and costs nothing.

Example(A field at four grid points, worked through exactly)Original

A one-dimensional "field" is measured at $p = 4$ grid points on $n = 5$ occasions. The data, already centred (each column sums to zero), are

X \;=\; \begin{pmatrix} 4 & 2 & -4 & -2\\ 0 & 2 & 0 & -2\\ -2 & 0 & 2 & 0\\ -2 & -4 & 2 & 4\\ 0 & 0 & 0 & 0 \end{pmatrix} .

Step 1: the sample covariance matrix. With $n - 1 = 4$ ,

S \;=\; \tfrac14 X^{\mathsf T}X \;=\; \begin{pmatrix} 6 & 4 & -6 & -4\\ 4 & 6 & -4 & -6\\ -6 & -4 & 6 & 4\\ -4 & -6 & 4 & 6 \end{pmatrix}, \qquad \operatorname{tr}S = 24 .

For instance the $(1,2)$ entry is $\tfrac14\bigl(4\cdot2 + 0\cdot2 + (-2)\cdot0 + (-2)\cdot(-4) + 0\bigr) = \tfrac14(8+0+0+8) = 4$ .

Step 2: guess and verify two eigenvectors. Rows $1$ and $3$ of $S$ are negatives of each other, and so are rows $2$ and $4$ ; the matrix therefore has rank at most $2$ , and the two sign patterns $f = (1,1,-1,-1)^{\mathsf T}$ and $g = (1,-1,-1,1)^{\mathsf T}$ suggest themselves. Multiplying out,

Sf = (6+4+6+4,\; 4+6+4+6,\; -6-4-6-4,\; -4-6-4-6)^{\mathsf T} = 20\,f ,

Sg = (6-4+6-4,\; 4-6+4-6,\; -6+4-6+4,\; -4+6-4+6)^{\mathsf T} = 4\,g .

So $\lambda_1 = 20$ and $\lambda_2 = 4$ , with $f \cdot g = 1 - 1 + 1 - 1 = 0$ . Since $\lambda_1 + \lambda_2 = 24 = \operatorname{tr}S$ and the remaining eigenvalues are non-negative, they must both be zero: the data lie exactly in a two-dimensional subspace, and $r = 2$ .

Step 3: normalise, and compute the principal components. $\lVert f\rVert = \lVert g\rVert = 2$ , so the EOFs are

e_1 = \tfrac12(1,1,-1,-1)^{\mathsf T}, \qquad e_2 = \tfrac12(1,-1,-1,1)^{\mathsf T} .

EOF 1 is a dipole: the left half of the domain up, the right half down. EOF 2 opposes the two outer grid points to the two inner ones — its signs across grid points $1,2,3,4$ are $+,-,-,+$ , so it changes sign between points $1$ and $2$ and between points $3$ and $4$ , but not between $2$ and $3$ . It is not the grid-scale pattern that alternates at every step: that vector, $(1,-1,1,-1)^{\mathsf T}$ , satisfies $S(1,-1,1,-1)^{\mathsf T} = 0$ here, so it carries no sample variance at all and spans one of the two null directions whose zero eigenvalues step 2 accounted for (the other is the uniform pattern $(1,1,1,1)^{\mathsf T}$ , also annihilated by $S$ ). The PC time series are $z_k = Xe_k$ :

z_1 = (6,\,2,\,-2,\,-6,\,0)^{\mathsf T}, \qquad z_2 = (2,\,-2,\,-2,\,2,\,0)^{\mathsf T} .

(First entry of $z_1$ : $\tfrac12(4 + 2 + 4 + 2) = 6$ .) Their sample variances are $\lVert z_1\rVert^2/4 = (36+4+4+36)/4 = 20 = \lambda_1$ and $\lVert z_2\rVert^2/4 = 16/4 = 4 = \lambda_2$ , as part 1 promises, and $z_1 \cdot z_2 = 12 - 4 + 4 - 12 + 0 = 0$ , as part 2 promises.

Step 4: explained variance.

\frac{\lambda_1}{\operatorname{tr}S} = \frac{20}{24} = \frac56 = 83.33\%, \qquad \frac{\lambda_2}{\operatorname{tr}S} = \frac{4}{24} = \frac16 = 16.67\% .

The dipole alone accounts for five sixths of the variance of the field; the two modes together account for all of it.

Step 5: cross-check against the SVD. The singular values must be $\sigma_k = \sqrt{(n-1)\lambda_k}$ , that is $\sigma_1 = \sqrt{80} = 4\sqrt5 \approx 8.9443$ and $\sigma_2 = \sqrt{16} = 4$ , with $\sigma_3 = \sigma_4 = 0$ ; and indeed $\lVert X\rVert_F^2 = \sum_{t,j}X_{tj}^2 = 96 = 80 + 16$ . Truncating after one mode leaves $\lVert X - \widehat X^{(1)}\rVert_F^2 = \sigma_2^2 = 16$ , and no rank-one matrix does better.

ApplicationOriginal

Climate fields. Applied to monthly anomalies of sea-level pressure over the northern hemisphere, the leading EOF is the annular seesaw between the polar cap and the midlatitude belt whose PC time series is published as the Arctic Oscillation, or northern annular mode, index; the same analysis restricted to the Atlantic sector gives the North Atlantic Oscillation, the seesaw between the Icelandic low and the Azores high; applied to tropical Pacific sea-surface temperature, the leading EOF is the tongue of warm anomaly whose PC is an El Nino index. The attraction is exactly part 4: a handful of numbers per month reproduces a field of tens of thousands of grid points to within a few per cent of its variance, so a century of maps compresses to a few time series that can be correlated with everything else. The attendant danger is the degeneracy caveat above — the leading EOF is a statistically optimal pattern, and only sometimes also a dynamical mode.

Model reduction and data assimilation. Truncating at $m$ modes replaces a state vector of length $p$ by one of length $m$ in the subspace that carries the most variance, which is how a reduced-rank background covariance is built for Optimal interpolation and variational data assimilation, and how the "essential dynamics" of a protein is extracted from a molecular-dynamics trajectory: the atomic coordinate matrix has $p = 3N_{\text{atoms}}$ columns, and a handful of collective modes carries most of the fluctuation.

Measurement and instrumentation. The same decomposition applied to a set of spectra separates the number of independent chemical species present (the number of eigenvalues above the noise floor) from the noise, and applied to an array of detectors it identifies common-mode drift as a leading mode to be regressed out. In both cases what is being used is that the noise, being uncorrelated between channels, contributes a nearly flat tail to the eigenvalue spectrum, while any shared signal contributes a few large eigenvalues.

Method

The periodogram and spectral estimation

Adapted

Add to PDF

Defines the periodogram of a finite record as the squared modulus of its discrete Fourier transform, shows that it is exactly the transform of the sample autocovariance and that its expectation is the true spectral density smoothed by the Fejer kernel of the data window, computes its distribution for Gaussian white noise to show that it never becomes consistent however long the record, and sets out the two repairs — tapering against leakage and segment averaging against variance.

Prerequisites: The discrete Fourier transform and the FFT, Stationary processes and spectral densities, Cesàro means and Fejér's theorem, The multivariate normal distribution, Expectation, variance and covariance, Sampling distributions

Note(Conventions, and two densities that differ by a factor)Original

Stationary processes and spectral densities defines the spectral density of a discrete-time stationary process through its autocorrelation $\rho(n) = c(n)/c(0)$ , so that the density

f(\lambda) = \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}\rho(n)\,\mathrm{e}^{-in\lambda}, \qquad \lambda\in(-\pi,\pi],

is a probability density: it integrates to one. In data analysis one wants the spectrum in the units of the data, so throughout this unit we use the power spectral density

g(\lambda) \;=\; c(0)\,f(\lambda) \;=\; \frac{1}{2\pi}\sum_{n\in\mathbb{Z}}c(n)\,\mathrm{e}^{-in\lambda}, \qquad\text{so that}\qquad c(n) = \int_{-\pi}^{\pi}\mathrm{e}^{in\lambda}g(\lambda)\,d\lambda \quad\text{and}\quad \int_{-\pi}^{\pi}g = c(0) = \operatorname{var}X_t .

The two differ only by the constant factor $c(0)$ , and the summability hypothesis is the same one, Grimmett and Stirzaker's (Grimmett & Stirzaker (2020), (9.3.18), p. 442) condition $\sum_n\lvert\rho(n)\rvert<\infty$ . Every formula below can be read in either normalisation by inserting or removing one factor of $\operatorname{var}X_t$ ; mixing them silently is the commonest unit error in this subject.

Time is measured in sampling intervals, so the frequency variable $\lambda$ is in radians per sample and the Nyquist frequency is $\lambda = \pi$ . For data sampled every $\Delta$ seconds, the physical angular frequency is $\omega = \lambda/\Delta$ and the density per unit $\omega$ is $\Delta\,g(\omega\Delta)$ ; this is the factor Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Exercise 13.10(a), p. 462) carry as $\tau$ in their Exercise 13.10.

Definition(Fourier frequencies, the periodogram, tapers)AdaptedAdapted from Arfken, Weber & Harris (2013), §19.1, Exercise 19.1.2, p. 945; Riley, Hobson & Bence (2006), §13.4, Exercise 13.10, p. 462

Let $x_0, x_1, \dots, x_{N-1}$ be a finite record. The Fourier frequencies are

\lambda_j = \frac{2\pi j}{N}, \qquad j = 0, 1, \dots, N-1 ,

the frequencies for which the record contains a whole number of cycles. Values of $j$ above $N/2$ are read as the negative frequencies $\lambda_j - 2\pi$ , so the $\lambda_j$ fill $(-\pi,\pi]$ .

The discrete Fourier sum of the record is

J_N(\lambda) \;=\; \sum_{t=0}^{N-1} x_t\,\mathrm{e}^{-i\lambda t} ,

which at $\lambda = \lambda_j$ is $N$ times the discrete Fourier transform $\hat x(j)$ of The discrete Fourier transform and the FFT, whose convention carries the factor $1/N$ in front. The periodogram is

I_N(\lambda) \;=\; \frac{1}{2\pi N}\,\bigl\lvert J_N(\lambda)\bigr\rvert^{2} .

Squaring discards the phase: this is the property Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 19.1.2, p. 945) single out when they call $a_n^2+b_n^2$ the power spectrum and note its invariance under a phase shift. The $2\pi$ and the $N$ are chosen so that $I_N$ estimates $g$ and not a multiple of it; several other normalisations are in use.

A taper (or data window) is a sequence of weights $w_0,\dots,w_{N-1}$ , not all zero, applied before transforming. The tapered periodogram is

I^{w}_N(\lambda) \;=\; \frac{1}{2\pi\sum_{t}w_t^{2}}\;\Bigl\lvert \sum_{t=0}^{N-1} w_t\,x_t\,\mathrm{e}^{-i\lambda t}\Bigr\rvert^{2} ,

the denominator being what keeps the scale right; $w_t \equiv 1$ gives $I_N$ back. Its spectral window is

\Phi^{w}_N(\theta) \;=\; \frac{1}{2\pi\sum_{t}w_t^{2}}\;\Bigl\lvert \sum_{t=0}^{N-1} w_t\,\mathrm{e}^{-i\theta t}\Bigr\rvert^{2} ,

a non-negative even function of $\theta$ with $\int_{-\pi}^{\pi}\Phi^{w}_N = 1$ . For $w_t\equiv1$ the geometric sum gives the closed form

\Phi_N(\theta) \;=\; \frac{1}{2\pi N}\,\frac{\sin^{2}(N\theta/2)}{\sin^{2}(\theta/2)} ,

which is the Fejer kernel of Cesàro means and Fejér's theorem in the normalisation Stakgold and Holst (Stakgold & Holst (2011), Exercise 2.3.6, p. 144) use.

Finally, the sample autocovariance of the record at lag $h$ , $\lvert h\rvert \le N-1$ , is

\hat c(h) \;=\; \frac{1}{N}\sum_{t=0}^{N-1-\lvert h\rvert}(x_t - \bar x)\,(x_{t+\lvert h\rvert} - \bar x), \qquad \bar x = \frac1N\sum_t x_t .

Note the divisor $N$ , not $N - \lvert h\rvert$ : the reason is part 1 below.

Method(What the periodogram is, and what it is not)Original

Let $X_0,\dots,X_{N-1}$ be a stretch of a real weakly stationary process with autocovariance $c$ and power spectral density $g$ , and let $I_N$ , $I^w_N$ , $\Phi^w_N$ and $\hat c$ be as defined above.

1. It is exactly the transform of the sample autocovariance, and it carries all of the mean square of the record. At every Fourier frequency $\lambda_j$ with $j \ne 0$ ,

I_N(\lambda_j) \;=\; \frac{1}{2\pi}\sum_{h=-(N-1)}^{N-1}\hat c(h)\,\mathrm{e}^{-ih\lambda_j} ,

and summing over the grid,

\frac{2\pi}{N}\sum_{j=0}^{N-1} I_N(\lambda_j) \;=\; \frac1N\sum_{t=0}^{N-1}x_t^{2} .

So the periodogram is the finite-sample copy of the inversion formula $g(\lambda) = (2\pi)^{-1}\sum_n c(n)\mathrm{e}^{-in\lambda}$ , and it partitions the mean square of the record over the frequency grid. The divisor $N$ in $\hat c$ is what makes the first identity exact; with $N - \lvert h\rvert$ it fails, and the right-hand side need not even be non-negative.

2. Its expectation is the true density smoothed by the spectral window. If $\sum_n\lvert c(n)\rvert < \infty$ and $\mathbb{E}X_t = 0$ , then for every $\lambda$

\mathbb{E}\,I^{w}_N(\lambda) \;=\; \int_{-\pi}^{\pi}\Phi^{w}_N(\lambda - \mu)\,g(\mu)\,d\mu ,

and in the untapered case this is

\mathbb{E}\,I_N(\lambda) \;=\; \frac{1}{2\pi}\sum_{\lvert h\rvert < N}\Bigl(1 - \frac{\lvert h\rvert}{N}\Bigr)c(h)\,\mathrm{e}^{-ih\lambda} \;\xrightarrow[N\to\infty]{}\; g(\lambda) ,

the convergence being uniform in $\lambda$ . The periodogram is therefore asymptotically unbiased but biased at any finite $N$ , and the bias has two distinguishable parts: the main lobe of $\Phi^w_N$ blurs sharp features (limited resolution), and the sidelobes of $\Phi^w_N$ transport power from strong frequencies to distant weak ones (leakage). Among tapers that are heaviest at the centre the main lobe is narrowest with no taper at all: for $w_t\equiv1$ the full width of the main lobe of $\Phi_N$ between its first zeros is exactly $4\pi/N$ , and every taper that lowers the sidelobes widens it — to $8\pi/N$ for the Hann taper of the method below. Narrow main lobe and low sidelobes pull against each other: for $w_t\equiv1$ the first sidelobe of $\Phi_N$ sits about $1.43$ frequency bins from the centre at about $4.7\%$ of the peak, that is $-13.3$ decibels, and a taper buys a lower figure than that with a wider main lobe.

3. It is not consistent. Let $X_t$ be independent $N(0,\sigma^2)$ , so that $g \equiv \sigma^2/2\pi$ . Then for each Fourier frequency $\lambda_j$ with $0 < j < N/2$ ,

\frac{I_N(\lambda_j)}{g} \;\sim\; \tfrac12\chi^{2}_{2} , \qquad\text{equivalently}\qquad I_N(\lambda_j) \sim \text{Exponential with mean } g ,

and these variables are independent for different such $j$ ; at the two self-conjugate frequencies $\lambda = 0$ and (for even $N$ ) $\lambda = \pi$ one gets $I_N/g \sim \chi^2_1$ instead. Consequently

\mathbb{E}\,I_N(\lambda_j) = g, \qquad \operatorname{var}I_N(\lambda_j) = g^{2}, \qquad \frac{\text{standard deviation}}{\text{mean}} = 1 \quad\text{for every }N .

Lengthening the record adds more frequencies; it does not improve any of them. A raw periodogram of a long record is a dense forest of spikes whose individual heights are no better determined than those from a short one.

4. Averaging is the repair. Split the record into $K$ segments of length $M = N/K$ , form the tapered periodogram of each, and average:

\bar I(\lambda) \;=\; \frac1K\sum_{k=1}^{K} I^{w,(k)}_M(\lambda) .

For Gaussian white noise and non-overlapping segments the $K$ terms are independent, and each has mean exactly $g$ whatever the taper, so $\mathbb{E}\,\bar I = g$ . The law, however, depends on the taper. Write the condition

\textbf{(T)}\qquad \sum_{t=0}^{M-1} w_t^{2}\,\mathrm{e}^{-2i\lambda_j t} \;=\; 0 .

At an interior Fourier frequency $\lambda_j$ of the segment grid, $0 < j < M/2$ , at which (T) holds,

\frac{\bar I(\lambda_j)}{g} \;\sim\; \frac{\chi^{2}_{2K}}{2K}, \qquad \mathbb{E}\,\bar I = g, \qquad \frac{\text{standard deviation}}{\text{mean}} = \frac{1}{\sqrt{K}} .

For $w_t\equiv1$ condition (T) holds at every interior $\lambda_j$ , so the untapered statement is exact across the grid. For a general taper it can fail, and where it fails the mean is still exactly $g$ but the chi-squared law and the figure $K^{-1/2}$ are only approximations. The failures are few and sit at the two ends of the band: for the Hann taper of the method below, with $M$ even, the only interior exceptions are $j = 1$ and $j = M/2-1$ , where the relative standard deviation is $\sqrt{37}/(6\sqrt{K}\,) = 1.0138\,K^{-1/2}$ instead of $K^{-1/2}$ . The price is resolution: the frequency grid is now spaced $2\pi/M$ instead of $2\pi/N$ , and the main lobe is $K$ times wider. This trade — variance down by $K$ , resolution down by $K$ — is the whole of the design problem. Welch's method is this procedure with a taper on each segment and the segments overlapped by half their length, which recovers part of the data that the tapers suppress at the cost of making the segments no longer independent.

Idea of the proofOriginal

The periodogram is a squared Fourier sum, so everything about it follows from writing the square as a double sum over pairs of times and collecting terms by their separation.

Do that to the data and you get the sample autocovariance: that is part 1, and it is an identity, with no probability in it. Do the same to the expectation and the pair $(s,t)$ contributes $c(s-t)$ , which depends only on $h = s-t$ by stationarity; there are $N - \lvert h\rvert$ pairs at separation $h$ , so the triangular factor $1 - \lvert h\rvert/N$ appears, and that factor is precisely what converts the exact inversion sum into a convolution with the Fejer kernel. The record length enters only through the width of that kernel, which is why the bias disappears as $N$ grows.

The variance does not, and the reason is a counting argument. A record of length $N$ gives $N$ numbers, the periodogram evaluates them at about $N/2$ frequencies, and each value is built out of just two of them — the sine and cosine amplitudes at that frequency. The amount of data behind each estimate is two numbers however long the record is, and two numbers cannot have small relative error. That is part 3, and it also tells us the only possible cure: pool neighbouring or repeated estimates until the count is large. Part 4 pools repetitions in time; smoothing the periodogram across adjacent frequencies pools neighbours and does the same arithmetic.

Proof(Parts 1-4)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.6, Eq. (20.120) and Exercises 20.6.1 and 20.6.3(a), pp. 1004-1007; Grimmett & Stirzaker (2020), §9.3, the discrete-time spectral representation and the inversion formula (9.3.16), pp. 441-442

Part 1. Fix $j \ne 0$ . Then $\sum_{t=0}^{N-1}\mathrm{e}^{-i\lambda_j t} = 0$ by the discrete orthogonality relation of The discrete Fourier transform and the FFT, so replacing $x_t$ by $y_t = x_t - \bar x$ changes $J_N(\lambda_j)$ by $-\bar x\sum_t\mathrm{e}^{-i\lambda_j t} = 0$ , that is, not at all. Expanding the modulus squared as a double sum,

\bigl\lvert J_N(\lambda_j)\bigr\rvert^{2} = \sum_{s=0}^{N-1}\sum_{t=0}^{N-1} y_sy_t\,\mathrm{e}^{-i(s-t)\lambda_j} .

Group the pairs by $h = s - t$ , which runs over $-(N-1),\dots,N-1$ ; for fixed $h \ge 0$ the admissible $t$ are $0,\dots,N-1-h$ , and the inner sum is $\sum_t y_{t+h}y_t = N\hat c(h)$ by the definition of $\hat c$ , with the divisor $N$ cancelling. For $h<0$ the same with $\lvert h\rvert$ , since $\hat c(-h) = \hat c(h)$ . Hence $\lvert J_N(\lambda_j)\rvert^2 = N\sum_{\lvert h\rvert<N}\hat c(h)\mathrm{e}^{-ih\lambda_j}$ , and dividing by $2\pi N$ gives the identity. Had we divided by $N - \lvert h\rvert$ in the definition of $\hat c$ , the cancellation would not occur and each lag would be reweighted by $N/(N-\lvert h\rvert)$ .

For the sum rule, the discrete Parseval identity of The discrete Fourier transform and the FFT, written for $J_N(\lambda_j) = N\hat x(j)$ , reads $\sum_{j}\lvert J_N(\lambda_j)\rvert^2 = N\sum_t\lvert x_t\rvert^2$ . Dividing by $2\pi N$ and then multiplying by $2\pi/N$ gives the stated form.

Part 2. Write $m_t = w_t\mathrm{e}^{-i\lambda t}$ . Since $\mathbb{E}X_t = 0$ ,

\mathbb{E}\Bigl\lvert\sum_t w_tX_t\mathrm{e}^{-i\lambda t}\Bigr\rvert^{2} = \sum_{s}\sum_{t} w_sw_t\,\mathbb{E}(X_sX_t)\,\mathrm{e}^{-i(s-t)\lambda} = \sum_{s}\sum_{t} w_sw_t\,c(s-t)\,\mathrm{e}^{-i(s-t)\lambda} ,

using stationarity in the last step to write $\mathbb{E}(X_sX_t) = c(s-t)$ . Absolute summability of $c$ lets us substitute the spectral representation $c(s-t) = \int_{-\pi}^{\pi}\mathrm{e}^{i(s-t)\mu}g(\mu)\,d\mu$ of Stationary processes and spectral densities and exchange the finite double sum with the integral:

= \int_{-\pi}^{\pi}\Bigl(\sum_{s}\sum_{t}w_sw_t\,\mathrm{e}^{i(s-t)(\mu-\lambda)}\Bigr)g(\mu)\,d\mu = \int_{-\pi}^{\pi}\Bigl\lvert\sum_t w_t\,\mathrm{e}^{it(\mu-\lambda)}\Bigr\rvert^{2}g(\mu)\,d\mu ,

the bracket having factorised into a sum times its own conjugate because $w$ is real. Dividing by $2\pi\sum_tw_t^2$ turns the integrand into $\Phi^w_N(\mu-\lambda)g(\mu)$ , and $\Phi^w_N$ is even, which gives the stated form. That $\int\Phi^w_N = 1$ is the same computation with $g$ replaced by the constant $1/2\pi$ and $c$ by $\delta_{h0}$ , or directly from orthogonality of the exponentials.

For the untapered case set $w_t\equiv1$ in the double sum and group by $h = s-t$ as in part 1: the number of pairs with $s - t = h$ is $N - \lvert h\rvert$ , so

\mathbb{E}\,I_N(\lambda) = \frac{1}{2\pi N}\sum_{\lvert h\rvert<N}(N - \lvert h\rvert)\,c(h)\mathrm{e}^{-ih\lambda} = \frac{1}{2\pi}\sum_{\lvert h\rvert<N}\Bigl(1-\frac{\lvert h\rvert}{N}\Bigr)c(h)\mathrm{e}^{-ih\lambda} .

Subtracting this from $g(\lambda) = (2\pi)^{-1}\sum_{h}c(h)\mathrm{e}^{-ih\lambda}$ and taking absolute values,

\bigl\lvert\mathbb{E}I_N(\lambda) - g(\lambda)\bigr\rvert \;\le\; \frac{1}{2\pi}\sum_{\lvert h\rvert<N}\frac{\lvert h\rvert}{N}\lvert c(h)\rvert + \frac{1}{2\pi}\sum_{\lvert h\rvert\ge N}\lvert c(h)\rvert ,

with no $\lambda$ on the right. The second sum is the tail of a convergent series and goes to $0$ . For the first, fix $\varepsilon>0$ , choose $H$ with $\sum_{\lvert h\rvert\ge H}\lvert c(h)\rvert<\varepsilon$ , and split at $H$ : the part with $\lvert h\rvert<H$ is at most $(H/N)\sum_h\lvert c(h)\rvert \to 0$ , and the part with $\lvert h\rvert \ge H$ is at most $\varepsilon$ because $\lvert h\rvert/N \le 1$ . Hence the bound is eventually below $2\varepsilon/2\pi$ , which is uniform convergence. The closed form of $\Phi_N$ is the geometric sum $\sum_{t<N}\mathrm{e}^{-i\theta t} = \mathrm{e}^{-i(N-1)\theta/2}\sin(N\theta/2)/\sin(\theta/2)$ , whose modulus squared is $\sin^2(N\theta/2)/\sin^2(\theta/2)$ ; the sidelobe figures are the local maxima of that function, which for large $N$ occur near $\theta = 2\pi\cdot1.43/N$ at $4.7\%$ of the central value $N^2$ .

Part 3. Write $J_N(\lambda_j) = A_j - iB_j$ with

A_j = \sum_{t=0}^{N-1}X_t\cos(\lambda_j t), \qquad B_j = \sum_{t=0}^{N-1}X_t\sin(\lambda_j t) .

Both are linear combinations of the independent Gaussians $X_t$ , so the whole family $\{A_j, B_j\}$ is jointly Gaussian with mean zero (The multivariate normal distribution), and its law is fixed by the covariances. With $0<j<N/2$ and $0<k<N/2$ , independence of the $X_t$ gives

\operatorname{cov}(A_j,A_k) = \sigma^2\sum_t\cos(\lambda_jt)\cos(\lambda_kt), \qquad \operatorname{cov}(A_j,B_k) = \sigma^2\sum_t\cos(\lambda_jt)\sin(\lambda_kt),

and similarly for $\operatorname{cov}(B_j,B_k)$ . These are the discrete trigonometric orthogonality sums of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 20.6.1, p. 1007): the cosine-sine sum vanishes always, and the cosine-cosine and sine-sine sums vanish for $j\ne k$ and equal $N/2$ for $j = k$ in the stated range, the exceptional values $N$ and $0$ occurring only at $j = 0$ and $j = N/2$ .

So $A_j$ and $B_j$ are independent $N(0,\sigma^2N/2)$ variables, independent also of $A_k, B_k$ for $k \ne j$ — for jointly Gaussian variables zero covariance is independence. Therefore

\frac{A_j^2+B_j^2}{\sigma^2N/2}\;\sim\;\chi^2_2 , \qquad I_N(\lambda_j) = \frac{A_j^2+B_j^2}{2\pi N} = \frac{\sigma^2}{4\pi}\cdot\frac{A_j^2+B_j^2}{\sigma^2N/2} \;\sim\; \frac{\sigma^2}{4\pi}\chi^2_2 ,

by the definition of the chi-squared distribution as a sum of squares of independent standard normals (Sampling distributions). Since $g = \sigma^2/2\pi$ here, this is $I_N/g \sim \tfrac12\chi^2_2$ . A $\chi^2_2$ variable has mean $2$ and variance $4$ , so $\mathbb{E}I_N = g$ and $\operatorname{var}I_N = (g/2)^2\cdot4 = g^2$ ; and $\tfrac12\chi^2_2$ is the exponential distribution of mean $1$ , because the chi-squared density with two degrees of freedom is $\tfrac12\mathrm{e}^{-u/2}$ .

At $j=0$ , $B_0 = 0$ and $A_0 = \sum_tX_t \sim N(0,N\sigma^2)$ , so $I_N(0) = A_0^2/(2\pi N) \sim (\sigma^2/2\pi)\chi^2_1$ ; the same happens at $j = N/2$ for even $N$ , where $\sin(\pi t) = 0$ . This is the Hermitian symmetry $F_p = F^*_{N-p}$ of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 20.6.3(a), p. 1007) seen from the probabilistic side: a real record of length $N$ has only $N$ independent real numbers in its transform, two at each of the $\lfloor (N-1)/2\rfloor$ interior frequencies and one at each self-conjugate end.

Part 4. The $K$ segments are disjoint blocks of independent variables, hence independent, so it is enough to find the law of one tapered segment periodogram. Repeat the decomposition of part 3 with the taper carried along: for a single segment $X_0,\dots,X_{M-1}$ put

A_j = \sum_{t=0}^{M-1}w_tX_t\cos(\lambda_jt), \quad B_j = \sum_{t=0}^{M-1}w_tX_t\sin(\lambda_jt), \quad I^{w}_M(\lambda_j) = \frac{A_j^{2}+B_j^{2}}{2\pi\sum_tw_t^{2}} ,

with $\lambda_j = 2\pi j/M$ . The pair $(A_j,B_j)$ is centred Gaussian, and independence of the $X_t$ gives

\operatorname{var}A_j = \sigma^2\!\sum_t w_t^2\cos^2(\lambda_jt), \quad \operatorname{var}B_j = \sigma^2\!\sum_t w_t^2\sin^2(\lambda_jt), \quad \operatorname{cov}(A_j,B_j) = \sigma^2\!\sum_t w_t^2\cos(\lambda_jt)\sin(\lambda_jt) .

Adding the first two, $\operatorname{var}A_j + \operatorname{var}B_j = \sigma^2\sum_tw_t^2$ for every taper and every $j$ , so

\mathbb{E}\,I^{w}_M(\lambda_j) = \frac{\sigma^2\sum_tw_t^2}{2\pi\sum_tw_t^2} = \frac{\sigma^2}{2\pi} = g

exactly: the mean never depends on the taper. The difference, on the other hand, does. By the double-angle identities,

\operatorname{var}A_j - \operatorname{var}B_j = \sigma^2\sum_t w_t^2\cos(2\lambda_jt), \qquad 2\operatorname{cov}(A_j,B_j) = \sigma^2\sum_t w_t^2\sin(2\lambda_jt) ,

and these are the real and imaginary parts of $\sigma^2\sum_tw_t^2\mathrm{e}^{-2i\lambda_jt}$ up to a sign. So $A_j$ and $B_j$ are uncorrelated with the same variance $\tfrac12\sigma^2\sum_tw_t^2$ precisely when condition (T) holds. For $w_t\equiv1$ the two sums are $\sum_t\cos(2\lambda_jt)$ and $\sum_t\sin(2\lambda_jt)$ , which vanish unless $2\lambda_j$ is $0$ modulo $2\pi$ , that is unless $2j \equiv 0 \pmod M$ ; for $0<j<M/2$ we have $0 < 2j < M$ , so (T) holds at every interior frequency and part 3 is recovered.

Where (T) holds, $A_j$ and $B_j$ are independent $N(0,\tfrac12\sigma^2\sum_tw_t^2)$ — zero covariance is independence for a jointly Gaussian pair — so $(A_j^2+B_j^2)\big/\bigl(\tfrac12\sigma^2\sum_tw_t^2\bigr)\sim\chi^2_2$ and $I^w_M(\lambda_j)\sim(\sigma^2/4\pi)\chi^2_2 = \tfrac12 g\chi^2_2$ , as in part 3. A sum of $K$ independent $\chi^2_2$ variables is $\chi^2_{2K}$ , so $K\bar I/g\sim\tfrac12\chi^2_{2K}$ , that is $\bar I/g \sim \chi^2_{2K}/(2K)$ . Its mean is $2K/(2K)=1$ and its variance is $4K/(2K)^2 = 1/K$ , giving relative standard deviation $K^{-1/2}$ .

Where (T) fails, let $\Sigma$ be the covariance matrix of $(A_j,B_j)$ and $\sigma_1^2\ge\sigma_2^2$ its eigenvalues. Rotating the plane by an orthogonal matrix leaves $A_j^2+B_j^2$ unchanged and diagonalises $\Sigma$ , so $A_j^2+B_j^2$ has the law of $\sigma_1^2Z_1^2+\sigma_2^2Z_2^2$ with $Z_1,Z_2$ independent standard normal. Hence its mean is $\sigma_1^2+\sigma_2^2 = \operatorname{tr}\Sigma = \sigma^2\sum_tw_t^2$ , confirming the mean once more, while its variance is $2\sigma_1^4+2\sigma_2^4$ , so the relative standard deviation of one segment is

\frac{\sqrt{2(\sigma_1^4+\sigma_2^4)}}{\sigma_1^2+\sigma_2^2} \;\ge\; 1 ,

with equality only when $\sigma_1 = \sigma_2$ , that is only under (T); and a weighted sum of two chi-squares with unequal weights is not a multiple of $\chi^2_2$ . Averaging $K$ independent segments divides this ratio by $\sqrt K$ .

For the Hann taper $w_t = \tfrac12(1-\cos(2\pi t/M))$ one has $w_t^2 = \tfrac18\bigl(3 - 4\cos(2\pi t/M) + \cos(4\pi t/M)\bigr)$ , a trigonometric polynomial carrying only the frequencies $0,\pm2\pi/M,\pm4\pi/M$ . By discrete orthogonality, $\sum_tw_t^2\mathrm{e}^{-2i\lambda_jt}$ therefore vanishes unless $2\lambda_j$ is one of those frequencies modulo $2\pi$ , that is unless $2j\equiv0,\pm1,\pm2\pmod M$ ; for even $M$ and $0<j<M/2$ the only solutions are $j=1$ and $j=M/2-1$ . At those two, $\sum_tw_t^2 = 3M/8$ and $\sum_tw_t^2\cos(2\lambda_jt) = M/16$ with the sine sum zero, so $\sigma_1^2 = \tfrac{7}{32}\sigma^2M$ and $\sigma_2^2 = \tfrac{5}{32}\sigma^2M$ and the ratio above is $\sqrt{2(7^2+5^2)}/(7+5) = \sqrt{148}/12 = \sqrt{37}/6 = 1.0138$ , independently of $M$ . $\square$

Remark(Where the hypotheses are used)Original

Stationarity is used once, in part 2, to replace $\mathbb{E}(X_sX_t)$ by a function of $s-t$ . Nothing else needs it: part 1 is an algebraic identity about numbers. This is why a periodogram of a record with a trend in it is not wrong so much as meaningless — the trend is a huge low-frequency feature whose sidelobes contaminate the whole spectrum, and no amount of averaging removes it. Detrend first.
Zero mean is assumed in part 2 and matters for the same reason. Adding a constant $\mu$ to every $x_t$ adds $\mu^2\Phi_N(\lambda)$ to the expectation of the periodogram — the cross term between the constant and the process averages to zero — and $\mu^2\Phi_N$ is, by the definition of $\Phi_N$ , a copy of the whole spectral window, peak and sidelobes alike, centred on zero frequency. That is a term of size $\mu^2N/2\pi$ at $\lambda=0$ with sidelobes stretching across the band. Subtracting $\bar x$ removes it exactly; part 1 is stated for $j\ne0$ because that is where subtracting the mean leaves $J_N$ unchanged, and the mean-removed periodogram is identically zero at $j=0$ .
Absolute summability of $c$ is used twice in part 2: to justify exchanging the double sum with the integral, and to make the tail estimate that gives uniform convergence. It is a sufficient condition, not a necessary one — Grimmett and Stirzaker (Grimmett & Stirzaker (2020), (9.3.18), p. 442) say so where they state it — and it excludes in particular every process with a spectral atom, such as a deterministic sinusoid of random phase, whose spectral distribution has a jump and no density. For those the periodogram does something different and better: at the matching Fourier frequency its value grows linearly in $N$ , which is how a line is detected rather than estimated. But an atom is not the only obstruction, and the summability is genuinely stronger than having a density: it also fails for processes whose spectral measure is absolutely continuous and whose autocovariance merely decays too slowly, the long-memory spectra with a singularity at the origin. The density $g(\lambda) = \log 2 - \log\lvert 2\sin(\lambda/2)\rvert$ is non-negative, with minimum $g(\pi)=0$ , and integrable, with total mass $2\pi\log2 = 4.35517$ , so it carries no atom whatever; yet its autocovariances are $c(n) = \pi/\lvert n\rvert$ for $n\ne0$ , and $\sum_n\lvert c(n)\rvert$ diverges harmonically. More generally $c(n)\sim\lvert n\rvert^{-\alpha}$ with $0<\alpha\le1$ goes with an integrable $\lvert\lambda\rvert^{\alpha-1}$ singularity of $g$ at $\lambda=0$ and is never summable. Part 2 therefore says nothing about exactly the case in which periodogram bias is worst — a spectrum that climbs without bound at low frequency, where the window's sidelobes have an unbounded feature to draw power from.
Gaussianity is used only in part 3, to pass from uncorrelated to independent and to name the law. Without it, $\mathbb{E}I_N$ and $\operatorname{var}I_N$ are unchanged to leading order but the distribution is not chi-squared at finite $N$ . The qualitative conclusion — relative error of order one, not shrinking with $N$ — is not special to the Gaussian case.
Whiteness is used in part 3 to make the covariance sums collapse to the orthogonality relations. The corresponding statement for a general Gaussian stationary process with a continuous, strictly positive spectral density is that $I_N(\lambda_j)/g(\lambda_j) \to \tfrac12\chi^2_2$ in distribution, independently across fixed distinct frequencies; we state this without proof, and use only the white-noise case, which is exact at every $N$ and is also the case that matters for setting detection thresholds.
Condition (T) is what part 4 needs of the taper, and only for the distributional claim: the mean $\mathbb{E}\bar I = g$ holds for every taper at every frequency, because it depends on $\sum_tw_t^2$ alone. What (T) buys is that the cosine and sine amplitudes be uncorrelated and of equal variance, which is what makes the two of them a $\chi^2_2$ rather than a weighted sum of two $\chi^2_1$ s. Dropping it costs little in practice — for the Hann taper it fails at two frequencies out of $M/2$ , and inflates the standard deviation there by $1.4\%$ — but the law is then not exactly chi-squared, and the error bars of step 8 are correspondingly approximate at those frequencies.
Independence of segments is used in part 4, both to add the chi-squared variables and to divide the variance by $K$ . Welch's overlapped segments are not independent, so the effective count of degrees of freedom is not $2K$ and has to be worked out for the particular taper and overlap; quoting $2K$ for overlapped segments is a genuine, if usually small, misstatement of the error bars.
Equispaced samples are assumed everywhere. For a gapped or irregular record the exponentials are no longer orthogonal, $\Phi_N$ acquires large sidelobes at the gap frequencies, and the clean statements above all fail.

Method(Estimating a spectrum from a record)Original

Check the sampling. Frequencies above $\pi/\Delta$ cannot be represented and appear folded back onto lower ones; this is the aliasing of Arfken, Weber and Harris's (Arfken, Weber & Harris (2013), Example 20.6.1, p. 1006) Example 20.6.1 and of The sampling theorem. Aliasing is a property of the measurement, not of the analysis, and must be prevented by an anti-alias filter (Linear filters) before sampling, not after.
Remove the mean, then any trend, by least squares (Least-squares fitting); remove a known periodic cycle too if it is not the object of study. Record how many parameters you removed.
Choose the trade-off. Decide the frequency resolution $\delta\lambda$ you need; that fixes the segment length $M \approx 2\pi/\delta\lambda$ , and then $K = \lfloor N/M\rfloor$ segments fix the relative error $K^{-1/2}$ . Doing this consciously, in this order, is the entire design.
Taper each segment. A Hann taper $w_t = \tfrac12\bigl(1-\cos(2\pi t/M)\bigr)$ is the default: it lowers the first sidelobe from the $-13.3$ decibels of the untapered window to about $-31.5$ decibels, at the cost of a main lobe twice as wide. Use no taper only when the spectrum is known to be nearly flat, and a more aggressive taper when a strong line sits near a weak feature.
Transform and square. Compute $\sum_t w_tx_t\mathrm{e}^{-i\lambda_jt}$ by the FFT (The discrete Fourier transform and the FFT), take the squared modulus, divide by $2\pi\sum_tw_t^2$ .
Average the segments and, if more smoothness is wanted, additionally average over adjacent frequencies; $L$ neighbouring bins multiply the degrees of freedom by about $L$ and widen the resolution by the same factor.
Check the scale. Verify $\frac{2\pi}{M}\sum_j\bar I(\lambda_j) \approx$ the sample variance of the record. If it does not, a factor of $2\pi$ , of $N$ , or of $\Delta$ is wrong. Convert to physical units at the very end: density per unit angular frequency is $\Delta\bar I$ , per unit ordinary frequency $2\pi\Delta\bar I$ .
Put error bars on, and compare with a null. With $\nu$ degrees of freedom the interval $[\nu\bar I/q_{1-\alpha/2},\ \nu\bar I/q_{\alpha/2}]$ , with $q_p$ the $p$ -quantile of $\chi^2_\nu$ , covers $g(\lambda)$ with probability $1-\alpha$ (Confidence intervals); plot it on a logarithmic axis, where its width is the same at every frequency. A peak counts as detected only against a stated null spectrum, usually the red-noise spectrum of Autoregressive models and red noise rather than a flat one, and only after allowing for the number of frequencies searched (Hypothesis testing).

Example(A line, its leakage, and the bias of a smooth spectrum)Original

Part A: a line exactly on the grid. Take $N = 8$ and $x_t = 3 + 2\cos(\pi t/2)$ for $t = 0,\dots,7$ , that is the record $(5,3,1,3,5,3,1,3)$ . Its mean is $3$ ; subtracting it leaves $y_t = (2,0,-2,0,2,0,-2,0)$ . The frequency $\pi/2$ is $\lambda_2 = 2\pi\cdot2/8$ , a Fourier frequency, so write $y_t = \mathrm{e}^{i\lambda_2 t} + \mathrm{e}^{-i\lambda_2t}$ . Then

J_8(\lambda_2) = \sum_{t=0}^{7}\bigl(\mathrm{e}^{i\lambda_2t}+\mathrm{e}^{-i\lambda_2t}\bigr)\mathrm{e}^{-i\lambda_2t} = 8 + 0 = 8 ,

the second sum vanishing by orthogonality, and likewise $J_8(\lambda_6) = 8$ while $J_8(\lambda_j) = 0$ for every other $j$ . Hence

I_8(\lambda_2) = I_8(\lambda_6) = \frac{8^2}{2\pi\cdot8} = \frac{4}{\pi} = 1.27324 ,

and zero elsewhere. Check with part 1: $\frac{2\pi}{8}\bigl(2\cdot\frac4\pi\bigr) = 2$ , and the mean square of $y$ is $(4+4+4+4)/8 = 2$ . The two agree. Note what happens if the mean is not removed: $J_8(0) = 24$ , so $I_8(0) = 576/(16\pi) = 11.459$ , nine times the signal peak, sitting where no signal is.

Part B: the same line half a bin off, with and without a taper. Now take $N = 16$ and $x_t = \cos(2\pi\cdot3.5\,t/16)$ , a frequency exactly between $\lambda_3$ and $\lambda_4$ . The record no longer contains a whole number of cycles, so no single $J_{16}(\lambda_j)$ captures it. Computing the sixteen values gives, for the untapered periodogram,

I_{16}(\lambda_3) = 0.25121, \quad I_{16}(\lambda_4) = 0.27130, \quad I_{16}(\lambda_2) = 0.02892, \quad I_{16}(\lambda_5) = 0.03538, \ \dots

with total $\sum_jI_{16}(\lambda_j) = 4/\pi = 1.27324$ . The two bins flanking the true frequency hold

\frac{2(0.25121+0.27130)}{1.27324} = 82.08\%

of the power (the factor $2$ counts the mirror bins $\lambda_{12},\lambda_{13}$ ), the nearest four hold $92.18\%$ , and the remaining $7.8\%$ is spread over the rest of the spectrum. That residue is leakage, and it is what a weak feature elsewhere in the spectrum has to compete with. Repeating with the Hann taper $w_t = \tfrac12(1-\cos(2\pi t/16))$ gives $96.14\%$ in the two nearest bins and $99.93\%$ in the nearest four: leakage falls by more than two orders of magnitude. The taper concentrates the peak better here, which is a feature of the half-bin offset: it is the offset that excites the untapered sidelobes maximally, and those sidelobes, $18$ decibels above the Hann ones, carry away more than the wider Hann main lobe does. The price of the taper is paid not here but in resolution, because the main lobe is twice as wide, $8\pi/N$ against $4\pi/N$ : two lines between one and two bins apart, which the untapered window still shows as two peaks, the Hann window merges into one.

Part C: the bias on a smooth spectrum. Let $X_n = Z_n + \theta Z_{n-1}$ with the $Z_n$ uncorrelated of mean $0$ and variance $\sigma^2$ . This is the moving average process of Grimmett and Stirzaker's (Grimmett & Stirzaker (2020), §9.7, Problem 7(a), p. 468) Problem 7 in §9.7, who write the coefficient $\alpha$ for our $\theta$ and take the $Z_n$ to have unit variance; we keep a general $\sigma^2$ , which only multiplies everything below. Its autocovariance is $c(0) = (1+\theta^2)\sigma^2$ , $c(\pm1) = \theta\sigma^2$ and $c(h) = 0$ for $\lvert h\rvert\ge2$ , so by the inversion formula

g(\lambda) = \frac{\sigma^2}{2\pi}\bigl(1 + \theta^2 + 2\theta\cos\lambda\bigr) .

Because only the lags $0$ and $\pm1$ are present, part 2 gives the expectation of the untapered periodogram in closed form at every $N$ :

\mathbb{E}\,I_N(\lambda) = \frac{\sigma^2}{2\pi}\Bigl(1+\theta^2 + 2\theta\Bigl(1-\frac1N\Bigr)\cos\lambda\Bigr) = g(\lambda) - \frac{\theta\sigma^2\cos\lambda}{\pi N} .

Put $\theta = \tfrac12$ , $\sigma^2 = 1$ , $N = 8$ . At $\lambda = 0$ the truth is $g(0) = 2.25/2\pi = 0.358099$ and the expectation is $2.125/2\pi = 0.338204$ , low by $5.56\%$ . At $\lambda = \pi$ the truth is $g(\pi) = 0.25/2\pi = 0.039789$ and the expectation is $0.375/2\pi = 0.059683$ , high by exactly $50\%$ . The absolute bias is the same at both ends, $0.019894 = 1/(16\pi)$ ; it is the relative bias that explodes where the spectrum is small, which is the general rule and the reason spectra are plotted on a logarithmic axis. At $\lambda = \pi/2$ , where $\cos\lambda = 0$ , the periodogram is exactly unbiased at every $N$ .

ApplicationOriginal

Detecting periodicity against a coloured background. The periodogram of a palaeoclimate record is how the Milankovitch periods near $19$ , $23$ , $41$ and $100$ thousand years were established; the periodogram of a stellar light curve is how a pulsation mode is identified. In both cases the finding is not the peak but the peak's height relative to a fitted null spectrum and to the number of frequencies searched. Part 3 is the reason this must be done carefully: a raw periodogram of pure red noise produces, at each of $N/2$ frequencies independently, an exponential variable, so the largest of them exceeds its mean by a factor $\ln(N/2)$ as a matter of course — about $6$ for a record of a thousand points. Peaks of that size are the null hypothesis, not evidence against it.

Instrument noise characterisation. The specification of an amplifier, an accelerometer or a gravitational-wave detector is quoted as an amplitude spectral density, $\sqrt{g(\lambda)}$ in units of signal per root hertz, and is measured by exactly the Welch procedure of part 4: minutes of data, split into hundreds of overlapping tapered segments, averaged. Here the smooth broadband shape is the quantity of interest and resolution is worth little, so $K$ is taken as large as the record allows. The narrow lines that stick out of such a spectrum — mains hum and its harmonics, mechanical resonances — are diagnosed by the opposite choice, one long untapered transform.

Turbulence and the inertial range. Measuring the exponent in a power-law spectrum $g(\lambda)\propto\lambda^{-5/3}$ is a case where leakage, not variance, is the enemy: the spectrum spans many decades, so the sidelobes of the untapered window at the energetic large scales sit above the true signal at the small scales and flatten the measured slope. Part 2 quantifies this, and the remedy is the taper of step 4 rather than more averaging.

Method

Linear filters

Adapted

Add to PDF

Shows that a linear time-invariant operation on a signal is a convolution, that complex exponentials pass through it unchanged except for multiplication by its transfer function, and that it multiplies the power spectral density of a stationary process by the squared gain; then turns those facts into a design recipe for moving averages and for low-pass, high-pass and band-pass filters, with the Gibbs ripple and the causality constraint that limit what can be built.

Prerequisites: The convolution theorem, The Fourier transform, Fourier inversion, The discrete Fourier transform and the FFT, Mean-square convergence and Parseval's identity, The Gibbs phenomenon, Cesàro means and Fejér's theorem, Stationary processes and spectral densities

Definition(Filters, impulse response, transfer function)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.5, Eqs. (20.100)-(20.106), pp. 997-999

A filter is a map $L$ taking an input signal to an output signal which is

linear: $L(\alpha x_1 + \beta x_2) = \alpha Lx_1 + \beta Lx_2$ , and
time-invariant: delaying the input by any amount delays the output by the same amount and changes nothing else.

In continuous time the input is $x:\mathbb{R}\to\mathbb{C}$ and the output is

y(t) \;=\; (\Phi * x)(t) \;=\; \int_{-\infty}^{\infty}\Phi(t-t')\,x(t')\,dt' ,

with $\Phi$ the impulse response, the output produced by a unit impulse at the origin. Its Fourier transform

\varphi(\omega) \;=\; \int_{-\infty}^{\infty}\Phi(t)\,\mathrm{e}^{-i\omega t}\,dt

is the transfer function; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.5, Eq. (20.102), p. 998) introduce it the other way round, as the factor by which the device multiplies a pure input frequency, and draw it as a labelled box between an input and an output wire.

In discrete time the input is a sequence $(x_n)$ and the filter is given by its weights $(a_k)_{k\in\mathbb{Z}}$ ,

y_n \;=\; \sum_{k\in\mathbb{Z}} a_k\,x_{n-k} ,

assumed throughout to satisfy $\sum_k\lvert a_k\rvert<\infty$ . Its transfer function is the $2\pi$ -periodic function

A(\lambda) \;=\; \sum_{k\in\mathbb{Z}} a_k\,\mathrm{e}^{-ik\lambda}, \qquad \lambda\in(-\pi,\pi] ,

whose modulus $\lvert A(\lambda)\rvert$ is the gain and whose argument $-\arg A(\lambda)$ is the phase lag at frequency $\lambda$ . The power gain is $\lvert A\rvert^2$ , and a frequency at which the power gain is half its maximum is a half-power point.

A filter is causal when the output at $t$ uses no input from after $t$ : $\Phi(t) = 0$ for $t<0$ , or $a_k = 0$ for $k<0$ . It is finite (FIR) when only finitely many $a_k$ are non-zero, and recursive (IIR) when $y_n$ is computed from earlier outputs as well as inputs; the autoregressive schemes of Autoregressive models and red noise are the standard recursive filters, and are not treated here.

Three special cases name themselves. A filter is a low-pass when $\lvert A\rvert$ is near $1$ for small $\lvert\lambda\rvert$ and near $0$ near $\lambda = \pi$ , a high-pass when the roles are exchanged, and a band-pass when $\lvert A\rvert$ is near $1$ only on an interval $0<\lambda_1<\lvert\lambda\rvert<\lambda_2<\pi$ . A moving average is a low-pass filter with non-negative weights summing to $1$ .

Method(What a filter does to a signal and to a spectrum)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.5, Eqs. (20.102)-(20.106) and (20.110)-(20.113), pp. 998-1001; Grimmett & Stirzaker (2020), §9.7, Problem 2, p. 468

Let $L$ be a filter as above.

1. Complex exponentials are its eigenfunctions. In discrete time, feeding $x_n = \mathrm{e}^{in\lambda}$ into the filter gives

y_n \;=\; A(\lambda)\,\mathrm{e}^{in\lambda} ,

and in continuous time $x(t) = \mathrm{e}^{i\omega t}$ gives $y(t) = \varphi(\omega)\mathrm{e}^{i\omega t}$ . A filter changes the amplitude and the phase of each frequency and nothing else; it creates no new frequency. This is the assumption Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.5, Eq. (20.102), p. 998) make explicitly as the definition of a linear response, listing its three parts: same frequency out as in, output proportional to input, and no interference between frequencies present at once.

2. Equivalently, the transform multiplies. $\widehat{Lx}(\omega) = \varphi(\omega)\hat x(\omega)$ , by the convolution theorem (The convolution theorem). Running the identification backwards recovers the convolution: because the output $\frac{1}{2\pi}\int\varphi(\omega)\hat x(\omega)\mathrm{e}^{i\omega t}d\omega$ is the inverse transform of a product, it is a convolution, which is how Arfken, Weber and Harris reach their Eq. (20.104). For a cascade of two filters the transfer functions multiply, $\varphi = \varphi_2\varphi_1$ , and the impulse responses convolve.

3. Real weights make the gain even and the phase odd; symmetric weights make the phase zero. If every $a_k$ is real then $A(-\lambda) = \overline{A(\lambda)}$ , so $\lvert A\rvert$ is even and $\arg A$ is odd, and a real input gives a real output. If in addition $a_{-k} = a_k$ , then

A(\lambda) = a_0 + 2\sum_{k\ge1}a_k\cos k\lambda

is real, so the filter delays nothing: the only phase it can apply is a sign, $0$ where $A>0$ and $\pi$ where $A<0$ . A symmetric filter is never causal, apart from the trivial multiple of the identity, and a causal filter is never zero-phase; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.5, Eqs. (20.110)-(20.113), pp. 1000-1001) show the continuous-time form of the obstruction, that causality ties the real and imaginary parts of $\varphi$ together through the relation

\int_0^{\infty}u(\omega)\cos\omega t\,d\omega = -\int_0^{\infty}v(\omega)\sin\omega t\,d\omega \qquad (t>0), \qquad \varphi = u+iv .

A filter run forwards and then backwards over a stored record has transfer function $\lvert A\rvert^2$ , which is real and non-negative: that is how zero phase is obtained in practice from a causal design.

4. Spectral transfer. Let $X = \{X_n\}$ be real, weakly stationary, with autocovariance $c_X$ and power spectral density $g_X(\lambda) = (2\pi)^{-1}\sum_n c_X(n)\mathrm{e}^{-in\lambda}$ , and let $Y_n = \sum_k a_kX_{n-k}$ with real weights and $\sum_k\lvert a_k\rvert<\infty$ . Then $Y$ is weakly stationary, with

c_Y(m) \;=\; \sum_{j}\sum_{k} a_ja_k\,c_X(m-k+j) \qquad\text{and}\qquad \boxed{\;g_Y(\lambda) \;=\; \bigl\lvert A(\lambda)\bigr\rvert^{2}\,g_X(\lambda)\;}

for every $\lambda\in(-\pi,\pi]$ . In particular, for white noise input of variance $\sigma^2$ , where $g_X \equiv \sigma^2/2\pi$ ,

g_Y(\lambda) = \frac{\sigma^2}{2\pi}\lvert A(\lambda)\rvert^{2}, \qquad \operatorname{var}Y_n = \sigma^{2}\sum_k a_k^{2} .

Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §9.7, Problem 2, p. 468) state this as Problem 2 of §9.7, for a one-sided finite filter $X_n = \sum_{k=0}^{r}a_kY_{n-k}$ and in terms of the normalised spectral density of the autocorrelation function, in which form a ratio of variances appears:

f_{\text{out}}(\lambda) = \frac{\sigma_{\text{in}}^2}{\sigma_{\text{out}}^2}\,f_{\text{in}}(\lambda)\,\bigl\lvert G_a(\mathrm{e}^{i\lambda})\bigr\rvert^{2}, \qquad G_a(z) = \sum_{k=0}^{r}a_kz^{k} .

Multiplying both sides by $\sigma^2_{\text{out}}$ turns their statement into the boxed one, since $g = \sigma^2f$ and $\lvert G_a(\mathrm{e}^{i\lambda})\rvert = \lvert A(\lambda)\rvert$ for real weights. The boxed form is the one used below.

Idea of the proofOriginal

A filter commutes with time shifts, and the complex exponentials are exactly the signals that a time shift multiplies by a constant. So in the basis of exponentials the filter must be diagonal, and the number sitting on the diagonal at frequency $\lambda$ is the transfer function. Everything in the statement is a restatement of that one sentence: the convolution form is the diagonal operator written in the time domain, the cascade rule is that diagonal matrices commute and multiply entrywise, and the spectral transfer is the same diagonalisation applied to the second moments rather than to the signal.

For part 4 the only extra ingredient is that the autocovariance is itself a superposition of exponentials, $c_X(h) = \int\mathrm{e}^{ih\lambda}g_X(\lambda)d\lambda$ . Expanding $c_Y$ as a double sum over the weights and substituting that representation lets the two sums be done separately, one producing $\sum_ja_j\mathrm{e}^{ij\lambda}$ and the other its complex conjugate; their product is $\lvert A\rvert^2$ . The squaring is not an accident of the algebra: covariance is quadratic in the signal, so a filter that halves an amplitude quarters the power at that frequency.

Proof(Parts 1-4)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.5, pp. 998-999; Grimmett & Stirzaker (2020), §9.7, Problem 2, p. 468

Part 1. Substitute $x_n = \mathrm{e}^{in\lambda}$ into the definition:

y_n = \sum_k a_k\,\mathrm{e}^{i(n-k)\lambda} = \Bigl(\sum_k a_k\mathrm{e}^{-ik\lambda}\Bigr)\mathrm{e}^{in\lambda} = A(\lambda)\,\mathrm{e}^{in\lambda} ,

the factorisation being legitimate because $\mathrm{e}^{i(n-k)\lambda} = \mathrm{e}^{in\lambda}\mathrm{e}^{-ik\lambda}$ and the series converges absolutely, by the hypothesis $\sum_k\lvert a_k\rvert<\infty$ and $\lvert\mathrm{e}^{-ik\lambda}\rvert = 1$ . The continuous-time computation is the same with an integral: $\int\Phi(t-t')\mathrm{e}^{i\omega t'}dt' = \mathrm{e}^{i\omega t}\int\Phi(u)\mathrm{e}^{-i\omega u}du$ after the substitution $u = t-t'$ .

Part 2. The identity $\widehat{\Phi*x} = \hat\Phi\hat x = \varphi\hat x$ is the convolution theorem of The convolution theorem, which needs only $\Phi, x\in L^1(\mathbb{R})$ . For the cascade, apply it twice: $\widehat{L_2L_1x} = \varphi_2\widehat{L_1x} = \varphi_2\varphi_1\hat x$ , and by uniqueness of the transform (Fourier inversion) the composite filter is the one with transfer function $\varphi_2\varphi_1$ , whose impulse response is $\Phi_2*\Phi_1$ .

Part 3. If the $a_k$ are real then $\overline{A(\lambda)} = \sum_k a_k\mathrm{e}^{ik\lambda} = A(-\lambda)$ , from which $\lvert A(-\lambda)\rvert = \lvert A(\lambda)\rvert$ and $\arg A(-\lambda) = -\arg A(\lambda)$ . If moreover $a_{-k} = a_k$ , pair the terms $k$ and $-k$ :

A(\lambda) = a_0 + \sum_{k\ge1}a_k\bigl(\mathrm{e}^{-ik\lambda}+\mathrm{e}^{ik\lambda}\bigr) = a_0 + 2\sum_{k\ge1}a_k\cos k\lambda ,

a real number. A real $A$ has argument $0$ or $\pi$ , which is the claim. Symmetry forces $a_{-k} = a_k \ne 0$ for some $k>0$ unless the filter is a multiple of the identity, so a symmetric filter is not causal. For the forward-backward trick, running the filter backwards over a stored record applies the weights $a_{-k}$ , whose transfer function is $\sum_ka_{-k}\mathrm{e}^{-ik\lambda} = \overline{A(\lambda)}$ for real weights, and by the cascade rule the two passes compose to $A\overline{A} = \lvert A\rvert^2$ . The causality relations are Arfken, Weber and Harris's (Arfken, Weber & Harris (2013), §20.5, Eqs. (20.111)-(20.113), pp. 1000-1001) Eqs. (20.111)-(20.113), obtained by splitting $\varphi(\omega) = \int_0^{\infty}\Phi(t)\mathrm{e}^{-i\omega t}dt$ into real and imaginary parts, observing that the first is even and the second odd in $\omega$ , and inverting each as a cosine and as a sine transform of the same $\Phi$ .

Part 4. $Y_n$ has constant mean $\mathbb{E}Y_n = (\sum_ka_k)\mathbb{E}X_0$ , finite because $\sum_k\lvert a_k\rvert<\infty$ . Take $\mathbb{E}X_n = 0$ without loss of generality, by subtracting the mean from both sides. Then

c_Y(m) = \mathbb{E}\bigl(Y_nY_{n+m}\bigr) = \mathbb{E}\Bigl(\sum_ja_jX_{n-j}\sum_ka_kX_{n+m-k}\Bigr) = \sum_j\sum_k a_ja_k\,c_X\bigl(m-k+j\bigr) ,

where the expectation has been taken inside the double sum — justified because $\sum_{j,k}\lvert a_ja_k\rvert\,\mathbb{E}\lvert X_{n-j}X_{n+m-k}\rvert \le c_X(0)\bigl(\sum_k\lvert a_k\rvert\bigr)^2 <\infty$ by the Cauchy-Schwarz inequality applied to each term — and stationarity of $X$ has been used to write $\mathbb{E}(X_{n-j}X_{n+m-k}) = c_X((n+m-k)-(n-j))$ . The right-hand side depends on $m$ alone, so $Y$ is weakly stationary.

Now substitute the spectral representation $c_X(h) = \int_{-\pi}^{\pi}\mathrm{e}^{ih\lambda}g_X(\lambda)\,d\lambda$ of Stationary processes and spectral densities and exchange the sums with the integral, which is permitted because $\sum_{j,k}\lvert a_ja_k\rvert\int\lvert g_X\rvert = c_X(0)(\sum\lvert a_k\rvert)^2<\infty$ and $g_X\ge0$ :

c_Y(m) = \int_{-\pi}^{\pi}\mathrm{e}^{im\lambda}\Bigl(\sum_j a_j\mathrm{e}^{ij\lambda}\Bigr)\Bigl(\sum_k a_k\mathrm{e}^{-ik\lambda}\Bigr)g_X(\lambda)\,d\lambda = \int_{-\pi}^{\pi}\mathrm{e}^{im\lambda}\,\bigl\lvert A(\lambda)\bigr\rvert^{2}g_X(\lambda)\,d\lambda ,

the last step because for real weights $\sum_ja_j\mathrm{e}^{ij\lambda} = \overline{A(\lambda)}$ while $\sum_ka_k\mathrm{e}^{-ik\lambda} = A(\lambda)$ . So $c_Y$ is the Fourier coefficient sequence of the function $\lvert A\rvert^2g_X$ , which is non-negative and integrable; by the uniqueness of the spectral distribution (Stationary processes and spectral densities) this function is the power spectral density of $Y$ , which is the boxed identity.

For the white-noise specialisation, $g_X = \sigma^2/2\pi$ is constant, so $g_Y = \sigma^2\lvert A\rvert^2/2\pi$ and

\operatorname{var}Y_n = c_Y(0) = \int_{-\pi}^{\pi}\frac{\sigma^2}{2\pi}\lvert A(\lambda)\rvert^2\,d\lambda = \sigma^2\sum_k a_k^2 ,

the last equality being Parseval's identity for the Fourier series $A$ (Mean-square convergence and Parseval's identity); it also follows directly from the double sum with $c_X(h) = \sigma^2\delta_{h0}$ . $\square$

Remark(Where the hypotheses are used)Original

Absolute summability of the weights is used three times: to factor the series in part 1, to take the expectation inside the double sum in part 4, and to make $A$ a continuous function. It fails for the one filter everybody wants. The ideal low-pass, $A(\lambda) = 1$ for $\lvert\lambda\rvert\le\lambda_c$ and $0$ beyond, has weights $a_0 = \lambda_c/\pi$ and $a_k = \sin(k\lambda_c)/(\pi k)$ , which decay like $1/k$ and are not absolutely summable. Truncating them to $\lvert k\rvert\le K$ gives a filter whose gain overshoots to about $1.09$ just inside the passband and dips to about $-0.09$ just inside the stopband, however large $K$ is: this is the Gibbs phenomenon of The Gibbs phenomenon, and a negative gain in what was supposed to be a stopband is a real defect, not a cosmetic one. Multiplying the truncated weights by a taper (the Bartlett weights $1-\lvert k\rvert/K$ give the Fejer means of Cesàro means and Fejér's theorem) removes the overshoot entirely, at the price of a transition band about twice as wide. A brick wall is not available at any finite cost.
Time invariance is what makes exponentials eigenfunctions, and it is the hypothesis quietly violated whenever the weights are allowed to depend on position — an edge-preserving or adaptive smoother is not a filter in this sense and none of part 1 to part 4 applies to it.
Stationarity of the input is used only in part 4, where it turns $\mathbb{E}(X_sX_t)$ into a function of $s-t$ . Parts 1 to 3 are statements about one signal and need no probability at all.
Real weights are used in part 3 and part 4. For complex weights the correct general statement is Grimmett and Stirzaker's (Grimmett & Stirzaker (2020), §9.7, Problem 2, p. 468), with $\lvert G_a(\mathrm{e}^{i\lambda})\rvert^2 = \lvert\sum_ka_k\mathrm{e}^{ik\lambda}\rvert^2$ , which is $\lvert A(-\lambda)\rvert^2$ and not $\lvert A(\lambda)\rvert^2$ . The two agree for real weights, which is every filter applied to real data, and the discrepancy is traceable to the conjugate in the definition of the covariance of complex random variables.
A filter cannot create power where there is none. Since $g_Y = \lvert A\rvert^2g_X$ pointwise, $g_X(\lambda) = 0$ forces $g_Y(\lambda) = 0$ . In particular no filter can undo aliasing: once two frequencies have been folded onto one by sampling, no operation on the samples separates them, which is why the anti-alias filter of The sampling theorem must be analogue and must act before the sampler.
Deconvolution is division, and division is unstable. Inverting a filter means multiplying by $1/A(\lambda)$ , which is unbounded wherever $A$ is small — and a low-pass filter is designed to make $A$ small. The inverse of a three-point running mean does not exist as a bounded filter at all, because $A(2\pi/3) = 0$ exactly. This is the ill-posedness already noted in The convolution theorem.
Smoothing manufactures autocorrelation. White noise through any filter comes out correlated, with $\rho_Y(m) = \sum_ka_ka_{k+m}\big/\sum_ka_k^2$ ; for both three-point filters of the example below this is $2/3$ at lag $1$ . A smoothed series therefore has far fewer independent values than it has points, and any test applied to it afterwards must use the reduced count, not $N$ .
Causality is a constraint, not a property of the mathematics. Part 3 says a zero-phase filter needs future values. That is unavailable in real time and perfectly available when post-processing a stored record, which is why real-time and offline filtering are different design problems with different answers.

Method(Designing and applying a filter)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.5, Eq. (20.102) and Exercises 20.5.3-20.5.4, pp. 998, 1002; Riley, Hobson & Bence (2006), §13.4, Exercise 13.12, p. 462

Write down the gain you want, as a function $\lvert A_{\text{target}}(\lambda)\rvert$ on $[0,\pi]$ : where it should be $1$ , where it should be $0$ , and how wide a transition you will tolerate between them. This is a specification, not yet a filter, and the width of the transition is the quantity that will cost you coefficients.
Decide the phase. Zero phase (symmetric weights, offline) or causal (one-sided weights, real time). If zero phase is wanted from a causal design, plan to run the filter forwards and then backwards, which squares the gain, so specify $\sqrt{\lvert A_{\text{target}}\rvert}$ for the single pass.
Invert the Fourier series. The weights are $a_k = \frac{1}{2\pi}\int_{-\pi}^{\pi}A_{\text{target}}(\lambda)\mathrm{e}^{ik\lambda}\,d\lambda$ . For a low-pass with cutoff $\lambda_c$ this gives $a_k = \sin(k\lambda_c)/(\pi k)$ ; for a high-pass, subtract a low-pass from the identity, $a_k^{\text{hp}} = \delta_{k0}-a_k^{\text{lp}}$ ; for a band-pass with edges $\lambda_1<\lambda_2$ , subtract two low-passes, $a_k = a_k^{(\lambda_2)}-a_k^{(\lambda_1)}$ . Cascading two filters multiplies their transfer functions, which is the other way to build a band-pass, and is what Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercises 20.5.3-20.5.4 and Fig. 20.12, p. 1002) ask for in Exercises 20.5.3 and 20.5.4 — with the warning that a physical cascade only multiplies if the second stage does not load the first.
Truncate and taper. Keep $\lvert k\rvert\le K$ and multiply by a taper such as $1-\lvert k\rvert/(K+1)$ , or use a smooth target in step 1 instead of a brick wall. Without this step, see the Gibbs remark above.
Renormalise. For a low-pass, divide by $\sum_ka_k$ so that $A(0) = 1$ exactly and constants pass untouched; for a high-pass or band-pass, check instead that $\sum_ka_k = 0$ exactly, so that a constant is annihilated and no trend leaks through.
Plot the realised gain $A(\lambda) = a_0+2\sum_{k\ge1}a_k\cos k\lambda$ before using the filter, on a logarithmic scale, and read off the actual half-power point and the stopband level. The realised gain, not the target, is what the data will see.
Apply it, and deal with the ends. A filter of half-width $K$ loses $K$ points at each end of a finite record. Either accept the loss, or state what you padded with; padding with zeros is equivalent to multiplying the record by a rectangular window and brings back the leakage of The periodogram and spectral estimation.
Check against the spectrum. Verify that the periodogram of the output is the periodogram of the input times $\lvert A\rvert^2$ , to within the sampling error of The periodogram and spectral estimation, and that the variance ratio is $\sum_ka_k^2$ for a white input. For an analogue filter the corresponding check is on the impulse response: Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §13.4, Exercise 13.12, p. 462) Exercise 13.12 is the model case, where a response shaped as a triangle of base $2T$ — two rectangles convolved — has for its transfer function the square of a sinc, decaying like $\omega^{-2}$ because the response is continuous but has a corner.

Example(Two three-point smoothers, and an

RC

low-pass)Original

Part A: gains. Consider the three-point running mean $a = (\tfrac13,\tfrac13,\tfrac13)$ and the binomial smoother $a = (\tfrac14,\tfrac12,\tfrac14)$ , both symmetric about $k=0$ . By part 3,

A_{\text{mean}}(\lambda) = \tfrac13\bigl(1+2\cos\lambda\bigr), \qquad A_{121}(\lambda) = \tfrac12\bigl(1+\cos\lambda\bigr) = \cos^{2}\tfrac{\lambda}{2} .

Both have $A(0) = 1$ . But $A_{\text{mean}}$ vanishes at $\cos\lambda = -\tfrac12$ , that is at $\lambda = 2\pi/3$ , and is negative beyond: at the Nyquist frequency $A_{\text{mean}}(\pi) = \tfrac13(1-2) = -\tfrac13$ . The running mean therefore inverts every oscillation faster than three samples per cycle, while $A_{121}\ge0$ everywhere and decreases monotonically to $A_{121}(\pi) = 0$ . The half-power point of the binomial smoother solves $\cos^4(\lambda/2) = \tfrac12$ , giving

\lambda_{1/2} = 2\arccos\bigl(2^{-1/4}\bigr) = 1.14372\ \text{rad} = 0.36406\,\pi .

Part B: the sign reversal, seen in integers. Take the seven-point record

x_t = t + (-1)^{t}, \qquad t = 0,\dots,6, \qquad\text{that is}\qquad x = (1,0,3,2,5,4,7) ,

a straight line with a Nyquist-frequency wobble of unit amplitude on top. Both filters pass the line exactly, because both have $\sum_ka_k = 1$ and $\sum_kka_k = 0$ by symmetry. The wobble is $\mathrm{e}^{i\pi t}$ , so by part 1 it comes out multiplied by $A(\pi)$ .

For the binomial smoother, $A_{121}(\pi) = 0$ and the wobble should vanish completely. Computing $y_t = (x_{t-1}+2x_t+x_{t+1})/4$ for $t = 1,\dots,5$ :

\tfrac{1+0+3}{4} = 1,\quad \tfrac{0+6+2}{4} = 2,\quad \tfrac{3+4+5}{4} = 3,\quad \tfrac{2+10+4}{4} = 4,\quad \tfrac{5+8+7}{4} = 5 ,

exactly $(1,2,3,4,5)$ — the line, with the wobble gone.

For the running mean, $A_{\text{mean}}(\pi) = -\tfrac13$ , so the wobble should survive at a third of its amplitude and with its sign flipped. Computing $y_t = (x_{t-1}+x_t+x_{t+1})/3$ :

\tfrac43,\quad \tfrac53,\quad \tfrac{10}{3},\quad \tfrac{11}{3},\quad \tfrac{16}{3} ,

whose residuals about the line $(1,2,3,4,5)$ are $+\tfrac13,-\tfrac13,+\tfrac13,-\tfrac13,+\tfrac13$ . The prediction is $A_{\text{mean}}(\pi)\cdot(-1)^t = -\tfrac13(-1)^t$ , which at $t=1,2,3,4,5$ is $+\tfrac13,-\tfrac13,+\tfrac13,-\tfrac13,+\tfrac13$ . The two agree exactly, including the sign.

Part C: variance reduction on white noise. By part 4, $\operatorname{var}Y/\operatorname{var}X = \sum_ka_k^2$ , which is $3\cdot\tfrac19 = \tfrac13 = 0.3333$ for the running mean and $\tfrac1{16}+\tfrac14+\tfrac1{16} = \tfrac38 = 0.3750$ for the binomial smoother. The cross-check through the spectrum is Parseval: for the binomial filter,

\frac{1}{2\pi}\int_{-\pi}^{\pi}\cos^{4}\tfrac{\lambda}{2}\,d\lambda = \frac38 ,

since the mean value of $\cos^4$ over a whole period is $3/8$ . So the running mean removes slightly more variance — and pays for it by reversing the sign of the highest frequencies. That is the trade the two filters illustrate, and it is why the binomial weights, not the equal weights, are the default three-point smoother.

Part D: the continuous analogue. A resistance $R$ in series with a capacitance $C$ to ground, with the output read across the capacitor, is the low-pass companion of the circuit in Arfken, Weber and Harris's (Arfken, Weber & Harris (2013), Example 20.5.1 and Exercise 20.5.1, pp. 999-1001) Example 20.5.1. Driving it at a single frequency, with $V_{\text{in}}\mathrm{e}^{i\omega t}$ , $V_{\text{out}}\mathrm{e}^{i\omega t}$ and $I\mathrm{e}^{i\omega t}$ the complex amplitudes, Kirchhoff's loop equation and its derivative give

V_{\text{in}}\mathrm{e}^{i\omega t} = RI\mathrm{e}^{i\omega t} + \frac1C\int\! I\mathrm{e}^{i\omega t}dt \;\Longrightarrow\; i\omega V_{\text{in}} = i\omega RI + \frac{I}{C} \;\Longrightarrow\; I = \frac{i\omega C\,V_{\text{in}}}{1+i\omega RC} .

The output is the capacitor voltage $V_{\text{out}} = I/(i\omega C)$ , so

\varphi(\omega) = \frac{V_{\text{out}}}{V_{\text{in}}} = \frac{1}{1+i\omega RC}, \qquad \lvert\varphi(\omega)\rvert = \frac{1}{\sqrt{1+(\omega RC)^2}}, \qquad \arg\varphi(\omega) = -\arctan(\omega RC) .

The power gain is $1/2$ exactly at $\omega = 1/(RC)$ , and $\lvert\varphi\rvert$ falls like $1/(\omega RC)$ above it — a roll-off of six decibels per octave. The impulse response is the causal exponential $\Phi(t) = (RC)^{-1}\mathrm{e}^{-t/RC}$ for $t>0$ and $0$ for $t<0$ , since

\int_0^{\infty}\frac{1}{RC}\mathrm{e}^{-t/RC}\mathrm{e}^{-i\omega t}\,dt = \frac{1}{RC}\cdot\frac{1}{(1/RC)+i\omega} = \frac{1}{1+i\omega RC} ,

which both confirms the transfer function and shows the filter is causal, as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.5, Eqs. (20.105)-(20.106), p. 999) require of any physical device. With $R = 1\,\mathrm{k}\Omega$ and $C = 1\,\mu\mathrm{F}$ the half-power point is at $\omega = 10^{3}\,\mathrm{s}^{-1}$ , that is at an ordinary frequency of $159.15\ \mathrm{Hz}$ .

ApplicationOriginal

Anti-aliasing, and why it must be analogue. Every digitiser is preceded by a low-pass filter with its cutoff below the Nyquist frequency of the sampler. Part 4 says why no digital alternative exists: filtering after sampling multiplies the already-folded spectrum by $\lvert A\rvert^2$ and cannot unfold it. The $RC$ section of the example is the cheapest such filter, and its six decibels per octave is usually not enough; cascading $n$ sections gives $6n$ decibels per octave, with the transfer functions multiplying as in part 2.

Climate and geophysical time series. A monthly series is routinely passed through a low-pass before the decadal signal is examined, and through a band-pass to isolate, say, the two-to-seven-year band in which El Nino lives. Both the remarks above bite here. The negative side lobes of a badly designed filter put spurious anticorrelations into the output; and because smoothing manufactures autocorrelation, a significance test applied to a filtered series using the raw number of months will reject the null far too often. The standard practice of quoting an effective sample size $N\big/\bigl(1+2\sum_{m\ge1}\rho_Y(m)\bigr)$ is a direct consequence of part 4.

Matched filtering. When the signal shape $s$ is known and the noise is white, the filter maximising the output signal-to-noise ratio at the arrival time is the time-reverse of the signal, $a_k = s_{-k}$ , with transfer function $\overline{\hat s}$ . The output is then the cross-correlation of data with template, and the detection statistic is its peak. This is how a gravitational-wave chirp, a radar return or a seismic arrival is pulled out of noise; for coloured noise of density $g$ , the weights are those of $\overline{\hat s}/g$ , which is a filter that whitens first and then matches.

Instrument response and its removal. A seismometer, a thermometer in a housing and an amplifier each impose their own $\varphi(\omega)$ on what they measure, and recovering the true input means dividing by it. Part 2 makes this division trivial in principle; the deconvolution remark makes it delicate in practice, and the standard compromise is to divide only over the band where $\lvert\varphi\rvert$ is comfortably above the noise and to taper to zero outside it.

Method

Autoregressive models and red noise

Adapted

Add to PDF

Builds the AR(1) process from uncorrelated noise, derives its autocorrelation $\alpha^{|m|}$ and its power spectrum, obtains the Yule-Walker equations that fit an AR( $p$ ) scheme to data, and uses the fitted AR(1) spectrum as the red-noise null against which a spectral peak is tested, with a worked twelve-point series whose apparent cycle survives a white-noise test and fails a red-noise one.

Prerequisites: Stationary processes and spectral densities, Expectation, variance and covariance, Power series

Note(Why a null hypothesis has to be red)Original

A geophysical record almost never looks like white noise. Yesterday's temperature anomaly, last month's sea level, last year's sunspot count all carry over into the next value, because the system has inertia: an ocean mixed layer, a reservoir, a population. Persistence alone, with no oscillator anywhere in the system, already produces a spectrum that rises towards low frequency — more power in the slow fluctuations than in the fast ones. Such a spectrum is called red, by analogy with light.

That is a nuisance, because it is also what a long-period cycle looks like. A peak in an estimated spectrum is evidence of a cycle only if it stands above what mere persistence would have produced, and to say how much is "above" one needs a quantitative null hypothesis. The AR(1) process supplies it with a single fitted number. This unit builds the process, computes its spectrum exactly, and shows how to fit it and test against it. The example below is a series whose apparent twelve-point cycle is comfortably significant against white noise and not significant at all against red noise; that reversal is the whole point of the method.

Definition(Autoregressive schemes)AdaptedAdapted from Grimmett & Stirzaker (2020), §8.2, the definition of an autoregressive scheme, p. 410, and §9.3, Example (9.3.23), p. 444

Let $\{Z_n : n \in \mathbb{Z}\}$ be uncorrelated real random variables with $\mathbb{E}(Z_n) = 0$ and $\operatorname{var}(Z_n) = \sigma_Z^{2}$ , the same finite value for every $n$ . Such a sequence is called discrete white noise, or the innovations.

A real sequence $\{X_n : n \in \mathbb{Z}\}$ is an autoregressive scheme of order $p$ , written AR( $p$ ), if

X_n \;=\; \sum_{i=1}^{p}\alpha_i X_{n-i} \;+\; Z_n , \qquad n \in \mathbb{Z} ,

for real constants $\alpha_1,\dots,\alpha_p$ . The case $p = 1$ ,

X_n \;=\; \alpha X_{n-1} + Z_n ,

is the one this unit is mostly about. Grimmett and Stirzaker introduce these schemes, alongside moving averages $Y_n = \sum_{i=0}^{r}\alpha_i X_{n-i}$ , as the two families a statistician fits to an observed series; they take $\sigma_Z^{2} = 1$ throughout, and we keep $\sigma_Z^{2}$ explicit, which multiplies every covariance below by $\sigma_Z^{2}$ and leaves every correlation and every normalised spectral density unchanged.

Two conventions. Frequency $\lambda$ runs over $(-\pi,\pi]$ , one unit of $n$ being one sampling interval, as in Stationary processes and spectral densities; a frequency $\lambda$ corresponds to a period of $2\pi/\lambda$ sampling intervals, and $\lambda = \pi$ is the Nyquist frequency. And beside the normalised spectral density $f$ of that unit, with $\int_{-\pi}^{\pi} f = 1$ , we use the unnormalised power spectrum

P(\lambda) \;=\; \sum_{m=-\infty}^{\infty} c(m)\,e^{-im\lambda} \;=\; 2\pi c(0)\,f(\lambda) , \qquad \frac{1}{2\pi}\int_{-\pi}^{\pi} P(\lambda)\,d\lambda \;=\; c(0) ,

which is the discrete-time form of the physicists' $S(\omega)$ in the conventions table of that unit. $P$ is what a spectrum estimate is compared against, because it carries the units of variance.

Method(The stationary AR(1) process, its correlations and its spectrum)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.2, Example (9.2.5), p. 437, and §9.3, Example (9.3.23) including the general AR($r$) criterion, p. 444

Let $\lvert\alpha\rvert < 1$ . Then:

Existence and uniqueness. The series $X_n = \sum_{k=0}^{\infty}\alpha^{k}Z_{n-k}$ converges in mean square, and the sequence so defined is the unique weakly stationary solution of $X_n = \alpha X_{n-1}+Z_n$ .
Moments. $\mathbb{E}(X_n) = 0$ and $c(m) \;=\; \operatorname{cov}(X_n, X_{n+m}) \;=\; \frac{\sigma_Z^{2}\,\alpha^{\lvert m\rvert}}{1-\alpha^{2}} , \qquad \rho(m) \;=\; \alpha^{\lvert m\rvert} , \qquad m \in \mathbb{Z} .$ In particular $c(0) = \sigma_Z^{2}/(1-\alpha^{2})$ and $\alpha = \rho(1)$ .
Spectrum. The spectral density and power spectrum are $f(\lambda) \;=\; \frac{1-\alpha^{2}}{2\pi\bigl(1 - 2\alpha\cos\lambda + \alpha^{2}\bigr)} , \qquad P(\lambda) \;=\; \frac{\sigma_Z^{2}}{\bigl\lvert 1-\alpha e^{-i\lambda}\bigr\rvert^{2}} \;=\; \frac{\sigma_Z^{2}}{1-2\alpha\cos\lambda+\alpha^{2}} ,$ for $-\pi \le \lambda \le \pi$ . For $0 < \alpha < 1$ the spectrum decreases monotonically from $P(0) = \frac{\sigma_Z^{2}}{(1-\alpha)^{2}} = c(0)\,\frac{1+\alpha}{1-\alpha} \qquad\text{to}\qquad P(\pi) = \frac{\sigma_Z^{2}}{(1+\alpha)^{2}} = c(0)\,\frac{1-\alpha}{1+\alpha} ,$ a ratio of $\bigl((1+\alpha)/(1-\alpha)\bigr)^{2}$ between the two ends of the band: this is the red shape. For $-1<\alpha<0$ the same formula increases monotonically (a blue spectrum), and $\alpha = 0$ returns the flat spectrum $P \equiv \sigma_Z^{2}$ of white noise.
Prediction. The best linear predictor of $X_{r+k}$ , $k\ge1$ , from $X_r, X_{r-1},\dots,X_{r-s}$ is $\alpha^{k}X_r$ , with mean squared error $\sigma_Z^{2}(1-\alpha^{2k})/(1-\alpha^{2})$ . Only the most recent value is used, whatever $s$ is.

General order. Grimmett and Stirzaker state, for the AR( $p$ ) scheme with unit-variance innovations, that $X$ is stationary if the complex roots $\theta_1,\dots,\theta_p$ of

A(z) \;=\; z^{p} - \alpha_1 z^{p-1} - \cdots - \alpha_p

all satisfy $\lvert\theta_j\rvert < 1$ , and that the spectral density is then $f(\lambda) = 1/\bigl(2\pi\sigma^{2}\lvert A(e^{-i\lambda})\rvert^{2}\bigr)$ with $\sigma^{2} = \operatorname{var}(X_0)$ . With a general innovation variance this reads

P(\lambda) \;=\; \frac{\sigma_Z^{2}}{\bigl\lvert A(e^{-i\lambda})\bigr\rvert^{2}} \;=\; \frac{\sigma_Z^{2}}{\bigl\lvert 1 - \alpha_1e^{-i\lambda} - \cdots - \alpha_pe^{-ip\lambda}\bigr\rvert^{2}} ,

the two forms of the denominator agreeing because $\lvert e^{-ip\lambda}\rvert = 1$ and because the $\alpha_i$ are real, so that conjugating the bracket leaves its modulus unchanged. Part 3 is the case $p=1$ .

Proposition(The Yule-Walker equations)Original

Let $\{X_n\}$ be a weakly stationary AR( $p$ ) scheme with autocovariance function $c$ , driven by innovations $\{Z_n\}$ that are uncorrelated with zero mean and variance $\sigma_Z^{2}$ , and suppose each $Z_n$ is uncorrelated with every earlier value $X_{n-1}, X_{n-2},\dots$ of the process. Then

c(j) \;=\; \sum_{i=1}^{p}\alpha_i\,c(j-i) , \qquad j = 1, 2, \dots, p , \tag{YW}

and

\sigma_Z^{2} \;=\; c(0) - \sum_{i=1}^{p}\alpha_i\,c(i) .

Dividing (YW) by $c(0)$ gives the same system in the correlations, $\rho(j) = \sum_i \alpha_i\rho(j-i)$ , which in matrix form is

\begin{pmatrix} \rho(0) & \rho(1) & \cdots & \rho(p-1)\\ \rho(1) & \rho(0) & \cdots & \rho(p-2)\\ \vdots & & \ddots & \vdots\\ \rho(p-1) & \rho(p-2) & \cdots & \rho(0) \end{pmatrix} \begin{pmatrix}\alpha_1\\ \alpha_2\\ \vdots\\ \alpha_p\end{pmatrix} = \begin{pmatrix}\rho(1)\\ \rho(2)\\ \vdots\\ \rho(p)\end{pmatrix} .

The matrix is the symmetric Toeplitz matrix of correlations; it is non-singular whenever the first $p$ values of the process are not almost surely linearly dependent. For $p=1$ the system is the single equation $\alpha = \rho(1)$ .

These are the Yule-Walker equations. The name is standard in the time-series literature but is not used by our source; what Grimmett and Stirzaker give, as their prediction equations (9.2.3), is the system $\sum_{i=0}^{s}a_i c(\lvert i-j\rvert) = c(k+j)$ , $0\le j\le s$ , satisfied by the coefficients of the best linear predictor of $X_{r+k}$ from $X_r,\dots,X_{r-s}$ . Taking $k=1$ and $s=p-1$ there, and relabelling $a_{i-1}$ as $\alpha_i$ , reproduces (YW) exactly. So the Yule-Walker system has two readings: it fits an AR( $p$ ) model, and it finds the best one-step linear predictor from $p$ past values. They coincide because an AR( $p$ ) scheme is precisely a process whose one-step prediction error is uncorrelated with the past.

Idea of the proofOriginal

Why the AR(1) process is a geometric smoothing of the noise. Substituting the recursion into itself $k$ times gives $X_n = \alpha^{k}X_{n-k} + \sum_{j<k}\alpha^{j}Z_{n-j}$ . If $\lvert\alpha\rvert<1$ the memory term $\alpha^{k}X_{n-k}$ dies away, and what is left is a geometrically weighted sum of all past shocks. Every statement about the process is then a geometric series: the variance is $\sigma_Z^{2}\sum\alpha^{2k}$ , the lag- $m$ covariance is the same sum shifted, and the spectrum is the squared modulus of the geometric series $\sum\alpha^{k}e^{-ik\lambda} = (1-\alpha e^{-i\lambda})^{-1}$ .

Why the spectrum is red. That last expression is largest when $\alpha e^{-i\lambda}$ points along the positive real axis, which for $\alpha>0$ happens at $\lambda = 0$ , and smallest when it points the other way, at $\lambda=\pi$ . Positive $\alpha$ therefore piles power at low frequency. Geometrically, $\lvert 1-\alpha e^{-i\lambda}\rvert$ is the distance from the point $\alpha e^{-i\lambda}$ on a circle of radius $\alpha$ to the point $1$ , and that distance is smallest when the circle point is nearest $1$ .

Why the Yule-Walker equations are an orthogonality statement. The defining feature of an AR( $p$ ) scheme is that the residual $Z_n = X_n - \sum\alpha_iX_{n-i}$ carries no linear information about the past. Multiplying the recursion by a past value $X_{n-j}$ and taking expectations records exactly that: the residual term drops out, and what survives is a linear relation among covariances. This is the projection theorem of Conditional expectation in coordinates, with "orthogonal" meaning "uncorrelated".

Why red noise is the right null. A test needs a null that is sceptical in the right direction. Against white noise, any persistent record shows a spuriously significant low-frequency peak, because the null underestimates the low-frequency power that persistence alone supplies. The AR(1) null absorbs persistence into one parameter, estimated from the same record, and asks whether a peak beats what that persistence would have delivered. It is the cheapest null that is not obviously wrong.

Proof(The AR(1) process: existence, moments and spectrum)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.2, Example (9.2.5), p. 437, and §9.3, Example (9.3.23), p. 444

Throughout $\lvert\alpha\rvert<1$ .

Step 1: the series converges. Put $S_K = \sum_{k=0}^{K}\alpha^{k}Z_{n-k}$ . For $L > K$ , the increments $S_L - S_K$ involve distinct, hence uncorrelated, innovations, so

\mathbb{E}\bigl((S_L-S_K)^{2}\bigr) = \sigma_Z^{2}\sum_{k=K+1}^{L}\alpha^{2k} \;\le\; \sigma_Z^{2}\,\frac{\alpha^{2(K+1)}}{1-\alpha^{2}} \;\xrightarrow[K\to\infty]{}\; 0 ,

using $\lvert\alpha\rvert<1$ to sum the geometric series. So $\{S_K\}$ is Cauchy in mean square and converges to a limit $X_n$ with $\mathbb{E}(X_n)=0$ .

Step 2: it solves the recursion. Splitting off the $k=0$ term, $S_K = Z_n + \alpha\sum_{k=0}^{K-1}\alpha^{k}Z_{n-1-k}$ . Letting $K\to\infty$ in mean square gives $X_n = Z_n + \alpha X_{n-1}$ .

Step 3: the covariances. For $m \ge 0$ , the products $Z_{n-j}Z_{n+m-k}$ have zero expectation unless $n-j = n+m-k$ , that is $k = j+m$ ; the surviving terms contribute $\sigma_Z^{2}$ each. Hence

c(m) = \mathbb{E}(X_nX_{n+m}) = \sum_{j,k\ge0}\alpha^{j+k}\,\mathbb{E}(Z_{n-j}Z_{n+m-k}) = \sigma_Z^{2}\sum_{j\ge0}\alpha^{\,j}\alpha^{\,j+m} = \frac{\sigma_Z^{2}\alpha^{m}}{1-\alpha^{2}} .

This does not depend on $n$ , so the process is weakly stationary; and $c(-m)=c(m)$ because the process is real, giving $c(m) = \sigma_Z^{2}\alpha^{\lvert m\rvert}/(1-\alpha^{2})$ and $\rho(m)=\alpha^{\lvert m\rvert}$ .

Step 4: uniqueness. Let $Y$ be any weakly stationary solution and set $D_n = Y_n - X_n$ . Subtracting the two recursions gives $D_n = \alpha D_{n-1}$ , hence $D_n = \alpha^{k}D_{n-k}$ for every $k\ge1$ . Both $X$ and $Y$ are stationary, so $\mathbb{E}(D_n^{2})$ is a constant $\delta$ independent of $n$ , and $\delta = \alpha^{2k}\delta$ for all $k$ . Since $\lvert\alpha\rvert<1$ , letting $k\to\infty$ forces $\delta = 0$ , so $Y_n = X_n$ almost surely.

Step 5: the spectrum. Since $\sum_m\lvert\rho(m)\rvert = 1 + 2\sum_{m\ge1}\lvert\alpha\rvert^{m} < \infty$ , the proposition on absolutely summable autocorrelations in Stationary processes and spectral densities applies, and $P(\lambda) = \sum_m c(m)e^{-im\lambda}$ . Split the sum at $m=0$ and sum two geometric series:

\frac{P(\lambda)}{c(0)} = 1 + \sum_{m\ge1}\alpha^{m}e^{-im\lambda} + \sum_{m\ge1}\alpha^{m}e^{im\lambda} = 1 + \frac{\alpha e^{-i\lambda}}{1-\alpha e^{-i\lambda}} + \frac{\alpha e^{i\lambda}}{1-\alpha e^{i\lambda}} .

Put the three terms over the common denominator $\lvert 1-\alpha e^{-i\lambda}\rvert^{2} = 1-2\alpha\cos\lambda+\alpha^{2}$ . The two numerators from the fractions are $\alpha e^{-i\lambda}(1-\alpha e^{i\lambda}) + \alpha e^{i\lambda}(1-\alpha e^{-i\lambda}) = 2\alpha\cos\lambda - 2\alpha^{2}$ , so

\frac{P(\lambda)}{c(0)} = \frac{\bigl(1-2\alpha\cos\lambda+\alpha^{2}\bigr) + \bigl(2\alpha\cos\lambda-2\alpha^{2}\bigr)}{1-2\alpha\cos\lambda+\alpha^{2}} = \frac{1-\alpha^{2}}{1-2\alpha\cos\lambda+\alpha^{2}} .

Multiplying by $c(0) = \sigma_Z^{2}/(1-\alpha^{2})$ gives $P(\lambda) = \sigma_Z^{2}/(1-2\alpha\cos\lambda+\alpha^{2})$ , and dividing $P$ by $2\pi c(0)$ gives $f$ . Monotonicity in $\lambda$ on $[0,\pi]$ is immediate, since $\cos\lambda$ decreases there and the denominator is $1+\alpha^{2}-2\alpha\cos\lambda$ .

Step 6: prediction. By Grimmett and Stirzaker's Theorem (9.2.1) the coefficients $a_0,\dots,a_s$ of the best linear predictor of $X_{r+k}$ satisfy $\sum_{i}a_ic(\lvert i-j\rvert) = c(k+j)$ for $0\le j\le s$ . Try $a_0 = \alpha^{k}$ and $a_i = 0$ otherwise: the left side becomes $\alpha^{k}c(j) = \sigma_Z^{2}\alpha^{k+j}/(1-\alpha^{2})$ , which is $c(k+j)$ because $k+j\ge0$ . So the trial solves the system and $\widehat X_{r+k} = \alpha^{k}X_r$ . Its error is

\operatorname{var}(X_{r+k} - \alpha^{k}X_r) = c(0) - 2\alpha^{k}c(k) + \alpha^{2k}c(0) = \frac{\sigma_Z^{2}}{1-\alpha^{2}}\bigl(1 - 2\alpha^{2k} + \alpha^{2k}\bigr) = \frac{\sigma_Z^{2}\bigl(1-\alpha^{2k}\bigr)}{1-\alpha^{2}} . \qquad\square

Proof(The Yule-Walker equations)Original

Assume the hypotheses of the proposition, so $\mathbb{E}(X_n) = 0$ and $c(j) = \mathbb{E}(X_nX_{n+j})$ .

Step 1: multiply by a past value. Fix $j$ with $1\le j\le p$ . Multiply $X_n = \sum_{i=1}^{p}\alpha_iX_{n-i} + Z_n$ by $X_{n-j}$ and take expectations:

\mathbb{E}(X_nX_{n-j}) = \sum_{i=1}^{p}\alpha_i\,\mathbb{E}(X_{n-i}X_{n-j}) + \mathbb{E}(Z_nX_{n-j}) .

Step 2: kill the noise term. Here the hypothesis that $Z_n$ is uncorrelated with $X_{n-1},X_{n-2},\dots$ is used: since $j\ge1$ and both factors have zero mean, $\mathbb{E}(Z_nX_{n-j}) = 0$ .

Step 3: read off the covariances. By stationarity $\mathbb{E}(X_nX_{n-j}) = c(-j) = c(j)$ and $\mathbb{E}(X_{n-i}X_{n-j}) = c(j-i)$ , which gives (YW).

Step 4: the variance equation. Now take $j=0$ , that is multiply by $X_n$ itself:

c(0) = \sum_{i=1}^{p}\alpha_i c(i) + \mathbb{E}(Z_nX_n) .

Substituting the recursion once more inside the last expectation,

\mathbb{E}(Z_nX_n) = \sum_{i=1}^{p}\alpha_i\,\mathbb{E}(Z_nX_{n-i}) + \mathbb{E}(Z_n^{2}) = 0 + \sigma_Z^{2} ,

again by step 2 and by $\operatorname{var}(Z_n)=\sigma_Z^{2}$ . Rearranging gives the stated expression for $\sigma_Z^{2}$ . $\square$

For $p=1$ the two equations are $c(1) = \alpha c(0)$ and $\sigma_Z^{2} = c(0)-\alpha c(1) = c(0)(1-\alpha^{2})$ , consistent with part 2 of the statement above.

Remark(Where the hypotheses are used, and what fails without them)Original

$\lvert\alpha\rvert<1$ is used in step 1 (convergence), step 4 (uniqueness) and step 5 (summing the series) of the first proof, and it cannot be dropped.

At $\alpha = 1$ the recursion is a random walk (Random walks). It has no stationary solution: starting from $X_0 = 0$ , $\operatorname{var}(X_n) = n\sigma_Z^{2}$ grows without bound. Formally the spectrum $\sigma_Z^{2}/(2-2\cos\lambda) \sim \sigma_Z^{2}/\lambda^{2}$ diverges at the origin and is not integrable, which is the frequency-domain image of the infinite variance. A red-noise test is therefore meaningless on an undifferenced random walk, and a fitted $\hat\alpha$ very close to $1$ is a warning that one is looking at a trend rather than at a stationary process.
At $\lvert\alpha\rvert>1$ the backwards series diverges, but there is still a stationary solution, the non-causal one $X_n = -\sum_{k\ge1}\alpha^{-k}Z_{n+k}$ built from future innovations; substituting it into the recursion verifies it directly. Its autocorrelation is $\alpha^{-\lvert m\rvert}$ , the same as that of an AR(1) with parameter $1/\alpha$ . So correlations alone cannot distinguish $\alpha$ from $1/\alpha$ , and the condition $\lvert\alpha\rvert<1$ is what picks the one solution that depends only on the past. This paragraph is our own observation, not something our source states.

Uncorrelated innovations with constant variance are used in step 3 of the first proof (only matching indices survive) and in step 2 of the second (the noise term vanishes). If the innovations are themselves correlated, (YW) acquires extra terms and fitting an AR( $p$ ) by Yule-Walker returns biased coefficients; the usual symptom is that the residuals of the fit are not white, which is why a residual autocorrelation check belongs in the recipe below.

$Z_n$ uncorrelated with the past of $X$ , not with all of it. This is automatic for the causal solution of step 1, since $X_{n-j}$ is built from $Z_{n-j},Z_{n-j-1},\dots$ . It is false for the non-causal solution above, and that is exactly why Yule-Walker estimation always returns a causal model.

Uncorrelated is weaker than independent. Nothing above uses independence or normality, so the autocovariance and spectrum hold for any white-noise innovations. The red-noise test below is a different matter: its distributional statement needs the process to be Gaussian, and the chi-squared law for a periodogram ordinate is in any case only asymptotic in the record length.

Conventions to check before comparing numbers. Three are easy to get wrong. Some authors write the AR(1) recursion with a minus sign, $X_n + \phi X_{n-1} = Z_n$ , which flips the sign of the correlation. Some define the spectrum over $[0,\pi]$ only, doubling it. And the "lag-1 autocorrelation" quoted with a record is sometimes computed dividing by $N$ and sometimes by $N-m$ ; the two differ by a factor $1-m/N$ , which matters at the long lags and is why the $N$ divisor, used below, is preferred — it is the one that always yields a non-negative definite estimated covariance sequence.

Method(Fitting an AR(1) null and testing a spectral peak against it)Original

Given a record $x_1,\dots,x_N$ sampled at equal intervals:

Remove the mean, and any trend you do not want tested. Put $\bar x = N^{-1}\sum_n x_n$ and $d_n = x_n - \bar x$ . A linear trend left in the record inflates the lowest frequencies and will inflate $\hat\alpha$ with them; if a trend is physically expected, fit and subtract it first, and say so.
Estimate the autocovariances at lags $0$ and $1$ : $\hat c_0 = \frac1N\sum_{n=1}^{N}d_n^{2} , \qquad \hat c_1 = \frac1N\sum_{n=1}^{N-1}d_nd_{n+1} .$
Solve Yule-Walker. For $p=1$ this is just $\hat\alpha = \hat c_1/\hat c_0$ , and the innovation variance is $\hat\sigma_Z^{2} = \hat c_0\bigl(1-\hat\alpha^{2}\bigr)$ . For order $p$ , solve the Toeplitz system of the proposition, for example by the Levinson-Durbin recursion.
Check the fit. Require $\lvert\hat\alpha\rvert<1$ (for order $p$ , that the roots of $A$ lie inside the unit circle). Compare $\hat\rho(2)$ with $\hat\alpha^{2}$ : a serious mismatch means AR(1) is the wrong null and a higher order, or a different model, is needed.
Form the null spectrum. $P_{\text{null}}(\lambda) \;=\; \frac{\hat\sigma_Z^{2}}{1 - 2\hat\alpha\cos\lambda + \hat\alpha^{2}} , \qquad \lambda \in [0,\pi] .$
Compute the spectral estimate on the same scale. At the Fourier frequencies $\lambda_j = 2\pi j/N$ use the periodogram $I(\lambda_j) \;=\; \frac1N\Bigl\lvert\sum_{n=1}^{N}d_n e^{-i\lambda_j n}\Bigr\rvert^{2} ,$ whose properties, and whose smoothed and windowed variants, are the subject of The periodogram and spectral estimation. The normalisation is fixed by Parseval: $N^{-1}\sum_{j=0}^{N-1}I(\lambda_j) = \hat c_0$ , the same total variance that $\frac{1}{2\pi}\int P$ carries.
Set the threshold. For a Gaussian stationary process, and $\lambda_j$ neither $0$ nor $\pi$ , the ratio $2I(\lambda_j)/P(\lambda_j)$ is approximately $\chi^{2}$ with two degrees of freedom. A one-sided test at level $q$ therefore rejects the null when $I(\lambda_j) \;>\; \tfrac12\,\chi^{2}_{2,\,q}\;P_{\text{null}}(\lambda_j) , \qquad \tfrac12\chi^{2}_{2,\,0.95} = 2.995732 , \quad \tfrac12\chi^{2}_{2,\,0.99} = 4.605170 .$ (If the estimate is a smoothed or Welch-averaged one with $\nu$ effective degrees of freedom, replace this factor by $\chi^{2}_{\nu,q}/\nu$ ; see Hypothesis testing for the vocabulary.)
Allow for multiple testing. You are usually scanning all $\lfloor N/2\rfloor$ frequencies at once, so at the nominal $5\%$ level about one ordinate in twenty will exceed the threshold under the null. Either nominate the frequency of interest in advance, or apply a correction: the probability that no ordinate among $M$ roughly independent ones exceeds the level- $q$ threshold is about $q^{M}$ , so testing at $q = 0.95^{1/M}$ controls the chance of any false peak at $5\%$ .
Quote the effective sample size. For other statistics computed from the same record, persistence costs degrees of freedom. Since $\sum_{m}\rho(m) = (1+\alpha)/(1-\alpha)$ for AR(1), the variance of the sample mean is approximately $\hat c_0(1+\hat\alpha)/\{N(1-\hat\alpha)\}$ , so the record is worth $N_{\text{eff}} \;=\; N\,\frac{1-\hat\alpha}{1+\hat\alpha}$ independent observations. Report it: it is usually smaller than people expect.

Example(A twelve-point record: significant against white noise, not against red)Original

A twelve-year index, one value per year, reads

x = (11,\,10,\,13,\,13,\,13,\,10,\,7,\,8,\,6,\,9,\,9,\,11) .

The eye sees a cycle: high for the first five years, low for the next four, recovering at the end — one oscillation in twelve years. Is it real?

Step 1: mean and deviations. The twelve values sum to $120$ , so $\bar x = 10$ exactly and

d = (1,\,0,\,3,\,3,\,3,\,0,\,-3,\,-2,\,-4,\,-1,\,-1,\,1) , \qquad \sum_n d_n = 0 .

Step 2: autocovariances. Summing squares, $\sum d_n^{2} = 1+0+9+9+9+0+9+4+16+1+1+1 = 60$ , and summing the eleven lagged products,

\sum_{n=1}^{11} d_nd_{n+1} = 0+0+9+9+0+0+6+8+4+1-1 = 36 .

Hence $\hat c_0 = 60/12 = 5$ and $\hat c_1 = 36/12 = 3$ .

Step 3: the fitted AR(1).

\hat\alpha = \frac{\hat c_1}{\hat c_0} = \frac{3}{5} = 0.6 , \qquad \hat\sigma_Z^{2} = \hat c_0\bigl(1-\hat\alpha^{2}\bigr) = 5\times 0.64 = 3.2 .

A check on the order: an AR(1) predicts $\rho(2) = \alpha^{2} = 0.36$ , while the lag-2 sample correlation is $\hat c_2/\hat c_0$ with $\sum_{n=1}^{10}d_nd_{n+2} = 3+0+9+0-9+0+12+2+4-1 = 20$ , giving $\hat c_2 = 20/12 = 5/3$ and $\hat\rho(2) = (5/3)/5 = 1/3 = 0.3333$ . That is within a few hundredths of $0.36$ , so AR(1) is an acceptable null.

Step 4: the null spectrum. With $\hat\alpha = 0.6$ ,

P_{\text{null}}(\lambda) = \frac{3.2}{1.36 - 1.2\cos\lambda} , \qquad P_{\text{null}}(0) = \frac{3.2}{0.16} = 20 , \qquad P_{\text{null}}(\pi) = \frac{3.2}{2.56} = 1.25 .

The spectrum is red by a factor $16 = (1.6/0.4)^{2}$ across the band, as part 3 of the statement predicts.

Step 5: the periodogram at the frequency of interest. The suspected cycle has period $12$ , which is $j=1$ , $\lambda_1 = \pi/6$ . Computing the two sums,

\sum_{n=1}^{12} d_n\cos(n\pi/6) = 0 , \qquad \sum_{n=1}^{12} d_n\sin(n\pi/6) = 11 + 3\sqrt3 ,

the first vanishing by an exact cancellation. Hence

I(\lambda_1) = \frac{\bigl(11+3\sqrt3\bigr)^{2}}{12} = \frac{37}{3} + \frac{11\sqrt3}{2} = 21.859613 .

As an arithmetic check, the twelve periodogram ordinates average to $\hat c_0$ : $N^{-1}\sum_j I(\lambda_j) = 5$ , which they do.

Step 6: the two tests.

Against white noise. The white null has $P_{\text{white}}\equiv\hat c_0 = 5$ at every frequency, so the $95\%$ threshold is $2.995732\times 5 = 14.978661$ . Since $21.859613 > 14.978661$ , the peak is significant at $5\%$ : on a white-noise null one would announce a twelve-year cycle.

Against red noise. Here

P_{\text{null}}(\pi/6) = \frac{3.2}{1.36 - 0.6\sqrt3} = 9.976010 ,

so the $95\%$ threshold is $2.995732\times 9.976010 = 29.885456$ . Now $21.859613 < 29.885456$ : the peak is not significant. Its ratio to the null is $21.859613/9.976010 = 2.191$ , which under the null would be exceeded with probability $e^{-2.191} = 0.112$ — about one chance in nine, before any correction for having scanned the five interior frequencies of a twelve-point record.

Step 7: the effective sample size. Finally,

N_{\text{eff}} = N\,\frac{1-\hat\alpha}{1+\hat\alpha} = 12\times\frac{0.4}{1.6} = 3 .

Twelve years of this index carry the information of three independent observations. That single number explains the whole result: one is not going to establish a cycle from three effective points, and the red-noise test is the formal statement of that intuition.

Application(Where red-noise tests are used)Original

Climate variability. Almost every climate index — the North Atlantic Oscillation, the Southern Oscillation Index, a basin-mean sea surface temperature, a borehole temperature profile — is first characterised by its lag-1 autocorrelation and its AR(1) spectrum, and claims of periodicity in them are tested against that spectrum rather than against white noise. The physical justification is Hasselmann's model of climate variability: a slow component (the ocean mixed layer) integrates fast weather forcing, and the resulting first-order stochastic equation is exactly the continuous-time limit of AR(1), the Ornstein-Uhlenbeck process of Stochastic differential equations, whose spectrum is the Lorentzian that the discrete formula above approximates for $\lambda\ll1$ . So the null is not an arbitrary convenience: it is the spectrum of a system with one reservoir and no oscillator, which is the simplest hypothesis a claimed climate cycle has to beat.

Degrees of freedom in correlation and trend tests. The step-9 figure $N_{\text{eff}} = N(1-\alpha)/(1+\alpha)$ is used far more often than the spectral test. A correlation between two autocorrelated records, or the significance of a fitted trend, must be assessed with $N_{\text{eff}}$ in place of $N$ ; using $N$ is the commonest way to manufacture a significant result from two unrelated red series. For $\alpha = 0.9$ the factor is $1/19$ , so a century of monthly data can be worth only a handful of independent numbers.

Instrument and detector noise. Amplifier drift, detector dark current and laser frequency noise all show spectra rising at low frequency, and an AR(1) fit gives a quick corner frequency: the spectrum of the fitted process falls to half its zero-frequency value at

\cos\lambda_{1/2} \;=\; \frac{1+\alpha^{2}-2(1-\alpha)^{2}}{2\alpha} ,

which for $\alpha$ near $1$ gives $\lambda_{1/2}\approx 1-\alpha$ (at $\alpha = 0.9$ the exact value is $0.105430$ against the estimate $0.1$ ), an inverse correlation time in units of the sampling interval. Deciding how long to integrate before a reference measurement is, in effect, choosing a sampling rate against that corner.

Prefiltering and prewhitening. Part 4 of the statement says the best one-step linear predictor of an AR(1) is $\alpha X_{n-1}$ , so the residual series $x_n - \hat\alpha x_{n-1}$ is an estimate of the innovations and should be white. "Prewhitening" a record this way before a spectral analysis makes the null flat and the test easier to read, at the cost of also flattening any genuine low-frequency signal; the transfer function $1-\hat\alpha e^{-i\lambda}$ that does it is a high-pass filter in the sense of Linear filters. Whether to prewhiten, or to test against a red null without filtering, is a question of which distortion you prefer; the two are algebraically the same test when the null is correct.

Model selection in econometrics and geophysics. The Yule-Walker equations are the entry point to the whole ARMA apparatus: fit increasing $p$ , watch the estimated innovation variance $\hat\sigma_Z^{2}$ stop falling, and penalise the extra parameters with an information criterion. For many physical records $p=1$ is enough, which is why the AR(1) null is the default and not merely the first case.

Method · advanced

The Kalman filter

Adapted

Add to PDF

Sets up the linear Gaussian state-space model, proves by induction on the Gaussian conditioning formula that the filtering distribution stays normal and obeys a forecast step and an update step with gain $K_k = P_k^{f}H_k^{\mathsf T}S_k^{-1}$ , identifies the update with a Bayesian posterior and with weighted least squares, and runs the recursion on a four-point scalar record where the forecast variance falls $1.875 \to 1.0380 \to 1.0023 \to 1.00015$ towards the fixed point $1$ and the final estimate $-0.000376 \pm 0.7071$ reproduces the exact nine-dimensional Gaussian conditional.

Prerequisites: The multivariate normal distribution, Least-squares fitting, Conditional expectation, Discrete-time Markov chains, Bayesian inference, Autoregressive models and red noise, Quadratic forms and the Rayleigh quotient, Linear systems and rank

Note(What this unit is for)Original

A satellite is tracked, a reactor's core temperature is inferred from thermocouples, a tide gauge is read every hour: in each case a system evolves according to a model that is only approximately right, and is measured by instruments that are only approximately accurate. Neither the model nor the data alone answers the question "where is the system now?".

The Kalman filter answers it, and answers it recursively. At every step it carries one mean vector and one covariance matrix, and the arrival of a new observation updates both in closed form. Nothing earlier than the current summary is ever revisited, so the cost per observation does not grow with the length of the record. The reason this is possible is a single structural fact: within the linear Gaussian model, the conditional distribution of the state given all data so far is normal, and a normal distribution is exactly described by a mean and a covariance.

Everything here is a consequence of two results already proved in The multivariate normal distribution: an affine image of a normal vector is normal, and the conditional law of one block of a normal vector given another is normal with a mean affine in the observed block and a covariance that does not depend on it. The first gives the forecast step, the second gives the update step, and induction does the rest.

Definition(The linear Gaussian state-space model)Original

Fix integers $n, p \ge 1$ . A linear Gaussian state-space model consists of an unobserved state sequence $x_0, x_1, x_2, \dots$ in $\mathbb{R}^{n}$ and an observed sequence $y_1, y_2, \dots$ in $\mathbb{R}^{p}$ satisfying, for $k \ge 1$ ,

x_k \;=\; F_k\,x_{k-1} + w_k , \qquad y_k \;=\; H_k\,x_k + v_k ,

where

$F_k$ (the transition matrix, $n \times n$ ) and $H_k$ (the observation matrix, $p \times n$ ) are known and non-random;
the process noise $w_k \sim N(0, Q_k)$ on $\mathbb{R}^{n}$ and the observation noise $v_k \sim N(0, R_k)$ on $\mathbb{R}^{p}$ have known covariances, with each $R_k$ symmetric positive definite and each $Q_k$ symmetric positive semi-definite;
the initial state is $x_0 \sim N(m_0, P_0)$ with $m_0$ and $P_0$ known, $P_0$ symmetric positive semi-definite;
the whole collection $x_0, w_1, w_2, \dots, v_1, v_2, \dots$ is independent.

Normal vectors are understood in the sense of The multivariate normal distribution, the singular case included: $N(\mu, V)$ with $V$ singular means a vector every linear combination of whose components is univariate normal, with mean vector $\mu$ and covariance matrix $V$ . This is why $Q_k$ and $P_0$ are allowed to be semi-definite — a noiseless component of the dynamics is a legitimate model — while $R_k$ must be definite, because an observation with no error at all would determine a linear combination of the state exactly and is excluded here.

Write $Y_k = (y_1, \dots, y_k)$ for the data up to time $k$ , with $Y_0$ empty. The three distributions one wants are

\underbrace{p(x_k \mid Y_{k-1})}_{\text{forecast}}, \qquad \underbrace{p(x_k \mid Y_{k})}_{\text{filtering}}, \qquad \underbrace{p(y_k \mid Y_{k-1})}_{\text{predictive}} .

Computing the filtering distribution sequentially in $k$ is filtering; computing $p(x_j \mid Y_k)$ for $j < k$ is smoothing and is not treated here.

The state sequence is a Markov chain in the sense of Discrete-time Markov chains with a continuous state space: $x_k$ depends on the past only through $x_{k-1}$ . The pair (hidden Markov state, noisy observation) is the structure the filter exploits.

Method(The Kalman filter)Original

In the model above, for every $k \ge 1$ :

(i) The filtering distribution is normal. Conditional on $Y_k = (y_1,\dots,y_k)$ ,

x_k \;\sim\; N\bigl(m_k,\, P_k\bigr) ,

where $m_k$ depends on the data and $P_k$ does not: $P_k$ is a deterministic function of $F_1,\dots,F_k$ , $H_1,\dots,H_k$ , $Q_1,\dots,Q_k$ , $R_1,\dots,R_k$ and $P_0$ alone.

(ii) The recursion. Starting from $m_0, P_0$ , the pairs $(m_k, P_k)$ are generated by two steps.

Forecast (time update).

m_k^{f} = F_k\,m_{k-1} , \qquad P_k^{f} = F_k\,P_{k-1}\,F_k^{\mathsf T} + Q_k .

Conditional on $Y_{k-1}$ , $x_k \sim N(m_k^{f}, P_k^{f})$ .

Update (measurement update). Form the innovation and its covariance,

d_k = y_k - H_k\,m_k^{f} , \qquad S_k = H_k\,P_k^{f}H_k^{\mathsf T} + R_k ,

then the Kalman gain

K_k \;=\; P_k^{f}\,H_k^{\mathsf T}\,S_k^{-1} ,

and set

m_k \;=\; m_k^{f} + K_k\,d_k , \qquad P_k \;=\; \bigl(I - K_kH_k\bigr)P_k^{f} .

The matrix $S_k$ is invertible, so $K_k$ is well defined.

(iii) The predictive distribution. Conditional on $Y_{k-1}$ , $y_k \sim N(H_km_k^{f}, S_k)$ . Hence the innovations $d_1, d_2, \dots$ are independent, with $d_k \sim N(0, S_k)$ unconditionally, and the log-likelihood of the whole record is

\log p(y_1,\dots,y_N) \;=\; -\tfrac12\sum_{k=1}^{N}\Bigl( d_k^{\mathsf T}S_k^{-1}d_k + \log\det S_k + p\log 2\pi \Bigr) .

(iv) Equivalent forms of the update. If in addition $P_k^{f}$ is invertible, then $P_k$ is invertible and

P_k^{-1} \;=\; \bigl(P_k^{f}\bigr)^{-1} + H_k^{\mathsf T}R_k^{-1}H_k , \qquad K_k \;=\; P_k\,H_k^{\mathsf T}R_k^{-1} ,

the information form: precisions add. Whether or not $P_k^{f}$ is invertible,

P_k \;=\; \bigl(I-K_kH_k\bigr)P_k^{f}\bigl(I-K_kH_k\bigr)^{\mathsf T} + K_kR_kK_k^{\mathsf T} ,

the Joseph form, which is manifestly symmetric and positive semi-definite.

Idea of the proofOriginal

Only two facts about normal vectors are needed, and the filter is the alternation of them.

Affine images of normal vectors are normal. So if the state at time $k-1$ is normal given the data so far, then $x_k = F_kx_{k-1} + w_k$ — an affine function of that state and of fresh noise independent of everything already seen — is normal too, and its mean and covariance follow from the usual rules. That is the forecast step; no new idea is involved.

Conditioning a normal vector on part of itself gives a normal vector. So if, given the data so far, the pair $(x_k, y_k)$ is jointly normal, then conditioning on the newly observed $y_k$ leaves $x_k$ normal, with mean shifted by $\operatorname{cov}(x_k,y_k)\operatorname{var}(y_k)^{-1}$ times the surprise $y_k - \mathbb{E}(y_k)$ and with covariance reduced by the Schur complement, which is where $K_k$ comes from. That is the update step.

The only thing to watch is that both steps must be carried out conditionally on the earlier data, and that this is legitimate: conditioning on $Y_{k-1}$ and then on $y_k$ is the same as conditioning on $Y_k$ at one stroke. The induction then runs, and the reason the covariance recursion decouples from the data is that the conditional covariance in the Gaussian conditioning formula does not involve the value conditioned on.

Proof(The filter, by induction on

k

)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, Theorem (4.9.6) and Definition (4.9.7), p. 133, and Exercises (4.9.7) and (4.9.9), pp. 134-135

Throughout, "normal" includes the singular case, in the linear-combination sense of Definition (4.9.7) of Grimmett and Stirzaker, as in The multivariate normal distribution.

Step 0: the induction hypothesis. Let $\mathcal{H}_k$ be the statement: conditional on $Y_k = (y_1,\dots,y_k)$ , the state $x_k$ is $N(m_k, P_k)$ , with $P_k$ non-random. For $k=0$ this is the assumption $x_0 \sim N(m_0,P_0)$ , since $Y_0$ is empty. Assume $\mathcal{H}_{k-1}$ and fix an observed value $Y_{k-1} = \eta$ .

Step 1: the new noises are untouched by the old data. Write $\mathcal{G}_{k-1}$ for the $\sigma$ -field generated by $x_0, w_1,\dots,w_{k-1}, v_1,\dots,v_{k-1}$ . By the model equations every $x_j$ with $j \le k-1$ and every $y_j$ with $j \le k-1$ is a function of these variables, so both $x_{k-1}$ and $Y_{k-1}$ are $\mathcal{G}_{k-1}$ -measurable. The independence assumption says that the pair $(w_k, v_k)$ is independent of $\mathcal{G}_{k-1}$ . Hence, conditionally on $Y_{k-1}=\eta$ :

$x_{k-1} \sim N(m_{k-1}, P_{k-1})$ , by $\mathcal{H}_{k-1}$ (here $m_{k-1}$ is the number obtained from the observed $\eta$ );
$w_k \sim N(0,Q_k)$ and $v_k \sim N(0,R_k)$ , with the same laws as unconditionally;
$x_{k-1}$ , $w_k$ and $v_k$ are mutually independent.

Only the independence hypothesis of the model is used here; without it the conditional law of $w_k$ would depend on $\eta$ and nothing below would survive.

Step 2: the forecast. Conditionally on $Y_{k-1}=\eta$ , the vector $(x_{k-1}, w_k, v_k)$ is normal with block-diagonal covariance $\operatorname{diag}(P_{k-1}, Q_k, R_k)$ : each block is normal and the three are independent, so every linear combination of all the components is a sum of three independent univariate normals and is therefore univariate normal. Now

\begin{pmatrix} x_k \\ y_k \end{pmatrix} = \begin{pmatrix} F_k & I & 0 \\ H_kF_k & H_k & I \end{pmatrix} \begin{pmatrix} x_{k-1} \\ w_k \\ v_k\end{pmatrix}

is an affine (here linear) image of it, substituting $y_k = H_k(F_kx_{k-1}+w_k) + v_k$ . By Theorem (4.9.6) of Grimmett and Stirzaker in the form stated in The multivariate normal distribution, the pair $(x_k,y_k)$ is jointly normal conditionally on $Y_{k-1}=\eta$ . Its moments come from linearity of expectation and the covariance law $V \mapsto AVA^{\mathsf T}$ :

\mathbb{E}(x_k \mid Y_{k-1}) = F_km_{k-1} =: m_k^{f} , \qquad \operatorname{var}(x_k\mid Y_{k-1}) = F_kP_{k-1}F_k^{\mathsf T} + Q_k =: P_k^{f} ,

the cross-term vanishing because $x_{k-1}$ and $w_k$ are conditionally independent. Likewise

\mathbb{E}(y_k\mid Y_{k-1}) = H_km_k^{f} , \qquad \operatorname{var}(y_k\mid Y_{k-1}) = H_kP_k^{f}H_k^{\mathsf T} + R_k = S_k ,

\operatorname{cov}(x_k, y_k \mid Y_{k-1}) = \operatorname{cov}\bigl(x_k,\, H_kx_k + v_k\bigr) = P_k^{f}H_k^{\mathsf T} ,

again because $v_k$ is conditionally independent of $x_k$ . This proves the forecast step and part (iii) of the statement, the predictive law being the second block of a jointly normal pair and hence normal by the marginal rule.

Step 3: $S_k$ is invertible. For $u \in \mathbb{R}^{p}$ with $u \ne 0$ ,

u^{\mathsf T}S_ku \;=\; \bigl(H_k^{\mathsf T}u\bigr)^{\mathsf T}P_k^{f}\bigl(H_k^{\mathsf T}u\bigr) + u^{\mathsf T}R_ku \;>\; 0 ,

the first term being $\ge 0$ because $P_k^{f} = F_kP_{k-1}F_k^{\mathsf T}+Q_k$ is positive semi-definite as a sum of two such, and the second being $>0$ because $R_k$ is positive definite. So $S_k$ is positive definite, hence invertible (Quadratic forms and the Rayleigh quotient). This is the one place the hypothesis $R_k \succ 0$ is used.

Step 4: the update. Conditioning in stages is legitimate: for any integrable $g$ , the tower property of Conditional expectation applied inside the conditional law gives $\mathbb{E}(g(x_k) \mid Y_{k-1}, y_k) = \mathbb{E}(g(x_k)\mid Y_k)$ , since $Y_k$ is exactly the pair $(Y_{k-1},y_k)$ . So the conditional law of $x_k$ given $Y_k$ may be computed by conditioning the jointly normal pair of Step 2 on its second block.

Apply the conditional formula of The multivariate normal distribution, statement (c), with $V_{11} = P_k^{f}$ , $V_{12} = P_k^{f}H_k^{\mathsf T}$ , $V_{22} = S_k$ . It gives a normal law with

\text{mean} = m_k^{f} + P_k^{f}H_k^{\mathsf T}S_k^{-1}\bigl(y_k - H_km_k^{f}\bigr) = m_k^{f} + K_kd_k ,

\text{covariance} = P_k^{f} - P_k^{f}H_k^{\mathsf T}S_k^{-1}H_kP_k^{f} = P_k^{f} - K_kH_kP_k^{f} = (I-K_kH_k)P_k^{f} .

The covariance is the Schur complement and does not contain $y_k$ ; since $P_{k-1}$ was non-random by $\mathcal{H}_{k-1}$ , so are $P_k^{f}$ , $S_k$ , $K_k$ and $P_k$ . This establishes $\mathcal{H}_k$ and completes the induction.

Step 5: independence of the innovations. By Step 2, $d_k = y_k - H_km_k^{f}$ has conditional law $N(0,S_k)$ given $Y_{k-1}$ , and $S_k$ is non-random. The conditional law therefore does not depend on $Y_{k-1}$ , so $d_k$ is independent of $Y_{k-1}$ with unconditional law $N(0,S_k)$ . Since $d_1,\dots,d_{k-1}$ are functions of $Y_{k-1}$ , the whole sequence is independent. The likelihood identity is then the factorisation $p(y_1,\dots,y_N) = \prod_k p(y_k\mid Y_{k-1})$ with each factor the $N(H_km_k^{f},S_k)$ density evaluated at $y_k$ , whose logarithm is $-\tfrac12(d_k^{\mathsf T}S_k^{-1}d_k + \log\det S_k + p\log2\pi)$ .

Step 6: the Joseph form. Expand, abbreviating $K = K_k$ , $H = H_k$ , $P^{f} = P_k^{f}$ , $R = R_k$ , $S = S_k$ :

(I-KH)P^{f}(I-KH)^{\mathsf T} + KRK^{\mathsf T} = P^{f} - KHP^{f} - P^{f}H^{\mathsf T}K^{\mathsf T} + K\bigl(HP^{f}H^{\mathsf T}+R\bigr)K^{\mathsf T} .

The last bracket is $S$ , and $KSK^{\mathsf T} = P^{f}H^{\mathsf T}S^{-1}SS^{-1}HP^{f} = P^{f}H^{\mathsf T}S^{-1}HP^{f} = KHP^{f}$ , a symmetric matrix; hence also $P^{f}H^{\mathsf T}K^{\mathsf T} = (KHP^{f})^{\mathsf T} = KHP^{f}$ . The right-hand side collapses to $P^{f} - KHP^{f} = (I-KH)P^{f} = P_k$ .

Step 7: the information form. Suppose $P^{f}$ is invertible. First,

KHP^{f}H^{\mathsf T} = P^{f}H^{\mathsf T}S^{-1}\bigl(S - R\bigr) = P^{f}H^{\mathsf T} - KR , \qquad\text{so}\qquad P^{f}H^{\mathsf T} - KHP^{f}H^{\mathsf T} = KR .

Now multiply out, using this identity at the marked step:

(I-KH)P^{f}\Bigl(\bigl(P^{f}\bigr)^{-1} + H^{\mathsf T}R^{-1}H\Bigr) = (I-KH) + \bigl(P^{f}H^{\mathsf T} - KHP^{f}H^{\mathsf T}\bigr)R^{-1}H \overset{(\ast)}{=} I - KH + KRR^{-1}H = I .

So $P_k = (I-KH)P^{f}$ is invertible with inverse $(P^{f})^{-1} + H^{\mathsf T}R^{-1}H$ . Finally

P_kH^{\mathsf T}R^{-1} = (I-KH)P^{f}H^{\mathsf T}R^{-1} = \bigl(P^{f}H^{\mathsf T} - KHP^{f}H^{\mathsf T}\bigr)R^{-1} = KRR^{-1} = K ,

by the same identity. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Linearity of the dynamics and of the observation operator is used in Step 2, and only there: it is what makes $(x_k,y_k)$ an affine image of a normal vector. Replace $H_k x_k$ by a non-linear $h(x_k)$ and the pair is no longer jointly normal; the filtering distribution is then generally not normal and is not determined by two moments. The standard repairs — linearising $h$ about the forecast (the extended filter) or propagating a sample (ensemble and sigma-point filters) — are approximations, not identities.
Gaussianity is used only to go from "the conditional mean and covariance are these" to "the conditional distribution is this". Drop it and keep the first and second moments, and the recursion still computes something meaningful: $m_k$ remains the best linear unbiased estimate of $x_k$ from $Y_k$ , by the Projection theorem (7.9.14) of Grimmett and Stirzaker quoted in Conditional expectation, because the orthogonality condition $\mathbb{E}\bigl((x_k-m_k)d_j^{\mathsf T}\bigr) = 0$ holds for $j \le k$ by construction. What is lost is that it is the best predictor of any kind: for a non-normal law the conditional mean is a non-linear function of the data and beats every linear rule.
Independence of the noises from each other and from the past is Step 1. If $\operatorname{cov}(w_k,v_k) = C_k \ne 0$ — common when the "observation" is a derived quantity contaminated by the same error as the forecast — then $\operatorname{cov}(x_k,y_k\mid Y_{k-1}) = P_k^{f}H_k^{\mathsf T} + C_k$ and the gain acquires that extra term. Using the stated formula anyway produces a $P_k$ that understates the true error.
$R_k$ positive definite is Step 3, and nothing else. A perfect observation, $R_k = 0$ , makes $S_k$ singular whenever $H_kP_k^{f}H_k^{\mathsf T}$ is, and the conditional law degenerates onto an affine subspace. $Q_k$ and $P_0$ , by contrast, may be singular, which is why the model is stated with them only semi-definite.
Known $F_k,H_k,Q_k,R_k$ . In practice $Q$ and $R$ are not known. Part (iii) is the escape: the log-likelihood of the record is available as a by-product of the filter, so the covariances can be estimated by maximising it, which is the method of Maximum likelihood applied to a dependent sample.
Numerical caution. The form $P_k = (I-K_kH_k)P_k^{f}$ is the cheapest but is not symmetric as written, and rounding can make it indefinite after many steps. The Joseph form costs more and keeps symmetry and positivity, because it is a sum of two expressions of the form $MVM^{\mathsf T}$ with $V$ positive semi-definite; square-root (Cholesky) implementations go further and propagate a factor of $P_k$ instead of $P_k$ itself (Linear systems and rank).

Remark(Three readings of the update step)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.6.1, eqs. (31.95)-(31.99), pp. 1273-1274; Riley, Hobson & Bence (2006), §31.5.5, eqs. (31.81)-(31.83), pp. 1264-1265

The same two formulas can be reached from three directions, and it is worth seeing that they agree.

As a conditional law. The route taken above.

As a Bayesian posterior. Treat the forecast $N(m_k^{f},P_k^{f})$ as the prior for $x_k$ and $y_k = H_kx_k+v_k$ as the likelihood. Posterior $\propto$ prior $\times$ likelihood (Bayesian inference) gives an exponent

-\tfrac12\Bigl[(x-m_k^{f})^{\mathsf T}(P_k^{f})^{-1}(x-m_k^{f}) + (y_k-H_kx)^{\mathsf T}R_k^{-1}(y_k-H_kx)\Bigr] ,

a quadratic in $x$ whose matrix is $(P_k^{f})^{-1}+H_k^{\mathsf T}R_k^{-1}H_k$ — the information form of part (iv). The scalar normal-normal conjugacy recorded in Bayesian inference, "precisions add and the mean is the precision-weighted average", is the case $n=p=1$ , $H=1$ .

As weighted least squares. Regard the forecast as an extra block of data. In the notation of Least-squares fitting, take the stacked observation vector, design (response) matrix and noise covariance

z = \begin{pmatrix} m_k^{f} \\ y_k\end{pmatrix}, \qquad R_{\text{des}} = \begin{pmatrix} I \\ H_k \end{pmatrix}, \qquad \mathsf{N} = \begin{pmatrix} P_k^{f} & 0 \\ 0 & R_k\end{pmatrix} ,

and estimate the parameter vector $x_k$ . Riley, Hobson and Bence's normal equations give the estimator $(R_{\text{des}}^{\mathsf T}\mathsf{N}^{-1}R_{\text{des}})^{-1}R_{\text{des}}^{\mathsf T}\mathsf{N}^{-1}z$ with covariance $(R_{\text{des}}^{\mathsf T}\mathsf{N}^{-1}R_{\text{des}})^{-1}$ ; here

R_{\text{des}}^{\mathsf T}\mathsf{N}^{-1}R_{\text{des}} = (P_k^{f})^{-1}+H_k^{\mathsf T}R_k^{-1}H_k , \qquad R_{\text{des}}^{\mathsf T}\mathsf{N}^{-1}z = (P_k^{f})^{-1}m_k^{f} + H_k^{\mathsf T}R_k^{-1}y_k ,

which is again the information form. So $m_k$ is a least-squares estimate and $P_k$ is the least-squares covariance matrix, and by the Gauss-Markov theorem of that unit it is the best linear unbiased combination of forecast and observation. This reading is the one the next unit, Optimal interpolation and variational data assimilation, develops.

Method(Running the filter)Original

Write the model. Choose the state $x$ so that the dynamics are linear and the observations are a linear function of it. Fix $F_k, H_k, Q_k, R_k$ and the prior $m_0, P_0$ . A deliberately inflated $P_0$ costs little: its influence decays.
Forecast. $m^{f} \leftarrow F_km_{k-1}$ and $P^{f} \leftarrow F_kP_{k-1}F_k^{\mathsf T}+Q_k$ .
Innovation. $d \leftarrow y_k - H_km^{f}$ and $S \leftarrow H_kP^{f}H_k^{\mathsf T}+R_k$ .
Gain. Solve $S\,K^{\mathsf T} = H_kP^{f}$ for $K^{\mathsf T}$ by Cholesky factorisation of $S$ rather than forming $S^{-1}$ (Linear systems and rank); $S$ is positive definite, so the factorisation exists. Cost $O(p^3 + np^2)$ per step, independent of $k$ .
Update. $m_k \leftarrow m^{f}+Kd$ and $P_k \leftarrow (I-KH_k)P^{f}$ , or the Joseph form if the record is long.
Monitor. Accumulate $d^{\mathsf T}S^{-1}d$ and $\log\det S$ for the likelihood of part (iii). The standardised innovations should look like independent standard normal variables; if $d^{\mathsf T}S^{-1}d$ is persistently larger than $p$ , the assumed $Q$ or $R$ is too small or the model is wrong, and the test of Hypothesis testing on $\sum_k d_k^{\mathsf T}S_k^{-1}d_k \sim \chi^2(Np)$ makes this quantitative.
If the model is time-invariant, iterate step 2-5 covariances alone until $P^{f}$ stops changing, and use the limiting gain throughout; see the example.

Example(A scalar red-noise record: four steps, by hand and by brute force)Original

The model. A quantity follows the AR(1) scheme of Autoregressive models and red noise and is measured with independent error:

x_k = \tfrac12 x_{k-1} + w_k , \quad w_k \sim N(0,\tfrac78) , \qquad y_k = x_k + v_k , \quad v_k \sim N(0,1) ,

so $n=p=1$ , $F=\tfrac12$ , $H=1$ , $Q=\tfrac78$ , $R=1$ . The prior is $x_0 \sim N(0,4)$ : deliberately vague, four times the observation variance. The four observations are

y_1 = 1.2, \quad y_2 = 0.4, \quad y_3 = -0.6, \quad y_4 = 0.1 .

Step 1: the fixed point first. For a time-invariant scalar model the forecast variance obeys $P^{f} \mapsto F^2\bigl(P^{f} - (P^{f})^2/(P^{f}+R)\bigr)+Q$ , whose fixed point solves

\bigl(P^{f}\bigr)^2 + P^{f}\bigl(R - F^2R - Q\bigr) - QR = 0 .

With $R=1$ , $F^2 = \tfrac14$ and $Q = \tfrac78$ the coefficient of $P^{f}$ is $1-\tfrac14-\tfrac78 = -\tfrac18$ and the constant term is $-\tfrac78$ , so

\bigl(P^{f}\bigr)^2 - \tfrac18P^{f} - \tfrac78 = 0 \quad\Longrightarrow\quad P^{f}_{\infty} = 1 ,

the other root $-\tfrac78$ being negative and so inadmissible for a variance.

The limiting gain and filtered variance are then $K_\infty = 1/(1+1) = \tfrac12$ and $P_\infty = (1-\tfrac12)\cdot 1 = \tfrac12$ . Consistency check of the fixed point: $F^2P_\infty + Q = \tfrac14\cdot\tfrac12+\tfrac78 = \tfrac18+\tfrac78 = 1 = P^{f}_\infty$ .

Step 2: the four steps, in exact arithmetic. With $m_0=0$ , $P_0=4$ :

$k$	$m_k^{f}$	$P_k^{f}$	$S_k$	$K_k$	$d_k$	$m_k$	$P_k$
1	$0$	$\tfrac{15}{8}=1.875$	$\tfrac{23}{8}$	$\tfrac{15}{23}=0.652174$	$1.2$	$\tfrac{18}{23}=0.782609$	$0.652174$
2	$0.391304$	$\tfrac{191}{184}=1.038043$	$\tfrac{375}{184}$	$\tfrac{191}{375}=0.509333$	$0.008696$	$0.395733$	$0.509333$
3	$0.197867$	$\tfrac{3007}{3000}=1.002333$	$\tfrac{6007}{3000}$	$\tfrac{3007}{6007}=0.500583$	$-0.797867$	$-0.201532$	$0.500583$
4	$-0.100766$	$\tfrac{48063}{48056}=1.000146$	$\tfrac{96119}{48056}$	$\tfrac{48063}{96119}=0.500036$	$0.200766$	$-0.000376$	$0.500036$

One line in full, to show the arithmetic. At $k=2$ : $m_2^{f} = \tfrac12\cdot\tfrac{18}{23} = \tfrac{9}{23} = 0.391304$ and $P_2^{f} = \tfrac14\cdot\tfrac{15}{23}+\tfrac78 = \tfrac{15}{92}+\tfrac78 = \tfrac{30+161}{184} = \tfrac{191}{184}$ . Then $S_2 = \tfrac{191}{184}+1 = \tfrac{375}{184}$ , $K_2 = \tfrac{191}{375}$ , $d_2 = 0.4 - \tfrac{9}{23} = \tfrac{1}{115} = 0.008696$ , and $m_2 = \tfrac{9}{23}+\tfrac{191}{375}\cdot\tfrac{1}{115} = \tfrac{742}{1875} = 0.395733$ , $P_2 = \bigl(1-\tfrac{191}{375}\bigr)\tfrac{191}{184} = \tfrac{191}{375}$ . (With $R=1$ one always has $P_k = K_k$ here, since $P_k = P_k^{f}R/(P_k^{f}+R)$ .)

Step 3: how fast the variance converges. The departures from the fixed point are

P_k^{f}-1 : \quad 0.875,\ \ 0.0380435,\ \ 0.00233333,\ \ 0.000145663 ,

with successive ratios $0.04348$ , $0.06133$ , $0.06243$ , approaching $F^2(1-K_\infty)^2 = \tfrac14\cdot\tfrac14 = \tfrac{1}{16} = 0.0625$ . Three observations have already wiped out the vague prior: by $k=3$ the variance is within $0.2\%$ of its limit.

Step 4: the answer, and an independent check. The filtered estimate after four observations is

\boxed{\,x_4 \mid Y_4 \;\sim\; N\bigl(-0.000376,\; 0.500036\bigr)\,}, \qquad \text{i.e. } x_4 = -0.0004 \pm 0.7071 .

As a check that owes nothing to the recursion, note that $(x_0,x_1,\dots,x_4,y_1,\dots,y_4)$ is one nine-dimensional normal vector: it is the image of the nine independent variables $(x_0,w_1,\dots,w_4,v_1,\dots,v_4)$ under an explicit $9\times 9$ matrix $M$ , so its covariance is $M\operatorname{diag}(4,\tfrac78,\tfrac78,\tfrac78,\tfrac78,1,1,1,1)M^{\mathsf T}$ . Conditioning the $x_4$ component on the four $y$ components by the single Schur-complement formula of The multivariate normal distribution — one $4\times4$ inversion, no recursion — was carried out in exact rational arithmetic and gave

\text{mean} = -\tfrac{361}{961190} = -0.00037558 , \qquad \text{variance} = \tfrac{48063}{96119} = 0.50003641 ,

the same rationals the table produces. The point of the filter is that the direct route needs a $k \times k$ inversion that grows with the record, while the recursion needs one $1\times1$ division per step.

Reading the numbers. The data are noisy ( $R=1$ ) and the signal is modestly persistent ( $F=\tfrac12$ ), so the filter settles on trusting each new observation exactly half as much as it trusts its own forecast, and halves the variance from $1$ to $\tfrac12$ at each update. The final estimate is near zero although the last observation was $+0.1$ : the previous value was estimated at $-0.20$ , the forecast from it was $-0.10$ , and the innovation $+0.20$ was only half believed.

Application(Where the recursion earns its keep)Original

Navigation and tracking. The original use. The state is position, velocity and attitude; $F$ encodes the equations of motion over one time step, $Q$ the unmodelled accelerations, and the observations are ranges, Doppler shifts or inertial readings. Every satellite navigation receiver runs a filter of this kind, and the reported accuracy is the diagonal of $P_k$ .

Estimating a trend as it happens. A trend plus red noise is a two-component state, $x_k = (\text{level},\text{slope})$ , with $F = \begin{pmatrix}1 & 1\\ 0 & 1\end{pmatrix}$ . The filter then produces an estimate of the current level and slope from all data so far, with uncertainties, and does not need to be refitted when a new point arrives. This is the recursive counterpart of Least-squares fitting, and the two agree when the model is static ( $F=I$ , $Q=0$ ).

Calibration and fault detection. Part (iii) turns the filter into a test. The standardised innovations are independent standard normals when the model is right, so a run of same-signed innovations, or a $d_k^{\mathsf T}S_k^{-1}d_k$ far above $p$ , is evidence of a drifting instrument or a wrong $F$ . Monitoring the innovations is cheaper and more sensitive than monitoring the state estimate itself.

Geophysical and climate analysis. When the state is a discretised field with $n$ in the millions, $P_k$ cannot be stored, let alone propagated, and the exact recursion is replaced by approximations that keep its structure: the gain is built from a prescribed or ensemble-estimated covariance. That is the subject of Optimal interpolation and variational data assimilation, where the forecast covariance $P_k^{f}$ becomes the prescribed background covariance $B$ and the update step becomes optimal interpolation.

Method · advanced

Optimal interpolation and variational data assimilation

Adapted

Add to PDF

Derives the best linear unbiased analysis $x_a = x_b + K(y-Hx_b)$ with $K = BH^{\mathsf T}(HBH^{\mathsf T}+R)^{-1}$ by completing the square in the gain, proves that it also minimises the 3D-Var cost function and that the two routes agree exactly when $B$ is invertible, and interpolates a temperature analysis from two stations, obtaining weights $0.740445$ and $-0.135366$ , an analysis of $1.9490\,\mathrm{K}$ and an analysis error variance $0.425475\,\mathrm{K^2}$ against a background variance of $1$ .

Prerequisites: The Kalman filter, Least-squares fitting, The multivariate normal distribution, Stationary processes and spectral densities, Bayesian inference, Quadratic forms and the Rayleigh quotient, Linear systems and rank, Propagation of uncertainty

Note(The problem this solves)Original

A weather forecast starts from a state of the atmosphere: wind, temperature and humidity on a grid with far more unknowns than there are measurements. Observations are sparse, irregularly placed, of mixed type, and noisy. The only reason the problem is solvable at all is that a previous forecast, valid now, supplies a first guess everywhere. Data assimilation is the business of correcting that guess with the observations, and reporting the uncertainty of the result.

Two formulations are in use, and they look quite different. Optimal interpolation writes the correction as a matrix times the observed discrepancies and chooses the matrix to minimise the error variance. Three-dimensional variational assimilation (3D-Var) writes down a cost function measuring disagreement with the guess and with the observations, and minimises it. This unit derives the first, derives the second, and shows that under the stated hypotheses they produce the same analysis and the same error covariance. The difference is computational, not statistical: one inverts a matrix of the size of the observation vector, the other descends a function of the size of the state vector, and in operational practice the sizes differ by several orders of magnitude.

The update step of The Kalman filter is the same calculation. There the prior covariance was the forecast covariance $P_k^{f}$ produced by the filter's own recursion; here it is a prescribed background covariance $B$ , because a state of dimension $10^{8}$ - $10^{9}$ , depending on the model resolution, makes the recursion impossible to run.

Definition(Background, observations and the analysis)Original

Let the unknown true state be a fixed (non-random) vector $x \in \mathbb{R}^{n}$ . Two sources of information about it are given.

A background (or first guess) $x_b \in \mathbb{R}^{n}$ , with background error $e_b = x_b - x$ satisfying $\mathbb{E}[e_b] = 0 , \qquad \operatorname{Cov}[e_b] = \mathbb{E}\bigl[e_be_b^{\mathsf T}\bigr] = B ,$ with $B$ symmetric positive semi-definite and known.
An observation vector $y \in \mathbb{R}^{p}$ , related to the state by a known $p\times n$ observation operator $H$ (linear here), $y = Hx + \varepsilon , \qquad \mathbb{E}[\varepsilon] = 0 , \qquad \operatorname{Cov}[\varepsilon] = R ,$ with $R$ symmetric positive definite and known.

Assume further that the two error vectors are uncorrelated, $\mathbb{E}\bigl[e_b\varepsilon^{\mathsf T}\bigr] = 0$ . No distributional assumption is made: only first and second moments enter until the final remark.

The quantities built from these are

d \;=\; y - Hx_b \quad (\text{the \textbf{innovation}, or observation-minus-background}), \qquad S \;=\; HBH^{\mathsf T} + R \quad (\text{its covariance}) ,

the second because $d = y - Hx_b = \varepsilon - He_b$ , whence $\operatorname{Cov}[d] = HBH^{\mathsf T}+R$ by the covariance law of Propagation of uncertainty and the assumed lack of correlation.

An analysis is any estimate $x_a$ of $x$ built from $x_b$ and $y$ ; it is linear if $x_a = Lx_b + Ky$ for fixed matrices $L, K$ , and unbiased if $\mathbb{E}[x_a - x] = 0$ whatever the true $x$ . Its quality is measured by the analysis error covariance $A = \mathbb{E}\bigl[(x_a-x)(x_a-x)^{\mathsf T}\bigr]$ .

Method(The best linear unbiased analysis)Original

Under the hypotheses above:

(i) Unbiasedness fixes $L$ . A linear analysis $x_a = Lx_b + Ky$ is unbiased for every true $x$ exactly when $L = I - KH$ , equivalently exactly when it takes the form

x_a \;=\; x_b + K\,(y - Hx_b) \;=\; x_b + Kd .

The correction is thus always a matrix times the innovation; only $K$ , the gain, is free.

(ii) The error covariance as a function of the gain. For such an analysis,

A(K) \;=\; (I-KH)\,B\,(I-KH)^{\mathsf T} + K\,R\,K^{\mathsf T} .

(iii) The optimal gain. With $S = HBH^{\mathsf T}+R$ (invertible, since $R$ is positive definite), put

K_\ast \;=\; B\,H^{\mathsf T}S^{-1} .

Then for every $K$ ,

A(K) \;=\; A(K_\ast) + (K-K_\ast)\,S\,(K-K_\ast)^{\mathsf T} , \qquad A(K_\ast) = (I - K_\ast H)B = B - BH^{\mathsf T}S^{-1}HB ,

so $A(K) - A(K_\ast)$ is positive semi-definite for every $K$ , and is zero only for $K = K_\ast$ . The analysis $x_a = x_b + K_\ast d$ is therefore the best linear unbiased estimate (BLUE) of $x$ : it has the smallest error variance for every component, and indeed for every linear functional $u^{\mathsf T}x$ of the state, among all linear unbiased analyses. Computing $x_a$ and $A$ from these formulas is optimal interpolation.

(iv) Orthogonality. The optimal analysis error is uncorrelated with the innovation: $\mathbb{E}\bigl[(x_a-x)\,d^{\mathsf T}\bigr] = 0$ . Conversely this property determines $K_\ast$ .

(v) The information form. If in addition $B$ is invertible, then $A = A(K_\ast)$ is invertible and

A^{-1} \;=\; B^{-1} + H^{\mathsf T}R^{-1}H , \qquad K_\ast \;=\; A\,H^{\mathsf T}R^{-1} .

Method(3D-Var, and its equivalence with optimal interpolation)Original

Assume in addition that $B$ is invertible, and define the 3D-Var cost function $J : \mathbb{R}^{n}\to\mathbb{R}$ ,

J(x) \;=\; \tfrac12\,(x-x_b)^{\mathsf T}B^{-1}(x-x_b) \;+\; \tfrac12\,(y-Hx)^{\mathsf T}R^{-1}(y-Hx) .

Then:

$J$ is strictly convex with Hessian $\nabla^2 J = B^{-1}+H^{\mathsf T}R^{-1}H$ , constant in $x$ and positive definite, so $J$ has exactly one minimiser.
That minimiser is the best linear unbiased analysis of the previous statement: $\operatorname*{arg\,min}_{x} J(x) \;=\; x_b + K_\ast(y - Hx_b) \;=\; x_a .$
The inverse Hessian is the analysis error covariance: $\bigl(\nabla^2J\bigr)^{-1} = A$ .
The minimum value is $J(x_a) = \tfrac12\,d^{\mathsf T}S^{-1}d$ , a quadratic form in the innovation alone.

So optimal interpolation and 3D-Var are two ways of computing one object. If the errors are in addition jointly normal, both coincide with the Bayesian posterior: $x \mid y \sim N(x_a, A)$ , since $-\log$ of (prior $\times$ likelihood) is $J$ up to a constant, by Bayesian inference and The multivariate normal distribution.

Idea of the proofOriginal

Three short arguments, one for each part.

Unbiasedness. The analysis error of $Lx_b+Ky$ contains a term $(L+KH-I)x$ in the unknown true state. Since the estimate must be unbiased for every $x$ , and $x$ is not random, that term has to vanish identically, which pins $L$ down with no optimisation at all.

Optimality. $A(K)$ is a quadratic function of the matrix $K$ with positive-definite "second derivative" $S$ . Completing the square in $K$ — exactly as one minimises $aK^2+bK+c$ — exhibits $A(K)$ as a constant matrix plus $(K-K_\ast)S(K-K_\ast)^{\mathsf T}$ , and the second term is positive semi-definite and vanishes only at $K_\ast$ . Because the comparison is between whole matrices, this gives optimality for every linear functional of the state at once, not merely for the trace.

Equivalence. Setting the gradient of $J$ to zero gives $(B^{-1}+H^{\mathsf T}R^{-1}H)x = B^{-1}x_b + H^{\mathsf T}R^{-1}y$ . One algebraic identity, verified by multiplying out, says that the matrix on the left is $A^{-1}$ and that $AH^{\mathsf T}R^{-1} = K_\ast$ ; with it the solution rearranges into $x_b + K_\ast d$ in two lines. The same identity shows why $2J$ is Riley, Hobson and Bence's weighted $\chi^2$ for the design matrix that stacks the background on top of the observations: the factor $\tfrac12$ is a convention of the variational literature, and it is the same factor that turns their $\tfrac12\,\partial^2\chi^2/\partial a_i\partial a_j$ into our $\nabla^2 J$ .

Proof(Unbiasedness, optimality, the identity, and 3D-Var)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.8, Exercise 31.6, p. 1299; Riley, Hobson & Bence (2006), §31.6.1, eqs. (31.92)-(31.99), pp. 1272-1274

(i) Unbiasedness. Substituting $x_b = x + e_b$ and $y = Hx+\varepsilon$ ,

x_a - x = Lx_b + Ky - x = (L + KH - I)\,x + L\,e_b + K\,\varepsilon .

Taking expectations and using $\mathbb{E}[e_b] = \mathbb{E}[\varepsilon] = 0$ leaves $\mathbb{E}[x_a-x] = (L+KH-I)x$ . This vanishes for every $x \in \mathbb{R}^{n}$ exactly when the matrix $L+KH-I$ is zero, since a matrix annihilating all of $\mathbb{R}^{n}$ is the zero matrix. Hence $L = I-KH$ and $x_a = (I-KH)x_b + Ky = x_b + K(y-Hx_b)$ . This is the matrix version of the step in Riley, Hobson and Bence's Exercise 31.6 where unbiasedness of a linear estimator is imposed before efficiency is optimised. Note that the hypothesis that $x$ is a fixed unknown vector, not a random one, is what forces the condition to hold for all $x$ .

(ii) The error covariance. With $L = I-KH$ the display above gives

x_a - x = (I-KH)\,e_b + K\,\varepsilon .

Since $\mathbb{E}\bigl[e_b\varepsilon^{\mathsf T}\bigr] = 0$ , the cross terms in $\mathbb{E}\bigl[(x_a-x)(x_a-x)^{\mathsf T}\bigr]$ vanish and

A(K) = (I-KH)\,\mathbb{E}\bigl[e_be_b^{\mathsf T}\bigr]\,(I-KH)^{\mathsf T} + K\,\mathbb{E}\bigl[\varepsilon\varepsilon^{\mathsf T}\bigr]\,K^{\mathsf T} = (I-KH)B(I-KH)^{\mathsf T} + KRK^{\mathsf T} .

It is the lack of correlation between background and observation errors, and nothing else, that is used here.

(iii) Completing the square. $S = HBH^{\mathsf T}+R$ is positive definite: for $u\ne0$ , $u^{\mathsf T}Su = (H^{\mathsf T}u)^{\mathsf T}B(H^{\mathsf T}u) + u^{\mathsf T}Ru > 0$ , the first term being $\ge 0$ because $B$ is positive semi-definite and the second $>0$ because $R$ is positive definite (Quadratic forms and the Rayleigh quotient). So $S^{-1}$ exists and is positive definite, and $K_\ast = BH^{\mathsf T}S^{-1}$ is defined.

Expand (ii):

A(K) = B - KHB - BH^{\mathsf T}K^{\mathsf T} + K\bigl(HBH^{\mathsf T}+R\bigr)K^{\mathsf T} = B - KHB - BH^{\mathsf T}K^{\mathsf T} + KSK^{\mathsf T} .

Separately, since $SK_\ast^{\mathsf T} = SS^{-1}HB = HB$ and $K_\ast S K_\ast^{\mathsf T} = BH^{\mathsf T}S^{-1}HB$ ,

(K-K_\ast)S(K-K_\ast)^{\mathsf T} = KSK^{\mathsf T} - KHB - BH^{\mathsf T}K^{\mathsf T} + BH^{\mathsf T}S^{-1}HB .

Subtracting, every $K$ -dependent term cancels:

A(K) - (K-K_\ast)S(K-K_\ast)^{\mathsf T} = B - BH^{\mathsf T}S^{-1}HB ,

a matrix independent of $K$ . Putting $K = K_\ast$ identifies it as $A(K_\ast) = B - K_\ast HB = (I-K_\ast H)B$ , and the stated decomposition follows.

For any $u \in \mathbb{R}^{n}$ , $u^{\mathsf T}\bigl(A(K)-A(K_\ast)\bigr)u = v^{\mathsf T}Sv \ge 0$ with $v = (K-K_\ast)^{\mathsf T}u$ , so the difference is positive semi-definite; taking $u$ to be a coordinate vector gives the componentwise variance statement. If the difference is the zero matrix then $v^{\mathsf T}Sv = 0$ for every $u$ , and positive definiteness of $S$ forces $v = 0$ for every $u$ , that is $K = K_\ast$ .

(iv) Orthogonality. From (ii), $x_a - x = (I-KH)e_b + K\varepsilon$ and $d = \varepsilon - He_b$ , so

\mathbb{E}\bigl[(x_a-x)d^{\mathsf T}\bigr] = -(I-KH)BH^{\mathsf T} + KR = -BH^{\mathsf T} + K\bigl(HBH^{\mathsf T}+R\bigr) = KS - BH^{\mathsf T} ,

using $\mathbb{E}[e_b\varepsilon^{\mathsf T}]=0$ twice. This vanishes precisely when $K = BH^{\mathsf T}S^{-1} = K_\ast$ . Geometrically this is the Projection theorem (7.9.14) of Grimmett and Stirzaker, applied not to $x$ itself — which is a fixed vector, not a member of any $L^2$ space — but to the random part of the problem: the correction $Kd$ is the orthogonal projection of the background error $-e_b$ onto the span of the components of $d$ , since $x_a - x = (x_b - x) + Kd = e_b + Kd$ is exactly minus the residual $(-e_b) - Kd$ of that projection. Vanishing of the display above is, up to sign, the orthogonality condition $\mathbb{E}\bigl[\bigl((-e_b) - Kd\bigr)d^{\mathsf T}\bigr] = 0$ of that theorem, and it says that the analysis error is uncorrelated with every linear function of the data: $d$ carries all the information the data hold about $e_b$ .

(v) The information identity. Assume $B$ invertible. First note

K_\ast HBH^{\mathsf T} = BH^{\mathsf T}S^{-1}\bigl(S-R\bigr) = BH^{\mathsf T} - K_\ast R , \qquad\text{so}\qquad BH^{\mathsf T} - K_\ast HBH^{\mathsf T} = K_\ast R . \tag{$\ast$}

Then

A\bigl(B^{-1}+H^{\mathsf T}R^{-1}H\bigr) = (I-K_\ast H)B\bigl(B^{-1}+H^{\mathsf T}R^{-1}H\bigr) = (I-K_\ast H) + \bigl(BH^{\mathsf T}-K_\ast HBH^{\mathsf T}\bigr)R^{-1}H ,

and by $(\ast)$ the last bracket is $K_\ast R$ , so the whole expression is $I - K_\ast H + K_\ast RR^{-1}H = I$ . Hence $A$ is invertible with $A^{-1} = B^{-1}+H^{\mathsf T}R^{-1}H$ . The same identity gives

AH^{\mathsf T}R^{-1} = (I-K_\ast H)BH^{\mathsf T}R^{-1} = \bigl(BH^{\mathsf T}-K_\ast HBH^{\mathsf T}\bigr)R^{-1} = K_\ast RR^{-1} = K_\ast .

(vi) 3D-Var. Differentiating the two quadratic forms,

\nabla J(x) = B^{-1}(x-x_b) - H^{\mathsf T}R^{-1}(y-Hx) , \qquad \nabla^2 J = B^{-1}+H^{\mathsf T}R^{-1}H ,

the Hessian being constant because $J$ is quadratic. It is positive definite ( $B^{-1}$ is, and $H^{\mathsf T}R^{-1}H$ is positive semi-definite), so $J$ is strictly convex and its unique stationary point is its unique minimum. Setting $\nabla J = 0$ ,

\bigl(B^{-1}+H^{\mathsf T}R^{-1}H\bigr)x = B^{-1}x_b + H^{\mathsf T}R^{-1}y .

By (v) the matrix on the left is $A^{-1}$ . Adding and subtracting $H^{\mathsf T}R^{-1}Hx_b$ on the right,

x = A\Bigl(\bigl(B^{-1}+H^{\mathsf T}R^{-1}H\bigr)x_b + H^{\mathsf T}R^{-1}(y-Hx_b)\Bigr) = x_b + AH^{\mathsf T}R^{-1}d = x_b + K_\ast d = x_a ,

using (v) once more for $AH^{\mathsf T}R^{-1}=K_\ast$ . Claim 3 is the first sentence of (v). For claim 4, write $x_a - x_b = K_\ast d$ and $y - Hx_a = d - HK_\ast d = (I - HK_\ast)d$ . Since $HK_\ast = HBH^{\mathsf T}S^{-1} = (S-R)S^{-1} = I - RS^{-1}$ , this is $RS^{-1}d$ , so

2J(x_a) = d^{\mathsf T}K_\ast^{\mathsf T}B^{-1}K_\ast d + d^{\mathsf T}S^{-1}R\,R^{-1}R\,S^{-1}d = d^{\mathsf T}S^{-1}\bigl(HBH^{\mathsf T}\bigr)S^{-1}d + d^{\mathsf T}S^{-1}RS^{-1}d = d^{\mathsf T}S^{-1}d ,

where $K_\ast^{\mathsf T}B^{-1}K_\ast = S^{-1}HB\,B^{-1}BH^{\mathsf T}S^{-1} = S^{-1}HBH^{\mathsf T}S^{-1}$ and the two terms recombine through $HBH^{\mathsf T}+R = S$ . $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Uncorrelated background and observation errors is used once, in part (ii), and it is the hypothesis most often violated in practice. If $\mathbb{E}[e_b\varepsilon^{\mathsf T}] = C \ne 0$ — for instance when an observation has already been used to make the background, as happens if the same radiosonde enters two successive analyses — then $A(K)$ acquires the cross terms $(I-KH)CK^{\mathsf T} + KC^{\mathsf T}(I-KH)^{\mathsf T}$ , and repeating the derivative of part (ii) gives the optimal gain

K_\ast^{C} = \bigl(BH^{\mathsf T} - C\bigr)\bigl(HBH^{\mathsf T} + R - HC - C^{\mathsf T}H^{\mathsf T}\bigr)^{-1} .

The minus signs are forced by our conventions, and are worth tracing: with $e_b = x_b - x$ and $y = Hx + \varepsilon$ the innovation is $d = \varepsilon - He_b$ , so

\operatorname{Cov}[d] = HBH^{\mathsf T} + R - HC - C^{\mathsf T}H^{\mathsf T} , \qquad \mathbb{E}\bigl[e_b d^{\mathsf T}\bigr] = C - BH^{\mathsf T} ,

and $K_\ast^{C} = -\mathbb{E}[e_bd^{\mathsf T}]\operatorname{Cov}[d]^{-1}$ exactly as in the uncorrelated case; writing a plus in the denominator would not even give a covariance. A scalar check, $H = 1$ , $B = 2$ , $R = 0.5$ , $C = 0.6$ : the analysis variance is $A(K) = 2(1-K)^2 + \tfrac12 K^2 + 1.2\,K(1-K)$ , minimised at $K = \tfrac{14}{13} = 1.076923$ with $A = \tfrac{32}{65} = 0.492308$ , and the formula above returns $1.4/1.3 = \tfrac{14}{13}$ exactly. The plus-sign version would give $K = 2.6/3.7 = 0.702703$ , where $A = 0.674361$ — worse by $37\%$ .

Using the uncorrelated $K_\ast$ anyway leaves the analysis unbiased, but its true error no longer matches the reported $A$ . The cross terms above give the discrepancy exactly:

A_{\text{true}} - A \;=\; (I-K_\ast H)\,C\,K_\ast^{\mathsf T} + K_\ast C^{\mathsf T}(I-K_\ast H)^{\mathsf T} .

A matrix of this form is in general indefinite, whatever the sign of $C$ , so the reported covariance can err in either direction, and in different directions in different components. Take $n = p = 2$ , $H = I$ and

B = \begin{pmatrix}5&4\\4&5\end{pmatrix} , \qquad R = \begin{pmatrix}4&-3\\-3&5\end{pmatrix} , \qquad C = \begin{pmatrix}2&\tfrac12\\\tfrac12&\tfrac14\end{pmatrix} ,

all three positive definite, and the joint error covariance $\left(\begin{smallmatrix}B&C\\C^{\mathsf T}&R\end{smallmatrix}\right)$ positive definite as well (its eigenvalues are $0.108328$ , $1.938988$ , $7.323769$ , $9.628915$ ). Then $S = \left(\begin{smallmatrix}9&1\\1&10\end{smallmatrix}\right)$ , the reported variances are $\operatorname{diag}A = (91/89,\,100/89) = (1.022472,\,1.123596)$ and the true ones are $(1.948933,\,0.660144)$ : the first component is understated, the second overstated, although $C$ is positive definite. Replacing $C$ by $-C$ keeps the joint covariance positive definite — it is congruent to the original by $\operatorname{diag}(I,-I)$ — and simply exchanges the two directions of error.

Only in the scalar case is there a clean rule, and only for $H > 0$ . There $0 < K_\ast H < 1$ , the discrepancy is $2(1-K_\ast H)K_\ast C$ , and it carries the sign of $C$ : positive correlation — the radiosonde case, where the shared observation has pushed background and observation errors the same way — makes the true error exceed the reported $A$ , so that the system is overconfident, which is worse than being wrong, while negative correlation leaves the reported $A$ merely pessimistic. The check above, $H=1$ , $B=2$ , $R=\tfrac12$ , $C=0.6$ , is of this kind: the reported $A = (1-K_\ast)B = 0.4$ against a true $0.592$ . Even here the sign of $H$ matters: $H = -1$ , $B = 2$ , $R = \tfrac12$ , $C = 0.4$ gives $K_\ast = -0.8$ , a reported $A = 0.4$ and a true $0.272$ . Overconfidence is the danger to watch for, but which way the error goes has to be computed, not inferred from the sign of $C$ .

Unbiased errors. Part (i) needs $\mathbb{E}[e_b] = 0$ . A systematic model bias $\beta = \mathbb{E}[e_b] \ne 0$ propagates straight through: the analysis inherits the bias $(I-KH)\beta$ , and no choice of gain removes it, because the derivation of (i) then yields $\mathbb{E}[x_a-x] = (I-KH)\beta$ for every admissible $K$ . Bias must be estimated and removed separately; this is why operational systems carry a bias-correction stage ahead of the analysis.
$R$ positive definite is used to invert $S$ in parts (iii)-(iv), and again as $R^{-1}$ itself from part (v) onwards, in the cost function and throughout (vi). $B$ may be singular — a background known exactly along some direction is allowed — and parts (i)-(iv) survive. Part (v) and the whole of 3D-Var need $B^{-1}$ , so the variational form is strictly less general than optimal interpolation. This matters: ensemble-based background covariances have rank below $n$ and are singular, which is why they are used through the gain formula, or regularised before $B^{-1}$ is formed.
Linearity of $H$ . With a non-linear observation operator $h$ the cost function $J(x) = \tfrac12(x-x_b)^{\mathsf T}B^{-1}(x-x_b)+\tfrac12(y-h(x))^{\mathsf T}R^{-1}(y-h(x))$ is still minimised numerically, but it need not be convex and the minimiser is no longer a linear function of $y$ . Linearising $h$ about the background and iterating (incremental 3D-Var) recovers a sequence of problems of the kind treated here; it converges to a stationary point of $J$ , not necessarily the global one.
No normality is needed for the BLUE statement, exactly as in the Gauss-Markov theorem of Least-squares fitting: only the first two moments enter. Normality buys one extra thing, stated in the last line of the 3D-Var statement, namely that $x_a$ is the conditional mean and therefore the best estimate of any form, linear or not.
Which costs less. Optimal interpolation inverts $S$ , of size $p\times p$ ; 3D-Var minimises over $\mathbb{R}^{n}$ , usually by a conjugate-gradient method whose cost per iteration is one multiplication by $B$ and one by $H$ . For a modern forecasting system $n \sim 10^{8}$ - $10^{9}$ and $p \sim 10^{7}$ , so neither matrix can be formed and both routes are run in a reduced space; the variational form is preferred because $B$ can be applied as a sequence of filters without ever being stored.

Method(Optimal interpolation at a point)Original

The classical use is to analyse one grid point at a time from nearby observations, which keeps every matrix small. Let the analysis point be $P$ and let $p$ observations be selected within some radius of it.

Specify the background error statistics. Write $B_{ij} = \sigma_b(i)\,\sigma_b(j)\,\rho(r_{ij})$ , with $\sigma_b$ the background error standard deviation and $\rho$ a correlation function of the separation $r_{ij}$ , decreasing from $\rho(0)=1$ . The choice of $\rho$ must be a valid correlation function, i.e. non-negative definite (Stationary processes and spectral densities); the Gaussian $\rho(r) = e^{-r^2/(2L^2)}$ is the usual default, $L$ being the correlation length.
Form $S = HBH^{\mathsf T}+R$ , the $p\times p$ matrix of background error covariances between the observation locations, plus the observation error covariance $R$ (diagonal if instrument errors are independent).
Form $b =$ the column of background error covariances between $P$ and the observation locations, that is the row of $BH^{\mathsf T}$ belonging to $P$ .
Solve $S\,w = b$ for the weight vector $w$ , by Cholesky factorisation (Linear systems and rank); do not invert $S$ .
Analyse. $x_a(P) = x_b(P) + w^{\mathsf T}d$ with $d$ the vector of innovations at the observation locations.
Report the uncertainty. The analysis error variance at $P$ is $\sigma_b(P)^2 - w^{\mathsf T}b$ , the diagonal entry of $A = B - BH^{\mathsf T}S^{-1}HB$ . It is never larger than $\sigma_b(P)^2$ , and the reduction $w^{\mathsf T}b$ measures what the observations were worth.
Check. Verify $2J(x_a) = d^{\mathsf T}S^{-1}d$ against its expectation $p$ : the innovations have covariance $S$ , so $\mathbb{E}\bigl[d^{\mathsf T}S^{-1}d\bigr] = \operatorname{tr}\bigl(S^{-1}S\bigr) = p$ . A value far above $p$ means $B$ or $R$ is too small, far below means they are too large.

Example(A temperature analysis from two stations)Original

Set-up. Three points lie on a line: the analysis point $P$ at $0$ , station $1$ at $60\,\mathrm{km}$ and station $2$ at $150\,\mathrm{km}$ . We analyse a temperature anomaly. The background field is uniform, $x_b = 1.00\,\mathrm{K}$ at all three points, with background error standard deviation $\sigma_b = 1\,\mathrm{K}$ everywhere and correlation

\rho(r) = \exp\bigl(-r^2/(2L^2)\bigr) , \qquad L = 100\,\mathrm{km} .

Both stations measure with independent error of standard deviation $0.5\,\mathrm{K}$ , so $R = 0.25\,I_2$ in $\mathrm{K^2}$ and $H$ simply picks out the two station values. The readings are $y_1 = 2.30\,\mathrm{K}$ and $y_2 = 1.10\,\mathrm{K}$ .

Step 1: the covariances. The three separations needed are $60$ , $150$ and $90$ km:

\rho_{P1} = e^{-0.18} = 0.835270 , \qquad \rho_{P2} = e^{-1.125} = 0.324652 , \qquad \rho_{12} = e^{-0.405} = 0.666977 .

Hence, in $\mathrm{K^2}$ ,

b = \begin{pmatrix}0.835270\\ 0.324652\end{pmatrix} , \qquad S = HBH^{\mathsf T}+R = \begin{pmatrix}1 & 0.666977\\ 0.666977 & 1\end{pmatrix} + \begin{pmatrix}0.25&0\\0&0.25\end{pmatrix} = \begin{pmatrix}1.25 & 0.666977\\ 0.666977 & 1.25\end{pmatrix} .

Step 2: the weights. $\det S = 1.25^2 - 0.666977^2 = 1.5625 - 0.444858 = 1.117642$ , so solving $Sw = b$ ,

w_1 = \frac{1.25\,(0.835270) - 0.666977\,(0.324652)}{1.117642} = \frac{0.827552}{1.117642} = 0.740445 ,

w_2 = \frac{1.25\,(0.324652) - 0.666977\,(0.835270)}{1.117642} = \frac{-0.151291}{1.117642} = -0.135366 .

The second weight is negative. This is not an error and it is the most instructive number here. Station $1$ alone would receive weight $\rho_{P1}/1.25 = 0.668216$ ; it is also strongly correlated with station $2$ , and through station $1$ the scheme already predicts an innovation of $\rho_{12} \times$ (station 1's innovation) at station $2$ . Since $\rho_{P1}\rho_{12} = 0.557106$ exceeds $\rho_{P2} = 0.324652$ , that prediction overshoots what station $2$ is worth at $P$ , and the optimal scheme subtracts the excess. Observations are never weighted in isolation: $S^{-1}$ decorrelates them first. This is the same screening effect that makes the far station's residual informative mainly as a correction to the near one.

Step 3: the analysis. The innovations are $d_1 = 2.30-1.00 = 1.30\,\mathrm{K}$ and $d_2 = 1.10-1.00 = 0.10\,\mathrm{K}$ , so

x_a(P) = 1.00 + (0.740445)(1.30) + (-0.135366)(0.10) = 1.00 + 0.9625785 - 0.0135366 = 1.949042\,\mathrm{K} .

Step 4: the error bar. By step 6 of the method,

A_{PP} = \sigma_b^2 - w^{\mathsf T}b = 1 - \bigl[(0.740445)(0.835270) + (-0.135366)(0.324652)\bigr] = 1 - 0.574525 = 0.425475\,\mathrm{K^2} ,

that is $\pm 0.652285\,\mathrm{K}$ , against the background's $\pm 1\,\mathrm{K}$ . Station $1$ alone would have given $1 - (0.835270)^2/1.25 = 0.441859\,\mathrm{K^2}$ : the second station, negative weight and all, is worth a further reduction of only $0.016384\,\mathrm{K^2}$ , because almost everything it knows about $P$ is already known through station $1$ .

Step 5: the same answer from 3D-Var. Take the state to be the three-vector $(x_P, x_1, x_2)$ with the full $3\times3$ background covariance built from the same $\rho$ (its eigenvalues are $0.074127$ , $0.685117$ , $2.240756$ , all positive, so $B^{-1}$ exists) and $H = \begin{pmatrix}0&1&0\\0&0&1\end{pmatrix}$ . Minimising

J(x) = \tfrac12 (x-x_b)^{\mathsf T}B^{-1}(x-x_b) + \tfrac12 (y-Hx)^{\mathsf T}R^{-1}(y-Hx)

numerically from the background as starting point returned

x_a = (1.949042,\; 1.951431,\; 1.265990) , \qquad J(x_a) = 0.873082 ,

whose first component is the analysis of Step 3, to all digits shown; the largest discrepancy between the whole minimiser and $x_b + K_\ast d$ computed from the gain formula was $1.1\times10^{-15}$ . The inverse Hessian had diagonal $(0.425475,\,0.180098,\,0.180098)$ , the first entry reproducing Step 4 and agreeing with $(I-K_\ast H)B$ to $1.7\times10^{-16}$ .

Step 6: the consistency check. Independently, $\tfrac12 d^{\mathsf T}S^{-1}d = \tfrac12\bigl[(1.30)(1.394277)+(0.10)(-0.663960)\bigr] = 0.873082$ , matching the minimum value above, as claim 4 of the 3D-Var statement requires. Its expectation is $\tfrac12 p = 1$ , and the observed $0.873$ is unremarkable for a $\tfrac12\chi^2(2)$ variable, so the assumed $B$ and $R$ are not contradicted by these two innovations.

Application(Where this is done)Original

Numerical weather prediction. Every operational forecast centre runs an analysis of this type every few hours. Optimal interpolation was the standard method from the 1970s until the 1990s, when 3D-Var replaced it, chiefly because the variational form handles observations that are non-local functions of the state — satellite radiances, which are weighted vertical integrals, and radar reflectivities — simply by putting the right $H$ into $J$ , whereas optimal interpolation needs the covariance between the analysis point and each such integral. The 4D-Var now used in several centres is the same cost function with the model's time evolution added as a constraint over an assimilation window, so that observations are compared with the forecast at their own observation times.

Reanalysis. The long gridded records used in climate research are produced by running a fixed forecast model and a fixed assimilation scheme over decades of archived observations. The point of freezing the scheme is that trends in the analysis should then reflect trends in the atmosphere rather than changes in the method, though changes in the observing network still leak in through $K_\ast$ , which depends on which observations exist.

Ocean, land surface and atmospheric composition. The same algebra analyses sea-surface temperature from satellite and buoy data, soil moisture, and ozone fields. In each case $B$ carries the physics: its correlation length sets how far a single observation is allowed to spread, and anisotropy in $B$ is how flow-dependent structure is injected into an otherwise static scheme.

Back to filtering. Cycling the analysis - use the forecast model to propagate $x_a$ forward to the next analysis time, call the result the new background, repeat - is exactly the loop of The Kalman filter, with the one difference that $B$ is prescribed rather than propagated as $P_k^{f} = F_kP_{k-1}F_k^{\mathsf T}+Q_k$ . Ensemble Kalman filters close that gap by estimating $B$ from the spread of a small ensemble of forecasts, and so sit between the two units: the gain formula of this unit, with a $B$ obtained from the recursion of the previous one.