--- title: Factor Models, Linear State Space Models, and the Kalman Filter date: 17 September 2026 (Lecture 8) author: 36-740 output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{r, include=FALSE} # General set-up options library(knitr) opts_chunk$set(size="small", background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \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 {.smaller} - $\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 {.smaller} - 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-points-of-mathematical-technique] ## 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 {.smaller} - 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-points-of-mathematical-technique] ## 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 {.smaller} - 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.) {.smaller} - 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-optimal-estimation, 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 {.smaller} #### 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; @Kalman-Bucy - 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-on-Kalman] + Though the "state" and "observable" terminology was older in physics, e.g., @Dirac-on-qm - KF was very quickly adopted by the US aerospace industry [@McGee-Schmidt-discovery-of-kalman-filter] + "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! [@Statistical-significance-of-netflix] - 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-latent-variable-models] - Sampling: $q$ latent variables but now $p \ll q$, and each observable is the sum of a random sample of latents [@Thomson-hierarchy-without-factor; @Thomson-causes-of-hierarchy; @Thomson-factorial-analysis; @Bartholomew-Deary-Lawn-on-Thomson] - Many others [@van-der-Maas-et-al-mutualism] - 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