The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
80,282 characters
\thispagestyle{empty}
\enlargethispage{1in}
\newlength{\oldparindent}
\oldparindent=\parindent
\parindent=0.0in
\vspace*{0.1in}
\noindent{\Large\bf Conditionally linear, matrix normal state space models}
\vspace{0.05in}
\noindent{
Drew D. Creal\footnote{
Department of Economics, University of Illinois Urbana-Champaign;
[email removed].},
Marcelo C. Medeiros,\footnote{
Department of Economics, University of Illinois Urbana-Champaign;
[email removed].}
and
Rodrigo Sarlo,\footnote{
Department of Electrical Engineering, Pontifical Catholic University of Rio de Janeiro;
[email removed].}
}
\vspace{0.05in}
\noindent{This version: \today}
{
{\bf Abstract}
\medskip
We develop a class of linear state space models for matrix-valued time series data where the state is a latent matrix normal process. We derive matrix versions of the Kalman filter, log-likelihood, and smoother enabling estimation of the latent state matrix as well as the model's parameters. To conduct Bayesian inference, we provide algorithms that draw from the joint posterior distribution of the latent state matrices conditional on the observed data and parameters. We apply these methods to a large panel of U.S. macroeconomic time series across the 50 U.S.\ states. The proposed framework accommodates mixed-frequency data, heteroskedasticity, and outliers within a unified matrix-valued structure. Empirically, we find that a small number of latent factors captures the joint dynamics across states and variables, providing a parsimonious and scalable approach to modeling high-dimensional macroeconomic systems.
\medskip
\medskip
\noindent{\bf Keywords: } matrix normal distribution, Kalman filter and smoother, simulation smoothing.
\medskip
}
\setcounter{footnote}{0}
\noindent \thispagestyle{empty}\addtocounter{page}{-1}\newpage{}
\setcounter{equation}{0}
\section{Introduction}
Conditionally linear Gaussian state space models and the Kalman filter are central tools in modern time series analysis. They enable estimation and forecasting in a wide range of empirically important models, including vector autoregressions, structural time series models, and dynamic factor models; see, e.g., \cite{CappeMoulinesRyden(05)}, \cite{DurbinKoopman(12)}, and \cite{ShumwayStoffer(25)}.
At the same time, matrix-valued time series data have become increasingly common in economics, engineering, and the biological sciences; see \cite{Tsay(24)} for a recent survey. In many applications, the data naturally take the form of an $m \times n$ matrix $\mathbf{Y}_t$, where rows correspond to variables and columns correspond to related cross-sectional units. For example, in macroeconomics, one may observe multiple economic series across U.S. states, where each column corresponds to a state and each row corresponds to a variable such as employment or income.
Despite the growing importance of matrix-valued data, the matrix analogue of conditionally linear Gaussian state space models and their associated estimation methods remain only partially developed in the literature. Matrix-variate dynamic linear models originate with \cite{QuintanaWest(87)}, are treated within the general Bayesian forecasting framework of \cite{WestHarrison(97)}, and are further developed by \cite{CarvalhoWest(07)} and \cite{WangWest(09)}. \cite{CarvalhoWest(07)} derive the posterior mean and covariance of the smoothed state matrix, but only for a model in which every column shares the same known regression vector and evolution matrix -- there is no estimated loading matrix, and no row/column dimension mismatch for a companion-form representation to resolve. None of these papers derive the matrix-variate log-likelihood, a simulation smoother for drawing the state matrices jointly from their posterior distribution, or a precision sampler. Existing work therefore does not provide a unified framework for likelihood evaluation, state smoothing, simulation smoothing, and posterior inference for the more general class of matrix-valued state space models we consider.
Our main contribution is twofold. First, we introduce a companion-form representation for matrix normal processes that substantially generalizes existing matrix state space formulations (e.g., \cite{WangWest(09)}), accommodating a broad class of models, including bilinear dynamic factor models. Second, building on this representation, we derive the matrix Kalman filter, log-likelihood, and smoothing distributions. We also develop simulation smoothing and precision sampling methods for drawing the latent state matrix jointly from its posterior distribution. These extend the simulation smoother of \cite{DurbinKoopman(02)} and the precision-based sampler of \cite{ChanJeliazkov(09)} to the matrix setting. By preserving the matrix structure of the data throughout, the framework achieves substantial dimension reduction relative to vectorized representations while remaining fully likelihood-based. Vectorizing the observation and transition equations recovers an ordinary vector state space model, so the filtering and smoothing distributions we derive coincide exactly with those of that vector model. Our contribution is computing them without ever forming or factoring the resulting large Kronecker covariance, instead working directly with much smaller matrices tied to the row and column dimensions of the data and the state.
The framework readily accommodates additional features commonly used in empirical applications, such as time-varying parameters, heteroskedasticity, and mixture models. In particular, we extend the methods of \cite{GerlachCarterKohn(00)} and \cite{DoucetAndrieu(01)} to draw discrete latent states without conditioning on the continuous state matrix, enabling efficient Bayesian inference in matrix-valued state space models.
\subsection{Relation to the Literature}
\cite{CarvalhoWest(07)} introduce a matrix-normal dynamic linear model with graphically structured cross-sectional covariance, deriving both a forward filter and a backward smoother for the posterior moments of the state matrix. \cite{WangWest(09)} extend this to genuine matrix observations with a matrix state space representation and filtering algorithms. In both cases, however, every column shares the same known regression vector and evolution matrix, so there is no estimated loading matrix and no row/column dimension mismatch to resolve, and their smoother gives only the smoothed mean and covariance, not a joint draw of the state path. We generalize this to a companion-form representation that accommodates estimated, series-specific loading matrices and a general row/column dimension mismatch, and add likelihood evaluation, smoothing, simulation smoothing, and precision-based posterior simulation.
\cite{ChoukrounWeissBarItzhackOshman(06)} derive a more general matrix Kalman filter for a sum of bilinear terms with fully unrestricted, non-separable covariance matrices, obtained by vectorizing the system and applying the standard vector Kalman filter as a minimum-variance estimator. Their approach abandons the Kronecker-separable, matrix normal structure our framework relies on.
A growing literature develops matrix autoregressive models for matrix-valued time series (\cite{ChenXiaoYang(21)}, \cite{WangLiuChen2019}), typically estimated by moment-based or least-squares methods, with a large-scale Bayesian version in \cite{ChanQi(25)}. A related strand imposes low-rank, two-way factor structure on matrix data (\cite{YuanGaoHeHuangGuo(23)}, \cite{QinWangZhuShia(25)}). None of this casts the problem as a state space model; doing so is what gives us likelihood-based filtering, smoothing, and joint posterior simulation of the latent state.
Finally, our framework also builds on the classic Bayesian dynamic factor model literature (e.g., \cite{AguilarWest(00)}, \cite{LopesWest2004}), which treats the data as vector-valued. We extend this to matrix-valued data by keeping the matrix structure explicit throughout: the matrix normal distribution with separable covariance is what buys the dimension reduction and computational gains over a vectorized factor model.
\section{Matrix normal state space models} \label{model}
\subsection{Model}
Let $\mathbf{Y}_{t}$ denote an $m \times n$ matrix time series observed for $t=1,\ldots,T$. We study models that can be fit into the state space representation
\begin{eqnarray}
\mathbf{Y}_{t} & = & \mathbf{D}_{t} + \mathbf{Z}_{t}\mathbf{A}_{t}\mathbf{W}_{t}^{\top} + \mathbf{E}_{1t}, \qquad \mathbf{E}_{1t} \sim \text{MN}\left( \pmb{0},\mathbf{H}_{t},\mathbf{U}_{\mathbf{Y},t}\right), \label{obsrvation equation} \\
\mathbf{A}_{t+1} &=& \mathbf{C}_{t} + \mathbf{T}_{t}\mathbf{A}_{t} + \mathbf{R}_{t}\mathbf{E}_{2t}, \qquad \mathbf{E}_{2t} \sim \text{MN}\left( \pmb{0},\mathbf{Q}_{t},\mathbf{U}\right), \label{transition equation} \\
\mathbf{A}_{1} & \sim & \text{MN}\left( \mathbf{A}_{1|0},\mathbf{P}_{1|0},\mathbf{U}\right). \label{initial condition}
\end{eqnarray}
The state matrix $\mathbf{A}_{t}$ is $s \times r$ and is potentially unobserved or latent. The model requires that $n \geq r$. The matrices of shocks $\mathbf{E}_{1t} $ and $\mathbf{E}_{2t}$ have dimensions $m \times n$ and $q \times r$, respectively. Both have matrix normal (MN) distributions with mean zero and covariance matrices $\mathds{V}\left[ \text{vec}\left(\mathbf{E}_{1t} \right)\right] \ = \ \mathbf{U}_{\mathbf{Y},t} \otimes \mathbf{H}_{t}$ and $\mathds{V}\left[ \text{vec}\left(\mathbf{E}_{2t} \right)\right] \ = \ \mathbf{U} \otimes \mathbf{Q}_{t}$.
The system matrices are functions of a $k \times 1$ vector of parameters $\pmb{\theta}$ that need to be estimated. We take $\pmb{\theta}$ as known in Sections \ref{model}-\ref{MKFS} and discuss their estimation below.
The state space representation (\ref{obsrvation equation})-(\ref{initial condition}) is written in companion form enabling it to encompass a wide range of models, which we illustrate through several examples.
\textbf{Example \#1:} Consider a matrix version of the local level model
\begin{eqnarray*}
\mathbf{Y}_{t} &=& \mathbf{A}_{t} + \mathbf{E}_{1t}, \qquad \mathbf{E}_{1t} \sim \text{MN}\left( \pmb{0},\pmb{\Sigma},\mathbf{U}\right),\\
\mathbf{A}_{t+1} &=& \mathbf{A}_{t} + \mathbf{E}_{2t}, \qquad \mathbf{E}_{2t} \sim \text{MN}\left( \pmb{0},\pmb{\Omega},\mathbf{U}\right),
\end{eqnarray*}
which can be placed in state space form (\ref{obsrvation equation})-(\ref{transition equation}) by defining $\mathbf{D}_{t} = \pmb{0}$, $\mathbf{Z}_{t} = \mathbf{I}$, $\mathbf{W}_{t} = \mathbf{I}$, $\mathbf{H}_{t} = \pmb{\Sigma}$, $\mathbf{C}_{t} = \pmb{0}$, $\mathbf{T}_{t} = \mathbf{I}$, $\mathbf{R}_{t} = \mathbf{I}$, and $\mathbf{Q}_{t} = \pmb{\Omega}$. This model can be extended to include trends, seasonals, and cycles as in the literature on structural time series models; see e.g. \cite{DurbinKoopman(12)}.
\textbf{Example \#2:} Consider a matrix bi-linear dynamic factor model
\begin{eqnarray}
\mathbf{Y}_{t} &=& \pmb{\Lambda} \pmb{\mathcal{F}}_{t}\mathbf{W}^{\top} + \mathbf{E}_{1t}, \qquad \mathbf{E}_{1t} \sim \text{MN}\left( \pmb{0},\pmb{\Sigma}_{t},\mathbf{U}_{\mathbf{Y}}\right), \label{Example 2a} \\
\pmb{\mathcal{F}}_{t+1} &=& \pmb{\Phi}_{1}\pmb{\mathcal{F}}_{t} + \pmb{\Phi}_{2}\pmb{\mathcal{F}}_{t-1} + \mathbf{E}_{2t}, \qquad \mathbf{E}_{2t} \sim \text{MN}\left( \pmb{0},\pmb{\Omega}_{t},\mathbf{U}\right). \label{Example 2c}
\end{eqnarray}
The left factor loadings $\pmb{\Lambda}$ control the common factors across series while the right factor loadings $\mathbf{W}$ control the common factors across spatial units. This model can be placed in state space form by defining the state and its matrices as
\[ \mathbf{A}_{t} \ = \ \left(\begin{matrix}
\pmb{\mathcal{F}}_{t} \\
\pmb{\mathcal{F}}_{t-1}
\end{matrix}\right) \ \ \mathbf{T}_{t} \ = \ \left(\begin{matrix}
\pmb{\Phi}_{1} & \pmb{\Phi}_{2} \\
\mathbf{I} & \pmb{0}
\end{matrix}\right) \ \ \mathbf{C}_{t} \ = \ \left(\begin{matrix}
\pmb{0} \\
\pmb{0}
\end{matrix}\right) \ \ \mathbf{R}_{t} \ = \ \left(\begin{matrix}
\mathbf{I} \\
\pmb{0}
\end{matrix}\right) \]
and
$\mathbf{D}_{t} = \pmb{0}$, $\mathbf{Z}_{t} = \left( \pmb{\Lambda} \ \pmb{0}\right) $, $\mathbf{W}_{t} = \mathbf{W}$, $\mathbf{H}_{t} = \pmb{\Sigma}_{t}$, and $\mathbf{Q}_{t} = \pmb{\Omega}_{t}$. We estimate a version of this model in our empirical application.
\textbf{Example \#3:} Consider the matrix autoregressive moving average process model
\begin{eqnarray*}
\mathbf{Y}_{t} &=& \pmb{\Phi}_{1}\mathbf{Y}_{t-1}\pmb{\Upsilon}_{1}^{\top} + \ldots + \pmb{\Phi}_{\overline{p}}\mathbf{Y}_{t-\overline{p}}\pmb{\Upsilon}_{\overline{p}}^{\top} + \mathbf{E}_{t} + \pmb{\Theta}_{1}\mathbf{E}_{t-1}\pmb{\Gamma}_{1}^{\top} + \ldots + \pmb{\Theta}_{\overline{q}}\mathbf{E}_{t-\overline{q}}\pmb{\Gamma}_{\overline{q}}^{\top},
\end{eqnarray*}
with $\mathbf{E}_{t} \sim \text{MN}\left( \pmb{0},\pmb{\Sigma},\mathbf{U}\right)$.
Under the restrictions that $\pmb{\Gamma}_{\ell} = \mathbf{I}$ for all lags $\ell$, this model can be fit into the companion form (\ref{obsrvation equation})-(\ref{transition equation}). For example, if $\overline{p} = 2$ and $\overline{q} = 2$, we define $\mathbf{D}_{t} \ = \ \pmb{\Phi}_{1}\mathbf{Y}_{t-1}\pmb{\Upsilon}_{1}^{\top} + \pmb{\Phi}_{2}\mathbf{Y}_{t-2}\pmb{\Upsilon}_{2}^{\top}$
and
\[ \mathbf{A}_{t} \ = \ \left(\begin{matrix}
\mathbf{E}_{t} \\
\mathbf{E}_{t-1} \\
\mathbf{E}_{t-2}
\end{matrix}\right) \ \ \mathbf{T}_{t} \ = \ \left(\begin{matrix}
\pmb{0} & \pmb{0} & \pmb{0} \\
\mathbf{I} & \pmb{0} & \pmb{0} \\
\pmb{0} & \mathbf{I} & \pmb{0}
\end{matrix}\right) \ \ \mathbf{C}_{t} \ = \ \left(\begin{matrix}
\pmb{0} \\
\pmb{0} \\
\pmb{0}
\end{matrix}\right) \ \ \mathbf{R}_{t} \ = \ \left(\begin{matrix}
\mathbf{I} \\
\pmb{0} \\
\pmb{0}
\end{matrix}\right). \]
The remaining matrices are $\mathbf{Z}_{t} \ = \ \left(\begin{matrix}
\mathbf{I} & \pmb{\Theta}_{1} & \pmb{\Theta}_{2}
\end{matrix}\right)$, $\mathbf{H}_{t} = \pmb{0}$, $\mathbf{W}_{t} = \mathbf{I}$, and $\mathbf{Q}_{t} = \pmb{\Sigma}$. If the dynamics of $\mathbf{Y}_{t}$ follow a general MARMA$(\overline{p},\overline{q})$, the observations $\mathbf{Y}_{t}$ cannot be placed in the state matrix unless the additional restrictions that $\pmb{\Upsilon}_{j} = \mathbf{I}$ are imposed for all lags $j$. Right multiplication destroys the matrix normal structure. The transition equation propagates the state by left multiplication only, so every lag must share the same right scale matrix, ruling out lag-specific $\pmb{\Gamma}_{\ell}$ or $\pmb{\Upsilon}_{j}$.
An important special case of the model is when the system matrices of the transition equation (\ref{transition equation}) are time-invariant and the transition matrix $\mathbf{T}$ has all eigenvalues less than one in modulus. Then, the stochastic process $ p\left( \mathbf{A}_{t+1}|\mathbf{A}_{t},\pmb{\theta}\right)$ has a stationary distribution that is matrix normal $p\left( \mathbf{A}_{t}|\pmb{\theta}\right) \ = \ \text{MN}\left( \overline{\mathbf{A}},\overline{\mathbf{P}},\mathbf{U}\right) $ with mean matrix $\overline{\mathbf{A}} \ = \ \left( \mathbf{I}_{s} - \mathbf{T}\right)^{-1}\mathbf{C}$ and left scale matrix $\text{vec}\left(\overline{\mathbf{P}}\right) \ = \ \left( \mathbf{I}_{s^{2}} - \mathbf{T}\otimes\mathbf{T}\right)^{-1}\text{vec}\left( \mathbf{R}\mathbf{Q}\mathbf{R}^{\top}\right)$.
The stationary distribution is often chosen for the initial condition of the model.
In order for the filtering and smoothing distributions to be known in closed form, the measurement and transition equations must have a matrix normal distribution with a common right scale matrix $\mathbf{U}$ at each date, for two reasons. First, a sum of matrix normals is itself matrix normal only when one pair of scale matrices is proportional. The Kronecker factorization $\mathbf{P}\otimes\mathbf{U}$ is identified only up to such a rescaling. Therefore, we hold $\mathbf{U}$ fixed and shared across dates to keep the joint distribution $p\left( \mathbf{A}_{1},\ldots,\mathbf{A}_{T}|\pmb{\theta}\right)$ matrix normal.
Secondly, the proof of the matrix Kalman filter and related algorithms relies on a key property of the matrix normal distribution, which is Theorem 2.3.12 in \cite{GuptaNagar(00)}. We re-state it using our notation.
\begin{lemma}\label{Theorem 1} Let $\mathbf{X}$ be an $m \times n$ matrix with distribution $\mathbf{X} \ \sim \ \text{MN}\left( \mathbf{M},\mathbf{P},\mathbf{U}\right)$, where $\mathbf{P}$ is $m\times m$ and $\mathbf{U}$ is $n\times n$. Consider a partition of the matrices as
\[ \mathbf{X} \ = \ \left( \begin{matrix}
\mathbf{X}_{1} \\
\mathbf{X}_{2}
\end{matrix}\right) \qquad \mathbf{M} \ = \ \left(\begin{matrix}
\mathbf{M}_{1} \\
\mathbf{M}_{2}
\end{matrix}\right) \qquad \mathbf{P} \ = \ \left(\begin{matrix}
\mathbf{P}_{11} & \mathbf{P}_{12} \\
\mathbf{P}_{21} & \mathbf{P}_{22}
\end{matrix}\right) \]
such that $\mathbf{X}_{1}$ is $m_{1} \times n$ and $\mathbf{X}_{2}$ is $m_{2} \times n$ with $m = m_{1} + m_{2}$.
Then, the conditional distribution $p\left(\mathbf{X}_{2}|\mathbf{X}_{1}\right)$ is matrix normal $\text{MN}\left( \mathbf{M}_{2|1},\mathbf{P}_{2|1},\mathbf{U}\right)$ with
\[\mathbf{M}_{2|1} \ = \ \mathbf{M}_{2} + \mathbf{P}_{21}\mathbf{P}_{11}^{-1}\left( \mathbf{X}_{1} - \mathbf{M}_{1}\right), \qquad
\mathbf{P}_{2|1} \ = \ \mathbf{P}_{22} - \mathbf{P}_{21}\mathbf{P}_{11}^{-1}\mathbf{P}_{12}. \]
\end{lemma}
\vskip -0.2cm
The Kalman filter and related algorithms are a recursive application of this lemma.
\subsection{Collapsing the observation equation when $n > r$}
When the column dimension of $\mathbf{Y}_t$ exceeds that of the state matrix $\mathbf{A}_t$, i.e.\ $n > r$, the measurement and transition equations in (\ref{obsrvation equation})-(\ref{transition equation}) do not share a common right scale matrix. The measurement disturbances have right scale $\mathbf{U}_{\mathbf{Y},t}$ while the state disturbances have right scale
$\mathbf{U}$. This prevents direct application of Lemma \ref{Theorem 1}.
Our solution is to right multiply the original observation equation (\ref{obsrvation equation}) by a matrix $\mathbf{J}_{t} = \left( \mathbf{J}_{t}^{*,\top} \ \mathbf{J}_{t}^{+,\top} \right)^{\top}$ that is a function of $\mathbf{W}_{t}$ and $\mathbf{U}_{\mathbf{Y},t}$. The transformation $\mathbf{J}_{t}$ and the scale matrices $\mathbf{U}_{\mathbf{Y},t}$ and $\mathbf{U}$ must be jointly specified
to satisfy the conditions
\begin{eqnarray}
\mathbf{J}_{t}^{*}\mathbf{W}_{t} &=& \mathbf{I}_{r}, \label{condition 1} \\
\mathbf{J}_{t}^{+}\mathbf{W}_{t} &=& \pmb{0}, \label{condition 1b} \\
\mathbf{J}_t^* \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{*,\top} &=& \mathbf{U}, \label{condition 2} \\
\mathbf{J}_t^* \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{+,\top} &=& \pmb{0}. \label{condition 3}
\end{eqnarray}
and $\mathbf{J}_{t}$ must be full rank. Condition (\ref{condition 1b}) ensures that the transformed data $\mathbf{Y}_{t}^{+} \ \equiv \ \mathbf{Y}_{t}\mathbf{J}_{t}^{+,\top}$ does not depend on the state $\mathbf{A}_t$. Condition (\ref{condition 2}) fixes the right scale matrix of $\mathbf{Y}_{t}^{*} \ \equiv \ \mathbf{Y}_{t}\mathbf{J}_{t}^{*,\top}$ to equal $\mathbf{U}$. Condition (\ref{condition 3}) ensures that $\mathbf{Y}_t^*$ and $\mathbf{Y}_t^+$ are independent.
Under these conditions, the transformations
create a new set of observation equations
\begin{eqnarray}
\mathbf{Y}_{t}^{*} & = & \mathbf{D}_{t}^{*} + \mathbf{Z}_{t}\mathbf{A}_{t} + \mathbf{E}_{1t}^{*}, \qquad \mathbf{E}_{1t}^{*} \sim \text{MN}\left( \pmb{0},\mathbf{H}_{t},\mathbf{U}\right), \label{obs equation 1} \\
\mathbf{Y}_{t}^{+} & = & \mathbf{D}_{t}^{+} + \mathbf{E}_{1t}^{+}, , \qquad \mathbf{E}_{1t}^{+} \sim \text{MN}\left( \pmb{0},\mathbf{H}_{t},\pmb{\Psi}_t\right), \label{obs equation 2}
\end{eqnarray}
where $\mathbf{D}_{t}^{*} = \mathbf{D}_{t}\mathbf{J}_{t}^{*,\top}$, $\mathbf{D}_{t}^{+} = \mathbf{D}_{t}\mathbf{J}_{t}^{+,\top}$, and $\pmb{\Psi}_t = \mathbf{J}_t^+ \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{+,\top}$ is a $(n - r) \times (n-r)$ positive definite matrix. The matrix $\pmb{\Psi}_t$ will typically contain unknown parameters. By construction, $\mathbf{Y}_t^+$ contains the information about $\pmb{\Psi}_t$.
In order to satisfy conditions (\ref{condition 1})-(\ref{condition 3}), a researcher must first decide how to model $\mathbf{U}_{\mathbf{Y},t}$ and the identifying restrictions they want to impose on the model. Once these are chosen, the matrix $\mathbf{J}_t^*$ is pinned down uniquely, as the generalized least squares projection onto the column space of $\mathbf{W}_t$,
\begin{eqnarray}
\mathbf{J}_t^* &=& \left( \mathbf{W}_t^\top \mathbf{U}_{\mathbf{Y},t}^{-1} \mathbf{W}_t \right)^{-1} \mathbf{W}_t^\top \mathbf{U}_{\mathbf{Y},t}^{-1}. \label{Jstar GLS}
\end{eqnarray}
$\mathbf{J}_t^+$ can be any $(n-r)\times n$ matrix whose rows are orthogonal to the columns of $\mathbf{W}_t$, so that $\mathbf{Y}_t^+\equiv\mathbf{Y}_t\mathbf{J}_t^{+,\top}$ captures the part of $\mathbf{Y}_t$ that the state matrix does not explain. Unlike $\mathbf{J}_t^*$, the matrix $\mathbf{J}_t^+$ is not pinned down by the conditions. Any basis for this leftover space works, and different choices only relabel the coordinates of $\mathbf{Y}_t^+$, leaving $\mathbf{J}_t^*$ and $\pmb{\Psi}_t$ untouched. After discussing identification, we give concrete examples of how to parameterize $\mathbf{U}_{\mathbf{Y},t}$.
The idea of collapsing the measurement equation in a linear state space model into a lower dimensional representation was introduced by \cite{JungbackerKoopman(15)}. In their setting, the observation vector is collapsed (by left multiplication) to reduce the computational burden.
In our setting, collapsing the column dimension of the state space model is required in order to apply Lemma \ref{Theorem 1}. We show in the online appendix that their left collapse can be applied simultaneously with our column collapse, further reducing the dimension of the filter.
\subsection{Identification and choice of parameterization}
The model (\ref{obsrvation equation})--(\ref{transition equation}) is invariant to rotations $\mathbf{A}_t \to \mathbf{R}_{\ell} \mathbf{A}_t \mathbf{R}_{r}^{-1}$, $\mathbf{Z}_t \to \mathbf{Z}_t \mathbf{R}_{\ell}^{-1}$, and $\mathbf{W}_t \to \mathbf{W}_t \mathbf{R}_{r}^{\top}$ for any nonsingular $\mathbf{R}_{\ell},\mathbf{R}_{r}$, leaving the likelihood unchanged:
\begin{eqnarray*}
\left(\mathbf{Z}_t\mathbf{R}_\ell^{-1}\right)\left(\mathbf{R}_\ell\mathbf{A}_t\mathbf{R}_r^{-1}\right)\left(\mathbf{W}_t\mathbf{R}_r^\top\right)^\top &=& \mathbf{Z}_t\mathbf{A}_t\mathbf{W}_t^\top.
\end{eqnarray*}
The transition equation and initial condition are invariant under the same reparameterization provided $\mathbf{T}_t \to \mathbf{R}_\ell\mathbf{T}_t\mathbf{R}_\ell^{-1}$, $\mathbf{C}_t \to \mathbf{R}_\ell\mathbf{C}_t\mathbf{R}_r^{-1}$, $\mathbf{R}_t \to \mathbf{R}_\ell\mathbf{R}_t$, $\mathbf{U} \to \left(\mathbf{R}_r^{-1}\right)^\top\mathbf{U}\mathbf{R}_r^{-1}$, $\mathbf{P}_{1|0}\to\mathbf{R}_\ell\mathbf{P}_{1|0}\mathbf{R}_\ell^\top$, and $\mathbf{A}_{1|0}\to\mathbf{R}_\ell\mathbf{A}_{1|0}\mathbf{R}_r^{-1}$. Both the conditional mean $\mathbf{T}_t\mathbf{A}_t+\mathbf{C}_t$ and the transition equation's covariance matrix $\mathbf{U}\otimes\mathbf{R}_t\mathbf{Q}_t\mathbf{R}_t^\top$ transform consistently under this substitution. Identification therefore requires restrictions on both the row and column spaces of the state and on the scale matrices. A convenient choice imposes orthonormality, $\mathbf{Z}_t^\top \mathbf{Z}_t = \mathbf{I}$ and $\mathbf{W}_t^\top \mathbf{W}_t = \mathbf{I}$, and restricts $\pmb{\Omega}$ and $\mathbf{U}$ to be diagonal with ordered eigenvalues. These restrictions eliminate ordering and scaling indeterminacies, though not all rotational indeterminacy.
Two further normalizations are standard here, as in any factor model identified up to rotation. The sign of each factor is fixed by convention rather than by the model. Flipping the sign of a column of $\mathbf{Z}_t$ (or $\mathbf{W}_t$) together with the corresponding row (or column) of $\mathbf{A}_t$ leaves the orthonormality restrictions and ordered scale matrices unchanged. A researcher must designate one observed series as an anchor for each factor and flip the sign of the entire column whenever that anchor's loading comes out negative. A second normalization applies only if two eigenvalues of $\pmb{\Omega}$ or $\mathbf{U}$ are exactly tied, in which case any rotation within that pair leaves the restrictions unchanged. This does not bind once the corresponding eigenvalues are well separated, which is easy to check directly.
An alternative identifying assumption is to impose a triangular structure on the factor loadings by requiring that $\mathbf{Z}_t$ or $\mathbf{W}_t$ are each lower triangular with positive diagonal elements as in \cite{GewekeZhou(96)}, while allowing the scale matrices in the transition equation to be unrestricted. Although this achieves identification, the resulting factors depend on which variables are ordered first in $\mathbf{Y}_t$. Relabeling the rows changes which loadings are restricted to zero, an arbitrary dependence absent a natural ordering among the series.
Finally, the Kronecker structure implies a scalar identification problem. For any scalar $c > 0$, rescaling $\mathbf{U} \to c\mathbf{U}$, $\mathbf{Q}_t \to c^{-1}\mathbf{Q}_t$, $\mathbf{H}_t \to c^{-1}\mathbf{H}_t$, and $\pmb{\Psi}_t \to c\pmb{\Psi}_t$ leaves $\mathbf{U}\otimes\mathbf{Q}_t$, $\mathbf{U}\otimes\mathbf{H}_t$, and $\pmb{\Psi}_t\otimes\mathbf{H}_t$ unchanged. Since $\mathbf{H}_t$ is shared between $\mathbf{Y}_t^*$ and $\mathbf{Y}_t^+$, this is a single redundancy. Fixing the scale of a single parameter in any of the scale matrices solves the problem. For example, consider the matrix $\mathbf{U}$, we could set $U_{1,1}=1$, $\left| \mathbf{U}\right| = 1$, or $\text{tr}\left( \mathbf{U}\right)/r =1$, which then pins down $\mathbf{Q}_t$, $\mathbf{H}_t$, and $\pmb{\Psi}_t$ as well.
\paragraph{Example 1i: Generalized least squares parameterization.}
Here, a researcher specifies the scale matrix $\mathbf{U}_{\mathbf{Y},t}$ as a set of estimable parameters, and $\mathbf{J}_t^*$ is simply (\ref{Jstar GLS}), following \cite{JungbackerKoopman(15)}. To complete the transformation we need a complement $\mathbf{J}_t^+$ satisfying $\mathbf{J}_t^*\mathbf{U}_{\mathbf{Y},t}\mathbf{J}_t^{+,\top}=\pmb{0}$ with $\left( \mathbf{J}_t^{*,\top}, \mathbf{J}_t^{+,\top} \right)$ full rank. Any such choice works, and we fix the scale by imposing $\left| \mathbf{J}_t^{+} \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{+,\top} \right| = 1$. The state's own scale matrix $\mathbf{U}$ is not estimated separately. It is implied by $\mathbf{W}_t$ and $\mathbf{U}_{\mathbf{Y},t}$ through $\mathbf{U} = \left(\mathbf{W}_t^\top \mathbf{U}_{\mathbf{Y},t}^{-1} \mathbf{W}_t\right)^{-1}$, as is $\pmb{\Psi}_{t} = \mathbf{J}_t^{+} \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{+,\top}$. Identification requires restricting $\mathbf{W}_t^\top \mathbf{U}_{\mathbf{Y},t}^{-1} \mathbf{W}_t$, typically to be diagonal with ordered entries, which removes rotational indeterminacy in the column space.
\paragraph{Example 2i: Orthogonal parameterization.}
A simpler alternative is to restrict $\mathbf{W}_t$ itself to have orthonormal columns, $\mathbf{W}_t^\top \mathbf{W}_t = \mathbf{I}_r$, at the cost of imposing structure on $\mathbf{U}_{\mathbf{Y},t}$. Rather than estimating $\mathbf{U}_{\mathbf{Y},t}$ directly, a researcher estimates $\mathbf{U}$ and $\pmb{\Psi}_t$ separately, and $\mathbf{U}_{\mathbf{Y},t}$ is built from them. Let $\mathbf{W}_{t,\perp}$ be an orthonormal basis for the complement of $\mathbf{W}_t$. Setting $\mathbf{J}_t^{*,\top} = \mathbf{W}_t$ and $\mathbf{J}_t^{+,\top} = \mathbf{W}_{t,\perp}$ satisfies conditions (\ref{condition 1})--(\ref{condition 3}) whenever the right scale matrix takes the form
\begin{eqnarray}
\mathbf{U}_{\mathbf{Y},t} &=& \mathbf{W}_t \mathbf{U} \mathbf{W}_t^\top + \mathbf{W}_{t,\perp} \pmb{\Psi}_t \mathbf{W}_{t,\perp}^\top, \label{U param}
\end{eqnarray}
with both $\mathbf{U}$ and $\mathbf{W}_{t,\perp} \pmb{\Psi}_t \mathbf{W}_{t,\perp}^\top$ estimated freely. Identification is achieved by imposing that $\mathbf{U}$ is diagonal with ordered entries, eliminating rotational indeterminacy.
Example 2i isn't a different rule for building $\mathbf{J}_t$. It is (\ref{Jstar GLS}) evaluated at a value of $\mathbf{U}_{\mathbf{Y},t}$ that happens to be block diagonal in the $\left(\mathbf{W}_t,\mathbf{W}_{t,\perp}\right)$ basis.
\section{Algorithms} \label{MKFS}
Throughout this section, we assume $\mathbf{J}_t$ is full rank and satisfies conditions (\ref{condition 1})-(\ref{condition 3}). All proofs are contained in the online appendix.
\subsection{Filtering and one-step ahead prediction}
\begin{proposition}
Consider the state space model (\ref{obsrvation equation})-(\ref{initial condition}). For each date $t=1,\ldots,T$, the filtering and one-step ahead predictive distributions are matrix normal
\[ p\left( \mathbf{A}_{t}|\mathbf{Y}_{1:t},\pmb{\theta}\right) \ = \ \text{MN}\left( \mathbf{A}_{t|t},\mathbf{P}_{t|t},\mathbf{U}\right) \qquad p\left( \mathbf{A}_{t+1}|\mathbf{Y}_{1:t},\pmb{\theta}\right) \ = \ \text{MN}\left( \mathbf{A}_{t+1|t},\mathbf{P}_{t+1|t},\mathbf{U}\right)\]
with parameter matrices that can be calculated recursively forwards in time as
\begin{eqnarray}
\mathbf{V}_{t} &=& \mathbf{Y}_{t}\mathbf{J}_{t}^{*,\top} - \mathbf{D}_{t}\mathbf{J}_{t}^{*,\top} - \mathbf{Z}_{t}\mathbf{A}_{t|t-1}, \label{KF: V} \\
\mathbf{F}_{t} &=& \mathbf{Z}_{t}\mathbf{P}_{t|t-1} \mathbf{Z}_{t}^{\top} + \mathbf{H}_{t}, \\
\mathbf{K}_{t} &=& \mathbf{P}_{t|t-1}\mathbf{Z}_{t}^{\top}\mathbf{F}_{t}^{-1}, \label{Kalman gain} \\
\mathbf{A}_{t|t} & = & \mathbf{A}_{t|t-1} + \mathbf{K}_{t}\mathbf{V}_{t}, \\
\mathbf{P}_{t|t} & = & \mathbf{P}_{t|t-1} - \mathbf{K}_{t} \mathbf{Z}_{t}\mathbf{P}_{t|t-1}, \\
\mathbf{A}_{t+1|t} &=& \mathbf{T}_{t}\mathbf{A}_{t|t} + \mathbf{C}_{t}, \\
\mathbf{P}_{t+1|t} &=& \mathbf{T}_{t}\mathbf{P}_{t|t}\mathbf{T}_{t}^{\top} + \mathbf{R}_{t}\mathbf{Q}_{t}\mathbf{R}_{t}^{\top}. \label{Pred}
\end{eqnarray}
\end{proposition}
\vskip -0.5cm
These recursions are the same as the standard Kalman filter but where $\mathbf{V}_{t}, \mathbf{A}_{t|t}$ and $\mathbf{A}_{t|t-1}$ are matrix valued instead of vectors.
The matrix $\mathbf{U}$ does not enter the recursions but it does enter the likelihood function, see sub-section \ref{loglikelihood}.
In many empirical studies, the Kalman filter is primarily used to evaluate the likelihood or for simulation smoothing. In this case, the prediction form of the Kalman filter is more computationally efficient because it avoids calculating the filtered values. It redefines the matrix $\mathbf{K}_{t} \ = \ \mathbf{T}_{t}\mathbf{P}_{t|t-1}\mathbf{Z}_{t}^{\top}\mathbf{F}_{t}^{-1}$ and replaces (\ref{Kalman gain})-(\ref{Pred}) with
\[ \mathbf{A}_{t+1|t} \ = \ \mathbf{T}_{t}\mathbf{A}_{t|t-1} + \mathbf{C}_{t} + \mathbf{K}_{t}\mathbf{V}_{t}, \qquad
\mathbf{P}_{t+1|t} \ = \ \mathbf{T}_{t}\mathbf{P}_{t|t-1}\mathbf{L}_{t}^{\top} + \mathbf{R}_{t}\mathbf{Q}_{t}\mathbf{R}_{t}^{\top}, \]
where $\mathbf{L}_{t} = \mathbf{T}_{t} - \mathbf{K}_{t}\mathbf{Z}_{t}$.
\subsection{Computational savings: Kronecker product assumption}
Time series models that fit into the state space representation (\ref{obsrvation equation})-(\ref{initial condition}) are admittedly restricted because they assume a (conditional) Kronecker structure on their covariance matrix. This structure yields substantial computational savings. A naive vectorized Kalman filter must invert an $mn\times mn$ matrix at every date, which is $O\left((mn)^{3}\right)$. A sufficiently sophisticated vector implementation could in principle recover some of this cost via Woodbury-type updates, but existing Bayesian matrix-factor implementations do not attempt this and instead vectorize directly (\cite{ChanZhang(24)}, \cite{BarigozziTrapin(25)}). Because the recursions $\mathbf{A}_{t|t},\mathbf{P}_{t|t}$ of Proposition 1 depend only on $\mathbf{Z}_t,\mathbf{H}_t,\mathbf{F}_t$ and never $\mathbf{U}$, the matrix Kalman filter instead inverts only the $m\times m$ matrix $\mathbf{F}_t$ which is $O\left( m^{3}\right)$ and independent of $n$ and $r$. The matrix $\mathbf{U}$ enters only the log-likelihood, requiring a single, cheap $O\left( r^{3}\right)$ inversion. In our application ($m=12$, $n=50$, $r=11$), this is $1{,}728$ versus $2.16\times10^{8}$ operations per date. Using the approach of \cite{JungbackerKoopman(15)} to collapse the row dimension as well (see the online appendix) reduces the cost further to $O\left( s^{3}\right)$ when $\mathbf{Z}_t$ has rank $s<m$. The same principle - exploiting a Kronecker-separable covariance matrix to avoid a large joint matrix inversion -- underlies the natural-conjugate priors used to make large Bayesian VARs computationally tractable; see \cite{CarrieroClarkMarcellino(16)}, \cite{CarrieroClarkMarcellino(19)}, and \cite{Chan(20)}.
\subsection{Log-likelihood function} \label{loglikelihood}
We now derive the log-likelihood for the model, when $\mathbf{J}_{t}$ rotates the data from $\mathbf{Y}_{t}$ to $\left[ \mathbf{Y}_{t}^{*} \ \mathbf{Y}_{t}^{+}\right]$ leaving observation equations (\ref{obs equation 1})-(\ref{obs equation 2}).
\begin{proposition} \label{proposition loglikelihood}
The log-likelihood of the state space model (\ref{obsrvation equation})-(\ref{initial condition}) is
\begin{eqnarray*}
\log p\left( \mathbf{Y}_{1:T}|\pmb{\theta}\right) &=& \log p\left( \mathbf{Y}_{1:T}^{*}|\pmb{\theta}\right) + \log p\left( \mathbf{Y}_{1:T}^{+}|\pmb{\theta}\right) + m\sum_{t=1}^{T}\log\left| \mathbf{J}_{t}\right|,
\end{eqnarray*}
where $\left| \mathbf{J}_{t}\right|$ is the Jacobian of the transformation from $\mathbf{Y}_{t}$ to $\left( \mathbf{Y}_{t}^{*,\top},\mathbf{Y}_{t}^{+,\top}\right)^{\top}$. The log-likelihood contribution for the transformed data $\mathbf{Y}_{t}^{*}$ is
\begin{eqnarray*}
\log p\left( \mathbf{Y}_{1:T}^* \mid \pmb{\theta}\right)
&=& -\frac{Tmr}{2} \log(2\pi)
- \frac{r}{2} \sum_{t=1}^T \log |\mathbf{F}_t|
- \frac{mT}{2} \log |\mathbf{U}| \\
& &
- \frac{1}{2} \sum_{t=1}^T \mathrm{tr}
\left(
\mathbf{F}_t^{-1} \mathbf{V}_t \mathbf{U}^{-1} \mathbf{V}_t^{\top}
\right),
\end{eqnarray*}
with $\mathbf{V}_t$ and $\mathbf{F}_t$ obtained from the matrix Kalman filter applied to $\mathbf{Y}_t^*$, and the log-likelihood contribution of $\mathbf{Y}_{t}^+$ is
\begin{eqnarray*}
\log p\left( \mathbf{Y}_{1:T}^+ \mid \pmb{\theta}\right)
&=& -\frac{Tm(n-r)}{2} \log(2\pi)
- \frac{n-r}{2} \sum_{t=1}^T \log \left| \mathbf{H}_t\right|
- \frac{m}{2}\sum_{t=1}^{T}\log |\pmb{\Psi}_{t}|\\
& & - \frac{1}{2} \sum_{t=1}^T \mathrm{tr}
\left(
\mathbf{H}_t^{-1} \left[\mathbf{Y}_t^+ - \mathbf{D}_{t}^{+}\right]\pmb{\Psi}_{t}^{-1}\left[\mathbf{Y}_t^{+}-\mathbf{D}_{t}^{+}\right]^{\top}
\right).
\end{eqnarray*}
The decomposition is also invariant to the choice of basis for the complement $\mathbf{J}_t^+$. The log-likelihood does not depend on which valid $\mathbf{J}_t^+$ is used.
\end{proposition}
Although the log-likelihood is invariant, $\pmb{\Psi}_t = \mathbf{J}_t^+\mathbf{U}_{\mathbf{Y},t}\mathbf{J}_t^{+,\top}$ itself is not. Replacing $\mathbf{J}_t^+$ with $\mathbf{B}\mathbf{J}_t^+$ for any nonsingular $\mathbf{B}$ changes $\pmb{\Psi}_t$ to $\mathbf{B}\pmb{\Psi}_t\mathbf{B}^\top$, with the resulting change in the likelihood contribution of $\mathbf{Y}_t^+$ exactly offset by the change in the Jacobian term (see the online appendix).
When $n=r$, no collapsing is needed, $\mathbf{W}_t=\mathbf{I}_n$, $\mathbf{J}_t=\mathbf{I}_n$, and $\mathbf{Y}_t^+$ is an empty matrix. Both the $\mathbf{Y}_t^+$ term and the Jacobian term vanish, and the log-likelihood reduces to the $\mathbf{Y}_t^*$ expression above with $r$ replaced by $n$.
The Jacobian $\left| \mathbf{J}_t \right|$ depends on the choice of transformation $\mathbf{J}_t$.
In the GLS parameterization (Example 1i), the transformation depends on $\mathbf{W}_t$ and $\mathbf{U}_{\mathbf{Y},t}$. Under the normalization $\left| \mathbf{J}_t^{+} \mathbf{U}_{\mathbf{Y},t} \mathbf{J}_t^{+,\top} \right| = 1$, the Jacobian is given by
$
\left| \mathbf{J}_t \right| = |\mathbf{U}_{\mathbf{Y},t}|^{-1/2} |\mathbf{U}|^{1/2},
$ where $\mathbf{U} = (\mathbf{W}_t^\top \mathbf{U}_{\mathbf{Y},t}^{-1} \mathbf{W}_t)^{-1}$. In the orthogonal parameterization (Example 2i), $\mathbf{J}_t$ is orthonormal so that $\left| \mathbf{J}_t \right| = 1$.
\subsection{Smoothing}
The marginal smoothing distribution $p\left( \mathbf{A}_{t}|\mathbf{Y}_{1:T},\pmb{\theta}\right)$ provides information about the state matrix $\mathbf{A}_{t}$ conditional on all of the data $\mathbf{Y}_{1:T}$. It can be calculated recursively backwards through time after running the Kalman prediction recursions forwards and storing the matrices $\mathbf{V}_{t}, \mathbf{Z}^{\top}\mathbf{F}_{t}^{-1}, \mathbf{L}_{t}, \mathbf{A}_{t|t-1}$ and $\mathbf{P}_{t|t-1}$ for each date $t$.
\begin{proposition}
Consider the state space model (\ref{obsrvation equation})-(\ref{initial condition}). For $t=T,\ldots,1$, the marginal smoothing distribution is matrix normal $p\left( \mathbf{A}_{t}|\mathbf{Y}_{1:T},\pmb{\theta}\right) \ = \ \text{MN}\left( \mathbf{A}_{t|T},\mathbf{P}_{t|T},\mathbf{U}\right)$.
The mean matrix $\mathbf{A}_{t|T} $ and scale matrix $\mathbf{P}_{t|T}$ can be calculated recursively backwards
\begin{eqnarray}
\mathbf{G}_{t-1} \ = \ \mathbf{Z}_{t}^{\top}\mathbf{F}_{t}^{-1}\mathbf{V}_{t} + \mathbf{L}_{t}^{\top}\mathbf{G}_{t}, & & \mathbf{A}_{t|T} \ = \ \mathbf{A}_{t|t-1} + \mathbf{P}_{t|t-1}\mathbf{G}_{t-1}, \label{state smooth 1} \\
\mathbf{N}_{t-1} \ = \ \mathbf{Z}_{t}^{\top}\mathbf{F}_{t}^{-1}\mathbf{Z}_{t} + \mathbf{L}_{t}^{\top}\mathbf{N}_{t}\mathbf{L}_{t}, & & \mathbf{P}_{t|T} \ = \ \mathbf{P}_{t|t-1} - \mathbf{P}_{t|t-1}\mathbf{N}_{t-1}\mathbf{P}_{t|t-1}, \label{state smooth 2}
\end{eqnarray}
for $t=T,\ldots,1$ with initial values $\mathbf{G}_{T} = \pmb{0}_{s \times r}$ and $\mathbf{N}_{T} = \pmb{0}_{s\times s}$.
\end{proposition}
\vskip -0.5cm
There are several different forms of the Kalman smoother in the literature. This is a matrix version of the method developed by \cite{deJong(89)}.
\subsection{Drawing from the joint smoothing distribution} \label{simulation smoothing}
\subsubsection{Simulation smoothing}
Simulation smoothing algorithms are methods for drawing from the joint posterior distribution $p\left( \mathbf{A}_{1:T}|\mathbf{Y}_{1:T},\pmb{\theta}\right)$. Simulation smoothing algorithms for linear, Gaussian state space models were originally developed by \cite{CarterKohn(94)}, \cite{FruhwirthSchnatter(94)}, \cite{deJongShephard(95)}, and \cite{DurbinKoopman(02)}.
The following algorithm is a matrix generalization of \cite{DurbinKoopman(02)}.
\begin{proposition}
Consider the state space model (\ref{obsrvation equation})-(\ref{initial condition}). A random draw $\mathbf{A}_{1}^{d}, \ldots,\mathbf{A}_{T}^{d}$ from $p\left( \mathbf{A}_{1:T}|\mathbf{Y}_{1:T},\pmb{\theta}\right)$
can be obtained using the following algorithm
\begin{enumerate}
\item For $t=1,\ldots,T$, simulate a new state matrix $\mathbf{A}_{t}^{\dagger}$ and data series $\mathbf{Y}_{t}^{*\dagger}$ from the model (\ref{transition equation})-(\ref{initial condition}) and (\ref{obs equation 1}). When running this step, all mean terms are set to zero, i.e. $\mathbf{D}_{t}^{*} = \pmb{0}$ and $\mathbf{C}_{t} = \pmb{0}$, including the initial state $\mathbf{A}_{1}^{\dagger} \sim \text{MN}\left( \pmb{0},\mathbf{P}_{1|0},\mathbf{U}\right)$.
\item Construct the artificial data series $\widehat{\mathbf{Y}}_{t}^{*} = \mathbf{Y}_{t}^{*} - \mathbf{Y}_{t}^{*\dagger}$. Run the matrix Kalman filter and smoother on the data series $\widehat{\mathbf{Y}}_{t}^{*}$, storing the smoothed estimates $\widehat{\mathbf{A}}_{t|T}$.
\item For $t=1,\ldots,T$, calculate $\mathbf{A}_{t}^{d} = \widehat{\mathbf{A}}_{t|T} + \mathbf{A}_{t}^{\dagger}$.
\end{enumerate}
\end{proposition}
In step 2, only equations (\ref{state smooth 1}) need calculated during the backwards pass of the Kalman smoother. This algorithm works well for most models and has the benefit that minimal changes are needed when going from one model to another.
\subsubsection{Precision sampler}
\cite{ChanJeliazkov(09)} proposed an algorithm for drawing from the joint smoothing distribution in a linear, Gaussian state space model that does not use the Kalman filter. It can be faster for some models because it takes the Cholesky decomposition of the posterior precision matrix which is a lower-triangular, banded matrix instead of the covariance matrix.
To implement the precision sampler, we need to re-define the state space model for the transformed data $\mathbf{Y}_{t}^{*} = \mathbf{Y}_{t}\mathbf{J}_{t}^{*,\top}$ without using the companion form. In this subsection, we use tilde's above the matrices to differentiate the notation. Let $\widetilde{\mathbf{A}}_{t}$ be a $\widetilde{s} \times r$ latent state with autoregressive dynamics of order $p$ given by
\begin{eqnarray*}
\mathbf{Y}_{t}^{*} &=& \widetilde{\mathbf{D}}_{t}^{*} + \widetilde{\mathbf{Z}}_{t}\widetilde{\mathbf{A}}_{t} + \widetilde{\mathbf{E}}_{1t}, \qquad \widetilde{\mathbf{E}}_{1t} \sim \text{MN}\left( \pmb{0},\mathbf{H}_{t},\mathbf{U}\right), \\
\widetilde{\mathbf{A}}_{t} &=& \widetilde{\mathbf{C}}_{t} + \widetilde{\mathbf{T}}_{1}\widetilde{\mathbf{A}}_{t-1} + \ldots + \widetilde{\mathbf{T}}_{p}\widetilde{\mathbf{A}}_{t-p} + \widetilde{\mathbf{E}}_{2t}, \qquad \widetilde{\mathbf{E}}_{2t} \sim \text{MN}\left( \pmb{0},\widetilde{\mathbf{Q}}_{t},\mathbf{U}\right),
\end{eqnarray*}
with initial condition $\widetilde{\mathbf{A}}_{1} \sim \text{MN}\left( \widetilde{\mathbf{A}}_{1|0},\widetilde{\mathbf{P}}_{1|0},\mathbf{U}\right)$. The $r \times r$ matrix $\mathbf{U}$ and the $m \times m$ matrix $\mathbf{H}_{t}$ are the same as before. The dimension of all other system matrices are adjusted. The transition matrices $\widetilde{\mathbf{T}}_{j}$ are also assumed to be constant over time.
Next, we stack the observation equations together
and define the following matrices
\[ \underset{Tm \times r}{\widetilde{\mathbf{Y}}^{*}} \ = \ \left[ \begin{matrix}
\mathbf{Y}_{1}^{*} \\
\vdots \\
\mathbf{Y}_{T}^{*}
\end{matrix}\right] \qquad \underset{Tm \times Tr}{\widetilde{\mathbf{Z}}} \ = \ \left[ \begin{matrix}
\widetilde{\mathbf{Z}}_{1} & \ldots & \pmb{0} \\
\vdots & \ddots & \vdots \\
\pmb{0} & \ldots &\widetilde{\mathbf{Z}}_{T}
\end{matrix}\right] \qquad \underset{Tm \times Tm}{\widetilde{\mathbf{H}}} \ = \ \left[ \begin{matrix}
\mathbf{H}_{1} & \ldots & \pmb{0} \\
\vdots & \ddots & \vdots \\
\pmb{0} & \ldots & \mathbf{H}_{T} \\
\end{matrix}\right] \]
with $\widetilde{\mathbf{D}}^{*} = \left( \widetilde{\mathbf{D}}_{1}^{*\top},\widetilde{\mathbf{D}}_{2}^{*\top},\ldots,\widetilde{\mathbf{D}}_{T}^{*\top}\right)^{\top}$.
We also stack the state matrices into a $T \tilde{s} \times r$ matrix $\widetilde{\mathbf{A}} = \left( \widetilde{\mathbf{A}}_{1}^{\top},\widetilde{\mathbf{A}}_{2}^{\top},\ldots,\widetilde{\mathbf{A}}_{T}^{\top}\right)^{\top}$ and define
\[ \underset{T\tilde{s} \times T\tilde{s}}{\widetilde{\mathbf{T}}} \ =\ \left[ \begin{matrix}
\mathbf{I} & \pmb{0} & \ldots & \ldots & \ldots & \pmb{0} \\
-\widetilde{\mathbf{T}}_{1} & \mathbf{I} & \pmb{0} & \ldots & \ldots & \vdots \\
-\widetilde{\mathbf{T}}_{2} & -\widetilde{\mathbf{T}}_{1} & \mathbf{I} & \ldots & \ldots & \vdots \\
\vdots & \vdots & \vdots & & \ddots & \pmb{0} \\
\pmb{0} & \ldots & -\widetilde{\mathbf{T}}_{p} & \ldots & -\widetilde{\mathbf{T}}_{1} & \mathbf{I} \\
\end{matrix}\right] \ \
\underset{T\tilde{s} \times T\tilde{s}}{\widetilde{\mathbf{Q}}} \ = \ \left[ \begin{matrix}
\widetilde{\mathbf{P}}_{1|0} & \pmb{0} & \ldots & \pmb{0} \\
\pmb{0} & \widetilde{\mathbf{Q}}_{2} & \ldots & \pmb{0} \\
\vdots & \vdots & \ddots & \vdots \\
\pmb{0} & \pmb{0} & \ldots & \widetilde{\mathbf{Q}}_{T} \\
\end{matrix}\right] \]
with $\widetilde{\mathbf{C}} = \left( \widetilde{\mathbf{A}}_{1|0}^{\top},\widetilde{\mathbf{C}}_{2}^{\top},\ldots,\widetilde{\mathbf{C}}_{T}^{\top}\right)^{\top}$. Moving average terms can be incorporated by making $\widetilde{\mathbf{Q}}$ a banded matrix, though the bandwidth grows with the MA order and erodes the sparsity advantage over the simulation smoother.
Using this notation, the model can be written in stacked form as
\begin{eqnarray*}
\widetilde{\mathbf{Y}}^{*} &=& \widetilde{\mathbf{D}}^{*} + \widetilde{\mathbf{Z}}\widetilde{\mathbf{A}} + \widetilde{\mathbf{E}}_{1}^{*}, \qquad \widetilde{\mathbf{E}}_{1}^{*} \sim \text{MN}\left( \pmb{0},\widetilde{\mathbf{H}},\mathbf{U}\right), \\
\widetilde{\mathbf{A}} &=& \widetilde{\mathbf{T}}^{-1}\widetilde{\mathbf{C}} + \widetilde{\mathbf{E}}_{2}, \qquad \widetilde{\mathbf{E}}_{2} \sim \text{MN}\left( \pmb{0},\widetilde{\mathbf{T}}^{-1}\widetilde{\mathbf{Q}}\widetilde{\mathbf{T}}^{-1,\top},\mathbf{U}\right).
\end{eqnarray*}
We can interpret the first equation as the likelihood and the second equation as the prior.
Using Bayes rule, the posterior distribution for the stacked state matrix is $\widetilde{\mathbf{A}} \sim \text{MN}\left( \widetilde{\mathbf{M}},\widetilde{\pmb{\Sigma}},\mathbf{U}\right)$ with $T \tilde{s} \times T\tilde{s}$ precision matrix $\widetilde{\pmb{\Sigma}}^{-1} = \widetilde{\mathbf{T}}^{\top}\widetilde{\mathbf{Q}}^{-1}\widetilde{\mathbf{T}} + \widetilde{\mathbf{Z}}^{\top}\widetilde{\mathbf{H}}^{-1}\widetilde{\mathbf{Z}}$.
and $T \tilde{s} \times r$ posterior mean matrix $\widetilde{\mathbf{M}} = \widetilde{\pmb{\Sigma}}\left( \widetilde{\mathbf{T}}^{\top}\widetilde{\mathbf{Q}}^{-1}\widetilde{\mathbf{C}} + \widetilde{\mathbf{Z}}^{\top}\widetilde{\mathbf{H}}^{-1}\left[\widetilde{\mathbf{Y}}^{*} - \widetilde{\mathbf{D}}^{*}\right]\right)$.
To draw from this distribution efficiently, we take the following steps
\begin{enumerate}
\item Calculate the matrix $\widetilde{\pmb{\Sigma}}^{-1}$ and its Cholesky decomposition $\widetilde{\pmb{\Sigma}}^{-1} = \mathbf{L}_{\sigma}\mathbf{L}_{\sigma}^{\top}$.
\item Calculate the Cholesky decomposition of $\mathbf{U} = \mathbf{L}_{\mathbf{U}}\mathbf{L}_{\mathbf{U}}^{\top}$.
\item Calculate the matrix $\widetilde{\mathbf{G}}= \mathbf{L}_{\sigma} \backslash \left( \widetilde{\mathbf{T}}^{\top}\widetilde{\mathbf{Q}}^{-1}\widetilde{\mathbf{C}} + \widetilde{\mathbf{Z}}^{\top}\widetilde{\mathbf{H}}^{-1}\left[\widetilde{\mathbf{Y}}^{*} - \widetilde{\mathbf{D}}^{*}\right]\right) / \mathbf{L}_{\mathbf{U}}^{\top}$ using forward substitution with $\mathbf{L}_{\sigma}$ and backward substitution with $\mathbf{L}_{\mathbf{U}}^{\top}$. Here, $\mathbf{A} \backslash \mathbf{B}$ and $\mathbf{B} / \mathbf{A}$ denote the solutions $\mathbf{X}$ of $\mathbf{A}\mathbf{X} = \mathbf{B}$ and $\mathbf{X}\mathbf{A} = \mathbf{B}$, respectively, computed by triangular solves rather than by forming $\mathbf{A}^{-1}$ explicitly.
\item A random $T\widetilde{s} \times r$ matrix $\mathbf{A}^{d}$ drawn from the joint distribution is then obtained by drawing $\widetilde{\mathbf{E}} \sim \text{MN}\left( \pmb{0},\mathbf{I}_{T\widetilde{s}},\mathbf{I}_{r}\right)$ and calculating $\mathbf{A}^{d} = \left( \mathbf{L}_{\sigma}^{\top}\backslash \left[ \widetilde{\mathbf{G}} + \widetilde{\mathbf{E}}\right] \right) \mathbf{L}_{\mathbf{U}}^{\top}.$
\end{enumerate}
In our experience, the precision sampler is typically faster than the simulation smoother for autogressive models and simple time series models. However, the simulation smoother is better for estimating unknown constant parameters in the matrices $\mathbf{C}_{t}$ and $\mathbf{D}_{t}$.
\subsection{Estimation of regression parameters}
Let $\mathbf{X}_t$ denote a matrix of observed covariates. Covariates can be incorporated by specifying the intercept matrices $\mathbf{D}_t$ and/or $\mathbf{C}_t$ in (\ref{obsrvation equation})--(\ref{transition equation}) as functions of $\mathbf{X}_t$. For example, one may set $\mathbf{C}_t = \mathbf{X}_t \pmb{\beta}$, where $\mathbf{X}_t$ is $s \times \ell$ and $\pmb{\beta}$ is an $\ell \times r$ matrix of unknown parameters. For Bayesian analysis, a conjugate prior can be specified as $\pmb{\beta} \sim \text{MN}(\underline{\mathbf{M}}, \underline{\mathbf{P}}, \mathbf{U})$ and $\pmb{\beta}$ can be drawn jointly with the latent states using the simulation smoother. For maximum likelihood estimation, a profile likelihood or concentrated estimator for $\pmb{\beta}$ can be derived analogously to \cite{deJong(91)}.
\subsection{Missing values} \label{Missing values}
In empirical work, time series data may have missing values. The state space model (\ref{obsrvation equation})-(\ref{initial condition}) can accommodate missing values as long as the entire row or column vector of $\mathbf{Y}_{t}$ is treated as missing. While potentially restrictive, this is still empirically relevant for researchers working with mixed frequency data (quarterly, monthly, etc.) whose release dates are the same for all spatial units; e.g. government agencies often release economic data with this pattern. Mixed frequency data is popular in economics; see, e.g. \cite{MarianoMurasawa(03)}, and \cite{SchorfheideSong(15)}.
To handle missing values, we define an $m_{t} \times m$ matrix $\mathbf{S}_{t}$ that selects out rows of the observed data at each date.\footnote{An entirely missing column can be handled symmetrically, by right-multiplying $\mathbf{Y}_t$ by an analogous $n_t\times n$ column-selection matrix $\mathbf{S}_t^{c}$; this transforms $\mathbf{W}_t\to\mathbf{S}_t^{c}\mathbf{W}_t$ and the column scale $\mathbf{U}_{\mathbf{Y},t}\to\mathbf{S}_t^{c}\mathbf{U}_{\mathbf{Y},t}\mathbf{S}_t^{c,\top}$, leaving $\mathbf{H}_t$ unchanged.} We then left-multiply the original observation equation (\ref{obsrvation equation}) by $\mathbf{S}_{t}$ to get $\mathbf{Y}_{t}^{-} \ = \ \mathbf{D}_{t}^{-} + \mathbf{Z}_{t}^{-}\mathbf{A}_{t}\mathbf{W}_{t}^{\top} + \mathbf{E}_{1t}^{-}$, with $\mathbf{E}_{1t}^{-} \sim \text{MN}\left( \pmb{0},\mathbf{H}_{t}^{-},\mathbf{U}_{\mathbf{Y}}\right)$.
The adjusted matrices are defined as $\mathbf{Y}_{t}^{-} = \mathbf{S}_{t}\mathbf{Y}_{t}$, $\mathbf{D}_{t}^{-} = \mathbf{S}_{t}\mathbf{D}_{t}$, $\mathbf{Z}_{t}^{-} = \mathbf{S}_{t}\mathbf{Z}_{t}$, and $\mathbf{H}_{t}^{-} \ = \ \mathbf{S}_{t}\mathbf{H}_{t}\mathbf{S}_{t}^{\top}$.
The Kalman filter and associated algorithms can then be applied as usual but with system matrices whose dimension are changing over time.
Entry-level (non-row-wise) missingness can also be handled by data augmentation, drawing the missing entries from their conditional Gaussian distribution within the Gibbs sampler and treating them as observed thereafter. This is the approach we use for the ragged-edge missing monthly observations in the empirical application (Section \ref{Gibbs sampler}). The frequentist analogue is the EM algorithm, which replaces each missing entry with its conditional expectation in the E-step. Both accommodate more general missingness patterns than row-selection.
\subsection{Efficient sampling of discrete indicators in Markov switching and mixture models} \label{Sampling discrete states}
To add flexibility to the model (\ref{obsrvation equation})-(\ref{initial condition}), researchers often allow the parameters within the system matrices to be time-varying. A common approach is to make the parameters a function of a finite mixture or Markov-switching variable; see, e.g. \cite{Kim(94)}, \cite{KimNelsonBook(99)} and \cite{GiordaniKohnvanDijk(07)}. Let $s_{t}$ denote an indicator variable that can take on a finite number of values $s_{t}=j$ for $j = 1,\ldots,J$. Each discrete state $s_{t}=j$ corresponds to a different set of parameters. The system matrices are functions of the discrete state
\begin{eqnarray}
\mathbf{Y}_{t} & = & \mathbf{D}\left( s_{t}\right) + \mathbf{Z}\left( s_{t}\right)\mathbf{A}_{t}\mathbf{W}_{t}^{\top} + \mathbf{E}_{1t}, \qquad \mathbf{E}_{1t} \sim \text{MN}\left( \pmb{0},\mathbf{H}\left( s_{t}\right),\mathbf{U}_{\mathbf{Y},t}\right), \label{obsrvation equation MS} \\
\mathbf{A}_{t+1} &=& \mathbf{C}\left( s_{t+1}\right) + \mathbf{T}\left( s_{t+1}\right)\mathbf{A}_{t} + \mathbf{R}\left( s_{t+1}\right)\mathbf{E}_{2t}, \qquad \mathbf{E}_{2t} \sim \text{MN}\left( 0,\mathbf{Q}\left( s_{t+1}\right),\mathbf{U}\right), \label{transition equation MS}
\end{eqnarray}
where $\mathbf{A}_{1} \sim \text{MN}\left( \mathbf{A}_{1|0}(s_{1}),\mathbf{P}_{1|0}\left( s_{1}\right),\mathbf{U}\right)$. One can generalize this further to make the system matrices a function of a vector of discrete states.\footnote{This extension is straightforward but notationally more complex. The algorithm for drawing the discrete states does not fundamentally change.}
\cite{GerlachCarterKohn(00)} and \cite{DoucetAndrieu(01)} develop algorithms for drawing from the conditional distribution $p \left( s_{t}|\mathbf{s}_{-t},\mathbf{Y}_{1:T},\pmb{\theta}\right) \propto p\left( \mathbf{Y}_{1:T}|\mathbf{s}_{-t},s_{t},\pmb{\theta}\right) p\left( s_{t}|\mathbf{s}_{-t},\pmb{\theta}\right)$ of a single discrete state but with the continuous states integrated out. The notation $\mathbf{s}_{-t}$ is a vector of all indicator variables but with $s_{t}$ omitted. For Bayesian estimation with MCMC, this improves mixing of the Markov chain because it is a collapsed Gibbs sampler. We extend their approach to the matrix state space model where the key result is the following proposition.
\begin{proposition} \label{proposition discrete state}
Consider the state space model (\ref{obsrvation equation MS})-(\ref{transition equation MS}). The likelihood of $\mathbf{Y}_{1:T}$ conditional on the indicators $s_{1:T}$ is given up to proportionality by
\begin{eqnarray*}
p\left( \mathbf{Y}_{1:T}|s_{t},\mathbf{s}_{-t},\pmb{\theta}\right)
&\propto & \left| \mathbf{F}_{t}\left( s_{1:t}\right)\right|^{-\frac{r}{2}} \left| \mathbf{P}_{t|t}\left( s_{1:t}\right)\right|^{-\frac{r}{2}} \left| \mathbf{P}_{t|T}\left( s_{1:T}\right)\right|^{\frac{r}{2}}\left|\mathbf{H}\left( s_{t}\right)\right|^{-\frac{(n-r)}{2}}\\
& & \exp\left( -\frac{1}{2}\text{tr}\left[ \mathbf{U}^{-1}\mathbf{V}_{t}\left( s_{1:t}\right)^{\top}\mathbf{F}_{t}\left( s_{1:t}\right)^{-1}\mathbf{V}_{t}\left( s_{1:t}\right)\right] \right) \\
& & \exp\left( -\frac{1}{2}\text{tr}\left[
\mathbf{H}\left( s_{t}\right)^{-1} \left[\mathbf{Y}_t^+ - \mathbf{D}\left( s_{t}\right)^{+}\right]\pmb{\Psi}_{t}^{-1}\left[\mathbf{Y}_t^{+}-\mathbf{D}\left( s_{t}\right)^{+}\right]^{\top}
\right]\right)
\\
& & \exp\left(-\frac{1}{2}\text{tr}\left[ \mathbf{U}^{-1}\mathbf{A}_{t|t}\left( s_{1:t}\right)^{\top}\mathbf{P}_{t|t}\left( s_{1:t}\right)^{-1}\mathbf{A}_{t|t}\left( s_{1:t}\right)\right]\right) \\
& & \exp\left( \frac{1}{2}\text{tr}\left[ \mathbf{U}^{-1}\mathbf{A}_{t|T}\left( s_{1:T}\right)^{\top}\mathbf{P}_{t|T}\left( s_{1:T}\right)^{-1}\mathbf{A}_{t|T}\left( s_{1:T}\right)\right]\right)
\end{eqnarray*}
where the constant of proportionality does not depend on $s_{t}$.
\end{proposition}
The key insight of \cite{GerlachCarterKohn(00)} and \cite{DoucetAndrieu(01)} is an algorithm for calculating $p\left( \mathbf{Y}_{1:T}|\mathbf{s}_{-t},s_{t}=j,\pmb{\theta}\right)$ for each state $s_{t} = j$ in a computationally efficient way. First, one runs a backwards pass of the Kalman information filter conditional on a previous MCMC draw of the discrete states $s_{1:T}$. Then, the algorithm iterates forward in time drawing the discrete states $s_{t}$ for $t=1,\ldots,T$ while evaluating $p\left( \mathbf{Y}_{1:T}|s_{t}=j,\mathbf{s}_{-t},\pmb{\theta}\right)$ for each state $s_{t}=j$ for $j=1,\ldots,J$.
The algorithm for drawing the discrete states $s_{t}$ from their conditional distribution $p \left( s_{t}|\mathbf{s}_{-t},\mathbf{Y}_{1:T},\pmb{\theta}\right)$ sequentially forward for $t=1,\ldots,T$ is
\begin{description}
\item[1.] Given the current set of indicators $s_{1:T}$ from a previous MCMC iteration, run the backwards information filter. Initialize $\mathbf{B}_{T|T} = \pmb{0}$ and $\pmb{\Pi}_{T|T} = \pmb{0}$. And, for $t=T-1
\ldots,1$, calculate
\begin{eqnarray*}
\underline{\pmb{\Pi}}_{t+1|T} &=& \pmb{\Pi}_{t+1|T} + \mathbf{Z}\left( s_{t+1}\right) ^{\top}\mathbf{H}\left( s_{t+1}\right)^{-1}\mathbf{Z}\left( s_{t+1}\right) \\
\underline{\mathbf{B}}_{t+1|T}& = & \mathbf{B}_{t+1|T} + \mathbf{Z}\left( s_{t+1}\right) ^{\top}\mathbf{H}\left( s_{t+1}\right) ^{-1}\left(\mathbf{Y}_{t+1} - \mathbf{D}\left( s_{t+1}\right) \right) \mathbf{J}_{t+1}^{*,\top} \\
\pmb{\Delta}_{t+1} &=& \mathbf{I}_{s} + \underline{\pmb{\Pi}}_{t+1|T}\mathbf{R}\left( s_{t+1}\right) \mathbf{Q}\left( s_{t+1}\right) \mathbf{R}\left( s_{t+1}\right) ^{\top} \\
\pmb{\Pi}_{t|T}\left( \mathbf{s}_{t+1:T}\right) &=& \mathbf{T}\left( s_{t+1}\right) ^{\top}\pmb{\Delta}_{t+1}^{-1}\underline{\pmb{\Pi}}_{t+1|T}\mathbf{T}\left( s_{t+1}\right) \\
\mathbf{B}_{t|T}\left( \mathbf{s}_{t+1:T}\right) &=& \mathbf{T}\left( s_{t+1}\right) ^{\top}\pmb{\Delta}_{t+1}^{-1}\left(\underline{\mathbf{B}}_{t+1|T}-\underline{\pmb{\Pi}}_{t+1|T}\mathbf{C}\left( s_{t+1}\right) \right)
\end{eqnarray*}
During the backwards pass, store $\mathbf{B}_{t|T}\left( \mathbf{s}_{t+1:T}\right) $ and $\pmb{\Pi}_{t|T}\left( \mathbf{s}_{t+1:T}\right)$.
\vskip 0.1cm
\item[2.] For $t=1,\ldots,T$ compute for each state $s_{t}=j$ for $j=1,\ldots,J$ at each date $t$
\begin{description}
\item[(2a.)] For $t=1$, the initial condition is $\mathbf{A}_{1|0}\left( s_{1}=j\right)$ and $\mathbf{P}_{1|0}\left( s_{1}=j\right)$. For $t=2,\ldots,T$, assume that we have already sampled $\mathbf{s}_{1:t-1}$ and calculated $\mathbf{A}_{t-1|t-1}\left(\mathbf{s}_{1:t-1}\right)$ and $\mathbf{P}_{t-1|t-1}\left( \mathbf{s}_{1:t-1}\right)$ at a previous iteration. Then, we calculate the one-step ahead predictive distribution for each $s_{t} = j$
\begin{eqnarray*}
\mathbf{A}_{t|t-1}\left( \mathbf{s}_{1:t}\right) &=& \mathbf{T}\left( s_{t}\right)\mathbf{A}_{t-1|t-1}\left( \mathbf{s}_{1:t-1}\right) + \mathbf{C}\left( s_{t}\right) \\
\mathbf{P}_{t|t-1}\left( \mathbf{s}_{1:t}\right) &=& \mathbf{T}\left( s_{t}\right)\mathbf{P}_{t-1|t-1}\left( \mathbf{s}_{1:t-1}\right)\mathbf{T}\left( s_{t}\right)^{\top} + \mathbf{R}\left( s_{t}\right)\mathbf{Q}\left( s_{t}\right)\mathbf{R}\left( s_{t}\right)^{\top}
\end{eqnarray*}
\item[(2b.)] Calculate the prediction error and prediction error variance
\begin{eqnarray*}
\mathbf{V}_{t}\left( \mathbf{s}_{1:t}\right) &=& \mathbf{Y}_{t}^{*}-\mathbf{D}_{t}\left( \mathbf{s}_{t}\right) \mathbf{J}_{t}^{*,\top} - \mathbf{Z}\left( s_{t}\right)\mathbf{A}_{t|t-1}\left( s_{1:t}\right)\\
\mathbf{F}_{t}\left( \mathbf{s}_{1:t}\right) &=& \mathbf{Z}\left( s_{t}\right)\mathbf{P}_{t|t-1}\left( s_{1:t}\right) \mathbf{Z}\left( s_{t}\right)^{\top} + \mathbf{H}\left( s_{t}\right)
\end{eqnarray*}
and the filter
\begin{eqnarray*}
\mathbf{K}_{t}\left( \mathbf{s}_{1:t}\right) &=& \mathbf{P}_{t|t-1}\left( \mathbf{s}_{1:t}\right)\mathbf{Z}\left( s_{t}\right)^{\top}\mathbf{F}_{t}\left( \mathbf{s}_{1:t}\right)^{-1} \\
\mathbf{A}_{t|t}\left( \mathbf{s}_{1:t}\right) & = & \mathbf{A}_{t|t-1}\left( \mathbf{s}_{1:t}\right) + \mathbf{K}_{t}\left( \mathbf{s}_{1:t}\right)\mathbf{V}_{t}\left( \mathbf{s}_{1:t}\right) \\
\mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right) & = & \mathbf{P}_{t|t-1}\left( \mathbf{s}_{1:t}\right) - \mathbf{K}_{t}\left( \mathbf{s}_{1:t}\right) \mathbf{Z}\left( s_{t}\right)\mathbf{P}_{t|t-1}\left( \mathbf{s}_{1:t}\right)
\end{eqnarray*}
\item[(2c.)] Combine the output of the backward information filter and the forward Kalman filter to calculate the smoothed estimate for each $s_{t} = j$
\begin{eqnarray*}
\pmb{\Gamma}_{t}\left( \mathbf{s}_{1:T}\right) &=& \mathbf{I}_{s} + \pmb{\Pi}_{t|T}\left( s_{t+1|T}\right)\mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right) \\
\mathbf{A}_{t|T}\left( \mathbf{s}_{1:T}\right) &=& \mathbf{A}_{t|t}\left( \mathbf{s}_{1:t}\right) + \mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right) \pmb{\Gamma}_{t}\left( \mathbf{s}_{1:T}\right)^{-1}\left[ \mathbf{B}_{t|T}\left( \mathbf{s}_{t+1:T}\right) \right. \\
& & \left. - \pmb{\Pi}_{t|T}\left( \mathbf{s}_{t+1:T}\right)\mathbf{A}_{t|t}\left( \mathbf{s}_{1:t}\right) \right] \\
\mathbf{P}_{t|T}\left( \mathbf{s}_{1:T}\right) &=& \mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right) - \mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right) \pmb{\Gamma}_{t}\left( \mathbf{s}_{1:T}\right)^{-1}\pmb{\Pi}_{t|T}\left( s_{t+1|T}\right)\mathbf{P}_{t|t}\left( \mathbf{s}_{1:t}\right)
\end{eqnarray*}
\item[(2d.)] Draw $s_{t}$ from the discrete distribution
\[ p\left( s_{t}=j|\mathbf{Y}_{1:T},\mathbf{s}_{-t},\pmb{\theta}\right) \propto p\left( \mathbf{Y}_{1:T}|s_{t}=j,\mathbf{s}_{-t},\pmb{\theta}\right) p\left( s_{t}=j|\mathbf{s}_{-t},\pmb{\theta}\right).\]
where $p\left( \mathbf{Y}_{1:T}|s_{t}=j,\mathbf{s}_{-t},\pmb{\theta}\right)$ is given in Proposition \ref{proposition discrete state}. Save $\mathbf{A}_{t|t}\left( s_{t}=j\right)$ and $\mathbf{P}_{t|t}\left( s_{t}=j\right)$ for the next iteration and return to step \textbf{(2a.)} if $t<T$.
\end{description}
\end{description}
While proving Proposition \ref{proposition discrete state}, we derived a second type of smoothing algorithm for the matrix state space model known as two-filter formula smoothing. This algorithm runs the matrix Kalman filter forward in time and the matrix backwards information filter recursively backwards, combining the two filters to calculate the smoothed distribution, see e.g \cite{Mayne(66)}.
\section{Application: state level macroeconomic data}
\subsection{Data}
In our empirical application, we estimate the matrix bilinear dynamic factor model given by (\ref{Example 2a})--(\ref{Example 2c}) on a mixed-frequency data set of $m = 12$ economic variables for $n = 50$ U.S. states. Using the FRED SD database (see \cite{BokunJacksonKliesenOwyang(24)}), we extract monthly series on state-level unemployment and employment across sectors, along with three quarterly series: nominal personal income, real gross domestic product, and the FHFA home price index. The sample spans February 1990 through August 2025, yielding $T = 427$ observations. Additional details on data construction and transformations are provided in the online appendix.
\subsection{Model} \label{Model application}
We allow for common factors across both variables and U.S. states through the loading matrices $\pmb{\Lambda} \in \mathbb{R}^{m \times r_{\ell}}$ and $\mathbf{W} \in \mathbb{R}^{n \times r_{r}}$. We adopt the orthogonal parameterization described in Example 2i, imposing $\pmb{\Lambda}^{\top}\pmb{\Lambda} = \mathbf{I}$ and $\mathbf{W}^{\top}\mathbf{W} = \mathbf{I}$ with $\mathbf{J}_{t}^{*} = \mathbf{W}$ and $\mathbf{J}_{t}^{+} = \mathbf{W}_{\perp}$. To address rotational indeterminacy, the covariance matrices $\mathbf{U}$ and $\pmb{\Omega}$ are restricted to be diagonal with eigenvalues ordered in descending order.
The observation right-scale matrix $\mathbf{U}_{\mathbf{Y},t}$ is defined in (\ref{U param}).
We allow $\pmb{\Psi}$ to be dense and estimate it directly.
We assume a diagonal baseline structure for the left scale matrix $\pmb{\Sigma}$
to reduce dimensionality. To accommodate heavy-tailed measurement errors, we multiply these row variances by series-specific Student's $t$ scale mixtures, so that the observation left scale is
$\pmb{\Sigma}_t=\mathrm{diag}(h_{1t}\sigma_1^2,\ldots,h_{mt}\sigma_m^2)$.
Here, $h_{it} \sim \text{I.G.}\left( \nu_i/2,\nu_i/2\right)$ are latent scale inverse gamma variables that represent the Student's $t$ distribution as a normal mixture.
We allow for additional time variation in the state equation by specifying the left scale matrix as $\pmb{\Omega}_t=\kappa_t\pmb{\Omega}$, where $\kappa_t$ is a scalar two-state Markov switching process. One regime sets $\kappa_t=1$ and the other sets $\kappa_t=5$, allowing the model to capture common extreme shocks such as those observed during COVID. To identify scale, we impose $\mathrm{tr}(\pmb{\Omega})/s = 1$ for the transition covariance. The autoregressive matrices $\pmb{\Phi}_1$ and $\pmb{\Phi}_2$ are restricted to lie in the stationary region.
We denote by $\mathbf{Y}_{1:T}^{o}$ the observed monthly and quarterly data. Following \cite{MarianoMurasawa(03)}, the model is specified at a monthly frequency. Quarterly observations are treated as averages of the latent monthly variables, and missing quarterly values $\mathbf{Y}_{1:T}^{q}$ are imputed within the MCMC algorithm.
\subsection{Bayesian estimation}
\subsubsection{Priors} \label{Priors}
We adopt weakly informative priors that respect the structural constraints of the model and facilitate efficient posterior simulation using a Gibbs sampler with parameter expanded data augmentation (PXDA), \cite{LiuWu(99)} and \cite{MengVanDyk(99)}. Full details of the prior specification and hyperparameter choices are provided in the online appendix.
The factor loadings $\pmb{\Lambda}$ and $\mathbf{W}$ are restricted to lie on the Stiefel manifold. We therefore place matrix von Mises--Fisher priors on each, following \cite{Hoff(09)}, which enforce orthonormality while allowing for flexible prior centering. We place a matrix normal prior on the autoregressive parameters
$\pmb{\Phi} = (\pmb{\Phi}_{1}, \pmb{\Phi}_{2}) \sim \mathrm{MN}(\underline{\pmb{\Phi}}, \underline{\mathbf{V}}, \pmb{\Omega}),
$ which centers the dynamics on a stationary VAR(2) process with diagonal shrinkage across lags and factors. The covariance matrices $\pmb{\Omega}$, $\mathbf{U}$, and $\pmb{\Psi}$ are assigned inverse Wishart priors in their parameter-expanded forms, which lead to conditionally conjugate updates within the PXDA framework.
The observation variances $\sigma_{i}^{2}$ are assigned inverse gamma priors, $
\sigma_{i}^{2} \sim \mathrm{IG}(\underline{a}_i, \underline{b}_i),
$ while the initial state follows
$
\mathbf{A}_1 \sim \mathrm{MN}(\pmb{0}, \underline{\mathbf{P}}_1, \mathbf{U}).
$
These priors are scaled using empirical moments of the data to ensure comparable magnitudes across series and states.
\subsubsection{Gibbs sampler} \label{Gibbs sampler}
The model in (\ref{Example 2a})--(\ref{Example 2c}) together with the priors leads to a partially collapsed Gibbs sampler that combines conjugate updates, PXDA, and Metropolis--Hastings steps on the Stiefel manifold. At a high level, each iteration proceeds as follows:
\begin{description}
\item[(1)] Draw the two-state Markov switching indicators $s_{t}$ for $t=1,\ldots,T$ that determine the transition-scale multiplier $\kappa_t$, with the latent states integrated out.
\item[(2)] Draw $(\pmb{\Psi},\mathbf{U},\mathbf{A}_{1:T},\mathbf{Y}_{1:T}^{q})$ jointly using a partially collapsed step. The covariance matrices $\pmb{\Psi}$ and $\mathbf{U}$ are drawn using PXDA with the latent states integrated out via the Kalman filter. Conditional on these draws, the latent states are sampled using a simulation smoother that skips the systematically missing quarterly data. Finally, the missing quarterly observations are drawn from their conditional distributions.
\item[(3)] Draw ragged-edge missing monthly observations from their conditional Gaussian distributions.
\item[(4)] Draw the left factor loadings $\pmb{\Lambda}$ using a Metropolis--Hastings algorithm on the Stiefel manifold, targeting the posterior implied by its matrix von Mises--Fisher prior.
\item[(5)] Draw the right factor loadings $\mathbf{W}$ using a Metropolis--Hastings algorithm on the Stiefel manifold that accounts for both the factor space and its orthogonal complement. During this step, the parameters in $\pmb{\Psi}$ are integrated out.
\item[(6)] Draw the state innovation covariance matrix $\pmb{\Omega}$ using a PXDA step, followed by a rotation and reparameterization of the state equation.
\item[(7)] Draw the autoregressive parameters $\pmb{\Phi} = (\pmb{\Phi}_1,\pmb{\Phi}_2)$ from a matrix normal distribution, imposing stationarity via rejection sampling.
\item[(8)] Draw the observation variances $\sigma^{2}_{i}$ and the Student's $t$ scale mixtures $h_{it}$ used for the measurement errors.
\end{description}
Several steps of the algorithm are worth highlighting. First, several blocks integrate out the latent states using the Kalman filter, leading to a partially collapsed sampler that improves mixing. Second, the factor loadings are constrained to lie on the Stiefel manifold and are therefore sampled using Metropolis--Hastings steps rather than direct draws from their full conditional distributions. The combination of PXDA and collapsed steps substantially improves mixing in high-dimensional settings. Finally, the discrete states are drawn in step (1) using the algorithm of Section 3.7. Full details of all algorithmic steps are provided in the online appendix.
\subsection{Empirical results}
\begin{figure}[!t]
\caption{Data and estimated conditional mean.} \label{fig:data_fit}
\begin{center}
\resizebox{0.80\textwidth}{!}{
\includegraphics{unemployment_rate.pdf}
\includegraphics{employment_data_CONS.pdf}}
\resizebox{0.80\textwidth}{!}{
\includegraphics{employment_data_PSERV.pdf}
\includegraphics{employment_data_OTOT.pdf}}
\end{center}
\textit{{\protect\footnotesize {Data and estimated conditional means for four series from Illinois. Top left: unemployment rate; Top right: change in employees, Construction; Bottom left: change in employees, Prof \& Business services; Bottom right: nominal personal income.}}}
\end{figure}
We compare WAIC across a grid of $(s,r)$ specifications, see Table~\ref{tab:waic}. WAIC keeps improving, with diminishing and increasingly noisy returns, as $r$ increases for each value of $s$. $p_{\text{waic}}$ fluctuates by several thousand across adjacent grid points, so WAIC alone does not identify a stable optimum. As the diagnostics below show, its improvement at large $r$ reflects the latent factors absorbing individual states' idiosyncratic dynamics rather than picking up common variation across states. We therefore do not select $(s,r)$ by minimizing WAIC directly, and instead combine it with two diagnostics calculated from the posterior draws. First, under $s=4$ the posterior mean of $\pmb{\Omega}$'s fourth diagonal element is small relative to the other three and shrinks further with $r$, indicating a spurious fourth row factor. Therefore, we restrict attention to $s=2,3$. Second, for both $s=2$ and $s=3$ we track each column of $\mathbf{W}$'s concentration ratio, which is the largest squared state loading divided by the column's total squared loading. For $r=1$ to $r=11$, every column's ratio stays below 0.2 for both $s$. From $r=14$ on, at least one column jumps above 0.4, dominated by a single state (Hawaii for $s=2$; persistently Texas for $s=3$ from $r=18$ on). This pattern replicates across independently reinitialized runs. We select $r=11$, the largest value before this appears. WAIC's real, non-noise improvement is concentrated in the move from $r=7$ to $r=9$--$11$ (roughly 5,000 points), while $r=9$ to $r=11$ is itself flat. Little genuine fit is sacrificed by stopping there. At the value of $r=11$, $s=2$ and $s=3$ are statistically tied on WAIC (613,273 versus 613,239). We focus on $s=3$ because it separates the panel into three economically distinct row factors -- unemployment, sector employment together with housing prices, and income/output. A model with $s=2$ cannot represent these factors without conflating two of them.
\begin{table}[!t]
\caption{WAIC across left-factor ($s$) and right-factor ($r$) specifications.} \label{tab:waic}
\begin{center}
\begin{tabular}{crrr}
\toprule
$r$ & $s=2$ & $s=3$ & $s=4$ \\
\midrule
7 & 620{,}066 & 618{,}033 & 618{,}527 \\
9 & 614{,}629 & 613{,}187 & 614{,}236 \\
11 & 613{,}273 & 613{,}239 & 613{,}354 \\
14 & 604{,}247 & 603{,}106 & 609{,}915 \\
16 & 597{,}263 & 597{,}510 & 606{,}927 \\
18 & 597{,}207 & 598{,}783 & 602{,}451 \\
20 & 594{,}151 & 594{,}732 & 599{,}441 \\
22 & 589{,}687 & 591{,}284 & 595{,}375 \\
\bottomrule
\end{tabular}
\end{center}
\end{table}
We estimate the model with a $3 \times 11$ latent factor matrix ($s=3$, $r=11$). Consequently, the
600 observed state-level series are summarized by 33 dynamic factors. The purpose
of the empirical exercise is to illustrate that the matrix state space structure
can capture the main common movements in a large mixed-frequency panel while
remaining computationally tractable.
\begin{figure}[!t]
\caption{Estimated factor loadings, most persistent state factor, and regime probabilities.} \label{fig:loadings_regime}
\begin{center}
\includegraphics[width=0.80\textwidth]{loadings_regime_combined.pdf}
\end{center}
\textit{{\protect\footnotesize {Estimated factor loadings, most persistent state factor, and regime probabilities. Top left: estimated left factor loadings $\pmb{\Lambda}$; Top right: estimated right factor loadings $\mathbf{W}$; Bottom left: posterior mean of the most persistent row factor, selected by an AR(1)-times-posterior-scale score; Bottom right: smoothed probability of a high or low variance state.}}}
\end{figure}
Figure \ref{fig:data_fit} reports four representative series from Illinois together with their estimated conditional means from the posterior draws. The examples include the unemployment rate, changes in the total employees in the construction sector, changes in the total employees in the professional and business services sector, and quarterly nominal personal income. The fitted conditional means track the persistent movements in the data while smoothing through high-frequency
idiosyncratic variation. This is especially clear in the employment series, where the observed monthly changes are volatile but the fitted component captures the lower-frequency movement. The model also captures the major business-cycle events in the sample, including the Great Recession and the sharp movements around the COVID period, while treating the largest transitory observations as measurement
noise or heavy-tailed shocks rather than forcing them entirely into the common state.
\begin{table}[!t]
\caption{Integrated autocorrelation times (IACT) and effective sample sizes (ESS), $s=3$, $r=11$, 50{,}000 post-burn-in draws.} \label{tab:iact}
\begin{center}
\begin{tabular}{lrrr}
\toprule
Quantity & \# params & IACT & ESS \\
\midrule
$\text{diag}(\pmb{\Omega})$, factor 1 & 1 & 65.6 & 762 \\
$\text{diag}(\pmb{\Omega})$, factor 2 & 1 & 65.9 & 759 \\
$\text{diag}(\pmb{\Omega})$, factor 3 & 1 & 5.4 & 9{,}316 \\
$\pmb{\Lambda}$ loading, factor 1 & 1 & 26.3 & 1{,}905 \\
$\pmb{\Lambda}$ loading, factor 2 & 1 & 29.8 & 1{,}677 \\
$\pmb{\Lambda}$ loading, factor 3 & 1 & 17.7 & 2{,}818 \\
$\text{tr}(\pmb{\Psi})$ & 1 & 53.3 & 938 \\
max eig$(\pmb{\Psi})$ & 1 & 39.8 & 1{,}258 \\
$\text{diag}(\mathbf{U})$, all factors & 11 & 49.8--63.4 & 789--1{,}005 \\
Idiosyncratic variances, all series & 12 & 9.6--41.3 & 1{,}209--5{,}222 \\
$\pmb{\Phi}_1,\pmb{\Phi}_2$ coefficients & 18 & 1.2--21.5 & 2{,}330--41{,}201 \\
$\mathbf{W}$ anchor loadings, all columns & 11 & 1.1--75.7 & 660--44{,}089 \\
\bottomrule
\end{tabular}
\end{center}
\end{table}
Figure \ref{fig:loadings_regime} summarizes the estimated factor structure and the probabilities for the two-state Markov switching variance. The posterior mean of the left loading matrix shows that the first row factor loads primarily on the two lower-frequency income and output series (nominal personal income and real gross state product), the second row factor loads broadly across the sector-level employment series together with the house price index, and the third row factor loads almost exclusively on the unemployment rate, with an offsetting loading on house prices. The posterior mean of the right loading matrix shows substantial heterogeneity across states, with the first column factor capturing a common national component (positive for all 50 states) and the remaining column factors capturing cross-state deviations from that common movement. The bottom-left panel plots the most persistent row factor by our AR(1)-times-scale ranking, the unemployment-loaded factor. It has a clear business-cycle interpretation, rising sharply during the Great Recession and again, briefly, during the COVID episode. The smoothed probability of the high-volatility transition state rises sharply and stays elevated for several months around exactly the 2008--09 financial crisis and the 2020 COVID episode. This is consistent with the role of $\kappa_t$ as a means of capturing common rare, large aggregate shocks.
Table \ref{tab:iact} reports integrated autocorrelation times (IACT) and effective sample sizes (ESS), computed via \cite{Geyer(92)}'s initial monotone sequence estimator, for a representative set of parameters. Individually interpretable, low-dimensional quantities -- the diagonal of $\pmb{\Omega}$, the $\pmb{\Lambda}$ anchor loadings, and the scale of $\pmb{\Psi}$ -- are reported separately. The remaining blocks, which have no individual economic interpretation and are too numerous to report one by one. Instead, we summarize them by the range of IACT/ESS values across each block. The diagonal entries of $\mathbf{U}$, the idiosyncratic variances, the autoregressive coefficients $\pmb{\Phi}_1,\pmb{\Phi}_2$, and the right factor loadings $\mathbf{W}$. The slowest-mixing quantities are the two largest diagonal elements of $\pmb{\Omega}$, with IACT $\approx 66$ and ESS $\approx 760$ out of 50{,}000 post-burn-in draws. Every other block mixes substantially faster. This confirms that the partially collapsed Gibbs sampler, combined with PXDA and the Metropolis--Hastings moves on the Stiefel manifold, delivers adequate mixing for this specification.
Overall, the empirical results indicate that a small matrix of latent factors can summarize a large panel of state-level macroeconomic series. The row factors capture comovement across economic variables, the column factors capture cross-state dependence, and the Markov switching transition scale identifies periods in which common shocks are unusually large.
\section{Conclusion}
We developed the Kalman filter, smoother, and posterior sampling algorithms for a large class of conditionally, linear matrix normal time series models.
There are additional algorithms that could be derived for this class of time series models that we have omitted due to space constraints. These include disturbance smoothing and disturbance simulation smoothing, see \cite{Koopman(93)} and \cite{deJongShephard(95)}. One can also develop methods for exact treatment of initial conditions in non-stationary models, see \cite{deJong(91)} and \cite{Koopman(97)}. Finally, the matrix Kalman filter can be used within a particle filter, see \cite{ChenLiu(00)}.
\subsection{Declaration of generative AI and AI-assisted technologies in the manuscript preparation process}
During the preparation of this work, the authors used Codex to proofread the paper and computer code. The authors reviewed and edited the output as needed and take full responsibility for the content of the published article.
\bibliographystyle{jf}
\bibliography{creal}