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.
173,877 characters
Latent Variable Modelling by Supervised Diffusion
\maketitle
\begin{abstract}
This paper proposes a new methodological framework for estimating
inferential models with latent variables. It also introduces a new
latent variable regression model called LARX: an extension of the
ubiquitous autoregressive model with exogenous inputs (ARX) in which any
or all input variables can be latent. In deriving the LARX model, a
minor contribution is also made to the field of matrix calculus: A new
matrix operator is defined and applied to solve a class of Lagrangian
optimisation problems with interactions between multiple coefficient
vectors subject to case-by-case constraints.
In the empirical section, the LARX model is used to re-examine the
relationship between stock market performance and real economic activity
in the United States. The LARX model attains an out-of-sample R-squared
of up to 79.7\% compared to 50.3\% for the baseline OLS approach. It also
reveals new information about the underlying drivers of the relationship
between stock returns and economic growth, including the predictive
power of sector rotations.
\end{abstract}
\section{Introduction}
\label{introduction}
A latent variable (LV) is an unobserved process which influences the
behaviour of an observed system. Many inferential models in economics
and finance use latent variables, which adds a layer of complexity to
their empirical estimation process. Prominent examples of such models
include the business cycle theory (\citep{Burns-Mitchell-1946}) and the
asset pricing theory (APT) (\citep{Ross-1976}), both of which are
characterised by a small number of unobserved processes driving a much
larger number of observed variables\footnote{The business cycle theory states that the unobserved state of the
business cycle drives a large number of observed macroeconomic
aggregates such as growth, inflation and unemployment. In APT, long-term
asset returns are determined by the respective assets' exposures to
unobserved risk premia.}.
One of the most popular methodologies for approximating latent variables
in both descriptive and inferential settings is principal component
analysis (PCA) attributed to \citep{Pearson-1901} and
\citep{Hotelling-1933}. PCA uses a group of observed variables to
construct a sequence of latent variables, known as Principal Components
(PC) or Diffusion Indices (DI)\footnote{Principal components are referred to as ``diffusion indices'' in
\citep{Stock-2002}, which is one of the first papers to demonstrate the
superior forecasting power of latent macroeconomic factors approximated
by principal components with respect to a wide range of US economic
aggregates. The term is not to be confused with the colloquial meaning
of a diffusion index as a summary statistic describing the share of
positive and negative values in a set of survey responses or
macroeconomic indicators, such as \href{https://www.ismworld.org/supply-management-news-and-reports/reports/ism-pmi-reports/}{the suite of US business activity
indices produced by the Institute of Supply Management}. This paper
adopts the terminology of \citep{Stock-2002} and uses the term
``diffusion index'', or DI, to refer to a linear combination of observed
variables constructed with the goal of approximating a latent variable
or factor.}, each capturing one independent
direction of variance. DIs have useful statistical properties, including
their ability to reduce the dimensionality of the input dataset and to
eliminate multicollinearity in the final regression problem. However,
PCA also has two major drawbacks: interpretability, and the restrictive
assumptions required to ensure the optimality of DI-based factor models
in inferential settings (e.g., see \citep{Bai-2006}).
This paper builds on the conceptual foundations of PCA and its
multivariate cousin, Canonical Correlation Analysis (CCA)
(\citep{Hotelling-1936}), to propose a new methodological framework for
estimating inferential models with latent variables. The new framework
is given the name `'Supervised Diffusion'' (SDF), following the
terminology of \citep{Stock-2002} which refers to principal components
as diffusion indices. SDF capitalises on the idea of latent variables as
linear combinations of observed variables. However, contrary to PCA, the
linear combination weights are optimised for inferential purposes.
The SDF formula makes it possible to transform any traditional
regression methodology into a latent variable model. As a practical
example, the ubiquitous autoregressive model with exogenous inputs (ARX)
is used to derive a new regression model called LARX: a superset of the
ARX model in which any or all variables can be latent. Like ARX, the
LARX model presents with several interesting special cases, including a
parsimonious lead-lag regression model referred to as ``Latent Shock
Regression'' (LSR), as well as a new canonical decomposition technique
best described as Canonical Autocorrelation Analysis (CAA).
In the broader ecosystem of latent variable research, SDF and the LARX
model contribute to a body of work extending beyond economics and
finance, where models based on CCA and its close relative Partial Least
Squares (PLS) (\citep{Wold-1982,Wold-1975}) are studied under the term
Latent Variable Regression (LVR). Over the past 90 years, LVR models
have been developed and applied in fields like industrial chemistry
(e.g., \citep{Burnham-MacGregor-1996,Burnham-McGregor-1999}), machine
learning (e.g., \citep{Wang-2020,Dai-2020,Vaerenberg-2018,Chi-2013}), and medicine, among other quantitative disciplines
(\citep{Uurtio-2017} provides a good overview). Traditional use cases
for these models include dimensionality reduction, identification of the
directions of correlation in multivariate data streams (e.g.,
\citep{Burnham-MacGregor-1996,Dong-2018,Qin-2021}), and multi-label
classification of images, videos, audio, and hand-written text (e.g.,
\citep{Wang-2020,Dai-2020,Vaerenberg-2018,Chi-2013}).
The rest of this paper is organised as follows. Section
\ref{blockwise_direct_sum} defines a new matrix operator which is
subsequently used to derive the LARX model. Section \ref{sdf} introduces the
supervised diffusion framework and examines its relationship to PCA,
least-squares regression and portfolio optimisation. Section \ref{larx}
derives the LARX methodology, making a minor contribution to the field
of matrix calculus. Section \ref{larx_special_cases} examines a few special
cases of the LARX model, including CCA, LSR and CAA. Section
\ref{empirical_results} applies the LARX model in an empirical setting by
examining the predictive power of the US stock market with respect to
real US economic activity. Section \ref{concluding_remarks} concludes and
discusses possible avenues for future work.
\section{Block-wise Direct Sum Operator}
\label{blockwise_direct_sum}
Some of the mathematical derivations in this paper involve block-wise
operations between matrices and vectors with blocks along one dimension.
These derivations become much easier with the help of a shorthand
notation which we can call a block-wise direct sum \(\mathbf{\bm{A}}^{\oplus}\) for an
arbitrary matrix or vector \(\mathbf{\bm{A}}\). This section briefly introduces the
block-wise direct sum operator and its relevant properties.
Let \(\left< {\mathbf{\bm{A}}_i | 1 \leq i \leq k} \right>\) be a sequence of \(k\) real
matrices with arbitrary dimensions. If all \(\mathbf{\bm{A}}_i\) have the same number
of columns, they can be concatenated vertically into a matrix with row
blocks. If all \(\mathbf{\bm{A}}_i\) have the same number of rows, they can be
concatenated horizontally into a matrix with column blocks. Formally:
\begin{flalign*}
& \left< {\mathbf{\bm{A}}} \right> \equiv \left< {\mathbf{\bm{A}}_i} \right> := \left< {\mathbf{\bm{A}}_i | 1 \leq i \leq k} \right>
& \text{ a sequence of matrices of arbitrary dimensions} \\[10pt]
& \mathbf{\bm{A}} \equiv \mathbf{\bm{A}}_v := \left[ \left< {\mathbf{\bm{A}}_i} \right> \right]_v
& \text{ a vertical concatenation of } \left< {\mathbf{\bm{A}}_i} \right> \text{ (default)} \\[10pt]
& \mathbf{\bm{A}}_h := \left[ \left< {\mathbf{\bm{A}}_i} \right> \right]_h
& \text{ a horizontal concatenation of } \left< {\mathbf{\bm{A}}_i} \right>
\end{flalign*}
The block-wise direct sum operator \(\mathbf{\bm{A}}^{\oplus}\) is then defined as:
\begin{flalign*}
& \mathbf{\bm{A}}^{\oplus} \equiv \left< {\mathbf{\bm{A}}} \right>^{\oplus} \equiv \left< {\mathbf{\bm{A}}_i} \right>^{\oplus}
:= \mathbf{\bm{A}}_1 \oplus \mathbf{\bm{A}}_2 \oplus \mathbf{\bm{A}}_3 \oplus \cdots \oplus \mathbf{\bm{A}}_k
= \begin{pmatrix}
\mathbf{\bm{A}}_1, & \0, & \cdots & \0 \\
\0, & \mathbf{\bm{A}}_2, & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0, & \0, & \cdots & \mathbf{\bm{A}}_k
\end{pmatrix} &
\end{flalign*}
Compatibility of dimensions is not a constraint for a direct sum between
matrices, so a block-wise direct sum is defined for any sequence of
matrices \(\left< {\mathbf{\bm{A}}_i} \right>\). The result is always a matrix whose number of
rows (columns) is the sum total number of rows (columns) across all
comprising \(\mathbf{\bm{A}}_i\).
The matrix \(\mathbf{\bm{A}}^{\oplus}\) can itself be mapped to either a sequence of \(k\)
row blocks or a sequence of \(k\) column blocks without slicing through
the original matrices in the sequence. For most use cases the block
structure of \(\mathbf{\bm{A}}^{\oplus}\) will not be relevant. For the remaining
scenarios let us apply the same shorthand notation as above:
\begin{flalign*}
& \mathbf{\bm{A}}^{\oplus} \equiv \mathbf{\bm{A}}^{\oplus}_v
:= \left( { \begin{array}{cccc}
\mathbf{\bm{A}}_1, & \0, & \cdots & \0 \\ \hline
\0, & \mathbf{\bm{A}}_2, & \cdots & \0 \\ \hline
\vdots & \vdots & \ddots & \vdots \\ \hline
\0, & \0, & \cdots & \mathbf{\bm{A}}_k
\end{array} } \right), \quad
\mathbf{\bm{A}}^{\oplus}_h
:= \left( { \begin{array}{c|c|c|c}
\mathbf{\bm{A}}_1 & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{A}}_2 & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{A}}_k \\
\end{array} } \right) &
\end{flalign*}
One of the main use cases for the block-wise direct sum operator lies in
defining block-wise inner products and block-wise quadratic forms between
matrices and vectors. For example, let \(\mathbf{\bm{B}}\) be a matrix with the same
number of rows and the same block structure as \(\mathbf{\bm{A}}\). The product
\(\left( \mathbf{\bm{A}}^{\oplus} \right)' \mathbf{\bm{B}}\) then yields:
\begin{flalign*}
& \left( \mathbf{\bm{A}}^{\oplus} \right)' \mathbf{\bm{B}}
= \begin{pmatrix}
\mathbf{\bm{A}}_1' & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{A}}_2' & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{A}}_k'
\end{pmatrix} \left( \begin{array}{c}
\mathbf{\bm{B}}_1 \\ \hline
\mathbf{\bm{B}}_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{B}}_k
\end{array} \right)
= \left( \begin{array}{c}
\mathbf{\bm{A}}_1' \mathbf{\bm{B}}_1 \\ \hline
\mathbf{\bm{A}}_2' \mathbf{\bm{B}}_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{A}}_k' \mathbf{\bm{B}}_k
\end{array} \right) &
\end{flalign*}
In our case, this operation will be particularly useful for solving
Lagrangian optimisation problems in the presence of piecemeal
constraints on the coefficient vector, e.g., when one portion of the
coefficient vector must have a unit length constraint, another a
zero-sum constraint, and so on. Several properties of the block-wise
direct sum operator will be relevant in this context. First of all, we
note that the transpose of \(\mathbf{\bm{A}}^{\oplus}\) is the same as the block-wise
direct sum of \(\mathbf{\bm{A}}'\):
\begin{proposition} \label{prop:bds_vs_transpose} The function
composition of the block-wise direct sum operator and the transpose
operator is commutative. In other words, \( \left( \mathbf{\bm{A}}^{\oplus} \right)' = \left( { \mathbf{\bm{A}}'
} \right)^{\oplus}\).
\end{proposition}
\begin{proof}
\begin{flalign*}
& {\left( \mathbf{\bm{A}}^{\oplus} \right)'} = \begin{pmatrix}
\mathbf{\bm{A}}_1' & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{A}}_2' & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{A}}_k'
\end{pmatrix}
= \mathbf{\bm{A}}_1' \oplus \mathbf{\bm{A}}_2' \oplus \mathbf{\bm{A}}_3' \oplus \cdots \oplus \mathbf{\bm{A}}_k'
= \left( {\mathbf{\bm{A}}'} \right)^{\oplus} &
\end{flalign*} \end{proof}
Second, the block-wise direct sum of a vector can be written as a
block-wise Kronecker product as defined in \citep{Khatri-Rao-1968},
which we denote by ``\(\odot\)'':
\begin{proposition} \label{prop:bds_vs_odot} Let \(\mathbf{\bm{a}}\) comprise \(k\)
blocks given by the sequence of vectors \(\left< {\mathbf{\bm{a}}_i | 1 \leq i \leq
k} \right>\). The matrix \(\mathbf{\bm{a}}^{\oplus}\) can be expressed as a block-wise Kronecker
product between \(\mathbf{\bm{a}}\) and an identity matrix \(I_k\) with vector
blocks along the same dimension as \(\mathbf{\bm{a}}\).
\end{proposition}
\begin{proof}For a column vector we have:
\begin{flalign*}
& \mathbf{\bm{a}}^{\oplus} = \begin{pmatrix}
\mathbf{\bm{a}}_1 & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{a}}_2 & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{a}}_k
\end{pmatrix} = \left( { \begin{array}{c}
\mathbf{\bm{a}}_1 \otimes \begin{bmatrix}1, & 0, & \cdots & 0 \end{bmatrix} \\ \hline
\mathbf{\bm{a}}_2 \otimes \begin{bmatrix}0, & 1, & \cdots & 0 \end{bmatrix} \\ \hline
\vdots \\ \hline
\mathbf{\bm{a}}_k \otimes \begin{bmatrix}0, & 0, & \cdots & 1 \end{bmatrix}
\end{array} } \right) = \left( { \begin{array}{c}
\mathbf{\bm{a}}_1 \\ \hline
\mathbf{\bm{a}}_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{a}}_k
\end{array} } \right) \odot \left( { \begin{array}{cccc}
1, & 0, & \cdots & 0 \\ \hline
0, & 1, & \cdots & 0 \\ \hline
\vdots & \vdots & \ddots & \vdots \\ \hline
0, & 0, & \cdots & 1
\end{array} } \right) = \mathbf{\bm{a}} \odot I_k &
\end{flalign*}
The proof for a row vector follows by symmetry. \end{proof}
Third, for two column vectors \(\mathbf{\bm{a}}\) and \(\bm{b}\) with the same length
and row block structure, the operation \(\left( \mathbf{\bm{a}}^{\oplus} \right)' \bm{b}\) is symmetric:
\begin{proposition} \label{prop:bds_vs_vectors}
For column vectors \(\mathbf{\bm{a}}\) and \(\bm{b}\) with identically sized row blocks
\(\mathbf{\bm{a}}_1,\mathbf{\bm{a}}_2,\mathbf{\bm{a}}_3,\ldots,\mathbf{\bm{a}}_k\) and \(\bm{b}_1,\bm{b}_2,\bm{b}_3,\ldots,\bm{b}_k\),
respectively, \(\left( \mathbf{\bm{a}}^{\oplus} \right)' \bm{b} = \left( \bm{b}^{\oplus} \right)' \mathbf{\bm{a}} \).
\end{proposition}
\begin{proof}
\begin{flalign*}
& \left( \mathbf{\bm{a}}^{\oplus} \right)' \bm{b} = \begin{pmatrix}
\mathbf{\bm{a}}_1' & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{a}}_2' & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{a}}_k'
\end{pmatrix} \left( \begin{array}{c}
\bm{b}_1 \\ \hline
\bm{b}_2 \\ \hline
\vdots \\ \hline
\bm{b}_k
\end{array} \right)
= \left( \begin{array}{c}
\mathbf{\bm{a}}_1' \bm{b}_1 \\ \hline
\mathbf{\bm{a}}_2' \bm{b}_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{a}}_k' \bm{b}_k
\end{array} \right)
= \left( \begin{array}{c}
\bm{b}_1' \mathbf{\bm{a}}_1 \\ \hline
\bm{b}_2' \mathbf{\bm{a}}_2 \\ \hline
\vdots \\ \hline
\bm{b}_k' \mathbf{\bm{a}}_k
\end{array} \right) = \left( \bm{b}^{\oplus} \right)' \mathbf{\bm{a}} &
\end{flalign*} \end{proof}
Lastly, the block-wise direct sum operator is commutative with respect to
a certain class of operations over matrix sequences. Specifically, for a
sequence of matrices \(\left< {\mathbf{\bm{A}}_i | 1 \leq i \leq k} \right>\) and an operation
\(f\) over matrix sequences of length \(k\), it can be shown that \(f
\left( { \left< {\mathbf{\bm{A}}} \right> } \right)^{\oplus} = \left[ f \left( { \left< {\mathbf{\bm{A}}^{\oplus}} \right> } \right) \right]\)
if \(f\) satisfies certain conditions. This, in turn, can be used to
prove that for two vectors \(\mathbf{\bm{a}}\) and \(\bm{b}\) with \(k\) row blocks
each, the block-wise Kronecker product \(\mathbf{\bm{a}} \odot \bm{b}\) can be factorised
in the same way as the traditional Kronecker product, namely:
\begin{equation*}
\mathbf{\bm{a}} \odot \bm{b} = \left( {\mathbf{\bm{a}} \odot I_{\bm{b}}} \right) \bm{b} = \left( { I_{\mathbf{\bm{a}}} \odot \bm{b} } \right) \mathbf{\bm{a}}
\end{equation*}
where \(I_{\mathbf{\bm{a}}}\) and \(I_{\bm{b}}\) are identity matrices with the same row
block structure as \(\mathbf{\bm{a}}\) and \(\bm{b}\), respectively. The corresponding
derivations are deferred to \ref{bds_vs_mmul} and \ref{block_kron_factorisation}.
\section{A Conceptual Case for a Supervised Diffusion Framework}
\label{sdf}
An inferential model with latent variables (from now, simply a ``latent
variable model'', or LVM for short) has two simultaneous objectives: to
estimate a functional relationship, and to approximate the latent
variables in it. Several types of latent variable models are widely used
in economics and finance\footnote{Two popular latent variable models in modern econometrics are
instrumental variable regression (e.g., \citep{Stock-2003-IV}) and
hidden Markov models (\citep{Baum-1966}). In the former, the causal
relationship between the explanatory and the dependent is obscured by an
unobserved common factor. In the latter, the variables and their
relationships are influenced by unobserved changes in regime such as the
stage of the business cycle or a bull/bear market.}, but this paper is concerned with a
special class of LVMs in which unobserved processes are approximated as
linear combinations of observed variables. Principal Component Analysis
(PCA) (\citep{Pearson-1901,Hotelling-1933}) is arguably the best-known
model of this kind.
This section proposes a new methodological framework for estimating
latent variable models. It is referred to as a ``Supervised Diffusion
Framework'' (SDF) following the terminology of \citep{Stock-2002} which
uses the term ``Diffusion Index'' (DI) to describe a composite indicator
obtained by PCA. Conceptually, SDF brings PCA into an inferential
setting: If the goal of PCA is to describe the latent sources of
variation in a single set of variables, the goal of SDF is to describe
the latent source(s) of functional dependency between two or more sets.
In this context, SDF also generalises the method of Canonical
Correlation Analysis (CCA) (\citep{Hotelling-1936}) whose goal is to
identify the latent sources of pairwise correlation -- a special type of
functional relationship (more on this in Section \ref{cca}).
Let us start with a couple of formal definitions:
\begin{definition}\label{def:lv}
A latent variable \(\tilde{y}\) is an unobserved process which can be
approximated as a linear combination over a vector of observed random
variables \(Y = \begin{bmatrix} y_1, & y_2, & y_3, & \ldots, & y_n
\end{bmatrix}\) with measurement error \(r_y\), such that \(\tilde{y} = Y\bm{\mathbf{w}} +
r_y\) where \(\bm{\mathbf{w}}\) has dimensions \(n \times 1\) and can't be a multiple
of a standard basis vector.
\end{definition}
\begin{definition}\label{def:lvm}
A latent variable model (LVM) is a regression model in which at least
one variable is latent.
\end{definition}
Now let us consider a generalised formula for a standard multiple
regression between a dependent variable \(y\) and a sequence of \(\mathbf{\kappa}\)
explanatory variables \(\left< {x_j} \right> := \left< {x_1,x_2,x_3,\ldots,x_{\mathbf{\kappa}}} \right> =
\left< {x_j} \right>_{j=1}^{\mathbf{\kappa}}\):
\begin{equation}\label{eqn:reg}
y = F_{\bm{\gamma}}(x_1,x_2,x_3,\ldots,x_{\mathbf{\kappa}}) + e \equiv F_{\bm{\gamma}}(\left< {x_j} \right>) + e
\end{equation}
Here, \(\bm{\gamma}\) denotes a vector of regression parameters, \(e\) is a
mean-zero error term, and \(F_{\bm{\gamma}}: \mathbb{R} \longrightarrow
\mathbb{R}\) is generally assumed to be at least once differentiable
with respect to \(\bm{\gamma}\). The standard least-squares solution to (\ref{eqn:reg})
is obtained by finding the vector \(\bm{\hat{\gamma}}\) which minimises the variance
of \(e\) using some sample observation vectors \(\mathbf{\bm{y}}\) and \(\mathbf{\bm{x}}_j\)
containing the data for \(y\) and \(x_j\), respectively.
In a slight abuse of notation, let \(F(\left< {\mathbf{\bm{a}}_j} \right>)\) represent a
row-wise operation over the sequence \(\left< {\mathbf{\bm{a}}_j} \right>\) when \(\mathbf{\bm{a}}_j\) are
sample vectors (or matrices). In other words, given a sequence of sample
vectors \(\left< {\mathbf{\bm{a}}_j} \right>_{j = 1}^{\mathbf{\kappa}}\) where \(\mathbf{\bm{a}}_j = \begin{bmatrix}
a_{1,j}, & a_{2,j}, & a_{3,j}, & \ldots, & a_{s,j} \end{bmatrix}'\) for
all \(j\) and \(s\) is the sample size, let:
\begin{equation*}
F(\left< {\mathbf{\bm{a}}_j} \right>) = \begin{bmatrix}
F(\left< {a_{1,j}} \right>_{j = 1}^{\mathbf{\kappa}}), & F(\left< {a_{2,j}} \right>_{j = 1}^{\mathbf{\kappa}}), &
F(\left< {a_{3,j}} \right>_{j = 1}^{\mathbf{\kappa}}), & \ldots, & F(\left< {a_{s,j}} \right>_{j = 1}^{\mathbf{\kappa}})
\end{bmatrix}'
\end{equation*}
We can then write the least squares optimisation problem for (\ref{eqn:reg})
in sample form as:
\begin{equation}\label{eqn:reg_solution}
\bm{\hat{\gamma}} = \underset{\bm{\gamma}}{\mathrm{argmin}} \| \mathbf{\bm{y}} - F_{\bm{\gamma}}(\left< {\mathbf{\bm{x}}_j} \right>) \|^2_2
\end{equation}
Now assume WLOG that all the variables in this model are latent. The
problem specified by (\ref{eqn:reg}) then becomes:
\begin{equation}\label{eqn:reg_lv_a}
\tilde{y} = F_{\bm{\gamma}}(\left< {\tilde{x}_j} \right>) + e
\end{equation}
According to Definition \(\ref{def:lv}\), we can approximate \(\tilde{y}\) and
all \(\tilde{x}_j\) using linear combinations of observed variables. Define
\(Y = \begin{bmatrix} y_1, & y_2, & y_3, & \ldots, & y_n \end{bmatrix}\)
and \(X_j = \begin{bmatrix} x_{1,j}, & x_{2,j}, & x_{3,j}, & \ldots, &
x_{m_j,j} \end{bmatrix}\) such that \(Y \bm{\mathbf{w}} = \tilde{y} + r_y\) and \(X_j \bm{\omega}_j
= \tilde{x}_j + r_j\) for some unknown weight vectors \(\bm{\mathbf{w}}\) and \(\bm{\omega}_j\) with
dimensions \(n \times 1\) and \(m_j \times 1\), respectively. Plugging
these expressions into (\ref{eqn:reg_lv_a}) and collecting the error terms
produces a regression problem of the form:
\begin{equation}\label{eqn:reg_lv_b}
\begin{split}
& Y \bm{\mathbf{w}} = F_{\bm{\gamma}}(\left< {X_j \bm{\omega}_j} \right>) + \epsilon \\
& \text{where } \epsilon = e + r_y + F_{\bm{\gamma}}(\left< {\tilde{x}_j} \right>) - F_{\bm{\gamma}}(\left< {X_j \bm{\omega}_j} \right>)
\end{split}
\end{equation}
Here, \(\epsilon\) collects all the possible sources of regression error
in (\ref{eqn:reg_lv_a}): \(e\) is the residual noise from the functional
relationship \(F_{\bm{\gamma}}: \left< {\tilde{x}_j} \right> \mapsto \tilde{y}\) under the true
regression parameters, \(r_y\) is the measurement noise from
approximating \(\tilde{y}\) based on \(Y\), and \(F_{\bm{\gamma}}(\left< {\tilde{x}_j} \right>) -
F_{\bm{\gamma}}(\left< {X_j \bm{\omega}_j} \right>)\) summarises the reduction in the quality of
regression fit caused by the measurement noise in all \(\tilde{x}_j\).
When we compare the traditional regression problem (\ref{eqn:reg}) to the
latent variable regression problem (\ref{eqn:reg_lv_b}), two important
differences jump out. First, to the extent that the statistical
properties of any regression model are defined by the properties of its
error term, (\ref{eqn:reg_lv_b}) looks different from (\ref{eqn:reg}) because of all
the additional components in \(\epsilon\). By extension, the
optimisation techniques which are appropriate for (\ref{eqn:reg}) may not
always be applicable to (\ref{eqn:reg_lv_b}). However, a good counter-example
to this argument is provided by any traditional regression model in
which the sample measurements of the observed variables \(y\) and
\(x_j\) contain measurement uncertainty. In these cases one could just
as well write out \(y\) and \(x_j\) as \(y = \tilde{y} + r_y\) and \(x_j =
\tilde{x}_j + r_j\) for some ``true'' intrinsic processes \(\tilde{y}\) and
\(\tilde{x}_j\) with the measurement noise removed\footnote{An argument can be made that this counter-example represents
many research settings prevalent in finance and economics. For example,
asset prices are susceptible to noise trading, while macroeconomic
aggregates will generally have non-negligible measurement error.}. Plugging these
expressions back into (\ref{eqn:reg}) would produce the same error term as in
(\ref{eqn:reg_lv_b})\footnote{Under this assumption both (\ref{eqn:reg}) and (\ref{eqn:reg_lv_a}) can also
be interpreted as errors-in-variables (EIV) models
(\citep{Griliches-1970}). A thorough discussion of the statistical
properties of EIV models is beyond the scope of this paper, but the
reader is referred to \citep{Wolter-1982} as a starting point.}. This allows us to draw a strong methodological
parallel between traditional regression models and LVMs with limited
loss of generality. Formally speaking, for any intrinsic random
processes \(\tilde{y}\) and \(\tilde{x}_j\) governed by a hypothesised functional
relationship \(F_{\bm{\gamma}}\), any viable approach to estimating \(F_{\bm{\gamma}}\)
as a traditional regression model should also be applicable in a latent
variable setting as long as all the available sample measurements for
\(\tilde{y}\) and \(\tilde{x}_j\) are prone to measurement uncertainty.
The second difference between (\ref{eqn:reg}) and (\ref{eqn:reg_lv_b}) stems from
the extra parameters required for estimating (\ref{eqn:reg_lv_b}) (recall that
the coefficient vectors \(\bm{\mathbf{w}}\) and \(\bm{\omega}_j\) are all initially unknown).
One could say that the asset pricing theory (APT) has largely succeeded
in estimating \(\bm{\omega}_j\) using a combination of economically motivated
heuristics\footnote{Arguably the most famous example of a heuristic approach to the
construction of asset pricing factors is \citep{Fama-French-1993} which
recombines a large cross-section of single stocks into 25 value-weighted
portfolios based on valuation and market capitalisation, and then uses
those portfolios as both the dependent variables and the building blocks
for the Fama-French three-factor model. A similar technique was
subsequently adopted in other prominent papers, such as
\citep{Asness-2013} which constructs 48 portfolios across asset classes
ranked on valuation and price momentum (the individual portfolios are
value-weighted in the case of single stocks and equal-weighted in other
cases).} and principal component analysis\footnote{For example, \citep{Bai-2006} uses PCA to show that the
Fama-French three factor model (FF3F) (\citep{Fama-French-1993}) is a
reasonable approximation of the latent factor space in the US equity
market; \citep{Lustig-2011} demonstrates that the first two principal
components in the currency market map to a theoretically motivated
two-factor model consisting of a ``dollar'' (market) factor and a
``slope'' (carry) factor; \citep{Diebold-2006} proposes a PCA-based
factor model for the term structure of interest rates; etc.}. On the other
hand, the APT methodology does not lend itself well to other research
settings. To illustrate, let us consider a generic asset pricing formula
for a representative test asset \(\tilde{y}_i\):
\begin{equation}\label{eqn:reg_lv_fm}
\tilde{y}_i = c_i + \sum_{j = 1}^{k}{\beta_{i,j} \tilde{x}_j} + e,
\quad i \in \mathbb{N}
\end{equation}
Here, \(\tilde{x}_j\) are the latent asset pricing factors, \(\beta_{i,j}\)
are asset \(i\)'s factor loadings, and \(c_i\) is the asset-specific
intercept term. This class of models is characterised by two unique
properties: a simple functional relationship, and a large cross-section
of investment instruments \(\tilde{y}_i\) which may or may not be latent in
the sense of being linear combinations\footnote{Historically, asset pricing tests have used single stocks,
bonds, currency pairs, regional benchmark indices, and even custom test
portfolios (\citep{Fama-French-1993}, \citep{Lustig-2011} and
\citep{Asness-2013} are some of the most famous examples).}. The key assumption is
that all \(\tilde{y}_i\) are driven by a small number of latent systematic
risk factors \(\tilde{x}_j\) and some idiosyncratic noise. A breakthrough in
estimating this class of models was achieved by \citep{Bai-2006} with a
proof that the space of the true unobserved factors \(\tilde{x}_j\) is
consistently spanned by the set of diffusion indices \(\tilde{y}^{di}_i\)
obtained from a PCA decomposition\footnote{It can be noted that several papers have subsequently applied
various extensions to the PCA methodology for the purpose of estimating
latent factors more effectively. For example, \citep{Boucher-2021} uses
a sparse (penalised) version of PCA to approximate latent factors for
the European equity market, \citep{Begusic-2020} uses a clustering
technique to identify globally and locally relevant latent factors, and
\citep{Lettau-2020} develops a procedure called Risk Premium PCA
(RP-PCA) for estimating latent factors that better capture
cross-sectional risk premia in an asset pricing context.} over all \(\tilde{y}_i\), as long as
the cross-section of test assets is large enough. This result allows the
researcher to reduce the validation process for any economically
motivated asset pricing model to a three-step procedure: calculate a set
of DI factors \(\{\tilde{y}^{di}_i\}\), propose a set of economically
motivated heuristic factors \(\{\tilde{x}^h_j\}\), and show that \(\tilde{x}^h_j\)
and \(\tilde{y}^{di}_i\) span the same space. Formal statistical tests for
this procedure have been proposed by \citep{Bai-2006},
\citep{Parker-2016} and \citep{Andreou-2025}, among others.
Unfortunately, this three-step approach becomes far less effective when
the number of test assets \(\tilde{y}_i\) is small and fixed\footnote{\citep{Ahn-2018} shows that CCA is preferable to PCA under these
circumstances.}, or when
the functional relationship no longer reduces to a simple linear model.
The latter is a big stumbling block in macroeconomics: Like APT,
macroeconomic theory has long entertained the notion of a limited number
of unobserved processes governing a larger set of observed macroeconomic
aggregates\footnote{One prominent example of this paradigm is the theory of business
cycles which goes back as far as \citep{Burns-Mitchell-1946}.}; however, the functional relationships in
macroeconomics can often be non-linear or have an unknown lead-lag
dimension. Empirically, \citep{Stock-2002} applies PCA to 215
macroeconomic variables in an attempt to approximate a set of latent
factors driving 8 key economic aggregates. The forecasting performance
of the resulting diffusion indices is compared to that of the
economically motivated observed predictors using an autoregressive model
with exogenous inputs (ARX) up to lag order \(q\). The underlying
regression model can be expressed as:
\begin{equation}\label{eqn:reg_arx_di}
y_{t+1} = c + \sum_{\tau = 0}^{q}{\phi_{\tau} y_{t-\tau}} +
\sum_{j=1}^{\mathbf{\kappa}} \sum_{\tau=0}^{q} {\beta_{j,\tau} \tilde{x}^{di}_{j,t-\tau}} + e
\end{equation}
The results obtained by \citep{Stock-2002} suggest that PCA-based
diffusion indices constructed from a large cross-section of
macroeconomic variables can produce more accurate forecasts for key
macroeconomic aggregates than the individual observed predictors.
However, unlike in APT, there is no longer any guarantee that these
diffusion indices produce an optimal factor model. This point is best
illustrated with same example that \citep{Stock-2002} uses in support of
PCA, namely, that of the empirical Phillips curve. Broadly speaking, the
Phillips curve postulates that inflation can be predicted by real
economic activity; however, there is no consensus about which measure of
economic activity is the most appropriate\footnote{Alternative measures of real economic activity include GDP,
employment, capacity utilisation, business surveys, and various activity
``nowcasts'', to name a few.}. One could say that
applying PCA to a broad set of real activity measures negates the
question of variable selection altogether because it allows the
researcher to extract all the relevant sources of variation from the
real activity data. However, PCA also reveals absolutely nothing about
the relevance of each of those sources of variation for inflation
specifically. Just like with the observed measures of real activity,
inflation can be influenced by any one of the diffusion indices of real
activity or, indeed, any linear combination thereof. What's more, if no
diffusion indices can be removed from the pool of candidate regressors
out of hand, PCA becomes somewhat redundant altogether (\ref{pca_redundancy}
expresses this argument more formally).
An effective estimation procedure for the broad class of latent variable
models defined by (\ref{eqn:reg_lv_b}) would need to address the shortcomings
of both PCA and domain-specific heuristic approaches. To this end, this
paper proposes a relatively simple solution: simultaneously estimate all
the coefficients in (\ref{eqn:reg_lv_b}) using the method of least squares.
Formally, define:
\begin{itemize}
\item \(n \geq 1\) proxy measurements for the dependent variable \(\tilde{y}\),
collected into a \(1 \times n\) row vector \(Y\)
\item \(m_j \geq 1\) proxy measurements for each explanatory variable
\(\tilde{x}_j\) where \(1 \leq j \leq \mathbf{\kappa}\), collected into a sequence of
row vectors \(\left< {X_j} \right> \equiv \left< {X_j} \right>_{j=1}^{\mathbf{\kappa}}\)
\item An \(s \times n\) matrix of sample observations \(\mathbf{\bm{Y}}\) for the
variables in \(Y\) and a sequence of sample observation matrices
\(\left< {\mathbf{\bm{X}}_j} \right>\), each with dimensions \(s \times m_j\), for the
explanatory variable vectors \(X_j\).
\end{itemize}
The least-squares solution to (\ref{eqn:reg_lv_b}) in sample form is then
obtained by the vectors \(\bm{\hat{\gamma}}\), \(\bm{\mathbf{\hat{w}}}\) and \(\left< {\bm{\hat{\omega}}_j} \right> \equiv
\left< {\bm{\hat{\omega}}_j} \right>_{j=1}^{\mathbf{\kappa}}\) which satisfy:
\begin{equation}\label{eqn:sdf_solution}
\bm{\hat{\gamma}},\bm{\mathbf{\hat{w}}},\left< {\bm{\hat{\omega}}_j} \right> = \underset{\bm{\gamma},\bm{\mathbf{w}},\left< {\bm{\omega}_j} \right>}{\mathrm{argmin}}
\| \mathbf{\bm{Y}} \bm{\mathbf{w}} - F_{\bm{\gamma}}(\left< {\mathbf{\bm{X}}_j \bm{\omega}_j} \right>) \|^2_2
\end{equation}
In practice, we would often need to impose constraints on \(\bm{\mathbf{w}}\), \(\bm{\omega}\)
and/or \(\bm{\gamma}\) because the economic interpretation of the model or the
underlying variables implies certain restrictions (e.g., an asset
pricing factor may have to be a zero-cost investment portfolio), or
because we wish to avoid trivial solutions such as \(\bm{\mathbf{w}} = \bm{0}\). For
example, PCA and partial least squares (PLS) models
(\citep{Wold-1975,Wold-1982}) impose a unit length constraint on \(\bm{\mathbf{w}}\)
and \(\bm{\omega}\), while CCA is implemented with unit variance constraints on
\(\mathbf{\bm{Y}}\bm{\mathbf{w}}\) and \(\mathbf{\bm{X}}\bm{\omega}\). For now, let us define an arbitrary set of
constraint functions \(\{ g_i \colon \mathbb{R}^{M+n+p} \longrightarrow
\mathbb{R}, \, i \in \mathbb{N} \}\), where \(M = \sum_{j =
1}^{\mathbf{\kappa}}{m_j}\) and \(p\) is the length of \(\bm{\gamma}\), and write the full
constrained optimisation problem as:
\begin{equation}\label{eqn:sdf_solution_full}
\begin{split}
\underset{\bm{\gamma},\bm{\mathbf{w}},\left< {\bm{\omega}_j} \right>}{\mathrm{min}} &
\| \mathbf{\bm{Y}} \bm{\mathbf{w}} - F_{\bm{\gamma}}(\left< {\mathbf{\bm{X}}_j \bm{\omega}_j} \right>) \|^2_2 \\
\text{s.t.} \quad &
\{g_i \left( \bm{\gamma},\bm{\mathbf{w}},\left< {\bm{\omega}_j} \right> \right) \geq 0: i \in \mathbb{N}\}
\end{split}
\end{equation}
This paper will henceforth refer to (\ref{eqn:sdf_solution_full}) as the
``Supervised Diffusion Framework'' (SDF) for latent variable modelling.
From a mathematical perspective, SDF can be viewed as a generalisation
over regression analysis on the one hand, and traditional variance
optimisation models on the other (Section \ref{sdf_const} further elaborates
on this point). From a practical perspective, the usefulness of SDF has
at least three facets. First, it address the problem of interpretability
inherent in PCA. Unlike PCA-based diffusion indices, supervised
diffusion indices (SDI) conform to a predefined functional relationship
in which the input variables would usually have an a priori economic
meaning\footnote{\citep{Ahn-2018} makes a similar argument with respect to latent
variables estimated by CCA. Section \ref{cca} shows that CCA is a special
type of a supervised diffusion model.}. Second, a supervised diffusion model (SDM) can be used
to improve the accuracy of empirical measurement for variables which are
not altogether latent yet not very accurately observed. For example, an
economic model describing the relationship between GDP and unemployment
could be used to better approximate both of these aggregates from sets
of alternative proxy measurements. Third, SDMs can produce synthetic
variables whose sole purpose is to track other random processes in
non-trivial ways. For example, an SDM can be used to build an investment
portfolio which tracks the market's expectation for a macroeconomic
index over multiple forecast horizons, or to devise a systematic trading
strategy which can be accurately predicted by a leading signal. The
concluding paragraph of Section \ref{sdf_const} further expands on this use
case by looking at SDMs through the lens of portfolio optimisation.
\subsection{SDF as a Variance Optimisation Framework}
\label{sdf_const}
The SDF formula is inherently a generalised variance optimisation
problem which subsumes a number of popular statistical techniques in
economics and finance as special cases. To demonstrate this, let us
explicitly add an intercept term \(c\) to equation (\ref{eqn:reg_lv_a}) and
rewrite it as:
\begin{equation}
\tilde{y} = c + F_{\bm{\gamma}} \left( \left< {\tilde{x}_j} \right> \right) + e
\end{equation}
The solution for \(c\) in a least squares setting is well documented in
prior literature, but it is replicated below for completeness. Defining
a shorthand \(\hat{\tilde{y}} := F_{\bm{\gamma}} \left( \left< {\tilde{x}_j} \right> \right)\), as well as
the sample counterparts for \(\tilde{y}\) and \(\hat{\tilde{y}}\) as \(\tilde{\mathbf{\bm{y}}}\) and
\(\hat{\tilde{\mathbf{\bm{y}}}}\), respectively, we can rewrite the optimisation problem
(\ref{eqn:sdf_solution}) as:
\begin{equation}\label{eqn:sdf_solution_short}
\min \| \tilde{\mathbf{\bm{y}}} - \left( \1_s c + \hat{\tilde{\mathbf{\bm{y}}}} \right)\|^2_2
= \min \left( \tilde{\mathbf{\bm{y}}}'\tilde{\mathbf{\bm{y}}} + \1_s'\1_s c + \hat{\tilde{\mathbf{\bm{y}}}}' \hat{\tilde{\mathbf{\bm{y}}}} - 2 c \1_s' \tilde{\mathbf{\bm{y}}}
- 2 \tilde{\mathbf{\bm{y}}}'\hat{\tilde{\mathbf{\bm{y}}}} + 2 c \1_s' \hat{\tilde{\mathbf{\bm{y}}}} \right)
\end{equation}
Here, \(s\) denotes the sample size of \(\tilde{\mathbf{\bm{y}}}\) and \(\hat{\tilde{\mathbf{\bm{y}}}}\), and
\(\1_s\) is a column vector of ones of sample length \(s\). Note that
this problem is convex with respect to \(c\), so an unconstrained
solution can be found by setting the partial derivative to zero:
\begin{equation}
\begin{split}
& \frac{\partial}{\partial \, c} \| \tilde{\mathbf{\bm{y}}} - \left( \1_s c + \hat{\tilde{\mathbf{\bm{y}}}} \right)\|^2_2
= 2 \1_s c - 2 \tilde{\mathbf{\bm{y}}} + 2 \hat{\tilde{\mathbf{\bm{y}}}} \\
& \1_s c = \tilde{\mathbf{\bm{y}}} - \hat{\tilde{\mathbf{\bm{y}}}}
\end{split}
\end{equation}
Pre-multiplying both sides by \(\1_s' / s\) yields:
\begin{equation}\label{eqn:sdf_solution_c}
\hat{c} = \frac{\1_s' \tilde{\mathbf{\bm{y}}}}{s} - \frac{\1_s' \hat{\tilde{\mathbf{\bm{y}}}}}{s} = \overline{\tilde{\mathbf{\bm{y}}}} - \overline{\hat{\tilde{\mathbf{\bm{y}}}}}
\end{equation}
where \(\overline{\tilde{\mathbf{\bm{y}}}}\) and \(\overline{\hat{\tilde{\mathbf{\bm{y}}}}}\) denote the sample means of
\(\tilde{\mathbf{\bm{y}}}\) and \(\hat{\tilde{\mathbf{\bm{y}}}}\), respectively. By plugging the unconstrained
solution for \(c\) back into (\ref{eqn:sdf_solution_short}) we arrive at a
simplified version of (\ref{eqn:sdf_solution_short}), namely:
\begin{equation}\label{eqn:sdf_solution_cov}
\min \| \left(\tilde{\mathbf{\bm{y}}} - \1_s\overline{\tilde{\mathbf{\bm{y}}}} \right) - \left(\hat{\tilde{\mathbf{\bm{y}}}} - \1_s\overline{\hat{\tilde{\mathbf{\bm{y}}}}} \right)\|^2_2
= \min \left( \mathbf{\bm{\Sigma}}_{\tilde{y}} - 2 \mathbf{\bm{\Sigma}}_{\tilde{y}\hat{\tilde{y}}} + \mathbf{\bm{\Sigma}}_{\hat{\tilde{y}}} \right)
\end{equation}
Here, \(\mathbf{\bm{\Sigma}}_A\) denotes the sample covariance matrix over random
variable vector \(A\), and \(\mathbf{\bm{\Sigma}}_{AB}\) denotes the sample covariance
matrix from \(A\) to random variable vector \(B\) with the elements of
\(A\) as rows and the elements of \(B\) as columns\footnote{In theory, a constraint could be imposed on \(c\) which makes it
deviate from the solution given by (\ref{eqn:sdf_solution_c}). However, even
with a constraint on \(c\) the simplification provided by
(\ref{eqn:sdf_solution_cov}) remains meaningful with the caveat that something
other than the sample mean would need to be subtracted from \(\tilde{\mathbf{\bm{y}}}\) or
\(\hat{\tilde{\mathbf{\bm{y}}}}\) in equation (\ref{eqn:sdf_solution_cov}). In other words, if a
constraint is imposed on \(c\), \(\mathbf{\bm{\Sigma}}\) will represent a biased sample
variance-covariance estimate for at least one of the variable vectors
involved. See \citep{Lettau-2020} for an example of when this may be
desirable.}. Note that in
this particular case, \(\mathbf{\bm{\Sigma}}_{\tilde{y}}\), \(\mathbf{\bm{\Sigma}}_{\tilde{y}\hat{\tilde{y}}}\) and \(\mathbf{\bm{\Sigma}}_{\hat{\tilde{y}}}\)
all resolve to scalar values because both \(\tilde{y}\) and \(\hat{\tilde{y}}\) are
univariate. Put simply, the SDF objective is nothing more than a
variance minimisation problem which subsumes a few popular models like
PCA and mean-variance portfolio optimisation as special cases. To
illustrate, consider a trivial SDM in which \(F_{\bm{\gamma}} \left( {
\left< {\tilde{x}_j} \right>} \right) = 0\). In this case, (\ref{eqn:sdf_solution_cov}) reduces to:
\begin{equation}\label{eqn:sdm_const}
\begin{split}
& \min \mathbf{\bm{\Sigma}}_{\tilde{y}} \equiv \underset{\bm{\mathbf{w}}}{\min} \, \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} \\
\text{s.t. } &
\{g_i \left( \bm{\mathbf{w}} \right) \geq 0: i \in \mathbb{N}\}
\end{split}
\end{equation}
where \(\{g_i \left( \bm{\mathbf{w}} \right) \geq 0: i \in \mathbb{N}\}\) is an
arbitrary set of constraints on \(\bm{\mathbf{w}}\) (without the inferential
component, \(\bm{\omega}\) and \(\bm{\gamma}\) become redundant). Here, the PCA objective
function\footnote{Technically speaking, the goal of PCA is to maximise variance
rather than to minimise it. However, objective (\ref{eqn:sdm_const}) is still
solved by PCA decomposition with the caveat that variance is minimised
by the last principal component and maximised by the first.} is obtained by applying a single constraint of the form
\(\bm{\mathbf{w}}' \bm{\mathbf{w}} = 1\). Similarly, if \(\mathbf{\bm{Y}}\) is a matrix of historical asset
returns, a minimum variance portfolio is obtained with the constraint \(\1_{n}' \bm{\mathbf{w}} = 1\), and a portfolio on the efficient frontier as in
\citep{Markowitz-1952} is produced with the help of a second constraint
of the form \(\overline{\mathbf{\bm{Y}}} \bm{\mathbf{w}} = r\) for a target return \(r\).
Conceptually, we can look at this result in three ways. First, SDF can
be described as a blend between traditional regression analysis and
variance optimisation models. SDF reduces a traditional regression
problem when all the input variables have one observed proxy each (i.e.,
when \(\bm{\mathbf{w}}\) and all \(\bm{\omega}_j\) are scalars or multiples of standard basis
vectors), and it reduces to a traditional variance minimisation problem
if there is no inferential component \(F_{\bm{\gamma}} \left( \left< {\tilde{x}_j} \right>
\right)\). Second, if we describe SDF as a ``supervised'' diffusion
methodology, we can conversely describe traditional variance
optimisation problems like PCA and mean-variance portfolio optimisation
as unsupervised diffusion models (once again, the ``supervised'' part
comes from the inferential component \(F_{\bm{\gamma}} \left( \left< {\tilde{x}_j} \right>
\right)\) which affects the estimates of the latent variables by
endowing them with inferential properties).
Third, if traditional mean-variance portfolio optimisation can be
described as a trivial SDM, then conversely, non-trivial SDMs have a
place in investment management. Notably, a latent variable constructed
from investible assets under portfolio-style constraints becomes a
viable investment strategy in its own right. This opens up a wide range
of potential applications for SDMs in the investment industry, with
examples including asset allocation, systematic trading, and the
construction of tracker funds. As a starting point, the empirical study
in Section \ref{empirical_results} produces a ``market SDI''which can be
interpreted as a tracker portfolio for growth expectations in the United
States.
\section{Latent Variable Autoregression with Exogenous Inputs}
\label{larx}
The SDF formula can be used to turn any traditional regression model
into a supervised diffusion model (SDM) by allowing the input variables
to be latent. This section provides a tangible example of this process
by deriving an SDM counterpart of the ubiquitous autoregressive model
with exogenous inputs (ARX).
The ARX model is one of the most popular regression models in economics
and finance and in time series analysis more broadly. In ARX models, the
dependent is expressed as a linear function of its own past values and
the (present and) past values of one or more explanatory variables.
Special cases of the ARX model include autoregressive models with no
exogenous inputs (AR), lead-lag regression models with no autoregressive
element, and multiple linear regression models with no lag structure.
To derive a latent variable ARX model (LARX), we can start from the
standard ARX specification and assume that the dependent and all the
explanatory variables are latent. First, let us define the LARX problem
in functional form. Let \(v\) be a ``version iterator'' representing
different versions\footnote{Note that in LARX models, versions need not be synonymous with
time series lags. For example, \(\tilde{y} \equiv Y \bm{\mathbf{w}}\) may represent the
return on an investment index and \(\tilde{x}_i \equiv X_i \bm{\mathbf{w}}\) may be some
other property of the same index, such as market capitalisation as in
the ``size'' factor of \cite{Fama-French-1993}.} of latent variable \(j\), i.e., a variable
identified by LV weight vector \(\bm{\omega}_j\), and let the total number of
versions for variable \(j\) be \(V_j \geq 1\). Let \(K\) denote the
total number of unique exogenous variables \(\tilde{x}_j\) excluding versions,
i.e., the total number of unique weight vectors \(\bm{\omega}_j\). We can refer
to the ``autoregressive'' versions of the dependent variable as
\(\tilde{a}_v\), with the number of autoregressive versions denoted by the
capital \(V_a \geq 0\). A generic formula for the LARX model can then be
written as:
\begin{equation}\label{eqn:larx_a}
\tilde{y} = c + \sum_{v = 1}^{V_a} \phi_{v} \tilde{a}_{v} +
\sum_{v = 1}^{V_1} \beta_{1,v} \tilde{x}_{1,v} +
\sum_{v = 1}^{V_2} \beta_{2,v} \tilde{x}_{2,v} +
\ldots + \sum_{v = 1}^{V_K} \beta_{K,v} \tilde{x}_{K,v} + e
\end{equation}
As a reminder, any input variable in this model becomes non-latent if
the corresponding weight vector only has one non-zero element, and the
entire problem reduces to a traditional ARX model when all the input
variables are non-latent.
Models of this kind are more conveniently solved in matrix form. The
matrix representation for (\ref{eqn:larx_a}) can be derived with the help of
the Kronecker product, traditionally denoted by ``\(\otimes\)'', and the
block-wise Kronecker product denoted by ``\(\odot\)'' as in
\citep{Khatri-Rao-1968}. The autoregressive terms can be written as:
\begin{equation*}
\sum_{v = 1}^{V_a} \phi_{v} \tilde{a}_{v} = A \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \mathbf{r}_a \bm{\phi}
\end{equation*}
where \(A = \begin{bmatrix}A_1, & A_2, & \cdots & A_{V_a}
\end{bmatrix}\) is a horizontal concatenation of the autoregressive
proxy vectors for the dependent, \(\bm{\phi} = \begin{bmatrix}\phi_1, & \phi_2,
& \cdots & \phi_{V_a} \end{bmatrix}'\) is a column vector containing the
respective autoregressive coefficients, and \(\mathbf{r}_a
= \begin{bmatrix}r_{a,1}, & r_{a,1}, & \cdots & r_{a,V_a}
\end{bmatrix}\) is a vector of approximation errors where \(r_{a,v} =
\tilde{a}_v - A_v \bm{\mathbf{w}}\) for \(1 \leq v \leq V_a\). Similarly, for each
individual explanatory variable \(\tilde{x}_j\) we have:
\begin{equation*}
\sum_{v = 1}^{V_j} \beta_{j,v} \tilde{x}_{j,v} =
X_j \left( \bm{\beta}_j \otimes \bm{\omega}_j \right) + \mathbf{r}_j \bm{\beta}_j
\end{equation*}
with \(X_j = \begin{bmatrix}X_{j,1}, & X_{j,2}, & \cdots & X_{j,V_j}
\end{bmatrix}\), \(\bm{\beta}_j = \begin{bmatrix}\beta_{j,1}, & \beta_{j,2}, &
\cdots & \beta_{j,V_j} \end{bmatrix}'\), and \(\mathbf{r}_j\) a (\(1
\times V_j\)) vector of approximation errors for the different versions
of \(\tilde{x}_j\). The explanatory terms can be written out as:
\begin{equation*}
\sum_{j = 1}^{K} \sum_{v = 1}^{V_j} \beta_{j,v} \tilde{x}_{j,v}
= \sum_{j = 1}^{K} X_j \left( \bm{\beta}_j \otimes \bm{\omega}_j \right) +
\sum_{j = 1}^{K} \mathbf{r}_j \bm{\beta}_j
= X \left( {\bm{\beta} \odot \bm{\omega}} \right) + \sum_{j = 1}^{K} \mathbf{r}_j \bm{\beta}_j
\end{equation*}
Here, \(X = \begin{bmatrix} X_1 | X_2 | \cdots | X_K\end{bmatrix}\) is a
row vector with \(K\) blocks corresponding to the individual \(X_j\),
\(\bm{\beta} = \begin{bmatrix} \bm{\beta}_1' | \bm{\beta}_2' | \cdots | \bm{\beta}_K'\end{bmatrix}'\) is
a column vector with \(K\) blocks for the individual \(\bm{\beta}_j\), and
\(\bm{\omega}\) has \(K\) blocks containing the individual \(\bm{\omega}_j\) such that
\(\bm{\omega} = \begin{bmatrix} \bm{\omega}_1' | \bm{\omega}_2' | \cdots | \bm{\omega}_K'\end{bmatrix}'\).
The complete LARX formula is then concisely defined in matrix form as:
\begin{equation} \label{eqn:larx_b}
Y \bm{\mathbf{w}} = c + A \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + X \left( {\bm{\beta} \odot \bm{\omega}} \right) + \epsilon
\end{equation}
where \(Y\), \(A\) and \(X\) comprise the underlying observed variable
space and \(\epsilon\) collects all the error terms.
Given a data sample of length \(s\) represented by observation matrices
\(\mathbf{\bm{Y}}\), \(\mathbf{\bm{A}}\) and \(\mathbf{\bm{X}}\), the least squares optimisation problem for
the LARX model can be defined as
\begin{equation}\label{eqn:clarx_optim_a}
\underset{c,\bm{\phi},\bm{\beta},\bm{\mathbf{w}},\bm{\omega}}{\min} \, \| \mathbf{\bm{Y}} \bm{\mathbf{w}} - \1_s c - \mathbf{\bm{A}} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) - \mathbf{\bm{X}} \left( {\bm{\beta} \odot \bm{\omega}} \right) \|^2_2
\end{equation}
In terms of optimisation constraints, the only requirement for a minimal
LARX implementation is a scaling constraint on \(\bm{\mathbf{w}}\) to avoid the
trivial solution given by \(\bm{\mathbf{w}} = \0\). However, in view of the
investment management use case discussed in the previous section, we can
start by deriving a constrained implementation of the LARX model (call
it CLARX) with full Markowitz-style constraints on the variance and sum
of weights of each latent variable. The ``unconstrained'' LARX model can
then be derived by setting the appropriate Lagrange multipliers to zero.
The variance constraint on the dependent takes the form of \(\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y}
\bm{\mathbf{w}} = \s2_y\). The variance constraint on explanatory variable \(j\)
takes the form of \(\bm{\omega}_j' \mathbf{\bm{\Sigma}}_{X_{j,c_j}} \bm{\omega}_j = \s2_j\), where \(c_j\)
is an arbitrarily chosen version (lag) of \(\tilde{x}_j\). The sum-of-weights
constraints take the form of \(\1_n' \bm{\mathbf{w}} = l_y\) and \(\1_{m_j}' \bm{\omega}_j =
l_j\), for some scalar values \(l_y\) and \(l_j\). The complete
constrained optimisation problem then becomes\footnote{To the author's best knowledge, the closest precedent for this
optimisation problem in prior literature is the EDACCA model defined in
\cite{Xu-2024}. The EDACCA model has a number of important differences,
including a different set of constraints and the use of a single
explanatory variable.}:
\begin{equation}\label{eqn:clarx_optim_b}
\begin{split}
\underset{c,\bm{\phi},\bm{\beta},\bm{\mathbf{w}},\bm{\omega}}{\min} \, & \| \mathbf{\bm{Y}} \bm{\mathbf{w}} - \1_s c - \mathbf{\bm{A}} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) - \mathbf{\bm{X}} \left( {\bm{\beta} \odot \bm{\omega}} \right) \|^2_2 \\[10pt]
\text{s.t. } & \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \s2_y, \quad
\bm{\omega}_j' \mathbf{\bm{\Sigma}}_{X_{j,c_j}} \bm{\omega}_j = \s2_j \text{ for } 1 \leq j \leq K, \\[10pt]
& \1_n' \bm{\mathbf{w}} = l_y, \quad \1_{m_j}' \bm{\omega}_j = l_j \text{ for } 1 \leq j \leq K
\end{split}
\end{equation}
The solution for \(c\) is covered by section \ref{sdf_const}. In this case it
reduces to:
\begin{equation}\label{eqn:clarx_solution_const}
\hat{c} = \overline{\mathbf{\bm{Y}}} \bm{\mathbf{w}} - \overline{\mathbf{\bm{A}}} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) - \overline{\mathbf{\bm{X}}} \left( {\bm{\beta} \odot \bm{\omega}} \right)
\end{equation}
where \(\overline{\mathbf{\bm{Y}}}\), \(\overline{\mathbf{\bm{A}}}\) and \(\overline{\mathbf{\bm{X}}}\) are row
vectors containing the column-wise sample means of \(Y\), \(A\) and
\(X\). Plugging this solution back into (\ref{eqn:clarx_optim_b}) and
expanding, we obtain a simplified optimisation problem of the form:
\begin{align} \label{eqn:clarx_optim_c}
\begin{split}
\underset{\bm{\phi},\bm{\beta},\bm{\mathbf{w}},\bm{\omega}}{\min} \, & \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} + \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \left( {\bm{\beta} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot \bm{\omega}} \right) \\[5pt]
& - 2 \left[ \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YX} \left( {\bm{\beta} \odot \bm{\omega}} \right) - \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right]
\end{split} \\[15pt]
\begin{split}
\text{s.t. } & \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \s2_y, \quad
\bm{\omega}_j' \mathbf{\bm{\Sigma}}_{X_{j,c_j}} \bm{\omega}_j = \s2_j \text{ for } 1 \leq j \leq K, \\[10pt]
& \1_n' \bm{\mathbf{w}} = l_y, \quad \1_{m_j}' \bm{\omega}_j = l_j \text{ for } 1 \leq j \leq K
\end{split} \nonumber
\end{align}
This is a convex problem with equality constraints, which means it can
be solved using the method of Lagrange multipliers (LM). As with the
rest of the model, we can represent the LM terms in matrix form with the
help of the following vectors:
\begin{align*}
\ts2_x & = \begin{bmatrix}
\s2_{x,1}, &
\s2_{x,2}, &
\s2_{x,3}, &
\ldots &
\s2_{x,K}
\end{bmatrix}', \quad
& \bm{\mathbf{\lambda}}_x = \begin{bmatrix}
\lambda_{x,1}, &
\lambda_{x,2}, &
\lambda_{x,3}, &
\ldots &
\lambda_{x,K}
\end{bmatrix}' \\[10pt]
\bm{l}_p & = \begin{bmatrix}
l_{p,1}, &
l_{p,2}, &
l_{p,3}, &
\ldots &
l_{p,K}
\end{bmatrix}', \quad
& \bm{\mathbf{\lambda}}_p = \begin{bmatrix}
\lambda_{p,1}, &
\lambda_{p,2}, &
\lambda_{p,3}, &
\ldots &
\lambda_{p,K}
\end{bmatrix}'
\end{align*}
Here, \(\ts2_x\) is a \(K \times 1\) column vector of variance targets
for the chosen versions of \(\tilde{x}_j\); \(\bm{l}_p\) is a \(K \times 1\) column
vector of sum-of-weights targets for the respective LV weight vectors
\(\bm{\omega}_j\); and \(\bm{\mathbf{\lambda}}_x\) and \(\bm{\mathbf{\lambda}}_p\) are the corresponding vectors of
LM coefficients. The full constraint terms for the Lagrangian function
can then be defined in matrix form with the help of the block-wise
direct sum operator introduced in Section \ref{blockwise_direct_sum}. The
sum-of-weights constraints become:
\begin{equation*}
\left( \bm{\omega}^{\oplus} \right)' \1_{\bm{\omega}} = \bm{l}_p
\end{equation*}
where \(\1_{\bm{\omega}}\) is a column vector of ones with the same block
structure as \(\bm{\omega}\).
For the variance constraints we require an indexer to identify which
version of each \(\tilde{x}_j\) has a variance constraint assigned to it. Let
\(\left< {\bm{\mathbf{u}}_j} \right>\) be a sequence of \(K\) logical vectors, each of size
\(V_j\) (i.e., the number of versions of \(\tilde{x}_j\)) with a value of 1 in
position of \(c_j\), i.e., the version of \(\tilde{x}_j\) that has a variance
constraint, and zeros elsewhere. The column vector \(\bm{\mathbf{u}} := \left[
\left< {\bm{\mathbf{u}}_j} \right> \right]\) is a vertical concatenation of \(\left< {\bm{\mathbf{u}}_j} \right>\) which
has the same size and block structure as the coefficient vector \(\bm{\beta}\).
We also need a block-diagonal matrix \(\mathbf{\bm{\Sigma}}_X^d\) containing the
covariance matrices of the individual \(X_{j,v_j}\) along the diagonal,
but no covariances across different variables or versions. Once again,
this can be done with the help of the block-wise direct sum operator:
\begin{equation*}
\begin{split}
& \mathbf{\bm{\Sigma}}^{d}_X = \left[ \left( \mathbf{\bm{X}} - \overline{\mathbf{\bm{X}}} \right)^{\oplus} \right]'
\left( \mathbf{\bm{X}} - \overline{\mathbf{\bm{X}}} \right)^{\oplus}
= \begin{pmatrix}
\mathbf{\bm{\Sigma}}^{d}_{X_1}, & \0, & \cdots & \0 \\
\0, & \mathbf{\bm{\Sigma}}^{d}_{X_2}, & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0, & \0, & \cdots & \mathbf{\bm{\Sigma}}^{d}_{X_K}
\end{pmatrix} \\
\text{with} \quad &
\mathbf{\bm{\Sigma}}^{d}_{X_j} = \left[ \left( \mathbf{\bm{X}}_j - \overline{\mathbf{\bm{X}}_j} \right)^{\oplus} \right]'
\left( \mathbf{\bm{X}}_j - \overline{\mathbf{\bm{X}}_j} \right)^{\oplus}
= \begin{pmatrix}
\mathbf{\bm{\Sigma}}_{X_{j,1}}, & \0, & \cdots & \0 \\
\0, & \mathbf{\bm{\Sigma}}_{X_{j,2}}, & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0, & \0, & \cdots & \mathbf{\bm{\Sigma}}_{X_{j,V_j}}
\end{pmatrix} \quad \text{for } 1 \leq j \leq K
\end{split}
\end{equation*}
The variance constraints on all \(\tilde{x}_j\) can then be expressed using a
single quadratic form:
\begin{equation*}
\left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)' \right]^{\oplus} \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) = \ts2_x
\end{equation*}
and the full Lagrangian function for (\ref{eqn:clarx_optim_c}) can be written
as:
\begin{equation} \label{eqn:clarx_lagr}
\begin{split}
\mathcal{L}( & \bm{\phi},\bm{\beta},\bm{\mathbf{w}},\bm{\omega},\lambda_y,\lambda_l,\bm{\mathbf{\lambda}}_x,\bm{\mathbf{\lambda}}_p) = \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} +
\left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \left( {\bm{\beta} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot \bm{\omega}} \right) \\[12pt]
& - 2 \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) - 2 \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YX} \left( {\bm{\beta} \odot \bm{\omega}} \right) + 2 \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right)
+ \lambda_y \left( \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} - \s2_y \right) \\[12pt]
& + \lambda_l \left( \1_n' \bm{\mathbf{w}} - l_y \right)
+ \bm{\mathbf{\lambda}}_x' \left\{ \left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)' \right]^{\oplus} \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) - \ts2_x \right\}
+ \bm{\mathbf{\lambda}}_p' \left[ \left( \bm{\omega}^{\oplus} \right)' \1_{\bm{\omega}} - \bm{l}_p \right]
\end{split}
\end{equation}
Note that in expanded form the terms under \(\bm{\mathbf{\lambda}}_x\) and \(\bm{\mathbf{\lambda}}_p\)
resolve to:
\begin{align*}
& \bm{\mathbf{\lambda}}_x' \left\{ \left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)' \right]^{\oplus} \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) - \ts2_x \right\} =
\sum_{j=1}^{K} \lambda_{x,j} \left( \, \bm{\omega}_j' \mathbf{\bm{\Sigma}}_{X_{j,c_j}} \bm{\omega}_j - \s2_{x,j} \right)\\
& \bm{\mathbf{\lambda}}_p' \left[ \left( \bm{\omega}^{\oplus} \right)' \1_{\bm{\omega}} - \bm{l}_p \right] =
\sum_{j = 1}^K \lambda_{p,j} \left(\bm{\omega}_j' \1_{m_j} - l_{p,j} \right)
\end{align*}
where \(m_j\) represents the number of observed variables used to
approximate \(\tilde{x}_j\).
Like both mean-variance portfolio optimisation and linear least-squares
regression, the Lagrangian function specified by (\ref{eqn:clarx_lagr}) is
convex and can be solved by setting various partial derivatives to zero.
The derivations for the individual coefficient vectors and Lagrange
multipliers are somewhat involved and deferred to \ref{derivation_w_o_p_g}
and \ref{derivation_ry_rl}. The solution is given by a fixed point over the
following system of equations:
\begin{subnumcases}{\hspace{-5em}\label{eqn:clarx_solution_full}}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \left[ \mathbf{\bm{\Sigma}}_{Y}^{-1} \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right)
- \rho_l \mathbf{\bm{\Sigma}}_{Y}^{-1} \1_n \right]
\label{eqn:clarx_fp_1} \\[7pt]
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) + \mathbf{\bm{M}}_2 \left( {\bm{\mathbf{\lambda}}_x \odot I_{\bm{\omega}}} \right) \right]^{-1}
\left[ \bm{\mathbf{v}}_3 - \frac{1}{2} \t1od \bm{\mathbf{\lambda}}_p \right]
\label{eqn:clarx_fp_2} \\[7pt]
\bm{\hat{\phi}} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1}
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} - \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right]
\label{eqn:clarx_fp_3} \\[7pt]
\bm{\hat{\beta}} = \left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \right]^{-1}
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) \right]
\label{eqn:clarx_fp_4} \\[7pt]
\rho_{y} = \frac{ \left( n \bm{\mathbf{w}} - l_y \1_n \right)' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) }
{n \s2_y - l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}}
\label{eqn:clarx_fp_5} \\[7pt]
\rho_l = \frac{1}{n} \left[ \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) -
\rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} \right]
\label{eqn:clarx_fp_6} \\[7pt]
\bm{\mathbf{\lambda}}_x = \left[ \mathbf{\bm{M}}_1 \bm{\mathbf{\Theta}} - \mathbf{\bm{L}} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \right]^{-1}
\left( \bm{\omega}^{\oplus} \mathbf{\bm{M}}_1 - \t1od \mathbf{\bm{L}} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right)
\label{eqn:clarx_fp_7} \\[7pt]
\bm{\mathbf{\lambda}}_p = 2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)'
\left[ \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x \right]
\label{eqn:clarx_fp_8}
\end{subnumcases}
with the shorthand notations:
\begin{flalign*}
& \1_{\bm{\omega}} \text{ a column vector of ones with the same length and block structure as } \bm{\omega} & \\
& \bm{\mathbf{v}}_1 = \left[
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right)
\right] \bm{\mathbf{w}} & \\
& \bm{\mathbf{v}}_2 = \left[\mathbf{\bm{\Sigma}}_{YX} - \left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AX} \right] \left( {\bm{\beta} \odot \bm{\omega}} \right) & \\
& \bm{\mathbf{v}}_3 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}} & \\
& \bm{\mathbf{v}}_4 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega} & \\
& \bm{\mathbf{\Theta}} = \operatorname{diag}(\ts2_x) & \\
& \mathbf{\bm{L}} = \operatorname{diag}(\bm{l}_p) & \\
& \mathbf{\bm{M}}_1 = \left( \t1od \right)' \t1od & \\
& \mathbf{\bm{M}}_2 = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) & \\
& \rho_{y} = 1 + \lambda_{y}, \quad \rho_l = \frac{\lambda_l}{2} &
\end{flalign*}
In the above expressions, \(I_a\) represents an identity matrix of
scalar size \(a\), while \(I_{\mathbf{\bm{a}}}\) represents an identity matrix with
the same size and block structure as that of the vector \(\mathbf{\bm{a}}\) (bold).
Here, it must be noted that the derivation steps required to arrive at
(\ref{eqn:clarx_solution_full}) follow directly from the properties of the
block-wise direct sum operator and the block-wise Kronecker product for
vectors introduced in Section \ref{blockwise_direct_sum}, \ref{bds_vs_mmul} and
\ref{block_kron_factorisation}. A thorough discussion of the relevance of
these results for the broader field of linear algebra is beyond the
scope of this paper; however, the derivations found in
\ref{derivation_w_o_p_g} and \ref{derivation_ry_rl} generalise rather easily to
other use cases in matrix calculus, including a larger class of
Lagrangian optimisation problems in which an arbitrary sequence of
coefficient vectors needs to be estimated in the presence of
interactions and case-by-case constraints.
In practice, the solution to (\ref{eqn:clarx_solution_full}) can be obtained
by fixed point iteration, which has sub-linear computational complexity
if the underlying regression problem is well-specified. Initial guesses
are required for \(\bm{\mathbf{w}}\), \(\bm{\omega}\) and \(\bm{\phi}\). Equations can be estimated
in the same order as shown above. The first iteration can start at
(\ref{eqn:clarx_fp_4}). Four matrix inversions need to be recalculated
at each iteration step (namely,
(\ref{eqn:clarx_fp_2})-(\ref{eqn:clarx_fp_4}) and
(\ref{eqn:clarx_fp_7})) because the corresponding matrices are derived
from the estimated coefficient vectors; however, equations
(\ref{eqn:clarx_fp_3}) and (\ref{eqn:clarx_fp_4}) can be reformulated in
terms of the Moore-Penrose inverses (\citep{Penrose-1955,Bjerhammar-1951,Moore-1920}) of the matrices \((\mathbf{\bm{A}} - \overline{\mathbf{\bm{A}}})
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)\) and \((\mathbf{\bm{X}} - \overline{\mathbf{\bm{X}}}) \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)\), respectively, which makes them
solvable by SVD.
\subsection{The Unconstrained Case and the Intuitive Interpretation}
\label{larx_ols_interpretation}
The constrained solution to the LARX model given by
(\ref{eqn:clarx_solution_full}) may look somewhat complex at first. However,
it simply says that the objective function is minimised when each of the
four coefficient vectors (\(\bm{\mathbf{w}}\), \(\bm{\omega}\), \(\bm{\phi}\), \(\bm{\beta}\)) comes as close
to a least-squares solution to (\ref{eqn:larx_b}) as possible without
violating a portfolio constraint. We can see this more clearly by
removing the portfolio-style constraints from the model and only keeping
the one constraint required for model convergence\footnote{A scaling constraint has to be imposed on the dependent weight
vector \(\bm{\mathbf{w}}\) to prevent the problem from converging to the trivial
solution given by \(\bm{\mathbf{w}} = \0\). Several possibilities exist for such a
constraint, the most popular ones being a length constraint \(\bm{\mathbf{w}}' \bm{\mathbf{w}} =
k\) and a unit variance constraint \(\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \s2\), where \(k\)
and \(\s2\) can be arbitrary scalars but are usually set to unity. A
unit length constraint is used in PCA and PLS; a unit variance
constraint is used in CCA. This paper chooses a unit variance constraint
for several reasons, including the ability to draw a parallel between
LARX and CCA, the usefulness of a variance constraint in portfolio
management settings, and the better interpretability of the concept of
variance in finance and economics more generally.}, namely, that
on the variance of \(\tilde{y}\). Setting the corresponding Lagrange
multipliers (\(\lambda_l\), \(\bm{\mathbf{\lambda}}_x\) and \(\bm{\mathbf{\lambda}}_p\)) to zero, we arrive at
the following (unconstrained) solution:
\begin{subnumcases}{\hspace{-5em}\label{eqn:larx_solution_full}}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right)
\label{eqn:larx_fp_1} \\[10pt]
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \right]^{-1} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}}
\label{eqn:larx_fp_2} \\[10pt]
\bm{\hat{\phi}} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1}
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} - \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right]
\label{eqn:larx_fp_3} \\[10pt]
\bm{\hat{\beta}} = \left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \right]^{-1}
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}}
\label{eqn:larx_fp_4} \\[10pt]
\rho_{y} = \frac{ \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) }{\s2_y}
\label{eqn:larx_fp_5}
\end{subnumcases}
with \(\bm{\mathbf{v}}_1\), \(\bm{\mathbf{v}}_2\) and \(\rho_{y}\) all defined as before. In this
model, each of the four coefficient vectors \(\bm{\mathbf{\hat{w}}}\), \(\bm{\hat{\omega}}\), \(\bm{\hat{\phi}}\),
\(\bm{\hat{\beta}}\) can be expressed as an OLS solution to a linear regression
problem of the form \(a = c + B \bm{\gamma}_i + e\), where both \(a\) and \(B\)
are derived using the other three vectors. For example, the coefficient
\(\bm{\hat{\beta}}\) solves:
\begin{equation*}
\begin{split}
& \tilde{r}_a = c + \widetilde{X} \bm{\beta} + e \\
\text{with } & \widetilde{X} := \begin{bmatrix} \tilde{x}_{1,1}, & \tilde{x}_{1,2}, & \ldots & \tilde{x}_{1, V_1}, & \ldots & \tilde{x}_{K, V_k} \end{bmatrix}
\quad \text{and} \quad \tilde{r}_a := \tilde{y} - \sum_{v = 1}^{V_a} \phi_v \tilde{a}_v
\end{split}
\end{equation*}
Here, \(\widetilde{X}\) is a vector of all explanatory variables in latent form,
which can be estimated by \(X \left( {I_{\bm{\beta}} \odot \bm{\hat{\omega}}} \right)\) if \(\bm{\hat{\omega}}\) is
known. Similarly, the variable \(\tilde{r}_a\) is the autoregressive residual
of \(\tilde{y}\) which can be estimated by \(Y \bm{\mathbf{\hat{w}}} - A \left( {\bm{\hat{\phi}} \otimes \bm{\mathbf{\hat{w}}}} \right)\). In other words, \(\bm{\hat{\beta}}\) solves a regression problem of the form \(a
= c + B \bm{\gamma}_i + e\) where \(a := Y \bm{\mathbf{\hat{w}}} - A \left( {\bm{\hat{\phi}} \otimes \bm{\mathbf{\hat{w}}}} \right)\) and
\(B := X \left( {I_{\bm{\beta}} \odot \bm{\hat{\omega}}} \right)\).
Similar transformations can be applied to derive the other coefficient
vectors. The corresponding formulas for \(a\) and \(B\) are provided in
Table \ref{table:larx_ols_matrices} below:
\begin{table}[ht]
\centering
\caption{OLS Inputs for LARX Coefficients}
\label{table:larx_ols_matrices}
\begin{threeparttable}
\begin{tabular}{CCC}
\toprule
\text{Coefficient vector } \bm{\gamma}_i & \text{Dependent variable } a & \text{Explanatory vector } B \\
\midrule
\bm{\mathbf{w}} & X \left( {\bm{\hat{\beta}} \odot \bm{\hat{\omega}}} \right) & Y - A \left( {\bm{\hat{\phi}} \otimes I_n} \right) \\
\bm{\phi} & Y \bm{\mathbf{\hat{w}}} - X \left( {\bm{\hat{\beta}} \odot \bm{\hat{\omega}}} \right) & A \left( { I_{V_a} \otimes \bm{\mathbf{\hat{w}}} } \right) \\
\bm{\omega} & Y \bm{\mathbf{\hat{w}}} - A \left( {\bm{\hat{\phi}} \otimes \bm{\mathbf{\hat{w}}}} \right) & X \left( { \bm{\hat{\beta}} \odot I_{\bm{\omega}} } \right) \\
\bm{\beta} & Y \bm{\mathbf{\hat{w}}} - A \left( {\bm{\hat{\phi}} \otimes \bm{\mathbf{\hat{w}}}} \right) & X \left( { I_{\bm{\beta}} \odot \bm{\hat{\omega}} } \right) \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\item \footnotesize Note: Formulas for a dependent variable \(a\) and
an explanatory vector \(B\) which conditionally produce a standard
linear regression model \(a = c + B \bm{\gamma}_i + e\) for which the given
coefficient vector \(\bm{\hat{\gamma}}_i\) provides the OLS solution.
\end{tablenotes}
\end{threeparttable}
\end{table}
This observation is useful in two respects. First, it provides an
alternative interpretation to the LARX model and its regression
parameters. Second, it creates a convenient shortcut for the purposes of
heuristic feature selection. If a coefficient vector can be estimated by
OLS, then so can its standard error, with the caveat that both are
obtained conditionally. Conditional standard errors fall short of
providing a full picture of model uncertainty, but they can shed some
light on the relative statistical significance of the individual
elements in each coefficient vector. This should allow the researcher to
perform heuristic feature selection for the LARX model in the same
manner as one would do with a standard least-squares regression.
\section{Special Cases of the LARX Model}
\label{larx_special_cases}
Like ARX, the LARX methodology subsumes several regression models as
special cases, and some of these special cases have interesting
mathematical properties in their own right. The main special cases of
the LARX model are presented in Table \ref{table:larx_summary}. Three of them
are examined more closely in this section.
\afterpage{
\begin{landscape}
\begin{table}[p]
\caption{Special Cases of the LARX model}
\label{table:larx_summary}
\begin{adjustbox}{width=\textwidth}
\begin{threeparttable}
\begin{tabular}{p{0.2\textwidth} p{0.55\textwidth} p{0.5\textwidth} p{0.45\textwidth}}
&
\thead{LARX: Latent Variable Autoregression with Exogenous Inputs} &
\thead{3LR: Latent Variable Lead-Lag Regression} &
\thead{LAR: Latent Variable Autoregression} \\
\toprule
Description &
The most generic implementation -- see Section \ref{larx}. &
No autoregressive term. An SDM counterpart to a multivariate
lead-lag regression. &
No exogenous inputs. An SDM equivalent to the traditional
autoregressive (AR) model. \\
\midrule
Regression formula &
\( \tilde{y} = c + \sum_{v = 1}^{V_a} \phi_{v} \tilde{a}_{v} +
\sum_{j = 1}^{K} \sum_{v = 1}^{V_j} \beta_{j,v} \tilde{x}_{j,v} + e \) &
\( \tilde{y} = c + \sum_{j = 1}^{K} \sum_{v = 1}^{V_j} \beta_{j,v} \tilde{x}_{j,v} + e \) &
\( \tilde{y} = c + \sum_{v = 1}^{V_a} \phi_{v} \tilde{a}_{v} + e \) \\
Formula in matrix form &
\(Y \bm{\mathbf{w}} = c + A \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + X \left( {\bm{\beta} \odot \bm{\omega}} \right) + \epsilon \) &
\(Y \bm{\mathbf{w}} = c + X \left( {\bm{\beta} \odot \bm{\omega}} \right) + \epsilon \) &
\(Y \bm{\mathbf{w}} = c + A \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \epsilon \) \\
\midrule
Solution (optional constraints in \color{red}{red}) &
\( \left\{ \begin{array}{l}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right)
\color{red}{- \frac{\rho_l }{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \1_n} \\
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \color{red}{+ \mathbf{\bm{M}}_2 \left( {\bm{\mathbf{\lambda}}_x \odot I_{\bm{\omega}}} \right)} \color{black}{} \right]^{-1}
\left[ \bm{\mathbf{v}}_3 \color{red}{- \frac{1}{2} \t1od \bm{\mathbf{\lambda}}_p} \color{black}{} \right] \\
\bm{\hat{\phi}} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1}
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} - \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right] \\
\bm{\hat{\beta}} = \left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \right]^{-1}
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) \right] \\
\rho_{y} = \frac{ \left( n \bm{\mathbf{w}} \color{red}{ - l_y \1_n } \color{black}{} \right)' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) }
{n \s2_y \color{red}{- l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}}} \\
\color{red}{\rho_l = \frac{1}{n} \left[ \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) -
\rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} \right] }\\
\color{red}{\bm{\mathbf{\lambda}}_x = \left[ \mathbf{\bm{M}}_1 \bm{\mathbf{\Theta}} - \mathbf{\bm{L}} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \right]^{-1}
\left( \bm{\omega}^{\oplus} \mathbf{\bm{M}}_1 - \t1od \mathbf{\bm{L}} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) } \\
\color{red}{\bm{\mathbf{\lambda}}_p = 2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)'
\left[ \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x \right] }
\end{array} \right.
\) &
\( \left\{ \begin{array}{l}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \left( {\bm{\beta} \odot \bm{\omega}} \right)
\color{red}{- \frac{\rho_l }{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \1_n} \\
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \color{red}{+ \mathbf{\bm{M}}_2 \left( {\bm{\mathbf{\lambda}}_x \odot I_{\bm{\omega}}} \right)} \color{black}{} \right]^{-1}
\left[ \bm{\mathbf{v}}_3 \color{red}{- \frac{1}{2} \t1od \bm{\mathbf{\lambda}}_p} \color{black}{} \right] \\
\bm{\hat{\beta}} = \left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \right]^{-1} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} \\
\rho_{y} = \frac{ \left( n \bm{\mathbf{w}} \color{red}{ - l_y \1_n } \color{black}{} \right)' \mathbf{\bm{\Sigma}}_{YX} \left( {\bm{\beta} \odot \bm{\omega}} \right) }
{n \s2_y \color{red}{- l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}}} \\
\color{red}{\rho_l = \frac{1}{n} \left[ \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) -
\rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} \right] }\\
\color{red}{\bm{\mathbf{\lambda}}_x = \left[ \mathbf{\bm{M}}_1 \bm{\mathbf{\Theta}} - \mathbf{\bm{L}} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \right]^{-1}
\left( \bm{\omega}^{\oplus} \mathbf{\bm{M}}_1 - \t1od \mathbf{\bm{L}} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) } \\
\color{red}{\bm{\mathbf{\lambda}}_p = 2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)'
\left[ \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x \right] }
\end{array} \right.
\) &
\( \left\{ \begin{array}{l}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_Y^{-1} \bm{\mathbf{v}}_1 \color{red}{- \frac{\rho_l}{\rho_{y}} \1_n } \\
\bm{\hat{\phi}} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} \\
\rho_{y} = \frac{ n \bm{\mathbf{w}}' \bm{\mathbf{v}}_1 \color{red}{- l_y \1_n' \bm{\mathbf{v}}_1} }
{n \s2_y \color{red}{- l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}} } \\
\color{red}{\rho_l = \frac{1}{n} \left[\1_n' \bm{\mathbf{v}}_1 - \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} \right]}
\end{array} \right.
\) \\
\midrule
Shorthand notations &
\( \begin{array}{l}
\1_{\bm{\omega}} \text{ a column vector of ones with the same length and block structure as } \bm{\omega} \\
\bm{\mathbf{v}}_1 = \left[
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right)
\right] \bm{\mathbf{w}} \\
\bm{\mathbf{v}}_2 = \left[\mathbf{\bm{\Sigma}}_{YX} - \left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AX} \right] \left( {\bm{\beta} \odot \bm{\omega}} \right) \\
\bm{\mathbf{v}}_3 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}} \\
\bm{\mathbf{v}}_4 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega} \\
\bm{\mathbf{\Theta}} = \operatorname{diag}(\ts2_x), \, \mathbf{\bm{L}} = \operatorname{diag}(\bm{l}_p), \,
\mathbf{\bm{M}}_1 = \left( \t1od \right)' \t1od, \, \mathbf{\bm{M}}_2 = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \\
\end{array}
\) &
\( \begin{array}{l}
\bm{\mathbf{v}}_3 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} \\
\bm{\mathbf{v}}_4 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega} \\
\bm{\mathbf{\Theta}} = \operatorname{diag}(\ts2_x), \quad \mathbf{\bm{L}} = \operatorname{diag}(\bm{l}_p) \\
\mathbf{\bm{M}}_1 = \left( \t1od \right)' \t1od, \quad \mathbf{\bm{M}}_2 = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \\
\rho_{y} = 1 + \lambda_{y}, \quad \rho_l = \frac{\lambda_l}{2}
\end{array}
\) &
\( \bm{\mathbf{v}}_1 \) as in the LARX model \\
\bottomrule
&
\rule{0pt}{40pt}\thead{LSR: Latent Shock Regression} &
\thead{LVMR/CCA: Latent Variable Multiple Regression} &
\thead{LAR(1)/CAA: First-order LAR Model} \\
\toprule
Description &
A lead-lag regression between an observed dependent and a single
latent explanatory. A parsimonious alternative to a traditional
multivariate lead-lag regression model. See Section \ref{lsr}. &
Linear regression between a latent dependent and one or more (latent
or non-latent) explanatory variables. Each latent variable enters with
one version (no duplicate weight vectors). See Section \ref{cca}. &
A latent variable autoregressive model of order LAR(1). Solvable by
Canonical Autocorrelation Analysis (CAA) when \(\mathbf{\bm{\Sigma}}_Y = \mathbf{\bm{\Sigma}}_A \)
(covariance stationary proxy vector). See Section \ref{caa}.
\\
\midrule
Regression formula &
\( y_t = c + \sum_{\tau = 1}^{F} \beta_{\tau} L^{\tau} \tilde{x}_t + e_t \)
\hspace{1.5em} where \(L\) denotes the lag operator &
\( \tilde{y} = c + \sum_{j = 1}^{K} \beta_{j} \tilde{x}_{j} + e \) &
\( \tilde{y}_t = c + \phi \tilde{y}_{t-1} + e_t \) \\
Formula in matrix form &
\( y_t = c + X \left( {\bm{\beta} \otimes \bm{\omega}} \right) + \epsilon_t \) \hspace{3em}
with \( X =
\begin{bmatrix} X_{t-1}, & X_{t-2}, & \ldots & X_{t-F} \end{bmatrix}
\) &
\( Y \bm{\mathbf{w}} = c + X \bm{\omega} + \epsilon \) \hspace{0.5em} where \(X\), \(\bm{\omega}\)
are vector concatenations of \(\left< {X_j} \right>\), \(\left< {\bm{\omega}_j} \right>\) &
\( Y \bm{\mathbf{w}} = c + \phi A \bm{\mathbf{w}} + \epsilon \) \\
\midrule
Solution &
\( \left\{ \begin{array}{l}
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \otimes I_m} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \otimes I_m} \right) \right]^{-1} \left( {\bm{\beta} \otimes I_m} \right)' \mathbf{\bm{\Sigma}}_{Xy} \\
\bm{\hat{\beta}} = \left[ \left( {I_F \otimes \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_F \otimes \bm{\omega}} \right) \right]^{-1} \left( {I_F \otimes \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{Xy}
\end{array} \right.
\) &
\makecell[l]{
\vspace{0.2em}
\( \left\{ \begin{array}{l}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \bm{\omega} \\
\bm{\hat{\omega}} = \mathbf{\bm{\Sigma}}_{X}^{-1} \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} \\
\rho_{y} = \frac{ \bm{\mathbf{w}}' \, \mathbf{\bm{\Sigma}}_{YX} \bm{\omega} }{\s2_y}
\end{array} \right. \vspace{0.2em}
\) \\or: \hspace{1em} \(
\rho_{y} \bm{\mathbf{\hat{w}}} = \mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \mathbf{\bm{\Sigma}}_{X}^{-1} \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}}
\)
\vspace{0.2em}
} &
\makecell[l]{
\vspace{0.2em}
\( \left\{ \begin{array}{l}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_Y^{-1} \phi
\left[\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} - \phi \mathbf{\bm{\Sigma}}_{A} \right] \bm{\mathbf{w}} \\
\hat{\phi} = \frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}}}{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{A} \bm{\mathbf{w}}} \\
\rho_{y} = \frac{\phi }{\s2_y} \bm{\mathbf{w}}' \left[\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} - \phi \mathbf{\bm{\Sigma}}_{A}
\right] \bm{\mathbf{w}}
\end{array} \right. \vspace{0.2em}
\) \\or: \hspace{1em} \(
\phi \bm{\mathbf{\hat{w}}} = \frac{1}{2} \mathbf{\bm{\Sigma}}_Y^{-1} \left( { \mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA}} \right) \bm{\mathbf{w}}
\) \hspace{1em} if \hspace{1em} \(\mathbf{\bm{\Sigma}}_Y = \mathbf{\bm{\Sigma}}_A \)
\vspace{0.2em}
}\\
\bottomrule
\end{tabular}
\end{threeparttable}
\end{adjustbox}
\end{table}
\end{landscape}
}
\subsection{LSR: A Parsimonious Alternative to a Traditional Lead-Lag Regression}
\label{lsr}
Just like the traditional ARX model, the LARX model is estimated by the
method of least squares applied to a collection of observed variables.
The key difference lies in how the regression coefficients are mapped to
the observed data. ARX models always assign a unique response
coefficient to each variable, whereas LARX models allow some
coefficients to enter the equation more than once as a side-effect of
the interaction between coefficient vectors and latent variable weight
vectors. As a result, LARX models will often be more parsimonious than
traditional regression models applied to the same dataset.
This difference is best exemplified by a class of models in which an
observed dependent variable \(y\) enters as a function of \(V\) versions
of a single latent explanatory variable \(\tilde{x}\) estimated using a \(\left( {1 \times m} \right)\) vector of observed proxy variables, \(X\). In a
time series context, the version iterator \(v\) is more intuitively
described as a lag iterator \(\tau\) relative to a time subscript \(t\),
so the corresponding regression problem becomes:
\begin{equation}\label{eqn:lsr_a}
\begin{split}
& y_t = c + \sum_{\tau = 1}^{F} \beta_{\tau} \tilde{x}_{t-\tau} + e_t
\quad \text{ or, equivalently:} \\
& y_t = c + \sum_{\tau = 1}^{F} \beta_{\tau} L^{\tau} \tilde{x}_{t} + e_t
\quad \text{where } L \text{ is the lag operator}
\end{split}
\end{equation}
As before, the latent variable \(\tilde{x}\) is approximated as \(\tilde{x}_t = X_t
\bm{\omega} - r_t\) using a time-invariant weight vector \(\bm{\omega}\) with dimensions
\(\left( {m \times 1} \right)\) and a measurement error term \(r_t\). In matrix
form, this reduces to:
\begin{equation} \label{eqn:lsr_b}
y_t = c + X \left( {\bm{\beta} \otimes \bm{\omega}} \right) + \epsilon_t
\end{equation}
where \(X := \begin{bmatrix} X_t, & X_{t-1}, & X_{t-2}, & \ldots &
X_{t-F} \end{bmatrix}\) is a \(\left( {1 \times mF} \right)\) vector of observed
variables which contains all the lag powers of the proxy vector \(X_t\),
and \(\bm{\beta}\) has dimensions \(F \times 1\). This paper will henceforth
refer to (\ref{eqn:lsr_b}) a Latent Shock Regression (LSR).
The unconstrained solution for the LSR model is given by:
\begin{subnumcases}{\hspace{-15em}\label{eqn:lsr_solution_full}}
\bm{\hat{\omega}} = \left[ \left( {\bm{\beta} \otimes I_m} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \otimes I_m} \right) \right]^{-1} \left( {\bm{\beta} \otimes I_m} \right)' \mathbf{\bm{\Sigma}}_{Xy}
\label{eqn:lsr_fp_1} \\[10pt]
\bm{\hat{\beta}} = \left[ \left( {I_F \otimes \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_F \otimes \bm{\omega}} \right) \right]^{-1} \left( {I_F \otimes \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{Xy}
\label{eqn:lsr_fp_2}
\end{subnumcases}
\begin{equation} \label{eqn:lsr_solution_const}
\hspace{-26em} \hat{c} = \overline{\mathbf{\bm{y}}} - \overline{\mathbf{\bm{X}}} \left( {\bm{\beta} \otimes \bm{\omega}} \right)
\end{equation}
One can look at the LSR model in two ways. On the one hand, it is an SDM
with \(F\) lags of a single latent factor. On the other hand, it is a
multivariate lead-lag regression with \(F\) lags of \(m\) observed
explanatory variables. The latter interpretation makes LSR directly
comparable to a traditional lead-lag regression model, which reveals one
key difference: In a standard model, each lag of each regressor would be
mapped to a unique regression coefficient, resulting in \(mF\) slope
parameters. Meanwhile, the LSR model only requires \(F\) coefficients
for the vector \(\bm{\beta}\) and \(m\) coefficients for the vector \(\bm{\omega}\),
resulting in \(F + m\) slope parameters. This difference can lead to a
substantial divergence in model complexity when both \(F\) and \(m\) are
large.
The gain in parsimony achieved by LSR models comes from one key
simplifying assumption. In LSR models, all the observed explanatory
variables are treated as facets (``proxies'') of a single unobserved
shock process \(\tilde{x}\). As a result, all observed explanatory variables
are assumed to affect the dependent with a shared lag profile: the lag
profile of \(\tilde{x}\) given by the \(\left( {F \times 1} \right)\) slope vector
\(\bm{\beta}\). Meanwhile, the weights of the observed variables in \(\tilde{x}\) are
time-invariant and given by the \(\left( {m \times 1} \right)\) vector \(\bm{\omega}\).
What would be an \(\left( {mF \times 1} \right)\) coefficient vector in a
traditional lead-lag regression is then replaced by an \(\left( {mF \times
1} \right)\) Kronecker product \(\left( {\bm{\beta} \otimes \bm{\omega}} \right)\), and the slope
coefficient for lag \(\tau\) of predictor \(i\) is given by the product
of the \(i\)'th element of \(\bm{\omega}\) and the \(\tau\)'th element of \(\bm{\beta}\).
This limiting assumption can lead to a reduction in the quality of
in-sample regression fit for LSR models compared to traditional lead-lag
regression models. On the other hand, it can also reduce the risk of
overfitting when the assumption about homogenous lag profiles is
realistic. Furthermore, the structural constraint of the LSR model can
be relaxed by slicing the proxy vector \(X_t\) into blocks and
constructing multiple latent shock processes \(\tilde{x}_j\) using a
block-wise Kronecker product \(\left( {\bm{\beta} \odot \bm{\omega}} \right)\). Multivariate LSR (MLSR)
specifications can provide additional flexibility in terms of the number
of unique lag profiles of transmission (and, accordingly, the number of
required regression parameters), serving as a middle ground between pure
(univariate) LSR and traditional lead-lag regression models\footnote{Given a \(\left( {1 \times m} \right)\) vector of observed regressors
\(X_t\), a multivariate LSR model would reduce to a univariate LSR model
when \(X_t\) has one block and to a traditional lead-lag regression
model when \(X_t\) has \(m\) blocks.}.
\subsection{LVMR: A Regression Model Solvable by Canonical Correlation Analysis}
\label{cca}
Section \ref{sdf} argues that the LARX model subsumes Canonical Correlation
Analysis (CCA). To demonstrate this, let us consider a subcategory of
LARX models in which each latent variable enters the equation with
exactly one version (no duplicate weight vectors). We can broadly
categorise these models as Latent Variable Multiple Regression (LVMR)
models because of their resemblance to the class of models examined in
\citep{Burnham-MacGregor-1996}. LVMR models have two simplifying
features compared to the full LARX specification: First, there are no
autoregressive lags and hence no vector \(\bm{\phi}\). Second, because each
latent explanatory variable only appears once, the entire vector \(\bm{\beta}\)
becomes redundant\footnote{Formally, if each \(\tilde{x}_j\) only has one version, then the
corresponding block \(\bm{\beta}_j\) of coefficient vector \(\bm{\beta}\) only has one
element. Each term \(\mathbf{\bm{X}}_j \left( { \bm{\beta}_j \otimes \bm{\omega}_j} \right)\) then reduces to
\(X_j \bm{\omega}_j \beta_j\) where \(\beta_j\) is a scalar. If there is no
scaling constraint imposed on \(\bm{\omega}_j\), \(\beta_j\) becomes redundant as
it can be ``absorbed'' into the solution for \(\bm{\omega}_j\). If there is a
scaling constraint on \(\bm{\omega}_j\), then it can be imposed by multiplying
the unconstrained solution by the appropriate scalar value after the
fact and \(\beta_j\) would represent that scalar value.}. With these simplifications, the regression
formula reduces to:
\begin{equation}\label{eqn:lvmr_b}
Y \bm{\mathbf{w}} = c + X \bm{\omega} + \epsilon
\end{equation}
The solution to this problem is given by:
\begin{subnumcases}{\hspace{-27em}\label{eqn:lvmr_solution}}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \bm{\omega}
\label{eqn:lvmr_fp_1} \\[10pt]
\bm{\hat{\omega}} = \mathbf{\bm{\Sigma}}_{X}^{-1} \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}}
\label{eqn:lvmr_fp_2} \\[10pt]
\rho_{y} = \frac{ \bm{\mathbf{w}}' \, \mathbf{\bm{\Sigma}}_{YX} \bm{\omega} }{\s2_y}
\label{eqn:lvmr_fp_3}
\end{subnumcases}
\begin{flalign} \label{eqn:lvmr_solution_const}
& \hat{c} = \overline{\mathbf{\bm{Y}}} \bm{\mathbf{w}} - \overline{\mathbf{\bm{X}}} \bm{\omega} &
\end{flalign}
LVMR models sit at the cusp between LARX, CCA and traditional multiple
linear regression models. Firstly, equation (\ref{eqn:lvmr_fp_2})
represents the least squares solution for a multiple linear regression
between \(\tilde{y}\) and the individual observed variables in \(X\)
(semantically, the vector \(\bm{\omega}\) could just as well be called \(\bm{\beta}\),
and whether any blocks of \(\bm{\omega}\) represent LV weight vectors for some
latent variables \(\tilde{x}_j\) is a question of interpretation only). As a
corollary, an LVMR model reduces to a standard multiple regression model
when \(\tilde{y}\) is non-latent. Secondly, solving (\ref{eqn:lvmr_solution}) is
equivalent to finding the canonical variates for \(Y\) and \(X\), as
previously observed by \citep{Dong-2018}. If we substitute
(\ref{eqn:lvmr_fp_2}) into (\ref{eqn:lvmr_fp_1}) and set \(\s2_y = 1\),
we obtain:
\begin{equation} \label{eqn:cca_fp}
\bm{\mathbf{\hat{w}}} = \frac{\mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \mathbf{\bm{\Sigma}}_{X}^{-1} \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}}}
{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y}^{-1} \mathbf{\bm{\Sigma}}_{YX} \mathbf{\bm{\Sigma}}_{X}^{-1} \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}}}
\end{equation}
which is the mathematical formula for CCA.
\subsection{Latent Variable Autoregression (LAR) and Canonical Autocorrelation Analysis (CAA)}
\label{caa}
Another interesting special case of the LARX model is a latent variable
autoregressive model with no exogenous inputs (LAR), which can be viewed
as an SDM counterpart of the traditional autoregressive model. This
section briefly discusses the LAR model and its own special case: a
latent autoregressive model of order LAR(1). With a few simplifying
assumptions, the LAR(1) model reduces to a new type of
eigendecomposition problem which is similar to PCA and CCA.
LAR models can be defined by stripping away the exogenous term from
equation (\ref{eqn:larx_b}), which results in the following formula:
\begin{equation} \label{eqn:clar}
Y \bm{\mathbf{w}} = c + A \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) + \epsilon
\end{equation}
The fixed point solution to this problem is given by:
\begin{subnumcases}{\hspace{-14em}\label{eqn:clar_solution_full}}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_Y^{-1} \left( { \bm{\mathbf{v}}_1 - \rho_l \1_n } \right)
\label{eqn:clar_fp_1} \\[10pt]
\bm{\hat{\phi}} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}}
\label{eqn:clar_fp_2} \\[10pt]
\rho_{y} = \frac{ \left( n \bm{\mathbf{w}} - l_y \1_n \right)' \bm{\mathbf{v}}_1 }
{n \s2_y - l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}}
\label{eqn:clar_fp_3} \\[10pt]
\rho_l = \frac{1}{n} \left[
\1_n' \bm{\mathbf{v}}_1 - \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}
\right] \label{eqn:clar_fp_4}
\end{subnumcases}
\begin{flalign*}
& \text{with } \bm{\mathbf{v}}_1 = \left[
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right)
\right] \bm{\mathbf{w}} &
\end{flalign*}
Similar methodologies have been examined by papers from other
disciplines\footnote{For example, \citep{Dong-2018} examines a similar type of model
under the name DiCCA (dynamic inner CCA), while \citep{Qin-2021}
proposes a LaVAR (latent vector autoregression) algorithm for achieving
a full canonical decomposition of the latent autoregressive structure in
\(Y\) allowing for interactions. First-order LAR models also bear a
strong resemblance to Min/Max Autocorrelation Factors (MAF) first
introduced in \citep{Switzer-1984} and since popularised in the
geosciences. A thorough comparison between these models is left to
future research.}. The main use case of these models outside of
economics of finance is the decomposition of multivariate sensor data
into time-persistent signals on the one hand, and serially uncorrelated
white noise on the other. In the context of investment management, the
same concept can be applied to derive trend following investment
strategies such as price momentum. Suppose the vector \(Y\) represents
the returns on a collection of investible assets, and \(Y \bm{\mathbf{w}}_i\)
represents the return on an investment strategy characterised by capital
allocation weights \(\bm{\mathbf{w}}_i\). The fixed point defined by
(\ref{eqn:clar_solution_full}) would then produce a vector of capital
allocation weights \(\bm{\mathbf{\hat{w}}}\) corresponding to an investment strategy with
the strongest possible price momentum or reversal signal, i.e., the
strongest possible correlation between past returns and future returns
in absolute terms, based on a sample of historical returns \(\mathbf{\bm{Y}}\).
A more refined result is obtained by examining the simplest type of LAR
problem, namely, a first-order autoregressive model of the form \(\tilde{y}_t
= c + \phi \tilde{y}_{t-1} + e_t\) with \(\tilde{y}_t = Y \bm{\mathbf{w}} + r_y\) and \(\tilde{y}_{t-1}
= A \bm{\mathbf{w}} + r_a\). In this case \(\bm{\phi}\) reduces to a scalar \(\phi\) and the
unconstrained solution becomes:
\begin{subnumcases}{\hspace{-19em}\label{eqn:lar1_solution_full}}
\bm{\mathbf{\hat{w}}} = \frac{1}{\rho_{y}} \mathbf{\bm{\Sigma}}_Y^{-1} \phi
\left[\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} - \phi \mathbf{\bm{\Sigma}}_{A} \right] \bm{\mathbf{w}}
\label{eqn:lar1_fp_1} \\[10pt]
\hat{\phi} = \frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}}}{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{A} \bm{\mathbf{w}}}
\label{eqn:lar1_fp_3} \\[10pt]
\rho_{y} = \frac{\phi }{\s2_y} \bm{\mathbf{w}}' \left[
\mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} -
\phi \mathbf{\bm{\Sigma}}_{A}
\right] \bm{\mathbf{w}}
\label{eqn:lar1_fp_2}
\end{subnumcases}
If we further assume that \(Y\) is covariance stationary, i.e., \(\mathbf{\bm{\Sigma}}_Y
= \mathbf{\bm{\Sigma}}_A\), \(\hat{\phi}\) reduces to \(\frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}}}{\s2_y}\) and
\(\rho_{y}\) can be refactored as:
\begin{flalign*}
& \rho_{y} = \frac{\phi }{\s2_y} \bm{\mathbf{w}}' \left[
\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} - \phi \mathbf{\bm{\Sigma}}_{A}
\right] \bm{\mathbf{w}}
= \frac{\phi }{\s2_y} \left[
\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} + \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YA} \bm{\mathbf{w}} - \phi \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{A} \bm{\mathbf{w}}
\right] & \\
& \quad = \phi \left[
\frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}}}{\s2_y} + \frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YA} \bm{\mathbf{w}}}{\s2_y}
- \phi \frac{\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{A} \bm{\mathbf{w}}}{\s2_y}
\right]
= \phi \left[ \phi + \phi - \phi \frac{\s2_y} {\s2_y} \right]
= \phi^2 &
\end{flalign*}
The solution for \(\bm{\mathbf{w}}\) then becomes:
\begin{equation} \label{eqn:w_eigenvec_lar1}
\hat{\phi} \bm{\mathbf{\hat{w}}} = \frac{1}{2} \mathbf{\bm{\Sigma}}_Y^{-1} \left[ \mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} \right] \bm{\mathbf{w}}
\end{equation}
or, equivalently:
\begin{equation} \label{eqn:w_eigenvec_lar1_b}
\hat{\phi} \bm{\mathbf{\hat{w}}} = \frac{1}{2} \left[ \mathbf{\bm{\Sigma}}_A^{-1} \mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_Y^{-1} \mathbf{\bm{\Sigma}}_{YA} \right] \bm{\mathbf{w}}
\end{equation}
In other words, if the proxy vector \(Y\) is covariance stationary, a
LAR(1) model specified by (\ref{eqn:lar1_solution_full}) can be estimated by
the eigendecomposition of the matrix \(\frac{1}{2} \mathbf{\bm{\Sigma}}_Y^{-1}
\left( {\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA}} \right)\), with each \(\bm{\mathbf{\hat{w}}}\) being its eigenvector
and \(\hat{\phi}\) the matching eigenvalue. An appropriate name for this
methodology would be ``Canonical Autocorrelation Analysis''\footnote{The term Canonical Autocorrelation Analysis was previously used
in \citep{Chen-2015} to describe a different type of statistical model.
However, given that the paper is in a different field, applying the same
term to equation (\ref{eqn:w_eigenvec_lar1}) should not invite any ambiguity.} (CAA)
because it produces a full canonical decomposition of the directions of
first-order autocorrelation in a group of covariance-stationary time
series variables \(Y\).
All components of equation (\ref{eqn:w_eigenvec_lar1}) have reasonably
intuitive interpretations. The matrices \(\mathbf{\bm{\Sigma}}_A^{-1} \mathbf{\bm{\Sigma}}_{AY}\) and
\(\mathbf{\bm{\Sigma}}_Y^{-1} \mathbf{\bm{\Sigma}}_{YA}\) are basically just matrices of least squares
regression coefficients: The first (second, third, \ldots{}) column of
\(\mathbf{\bm{\Sigma}}_A^{-1} \mathbf{\bm{\Sigma}}_{AY}\) is a vector of OLS coefficients from a
regression of the first (second, third, \ldots{}) variable in \(Y\) on all
the variables in \(A\). For the matrix \(\mathbf{\bm{\Sigma}}_Y^{-1} \mathbf{\bm{\Sigma}}_{YA}\) the roles
are reversed, i.e., its columns represent the OLS coefficient vectors
obtained from regressing each variable in \(A\) on all variables in
\(Y\). The coefficient \(\phi\) is both the autocorrelation coefficient
of \(\tilde{y}\) and the eigenvalue of \(\frac{1}{2} \mathbf{\bm{\Sigma}}_Y^{-1} \left( {
\mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA}} \right)\), which means that the first (last) eigenvector
produces the strongest (weakest) absolute autocorrelation signal.
In an investment management context, equation (\ref{eqn:w_eigenvec_lar1}) can
be used to extract momentum and reversal patterns from a collection of
investible assets. Let \(Y\) and \(A\) represent the returns on an
investment opportunity set at times \(t\) and \(t-1\), respectively.
Each eigenvector \(\bm{\mathbf{\hat{w}}}_i\) of the matrix \(\mathbf{\bm{\Sigma}}_Y^{-1} \left( {\mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA}} \right)\) would be a vector of capital allocation weights. The
investment strategies produced by all \(\bm{\mathbf{\hat{w}}}_i\) would capture all the
distinct directions of autocorrelation in \(Y\). The dominant
eigenvectors would produce investment strategies with the strongest
momentum or reversal signals, while the last eigenvectors would produce
strategies whose returns are closest to a random walk. Such a collection
of investment portfolios could be useful for a range of tasks, including
identifying market inefficiencies, predicting asset returns, and pricing
systematic investment risks.
\section{Empirical Application: Stock Markets and Economic Activity in the US}
\label{empirical_results}
As an example of how the LARX model can be used in the real world, let
us briefly examine the relationship between equity market performance
and real economic activity in the United States. A good starting point
for this analysis is provided by \citep{Ball-2021} who find that
de-trended levels of the S\&P 500 index have in-sample predictive power
over the de-trended levels of real US GDP.
The relationship between stock returns and real economic activity has a
strong foundation in economic theory. For example, the Consumption-based
Capital Asset Pricing Model (CCAPM) (\citep{Lucas-1978,Breeden-1979})
and the Investment CAPM (ICAPM) (see \citep{Zhang-2017} and the
references therein) put forward consumption and investment,
respectively, as the main sources of priced risk in the equity market.
With both consumption and investment also being key components in the
expenditure model of GDP, the relationship between stock returns and GDP
growth follows naturally.
In reality, however, market aggregates and macroeconomic aggregates are
designed to measure different things. Companies in the S\&P 500 are
weighted based on their market capitalisation rather than their relative
contributions to the real economy. US GDP is designed to measure
broad-based economic activity in the US -- not just that of the Fortune
500 companies. These differences suggest that the true strength of the
relationship between stock performance and real economic output would be
underestimated by a model linking the S\&P 500 to real US GDP.
The LARX model can be used to address this discrepancy. Both the S\&P 500
and US GDP are composite measures, which means they can be broken down
into their constituent parts and reassembled into latent measures of
market performance and economic growth, respectively. We can use the
five expenditure components of US GDP as proxy variables for a
supervised diffusion index (SDI) of US real economic activity (a ``real
activity SDI''), and 10 GICS\footnote{\href{https://www.msci.com/indexes/index-resources/gics}{Global Industry Classification Standard}} level 1 sector sub-indices\footnote{As of 2016, listed real estate (RE) was added as the eleventh
GICS level 1 sector of the S\&P 500. The RE sector is excluded from this
study for two reasons: Firstly, its data history only starts in Q4 2001
and would reduce the sample size from 138 quarterly observations to 90.
Secondly, the RE sector only has a 2.25\% weight in the S\&P 500 as of
April 2025 -- the second lowest weight after materials at 1.99\%.} of
the S\&P 500 as proxy variables for an SDI of market growth expectations
(a ``market SDI''). The component weights of the SDI measures can then
be compared to the component weights in the official aggregates to
explore whether the relationship between the S\&P 500 and real US GDP is
distorted by a misalignment in sector composition, the relative
importance of different sources of demand, or both.
\subsection{Data and Methodology}
\label{data_and_methodology}
Historical data for US GDP and its expenditure components are retrieved
from the Economic Database of the Federal Reserve Bank of St. Louis
(``FRED''). Historical index levels for the S\&P 500 and its GICS level 1
sector constituents are retrieved from Investing.com. A full data
reference can be found in Table \ref{table:data_reference}.
Empirical results from \citep{Ball-2021} suggest that the
best-performing model for real US GDP at time \(t\) contains the S\&P 500
at times \(t\) to \(t-3\) and two autoregressive lags, obtaining a
reported adjusted R-squared of 66.61\% for a quarterly data sample
between Q1 1999 and Q4 2020. Let us take the same specification as the
baseline but make three noteworthy changes to the experiment design:
First, model performance is going to be measured out of sample, rather
than in-sample. Second, log-percent changes will be used instead of
de-trended levels, which is a more common approach to measuring both
economic activity and stock market performance. Third, the sample period
will cover Q4 1989 to Q3 2025, capturing the full available data history
as of the time of writing\footnote{The main constraint on history length is the GICS classification
for the 10 original S\&P 500 sectors which was only introduced in 1999
with the data back-dated to 1989.}. Revised estimates will be used for all
economic aggregates following the tentative finding in \citep{Ball-2021}
that the link is stronger between equity performance and revised GDP
numbers as opposed to point-in-time (``vintage'') releases. A total of
four regression models will be compared:
\begin{subequations} \label{eqn:regression_models}
Baseline OLS/ARX model: Real GDP growth \(g\) vs S\&P 500 returns \(r\):
\begin{equation} \tag{\ref*{eqn:regression_models}}
g_t = c + \sum_{\tau = 1}^2 \phi_{t-\tau} g_{t-\tau} +
\sum_{\tau = 0}^3 \beta_{t-\tau} r_{t-\tau} + e
\end{equation}
LARX model a): Latent explanatory: Real GDP growth vs market SDI returns \(\tilde{r}\):
\begin{equation}\label{eqn:regression_g_nonlatent}
g_t = c + \sum_{\tau = 1}^2 \phi_{t-\tau} g_{t-\tau} +
\sum_{\tau = 0}^3 \beta_{t-\tau} \tilde{r}_{t-\tau} + e
\end{equation}
LARX model b): Latent dependent: Change in the real activity SDI \(\tilde{g}\) vs S\&P 500 returns:
\begin{equation}\label{eqn:regression_r_nonlatent}
\tilde{g}_t = c + \sum_{\tau = 1}^2 \phi_{t-\tau} \tilde{g}_{t-\tau} +
\sum_{\tau = 0}^3 \beta_{t-\tau} r_{t-\tau} + e
\end{equation}
LARX model c): Both latent: Change in the real activity SDI vs market SDI returns:
\begin{equation}\label{eqn:regression_all_latent}
\tilde{g}_t = c + \sum_{\tau = 1}^2 \phi_{t-\tau} \tilde{g}_{t-\tau} +
\sum_{\tau = 0}^3 \beta_{t-\tau} \tilde{r}_{t-\tau} + e
\end{equation}
\end{subequations}
For a cleaner comparison, variance constraints are imposed on the market
SDI and the real activity SDI to mimic the full-sample variance of the
S\&P 500 and real GDP growth, respectively. Although this information is
only available with the benefit of hindsight, the choice of variance
constraint does not affect the performance of the LARX model as it is
scale invariant. All regressions use exponentially decaying sample
weights with a half-life of 10 years, in order to allow for some drift
in the regression parameters. The COVID lockdown period of Q2 and Q3
2020 is treated as a statistical outlier and removed from the data
sample\footnote{US GDP shows a contraction of 8.2\% in Q2 2020 followed by a 7.5\%
rebound in Q3 2020 -- a -13.8 sigma event and a 12.6 sigma event,
respectively, based on the standard deviation of quarterly US GDP growth
excluding these two quarters.}. A minimum of 40 degrees of freedom is set as a requirement
for producing a forecast, which corresponds to 10 years of quarterly
data on top of one data point lost to each estimated coefficient,
including Lagrange multipliers. An additional three data points are lost
to the lag operator and one to the percent change calculation.
Ultimately, forecast coverage starts in Q3 2002 for the baseline model
(longest) and in Q3 2006 for the model with all latent variables
(shortest).
\subsection{Out-of-Sample Forecasting Performance}
\label{sec:org885f59d}
Figure \ref{NEW_fig-fc_ols_vs_clarx} plots the rolling out-of-sample (OOS)
predictions made by the baseline model (top left), as well as LARX
models a) (top right), b) (bottom left), and c) (bottom right). Each
plot overlays the actual values of the dependent (real GDP growth in the
top row, the real activity SDI in the bottom row), as well as the naïve
forecast for the dependent based a rolling sample mean (a.k.a. the
``benchmark''). The grey text boxes report each model's OOS R-squared
statistic as defined in \citep{Campbell-2007}.
\begin{figure}[ht]
\includegraphics[width=\textwidth]{NEW_fig-fc_ols_vs_clarx.pdf}
\caption{Out-of-sample forecasting performance} \label{NEW_fig-fc_ols_vs_clarx}
\floatfoot{ \footnotesize
Note: Rolling out-of-sample forecasting performance of the four
models specified by equations
(\ref{eqn:regression_models})-(\ref{eqn:regression_all_latent}).
The target variance for the real activity SDI is set to the
full-sample variance of quarterly US real GDP growth (approx.
0.000124). The variance of the market SDI is set to the
full-sample variance of the S\&P 500 quarterly log-return
(approx. 0.006348).
}
\end{figure}
The baseline model from \citep{Ball-2021} (top left) does well despite
the changes to the experiment design. It achieves an OOS R-squared of
50.3\% -- reasonably close to the in-sample adjusted R-squared of 66.61\%
reported in the original study.
All three LARX specifications further improve on the baseline model.
LARX model a) (market SDI and real GDP growth -- top right) yields an
OOS R-squared of 63.9\%. LARX model b) (the S\&P 500 and the real activity
SDI -- bottom left) produces an OOS R-squared of 67.0\%. LARX model c)
(marked SDI and real activity SDI -- bottom right) performs best with an
OOS R-squared of 79.7\%.
For completeness, it must be noted that the dependent variable in models
b) and c) is different from the dependent variable in models a) and b).
However, the comparison is still meaningful for at least three reasons.
First, the variance of the real activity SDI matches the variance of
real US GDP growth by design. Second, the goal of the study is to assess
the importance of a common unobserved factor in both the stock market
and the real economy, which is a bi-directional problem requiring
optimisation on both ends. Third, rearranging equations
(\ref{eqn:regression_models})-(\ref{eqn:regression_all_latent}) to use
the stock market return as the dependent yields a largely similar set of
results, as reported in \ref{empirical_results_rev}.
\subsection{Explaining the Outperformance of the LARX Models}
\label{sec:orgef6b772}
We can dive deeper into the empirical results for additional insight.
Overall, LARX models a) and b) both outperform the baseline model by
roughly the same amount, while LARX model c) yields approximately the
same improvement as LARX models a) and b) combined. This suggests three
things: First, the sector composition of the market SDI better reflects
the sector composition of the real economy. Second, the expenditure
composition of the real activity SDI better reflects the main sources of
demand for the output of Fortune 500 companies. Third, sector
composition and expenditure composition represent two equally important
yet distinct sources of discrepancy between the stock market and the
real economy, as fixing one has little to no bearing on the extent of
other.
To gain a better understanding of the underlying dynamics, let us look
at the weight vectors of the SDI measures obtained from LARX model c)
(the best-performing specification) side by side with their non-latent
counterparts.
\begin{figure}[ht]
\includegraphics[width=\textwidth]{NEW_fig-weight_evolution_spx.pdf}
\caption{Sector composition: Latent Measure of Market Growth Expectations vs the S\&P 500}
\label{NEW_fig-weight_evolution_spx}
\end{figure}
Figure \ref{NEW_fig-weight_evolution_spx} compares the sector weights of the
market SDI (left) with the sector weights of the S\&P 500 (right)
approximated using a rolling regression of S\&P 500 returns on the
coincident returns of the sector sub-indices. The weight vector of the
market SDI is scaled so that the positive weights add up to 100\%.
Here, two observations jump out. First, sector rotations, a.k.a.
relative sector performance, comprise between 40\% and 80\% of the market
SDI at any given point in time. This suggests that sector rotations are
about as useful for forecasting real economic activity as is the overall
direction of the equity market. Sector rotations have been previously
examined in the investment management literature for the purposes of
capturing systematic risk premia (see \citep{Lamponi-2014}). A link
between sector rotations and the business cycle has been examined for
the purposes of constructing systematic investment strategies by
\citep{Molchanov-2024}, albeit with limited success. To the author's
best knowledge, no previous study demonstrates that sector rotations
have out-of-sample predictive power with respect to real economic
growth.
Second, the composition of the market SDI fluctuates quite strongly over
time, and many of these fluctuations seem to correspond to major
economic trends or events. For example, the performance of the energy
(``En.'') sector changes its economic meaning after the introduction of
major green energy initiatives by the Obama administration in 2009, and
again around the commodity super-cycle of 2015-2016. Similarly, the
healthcare sector (``HC'') has a positive weight until the introduction
of the Affordable Care Act (a.k.a. ``Obamacare'') in 2010 and a negative
weight thereafter.
\begin{figure}[ht]
\includegraphics[width=\textwidth]{NEW_fig-weight_evolution_gdp.pdf}
\caption{Expenditure composition: Latent Economic Activity Measure vs US GDP}
\label{NEW_fig-weight_evolution_gdp}
\end{figure}
Figure \ref{NEW_fig-weight_evolution_gdp} plots the evolution of the expenditure
weights in the real activity SDI (left) against the evolution of the
expenditure weights in real US GDP (right).
Both measures assign the largest weight to consumer spending
(``Cons.''), which is an expected result for a consumption-led economy
like the United States. The role of private investment (``Inv.'') is
smaller and largely similar for both aggregates.
The biggest differences lie in the role of government spending
(``Govt.'') and the trade balance (``Imp.'', ``Exp.''). In the real
activity SDI, stronger government spending is associated with weaker
equity market performance, and vice versa, outside of the immediate
aftermath of the 2008 crisis. In other words, equity investors seem to
welcome fiscal expansion in a crisis, but view it as a negative in
normal times, a.k.a. the Keynesian notion of fiscal policy as a
counter-cyclical buffer. The effect becomes more pronounced after the
outbreak of COVID-19, perhaps owing to a ballooning US fiscal deficit in
2020 and a sharp rise in interest rates since 2021.
Lastly, the interplay between the trade balance and the equity market
seems to be much more complex than the role of imports and exports in
the national accounts. \citep{Reinbold-2019} offers a brief primer on
the history of the US trade balance across the value chain at different
stages of the country's industrialisation. From an accounting
perspective, imports subtract from GDP while exports add to it, both
with a multiplier of 1. However, today many large US companies have
manufacturing facilities in other countries, which means that their
products need to be imported before they are sold domestically, but they
don't need to exported to be sold abroad. As a result, the link between
stock returns and exports is weakened, whereas the link between stock
returns and imports becomes positive inasmuch as imports become a
leading indicator of domestic sales. The real activity SDI captures this
relationship by assigning a time-varying weight to US exports and a
consistently positive weight to US imports, with the exception of the
latest quarter in the sample (3Q 2025) characterised by major swings in
US tariff policy.
\subsection{Questions for Future Research}
\label{sec:org57fb151}
A more thorough examination of the relationship between stock returns
and real economic activity is beyond the scope of this paper; however,
several avenues of further investigation may be of interest to future
research. Three examples are listed below.
First, further insight may be gained by constructing the market SDI and
the real activity SDI at higher levels of granularity. The market SDI
can be estimated using GICS level 2 indices or even single stocks, while
the real activity SDI can go all the way down to the itemised national
accounts. Of course, the benefits of a more granular approach should be
weighed against a higher risk of overfitting.
Second, this paper's findings are not examined from an asset pricing
standpoint. No attempts are made to test whether the market SDI produced
by the LARX model is a valid asset pricing factor in the CCAPM or ICAPM
sense, a profitable investment strategy, or a consistent leading signal
for business cycle rotations.
Third, only the case of the United States is examined. It may be
informative to perform similar studies for countries with other economic
models and sector compositions, in order to see whether the relationship
between stock returns and domestic economic activity is a universal
phenomenon.
\section{Concluding remarks}
\label{concluding_remarks}
This paper proposes a new methodological framework for estimating latent
variable models called Supervised Diffusion (SDF). A fixed point
solution is derived for a new supervised diffusion model (SDM) called
LARX: a superset of the traditional autoregressive model with exogenous
inputs (ARX) in which any or all variables can be latent. The
derivations of the LARX model result in a minor contribution to the
field of matrix calculus: A block-wise direct sum operator is introduced
and applied to solve a class of Lagrangian optimisation problems with
interactions between multiple coefficient vectors in the presence of
piecemeal constraints. Several special cases of the LARX models are also
examined in more detail, including a parsimonious lead-lag regression
model called Latent Shock Regression (LSR), as well as a new canonical
decomposition technique called Canonical Autocorrelation Analysis (CAA).
In the empirical section, the LARX model is used to re-examine the link
between stock market performance and real economic activity in the
United States. The LARX model attains an out-of-sample (OOS) R-squared
of 79.7\%, compared to an OOS R-squared of 50.3\% for the baseline OLS
specification. The LARX model also provides novel insights about the
interplay between stock returns and the real economy, corroborating a
few notions which were previously best described as ``rules of thumb''.
For example, new evidence is found in support of the importance of
sector rotations in predicting economic growth, as well as the complex
and multifaceted roles of fiscal policy and international trade in the
economic value chain.
SDMs like LARX have many potential use cases in economics and finance.
For example, they can be used to obtain better estimates for variables
whose main goal is to track other processes, like diffusion indices of
business activity (e.g., \citep{Owens-2005}), surveys of consumer
sentiment (e.g., \citep{Curtin-2000}), or composite indicators of
financial stress (e.g., \citep{CISS-2012}). Supervised diffusion indices
can also replace various heuristic investment techniques, such as
``buying past winners and selling past losers'' as a means of capturing
asset price momentum as in \citep{Jegadeesh-1993}.
In the scientific process, SDMs have a place in research settings
characterised by noisy and unreliable empirical data and little scope
for controlled experiments -- a common problem in macroeconomics. SDMs
can enable the researcher to place relatively more trust in the research
hypothesis by expanding the search for its supporting evidence to
arbitrary linear combinations over the observed data series. Formally
speaking, the use SDMs reduces the risk of Type II error at the expense
of a somewhat higher chance of Type I error.
Sections \ref{sdf}-\ref{larx_special_cases} leave much scope for further
methodological work. For example, the topics of statistical significance
and feature selection\footnote{In the specific case of CCA, \citep{Bagozzi-1981} and
\citep{Ahn-2018} offer good starting points in the analysis of
statistical significance and feature selection, respectively.} were only briefly touched upon in Section
\ref{larx_ols_interpretation}. Furthermore, the LARX methodology lends itself
well to a number of adjustments which are common practice with the
traditional ARX model, including moving average errors (MA), conditional
heteroskedasticity (GARCH), seasonality patterns, various forms of
coefficient regularisation (e.g., LASSO, Ridge, Elastic Net - see, for
example, \citep{Vinod-1976}), and various covariance adjustment
techniques such as Generalised Least Squares (e.g., \citep{Aitken-1936})
and portfolio-style covariance shrinkage (e.g., see
\citep{Ledoit-Wolf-2022}).
Non-linear SDMs may warrant a closer look as well. For example, in the
context of the asset pricing theory, the LARX model can be used for
estimating risk factors such as price momentum (\citep{Jegadeesh-1993}),
earnings momentum (\citep{Barth-1999}) and company size
(\citep{Fama-French-1993}), which are governed by linear relationships.
However, it cannot be used as easily to estimate factors such as value
(e.g., \citep{Asness-2013,Fama-French-1993,Shiller-1988}) or quality
(e.g., see \citep{Hsu-2017}). Valuations are, generally speaking, ratios
of prices to fundamentals, e.g., a portfolio's price-to-earnings (PE)
valuation can be calculated as \(\frac{P \bm{\mathbf{w}}}{EPS \bm{\mathbf{w}}}\), where \(P\) and
\(EPS\) are vectors of company share prices and earnings per share, and
\(\bm{\mathbf{w}}\) is a vector of capital allocation weights. A value factor SDI
would then be governed by a functional relationship similar to:
\begin{equation} \label{eqn:reg_value_factor}
Y_t \bm{\mathbf{w}} = c + \frac{P_{t-1} \bm{\mathbf{w}}}{EPS_{t-1} \bm{\mathbf{w}}} + \epsilon_t
\end{equation}
where \(Y\) is a vector of constituent returns, \(c\) is the intercept,
\(\epsilon\) is the error term, and \(t\) is a time subscript. The
quality factor, on the other hand, can be based on earnings volatility
(e.g., see \citep{Dichev-2009}), which would necessitate a quadratic
SDM.
To summarise, SDMs such as LARX can be viewed as a rather natural
evolution of traditional regression analysis. At the same time, an
effective application of these models requires a slight paradigm shift
on the part of the researcher, because their input variables are
designed to change shape according to the functional relationship at
hand. This feature may not be universally useful, but it does create an
opportunity to re-examine a number of existing models and discover new
relationships in economics and finance, as well as in a number of other
data-intensive fields.
\begin{appendix}
\section{Commutativity of the blockwise direct sum operator}
\label{bds_vs_mmul}
\begin{proposition} \label{prop:bds_vs_mmul}
Let \(S\) be a set of all matrix sequences of length \(k\), and let the
sequence \(\left< {\mathbf{\bm{A}}} \right> \equiv \left< {\mathbf{\bm{A}}_i | 1 \leq i \leq k} \right>\) be an element
in \(S\). Let the matrix \(\mathbf{\bm{A}}^{\oplus}\) represent the direct sum over the
elements in \(\left< {\mathbf{\bm{A}}} \right>\). For any function \(f: S \rightarrow S\), if
\(f \left( {\left< {\mathbf{\bm{A}}} \right>} \right)\) can be expressed as a sequence \( \left< {\mathbf{\bm{M}}_i \mathbf{\bm{A}}_i} \right>
\equiv \left< {\mathbf{\bm{M}}_i \mathbf{\bm{A}}_i | 1 \leq i \leq k} \right>\) for some arbitrary sequence
of matrices \(\left< {\mathbf{\bm{M}}_i | 1 \leq i \leq k} \right>\), then \( f
\left(\left< {\mathbf{\bm{A}}} \right> \right)^{\oplus}_v = \left[ f \left( \left< {\mathbf{\bm{A}}^{\oplus}_v} \right> \right)
\right]_v \). If \( f \left(\left< {\mathbf{\bm{A}}} \right>\right) \) can be expressed as a
sequence \(\left< {\mathbf{\bm{A}}_i \mathbf{\bm{M}}_i} \right> \equiv \left< {\mathbf{\bm{A}}_i \mathbf{\bm{M}}_i | 1 \leq i \leq k} \right>\)
for an arbitrary sequence of matrices \(\left< {\mathbf{\bm{M}}_i | 1 \leq i \leq k} \right>\),
then \( f \left( { \left< {\mathbf{\bm{A}}} \right> } \right)^{\oplus}_h = \left[ f \left( { \left< {\mathbf{\bm{A}}^{\oplus}_h} \right>
} \right) \right]_h \). \end{proposition}
\begin{proof}
For a block matrix \(\mathbf{\bm{A}}\), take \(\left< {\mathbf{\bm{A}}} \right>\) to denote the sequence of
the blocks in \(\mathbf{\bm{A}}\). For a sequence of matrices \(\left< {\mathbf{\bm{A}}} \right>\), denote
its \(i\)'th element by \(\left< {\mathbf{\bm{A}}} \right>_i\). For the case of
left-multiplication we then have:
\begin{flalign*}
& f \left( { \left< {\mathbf{\bm{A}}} \right> } \right)^{\oplus}_v
= \left< {\mathbf{\bm{M}}_i \mathbf{\bm{A}}_i} \right>^{\oplus}_v
= \left( { \begin{array}{cccc}
\mathbf{\bm{M}}_1 \mathbf{\bm{A}}_1, & \0, & \cdots & \0 \\ \hline
\0, & \mathbf{\bm{M}}_2 \mathbf{\bm{A}}_2, & \cdots & \0 \\ \hline
\vdots & \vdots & \ddots & \vdots \\ \hline
\0, & \0, & \cdots & \mathbf{\bm{M}}_k \mathbf{\bm{A}}_k
\end{array} } \right) & \\[10pt]
& \left[ f \left( { \left< {\mathbf{\bm{A}}^{\oplus}_v} \right> } \right) \right]_v
= \left[ \left< {\mathbf{\bm{M}}_i \left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_i | 1 \leq i \leq k} \right>q i \leq k} \right]_v
= \left( { \begin{array}{c}
\mathbf{\bm{M}}_1 \left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_1 \\ \hline
\mathbf{\bm{M}}_2 \left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{M}}_k \left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_k
\end{array} } \right) = & \\[10pt]
& \hspace{1em} = \left( { \begin{array}{c}
\mathbf{\bm{M}}_1 \begin{pmatrix} \mathbf{\bm{A}}_1, & \0, & \cdots, & \0 \end{pmatrix} \\ \hline
\mathbf{\bm{M}}_2 \begin{pmatrix} \0, & \mathbf{\bm{A}}_2, & \cdots, & \0 \end{pmatrix} \\ \hline
\vdots \\ \hline
\mathbf{\bm{M}}_k \begin{pmatrix} \0, & \0, & \cdots, & \mathbf{\bm{A}}_k \end{pmatrix}
\end{array} } \right)
= \left( { \begin{array}{cccc}
\mathbf{\bm{M}}_1 \mathbf{\bm{A}}_1, & \0, & \cdots & \0 \\ \hline
\0, & \mathbf{\bm{M}}_2 \mathbf{\bm{A}}_2, & \cdots & \0 \\ \hline
\vdots & \vdots & \ddots & \vdots \\ \hline
\0, & \0, & \cdots & \mathbf{\bm{M}}_k \mathbf{\bm{A}}_k
\end{array} } \right) = f \left( {\left< {\mathbf{\bm{A}}} \right>} \right)^{\oplus}_v &
\end{flalign*}
For the case of right-multiplication we have:
\begin{flalign*}
& f \left( {\left< {\mathbf{\bm{A}}} \right>} \right)^{\oplus}_h
= \left< {\mathbf{\bm{A}}_i \mathbf{\bm{M}}_i} \right>^{\oplus}_h
= \left( { \begin{array}{c|c|c|c}
\mathbf{\bm{A}}_1 \mathbf{\bm{M}}_1 & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{A}}_2 \mathbf{\bm{M}}_2 & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{A}}_k \mathbf{\bm{M}}_k \\
\end{array} } \right) & \\[10pt]
& \left[ f \left( {\left< {\mathbf{\bm{A}}^{\oplus}_h} \right>} \right) \right]_h
= \left[ \left< {\left< {\mathbf{\bm{A}}^{\oplus}_h} \right>_i \mathbf{\bm{M}}_i | 1 \leq i \leq k} \right>q i \leq k} \right]_v
= \left[ \begin{array}{c|c|c|c}
\left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_1 \mathbf{\bm{M}}_1 &
\left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_2 \mathbf{\bm{M}}_2 &
\cdots &
\left< {\mathbf{\bm{A}}^{\oplus}_v} \right>_k \mathbf{\bm{M}}_k
\end{array} \right] = & \\[10pt]
& \hspace{1em} = \left( { \begin{array}{c|c|c|c}
\begin{bmatrix} \mathbf{\bm{A}}_1 \\ \0 \\ \vdots \\ \0 \end{bmatrix} \mathbf{\bm{M}}_1 &
\begin{bmatrix} \0 \\ \mathbf{\bm{A}}_2 \\ \vdots \\ \0 \end{bmatrix} \mathbf{\bm{M}}_2 &
\cdots &
\begin{bmatrix} \0 \\ \0 \\ \vdots \\ \mathbf{\bm{A}}_k \end{bmatrix} \mathbf{\bm{M}}_k
\end{array} } \right)
= \left( { \begin{array}{c|c|c|c}
\mathbf{\bm{A}}_1 \mathbf{\bm{M}}_1 & \0 & \cdots & \0 \\
\0 & \mathbf{\bm{A}}_2 \mathbf{\bm{M}}_2 & \cdots & \0 \\
\vdots & \vdots & \ddots & \vdots \\
\0 & \0 & \cdots & \mathbf{\bm{A}}_k \mathbf{\bm{M}}_k \\
\end{array} } \right) = & \\[10pt]
& \hspace{1em} = f \left( { \left< {\mathbf{\bm{A}}} \right> } \right)^{\oplus}_h &
\end{flalign*}
\end{proof}
\section{Blockwise Kronecker Product Factorisation for Vectors}
\label{block_kron_factorisation}
\begin{proposition} \label{prop:block_kron_factorisation}
Let \(\mathbf{\bm{a}}\) and \(\bm{b}\) be two column vectors, each comprised of \(k\)
row blocks of arbitrary lengths. The blockwise Kronecker product \(\mathbf{\bm{a}}
\odot \bm{b}\) can be factorised as \(\mathbf{\bm{a}} \odot \bm{b} = \left( {\mathbf{\bm{a}} \odot
I_{\bm{b}}} \right) \bm{b} = \left( { I_{\mathbf{\bm{a}}} \odot \bm{b} } \right) \mathbf{\bm{a}}\), where \(I_{\bm{b}}\) and
\(I_{\mathbf{\bm{a}}}\) are identity matrices with the same number of rows and row
block structure as \(\bm{b}\) and \(\mathbf{\bm{a}}\), respectively.
\end{proposition}
\begin{proof}
Let the vector \(\mathbf{\bm{a}}\) have dimensions \(M \times 1\) and \(\left< {\mathbf{\bm{a}}} \right>
\equiv \left< {\underset{m_i \times 1}{\mathbf{\bm{a}}_i} | 1 \leq i \leq k} \right>\) be the
sequence of vectors which represent the row blocks in \(\mathbf{\bm{a}}\) such that
\(\sum_{i = 1}^k{m_i} = M\). Similarly, let the vector \(\bm{b}\) have
dimensions \(N \times 1\) and the sequence of vectors \(\left< {\bm{b}} \right> \equiv
\left< {\underset{n_i \times 1}{\bm{b}_i} | 1 \leq i \leq k} \right>\) represent the
row blocks in \(\bm{b}\) such that \(\sum_{i = 1}^k{n_i} = N\).
The blockwise Kronecker product \(\mathbf{\bm{a}} \odot \bm{b}\) can then be defined as:
\begin{flalign*}
& \mathbf{\bm{a}} \odot \bm{b} = \left( {
\begin{array}{c}
\mathbf{\bm{a}}_1 \otimes \bm{b}_1 \\ \hline
\mathbf{\bm{a}}_2 \otimes \bm{b}_2 \\ \hline
\vdots \\ \hline
\mathbf{\bm{a}}_k \otimes \bm{b}_k
\end{array}
} \right) &
\end{flalign*}
Note that by the properties of the Kronecker product the following holds
for any two matrices \(\mathbf{\bm{A}}\) and \(\mathbf{\bm{B}}\):
\begin{flalign} \label{eqn:kron_factorisation_matrices}
& \underset{m \times p}{\mathbf{\bm{A}}} \otimes \underset{n \times q}{\mathbf{\bm{B}}}
= \left( {\mathbf{\bm{A}} \otimes I_n} \right) \left( {I_p \otimes \mathbf{\bm{B}}} \right)
= \left( {I_m \otimes \mathbf{\bm{B}}} \right) \left( {\mathbf{\bm{A}} \otimes I_q} \right) &
\end{flalign}
In the special case of a Kronecker product between two vectors, \(p\)
and \(q\) reduce to \(1\) and the identity matrices \(I_p\) and \(I_q\)
become \(1\) by association. As a result, for any given Kronecker
product \(\mathbf{\bm{a}}_i \otimes \bm{b}_i\), the following holds:
\begin{flalign} \label{eqn:kron_factorisation_vectors}
& \underset{m_i \times 1}{\mathbf{\bm{a}}_i} \otimes \underset{n_i \times 1}{\bm{b}_i}
= \left( {\mathbf{\bm{a}}_i \otimes I_p} \right) \bm{b}_i
= \left( {I_m \otimes \bm{b}_i} \right) \mathbf{\bm{a}}_i &
\end{flalign}
This allows us to rewrite the blockwise Kronecker product \(\mathbf{\bm{a}} \odot
\bm{b}\) in two alternative ways:
\begin{subequations}
\begin{flalign}
& \mathbf{\bm{a}} \odot \bm{b} = \left( {
\begin{array}{c}
\left( {\mathbf{\bm{a}}_1 \otimes I_{n_1}} \right) \bm{b}_1 \\ \hline
\left( {\mathbf{\bm{a}}_2 \otimes I_{n_2}} \right) \bm{b}_2 \\ \hline
\vdots \\ \hline
\left( {\mathbf{\bm{a}}_k \otimes I_{n_k}} \right) \bm{b}_k
\end{array}
} \right)right) \bm{b}_k
\end{array}
} & \label{eqn:bk_deriv_1a} \\[15pt]
& \mathbf{\bm{a}} \odot \bm{b} = \left( {
\begin{array}{c}
\left( {I_{m_1} \otimes \bm{b}_1} \right) \mathbf{\bm{a}}_1 \\ \hline
\left( {I_{m_2} \otimes \bm{b}_2} \right) \mathbf{\bm{a}}_2 \\ \hline
\vdots \\ \hline
\left( {I_{m_k} \otimes \bm{b}_k} \right) \mathbf{\bm{a}}_k
\end{array}
} \right)\right) \mathbf{\bm{a}}_k
\end{array}
} & \label{eqn:bk_deriv_1b}
\end{flalign}
\end{subequations}
Define a sequence of matrices \(\left< {\mathbf{\bm{a}}_i \otimes I_{n_i} | 1 \leq i
\leq k} \right> \equiv \left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right>\), and another
sequence \(\left< {I_{m_i} \otimes \bm{b}_i | 1 \leq i \leq k} \right> \equiv
\left< {I_{m_i} \otimes \bm{b}_i} \right> \). We can then rewrite
(\ref{eqn:bk_deriv_1a}) and (\ref{eqn:bk_deriv_1b}) using the blockwise
direct sum operator introduced in Section \ref{blockwise_direct_sum}:
\begin{subequations}
\begin{flalign}
& \mathbf{\bm{a}} \odot \bm{b} = \left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right>^{\oplus} \bm{b} & \label{eqn:bk_deriv_2a} \\[10pt]
& \mathbf{\bm{a}} \odot \bm{b} = \left< {I_{m_i} \otimes \bm{b}_i} \right>^{\oplus} \mathbf{\bm{a}} & \label{eqn:bk_deriv_2b}
\end{flalign}
\end{subequations}
It then remains to show that \(\left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right>^{\oplus}\) and
\(\left< {I_{m_i} \otimes \bm{b}_i} \right>^{\oplus}\) can be written as \( \mathbf{\bm{a}} \odot
I_{\bm{b}} \) and \( I_{\mathbf{\bm{a}}} \odot \bm{b} \), respectively. This can be done
with the help of Proposition \ref{prop:bds_vs_mmul}.
First of all, consider the sequences of identity matrices
\(\left< {I_{m_i}} \right> \equiv \left< {I_{m_i} | 1 \leq i \leq k} \right>\) and
\(\left< {I_{n_i}} \right> \equiv \left< {I_{n_i} | 1 \leq i \leq k} \right>\) on a standalone
basis. The sequences \(\left< {I_{m_i} \otimes \bm{b}_i} \right>\) and \(\left< {\mathbf{\bm{a}}_i
\otimes I_{n_i}} \right>\) can then be expressed as functions \(f: S \rightarrow
S\) and \(g: S \rightarrow S\) where \(S\) represents the set of all
matrix sequences of length \(k\). such that:
\begin{flalign*}
& f \left( {\left< {I_{m_i}} \right>} \right)
= \left< {\left< {I_{m_i}} \right>_i \otimes \left< {\bm{b}} \right>_i | 1 \leq i \leq k} \right>>_i | 1 \leq i \leq k}
\equiv \left< {I_{m_i} \otimes \bm{b}_i | 1 \leq i \leq k} \right>
\equiv \left< {I_{m_i} \otimes \bm{b}_i} \right> & \\[10pt]
& g \left( {\left< {I_{n_i}} \right>} \right)
= \left< {\left< {\mathbf{\bm{a}}} \right>_i \otimes \left< {I_{n_i}} \right>_i | 1 \leq i \leq k} \right>>_i | 1 \leq i \leq k}
\equiv \left< {\mathbf{\bm{a}}_i \otimes I_{n_i} | 1 \leq i \leq k} \right>
\equiv \left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right> &
\end{flalign*}
Proposition \ref{prop:bds_vs_mmul} applies because we can express
\(\left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right>\) and \(\left< {I_{m_i} \otimes \bm{b}_i} \right>\) as
left-multiplications over the sequences of identity matrcies
\(\left< {I_{m_i}} \right>\) and \(\left< {I_{n_i}} \right>\):
\begin{flalign*}
& \left< {I_{m_i} \otimes \bm{b}_i} \right> = \left< {\left( {I_{m_i} \otimes \bm{b}_i} \right) I_{m_i}} \right> & \\[10pt]
& \left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right> = \left< {\left( {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right) I_{n_i}} \right> &
\end{flalign*}
It then follows that:
\begin{flalign*}
& f \left( {\left< {I_{m_i}} \right>} \right)^{\oplus}_v
= \left[ f \left( {\left< {\left< {I_{m_i}} \right>^{\oplus}_v} \right>{\oplus}_v}} \right) \right]_v & \\[10pt]
& g \left( {\left< {I_{n_i}} \right>} \right)^{\oplus}_v
= \left[ g \left( {\left< {\left< {I_{n_i}} \right>^{\oplus}_v} \right>{\oplus}_v}} \right) \right]_v &
\end{flalign*}
Next, note that \(\left< {I_{m_i}} \right>^{\oplus}_v\) produces an identity matrix
of size \(M\) with a row block structure of \(\mathbf{\bm{a}}\), while
\(\left< {I_{n_i}} \right>^{\oplus}_v\) produces an identity matrix of size \(N\)
with a row block structure of \(\bm{b}\). In other words,
\(\left< {I_{m_i}} \right>^{\oplus}_v = I_{\mathbf{\bm{a}}}\) and \(\left< {I_{m_i}} \right>^{\oplus}_v =
I_{\bm{b}}\). This means:
\begin{flalign*}
& \left[ f \left( {\left< {\left< {I_{m_i}} \right>^{\oplus}_v} \right>{\oplus}_v}} \right) \right]_v
= \left[ f \left( {\left< {I_{\mathbf{\bm{a}}}} \right>} \right) \right]_v
= \left[ \left< { \left< { I_{\mathbf{\bm{a}}} } \right>_i \otimes \bm{b}_i } \right>mes \bm{b}_i } \right]_v & \\[10pt]
& \left[ g \left( {\left< {\left< {I_{n_i}} \right>^{\oplus}_v} \right>{\oplus}_v}} \right) \right]_v
= \left[ g \left( {\left< {I_{\bm{b}}} \right>} \right) \right]_v
= \left[ \left< { \mathbf{\bm{a}}_i \otimes \left< {I_{\bm{b}}} \right>_i } \right>\right>_i } \right]_v &
\end{flalign*}
Lastly, we note that the matrix representation of a pairwise Kronecker
product over two sequences of matrices is the same as a blockwise
Kronecker product if the sequences represent matrix blocks along the
same dimension. For example, \(\left[ \left< { \mathbf{\bm{a}}_i \otimes \left< {I_{\bm{b}}} \right>_i
} \right>\right>_i
} \right]_v = \mathbf{\bm{a}} \odot I_{\bm{b}}\) because \(\mathbf{\bm{a}}_i\) and
\(\left< {I_{\bm{b}}} \right>_i\) are row blocks of \(\mathbf{\bm{a}}\) and \(I_{\bm{b}}\),
respectively.
Putting all the steps together we have:
\begin{align}
& \mathbf{\bm{a}} \odot \bm{b} = \left< {\mathbf{\bm{a}}_i \otimes I_{n_i}} \right>^{\oplus} \bm{b}
= \left[ \left< {\mathbf{\bm{a}}_i \otimes \left< {\left< {I_{n_i}} \right>^{\oplus}_v} \right>{\oplus}_v}_i} \right> \right>{\oplus}_v}_i} \right]_v \bm{b}
= \left[ \left< { \mathbf{\bm{a}}_i \otimes \left< {I_{\bm{b}}} \right>_i } \right>\right>_i } \right]_v \bm{b}
= \left( { \mathbf{\bm{a}} \odot I_{\bm{b}} } \right) \bm{b} \label{eqn:bk_deriv_3a} \\[10pt]
& \mathbf{\bm{a}} \odot \bm{b} = \left< {I_{m_i} \otimes \bm{b}_i} \right>^{\oplus} \mathbf{\bm{a}}
= \left[ \left< {\left< {\left< {I_{m_i}} \right>^{\oplus}_v } \right>\oplus}_v }_i \otimes \bm{b}_i} \right>}_v }_i \otimes \bm{b}_i} \right]_v \mathbf{\bm{a}}
= \left[ \left< {\left< {I_{\bm{b}}} \right>_i \otimes \bm{b} } \right>times \bm{b} } \right]_v \mathbf{\bm{a}}
= \left( { I_{\mathbf{\bm{a}}} \odot \bm{b} } \right) \mathbf{\bm{a}} \label{eqn:bk_deriv_3b}
\end{align}
\end{proof}
\section{Conditions for the Redundancy of PCA Decomposition}
\label{pca_redundancy}
\begin{proposition} \label{prop:pca_redundancy}
Take a \(\left( {1 \times m} \right)\) vector of observed factors \(X\) whose
covariance matrix has full rank. Define the PCA-based diffusion indices
obtained from \(X\) as \(\widetilde{X}^{di} := X \mathbf{\bm{\Omega}}^{di}\), where \(\mathbf{\bm{\Omega}}^{di}\) is
an \(\left( {m \times m} \right)\) matrix of principal component weights. Assume
WLOG that the ``true'' latent factor \(\tilde{x}^*\) that drives the dependent
can be expressed as a linear combination of the diffusion indices such
that \(\tilde{x}^* := \widetilde{X}^{di} \bm{\omega}^*\) for some weight vector \(\bm{\omega}^*\). It then
follows from the properties of PCA that \(\tilde{x}^*\) can be expressed in
terms of \(X\) directly. \end{proposition}
\begin{proof}
By the properties of PCA, \(\mathbf{\bm{\Omega}}^{di}\) is a rotation matrix which means
it is invertible with \(\left( {\mathbf{\bm{\Omega}}^{di}} \right)^{-1} = \left( {\mathbf{\bm{\Omega}}^{di}} \right)'\). It then
follows that any valid solution for \(\bm{\omega}^*\) in the DI space can be
expressed as \(\bm{\omega}^* = \left( {\mathbf{\bm{\Omega}}^{di}} \right)' \bm{\omega}^{**}\) for an equally valid
weight vector \(\bm{\omega}^{**}\) in the observed variable space. Formally:
\begin{equation*}
\tilde{x}^* = \widetilde{X}^{di} \bm{\omega}^* = X \mathbf{\bm{\Omega}}^{di} \bm{\omega}^* = X \mathbf{\bm{\Omega}}^{di}
\left( {\mathbf{\bm{\Omega}}^{di}} \right)' \bm{\omega}^{**} = X \bm{\omega}^{**}
\end{equation*}
The same result extends to the multivariate case in which the dependent
is a function of more than one latent predictor because the derivations
apply after replacing the vectors \(\bm{\omega}^*\) and \(\bm{\omega}^{**}\) with matrices
\(\mathbf{\bm{\Omega}}^*\) and \(\mathbf{\bm{\Omega}}^{**}\) where each column is a weight vector for one
latent factor \(\tilde{x}^*_j\).
\end{proof}
\section{LARX: Derivation of the Coefficient and Weight Vectors}
\label{derivation_w_o_p_g}
This section covers the derivations of the sample estimates \(\bm{\mathbf{\hat{w}}}\),
\(\bm{\hat{\omega}}\), \(\bm{\hat{\phi}}\) and \(\bm{\hat{\beta}}\) for the coefficient vectors \(\bm{\mathbf{w}}\), \(\bm{\omega}\),
\(\bm{\phi}\) and \(\bm{\beta}\), respectively. The customary ``hat'' superscripts are
omitted from here on out because they are implied everywhere. The
derivations for the Lagrange multipliers are presented in
\ref{derivation_ry_rl} and \ref{derivation_tlx_tlp}.
First, we note that the Lagrangian optimisation problem (\ref{eqn:clarx_lagr})
is convex, which means that the solution is obtained by setting various
partial derivatives to zero. Second, we note that the properties of the
Kronecker product and the block-wise Kronecker product allow us to
factorise \(\left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)\) and \(\left( {\bm{\beta} \odot \bm{\omega}} \right)\) for compatibility with traditional matrix
calculus:
\begin{subequations}
\begin{align}
& \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) = \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \bm{\phi} = \left( {\bm{\phi} \otimes I_n} \right) \bm{\mathbf{w}} \label{eqn:factorisation_p_w} & \\
& \left( {\bm{\beta} \odot \bm{\omega}} \right) = \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \bm{\beta} = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega} \label{eqn:factorisation_b_o} &
\end{align}
\end{subequations}
Here, \(I_a\) represents an identity matrix of size \(a\) (where \(a\)
is a scalar value), while \(I_{\mathbf{\bm{a}}}\) represents an identity matrix with
the same size and row block structure as vector \(\mathbf{\bm{a}}\) (vectors are
conventionally represented by bold letters). The property of the
Kronecker product relevant for (\ref{eqn:factorisation_p_w}) is well
established but reiterated for completeness in equations
(\ref{eqn:kron_factorisation_matrices}-\ref{eqn:kron_factorisation_vectors}).
The factorisation of the block-wise Kronecker product used for
(\ref{eqn:factorisation_b_o}) is proved in \ref{block_kron_factorisation}.
Taking \(\mathbf{\bm{\Sigma}}_{BA}\) to denote the transpose of \(\mathbf{\bm{\Sigma}}_{AB}\), the
solution for \(\bm{\beta}\) becomes:
\begin{align} \label{eqn:clarx_solution_beta}
& \frac{\partial}{\partial \, \bm{\beta}} \, \mathcal{L} = 2 \left[
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \bm{\beta} -
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} +
\left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) \right] \nonumber \\[10pt]
& \bm{\beta} = \left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \right]^{-1}
\left[ \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} - \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right)' \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) \right]
\end{align}
Similarly, recalling that \(\left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) = \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \bm{\phi}\), the solution for \(\bm{\phi}\) is:
\begin{align} \label{eqn:clarx_solution_phi}
& \frac{\partial}{\partial \, \bm{\phi}} \, \mathcal{L} = 2 \left[
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \bm{\phi} -
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} +
\left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right] \nonumber \\[10pt]
& \bm{\phi} = \left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right) \right]^{-1}
\left[ \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AY} \bm{\mathbf{w}} - \left( {I_{V_a} \otimes \bm{\mathbf{w}}} \right)' \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right) \right]
\end{align}
For the dependent weight vector \(\bm{\mathbf{w}}\), we recall that \(\left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right) = \left( {\bm{\phi} \otimes I_n} \right) \bm{\mathbf{w}}\).
The partial derivative of \(\mathcal{L}\) with respect to \(\bm{\mathbf{w}}\) is then given
by:
\begin{equation*}
\begin{split}
\frac{\partial}{\partial \, \bm{\mathbf{w}}} \, \mathcal{L}
& = 2 \left[ \left( 1 + \lambda_{y} \right) \mathbf{\bm{\Sigma}}_{Y} +
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} -
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right)
\right] \bm{\mathbf{w}} \\
& - 2 \left[
\mathbf{\bm{\Sigma}}_{YX} \left( {\bm{\beta} \odot \bm{\omega}} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AX} \left( {\bm{\beta} \odot \bm{\omega}} \right)
\right] + \lambda_l \1_n
\end{split}
\end{equation*}
Note that \(\frac{\partial}{\partial \, \bm{\mathbf{w}}} \, \bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) \bm{\mathbf{w}}\)
resolves to \(\left[ \left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} + \mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}}\) because
the quadratic form is not symmetric. Setting the partial derivative to
zero and expressing in terms of \(\bm{\mathbf{w}}\), we get:
\begin{equation} \label{eqn:clarx_derivation_w}
\begin{split}
(1 + \lambda_{y} ) \, \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} & =
\left[
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right) \right] \, \bm{\mathbf{w}} \\[10pt]
& + \left[
\mathbf{\bm{\Sigma}}_{YX} -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AX}
\right] \left( {\bm{\beta} \odot \bm{\omega}} \right) - \frac{\lambda_l}{2} \1_n
\end{split}
\end{equation}
For ease of notation, define:
\begin{align*}
\underset{1 \times n}{\bm{\mathbf{v}}_1} = & \left[
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AY} +
\mathbf{\bm{\Sigma}}_{YA} \left( {\bm{\phi} \otimes I_n} \right) -
\left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{A} \left( {\bm{\phi} \otimes I_n} \right)
\right] \bm{\mathbf{w}} \\[10pt]
\underset{1 \times n}{\bm{\mathbf{v}}_2} = & \left[\mathbf{\bm{\Sigma}}_{YX} - \left( {\bm{\phi} \otimes I_n} \right)' \mathbf{\bm{\Sigma}}_{AX} \right] \left( {\bm{\beta} \odot \bm{\omega}} \right)
\end{align*}
Equation (\ref{eqn:clarx_derivation_w}) then becomes:
\begin{equation} \label{eqn:clarx_derivation_w_short}
(1 + \lambda_{y} ) \, \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 - \frac{\lambda_l}{2} \1_n
\end{equation}
Setting \(\rho_{y} = (1 + \lambda_{y})\), \(\rho_l = \frac{\lambda_l}{2}\) and
pre-multiplying both sides by \(\frac{1}{\rho_{y}} \, \mathbf{\bm{\Sigma}}_{Y}^{-1}\) we
get:
\begin{equation} \label{eqn:clarx_solution_w}
\bm{\mathbf{w}} = \frac{1}{\rho_{y}} \left[ \mathbf{\bm{\Sigma}}_{Y}^{-1} \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right)
- \rho_l \mathbf{\bm{\Sigma}}_{Y}^{-1} \1_n \right]
\end{equation}
For the weight vector \(\bm{\omega}\), we recall that \(\left( {\bm{\beta} \odot \bm{\omega}} \right) = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega}\).
Furthermore, we can refactor the constraint terms for compatibility with
traditional matrix calculus using the properties of the block-wise direct
sum operator. For the portfolio constraints on \(\bm{\omega}_j\) we can apply
Proposition \ref{prop:bds_vs_vectors} to obtain:
\begin{equation*}
\left( \bm{\omega}^{\oplus} \right)' \1_{\bm{\omega}} = \left( \t1od \right)' \bm{\omega}
\end{equation*}
For the variance constraints on \(\tilde{x}_j\), we note three things. First
of all, because \(\bm{\mathbf{u}}\) has the same length and row block structure as
\(\bm{\beta}\), the term \(\left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)\) can be factorised in the same way as \(\left( {\bm{\beta} \odot \bm{\omega}} \right)\),
namely:
\begin{equation*}
\left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega} = \left( {I_{\bm{\beta}} \odot \bm{\omega}} \right) \bm{\mathbf{u}}
\end{equation*}
Second, by applying Proposition \ref{prop:bds_vs_transpose} we have:
\begin{equation*}
\left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)' \right]^{\oplus} = \left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)^{\oplus} \right]'
\end{equation*}
Third, the operation \(\left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)\) can be expressed as a left-multiplication
over the sequence of blocks in \(\bm{\omega}\), i.e., \(\left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega} = \left[
\left< {\left( {\bm{\mathbf{u}}_j \otimes \left< {I_{\bm{\omega}}} \right>_j} \right) \bm{\omega}_j | 1 \leq j \leq K} \right>q j \leq K}
\right]_v\), which means that according to Proposition \ref{prop:bds_vs_mmul}
we have:
\begin{equation*}
\left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)^{\oplus} \equiv \left[ \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega} \right]^{\oplus} = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega}^{\oplus}
\end{equation*}
Putting these transformations together, we can rewrite the block-wise
quadratic form as:
\begin{equation*}
\begin{split}
\left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)' \right]^{\oplus} & \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)
= \left[ \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)^{\oplus} \right]' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)
= \left\{ \left[ \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega} \right]^{\oplus} \right\}' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right) \\
& = \left[ \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega}^{\oplus} \right]' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot \bm{\omega}} \right)
= \left( { \bm{\omega}^{\oplus} } \right)' \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_x \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega}
\end{split}
\end{equation*}
Applying these transformations, the partial derivative of
(\ref{eqn:clarx_lagr}) with respect to \(\bm{\omega}\) is:
\begin{equation*}
\frac{\partial}{\partial \, \bm{\omega} } \, \mathcal{L}
= 2 \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \left[
\mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) \bm{\omega} -
\mathbf{\bm{\Sigma}}_{XY} \bm{\mathbf{w}} +
\mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes \bm{\mathbf{w}}} \right)
\right] + 2 \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x + \t1od \bm{\mathbf{\lambda}}_p
\end{equation*}
where \(\mathbf{\bm{M}}_2 = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)\). Furthermore, Propositions
\ref{prop:bds_vs_odot} and \ref{prop:block_kron_factorisation} prove that:
\begin{equation*}
\bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x = \left( {\bm{\omega} \odot I_K} \right) \bm{\mathbf{\lambda}}_x = \left( {\bm{\mathbf{\lambda}}_x \odot I_{\bm{\omega}}} \right) \bm{\omega}
\end{equation*}
Re-arranging for \(\bm{\omega}\), we get:
\begin{equation} \label{eqn:clarx_solution_o}
\bm{\omega} = \left[ \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}_{X} \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right) + \mathbf{\bm{M}}_2 \left( {\bm{\mathbf{\lambda}}_x \odot I_{\bm{\omega}}} \right) \right]^{-1}
\left[ \bm{\mathbf{v}}_3 - \frac{1}{2} \t1od \bm{\mathbf{\lambda}}_p \right]
\end{equation}
with \(\bm{\mathbf{v}}_3 = \left( {\bm{\beta} \odot I_{\bm{\omega}}} \right)' \left[ \mathbf{\bm{\Sigma}}_{XY} - \mathbf{\bm{\Sigma}}_{XA} \left( {\bm{\phi} \otimes I_n} \right) \right] \bm{\mathbf{w}}\).
\section{LARX: Derivation of the Lagrange Multipliers for the Dependent}
\label{derivation_ry_rl}
Recalling that \(\rho_{y} = (1 + \lambda_{y})\) and \(\rho_l = \frac{\lambda_l}{2}\),
rewrite equation (\ref{eqn:clarx_derivation_w_short}) as:
\begin{flalign} \label{eqn:clarx_ry_1}
& \rho_{y} \, \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 - \rho_l \1_n &
\end{flalign}
To solve for \(\rho_l\), pre-multiply both sides of (\ref{eqn:clarx_ry_1}) by
\(\1_n'\) and re-arrange:
\begin{flalign*}
& \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} = \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - n \rho_l & \\[10pt]
& n \rho_l = \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} &
\end{flalign*}
Dividing both sides by \(n\) produces:
\begin{flalign} \label{eqn:clarx_solution_rl}
& \rho_l = \frac{1}{n} \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \frac{1}{n} \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} &
\end{flalign}
We can find an alternative solution for \(\rho_l\) by pre-multiplying both
sides of (\ref{eqn:clarx_ry_1}) with \(\bm{\mathbf{w}}'\) and recalling \(\bm{\mathbf{w}}' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}} =
\s2_y\) and \(\bm{\mathbf{w}}' \1_n = l_y\):
\begin{flalign*}
& \rho_{y} \s2_y = \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \rho_l l_y & \\[10pt]
& \rho_l l_y = \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \rho_{y} \s2_y &
\end{flalign*}
dividing both sides by \(l_y\) we get:
\begin{flalign} \label{eqn:clarx_solution_rl_alt}
& \rho_l = \frac{1}{l_y} \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) -
\frac{\s2_y}{l_y}\rho_{y} &
\end{flalign}
This alternative solution is not very practical because it does not
allow for the case of \(l_y = 0\). However, we can use it in conjunction
with (\ref{eqn:clarx_solution_rl}) to eliminate \(\rho_l\) and solve for
\(\rho_y\):
\begin{flalign*}
& \frac{1}{n} \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \frac{1}{n} \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}
= \frac{1}{l_y} \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \frac{\s2_y}{l_y}\rho_{y} & \\[10pt]
& \frac{\s2_y}{l_y}\rho_{y} - \frac{1}{n} \rho_{y} \, \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}
= \frac{1}{l_y} \bm{\mathbf{w}}' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) - \frac{1}{n} \1_n' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) & \\[10pt]
& \rho_{y} \frac{n \s2_y - l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}}{n l_y}
= \frac{\left( n \bm{\mathbf{w}} - l_y \1_n \right)' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) }{n l_y} &
\end{flalign*}
Rearranging for \(\rho_{y}\) we get:
\begin{flalign} \label{eqn:clarx_solution_ry}
& \rho_{y} = \frac{ \left( n \bm{\mathbf{w}} - l_y \1_n \right)' \left( \bm{\mathbf{v}}_1 + \bm{\mathbf{v}}_2 \right) }
{n \s2_y - l_y \1_n' \mathbf{\bm{\Sigma}}_{Y} \bm{\mathbf{w}}} &
\end{flalign}
\section{LARX: Derivation of the Lagrange Multipliers for the Explanatory}
\label{derivation_tlx_tlp}
Start from the first-order condition for \(\bm{\omega}\). Expressing in terms of
\(\bm{\mathbf{\lambda}}_p\), we get:
\begin{flalign} \label{eqn:deriv_tlpx_1}
& \t1od \bm{\mathbf{\lambda}}_p = 2 \bm{\mathbf{v}}_3 - 2 \bm{\mathbf{v}}_4 - 2 \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x &
\end{flalign}
Define the following shorthand notations for convenience:
\begin{flalign*}
& \underset{K \times K}{\bm{\mathbf{\Theta}}} = \operatorname{diag}(\ts2_x) = \begin{pmatrix}
\s2_{x,1}, & 0, & \cdots & 0 \\
0, & \s2_{x,2}, & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0, & 0, & \cdots & \s2_{x,K}
\end{pmatrix}, \underset{K \times K}{\mathbf{\bm{L}}} = \operatorname{diag}(\bm{l}_p) = \begin{pmatrix}
l_{p,1}, & 0, & \cdots & 0 \\
0, & l_{p,2}, & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0, & 0, & \cdots & l_{p,K}
\end{pmatrix}, & \\[10pt]
& \underset{K \times K}{\mathbf{\bm{M}}_1} = \left( \t1od \right)' \t1od = \begin{pmatrix}
m_1, & 0, & \cdots & 0 \\
0, & m_2, & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0, & 0, & \cdots & m_K
\end{pmatrix}, \quad \underset{K \times K}{\mathbf{\bm{M}}_2} = \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right)' \mathbf{\bm{\Sigma}}^{d}_X \left( {\bm{\mathbf{u}} \odot I_{\bm{\omega}}} \right) &
\end{flalign*}
A system of two vector equations, each with \(K\) rows and \(K\)
unknowns, is produced by pre-multiplying the first-order condition for
\(\bm{\omega}\) with \(\left( {\bm{\omega}^{\oplus}} \right)'\) and \(\left( {\t1od} \right)'\), respectively:
\begin{subnumcases}{\hspace{-15em}}
\mathbf{\bm{L}} \bm{\mathbf{\lambda}}_p = 2 \left( \bm{\omega}^{\oplus} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - 2 \bm{\mathbf{\Theta}} \bm{\mathbf{\lambda}}_x
\label{eqn:deriv_tlpx_cond_1} \\[10pt]
\mathbf{\bm{M}}_1 \bm{\mathbf{\lambda}}_p = 2 \left( \t1od \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right)
- 2 \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x
\label{eqn:deriv_tlpx_cond_2}
\end{subnumcases}
Pre-multiplying \ref{eqn:deriv_tlpx_cond_1} and
\ref{eqn:deriv_tlpx_cond_2} by \(\mathbf{\bm{L}}^{-1}\) and \(\mathbf{\bm{M}}_1^{-1}\),
respectively, we get:
\begin{subnumcases}{\hspace{-12em}}
\bm{\mathbf{\lambda}}_p = 2 \mathbf{\bm{L}}^{-1} \left( \bm{\omega}^{\oplus} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - 2 \mathbf{\bm{L}}^{-1} \bm{\mathbf{\Theta}} \bm{\mathbf{\lambda}}_x &
\label{eqn:clarx_solution_tlp_a} \\[10pt]
\bm{\mathbf{\lambda}}_p = 2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) -
2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x &
\label{eqn:clarx_solution_tlp_b}
\end{subnumcases}
The solution for \(\bm{\mathbf{\lambda}}_p\) given by \ref{eqn:clarx_solution_tlp_b} is
more practical because \ref{eqn:clarx_solution_tlp_a} does not allow for
the case of zero-sum weights (zeros on the diagonal of \(\mathbf{\bm{L}}\) would
make it uninvertible).
The solution for \(\bm{\mathbf{\lambda}}_x\) can be derived by setting the right-hand side
of \ref{eqn:clarx_solution_tlp_a} equal to the right-hand side of
\ref{eqn:clarx_solution_tlp_b}:
\begin{flalign*}
& 2 \mathbf{\bm{L}}^{-1} \left( \bm{\omega}^{\oplus} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - 2 \mathbf{\bm{L}}^{-1} \bm{\mathbf{\Theta}} \bm{\mathbf{\lambda}}_x =
2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) -
2 \mathbf{\bm{M}}_1^{-1} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x &
\end{flalign*}
We can pre-multiply both sides by \(\mathbf{\bm{M}}_1 \mathbf{\bm{L}}\) to avoid problems with
inverting \(\mathbf{\bm{L}}\) in the presence of zero-sum weight constraints (note
that \(\mathbf{\bm{M}}_1 \mathbf{\bm{L}} = \mathbf{\bm{L}} \mathbf{\bm{M}}_1\) because both are diagonal):
\begin{flalign*}
& 2 \mathbf{\bm{M}}_1 \left( \bm{\omega}^{\oplus} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - 2 \mathbf{\bm{M}}_1 \bm{\mathbf{\Theta}} \bm{\mathbf{\lambda}}_x =
2 \mathbf{\bm{L}} \left( \t1od \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right) - 2 \mathbf{\bm{L}} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \bm{\mathbf{\lambda}}_x &
\end{flalign*}
Re-arranging for \(\bm{\mathbf{\lambda}}_x\) we get:
\begin{equation} \label{eqn:clarx_solution_tlx}
\bm{\mathbf{\lambda}}_x = \left[ \mathbf{\bm{M}}_1 \bm{\mathbf{\Theta}} - \mathbf{\bm{L}} \left( \t1od \right)' \mathbf{\bm{M}}_2 \bm{\omega}^{\oplus} \right]^{-1}
\left( \bm{\omega}^{\oplus} \mathbf{\bm{M}}_1 - \t1od \mathbf{\bm{L}} \right)' \left( \bm{\mathbf{v}}_3 - \bm{\mathbf{v}}_4 \right)
\end{equation}
\section{Out-of-Sample Forecasting Performance of the Alternative Forecasting Models}
\label{empirical_results_rev}
\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{NEW_fig-fc_ols_vs_clarx_reversed.pdf}
\caption{Out-of-sample performance: Market returns as the dependent} \label{NEW_fig-fc_ols_vs_clarx_reversed}
\floatfoot{ \footnotesize
Note: Rolling out-of-sample forecasting performance of the
alternative forecasting models with stock market returns as the
dependent variable regressed on lags zero to 2 of real economic
activity and 3 autoregressive lags as specified by equations
(\ref{eqn:regression_models_rev})-(\ref{eqn:regression_all_latent_rev}).
The target variance for the real activity SDI is set to the
full-sample variance of quarterly US real GDP growth (approx.
0.000124). The variance of the market SDI is set to the
full-sample variance of the S\&P 500 quarterly log-return
(approx. 0.006348).
}
\end{figure}
\begin{subequations} \label{eqn:regression_models_rev}
Alternative baseline OLS/ARX model: The S\&P 500 \(r\) vs real GDP \(g\):
\begin{equation} \tag{\ref*{eqn:regression_models_rev}}
r_t = c + \sum_{\tau = 1}^3 \phi_{t-\tau} r_{t-\tau} +
\sum_{\tau = 0}^2 \beta_{t-\tau} g_{t-\tau} + e
\end{equation}
Alternative LARX model: Latent explanatory: The S\&P 500 vs the real activity SDI
\(\tilde{g}\):
\begin{equation}\label{eqn:regression_r_nonlatent_rev}
r_t = c + \sum_{\tau = 1}^3 \phi_{t-\tau} r_{t-\tau} +
\sum_{\tau = 0}^2 \beta_{t-\tau} \tilde{g}_{t-\tau} + e
\end{equation}
LARX model 2b: Latent dependent: Market SDI \(\tilde{r}\) vs real GDP:
\begin{equation}\label{eqn:regression_g_nonlatent_rev}
\tilde{r}_t = c + \sum_{\tau = 1}^3 \phi_{t-\tau} \tilde{r}_{t-\tau} +
\sum_{\tau = 0}^2 \beta_{t-\tau} g_{t-\tau} + e
\end{equation}
LARX model 2c: Both latent: Market SDI vs the real activity SDI:
\begin{equation}\label{eqn:regression_all_latent_rev}
\tilde{r}_t = c + \sum_{\tau = 1}^3 \phi_{t-\tau} \tilde{r}_{t-\tau} +
\sum_{\tau = 0}^2 \beta_{t-\tau} \tilde{g}_{t-\tau} + e
\end{equation}
\end{subequations}
\section{Data reference}
\label{sec:org483dd8c}
\begin{table}[H]
\centering
\caption{Data series used in the empirical study}
\label{table:data_reference}
\begin{adjustbox}{width=\textwidth, nofloat=table}
\begin{threeparttable}
\begin{tabular}{lllll}
\toprule
Dataset & Source & Ticker\tnote{1} & Frequency & History start \\
\midrule
Real GDP & U.S. Bureau of Economic Analysis & GDPC1 & Quarterly & 1947Q1 \\
\quad Personal Consumption Expenditure (Cons.) & U.S. Bureau of Economic Analysis & PCECC96 & Quarterly & 1947Q1 \\
\quad Gross Private Domestic Investment (Inv.) & U.S. Bureau of Economic Analysis & GPDIC1 & Quarterly & 1947Q1 \\
\quad Government Consumption and Investment (Govt.) & U.S. Bureau of Economic Analysis & GCEC1 & Quarterly & 1947Q1 \\
\quad Exports of Goods and Services (Exp.) & U.S. Bureau of Economic Analysis & EXPGSC1 & Quarterly & 1947Q1 \\
\quad Imports of Goods and Services (Imp.) & U.S. Bureau of Economic Analysis & IMPGSC1 & Quarterly & 1947Q1 \\
\hdashline S\&P 500 & Investing.com & US500 & Monthly & 1989-10 \\
\quad Energy (En.) & Investing.com & SPNY & Monthly & 1989-10 \\
\quad Materials (Mat.) & Investing.com & SPLRCM & Monthly & 1989-10 \\
\quad Industrials (Ind.) & Investing.com & SPLRCI & Monthly & 1989-10 \\
\quad Financials (Fin.) & Investing.com & SPSY & Monthly & 1989-10 \\
\quad Healthcare (HC) & Investing.com & SPXHC & Monthly & 1989-10 \\
\quad Consumer Discretionary (Disc.) & Investing.com & SPLRCD & Monthly & 1989-10 \\
\quad Consumer Staples (Stapl.) & Investing.com & SPLRCS & Monthly & 1989-10 \\
\quad Communication (Telco) & Investing.com & SPLRCL & Monthly & 1989-10 \\
\quad Technology (IT) & Investing.com & SPLRCT & Monthly & 1989-10 \\
\quad Utilities (Util.) & Investing.com & SPLRCU & Monthly & 1989-10 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\item[1] Data for U.S. GDP and its individual expenditure components was
retrieved from the St. Louis Federal Reserve economic database (FRED)
on 20 October 2025. The corresponding tickers are identifiers for the
FRED database. Data for the S\&P 500 and its sector sub-indices was
retrieved directly from Investing.com on 20 October 2025 using the
tickers above.
\end{tablenotes}
\end{threeparttable}
\end{adjustbox}
\end{table}
\end{appendix}
\bibliographystyle{unsrtnat}
\bibliography{/home/daniil/Research/bibliography/phd1}