Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
1818 factor-augmented tree ensembles
\begingroup
\footnote{ I thank Matteo Barigozzi and Kostas Kalogeropoulos for their valuable suggestions and supervision; Serena Lariccia, Lucrezia Reichlin, Esther Ruiz Ortega, Veronica Veggente, Qiwei Yao, the 2022 IMS Annual Meeting in Probability and Statistics and the 2023 Italian Congress of Econometrics and Empirical Economics participants for their helpful comments on a preliminary draft of this article.}
\addtocounter{footnote}{-1}
\endgroup
\begingroup
\footnote{ Disclaimer: Filippo Pellegrino is funded by J.P. Morgan Chase & Co under a J.P. Morgan A.I. Research Award. Any views or opinions expressed herein are solely those of the authors listed, and may differ from the views and opinions expressed by J.P. Morgan Chase & Co. or its affiliates. This material is not a product of the Research Department of J.P. Morgan Securities LLC. This material does not constitute a solicitation or offer in any jurisdiction.}
\addtocounter{footnote}{-1}
\endgroup
\thispagestyle{empty}
\null
center[center omitted — 82 chars of source]
abstractThis manuscript proposes to extend the information set of time-series regression trees with latent stationary factors extracted via state-space methods. In doing so, this approach generalises time-series regression trees on two dimensions. First, it allows to handle predictors that exhibit measurement error, non-stationary trends, seasonality and/or irregularities such as missing observations. Second, it gives a transparent way for using domain-specific theory to inform time-series regression trees. Empirically, ensembles of these factor-augmented trees provide a reliable approach for macro-finance problems. This article highlights it focussing on the lead-lag effect between equity volatility and the business cycle in the United States.
{Keywords: Ensemble learning, Factor models, Macro-finance, State-space models.} \\
{JEL: C32, C58, E32, E44, G17.}
\null
\setcounter{page}{1}
\let\footnote=\endnote
\patchcmd{\enoteformat}{1.8em}{0pt}
Introduction
In time series, the simplicity of regression trees morgan1963problems, breiman1984classification, quinlan1986induction comes at a cost: irregularities, complicated periodic patterns and non-stationary trends cannot be explicitly modelled, and this is unfortunate given that many real-world examples are subject to them. This is especially true in macroeconomics and finance as shown by a broad body of literature. This includes articles focussed on the output gap, the natural rate of interest and other structural components hamilton2016equilibrium, holston2017measuring, jarocinski2018inflation, fries2018national, del2019global, barigozzi2021measuring, hasenzagl2022a, hasenzagl2022b, non-stationary forecasting settings and impulse response functions of economic data in levels sims1980macroeconomics, doan1984forecasting, litterman1986forecasting, banbura2010large, giannone2015prior, barigozzi2021large, and problems with a prevalence of missing observations comprising nowcasting giannone2008nowcasting, banbura2014maximum, cimadomo2022nowcasting, brown2022nowcasting.
Following, in spirit, harvey1998messy, this paper proposes to pre-process problematic predictors using state-space representations general enough to deal with all these complexities at once. This operation can be thought as an automated feature engineering process that extracts stationary patterns hidden across multiple predictors, while handling problematic data characteristics. Besides, when the state-space representation is compatible with domain-specific theory, this becomes a transparent way for extracting signals with structural interpretation. The stationary common components recovered from the data, referred hereinbelow as stationary dynamic factors, are then employed as regular predictors for standard time-series regression trees. This manuscript calls them factor-augmented regression trees to stress their dependence on latent components.
For this article, I have built on a broad body of theoretical research on time series. Indeed, factor models originated in psychometrics lawley1962factor as a dimensionality reduction technique. They were later generalised to take into account the autocorrelation structure of time series with the work of geweke1977DFM on dynamic factor models. Over time, these methodologies have been further developed within the state-space literature pioneered by harvey1985trends to be compatible with data exhibiting peculiar patterns (e.g., non-stationary trends, seasonality) and missing observations. Relevant improvements include the generalized dynamic factor model forni2000generalized, forni2001generalized, forni2005generalized, forni2009opening, the factor-augmented vector autoregression bernanke2005measuring, quasi maximum likelihood estimation methods doz2012quasi, barigozzi2020quasi, the non-stationary dynamic factor model barigozzi2021large and functional analysis extensions li2020nonlinear.
As for standard regression trees, the factor-augmented version can suffer from over-fitting. Tree ensembles are an efficient way to reduce it without having to use complex vectors of hyperparameters. In order to do that, these methods generally fit a series of regression trees on a range of data subsamples and return aggregate forecasts. This article constructs the ensembles following breiman1996bagging. These factor-augmented ensembles are similar to the rotation forest proposed in rodriguez2006rotation and pardo2013rotation, but they take into account the autocorrelation structure in the data when estimating the factors and have the higher flexibility embedded in state-space modelling.
Factor-augmented regression trees and their ensembles are also strongly motivated by empirical results in economics and finance. Recent literature on semi-structural models, including hasenzagl2022a, hasenzagl2022b, proposed to enrich statistical trend-cycle decompositions by using a minimal set of economic-driven restrictions. The main advantage of these semi-structural models is that they are able to extract unobserved cyclical and persistent components with economic interpretation, while allowing the data to speak. However, it is often hard to determine reasonable restrictions for most high-dimensional problems. Indeed, in macroeconomics and finance, theory is unclear on the exact dynamics of classical aggregate variables and disaggregated indicators. Also, the literature is not mature enough to understand the precise drivers of new data (e.g., Google searches). Factor-augmented regression trees and their ensembles can be seen as a bridge between the output of small-dimensional semi-structural models (i.e., interpretable cyclical unobserved components) and time series that are not entirely understood from the theoretical standpoint and/or exhibit non-linear dynamics.
As a result, these factor-augmented models are well-suited for tackling macro-finance problems. Indeed, it is often unclear how to model specific drivers and non-linear links between macroeconomic and financial data. For instance, forecasting the yield curve while handling the zero lower bound kim2012term, swanson2014measuring, bauer2016monetary, bauer2020interest, studying exchange rate predictability meese1983empirical, engel2005exchange, rossi2013exchange, or the links between equity and economic aggregates nelson1976inflation, fama1977asset, cooper2009time, campbell2009stock. This article focusses on the latter and studies the lead-lag effect between equity volatility and a factor representing the business cycle in the United States. Results show that this is is a strong predictor of equity volatility, especially for small cap stocks. I find this in agreement with a large range of papers such as gertler1994monetary, sharpe1994financial and bernanke1996financial whereas it is argued that small firms are more vulnerable to recessions due to the lack of financing options.
Methodology
Regression trees
This subsection describes the population model implied by standard regression trees and their most common estimation method breiman1984classification.
assumption[Data]
Let $T \in \mathbb{N}$ and $n \in \mathbb{N}_{0}$. Assume that $Y_{t}$ and $Z_{j,t}$ are finite realisations of some real-valued mean-stationary stochastic processes observed at the time periods $t=1, \ldots, T$ and with $j=1, \ldots, n$.
assumption[Predictors of standard regression trees]
Let $\boldsymbol{\mathbf{X}}_t \vcentcolon= (Y_t \; Z_{1,t} \; \ldots \; Z_{n,t})'$ be $m \times 1$ dimensional and defined for any point in time $t \in \mathbb{Z}$.
remarkThroughout the manuscript, the dependency on $n$ and $T$ is highlighted only when strictly necessary to ease the reading experience. Furthermore, specific realisations at some integer point in time $t$ and their general value in the underlying stochastic processes are denoted with the same symbols. However, it should be clear from the context whether the manuscript is referring to the first or second category.
This article describes the regression trees as non-linear forecasting models for $Y_t$ based on the information included in $\boldsymbol{\mathbf{X}}_{t-1}, \ldots, \boldsymbol{\mathbf{X}}_{t-p}$, for some number of lags $p \in \mathbb{N}$. Without loss of generality, this manuscript focusses on one-step ahead forecasts. Long-run predictions can be generated by computing direct forecasts marcellino2006comparison.
assumption[Lags]
Let $0 < p \ll T-1$.
assumption[Regression tree model]
In a regression tree setting,
\begin{align*}
Y_{t} = \sum_{i=1}^{|\mathscr{F}|} b_{i} \mathbb{I} \left\{ ( \boldsymbol{\mathbf{X}}_{t-1} \, \ldots \, \boldsymbol{\mathbf{X}}_{t-p} ) \in \mathscr{F}_{i} \right\} + \epsilon_t,
\end{align*}
whereas $\mathscr{F}$ is an indexed family of disjoint sets of matrices, every $b_{i}$ is a finite constant, $\epsilon_t
\mathrel{\overset{i.i.d.}{\scalebox{1.5}[1]{$\sim$}}}
\left(0, \sigma^2\right)$ with $\sigma > 0$ and finite, for any integer $t$ and $i=1, \ldots, |\mathscr{F}|$. Besides, regression trees assume that
\begin{align*}
\mathbb{E} (Y_{t} | \boldsymbol{\mathbf{X}}_{t-1}, \ldots, \boldsymbol{\mathbf{X}}_{t-p}, \boldsymbol{\mathbf{b}}, \sigma, \mathscr{F}) = \sum_{i=1}^{|\mathscr{F}|} b_{i} \mathbb{I} \left\{ ( \boldsymbol{\mathbf{X}}_{t-1} \, \ldots \, \boldsymbol{\mathbf{X}}_{t-p} ) \in \mathscr{F}_{i} \right\},
\end{align*}
for any integer $t$.
Regression trees estimate $\boldsymbol{\mathbf{b}}$, $\sigma$ and $\mathscr{F}$ recursively partitioning the predictor space to find the best fit. There are several modelling choices to take when performing this operation. This article follows common practice by focussing on binary partitions and using the CART algorithm breiman1984classification. At its first iteration, this estimation method looks for the best possible way to split the predictor space into two regions. This assessment is performed fitting a constant model in each region and minimising the mean square forecast error. Moreover, the splits are computed by inspecting, in turn, each covariate separately. The algorithm iteratively repeats the same operation for each of the resulting regions until some stopping criteria is reached. This manuscript uses \href{https://github.com/JuliaAI/DecisionTree.jl}{DecisionTree.jl} to implement it and refers to the estimated parameters with $\boldsymbol{\mathbf{\hat{\theta}}}(\boldsymbol{\mathbf{\gamma}})$ where $\boldsymbol{\mathbf{\gamma}}$ is a vector of hyperparameters.
Factor-augmented regression trees
This subsection introduces the factor-augmented regression trees: a version of the model in (ref) able to handle predictors with irregularities such as structural breaks and missing observations, intricate periodic patterns and non-stationary trends. In order to deal with these complexities, this subsection introduces a series of changes to the model and estimation algorithm.
Factor-augmented regression trees allow for these complexities in the data redefining $\boldsymbol{\mathbf{Z}}_{t}$ and $\boldsymbol{\mathbf{X}}_{t}$ as detailed in \crefrange{ch2:assumption:data_predictors_f1}{ch2:assumption:data_predictors_f3}.
assumption[State-space representation: data]
Assume that $Z_{i,t}$ is finite realisations of some real-valued stochastic process observed at the time periods in the set $\mathscr{T}_i \subseteq \{t : t \in \mathbb{Z},\, 1 \leq t \leq T\}$ for every $i=1, \ldots, n$.
assumption[State-space representation: structure]
Let $\boldsymbol{\mathbf{Z}}_{t}$ be a $n \times 1$ real random vector that allows the state-space representation
\begin{align*}
Z_{i,t} &= g_{i,t}(\boldsymbol{\mathbf{\Phi}}_t, \xi_{i,t}) ,\\
\boldsymbol{\mathbf{\Phi}}_t &= \boldsymbol{\mathbf{h}}_t(\boldsymbol{\mathbf{\Phi}}_{t-1}, \boldsymbol{\mathbf{\zeta}}_t),
\end{align*}
where $g_{i,t}$ and $\boldsymbol{\mathbf{h}}_t$ are continuous and differentiable functions, $\boldsymbol{\mathbf{\Phi}}_t$ denotes a vector of $q > 0$ latent states, $\boldsymbol{\mathbf{\xi}}_t
\mathrel{\overset{i.i.d.}{\scalebox{1.5}[1]{$\sim$}}}
(\boldsymbol{\mathbf{0}}_{n \times 1}, \boldsymbol{\mathbf{R}}_t)$ and $\boldsymbol{\mathbf{\zeta}}_t
\mathrel{\overset{i.i.d.}{\scalebox{1.5}[1]{$\sim$}}}
(\boldsymbol{\mathbf{0}}_{q \times 1}, \boldsymbol{\mathbf{Q}}_t)$, for any integer $t$ and $i=1,\ldots, n$. Also, it is assumed that every $\boldsymbol{\mathbf{\Phi}}_t$ includes a $\bar{q} \times 1$ vector of stationary common factors $\boldsymbol{\mathbf{\phi}}_t$, with $0 < \bar{q} \ll n$ and $q \geq \bar{q}$. Since the observations start from the time period $t=1$, it is further assumed that $\boldsymbol{\mathbf{\Phi}}_1 = \boldsymbol{\mathbf{h}}_0(\boldsymbol{\mathbf{\zeta}}_0)$. This allows the evaluation of the state-space representation with the observed data.
remark[Non-stationarity]
Differently than with \crefrange{ch2:assumption:data_target}{ch2:assumption:data_predictors}, (ref) does not assume that the underlying stochastic process is stationary. As a result, (ref) is compatible with non-stationary trends and co-integrated relationships.
assumption[Predictors of factor-augmented trees]
Factor-augmented regression trees include the stationary common factors in the predictors (as is, transformed in a way that does not alter data ordering and preserves stationarity, or both). Formally, this is achieved including these common components in the predictor matrix $\boldsymbol{\mathbf{X}}_t$ jointly with $Y_t$ and updating $m$ accordingly.
remark[Information set]
Factor-augmented regression trees extend the information set of a tree autoregression for $Y_t$ with stationary common factors, while discarding idiosyncratic noise in the predictors and non-stationary trends, and handling data irregularities. The simplest case is when the predictor matrix is extended to include these stationary factors as they come out from the state-space. Formally, this is achieved by letting $\boldsymbol{\mathbf{X}}_t \vcentcolon= (Y_t \,\; \boldsymbol{\mathbf{\phi}}_t)$ be a $m \times 1$ vector of time series, with $m \vcentcolon= 1+\bar{q}$.
The structure of the model is exactly as described in (ref), but uses the newly defined predictor matrix and value for $m$. However, the estimation process is different and structured as a two-step method. In the first step, the state space in (ref) is estimated with any algorithm compatible with the data complexities described above, including, but not limited to, the EM dempster1977maximum, rubin1982algorithms, shumway1982approach, watson1983alternative, banbura2014maximum, barigozzi2020quasi, ECM meng1993maximum, pellegrino2020selecting and ECME algorithms liu1994ecme.\footnote{Bayesian techniques surveyed, for instance, in sarkka2013bayesian can also be used. In that case, $\boldsymbol{\mathbf{\phi}}_t$ would be a point estimate (e.g., mean or median) of the stationary dynamic factors distribution at time $t$.} In the second and final step, the predictor matrix is formed on the basis of the estimated states and the regression tree is trained with CART.
It is worth stressing that the main difference between factor-augmented regression trees and the individual base learners of rotation forests rodriguez2006rotation, pardo2013rotation lies in the technique used for reducing the dimensionality of the data. Instead of using Principal Component Analysis, factor-augmented regression tree models use a state-space. In doing so, this approach explicitly models the temporal factors dynamics\footnote{This is fundamentally the same difference between traditional and dynamic factor models barigozzi2020quasi.}, permits to pinpoint specific unobserved components and allows for data that exhibits peculiar patterns such as non-stationary trends. Note that factor-augmented regression trees could be extended to use selected idiosyncratic periodic patterns as additional predictors. This could be done by redefining $\boldsymbol{\mathbf{\phi}}$ into a vector of “selected cycles”, both common and idiosyncratic. However, this would increase the computational burden and, without limitations, the risk of generating spurious splits.
Tree ensembles
Tree ensembles are methods that combine multiple regression trees, in order to produce more efficient predictions than the individual base learners (i.e., the trees themselves).
For simplicity of illustration, this article focusses on ensemble averaging and, in particular, on bootstrap aggregating or bagging breiman1996bagging. This method obtains the increase in efficiency estimating a large number of regression trees on random data subsamples and combining their predictions taking a sample average. Intuitively, this reduces over-fitting since the base learners are not trained on the original data, but on random subsamples generated from it. The more heterogeneous and numerous the subsamples the better in terms of efficiency. This can be formalised following an approach equivalent to hastie2009elements.
This article follows common practice and uses the bootstrap version proposed in efron1983leisurely to generate the subsamples. This approach considers each covariate-response pair as a single datapoint and constructs data subsamples via independent bootstrap efron1979bootstrap, efron1979computers, efron1981nonparametric. In other words, it resamples covariate-response pairs from the original data to generate the subsamples. In particular, in the case of the factor-augmented trees this is done focussing on the factor-response pairs.
figure[figure omitted — 138 chars of source]
Monte Carlo study
Before delving into the empirical analysis, these factor-augmented ensembles are examined through the lens of the following Monte Carlo study.
I simulate data from
align*[align* omitted — 259 chars of source]
where $0 \leq w \leq 1$. The dynamics of the stationary cycle $\phi_{t}$ are set as in clark1987cyclical to simulate a component similar to the US business cycle. The coefficients linking the simulated target to the stationary cycle are such that $\rho_{1}
\mathrel{\overset{i.i.d.}{\scalebox{1.5}[1]{$\sim$}}}
U(-0.1, 0)$ and $\rho_{2}
\mathrel{\overset{i.i.d.}{\scalebox{1.5}[1]{$\sim$}}}
U(0, 0.1)$. For simplicity, $\rho_{3}=1$. The model allows for non-linear links between the target and cycle. These non-linearities are controlled by the scalar $w$. The constraint $\rho_{1} \leq \rho_{2}$ simulates data resembling procyclical asset returns that react differently to recessions than expansions.
I generate the model $1,000$ times for $w=0, 0.1, \ldots, 1$ and $T=100, 200$. For each replicate, I use the first half of the sample to regress the target onto the cycle and predict the remaining target data points out-of-sample. The regressions are performed both via OLS (as in standard factor regressions) and bagging (as in the factor-based ensembles proposed herein). Figure (ref) shows the out-of-sample mean squared error for all models, aggregated over the 1,000 simulations. The first row reports the results for the case in which the estimation and forecasts are based on the true cycle. The second row shows the equivalent results for the case in which the cycle is contaminated with random noise drawn from a normal distribution with the same variance as the cycle innovations themselves. The results indicate that traditional factor regressions are subpar, except when the data generating process is strongly linear. In fact, even with minor non-linearities, bagging produces more accurate out-of-sample predictions. More broadly, this suggests that factor-augmented ensembles are likely to outperform traditional approaches based on linear models when predicting data with suspected non-linearities.
table[table omitted — 1,956 chars of source]
Advantages for empirical macro-finance
Academic insights indicate that financial returns are linked to macroeconomic fundamentals. However, these relations are usually hard to find and may be subject to non-linearities.
In principle, a simple way to exploit this behaviour would be running a non-linear predictive regressions using the lagged business cycle as predictor and some function of a financial return of interest as a response. However, this is easier said than done. Indeed, the business cycle itself is an unobserved variable that reflects the cyclical co-movement between a series of non-stationary economic indicators (e.g., output, unemployment and inflation). Besides, the non-linear links with the financial returns have an unclear form and, thus, it is hard to model them in a parametric way.
Factor-augmented regression trees and their ensembles represent a simple approach to the problem, compatible with its complexities. In fact, the state-space in (ref) can be thought as a way for extracting the business cycle from a set of predictors and the regression tree as a model that does not require an a-priori parametrisation of the non-linear link between macroeconomic and financial data.
Empirical setting
This subsection illustrates the advantages of factor-augmented ensembles focussing on the United States. Specifically, it employs these techniques for predicting equity volatility -- measured as squared returns -- of the financial indices in (ref) as a function of lagged information on themselves and the business cycle.
In order to estimate the state of the economy in real time, this section uses a state-space representation similar, in spirit, to the one proposed in hasenzagl2022a, hasenzagl2022b. This modelling choice implies that each macroeconomic indicator in (ref) is considered as the sum of non-stationary trends and causal cycles, one of which can be interpreted as the US business cycle. The trends account for the persistence in the data and provide a view on a series of structural components such as the natural rate of unemployment and trend inflation. By linking together key variables such as the real personal consumption expenditures, unemployment rate and inflation through the business cycle, the model is compatible with economic relationships such as the Phillips curve and the Okun's law (interpreting consumption as a proxy for GDP). A complex lag structure in the coefficients associated with the business cycle allows to take into account frictions in the economy (for instance, in the labour market). Finally, the idiosyncratic cycles account for autocorrelation in the error terms (if any).
This approach can be formalised as follows.
assumption[State-space representation: trend-cycle model]
For any integer $t$, let $\boldsymbol{\mathbf{Z}}_{t}$ denote the macroeconomic indicators included in the first block of (ref). For simplicity, assume that the data in $\boldsymbol{\mathbf{Z}}_{t}$ is the order reported in the table and standardised such that each $i$-th series is divided for a given scaling factor $\eta_{i}$, for $i=1, \ldots, n$ with $n=9$. Hence, let
\begin{align*}
\left( \begin{array}{c} Z_{1,t} \\ Z_{2,t} \\ Z_{3,t} \\ Z_{4,t} \\ Z_{5,t} \\ Z_{6,t} \\ Z_{7,t} \\ Z_{8,t} \\ Z_{9,t} \end{array} \right) &= \left( \begin{array}{c} \tau_{1,t} \\ \tau_{2,t} \\ \tau_{3,t} \\ \tau_{4,t} \\ \tau_{5,t} \\ \tau_{6,t} \\ \tau_{7,t} \\ \frac{\tau_{8,t}}{\eta_{8}} \\ \frac{\tau_{8,t}}{\eta_{9}} \end{array} \right) + \left(
\begin{array}{cccc}
1 \\
\Upsilon_{1,1} + \Upsilon_{1,2} L + \ldots + \Upsilon_{1,p} L^{p-1} \\
\Upsilon_{2,1} + \Upsilon_{2,2} L + \ldots + \Upsilon_{2,p} L^{p-1} \\
\Upsilon_{3,1} + \Upsilon_{3,2} L + \ldots + \Upsilon_{3,p} L^{p-1} \\
\Upsilon_{4,1} + \Upsilon_{4,2} L + \ldots + \Upsilon_{4,p} L^{p-1} \\
\Upsilon_{5,1} + \Upsilon_{5,2} L + \ldots + \Upsilon_{5,p} L^{p-1} \\
\Upsilon_{6,1} + \Upsilon_{6,2} L + \ldots + \Upsilon_{6,p} L^{p-1} \\
\Upsilon_{7,1} + \Upsilon_{7,2} L + \ldots + \Upsilon_{7,p} L^{p-1} \\
\Upsilon_{8,1} + \Upsilon_{8,2} L + \ldots + \Upsilon_{8,p} L^{p-1}
\end{array} \right)
\psi_{1,t} + \left( \begin{array}{c} \psi_{2,t} \\ \psi_{3,t} \\ \psi_{4,t} \\ \psi_{5,t} \\ \psi_{6,t} \\ \psi_{7,t} \\ \psi_{8,t} \\ \psi_{9,t} \\ \psi_{10,t} \end{array} \right) + \boldsymbol{\mathbf{\xi}}_{t}
\end{align*}
where $\psi_{1,t}$ is a causal AR($p$) denoting the business cycle; $\tau_{1,t}, \ldots, \tau_{8,t}$ are second-order smooth trends kitagawa1996smoothness; headline and core inflation share a common trend (the so-called trend inflation); $\psi_{2, t}, \ldots, \psi_{10, t}$ are causal AR(1) idiosyncratic noises; $\boldsymbol{\mathbf{\xi}}_t
\mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}}
N \left(\boldsymbol{\mathbf{0}}_{9 \times 1}, \varepsilon \cdot \boldsymbol{\mathbf{I}}_{9} \right)$ for a small positive $\varepsilon$, similarly to banbura2014maximum.\footnote{In this empirical example, $\varepsilon = 10^{-4}$.} Going forward, the number of lags $p$ is set to $12$ months in order to allow for a relatively large memory in the business cycle. The dynamics for the states and estimation method for this trend-cycle model are further detailed in (ref).
remark[Trend inflation]
Post estimation, the standardisation is removed to use the original units. In doing so, the scaling factors associated with trend inflation are also removed. Hence, headline and core inflation have the same trend once the standardisation is lifted.
The factor-augmented ensembles extend the information set of traditional autoregression trees including with the estimated business cycle. The latter is used both in levels and with a selected range of transformations. Formally, in order to compute a prediction referring to a generic time $t+1$, the factor-augmented ensembles use a vector of predictors containing the target values referring to time $t, \ldots, t-11$ and the augmentation
align[align omitted — 466 chars of source]
where $\hat{\psi}_{1,t+j|t}$ denotes the estimate of the business cycle for a generic period $t+j$ computed with the information set available at time $t$. While the first block in the factor augmentation gives a direct view on the business cycle levels, the following ones are useful for computing splits directly on its turning points and making a better use of the data.
Each ensemble is regulated via a vector of hyperparameters that includes those specifics to the state-space estimation and the minimum number of observations per leaf. These tuning parameters are determined on a sample going from January 1984 to the end of January 2005. The ALFRED data vintage used for structuring the macroeconomic selection sample includes the information was available right before the end of January 2005. Since this article uses a two-step method, hyperparameters are selected first for the trend-cycle model and then for the factor-augmented ensembles. The trend-cycle model is tuned as illustrated in (ref). Next, the minimum number of observations per leaf of each ensemble is determined with a pseudo out-of-sample criterion and a grid search on the equally spaced $\mathscr{H}_{RT} \vcentcolon= \{0.01, 10, 15, \ldots, 0.5\}$ with $|\mathscr{H}_{RT}|=25$.\footnote{This difference in the selection method is determined by the higher computational complexity required to estimate and forecast with factor-augmented tree ensembles.} The minimum number of observations per leaf is expressed in percentage terms with respect to the number of time periods available. Both steps use the first half of the selection sample for the estimation and the second half to validate the results.
In-sample results
This subsection presents key in-sample findings to aid understanding each stage involved in constructing these factor-augmented ensembles
figure[figure omitted — 203 chars of source]
figure[figure omitted — 217 chars of source]
In the first step, the model extracts persistent and transitory components from macroeconomic data. \Crefrange{ch2:fig:tc_model_trends}{ch2:fig:tc_model_cycles} shows an in-sample snapshot of the economy captured by this trend-cycle decomposition based on monthly data from January 1984 to February 2020. This is the full sample available on the 28th February 2020 as it was recorded on the ALFRED database. In other words, the last pre COVID-19 vintage for the macroeconomic data in this analysis. (ref) compares the data in levels with the estimated trends, while (ref) decomposes the cycles into common and idiosyncratic fluctuations. The results are compatible with papers including hasenzagl2022a, hasenzagl2022b and barigozzi2021measuring, and show that, while there is a strong heterogeneity across macroeconomic indicators, the business cycle explains most of the cyclical fluctuations. Besides, these findings also show that the business cycle is synchronised with the NBER dating, which means that expansions (contractions) in this latent component are usually associated with economic growth (recessions).
figure[figure omitted — 449 chars of source]
The second step builds on the trend-cycle decomposition and uses the business cycle -- both in levels and transformed as in (ref) -- to extend the information set of tree ensembles for equity volatility. One of the advantages of tree-based models with respect to alternative machine learning techniques is that they are highly interpretable and allow to inspect the predictability drivers. (ref) and the additional figures in (ref) compare the factor-augmented ensembles estimated pre and post COVID-19 (i.e., estimating it first with data as at 28th February 2020 and then with the latest vintage available) by looking at the bagging importance weights: the number of times, in percentage points, that a predictor is selected to create a split in the underlying factor-augmented regression trees which is, essentially, an internal ranking. (ref) reports the total importance of the lagged target and factor augmentation. Mid to small cap shares are usually more vulnerable than blue chips to changes in economic conditions. Indeed, a broad range of papers such as gertler1994monetary, sharpe1994financial and bernanke1996financial argue that small firms do not have a broad range of financing options and mostly use intermediaries to access credit. This leaves them more at risk during a downturn when banks become more selective with respect to credit extensions. Therefore, it is not surprising that the factor augmentation is especially crucial for the Wilshire indices referring to these market capitalisations. In addition, (ref) highlights how the factor augmentation is even more relevant post COVID-19, a period of unprecedented high volatility and uncertainty. \Crefrange{ch2:fig:importance_disaggr_pre}{ch2:fig:importance_disaggr_post} highlight a high heterogeneity across targets.
Real-time evaluation
Influential studies orphanides2002unreliability have warned about the potential risks of utilising cyclical latent variables of economic significance, such as the business cycle or the output gap, without ensuring that the associated forecasting benefits are robust. As a result, this subsection assesses the real-time predictability of these factor-augmented tree ensembles both pre and post COVID-19.
I have started this analysis supposing to be right before the end of January 2005, the moment in time when the hyperparameters were calibrated (cf. (ref)). Based on the information available at that time, I estimated the ensembles, in turn, for each target and produced the corresponding one-step-ahead forecasts. Subsequently, I performed out-of-sample testing by re-estimating the models and generating fresh predictions every time a new ALFRED vintage became available. The macroeconomic test sample comprises 861 vintages, each with a minimum of 2277 observations, and, in adherence to the aforementioned guidelines, the evaluation does not look forward in the future.
table[table omitted — 1,851 chars of source]
(ref) provides a summary of the out-of-sample findings, showing the mean squared error associated with the bagging forecasts relative to the one for a naive benchmark constant at zero. In other words an metrics that, when lower than one, indicates the ensembles are more accurate than this simple model. Besides, (ref) compares standard autoregressive ensembles (i.e., without the factor-augmentation) with those that leverage on the business cycle and its transformation. These figures reveal two crucial findings. First, the factor-augmented models are more accurate than classic bagging for all targets and especially for the volatility of mid to small cap stocks, in accordance with (ref). Second, almost all ensembles outperform the corresponding naive benchmark, which may appear obvious initially, but not necessarily so for linear forecasting models such as those detailed in (ref). Hence, this provides further confirmation of the superiority of bootstrap aggregating.
Concluding remarks
Most econometric techniques rely on simplifying assumptions that are typically based on highly stylised theoretical models. While these assumptions enable researchers to obtain easily interpretable estimates, it is challenging to justify their use for complex real-world problems, which are frequently characterised by non-linear relationships with unclear parametric forms -- this is especially true when modelling the joint dynamics between macroeconomic and financial data. Conversely, machine learning techniques are unconstrained by the nature of the problem, but heavily reliant on observable data and often difficult to interpret. This manuscript proposes to bridge the gap between these families of techniques by developing a hybrid between structural time series harvey1985trends, harvey1990forecasting and machine learning. More specifically, the manuscript presents an approach where indicators with unclear or complex dynamics are connected to latent states of economic interest using tree ensembles breiman1996bagging, breiman2001random. This is a two-stage method since these latent variables are estimated separately and before the ensembles via state-space models informed by domain-specific theory.
From a methodological perspective, this is advantageous for two reasons. First, it generalises tree-based regressions to allow for predictors with non-stationary trends, seasonality, peculiar cyclical patterns and missing values. Second, this approach permits to inform ensembles with domain-specific theory extending their information set using latent variables with structural interpretation. Empirically, this is of particular relevance for macro-finance. Indeed, it is relatively straightforward to link returns and other financial variables of interest with unobserved components representing economic concepts such as the output gap. This article describes in greater detail the lead-lag effect between equity volatility and the US business cycle and finds that the current economic conditions are strong predictors of equity volatility when allowing for non-linearities. This is studied both in-sample and out-of-sample, and results are particularly strong for stocks with small market capitalisation. The latter is in agreement with the broad body of literature including gertler1994monetary, sharpe1994financial and bernanke1996financial that find similar firms more vulnerable to downturns due to a lack of financing options.
Factor-augmented ensembles could also be beneficial for related empirical problems such as the study of links between foreign exchange rates and fundamentals meese1983empirical, engel2005exchange, rossi2013exchange. Similar problems are left open for future research.
{
\titlespacing
0pt
{12pt plus 4pt minus 2pt}{0pt plus 2pt minus 2pt}
{1.5ex}
\theendnotes
}
{
\singlespacing
}
appendix\section{Business cycle estimation}
\subsection{Trend-cycle model}
Re-write the model in (ref) in the state-space form
\begin{align*}
\boldsymbol{\mathbf{Z}}_{t} &= \boldsymbol{\mathbf{B}} \boldsymbol{\mathbf{\Phi}}_{t} + \boldsymbol{\mathbf{\xi}}_{t}, \\
\boldsymbol{\mathbf{\Phi}}_{t} &= \boldsymbol{\mathbf{C}} \boldsymbol{\mathbf{\Phi}}_{t-1} + \boldsymbol{\mathbf{D}} \boldsymbol{\mathbf{\zeta}}_{t}.
\end{align*}
The innovation in the transition equation $\boldsymbol{\mathbf{\zeta}}_{t}
\mathrel{\overset{w.n}{\scalebox{1.5}[1]{$\sim$}}}
N \left(\boldsymbol{\mathbf{0}}_{r \times 1}, \boldsymbol{\mathbf{\Sigma}} \right)$ with $\boldsymbol{\mathbf{\Sigma}}$ being a $r \times r$ positive definite real diagonal matrix and $r=18$. The transition matrices are sparse and the non-zero entries are such that
\begin{align*}
&\boldsymbol{\mathbf{C}} \vcentcolon= \left(
\begin{array}{ccc | ccc | ccc | ccccc}
1 & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \ddots & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & 1 & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\hline
\cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\hline
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \pi_{1} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \pi_{n} & \cdot & \cdot & \cdot & \cdot & \cdot \\
\hline
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \pi_{n+1} & \pi_{n+2} & \ldots & \pi_{n+p-1} & \pi_{n+p} \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \ldots & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \ddots & \vdots & \vdots \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \vdots & \ddots & \ddots & \vdots & \vdots \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ldots & \ldots & 1 & \cdot \\
\end{array} \right),\\[-1em]
& \begin{array}{ccccccc ccc cccc ccccc}
\multicolumn{7}{c}{\underbrace_{q \times 8}} &
\multicolumn{3}{c}{\underbrace_{q \times 8}} &
\multicolumn{4}{c}{\underbrace_{q \times 9}} & \multicolumn{5}{c}{\underbrace_{q \times p}}
\end{array} \\[0.4em]
&\boldsymbol{\mathbf{D}} \vcentcolon= \left(
\begin{array}{ccc | ccc | ccc}
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\hline
1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot \\
\hline
\cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot \\
\hline
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot
\end{array} \right),\\[-1em]
& \begin{array}{ccc ccc ccc}
\multicolumn{3}{c}{\underbrace_{q \times 8}} &
\multicolumn{3}{c}{\underbrace_{q \times 9}} & \multicolumn{3}{c}{\underbrace_{q \times 1}}
\end{array}
\end{align*}
where $\boldsymbol{\mathbf{\pi}}$ is a $n+p \times 1$ vector of finite real parameters and $q=25+p$. The measurement coefficient matrix is also sparse and the non-zero entries are such that
\begin{align*}
&\boldsymbol{\mathbf{B}} \vcentcolon= \left(
\begin{array}{cccccccc | ccc | ccccccccc | cccc }
1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ldots & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \ldots & \cdot \\
\cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ldots & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \tilde{\Upsilon}_{1,1} & \tilde{\Upsilon}_{1,2} & \ldots & \tilde{\Upsilon}_{1,p} \\
\cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ldots & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \tilde{\Upsilon}_{2,1} & \tilde{\Upsilon}_{2,2} & \ldots & \tilde{\Upsilon}_{2,p} \\
\cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \tilde{\Upsilon}_{3,1} & \tilde{\Upsilon}_{3,2} & \ldots & \tilde{\Upsilon}_{3,p} \\
\cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \tilde{\Upsilon}_{4,1} & \tilde{\Upsilon}_{4,2} & \ldots & \tilde{\Upsilon}_{4,p} \\
\cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \tilde{\Upsilon}_{5,1} & \tilde{\Upsilon}_{5,2} & \ldots & \tilde{\Upsilon}_{5,p} \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \tilde{\Upsilon}_{6,1} & \tilde{\Upsilon}_{6,2} & \ldots & \tilde{\Upsilon}_{6,p} \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \frac{1}{\eta_8} & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \tilde{\Upsilon}_{7,1} & \tilde{\Upsilon}_{7,2} & \ldots & \tilde{\Upsilon}_{7,p} \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \frac{1}{\eta_9} & \cdot & \ldots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \tilde{\Upsilon}_{8,1} & \tilde{\Upsilon}_{8,2} & \ldots & \tilde{\Upsilon}_{8,p}
\end{array} \right),\\[-1em]
& \begin{array}{cccccccc ccc ccccccccc cccc}
\multicolumn{8}{c}{\underbrace_{9 \times 8}} &
\multicolumn{3}{c}{\underbrace_{9 \times 8}} &
\multicolumn{9}{c}{\underbrace_{9 \times 9}} & \multicolumn{4}{c}{\underbrace_{9 \times p}}
\end{array}
\end{align*}
where $\boldsymbol{\mathbf{\tilde{\Upsilon}}}$ is a $8 \times p$ matrix of finite real parameters and $\tilde{\Upsilon}_{i,j} \approx \Upsilon_{i,j}$, for any $i=1,\ldots,8$ and $j=1, \ldots, p$.\footnote{This numerical approximation is governed by $\varepsilon$ in (ref).} This representation implies that
\begin{align*}
&\boldsymbol{\mathbf{\Phi}}_{t} \vcentcolon= \left( \begin{array}{ccc | ccc | cccc | cccc} \tau_{1,t} & \ldots & \tau_{8,t} & \delta_{1,t} & \ldots & \delta_{8,t} & \psi_{2,t} & \psi_{3,t} & \ldots & \psi_{10, t} & \psi_{1,t} & \psi_{1,t-1} & \ldots & \psi_{1,t-p+1} \end{array} \right)'. \\[-1em]
& \begin{array}{ccc ccc cccc cccc}
\multicolumn{3}{c}{\underbrace_{8 \times 1}} & \multicolumn{3}{c}{\underbrace_{8 \times 1}} &
\multicolumn{4}{c}{\underbrace_{9 \times 1}} &
\multicolumn{4}{c}{\underbrace_{p \times 1}}
\end{array}
\end{align*}
The initial conditions for the states are such that $\boldsymbol{\mathbf{\Phi}}_{0}
\mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}}
N (\boldsymbol{\mathbf{\mu}}_{0}, \boldsymbol{\mathbf{\Omega}}_{0})$, where $\boldsymbol{\mathbf{\mu}}_{0}$ and $\boldsymbol{\mathbf{\Omega}}_{0}$ denote a $q \times 1$ real vector and a $q \times q$ positive definite real covariance matrix. The matrix $\boldsymbol{\mathbf{\Omega}}_{0}$ is sparse and the entries that can differ from zero are those with coordinates $(i,j) \in \{(i,j) : i=j \text{ and } 1 \leq i \leq 25\} \cup \{(i,j) : 25 < i \leq q \text{ and } 25 < j \leq q\}$.
\begin{remark}
The empirical application assumes that $\boldsymbol{\mathbf{\Sigma}}$ is diagonal. This implies an exact factor model (i.e., no cross-sectional dependence in the idiosyncratic components) and, in doing so, it simplifies the narrative. However, this assumption may be too restrictive for more general problems. For similar problems, the assumptions could be relaxed and the ECM algorithm described in this appendix could be presented as a penalised quasi maximum likelihood estimation method building on the theoretical results in barigozzi2020quasi.
\end{remark}
\subsection{The Expectation-Conditional Maximisation algorithm}
Denote the model free parameters with
\begin{align*}
&\boldsymbol{\mathbf{\vartheta}} \vcentcolon= \left( \begin{array}{cccccccc} \boldsymbol{\mathbf{\mu}}_0' & vech(\boldsymbol{\mathbf{\Omega}}_0)' & vec(\boldsymbol{\mathbf{\tilde{\Upsilon}}})' & \boldsymbol{\mathbf{\pi}}' & \Sigma_{1,1} & \Sigma_{2,2} & \ldots & \Sigma_{r, r} \end{array} \right)'.
\end{align*}
The ECM algorithm estimates these coefficients by repeating the optimisation process illustrated in (ref) until convergence.
\begin{definition}[ECM estimation routine]
At any $k+1 > 1$ iteration, the ECM algorithm computes the vector of coefficients
\begin{align*}
\boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) \vcentcolon= \operatorname*{arg\,max}_{\boldsymbol{\mathbf{\vartheta}} \, \in \, \mathscr{R}} \, \mathbb{E} \left [\mathcal{L}(\boldsymbol{\mathbf{\vartheta}} \,|\, \boldsymbol{\mathbf{Z}}_{1:s}, \boldsymbol{\mathbf{\Phi}}_{1:s}) \,|\, \mathscr{Z}(s), \, \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] - \mathbb{E} \left [\mathcal{P}(\boldsymbol{\mathbf{\vartheta}}, \boldsymbol{\mathbf{\gamma}}) \,|\, \mathscr{Z}(s), \, \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right],
\end{align*}
where $\mathscr{R}$ denotes the region in which the AR cycles (common and idiosyncratic) are causal, $\mathscr{Z}(s)$ is the information set available at time $s$,
\begin{align}
\mathcal{L}(\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \boldsymbol{\mathbf{Z}}_{1:s}, \boldsymbol{\mathbf{\Phi}}_{1:s}) \simeq &-\frac{1}{2}\ln|\underline{\boldsymbol{\mathbf{\Omega}}_0}| - \frac{1}{2}\operatorname{Tr}\Big[\underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} (\boldsymbol{\mathbf{\Phi}}_0 - \underline{\boldsymbol{\mathbf{\mu}}_0})(\boldsymbol{\mathbf{\Phi}}_0 - \underline{\boldsymbol{\mathbf{\mu}}_0})'\Big] \\
&-\frac{s}{2}\ln|\underline{\boldsymbol{\mathbf{\Sigma}}}| - \frac{1}{2} \operatorname{Tr} \Big[\sum_{t=1}^s \underline{\boldsymbol{\mathbf{\Sigma}}}^{-1} (\boldsymbol{\mathbf{\Phi}}_{*, t} - \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\Phi}}_{t-1})(\boldsymbol{\mathbf{\Phi}}_{*, t} - \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\Phi}}_{t-1})'\Big] \nonumber \\
&-\frac{s}{2}\ln|\underline{\boldsymbol{\mathbf{R}}}| - \frac{1}{2} \operatorname{Tr}\Big[\sum_{t=1}^s \underline{\boldsymbol{\mathbf{R}}}^{-1} (\boldsymbol{\mathbf{Z}}_t - \underline{\boldsymbol{\mathbf{B}}}\boldsymbol{\mathbf{\Phi}}_{t})(\boldsymbol{\mathbf{Z}}_t - \underline{\boldsymbol{\mathbf{B}}}\boldsymbol{\mathbf{\Phi}}_{t})'\Big], \nonumber
\end{align}
$\boldsymbol{\mathbf{\Phi}}_{*, t} \equiv \underline{\boldsymbol{\mathbf{D}}}' \boldsymbol{\mathbf{\Phi}}_{t}$, $\underline{\boldsymbol{\mathbf{C}}}_{\,*} \equiv \underline{\boldsymbol{\mathbf{D}}}' \underline{\boldsymbol{\mathbf{C}}}$, and the underlined matrices denote the state-space coefficients implied by $\underline{\boldsymbol{\mathbf{\vartheta}}}$. Besides,
\begin{align*}
\mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}) \vcentcolon= +\frac{1-\alpha}{2} &\left( \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,1:n} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1)^{\frac{1}{2}} \big\Vert_{\text{F}}^2 + \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,n+1:n+p}' \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1)^{\frac{1}{2}} \big\Vert_{\text{F}}^2 + \big\Vert \underline{\boldsymbol{\mathbf{\tilde{\Upsilon}}}} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p)^{\frac{1}{2}} \big\Vert_{\text{F}}^2 \right) \\
+ \frac{\alpha}{2} &\left( \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,1:n} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1) \big\Vert_{1,1} + \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,n+1:n+p}' \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1) \big\Vert_{1,1} + \big\Vert \underline{\boldsymbol{\mathbf{\tilde{\Upsilon}}}} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p) \big\Vert_{1,1} \right)
\end{align*}
where, for any $l \in \mathbb{N}$,
\begin{align*}
\boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, l) \vcentcolon= \lambda \begin{pmatrix}
1 & 0 &\ldots & 0 \\
0 &\beta &\ldots & 0 \\
\vdots &\ddots &\ddots &\vdots \\
0 &\ldots &\ldots & \beta^{l-1} \end{pmatrix},
\end{align*}
$\lambda \geq 0$, $0 \leq \alpha \leq 1$ and $\beta \geq 1$ are hyperparameters included in $\boldsymbol{\mathbf{\gamma}}$. The state-space coefficients for the first iteration are initialised as in (ref).
\end{definition}
\begin{remark}[Objective functions]
The function in (ref) is the so-called complete-data (i.e., fully observed data and known latent states) log-likelihood, while $\mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}})$ represents the generalised elastic-net penalty used in pellegrino2020selecting.
\end{remark}
\begin{remark}[Underlined coefficients]
Some of the underlined state-space coefficients are partially or fully fixed in accordance with the structure in (ref). For instance, $\boldsymbol{\mathbf{D}} = \underline{\boldsymbol{\mathbf{D}}}$ since $\boldsymbol{\mathbf{D}}$ does not contain free parameters.
\end{remark}
\begin{assumption}[Convergence]
The ECM algorithm is said to be converged when the criteria in (ref) are met.
\end{assumption}
The optimisation in (ref) is performed in two steps. The first one (E-step) involves the computation of the expectations in (ref). The second step (CM-step) conditionally maximises the resulting expected penalised log-likelihood with respect to the free parameters.
It is convenient to write down the E-step on the basis of the output of a Kalman smoother compatible with incomplete time series, as in shumway1982approach and watson1983alternative. The required output is introduced in (ref) and used in (ref) to compute the expected log-likelihood.
\begin{definition}[Kalman smoother output]
The hereinbefore mentioned Kalman smoother output is
\begin{align*}
&\boldsymbol{\mathbf{\hat{\Phi}}}_{t} \vcentcolon= \mathbb{E} \, \Big[\boldsymbol{\mathbf{\Phi}}_t \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big], \\
&\boldsymbol{\mathbf{\hat{P}}}_{t,t-j} \vcentcolon= \text{Cov} \, \Big[\boldsymbol{\mathbf{\Phi}}_t, \boldsymbol{\mathbf{\Phi}}_{t-j} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big],
\end{align*}
for any $k \geq 0$, $0 \leq j \leq t$ and $t \geq 0$. Let also $\boldsymbol{\mathbf{\hat{P}}}_{t} \equiv \boldsymbol{\mathbf{\hat{P}}}_{t,t}$.
\end{definition}
\begin{remark}
These estimates are computed as in pellegrino2020selecting.
\end{remark}
Furthermore, (ref) is useful to formalise which measurements are observed at every single point in time.
\begin{definition}[Observed measurements]
Let
\begin{align*}
&\mathscr{T} \vcentcolon= \bigcup_{i=1}^n \mathscr{T}_i, \\
&\mathscr{T}(s) \vcentcolon= \{t : t \in \mathscr{T}, \, 1 \leq t \leq s\},
\end{align*}
describe two sets representing the points in time (either over the full sample or up to time $s$) in which we observe at least one measurement, for $1 \leq s \leq T$. Let also
\begin{align*}
\mathscr{V}_{t} \vcentcolon= \{i : t \in \mathscr{T}_i, \, 1 \leq i \leq n\},
\end{align*}
for $1 \leq t \leq T$. Thus, let
\begin{align*}
&\boldsymbol{\mathbf{Z}}_{t}^{obs} \vcentcolon= \big( Z_{i,t} \big)_{i \in \mathscr{V}_t} \\
&\boldsymbol{\mathbf{B}}_{t}^{obs} \vcentcolon= \boldsymbol{\mathbf{A}}_{t} \boldsymbol{\mathbf{B}}
\end{align*}
be the vector of observed measurements at time $t$ and the corresponding $|\mathscr{V}_t| \times q$ matrix of coefficients, for any $t \in \mathscr{T}$. Every $\boldsymbol{\mathbf{A}}_t$ is indeed a selection matrix constituted by ones and zeros that permits to retrieve the appropriate rows of $\boldsymbol{\mathbf{B}}$ for every $t \in \mathscr{T}$.
\end{definition}
\begin{proposition}
Let
\begin{align*}
\mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \equiv \mathbb{E} \left [\mathcal{L}(\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \boldsymbol{\mathbf{Z}}_{1:s}, \boldsymbol{\mathbf{\Phi}}_{1:s}) \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right].
\end{align*}
Building on (ref), it follows that
\begin{align*}
\mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \simeq &-\frac{1}{2}\ln|\underline{\boldsymbol{\mathbf{\Omega}}_0}| - \frac{1}{2}\operatorname{Tr}\Big[\underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} (\boldsymbol{\mathbf{\hat{E}}} - \boldsymbol{\mathbf{\hat{\Phi}}}_0 \underline{\boldsymbol{\mathbf{\mu}}_0}' - \underline{\boldsymbol{\mathbf{\mu}}_0} \boldsymbol{\mathbf{\hat{\Phi}}}_0' + \underline{\boldsymbol{\mathbf{\mu}}_0} \, \underline{\boldsymbol{\mathbf{\mu}}_0}') \Big] \\
&-\frac{s}{2}\ln|\underline{\boldsymbol{\mathbf{\Sigma}}}| - \frac{1}{2} \operatorname{Tr} \Big[\underline{\boldsymbol{\mathbf{\Sigma}}}^{-1} (\boldsymbol{\mathbf{\hat{F}}}_s - \boldsymbol{\mathbf{\hat{G}}}_s \underline{\boldsymbol{\mathbf{C}}}_{\,*}' - \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\hat{G}}}_s' + \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\hat{H}}}_s \underline{\boldsymbol{\mathbf{C}}}_{\,*}') \Big] \\
&-\frac{1}{2 \varepsilon} \operatorname{Tr} \left\{ \sum_{t \in \mathscr{T}(s)} \left[ \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right) \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right)' + \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{P}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{{obs}'} \right] \right\},
\end{align*}
where
\begin{align*}
&\boldsymbol{\mathbf{\hat{E}}} \vcentcolon= \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_0 \boldsymbol{\mathbf{\Phi}}_0' \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] = \boldsymbol{\mathbf{\hat{\Phi}}}_0 \boldsymbol{\mathbf{\hat{\Phi}}}_0' + \boldsymbol{\mathbf{\hat{P}}}_0, \\
&\boldsymbol{\mathbf{\hat{F}}}_s \vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{*, t} \boldsymbol{\mathbf{\Phi}}_{*, t}' \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] = \sum_{t=1}^s \underline{\boldsymbol{\mathbf{D}}}' \left(\boldsymbol{\mathbf{\hat{\Phi}}}_t \boldsymbol{\mathbf{\hat{\Phi}}}_t' + \boldsymbol{\mathbf{\hat{P}}}_t\right) \underline{\boldsymbol{\mathbf{D}}}, \\
&\boldsymbol{\mathbf{\hat{G}}}_s \vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{*, t} \boldsymbol{\mathbf{\Phi}}_{t-1}' \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] = \sum_{t=1}^s \underline{\boldsymbol{\mathbf{D}}}' \left(\boldsymbol{\mathbf{\hat{\Phi}}}_t \boldsymbol{\mathbf{\hat{\Phi}}}_{t-1}' + \boldsymbol{\mathbf{\hat{P}}}_{t, t-1}\right), \\
&\boldsymbol{\mathbf{\hat{H}}}_s \vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{t-1} \boldsymbol{\mathbf{\Phi}}_{t-1}' \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] = \sum_{t=1}^s \left(\boldsymbol{\mathbf{\hat{\Phi}}}_{t-1} \boldsymbol{\mathbf{\hat{\Phi}}}_{t-1}' + \boldsymbol{\mathbf{\hat{P}}}_{t-1}\right).
\end{align*}
\end{proposition}
\begin{proof}
The proof is analogous to the one in pellegrino2020selecting.
\end{proof}
\begin{remark}
$\boldsymbol{\mathbf{D}} = \underline{\boldsymbol{\mathbf{D}}}$ plays the role of a selection matrix. In particular premultiplying by $\underline{\boldsymbol{\mathbf{D}}}'$ allows to select rows and postmultiplying by $\underline{\boldsymbol{\mathbf{D}}}$ columns.
\end{remark}
\begin{lemma}
The conditional expectation for the penalty in (ref) is
\begin{align*}
\mathbb{E} \left[ \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}) \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] = \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}).
\end{align*}
\end{lemma}
\begin{proof}
A formal proof is not reported since it is immediate. Indeed, the penalty function in this ECM algorithm depends only on the current vector of coefficients and hyperparameters.
\end{proof}
The CM-step conditionally maximises the expected penalised log-likelihood
\begin{align}
\mathcal{M}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \vcentcolon= \mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] - \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}})
\end{align}
to estimate the state-space parameters. The estimated coefficients are denoted with an “hat” symbol. Besides, an $s$ subscript is used for highlighting the sample size and a superscript for denoting the ECM iteration.
\begin{lemma}
The ECM estimator at a generic iteration $k+1 > 0$ for $\boldsymbol{\mathbf{\mu}}_0$ is
\begin{align*}
&\boldsymbol{\mathbf{\hat{\mu}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \boldsymbol{\mathbf{\hat{\Phi}}}_0
\end{align*}
and the estimator for $\boldsymbol{\mathbf{\Omega}}_0$ is a sparse covariance matrix whose non-zero entries are
\begin{align*}
&\Big[ \boldsymbol{\mathbf{\hat{\Omega}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) \Big]_{i,j} = \Big[\boldsymbol{\mathbf{\hat{P}}}_0\Big]_{i,j},
\end{align*}
for $(i,j) \in \{(i,j) : i=j \text{ and } 1 \leq i \leq 25\} \cup \{(i,j) : 25 < i \leq q \text{ and } 25 < j \leq q\}$.
\end{lemma}
\begin{proof}
The derivative of (ref) with respect to $\underline{\boldsymbol{\mathbf{\mu}}_0}$ is
\begin{align*}
&\frac{\partial \mathcal{M}_e \left[ \underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}} \,|\, \mathscr{Z}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right]}{\partial \underline{\boldsymbol{\mathbf{\mu}}_0}} = -\frac{1}{2} \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} \left(-2\boldsymbol{\mathbf{\hat{\Phi}}}_0 +2\underline{\boldsymbol{\mathbf{\mu}}_0} \right).
\end{align*}
It follows that the maximiser for the expected penalised log-likelihood is
\begin{align*}
\boldsymbol{\mathbf{\hat{\mu}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \boldsymbol{\mathbf{\hat{\Phi}}}_0.
\end{align*}
The derivative of (ref) with respect to $\underline{\boldsymbol{\mathbf{\Omega}}_0}$ and fixing $\underline{\boldsymbol{\mathbf{\mu}}_0} = \boldsymbol{\mathbf{\hat{\mu}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}})$ is
\begin{align*}
&-\frac{1}{2} \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} +\frac{1}{2} \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} \left[\boldsymbol{\mathbf{\hat{E}}} - \boldsymbol{\mathbf{\hat{\Phi}}}_0 \boldsymbol{\mathbf{\hat{\Phi}}}_0' \right] \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} = -\frac{1}{2} \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} +\frac{1}{2} \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1} \boldsymbol{\mathbf{\hat{P}}}_0 \, \underline{\boldsymbol{\mathbf{\Omega}}_0}^{-1}.
\end{align*}
as also shown in pellegrino2020selecting. Given the assumptions in (ref) on the structure of $\boldsymbol{\mathbf{\Omega}}_0$, it follows that $\boldsymbol{\mathbf{\hat{\Omega}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}})$ is a sparse matrix whose non-zero entries are
\begin{align*}
&\Big[ \boldsymbol{\mathbf{\hat{\Omega}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) \Big]_{i,j} = \Big[\boldsymbol{\mathbf{\hat{P}}}_0\Big]_{i,j},
\end{align*}
for $(i,j) \in \{(i,j) : i=j \text{ and } 1 \leq i \leq 25\} \cup \{(i,j) : 25 < i \leq q \text{ and } 25 < j \leq q\}$.
\end{proof}
\begin{definition}
Let $\boldsymbol{\mathbf{\tilde{\Gamma}}}(\boldsymbol{\mathbf{\gamma}})$ be a diagonal $q \times q$ matrix whose non-zero entries are such that
\begin{align*}
\boldsymbol{\mathbf{\tilde{\Gamma}}}(\boldsymbol{\mathbf{\gamma}}) \vcentcolon= \begin{pmatrix}
\cdot & \cdot & \cdot \\
\cdot & \lambda \, \boldsymbol{\mathbf{I}}_{9} & \cdot \\
\cdot & \cdot & \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p)
\end{pmatrix}.
\end{align*}
\end{definition}
\begin{lemma}
The ECM estimator at a generic iteration $k+1 > 0$ for $\boldsymbol{\mathbf{C}}$ is such that
\begin{align*}
\hat{C}_{i,j,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \frac{\mathcal{S} \, \left[ \hat{\Sigma}_{i, i, s}^{k^{-1}}(\boldsymbol{\mathbf{\gamma}}) \big( \hat{G}_{i, j, s} - \sum_{\substack{l=1,\, l \neq j}}^{q} \hat{C}_{i, l, s}^{k+\mathbb{I}_{l<j}}(\boldsymbol{\mathbf{\gamma}}) \, \hat{H}_{l, j, s} \big), \; \frac{\alpha}{2} \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \right]}{\hat{\Sigma}_{i,i,s}^{k^{-1}} (\boldsymbol{\mathbf{\gamma}}) \, \hat{H}_{j,j,s} + (1-\alpha) \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}})},
\end{align*}
for any $(i,j) \in \{(i,j) : i=j \text{ and } 17 \leq i \leq 25\} \cup \{(i,j) : i=26 \text{ and } 26 \leq j \leq q\}$, and constant to the values in (ref) for the remaining entries.
\end{lemma}
\begin{proof}
Given that the absolute value function in the penalty is not differentiable at zero, this part of the ECM algorithm estimates, in turn, the free entries of $\boldsymbol{\mathbf{C}}$ (i.e., $\pi_{1}, \ldots, \pi_{n+p}$) while fixing $\underline{\boldsymbol{\mathbf{\Sigma}}} = \boldsymbol{\mathbf{\hat{\Sigma}}}^k_s(\boldsymbol{\mathbf{\gamma}})$ and any other free entry of $\boldsymbol{\mathbf{C}}$ to their latest estimate. For any $\underline{C}_{\,i,j} \neq 0$ corresponding to a free parameter, the derivative of (ref) with respect to $\underline{C}_{\,i,j}$ having fixed the coefficients as described in the previous sentence is
\begin{align*}
&+ \hat{\Sigma}_{i, i, s}^{k^{-1}}(\boldsymbol{\mathbf{\gamma}}) \left( \hat{G}_{i, j, s} - \underline{C}_{\,i,j} \hat{H}_{j,j,s} - \sum_{\substack{l=1 \\ l \neq j}}^q \hat{C}_{i, l, s}^{k+\mathbb{I}_{l < j}}\,(\boldsymbol{\mathbf{\gamma}}) \, \hat{H}_{l, j, s} \right) -(1-\alpha) \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \, \underline{C}_{\,i,j} - \frac{\alpha}{2} \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \, \text{sign} (\underline{C}_{\,i,j}),
\end{align*}
since $\boldsymbol{\mathbf{\hat{\Sigma}}}^k_s(\boldsymbol{\mathbf{\gamma}})$ is diagonal. It follows that
\begin{align*}
\hat{C}_{i,j,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \frac{\mathcal{S} \, \left[ \hat{\Sigma}_{i, i, s}^{k^{-1}}(\boldsymbol{\mathbf{\gamma}}) \big( \hat{G}_{i, j, s} - \sum_{\substack{l=1,\, l \neq j}}^{q} \hat{C}_{i, l, s}^{k+\mathbb{I}_{l<j}}(\boldsymbol{\mathbf{\gamma}}) \, \hat{H}_{l, j, s} \big), \; \frac{\alpha}{2} \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \right]}{\hat{\Sigma}_{i,i,s}^{k^{-1}} (\boldsymbol{\mathbf{\gamma}}) \, \hat{H}_{j,j,s} + (1-\alpha) \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}})},
\end{align*}
for any $(i,j) \in \{(i,j) : i=j \text{ and } 17 \leq i \leq 25\} \cup \{(i,j) : i=26 \text{ and } 26 \leq j \leq q\}$, and constant to the values in (ref) for the remaining entries.
\end{proof}
\begin{lemma}
The ECM estimator at a generic iteration $k+1 > 0$ for $\boldsymbol{\mathbf{\Sigma}}$ is such that
\begin{align*}
&\hat{\Sigma}_{i,i,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \frac{1}{s} \left[\boldsymbol{\mathbf{\hat{F}}}_s - \boldsymbol{\mathbf{\hat{G}}}_s \boldsymbol{\mathbf{\hat{C}}}^{{k+1}'}_{s}(\boldsymbol{\mathbf{\gamma}}) - \boldsymbol{\mathbf{\hat{C}}}^{k+1}_{s}(\boldsymbol{\mathbf{\gamma}})\, \boldsymbol{\mathbf{\hat{G}}}_s' + \boldsymbol{\mathbf{\hat{C}}}^{k+1}_{s}(\boldsymbol{\mathbf{\gamma}})\, \boldsymbol{\mathbf{\hat{H}}}_s \boldsymbol{\mathbf{\hat{C}}}^{{k+1}'}_{s}(\boldsymbol{\mathbf{\gamma}})\right]_{i,i}
\end{align*}
for $i=1, \ldots, r$ and zero for the remaining entries.
\end{lemma}
\begin{proof}
The proof is equivalent to the one reported in pellegrino2020selecting. However, in this manuscript, $\boldsymbol{\mathbf{\hat{\Sigma}}}^{k+1}_{s} (\boldsymbol{\mathbf{\gamma}})$ is diagonal as indicated in (ref).
\end{proof}
\begin{lemma}
Let
\begin{alignat*}{2}
&\boldsymbol{\mathbf{\hat{M}}}_{s} &&\vcentcolon= \sum_{t \in \mathscr{T}(s)} \boldsymbol{\mathbf{A}}_{t}' \boldsymbol{\mathbf{Z}}_{t}^{obs} \,\boldsymbol{\mathbf{\hat{\Phi}}}_{t}', \\[0.5em]
&\boldsymbol{\mathbf{\hat{N}}}_{t} &&\vcentcolon= \boldsymbol{\mathbf{A}}_{t}' \boldsymbol{\mathbf{A}}_{t}, \\[0.5em]
&\boldsymbol{\mathbf{\hat{O}}}_{t} &&\vcentcolon= \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' + \boldsymbol{\mathbf{\hat{P}}}_{t}.
\end{alignat*}
The ECM estimator at a generic iteration $k+1 > 0$ for $\boldsymbol{\mathbf{B}}$ is such that
\begin{align*}
\hat{B}_{i,j,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \frac{\mathcal{S} \, \left[ \hat{M}_{i,j,s} - \sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \sum_{\substack{l=1, l \neq j}}^{q} \hat{B}^{k+\mathbb{I}_{l < j}}_{i,l,s} (\boldsymbol{\mathbf{\gamma}}) \, \hat{O}_{l,j,t}, \; \frac{\alpha}{2} \, \varepsilon \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \right]}{\sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \hat{O}_{j,j,t} + (1-\alpha) \, \varepsilon \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}})},
\end{align*}
for any $(i,j) \in \{(i,j) : 2 \leq i \leq n \text{ and } 26 \leq j \leq q\}$, and constant to the values in (ref) for the remaining entries.
\end{lemma}
\begin{proof}
Note that
\begin{align*}
&\sum_{t \in \mathscr{T}(s)} \left[ \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right) \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right)' + \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{P}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{{obs}'} \right] \\
&\quad = \sum_{t \in \mathscr{T}(s)} \left[ \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right) \left( \boldsymbol{\mathbf{Z}}_{t}^{obs} - \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right)' + \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \boldsymbol{\mathbf{\hat{P}}}_{t} \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' \right] \\
&\quad = \sum_{t \in \mathscr{T}(s)} \left[ \boldsymbol{\mathbf{Z}}_{t}^{obs} \, \boldsymbol{\mathbf{Z}}_{t}^{{obs}'} - \boldsymbol{\mathbf{Z}}_{t}^{obs} \,\boldsymbol{\mathbf{\hat{\Phi}}}_{t}' \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' - \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \boldsymbol{\mathbf{Z}}_{t}^{{obs}'} + \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \left( \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' + \boldsymbol{\mathbf{\hat{P}}}_{t} \right) \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' \right].
\end{align*}
Note also that all $\boldsymbol{\mathbf{\hat{N}}}_{t}$ are diagonal. Indeed, at any point in time $t$ when all series are observed $\boldsymbol{\mathbf{A}}_{t} = \boldsymbol{\mathbf{\hat{N}}}_{t} = \boldsymbol{\mathbf{I}}_{n}$. Besides, at any other $t \in \mathscr{T}(s)$,
\begin{align*}
\hat{N}_{i,i,t} = \begin{cases}
1 & \text{if the } i \text{-th series is observed at time } t, \\
0 & \text{otherwise},
\end{cases}
\end{align*}
for $i=1, \ldots, n$. Given that the absolute value function in the penalty is not differentiable at zero, this part of the ECM algorithm estimates, in turn, the free entries of $\boldsymbol{\mathbf{B}}$ (i.e., $\tilde{\Upsilon}_{1, 1}, \ldots, \tilde{\Upsilon}_{1, p}, \ldots, \tilde{\Upsilon}_{8, p}$) while fixing any other free entry of $\boldsymbol{\mathbf{B}}$ to their latest estimate. For any $\underline{B}_{\,i,j} \neq 0$ corresponding to a free parameter, the derivative of (ref) with respect to $\underline{B}_{\,i,j}$ having fixed the coefficients as described in the previous sentence is
\begin{align*}
&+ \varepsilon^{-1} \left( \hat{M}_{i,j,s} - \sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \sum_{\substack{l=1 \\ l \neq j}}^{q} \hat{B}^{k+\mathbb{I}_{l < j}}_{i,l,s} (\boldsymbol{\mathbf{\gamma}}) \, \hat{O}_{l,j,t} \right) - \underline{B}_{\,i,j} \left( \varepsilon^{-1} \sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \hat{O}_{j,j,t} + (1-\alpha) \, \tilde{\Gamma}_{j,j} (\boldsymbol{\mathbf{\gamma}}) \right) \\[0.5em]
&\qquad -\frac{\alpha}{2} \, \tilde{\Gamma}_{j,j} (\boldsymbol{\mathbf{\gamma}}) \, \text{sign} (\underline{B}_{\,i,j}),
\end{align*}
since all $\boldsymbol{\mathbf{\hat{N}}}_{t}$ are diagonal. It follows that
\begin{align*}
\hat{B}_{i,j,s}^{k+1}(\boldsymbol{\mathbf{\gamma}}) = \frac{\mathcal{S} \, \left[ \hat{M}_{i,j,s} - \sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \sum_{\substack{l=1, l \neq j}}^{q} \hat{B}^{k+\mathbb{I}_{l < j}}_{i,l,s} (\boldsymbol{\mathbf{\gamma}}) \, \hat{O}_{l,j,t}, \; \frac{\alpha}{2} \, \varepsilon \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \right]}{\sum_{t \in \mathscr{T}(s)} \hat{N}_{i,i,t} \hat{O}_{j,j,t} + (1-\alpha) \, \varepsilon \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}})},
\end{align*}
for any $(i,j) \in \{(i,j) : 2 \leq i \leq n \text{ and } 26 \leq j \leq q\}$, and constant to the values in (ref) for the remaining entries.
\end{proof}
\subsection{Initialisation of the Expectation-Maximisation algorithm}
The first step in the initialisation involves computing a first approximation for the trends. This is achieved via univariate trend-cycle decompositions. In the case of headline and core inflation, the initialisation of the trend involves a further operation. Trend inflation is initialised by taking the mean between the persistent components estimated for headline and core inflation, appropriately rescaled by $\eta_8$ and $\eta_9$. The variances of the innovations are calculated on the double differenced initial trends.
The second step involves the initialisation of the cycles, which is performed on the de-trended data. The business cycle is approximated by the first principal component of the de-trended data and a series of ridge regressions is used for computing the coefficients of the cycles. The restrictions described in (ref) are enforced on each regression. The variances of the innovations are computed on the sample residuals.
\subsection{Enforcing causality during the estimation}
The ECM algorithm used in this manuscript ensures that the AR states (i.e., the common cycle and idiosyncratic noise components) are causal at every iteration. This is achieved with the approach proposed in pellegrino2020selecting for vector autoregressions.
\subsection{Hyperparameter selection}
The hyperparameters are selected using the artificial jackknife selection method proposed in pellegrino2020selecting. In the empirical application in (ref), the grid of candidate hyperparameters $\mathscr{H} = \mathscr{H}_p \times \mathscr{H}_\lambda \times \mathscr{H}_\alpha \times \mathscr{H}_\beta$ is such that $\mathscr{H}_p \vcentcolon= \{12\}$, $\mathscr{H}_\lambda \vcentcolon= [10^{-2}, 2.5]$, $\mathscr{H}_\alpha \vcentcolon= [0, 1]$ and $\mathscr{H}_\beta \vcentcolon= [1, 1.2]$. The selection process returns the specification with the lowest expected forecast error for headline inflation, following a rational similar to jarocinski2018inflation.\footnote{The weights in pellegrino2020selecting are set to be equal to zero for all variables with the exception of headline inflation which has weight equal to one.}
\section{Additional material}
\begin{algorithm}[H]
\textsf{Initialization}\\
The ECM algorithm is initialised as described in (ref).
\textsf{Estimation}\\
\For{$k \gets 1$ to $max\_iter$}{
\For{$j \gets 1$ to $m$}{
Run the Kalman filter and smoother using $\boldsymbol{\mathbf{\hat{\vartheta}}}^{k-1}_s(\boldsymbol{\mathbf{\gamma}})$\;
\If{converged}{
Store the parameters and stop the loop.
}
Estimate $\boldsymbol{\mathbf{\hat{\mu}}}_{s,0}^{k}(\boldsymbol{\mathbf{\gamma}})$ and $\boldsymbol{\mathbf{\hat{\Omega}}}_{s,0}^{k}(\boldsymbol{\mathbf{\gamma}})$ as in (ref)\;
Estimate $\boldsymbol{\mathbf{\hat{C}}}^{k}_{s}(\boldsymbol{\mathbf{\gamma}})$, $\boldsymbol{\mathbf{\hat{\Sigma}}}^{k}_{s}(\boldsymbol{\mathbf{\gamma}})$ and $\boldsymbol{\mathbf{\hat{B}}}^{k}_{s}(\boldsymbol{\mathbf{\gamma}})$ as in \crefrange{ch2:lemma:ecm_c_star}{ch2:lemma:ecm_loadings}\;
Build $\boldsymbol{\mathbf{\hat{\vartheta}}}^{k}_s(\boldsymbol{\mathbf{\gamma}})$\;
}
}
\textbf{Notes} \begin{itemize}
• The results are computed fixing $max\_iter$ to 1000. This is a conservative number, since the algorithm generally requires substantially less iterations to converge.
• The ECM algorithm is considered to be converged when the estimated coefficients (all relevant parameters in \crefrange{ch2:lemma:ecm_c_star}{ch2:lemma:ecm_loadings}) do not significantly change in two subsequent iterations. This is done by computing the absolute relative change per parameters and comparing at the same time the median and $95^{th}$ quantile respectively with a fixed tolerance of $10^{-3}$ and $10^{-2}$. Intuitively, when the coefficients do not change much, the expected log-likelihood and the parameters in (ref) should also be stable.
• The scalar $\varepsilon$ is summed to the denominator of each relative change in order to ensure numerical stability.
\end{itemize}
\caption{ECM algorithm for the trend-cycle decomposition}
\end{algorithm}
The replication code for this paper is available on \href{https://github.com/fipelle/replication-pellegrino-2022-ensembles}{GitHub}.
\begin{figure}[!h]
\caption{Importance weights pre COVID-19: top 10 predictors. Pre COVID-19 weights are computed using the macroeconomic series available on the 28th February 2020 on ALFRED and the corresponding Wilshire data.}
\end{figure}
\begin{figure}[!h]
\caption{Importance weights post COVID-19: top 10 predictors.}
\end{figure}
\begin{table}[!h]
\caption{Glossary for the acronyms in (ref).}
\begin{tabularx}{\textwidth}{@ll@}
\toprule
Acronym & Description \\
\midrule
BEA & Bureau of Economic Analysis \\
BLS & Bureau of Labor Statistics \\
CPI & Consumer Price Index \\
FRB & Federal Reserve Board \\
FRBSL & Federal Reserve Bank of St. Louis \\
TMI & Total Market Index \\
WA & Wilshire Associates \\
\bottomrule
\end{tabularx}
\end{table}
\begin{table}[!h]
\caption{Mean squared errors associated to ridge forecasts, relative to those for naive predictions constant at zero. This is an exact replica of (ref), except for the forecasting model. Indeed, this table uses a ridge regression (with shrinkage selected in an equally spaced grid with cardinality 25 ranging from 0.1 to 100) instead of bootstrap aggregating. }
\begin{tabularx}{\textwidth}{@l ZZ c ZZ@}
\toprule
Target & \multicolumn{2}{c}{Pre COVID-19} & & \multicolumn{2}{c}{Post COVID-19} \\
\cmidrule(l){2-3} \cmidrule(l){5-6}
& Autoregressive & Augmented & & Autoregressive & Augmented \\
\midrule
WILL5000IND & 0.983 & 1.041 & & 0.960 & 3.378 \\
WILLLRGCAP & 0.986 & 1.064 & & 0.960 & 3.374 \\
WILLLRGCAPVAL & 0.944 & 0.916 & & 0.938 & 2.843 \\
WILLLRGCAPGR & 1.070 & 1.304 & & 1.015 & 3.623 \\
WILLMIDCAP & 0.996 & 1.062 & & 0.972 & 2.436 \\
WILLMIDCAPVAL & 0.964 & 0.936 & & 0.966 & 1.770 \\
WILLMIDCAPGR & 1.081 & 1.344 & & 1.040 & 3.375 \\
WILLSMLCAP & 0.988 & 1.065 & & 0.973 & 2.525 \\
WILLSMLCAPVAL & 0.946 & 0.884 & & 0.944 & 1.976 \\
WILLSMLCAPGR & 1.075 & 1.459 & & 1.033 & 3.070 \\
\bottomrule
\end{tabularx}
\end{table}