\[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\det}{det} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\Signal}{S} \newcommand{\SignalNoise}{N} \newcommand{\AutoCov}{\gamma} \newcommand{\NoiseAutoCov}{\xi} \newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]} \newcommand{\Periodicity}{\tau} \newcommand{\Loadings}{\mathbf{b}} \newcommand{\FactorNoise}{\epsilon} \newcommand{\Uniqueness}{\mathbf{\psi}} \newcommand{\AltLoadings}{\mathbf{c}} \newcommand{\AltFactorNoise}{\zeta} \newcommand{\AltUniqueness}{\mathbf{\rho}} \newcommand{\DynMatrix}{\mathbf{a}} \newcommand{\DynNoise}{\eta} \newcommand{\DynVar}{\xi} \newcommand{\KalmanGain}{\mathbf{K}} \]

Previously on:

Last few weeks:
- Optimal linear predictors
- Restrictions on covariance functions to help estimate them
Last time:
- Estimating an un-observed latent signal from the observed random process
This time: More about latent states, leading to restrictions on covariance functions
We’ll begin by an (apparent) complete change of topic

Factor models

\(\vec{F} \sim (0, \mathbf{I}_q)\), a \(q\)-dimensional vector
- If \(\vec{F}\) isn’t centered, center it
- If \(\Var{\vec{F}} \neq \mathbf{I}_q\), rotate and scale to a new coordinate system where that is true
\(\vec{X}\) is a \(p\)-dimensional vector \[ X_i = \sum_{j=1}^{q}{\Loadings_{ij} F_j} + \FactorNoise_i \] with \(\Expect{\FactorNoise_i|\vec{F}} = 0\), \(\Cov{\FactorNoise_i, \FactorNoise_j|\vec{F}} = \Uniqueness_i \delta_{ij}\)
In matrix form, \[\begin{eqnarray} \vec{X} & = & \Loadings \vec{F} + \vec{\FactorNoise}\\ \end{eqnarray}\]
- \(\Loadings\) is a \(p \times q\) matrix
Implies a very specific type of covariance matrix, “low rank plus noise”: \[\begin{eqnarray} \Var{\vec{X}} & = & \Loadings \Var{\vec{F}} \Loadings^T + \Var{\vec{\FactorNoise}}\\ & = & \Loadings \Loadings^T + \mathbf{\Uniqueness}\\ & = & \mathrm{low\ rank} + \mathrm{diagonal} \end{eqnarray}\]

Estimating the factors from the observables

We want \(\LinearExpect{\vec{F}|\vec{X}}\), we know \(\Var{\vec{X}}\), so we need \(\Cov{\vec{X}, \vec{F}}\) \[\begin{eqnarray} \Cov{X_i, F_j} & = & \Cov{ \sum_{j=1}^{q}{\Loadings_{ij} F_j} + \FactorNoise_i, F_j}\\ & = & \sum_{j=1}^{q}\Loadings_{ij} \Cov{F_i, F_j}\\ & = & \Loadings_{ij}\\ \Cov{\vec{X}, \vec{F}} & = & \Loadings \end{eqnarray}\]
\(\therefore\) \[ \beta = (\Loadings \Loadings^T + \mathbf{\Uniqueness})^{-1} \Loadings \] and \[ \LinearExpect{\vec{F}|\vec{X}} = \beta^T \vec{X} \]

(Thomson 1936)

Predicting some observables from others via the factor

Say that \(\vec{Y} = \AltLoadings\vec{F} + \AltFactorNoise\), \(\Var{\AltFactorNoise} = \mathbf{\AltUniqueness}\)
We use \(\vec{X}\) to estimate \(\vec{F}\) and get \[ \LinearExpect{\vec{F}|\vec{X}} = \Loadings^T (\Loadings \Loadings^T + \mathbf{\Uniqueness})^{-1} \vec{X} \]
We use the estimate to predict \(\vec{Y}\) and get \[ \LinearExpect{\vec{Y}|\LinearExpect{\vec{F}|\vec{X}}} = \AltLoadings \LinearExpect{\vec{F}|\vec{X}} = \AltLoadings \Loadings^T (\Loadings \Loadings^T + \mathbf{\Uniqueness})^{-1} \vec{X} \]
Hooray, the factor variables let us predict things!

Predicting some observables from others without the factor

What’s the covariance between the two sets of observables? \[\begin{eqnarray} \Cov{X_i, Y_k} & = & \Cov{\sum_{j=1}^{q}{\Loadings_{ij} F_j + \FactorNoise_i}, \sum_{j^{\prime}=1}^{q}{\AltLoadings_{k j^{\prime}}F_j + \AltFactorNoise_k}}\\ & = & \sum_{j, j^{\prime}}{\Loadings_{ij} \AltLoadings_{k j^{\prime}} \Cov{F_j, F_{j^{\prime}}}}\\ & = & \sum_{j}{\Loadings_{ij} \AltLoadings_{kj}}\\ \Cov{\vec{X}, \vec{Y}} = \Loadings \AltLoadings^T \end{eqnarray}\]
What predictions do get for \(\vec{Y}\) from \(\vec{X}\)? \[ \LinearExpect{\vec{Y}|\vec{X}} = \AltLoadings \Loadings^T (\Loadings \Loadings^T + \mathbf{\Uniqueness})^{-1} \vec{X} \]
Boo, the factor variables don’t actually matter for the predictions

(Thomson 1936)

Estimating \(\Loadings\)

The classic procedure is alternating least squares
We have \(n\) observations of \(\vec{X}\), in a matrix \(\mathbf{X}\)
Now the value of observable \(i\) for data point \(r\) is \[ X_{ri} = \sum_{j=1}^{q}{\Loadings_{ij} F_{rj}} + \FactorNoise_{ri} \]
Start with a guess about \(\Loadings\)
- A classic, but not mandatory, starting point: PCA
Estimate \(\mathbf{F}\) from \(\mathbf{X}\) and \(\Loadings\)
- Using our linear predictor, or just plain linear regression
  - Our linear predictor is more efficient
Linearly regress \(X_{ri}\) on \(F_{rj}\) to get a new estimate for \(\Loadings_{ij}\)
Cycle to convergence

Summing up on factor models (for our purposes)

Assume \(\vec{X} = \Loadings\vec{F} + \vec{\FactorNoise}\)
Then the covariances are very structured
- There are other ways to generate that structure though (see backup)
and we can use the covariances to estimate the latent factor variables
and predict other observables

Factor models over time

Assume a dynamic factor model \[\begin{eqnarray} \vec{X}(t) & = & \Loadings \vec{F}(t) + \FactorNoise(t)\\ \vec{F}(t+1) & = & \DynMatrix \vec{F}(t) + \DynNoise(t+1) \end{eqnarray}\]
- \(\DynNoise(t)\) is sometimes called the innovation at time \(t\)
- For simplicity, \(\Var{\FactorNoise(t)} = \Uniqueness\) and \(\Var{\DynNoise(t)}= \DynVar\) (constant in \(t\))
Consequences: \[\begin{eqnarray} \Var{\vec{X}(t)} & = & \Loadings \Var{\vec{F}(t)} \Loadings^T + \Uniqueness\\ \Cov{\vec{X}(t+1), \vec{X}(t)} & = & \Cov{\Loadings(\DynMatrix \vec{F}(t) + \DynNoise(t+1)) + \FactorNoise(t+1), \Loadings \vec{F}(t) + \FactorNoise(t)}\\ & = & \Loadings \DynMatrix \Var{\vec{F}(t)} \Loadings^T\\ \Cov{\vec{X}(t+h), \vec{X}(t)} & = & \Loadings \DynMatrix^h \Var{\vec{F}(t)} \Loadings^T \end{eqnarray}\]
- Obvious simplifications if \(\Var{\vec{F}(t)} = \mathbf{I}_q\)

Kalman filter

State vector \(S(t)\)
- could be \(\vec{F}(t)\) or we could deal with some other model
Linear time evolution: \(S(t+1) = \DynMatrix S(t) + \DynNoise(t+1)\)
- \(\Var{\DynNoise(t)} = \DynVar\), \(\Cov{\DynNoise(t), \DynNoise(s)} = 0\)
Linear observation: \(X(t) = \Loadings S(t) + \FactorNoise(t)\)
- \(\Var{\FactorNoise(t)} = \Uniqueness\) (no longer necessarily diagonal), \(\Cov{\FactorNoise(t), \FactorNoise(s)} = 0\), \(\Cov{\FactorNoise(t), \DynNoise(s)} = 0\)
- Because everything’s really linear, best linear prediction is the conditional expectation
Start with some guess \(\hat{S}(t|t)\) i.e. using all obs. from times \(\leq t\)
Conditional mean \(\Expect{S(t+1)|\hat{S}(t|t)} = \DynMatrix \hat{S}(t|t)\)
- Let’s say this is \(\hat{S}(t+1|t)\)
Conditional mean for \(X\) is \(\Expect{X(t+1)|\hat{S}(t+1|t)} = \Loadings \hat{S}(t+1|t) = \Loadings \DynMatrix \hat{S}(t|t)\)
Update: for some Kalman gain matrix \(\KalmanGain_t\) \[ \hat{S}(t+1|t+1) = \hat{S}(t+1|t) + \KalmanGain_{t+1} \left(X(t+1) - \Expect{X(t+1)|\hat{S}(t+1|t)}\right) \]
- This has the form of the linear predictor, centered on the means…

Kalman filter (cont’d.)

What would be the best gain matrix? Optimal linear prediction says: \[ \KalmanGain_{t+1}= \Var{X(t+1)|S(t+1|t)}^{-1} \Cov{X(t+1), S(t+1|t)} \]
We (basically) worked out every term when we talked about dynamic factor models: \[ \KalmanGain_{t+1} = (\Loadings\Var{S(t+1|t)} \Loadings^T + \Uniqueness)^{-1} \Loadings \Var{S(t+1|t)} \]
We know how big the prediction errors of linear prediction are, and here that’s the conditional variance: \[ \Var{S(t+1|t+1)} = \left(\mathbf{I}_q - \KalmanGain_{t+1} \Loadings\right) \Var{S(t+1|t)} \]
- Some character-building matrix algebra omitted here
We know how the variance propagates: \[ \Var{S(t+1|t)} = \DynMatrix \Var{S(t|t)} \DynMatrix^T + \DynVar \]

Kalman filter (cont’d.)

Summing up

Start with initial estimate \(\hat{S}(t|t)\) and \(\Var{S(t|t)}\)
Predict new mean state \(\hat{S}(t+1|t) = \DynMatrix \hat{S}(t|t)\), with variance \(\Var{S(t+1|t)} = \DynMatrix \Var{S(t|t)} \DynMatrix^T + \DynVar\)
Predict new mean observable \(\Loadings \hat{S}(t+1|t)\)
Measure new observable \(X(t+1)\)
Compute Kalman gain \(\KalmanGain_t = (\Loadings \Var{S(t+1|t)} \Loadings^T + \Uniqueness)^{-1} \Loadings \Var{S(t+1|t)}\)
Updates: \[\begin{eqnarray*} \hat{S}(t+1|t+1) & = & \hat{S}(t+1|t) + \KalmanGain_t (X(t+1) - \Loadings \hat{S}(t+1|t))\\ \Var{S}(t+1|t+1) & = & (\mathbf{I}_q - \KalmanGain_t \Loadings) \Var{S(t+1|t)} \end{eqnarray*}\]

Some comments on this presentation of the Kalman filter

Recursive: we don’t need the whole observed past, just the previous estimate
- Do need to track the conditional variance
The Kalman filter is the optimal recursive linear filter for a linear state-space model
- For LSSM, one can show that the optimal recursive filter is the optimal filter over all
  - This is not obvious
Nothing needs to be Gaussian
- This used to be fairly widely understood (e.g., Liebelt (1967), pp. 161–165) but has become slightly arcane lore

What if we do assume Gaussians?

Assume \(S(0)\), \(\DynNoise\) and \(\FactorNoise\) are all Gaussians
- So uncorrelated \(\Rightarrow\) independent, etc.
\(X(t)\) and \(S(t)\) are then (jointly) Gaussian
Even \((X(1:n), S(0:n))\) is jointly Gaussian
\(S(t)|X(1:t-1))\) and \(S(t)|X(1:t))\) are also Gaussian
The Kalman filter gives us the mean and variance of \(S(t)|X(1:t)\)
In fact, \(X(t+1)|X(1:t)\) is also Gaussian, with mean \(\Loadings S(t+1|t)\) and variance \(\Loadings \Var{S(t+1|t)} \Loadings^T + \Uniqueness\)
This lets us calculate the likelihood of \(X(1:n)\)

Finding the parameters of a linear state-space model

Crude approach

Guess \(\Loadings\)
Until converged:
- Estimate \(S(t)\) from \(X(t)\) using \(\Loadings\)
- Estimate \(\Loadings\) from \(S(t), X(t)\) pairs
Estimate \(\DynMatrix\) from \(S(t), S(t+1)\) pairs
Estimate noise variances from residuals

Slightly less crude

Guess \(\Loadings\), \(\DynMatrix\)
Until converged:
- Estimate \(S(t)\) using Kalman filter (i.e., using \(X(t)\), \(\Loadings\) and \(\DynMatrix\))
- Estimate \(\Loadings\) by regressing \(X(t)\) on \(S(t)\)
- Estimate \(\DynMatrix\) by regressing \(S(t+1)\) on \(S(t)\)
- Estimate noise variances from residuals

More principled, if you believe things are Gaussian

Maximize the likelihood

What about space?

Kalman filter really needs a bottleneck, all dependence on the past is channeled through the present state
- More or less defines “state”
This does not make sense for purely spatial data
For space-time data, treat the whole spatial field at time \(t\) as one big vector of observables
- We’ll come back to more “local” filters when we look at prediction processes in a week
Q: Couldn’t we do something like a hidden Markov random field? A: Yes, but the way to do that will be part of the general theory for hidden Markov models we’ll cover next week

Some history lessons

The Kalman filter was introduced by Kalman (1960); Kalman and Bucy (1961)
Conceptually, this was one of the first places people clearly distinguished between a latent state, which describes the dynamics, and the more-or-less noisy observables (Raginsky 2024)
- Though the “state” and “observable” terminology was older in physics, e.g., Dirac (1935)
KF was very quickly adopted by the US aerospace industry (McGee and Schmidt 1985)
- “Where is the rocket? Where is it going?”
- Recursive form saved memory and computing when those were very expensive

Summing up

Factor models: \(p\) observables are linear in \(q < p\) latents
- Implies low-rank-plus-diagonal covariances
- Linear predictors let us estimate latent factors
Dynamic factor models: link observables over time by linear evolution of latents
DFMs \(\subset\) linear state space models (LSSM)
- Linear dynamical system, plus noise, in the state space
- Observables are linear in the state, plus more noise
Kalman filter is a recursive optimal linear state-estimator
- bonus: predictions for the observables
- drawback: need to keep updating a variance matrix, but oh well
- Can work for temporal or spatiotemporal data, purely spatial is ???
Next week: non-linear state space models

Backup: Alternative estimators for factor models

The most common approaches for factor models these days, with not-too-big-data and not-too-high dimension, is to assume that everything’s Gaussian and use maximum likelihood
- More statistically efficient than, e.g., alternating least squares if everything’s Gaussian
- More annoying to deal with, e.g., missing data than the alternating-least-squares approach I sketched
  - Missing data is important for many applications! (Feuerverger, He, and Khatri 2012)
There are some tricks which begin with an initial PCA for larger data

Backup: Other ways to get low-rank-plus-noise covariances

Mixture models: Every \(q\)-factor model has the same variance matrix as some mixture model with \((q+1)\) components, and vice versa (Bartholomew 1987)
Sampling: \(q\) latent variables but now \(p \ll q\), and each observable is the sum of a random sample of latents (Thomson 1916, 1919, 1939; Bartholomew, Deary, and Lawn 2009)
Many others (Maas et al. 2006)
Many psychologists have a bad habit of leaping from “this covariance matrix is close to being rank \(q\) (plus diagonal)” to “there are \(q\) continuous latent factors”

Backup: factor models vs. PCA

Geometry

PCA approximates the data vectors by a \(q\)-dimensional linear subspace
Factor models expect the data to lie in \(q\)-dimensional subspace plus noise

Variance matrix

PCA: \(\Var{\vec{X}} \approx \mathbf{w}^{(q)} \mathbf{\Lambda}^{(q)} (\mathbf{w}^{(q)})^T\)
- A low-rank matrix (when \(q < p\))
Factor model: \(\Var{\vec{X}} = \Loadings \Loadings^T + \Uniqueness\)
- Low rank plus diagonal
If \(\Uniqueness \approx 0\) then \(\Loadings \approx \mathbf{w}^{(q)} \sqrt{\mathbf{\Lambda}^{(q)}}\), but not in general
- Much confusion on this point, especially in psychology

Backup: Linear dynamical systems in multiple dimensions

\[ \vec{x}(t) = \DynMatrix\vec{x}(t-1) \]

Suppose the eigenvectors \(\vec{v}_1, \ldots \vec{v}_p\) of \(\DynMatrix\) form a basis
Then with the right coefficients \(c_1, \ldots c_p\), \[ \vec{x}(0) = \sum_{j=1}^{p}{c_j \vec{v}_j} \]

Linear dynamical systems in multiple dimensions

Dynamics are just multiplying: \[\begin{eqnarray} \vec{x}(t) & = & \DynMatrix\vec{x}(t-1)\\ & = & \DynMatrix^t \vec{x}(0)\\ & = & \DynMatrix^t \left(\sum_{j=1}^{p}{c_j \vec{v}_j}\right)\\ & = & \sum_{j=1}^{p}{c_j \DynMatrix^t \vec{v}_j}\\ & = & \sum_{j=1}^{p}{c_j \lambda_j^t \vec{v}_j} \end{eqnarray}\]

Eigenvalues determine the dynamics of a linear system

Try the easy case first: all eigenvalues \(\lambda_1, \ldots \lambda_p\) are real

\(\lambda_j > 1\): grow along that direction
\(0 \leq \lambda_j < 1\): shrink along that direction towards the origin \(\vec{0}\)
\(\lambda_j < -1\): flip around the origin, grow in that direction
\(-1 < \lambda_j \leq 0\): flip around the origin and shrink

Eigenvalues determine the dynamics of a linear system

Some eigenvalues can be complex
- The corresponding coefficients are also complex
These always come in complex-conjugate pairs
- so the coefficients are complex-conjugate pairs
The formula \(\vec{x}(t) = \sum_{j=1}^{p}{\lambda^t_j c_j \vec{v}_j}\) still works
- the imaginary parts always cancel exactly
Complex eigenvalues \(\Leftrightarrow\) rotation

Morals on linear, deterministic dynamical systems

Find the eigenvalues and eigenvectors
- \(|\lambda_j| < 1 \Rightarrow\) exponential decay along \(\vec{v}_j\)
- \(|\lambda_j| > 1 \Rightarrow\) exponential growth
- \(|\lambda_j| = 1 \Rightarrow\) eternal recurrence
- \(\mathrm{Im}(\lambda_j) \neq 0 \Rightarrow\) rotations in the space spanned by those eigenvectors
Higher-order dependence on the past \(\Rightarrow\) first-order dependence with extra memory variables

Adding on noise

If the eigenvalues are less than \(1\), \(\DynMatrix\vec{X}(t-1)\) has less variance than \(\vec{X}(t-1)\)
- \(\Var{\DynMatrix \vec{X}(t-1)} = \DynMatrix\Var{\vec{X}(t-1)}\DynMatrix^T\)
The innovation contributes some extra variance
- namely \(\Var{\vec{\DynNoise}(t)} = \DynVar\)
- In sum, \(\Var{\vec{X}(t)} = \DynMatrix\Var{\vec{X}(t-1)}\DynMatrix^T + \DynVar\)
To be stationary, the shrinkage has to exactly balance the new variance
- Stationarity requires \(\Var{\vec{X}(t)} = \Var{\vec{X}(t-1)}\) which just won’t be possible if any eigenvalues of \(\DynMatrix\) are \(\geq 1\) in magnitude
- A stationary solution is always possible if all the eigenvalues of \(\DynMatrix\) are \(< 1\) in magnitude

References

Bartholomew, David J. 1987. Latent Variable Models and Factor Analysis. New York: Oxford University Press.

Bartholomew, David J., Ian J. Deary, and Martin Lawn. 2009. “A New Lease on Life for Thomson’s Bonds Model of Intelligence.” Psychological Review 116:567–79. https://doi.org/10.1037/a0016262.

Dirac, P. A. M. 1935. Principles of Quantum Mechanics. Oxford: Clarendon Press.

Feuerverger, Andrey, Yu He, and Shashi Khatri. 2012. “Statistical Significance of the Netflix Challenge.” Statistical Science 27:202–31. https://doi.org/10.1214/11-STS368.

Kalman, R. E., and R. S. Bucy. 1961. “New Results in Linear Filtering and Prediction.” ASME Transactions, Journal of Basic Engineering 83D:95–108.

Kalman, Rudolf E. 1960. “A New Approach to Linear Filtering and Prediction Problems.” ASME Transactions, Journal of Basic Engineering 82D:35–50.

Liebelt, Paul B. 1967. An Introduction to Optimal Estimation. Reading, Massachusetts: Addison-Wesley.

Maas, Han L. J. van der, Conor V. Dolan, Raoul P. P. P. Grasman, Jelte M. Wicherts, Hilde M. Huizenga, and Maartje E. J. Raijmakers. 2006. “A Dynamical Model of General Intelligence: The Positive Manifold of Intelligence by Mutualism.” Psychological Review 113:842–61. https://doi.org/10.1037/0033-295X.113.4.842.

McGee, Leonard A., and Stanley F. Schmidt. 1985. “Discovery of the Kalman Filter as a Practical Tool for Aerospace and Industry.” 86847. NASA Technical Memorandum. https://ntrs.nasa.gov/citations/19860003843.

Raginsky, Maxim. 2024. “The State-Space Revolution in the Study of Complex Systems.” In Foundational Papers in Complexity Science, edited by David C. Krakauer, 1:449–57. Santa Fe, New Mexico: Santa Fe Institute.

Thomson, Godfrey H. 1916. “A Hierarchy Without a General Factor.” British Journal of Psychology 8:271–81. https://doi.org/10.1111/j.2044-8295.1916.tb00133.x.

———. 1919. “On the Cause of Hierarchical Order Among the Correlation Coefficients of a Number of Variates Taken in Pairs.” Proceedings of the Royal Society of London A 95:400–408. http://www.jstor.org/stable/93637.

———. 1936. “Some Points of Mathematical Technique in the Factorial Analysis of Ability.” Journal of Educational Psychology 27:37–54. https://doi.org/10.1037/h0062007.

———. 1939. The Factorial Analysis of Human Ability. Boston: Houghton Mifflin Company. http://www.archive.org/details/factorialanalysi032965mbp.

Factor Models, Linear State Space Models, and the Kalman Filter

Previously on:

Factor models

Estimating the factors from the observables

Predicting some observables from others via the factor

Predicting some observables from others without the factor

Estimating \(\Loadings\)

Summing up on factor models (for our purposes)

Factor models over time

Kalman filter

Kalman filter (cont’d.)

Kalman filter (cont’d.)

Summing up

Some comments on this presentation of the Kalman filter

What if we do assume Gaussians?

Finding the parameters of a linear state-space model

Crude approach

Slightly less crude

More principled, if you believe things are Gaussian

What about space?

Some history lessons

Summing up

Backup: Alternative estimators for factor models

Backup: Other ways to get low-rank-plus-noise covariances

Backup: factor models vs. PCA

Geometry

Variance matrix

Backup: Linear dynamical systems in multiple dimensions

Linear dynamical systems in multiple dimensions

Eigenvalues determine the dynamics of a linear system

Eigenvalues determine the dynamics of a linear system

Morals on linear, deterministic dynamical systems

Adding on noise

References