Multidimensional dynamic factor models

abstractThis paper generalises dynamic factor models for multidimensional dependent data. In doing so, it develops an interpretable technique to study complex information sources ranging from repeated surveys with a varying number of respondents to panels of satellite images. We specialise our results to model microeconomic data on US households jointly with macroeconomic aggregates. This results in a powerful tool able to generate localised predictions, counterfactuals and impulse response functions for individual households, accounting for traditional time-series complexities depicted in the state-space literature. The model is also compatible with the growing focus of policymakers for real-time economic analysis as it is able to process observations online, while handling missing values and asynchronous data releases.

keywordDynamic Factor Models, Heterogeneous Agent Models, Multidimensional Time Series, Real-time Economic Analysis, State-space Models

frontmatter\address[add1]{University of Bologna} \address[add2]{Imperial College London} \address[add3]{London School of Economics and Political Science} \cortext[cor1]{Corresponding author: [email removed].}

assumption[Data] Let $\boldsymbol{\mathbf{H}}_t \in \mathbb{R}^{N_t \times K}$ be a data matrix with $N_{t} > 0$ and $K>0$, and denote with $\boldsymbol{\mathbf{h}}_{t}$ the $N_t K \times 1$ vectorisation of $\boldsymbol{\mathbf{H}}_{t}$, for every point in time $t$. Besides, assume that $\boldsymbol{\mathbf{H}}_t$ is a finite realisation of some stochastic process observed at any point in time $t \in \mathscr{T} \subseteq \{1, \ldots, T\}$ where $T \geq 1$.

remarkIn our notation, $N_{t}$ denotes the number of subjects at each point in time. It is important to stress that we talk about “subjects” figuratively. Indeed, our definition is not restricted to individuals, but extends to any abstract thing with a data structure compatible with (ref).

example[Time-series cross sections] In the case of time-series cross-sectional data \begin{align*} \boldsymbol{\mathbf{H}}_t = \begin{pmatrix} H_{1,1,t} & \ldots & H_{1,K,t} \\ \vdots & \ddots & \vdots \\ H_{N_t,1,t} & \ldots & H_{N_t,K,t} \end{pmatrix} \end{align*} represents a cross section referring to time $t$ and \begin{align*} \boldsymbol{\mathbf{h}}_t = \begin{pmatrix} H_{1, 1, t} & \ldots & H_{1, K, t} & \ldots & H_{N_t, 1, t} & \ldots & H_{N_t, K, t} \end{pmatrix}', \end{align*} where $N_t > 0$ is the number of cross-sectional observations for time $t$ and $K$ is the total number of covariates. In social sciences, similar datasets generally represent complex surveys with a varying number of respondents. However, $\boldsymbol{\mathbf{H}}_{t}$ could also represent more exotic data. For instance, the RGB representation of a satellite image taken at time $t$ with $N_{t}$ pixels.

example[Time series] The banal case in which $\boldsymbol{\mathbf{H}}_{t} \in \mathbb{R}^{N_{t}}$ gives a time series dataset. Indeed, $\boldsymbol{\mathbf{H}}_t = (H_{1,t} \;\, \ldots, H_{N_{t},t})'$ for any $t \in \mathscr{T}$. The value taken by $N_{t}$ over all $t \in \mathscr{T}$ controls whether this dataset represents a univariate or a multivariate time series, and if it is fully observed.

definition[Identifiers] In order to allow for this complexity, we let \begin{align*} \mathscr{S}_t \vcentcolon= \{ f(i,t) : 1 \leq i \leq N_t K\}, \end{align*} for any $t \in \mathscr{T}$. The function $f:\mathbb{N} \times \mathbb{N} \rightarrow \mathbb{N}$ is a convenient way for categorising different subjects. Indeed, we structure it to be equal to one when evaluated at $(1,1)$ and to have an incremental value for any pair referring to a new feature of the same subject, or to a new subject. This implies that $f(i_1, t_1)=f(i_2, t_2)$ if and only if $(i_1, t_1)$ and $(i_2, t_2)$ refer to the same subject-characteristic pair. Hence, we let \begin{align*} \mathscr{S} \vcentcolon= \bigcup_{t \in \mathscr{T}} \mathscr{S}_t \end{align*} be the set of all (observed) subject-characteristic pairs. For simplicity, we let $N \vcentcolon= \frac{|\mathscr{S}|}{K}$ be the number of unique subjects.

remarkUnder (ref), a minimum of one subject is observed across the whole sample and thus $N > 0$. Besides, note that $N$ is a natural number by construction.

definitionIndeed, we let $\boldsymbol{\mathbf{Y}}_{t} \vcentcolon= \boldsymbol{\mathbf{W}}_{t} \, \boldsymbol{\mathbf{h}}_{t}$ to be a $NK \times 1$ vector of time series where $\boldsymbol{\mathbf{W}}_{t}$ is a $NK \times N_{t} K$ matrix such that \begin{align*} W_{i, j, t} = \begin{cases} 1, & if f(j,t)=i and i \in \mathscr{S}_{t}, \\ 0, &otherwise, \end{cases} \end{align*} for any $t \in \mathscr{T}$, $1 \leq i \leq NK$ and $1 \leq j \leq N_{t} K$. In order to have a record of the non-missing entries of each $i$-th subject-characteristic pair, we also let $\mathscr{T}_{i} \subseteq \mathscr{T}$ to be the set of points in time in which it is observed. Interestingly, $\mathscr{T} = \bigcup_{i=1}^N \mathscr{T}_{i}$.

remarkNote that by taking track of the observed datapoints via the $\mathscr{T}_{i}$ sets, we can distinguish between real zeros and missing values. Indeed, we do so in the estimation method described below. Moreover, due to (ref), all $\mathscr{T}_{i}$ referring to a single subject are identical (i.e., when we observe a subject we measure all of its characteristics).

assumption[Generic MDFM] Going forward, we assume that the model for any $\boldsymbol{\mathbf{Y}}_{t}$ is the multidimensional dynamic factor model \begin{alignat*}{2} \boldsymbol{\mathbf{Y}}_{t} &= \boldsymbol{\mathbf{B}}(L) \boldsymbol{\mathbf{\Phi}}_{t} + \boldsymbol{\mathbf{e}}_{t}, \qquad &&\boldsymbol{\mathbf{e}}_{t} \mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}} N(\boldsymbol{\mathbf{0}}_{NK \times 1}, \boldsymbol{\mathbf{R}}), \\ \boldsymbol{\mathbf{\Phi}}_{t} &= \boldsymbol{\mathbf{C}}(L) \boldsymbol{\mathbf{\Phi}}_{t-1} + \boldsymbol{\mathbf{D}} \boldsymbol{\mathbf{u}}_{t}, \qquad &&\boldsymbol{\mathbf{u}}_{t} \mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}} N(\boldsymbol{\mathbf{0}}_{r \times 1}, \boldsymbol{\mathbf{\Sigma}}), \end{alignat*} where $\boldsymbol{\mathbf{\Phi}}_{t}$ denotes a vector of latent components (stationary and/or non-stationary) and for some positive definite covariance matrices $\boldsymbol{\mathbf{R}}$ and $\boldsymbol{\mathbf{\Sigma}}$. The vector $\boldsymbol{\mathbf{\Phi}}_{t}$ is $q$-dimensional with $1 \leq q \ll NK$, linked to the measurements via the matrix $\boldsymbol{\mathbf{B}}(L)$ and with dynamics determined by $\boldsymbol{\mathbf{C}}(L)$. Besides, $\boldsymbol{\mathbf{D}}$ is $q \times r$ with $1 \leq r \leq q$.

remarkIn this article, the common components are not restricted to be stationary.

definition[Groups] Define $\mathscr{G}$ as the Cartesian product of the household attributes on education and ethnicity in (ref): a set with cardinality four such that each member is a unique combination of characteristics that identifies a specific demographic group. For simplicity, we refer to these groups in the order: (0, 0), (0, 1), (1, 0), (1, 1) whereas zero and one refers to the values taken by the binary variables EDUC_HH and WHITE_HH. Finally, we also let $\boldsymbol{\mathbf{\omega}} \equiv \boldsymbol{\mathbf{\omega}}(\mathscr{G})$ be the vector of integers denoting the number of households per group observed across all periods.

definition[Empirical data] We then arrange the data to match the structure in (ref) and let \begin{align*} \boldsymbol{\mathbf{Y}}_{t} = \big(\boldsymbol{\mathbf{X}}_{t}' \;\, \boldsymbol{\mathbf{Z}}_{1,t}' \;\, \ldots \;\, \boldsymbol{\mathbf{Z}}_{4,t}' \big)', \end{align*} where $\boldsymbol{\mathbf{X}}_{t}$ denotes the vector of macroeconomic aggregates and each $\boldsymbol{\mathbf{Z}}_{i,t}$ represents the vector of real income per head for all households in group $1 \leq i \leq 4$.

remarkNote that every $\boldsymbol{\mathbf{Z}}_{i,t}$ is $\omega_i \times 1$ dimensional. Since the CE PUMD is structured to survey the same household for a maximum of 4 quarters, the $\boldsymbol{\mathbf{Z}}_{i,t}$ vectors are sparse. Besides, the missing observations in $\boldsymbol{\mathbf{Y}}_{t}$ are handled as in (ref). Finally, recall that this application focusses on the four demographic groups indicated in (ref).

assumptionFormally, we let \begin{align*} \left( \begin{array}{c} X_{1,t} \\ X_{2,t} \\ \vdots \\ X_{8,t} \\ \boldsymbol{\mathbf{Z}}_{1,t} \\ \boldsymbol{\mathbf{Z}}_{2,t} \\ \boldsymbol{\mathbf{Z}}_{3,t} \\ \boldsymbol{\mathbf{Z}}_{4,t} \end{array} \right) &= \left( \begin{array}{c} \tau_{1,t} \\ \tau_{2,t} \\ \vdots \\ \tau_{8,t} \\ \tau_{9,t} \, \boldsymbol{\mathbf{\iota}}_{\omega_1} \\ \tau_{10,t} \, \boldsymbol{\mathbf{\iota}}_{\omega_2} \\ (\tau_{9,t} + \tau_{11,t}) \, \boldsymbol{\mathbf{\iota}}_{\omega_3} \\ (\tau_{10, t} + \tau_{11, t}) \, \boldsymbol{\mathbf{\iota}}_{\omega_4} \end{array} \right) + \left( \begin{array}{cccc} 1 \\ \sum_{i=1}^{p} \Lambda_{1,i} L^{i-1} \\ \vdots \\ \sum_{i=1}^{p} \Lambda_{7,i} L^{i-1} \\ \sum_{i=1}^{p} \Lambda_{8,i} \, \boldsymbol{\mathbf{\iota}}_{\omega_1} L^{i-1} \\ \sum_{i=1}^{p} \Lambda_{9,i} \, \boldsymbol{\mathbf{\iota}}_{\omega_2} L^{i-1} \\ \sum_{i=1}^{p} \Lambda_{10,i} \, \boldsymbol{\mathbf{\iota}}_{\omega_3} L^{i-1} \\ \sum_{i=1}^{p} \Lambda_{11,i} \, \boldsymbol{\mathbf{\iota}}_{\omega_4} L^{i-1} \\ \end{array} \right) \psi_{t} + \left( \begin{array}{c} \xi_{1,t} \\ \xi_{2,t} \\ \vdots \\ \xi_{8,t} \\ \xi_{9,t} \\ \xi_{10,t} \\ \xi_{11,t} \\ \xi_{12, t} \end{array} \right) + \boldsymbol{\mathbf{e}}_{t} \end{align*} where $\psi_{t}$ is a causal AR($p$) cycle denoting the business cycle; the $\tau$ denote smooth trends of order two modelled as in kitagawa1996smoothness; $\xi_{1,t}, \ldots, \xi_{8+|\mathscr{G}|, t}$ are causal AR(1) latent components representing idiosyncratic noise; $\boldsymbol{\mathbf{\iota}}$ denotes a vector of ones with length indicated in the subscript. Hereinafter, the number of lags $p$ is assumed being equal to $4$ (quarters).

remark[Trends] Recall that a generic smooth trend $\underline{\tau}$ modelled as in kitagawa1996smoothness is of order $k$ if $(1-L)^k \, \underline{\tau}$ is a white noise. Besides, note that the income figures share common trends. In particular: $\boldsymbol{\mathbf{\tau}}_{9}$ models the persistent part of income for not college educated, not white households; $\boldsymbol{\mathbf{\tau}}_{10}$ models the persistent part of income for not college educated, white households; $\boldsymbol{\mathbf{\tau}}_{11}$ models the persistent offset of college educated households.

remark[CE PUMD data] (ref) implies that each household is modelled as a function of its own group and the dedicated parameters. In other words, all members of the $i$-th group are modelled via the same set of coefficients and latent factors, for every $1 \leq i \leq 4$. While the generic structure proposed in (ref) could allow for a more disaggregate model, we do not have enough observations in the CE PUMD to do it. That being said, the model in (ref) has quite a few advantages compared to these granular theoretical alternatives. Most importantly, it is less subject to idiosyncratic noise and due to the dimensionality reduction into group factors it is easier to interpret.

subappendices\section{CE PUMD} The demographic characteristics in (ref) are constructed at the household level. The following paragraphs give further details on each variable. \begin{itemize} • Census region: Census Bureau classification for US regions (1 Northeast, 2 Midwest, 3 South, 4 West). We use it for excluding CUs that moved across the US during the sampling period. • College educated households: describes the highest level of education of the reference person and spouse (if any). It is a dummy variable equal to 1 for CUs in which the highest education level is, at least, at an undergraduate level and 0 otherwise. • Family size: Number of family members. We use it for computing real household income per head (before tax). • Family type: Family categorisation. We use it for determine whether we there is a spouse to consider when constructing the other variables in this appendix. • Prime working age: dummy variable equal to 1 for CUs with average age between 25 and 55 years (excluded) and 0 otherwise. The average age is computed by taking the sample mean between the age of the reference person and spouse (if any). We use it for excluding non prime working age households. • Urban consumers: dummy variable equal to 1 for urban CUs and 0 otherwise. We use it for excluding rural CUs. • \emph{White household}: dummy variable equal to 1 for white CUs and 0 otherwise. \end{itemize} The nominal income per head (before tax) is computed by constructing total nominal income from the ITBI files and dividing it for the number of CUs members in the FMLI files. The identifiers or Universal Classification Code (UCC) for each single income component used for computing this total are summarised in (ref). \begin{table}[!t] \begin{tabularx}{\textwidth}{@lZ@} \toprule Releases & Universal Classification Codes (UCCs) \\ \midrule \multirow{3}{*}{1990 to 2003} & 900000, 900010, 900020, 900030, 900040, 900080, 900050, 900060, \\ & 900070, 900100, 900110, 900090, 900120, 900150, 900131, 900132, \\ & 800700I, 800710I, 900140 \\ \midrule \multirow{3}{*}{2004 to 2012} & 900000, 900010, 900020, 900030, 900040, 900080, 900050, 900060, \\ & 900070, 900100, 900110, 900090, 900120, 900150, 900131, 900132, \\ & 800700, 800710, 900140 \\ \midrule \multirow{2}{*}{2013 to 2020} & 900000, 900160, 900030, 900170, 900180, 900190, 900200, \\ & 900090, 900120, 900150, 900210, 800700, 800710, 900140 \\ \bottomrule \end{tabularx} \caption{Universal Classification Codes (UCCs) used for computing nominal income. \\ \textbf{Source:} Census Bureau for the Bureau of Labor Statistics.} \end{table} \section{ECM algorithm} This appendix develops an ECM algorithm meng1993maximum to estimate the MDFM in (ref). The design builds on pellegrino2020selecting and pellegrino2021factoraugmented. This manuscript uses the “hat” symbol to denote the estimated coefficients, an $s$ subscript to indicate the sample size and a $k$ superscript for the ECM iteration. \subsection{State-space representation} Recall that \begin{alignat*}{2} \boldsymbol{\mathbf{Y}}_{t} &= \boldsymbol{\mathbf{B}}(L) \boldsymbol{\mathbf{\Phi}}_{t} + \boldsymbol{\mathbf{e}}_{t}, \qquad &&\boldsymbol{\mathbf{e}}_{t} \mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}} N(\boldsymbol{\mathbf{0}}_{NK \times 1}, \boldsymbol{\mathbf{R}}), \\ \boldsymbol{\mathbf{\Phi}}_{t} &= \boldsymbol{\mathbf{C}}(L) \boldsymbol{\mathbf{\Phi}}_{t-1} + \boldsymbol{\mathbf{D}} \boldsymbol{\mathbf{u}}_{t}, \qquad &&\boldsymbol{\mathbf{u}}_{t} \mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}} N(\boldsymbol{\mathbf{0}}_{r \times 1}, \boldsymbol{\mathbf{\Sigma}}). \end{alignat*} The matrices $\boldsymbol{\mathbf{C}}(L)$ and $\boldsymbol{\mathbf{D}}$ are sparse and their non-zero entries are such that \begin{align*} \boldsymbol{\mathbf{C}}(L) = &\left( \begin{array}{c | ccc | ccccc | c} 2 \boldsymbol{\mathbf{I}}_{7+|\mathscr{G}|} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & - \boldsymbol{\mathbf{I}}_{7+|\mathscr{G}|} \\ \hline \cdot & \pi_{1} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \pi_{8+|\mathscr{G}|} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \hline \cdot & \cdot & \cdot & \cdot & \pi_{8+|\mathscr{G}|+1} & \pi_{8+|\mathscr{G}|+2} & \ldots & \pi_{8+|\mathscr{G}|+p-1} & \pi_{8+|\mathscr{G}|+p} & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot \\ \hline \boldsymbol{\mathbf{I}}_{7+|\mathscr{G}|} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \end{array} \right) \\[-1em] & \begin{array}{c ccc ccccc c} \multicolumn{1}{c}{\underbrace_{q \times 7+|\mathscr{G}|}} & \multicolumn{3}{c}{\underbrace_{q \times 8+|\mathscr{G}|}} & \multicolumn{5}{c}{\underbrace_{q \times p}} & \multicolumn{1}{c}{\underbrace_{q \times 7+|\mathscr{G}|}} \end{array} \end{align*} and \begin{align*} \boldsymbol{\mathbf{D}} = &\left( \begin{array}{c | c | c } \boldsymbol{\mathbf{I}}_{7+|\mathscr{G}|} & \cdot & \cdot \\ \hline \cdot & \boldsymbol{\mathbf{I}}_{8+|\mathscr{G}|} & \cdot \\ \hline \cdot & \cdot & \phantom{1} 1 \phantom{1} \\ \hline \cdot & \cdot & \cdot \end{array} \right) \\[-1em] & \begin{array}{c c c} \multicolumn{1}{c}{\underbrace_{q \times 7+|\mathscr{G}|}} & \multicolumn{1}{c}{\underbrace_{q \times 8+|\mathscr{G}|}} & \multicolumn{1}{c}{\underbrace_{q \times 1}} \end{array} \end{align*} where $\boldsymbol{\mathbf{\pi}}$ is a $8+|\mathscr{G}|+p \times 1$ vector of finite real parameters which ensures that the cyclical components are causal. Due to the structure of $\boldsymbol{\mathbf{C}}$ and $\boldsymbol{\mathbf{D}}$, it follows that $r = 16+2|\mathscr{G}|$ and $q=22+3|\mathscr{G}|+p$. The measurement coefficient matrix $\boldsymbol{\mathbf{B}}(L)$ is also sparse and its non-zero entries are \begin{footnotesize} \begin{align*} &\left( \begin{array}{ccccccc | ccccccc | cccc | c} 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \phantom \cdot \phantom & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \ldots & \cdot & \cdot \\ \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \Lambda_{1,1} & \Lambda_{1,2} & \ldots & \Lambda_{1,p} & \cdot \\ \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \vdots & \vdots & \vdots & \vdots & \cdot \\ \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \Lambda_{7,1} & \Lambda_{7,2} & \ldots & \Lambda_{7,p} & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \Lambda_{8,1} \, \boldsymbol{\mathbf{\iota}}_{\omega_{1}} & \Lambda_{8,2} \, \boldsymbol{\mathbf{\iota}}_{\omega_{1}} & \ldots & \Lambda_{8,p} \, \boldsymbol{\mathbf{\iota}}_{\omega_{1}} & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \ddots & \cdot & \vdots & \vdots & \vdots & \vdots & \phantom{1} \cdot \phantom{1} \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \Lambda_{7+|\mathscr{G}|,1} \, \boldsymbol{\mathbf{\iota}}_{\omega_{|\mathscr{G}|}}& \Lambda_{7+|\mathscr{G}|,2} \, \boldsymbol{\mathbf{\iota}}_{\omega_{|\mathscr{G}|}} & \ldots & \Lambda_{7+|\mathscr{G}|,p} \, \boldsymbol{\mathbf{\iota}}_{\omega_{|\mathscr{G}|}} & \cdot \end{array} \right).\\[-1em] & \begin{array}{ccccccc ccccccc cccc c} \multicolumn{7}{c}{\underbrace_{NK \times 7+|\mathscr{G}|}} & \multicolumn{7}{c}{\underbrace_{NK \times 8+|\mathscr{G}|}} & \multicolumn{4}{c}{\underbrace_{NK \times p}} & \multicolumn{1}{c}{\underbrace_{NK \times 7+|\mathscr{G}|}} \end{array} \end{align*} \end{footnotesize} As a result, \begin{align*} &\boldsymbol{\mathbf{\Phi}}_{t} \vcentcolon= \left( \begin{array}{ccc | ccc | cccc | ccc} \tau_{1,t} & \ldots & \tau_{7+|\mathscr{G}|,t} & \xi_{1,t} & \ldots & \xi_{8+|\mathscr{G}|, t} & \psi_{t} & \psi_{t-1} & \ldots & \psi_{t-p+1} & \tau_{1,t-1} & \ldots & \tau_{7+|\mathscr{G}|,t-1} \end{array} \right)'. \\[-1em] & \begin{array}{ccc ccc cccc ccc} \multicolumn{3}{c}{\underbrace_{7+|\mathscr{G}| \times 1}} & \multicolumn{3}{c}{\underbrace_{8+|\mathscr{G}| \times 1}} & \multicolumn{4}{c}{\underbrace_{p \times 1}} & \multicolumn{3}{c}{\underbrace_{7+|\mathscr{G}| \times 1}} \end{array} \end{align*} \begin{assumption}[Initial conditions] Given that we observe data at time $t=1$, we further assume that $\boldsymbol{\mathbf{\Phi}}_{0} \mathrel{\overset{w.n.}{\scalebox{1.5}[1]{$\sim$}}} N(\boldsymbol{\mathbf{\mu}}_{0}, \boldsymbol{\mathbf{\Omega}}_{0})$ for some finite vector of real parameters $\boldsymbol{\mathbf{\mu}}_{0}$ and a positive definite covariance matrix $\boldsymbol{\mathbf{\Omega}}_{0}$. The latter is assumed to be sparse and such that the only entries allowed to differ from zero are those denoting the initial auto-covariances of each state. \end{assumption} \begin{remark}[Non-zero entries of $\boldsymbol{\mathbf{\Omega}}_{0}$] In other words, the entries of $\boldsymbol{\mathbf{\Omega}}_{0}$ that are allowed to differ from zero are those with coordinates $(i,j)$ in the union of the following sets: \begin{itemize} • $\{(i,j) : i=j \text{ and } 1 \leq i < r\}$; • $\{(i,j) :r \leq i < r+p \text{ and } r \leq j < r+p\}$. \end{itemize} \end{remark} \subsection{Estimation} This manuscript builds on the theoretical results in barigozzi2020quasi and estimates the model via quasi penalised maximum likelihood estimation (PQMLE) by considering $\boldsymbol{\mathbf{\Sigma}}$ as a diagonal matrix and $\boldsymbol{\mathbf{R}} = \varepsilon \, \boldsymbol{\mathbf{I}}_{NK}$ for a small positive $\varepsilon$.\footnote{We set $\varepsilon = 10^{-2}$.} Formally, this implies that the free parameters to estimate are \begin{align*} &\boldsymbol{\mathbf{\vartheta}} \vcentcolon= \left( \begin{array}{cccccccc} \boldsymbol{\mathbf{\mu}}_0' & \text{vech}(\boldsymbol{\mathbf{\Omega}}_0)' & \text{vec}(\boldsymbol{\mathbf{\Lambda}})' & \boldsymbol{\mathbf{\pi}}' & \Sigma_{1,1} & \Sigma_{2,2} & \ldots & \Sigma_{r, r} \end{array} \right)'. \end{align*} The estimation is performed with an ECM algorithm: an optimisation method that repeats the operations in (ref) until it reaches 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}_{\underline{\boldsymbol{\mathbf{\vartheta}}} \, \in \, \mathscr{R}} \, \mathbb{E} \left [\mathcal{L}(\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \boldsymbol{\mathbf{Y}}_{1:s}, \boldsymbol{\mathbf{\Phi}}_{1:s}) \,|\, \mathscr{Y}(s), \, \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] - \mathbb{E} \left [\mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}) \,|\, \mathscr{Y}(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{Y}(s)$ is the information set available at time $s$, \begin{align} \mathcal{L}(\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \boldsymbol{\mathbf{Y}}_{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}}_{1:r, t} - \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\Phi}}_{t-1})(\boldsymbol{\mathbf{\Phi}}_{1:r, 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{Y}}_t - \underline{\boldsymbol{\mathbf{B}}}\boldsymbol{\mathbf{\Phi}}_{t})(\boldsymbol{\mathbf{Y}}_t - \underline{\boldsymbol{\mathbf{B}}}\boldsymbol{\mathbf{\Phi}}_{t})'\Big], \nonumber \end{align} $\underline{\boldsymbol{\mathbf{C}}}_{\,*} \equiv \underline{\boldsymbol{\mathbf{C}}}_{\,1:r, 1:q}$ and the underlined coefficients denote the parameters implied by $\underline{\boldsymbol{\mathbf{\vartheta}}}$. The function in (ref) is the so-called complete-data (i.e., fully observed data and known latent states) log-likelihood. 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:8+|\mathscr{G}|} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1)^{\frac{1}{2}} \big\Vert_{\text{F}}^2 + \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,8+|\mathscr{G}|+1:8+|\mathscr{G}|+p}' \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p)^{\frac{1}{2}} \big\Vert_{\text{F}}^2 + \big\Vert \underline{\boldsymbol{\mathbf{\Lambda}}} \, \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:8+|\mathscr{G}|} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1) \big\Vert_{1,1} + \big\Vert \underline{\boldsymbol{\mathbf{\pi}}}_{\,8+|\mathscr{G}|+1:8+|\mathscr{G}|+p}' \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p) \big\Vert_{1,1} + \big\Vert \underline{\boldsymbol{\mathbf{\Lambda}}} \, \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p) \big\Vert_{1,1} \right) \end{align*} is a version of the elastic-net penalty in pellegrino2020selecting in which, for any $l \in \mathbb{N}$, \begin{align*} \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, l) \vcentcolon= \rho \begin{pmatrix} 1 & 0 &\ldots & 0 \\ 0 &\beta &\ldots & 0 \\ \vdots &\ddots &\ddots &\vdots \\ 0 &\ldots &\ldots & \beta^{l-1} \end{pmatrix}, \end{align*} $\rho \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{assumption}[Convergence] The ECM algorithm is considered to be converged when the estimated coefficients 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 with a fixed tolerance of $10^{-3}$ and $10^{-2}$ respectively. \end{assumption} The operation in (ref) is performed in two steps: the so-called E-step and CM-step. The E-step computes the expectations in (ref), while the CM-step conditionally maximises them with respect to the free parameters. We write down the E-step on the basis of the output of a Kalman smoother compatible with incomplete data, as originally proposed in shumway1982approach and watson1983alternative. For that, we use the notation in (ref). \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{Y}(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{Y}(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} We also use the notation in (ref) to further deal with missing observations. \begin{definition}[Observed measurements] Recall that $\mathscr{T} = \bigcup_{i=1}^N \mathscr{T}_i$ and let \begin{align*} &\mathscr{T}(s) \vcentcolon= \{t : t \in \mathscr{T}, \, 1 \leq t \leq s\}, \end{align*} for $1 \leq s \leq T$. Let also \begin{align*} \mathscr{D}_{t} \vcentcolon= \{i : t \in \mathscr{T}_i, \, 1 \leq i \leq NK\}, \end{align*} for $1 \leq t \leq T$. Finally, let \begin{align*} &\boldsymbol{\mathbf{Y}}_{t}^{obs} \vcentcolon= \big( Y_{i,t} \big)_{i \in \mathscr{D}_t} \\ &\boldsymbol{\mathbf{B}}_{t}^{obs} \vcentcolon= \boldsymbol{\mathbf{A}}_{t} \boldsymbol{\mathbf{B}} \end{align*} be the $|\mathscr{D}_t| \times 1$ vector of observed measurements at time $t$ and the corresponding $|\mathscr{D}_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} Moreover, in order to simplify the notation, we let \begin{align*} \mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \equiv \mathbb{E} \left [\mathcal{L}(\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \boldsymbol{\mathbf{Y}}_{1:s}, \boldsymbol{\mathbf{\Phi}}_{1:s}) \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \end{align*} for $1 \leq s \leq T$. It then follows directly from pellegrino2021factoraugmented that \begin{align*} \mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Y}(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} \big( \boldsymbol{\mathbf{\hat{\Phi}}}_0 \boldsymbol{\mathbf{\hat{\Phi}}}_0' + \boldsymbol{\mathbf{\hat{P}}}_0 - \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) \Big] \\ &-\frac{s}{2}\ln|\underline{\boldsymbol{\mathbf{\Sigma}}}| - \frac{1}{2} \operatorname{Tr} \Big[\underline{\boldsymbol{\mathbf{\Sigma}}}^{-1} \big( \boldsymbol{\mathbf{\hat{E}}}^{(1)}_s - \boldsymbol{\mathbf{\hat{E}}}^{(2)}_s \underline{\boldsymbol{\mathbf{C}}}_{\,*}' - \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\hat{E}}}^{(2)'}_s + \underline{\boldsymbol{\mathbf{C}}}_{\,*} \boldsymbol{\mathbf{\hat{E}}}^{(3)}_s \underline{\boldsymbol{\mathbf{C}}}_{\,*}' \big) \Big] \\ &-\frac{1}{2 \varepsilon} \operatorname{Tr} \left\{ \sum_{t \in \mathscr{T}(s)} \left[ \big( \boldsymbol{\mathbf{Y}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \big) \big( \boldsymbol{\mathbf{Y}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \big)' + \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{P}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{{obs}'} \right] \right\}, \end{align*} where \begin{alignat*}{3} &\boldsymbol{\mathbf{\hat{E}}}^{(1)}_s &&\vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{1:r, t} \boldsymbol{\mathbf{\Phi}}_{1:r, t}' \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] &&= \sum_{t=1}^s \left(\boldsymbol{\mathbf{\hat{\Phi}}}_t \boldsymbol{\mathbf{\hat{\Phi}}}_t' + \boldsymbol{\mathbf{\hat{P}}}_t\right)_{1:r, 1:r}, \\ &\boldsymbol{\mathbf{\hat{E}}}^{(2)}_s &&\vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{1:r, t} \boldsymbol{\mathbf{\Phi}}_{t-1}' \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \Big] &&= \sum_{t=1}^s \left(\boldsymbol{\mathbf{\hat{\Phi}}}_t \boldsymbol{\mathbf{\hat{\Phi}}}_{t-1}' + \boldsymbol{\mathbf{\hat{P}}}_{t, t-1}\right)_{1:r, 1:q}, \\ &\boldsymbol{\mathbf{\hat{E}}}^{(3)}_s &&\vcentcolon= \sum_{t=1}^s \mathbb{E}\Big[\boldsymbol{\mathbf{\Phi}}_{t-1} \boldsymbol{\mathbf{\Phi}}_{t-1}' \,|\, \mathscr{Y}(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{alignat*} for $1 \leq s \leq T$. Furthermore, it follows from pellegrino2021factoraugmented that \begin{align*} \mathbb{E} \left[ \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}) \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] = \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}). \end{align*} The CM-step conditionally maximises the expected penalised log-likelihood \begin{align} \mathcal{M}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}} \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] \vcentcolon= \mathcal{L}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}} \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right] - \mathcal{P}(\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}}) \end{align} with respect to the free parameters. The CM-steps for all free parameters of the transition equation are reported in pellegrino2021factoraugmented. For clarity, we recall just the lemmas' statements in \crefrange{lemma:ecm_init}{lemma:mstep_sigma} with the adequate minimal notational changes. \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$, denoted with $\boldsymbol{\mathbf{\hat{\Omega}}}_{0,s}^{k+1}(\boldsymbol{\mathbf{\gamma}})$, is a sparse covariance matrix whose entries allowed to differ from zero 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*} and the coordinates $(i,j)$ are those described in (ref). \end{lemma} \begin{lemma} Note that $\boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, 1) = \rho$ and let $\boldsymbol{\mathbf{\tilde{\Gamma}}}(\boldsymbol{\mathbf{\gamma}})$ be a $q \times q$ sparse matrix whose non-zero elements are \begin{align*} \boldsymbol{\mathbf{\tilde{\Gamma}}}(\boldsymbol{\mathbf{\gamma}}) \vcentcolon= \begin{pmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \rho \, \boldsymbol{\mathbf{I}}_{8+|\mathscr{G}|} & \cdot & \cdot \\ \cdot & \cdot & \boldsymbol{\mathbf{\Gamma}}(\boldsymbol{\mathbf{\gamma}}, p) & \cdot \\ \cdot& \cdot & \cdot & \cdot \end{pmatrix}. \end{align*} Moreover, let \begin{align*} \mathscr{U}_{C} &\vcentcolon= \{1, \ldots, r\} \times \{1, \ldots, q\}, \\ \mathscr{U}_{\pi} &\vcentcolon= \{(i,j) : i=j \text{ and } 7 + |\mathscr{G}| < i < r\} \cup \{(i,j) : i=r \text{ and } r \leq j < r+p\}, \end{align*} wherein the latter can be partitioned as $\{\mathscr{C}(i,j), (i,j), \mathscr{C}''(i,j)\}$ for any $(i,j) \in \mathscr{U}_{\pi}$. Hence, 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}}) \, \hat{E}^{(2)}_{i, j, s} - \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{C} \\ (l_1, l_2) \neq (i,j)}} \hat{\Sigma}_{i, l_1, s}^{k^{-1}} (\boldsymbol{\mathbf{\gamma}}) \, \hat{C}_{l_1,l_2, s}^{k+\mathbb{I}_{(l_1,l_2) \in \mathscr{C}(i,j)}}(\boldsymbol{\mathbf{\gamma}})\, \hat{E}^{(3)}_{l_2, j, s}, \; \frac{\alpha}{2} \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}}) \right]}{\hat{\Sigma}_{i,i,s}^{k^{-1}} (\boldsymbol{\mathbf{\gamma}}) \, \hat{E}^{(3)}_{j,j,s} + (1-\alpha) \, \tilde{\Gamma}_{j,j}\,(\boldsymbol{\mathbf{\gamma}})}, \end{align*} for any $(i,j) \in \mathscr{U}_{\pi}$ and constant to the values in (ref) for the remaining entries, and with $\mathcal{S}$ being the soft-thresholding operator. \end{lemma} \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{E}}}^{(1)}_s - \boldsymbol{\mathbf{\hat{E}}}^{(2)}_s \boldsymbol{\mathbf{\hat{C}}}^{{k+1}'}_{s}(\boldsymbol{\mathbf{\gamma}}) - \boldsymbol{\mathbf{\hat{C}}}^{k+1}_{s}(\boldsymbol{\mathbf{\gamma}}) \, \boldsymbol{\mathbf{\hat{E}}}^{(2)'}_s + \boldsymbol{\mathbf{\hat{C}}}^{k+1}_{s} (\boldsymbol{\mathbf{\gamma}}) \, \boldsymbol{\mathbf{\hat{E}}}^{(3)}_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} The CM-step for the free parameters of the measurement equation requires an ad-hoc approach due to the implicit equality constraints for the households described in (ref) -- i.e., all households within a given group have the same factor loadings. While linear constraints have been handled before in EM-like algorithms for time-series models banbura2014maximum, our problem is a bit more complicated. Indeed, it is not advised to estimate an unconstrained version of the model and apply the restrictions ex-post, since each consumer unit is observed for very short periods of time. Hence, we handle this CM-step as the constrained optimisation problem in (ref).\footnote{We do not need to use Lagrangian multipliers, given that the constraints can be implemented by directly plugging them into the expected log-likelihood via $\boldsymbol{\mathbf{B}}$ as described in (ref).} \begin{proposition} Let \begin{alignat*}{2} &\boldsymbol{\mathbf{\ddot{A}}}_{t} &&\vcentcolon= \boldsymbol{\mathbf{A}}_{t}' \boldsymbol{\mathbf{A}}_{t}, \\[0.5em] &\boldsymbol{\mathbf{\hat{F}}}_{t} &&\vcentcolon= \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' + \boldsymbol{\mathbf{\hat{P}}}_{t}, \\[0.5em] &\boldsymbol{\mathbf{\hat{G}}}_{s} &&\vcentcolon= \sum_{t \in \mathscr{T}(s)} \boldsymbol{\mathbf{A}}_{t}' \boldsymbol{\mathbf{Y}}_{t}^{obs} \,\boldsymbol{\mathbf{\hat{\Phi}}}_{t}', \end{alignat*} and \begin{alignat*}{2} &\hat{\CyBe}_{i, j, s} &&\vcentcolon= \sum_{k=\omega_{i}^{\dag}}^{\omega_{i}^{\ddag}} \hat{G}_{k, j, s}, \\ &\CyEl_{i, j, t} &&\vcentcolon= \sum_{k=\omega_{i}^{\dag}}^{\omega_{i}^{\ddag}} \ddot{A}_{k, j, t}, \\ &\CySha_{i, j, t} &&\vcentcolon= \sum_{k=\omega_{i}^{\dag}}^{\omega_{i}^{\ddag}} \sum_{l=\omega_{j}^{\dag}}^{\omega_{j}^{\ddag}} \ddot{A}_{k, l, t}, \end{alignat*} where $\omega_{i}^{\dag} \vcentcolon= 9+\sum_{1 \leq j < i-7} \omega_{j}$ and $\omega_{i}^{\ddag} \vcentcolon= 8+\sum_{1 \leq j \leq i-7} \omega_{j}$, for $i=8, \ldots, 7+|\mathscr{G}|$. Let also \begin{align*} \mathscr{U}_{\Lambda} &\vcentcolon= \{1, \ldots, 7+|\mathscr{G}|\} \times \{0, \ldots, p-1\}. \end{align*} It follows that, when the penalty is not active and at a generic $k+1$ iteration of the ECM algorithm, the factor loadings \begin{align*} \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} &= \frac{1}{\sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t}} \, \Bigg\{ \hat{G}_{i+1,j+r,s} - \sum_{t \in \mathscr{T}(s)} \bigg[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \ddot{A}_{i+1,l_1,t} \\[0.5em] &\quad \qquad + \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1, l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \ddot{A}_{l_1+1, i+1, t} + \mathbb{I}_{l_1 > 7} \, \CyEl_{l_1,i+1,t} \big) \bigg] \Bigg\} \end{align*} for $i=1, \ldots, 7$, and \begin{align*} \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} &= \frac{1}{\sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t}} \, \Bigg\{\hat{\CyBe}_{i, j+r, s} - \sum_{t \in \mathscr{T}(s)} \bigg[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \CyEl_{i,l_1,t} \\ &\quad \qquad + \hat{F}_{r,j+r,t} \, \CyEl_{i,1,t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \CyEl_{i,l_1+1,t} + \mathbb{I}_{l_1 > 7} \, \CySha_{l_1,i,t} \big) \bigg] \Bigg\} \end{align*} for $i=8, \ldots, 7+|\mathscr{G}|$ where $\hat{\Lambda}_{l_{1}, l_{2}+1}^{\diamond}$ is the most up-to-date estimate for $\Lambda_{l_{1}, l_{2}+1}$ available when computing $\hat{\Lambda}^{k+1}_{i, j+1}$, for $j=0, \ldots, p-1$. More generally, it follows that, at a generic $k+1$ iteration of the ECM algorithm, the factor loadings \begin{align*} \hat{\Lambda}^{k+1}_{\,i,j+1}(\boldsymbol{\mathbf{\gamma}}) = \begin{cases} \frac{\mathcal{S} \left[ \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t} \, ,\; \frac{\varepsilon \alpha}{2} \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \right]}{\varepsilon (1-\alpha) \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \;+\; \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t}}, & 1 \leq i \leq 7, \\[1em] \frac{\mathcal{S} \left[ \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t} \, ,\; \frac{\varepsilon \alpha}{2} \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \right]}{\varepsilon (1-\alpha) \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \;+\; \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t}}, & 8 \leq i \leq 7+|\mathscr{G}|, \end{cases} \end{align*} for $j=0, \ldots, p-1$. \end{proposition} \begin{proofenv} We develop the proof in three steps. Step (i) derives the part of the expected log-likelihood that depends on the factor loadings. Step (ii) solves the maximisation problem assuming that the penalty is not active (i.e., $\rho=0$). This leads to a CM-step similar to the M-step usually employed in non-regularised EM-algorithms for dynamic factor models such as banbura2014maximum. Step (iii) builds on that to write down the formula for the final estimator. (i) Note that, for the linearity of the trace, \begin{align*} &\operatorname{Tr} \left\{ \sum_{t \in \mathscr{T}(s)} \left[ \left( \boldsymbol{\mathbf{Y}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right) \left( \boldsymbol{\mathbf{Y}}_{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\} \\ &\qquad = \sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left[ \left( \boldsymbol{\mathbf{Y}}_{t}^{obs} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \right) \left( \boldsymbol{\mathbf{Y}}_{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]. \end{align*} The part of this trace that depends on the measurement coefficients is \begin{align*} &\sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{{obs}'} - \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \, \boldsymbol{\mathbf{Y}}_{t}^{{obs}'} - \boldsymbol{\mathbf{Y}}_{t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' \, \underline{\boldsymbol{\mathbf{B}}}_{\,t}^{{obs}'} \right) \\ &\qquad = \sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' - \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \, \boldsymbol{\mathbf{Y}}_{t}^{{obs}'} - \boldsymbol{\mathbf{Y}}_{t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' \, \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' \right) \\ &\qquad = \sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \boldsymbol{\mathbf{A}}_{t} \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \boldsymbol{\mathbf{A}}_{t}' - \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t} \, \boldsymbol{\mathbf{Y}}_{t}^{{obs}'} \, \boldsymbol{\mathbf{A}}_{t} - \boldsymbol{\mathbf{A}}_{t}' \, \boldsymbol{\mathbf{Y}}_{t}^{obs} \, \boldsymbol{\mathbf{\hat{\Phi}}}_{t}' \, \underline{\boldsymbol{\mathbf{B}}}' \right) \\ &\qquad = \sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) - \operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right) - \operatorname{Tr} \left(\boldsymbol{\mathbf{\hat{G}}}_{s} \, \underline{\boldsymbol{\mathbf{B}}}' \right) \\ &\qquad = \sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) - 2 \operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right). \end{align*} In order to write down the CM-step for the factor loadings, we develop these traces as functions of $\underline{\boldsymbol{\mathbf{\Lambda}}}$.\footnote{Indeed, these are the only free parameters in $\underline{\boldsymbol{\mathbf{B}}}$.} We start with the simpler trace. Under the identification restrictions and constraints in (ref), \begin{align*} &\operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right) \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{q} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=1}^{NK} \sum_{j=r}^{q} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=1}^{8} \sum_{j=r}^{q} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=9}^{NK} \sum_{j=r}^{q} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=1}^{8} \sum_{j=r}^{r+p-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=9}^{NK} \sum_{j=r}^{r+p-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} + \sum_{i=1}^{8} \sum_{j=0}^{p-1} \underline{B}_{\,i,j+r} \, \hat{G}_{i,j+r,s} + \sum_{i=9}^{NK} \sum_{j=0}^{p-1} \underline{B}_{\,i,j+r} \, \hat{G}_{i,j+r,s} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \underline{B}_{\,i,j} \, \hat{G}_{i,j,s} \, + \, \hat{G}_{1,r,s} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \, \hat{G}_{i+1,j+r,s} + \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \hat{\CyBe}_{i, j+r, s}. \end{align*} Thus, \begin{align} \operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right) \propto \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \, \hat{G}_{i+1,j+r,s} + \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \hat{\CyBe}_{i, j+r, s}. \end{align} Next, we focus on the most complicated component of the expected log-likelihood. Since $\boldsymbol{\mathbf{\ddot{A}}}_{t}$ and $\boldsymbol{\mathbf{\hat{F}}}_{t}$ are symmetric, \begin{align*} &\operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{q} \sum_{k=1}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + \sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad \qquad + \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \sum_{k=r}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + \sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad= \sum_{i=1}^{NK} \sum_{j=1}^{r-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + 2\sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad \qquad + \sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t}. \end{align*} Thus, the only part of the latter trace that depends on the factor loadings is \begin{align} \sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + 2\sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t}. \end{align} The first term of (ref) is \begin{align*} &\sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad= \sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + \sum_{i=9}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad \qquad + \sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + \sum_{i=9}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad= \sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} + \sum_{i=9}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad \qquad + 2\sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t}. \end{align*} Under the identification restrictions and constraints in (ref), \begin{align*} &\sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=1}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \hat{F}_{r,r,t} \, \ddot{A}_{1,1,t} + 2 \sum_{i=2}^{8} \sum_{j=r}^{q} \underline{B}_{\,i,j} \, \hat{F}_{j,r,t} \, \ddot{A}_{1,i,t} + \sum_{i=2}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=2}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \hat{F}_{r,r,t} \, \ddot{A}_{1,1,t} + 2 \sum_{i=2}^{8} \sum_{j=r}^{r+p-1} \underline{B}_{\,i,j} \, \hat{F}_{j,r,t} \, \ddot{A}_{1,i,t} + \sum_{i=2}^{8} \sum_{j=r}^{r+p-1} \sum_{k=r}^{r+p-1} \sum_{l=2}^{8} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \hat{F}_{r,r,t} \, \ddot{A}_{1,1,t} + 2 \sum_{i=2}^{8} \sum_{j=0}^{p-1} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i,t} + \sum_{i=2}^{8} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=2}^{8} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,i,t} \\ &\qquad = \hat{F}_{r,r,t} \, \ddot{A}_{1,1,t} + 2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=1}^{7} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \ddot{A}_{l+1, i+1, t}. \end{align*} Also, \begin{align*} &\sum_{i=9}^{NK} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=9}^{NK} \sum_{j=r}^{r+p-1} \sum_{k=r}^{r+p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=9}^{NK} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CySha_{l,i,t}. \end{align*} Moreover, \begin{align*} &\sum_{i=1}^{8} \sum_{j=r}^{q} \sum_{k=r}^{q} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=1}^{8} \sum_{j=r}^{r+p-1} \sum_{k=r}^{r+p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=1}^{8} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{k=0}^{p-1} \sum_{l=9}^{NK} \hat{F}_{r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,1,t} + \sum_{i=2}^{8} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=9}^{NK} \sum_{j=0}^{p-1} \hat{F}_{r,j+r,t} \, \underline{B}_{\,i,j+r} \, \ddot{A}_{i,1,t} + \sum_{i=2}^{8} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=9}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k+r,t} \, \underline{B}_{\,l,k+r} \, \ddot{A}_{l,i,t} \\ &\qquad = \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \hat{F}_{r,j+r,t} \, \underline{\Lambda}_{\,i,j+1} \, \CyEl_{i,1,t} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CyEl_{l,i+1,t}. \end{align*} Thus, the first term of (ref) is proportional to \begin{align} &2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=1}^{7} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \ddot{A}_{l+1, i+1, t} \\ &\qquad + \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CySha_{l,i,t} + 2 \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \hat{F}_{r,j+r,t} \, \underline{\Lambda}_{\,i,j+1} \, \CyEl_{i,1,t} \nonumber \\ &\qquad \qquad + 2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CyEl_{l,i+1,t}. \nonumber \end{align} Finally, the second term of (ref) is \begin{align*} &2\sum_{i=1}^{NK} \sum_{j=r}^{q} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = 2\sum_{i=1}^{NK} \sum_{j=r}^{r+p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j} \, \hat{F}_{j,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = 2\sum_{i=1}^{NK} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{l,i,t} \\ &\qquad = 2\sum_{i=1}^{NK} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t} \\ &\qquad = 2\sum_{i=1}^{8} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t} + 2\sum_{i=9}^{NK} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t}, \end{align*} where \begin{align*} &\sum_{i=1}^{8} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t} \\ &\qquad = \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \hat{F}_{r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{1,l,t} + \sum_{i=2}^{8} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t} \\ &\qquad = \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \hat{F}_{r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{1,l,t} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i+1,l,t} \end{align*} and \begin{align*} &\sum_{i=9}^{NK} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{B}_{\,i,j+r} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i,l,t} = \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \CyEl_{i,l,t}. \end{align*} Hence, the second term of (ref) is proportional to \begin{align} 2\sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i+1,l,t} + 2\sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \CyEl_{i,l,t}. \end{align} Combining \crefrange{eq:trace_bobn_1}{eq:trace_bobn_3}, it follows that \begin{align} &\operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) \\ &\qquad \propto 2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=1}^{7} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \ddot{A}_{l+1, i+1, t} \nonumber \\ &\qquad \qquad + 2 \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \hat{F}_{r,j+r,t} \, \underline{\Lambda}_{\,i,j+1} \, \CyEl_{i,1,t} + \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CySha_{l,i,t} \nonumber \\ &\qquad \qquad \qquad + 2\sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i+1,l,t} + 2\sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \CyEl_{i,l,t} \nonumber \\ &\qquad \qquad \qquad \qquad + 2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CyEl_{l,i+1,t} \nonumber. \end{align} Finally, it follows from (ref) and (ref) that \begin{align} &\sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) - 2 \operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right) \\ &\qquad \propto \sum_{t \in \mathscr{T}(s)} \sum_{i=1}^{7} \sum_{j=0}^{p-1} \bigg( 2 \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{k=0}^{p-1} \sum_{l=1}^{7} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \ddot{A}_{l+1, i+1, t} \nonumber \\ &\qquad \qquad + 2 \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \ddot{A}_{i+1,l,t} \bigg) -2 \sum_{i=1}^{7} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \, \hat{G}_{i+1,j+r,s} \nonumber \\ &\qquad \qquad \qquad + \sum_{t \in \mathscr{T}(s)} \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \bigg( 2 \hat{F}_{r,j+r,t} \, \underline{\Lambda}_{\,i,j+1} \, \CyEl_{i,1,t} + \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CySha_{l,i,t} \nonumber \\ &\qquad \qquad \qquad \qquad + 2 \sum_{k=1}^{r-1} \sum_{l=1}^{NK} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k,t} \, \underline{B}_{\,l,k} \, \CyEl_{i,l,t} \bigg) -2 \sum_{i=8}^{7+|\mathscr{G}|} \sum_{j=0}^{p-1} \underline{\Lambda}_{\,i, j+1} \hat{\CyBe}_{i, j+r, s} \nonumber \\ &\qquad \qquad \qquad \qquad \qquad + 2 \sum_{t \in \mathscr{T}(s)} \sum_{i=1}^{7} \sum_{j=0}^{p-1} \sum_{k=0}^{p-1} \sum_{l=8}^{7+|\mathscr{G}|} \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,k+r,t} \, \underline{\Lambda}_{\,l,k+1} \, \CyEl_{l,i+1,t} \nonumber. \end{align} We have rearranged the terms in (ref) so that the first two rows refer to the factor loadings of the macroeconomic indicators, the third and fourth row refer to the ones of the households and the last row to both of them. (ii) When the penalty is not active, the CM-step is computed from (ref) since \begin{align*} \frac{\partial \mathcal{M}_e \left[\underline{\boldsymbol{\mathbf{\vartheta}}}, \boldsymbol{\mathbf{\gamma}} \,|\, \mathscr{Y}(s), \boldsymbol{\mathbf{\hat{\vartheta}}}_{s}^{k}(\boldsymbol{\mathbf{\gamma}}) \right]}{\partial \underline{\boldsymbol{\mathbf{\Lambda}}}} = -\frac{1}{2 \varepsilon} \frac{\partial \left[\sum_{t \in \mathscr{T}(s)} \operatorname{Tr} \left( \underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{F}}}_{t} \, \underline{\boldsymbol{\mathbf{B}}}' \, \boldsymbol{\mathbf{\ddot{A}}}_{t} \right) - 2 \operatorname{Tr} \left(\underline{\boldsymbol{\mathbf{B}}} \, \boldsymbol{\mathbf{\hat{G}}}_{s}' \right)\right]}{\partial \underline{\boldsymbol{\mathbf{\Lambda}}}}. \end{align*} We structure the CM-step by following analogous steps to those in pellegrino2021factoraugmented. Indeed, we estimate $\boldsymbol{\mathbf{\Lambda}}$ one entry at the time, starting from the $\Lambda_{1, 1}$ and in column-major order. In other words, the derivative of (ref) with respect to $\Lambda_{i,j+1}$ is taken having fixed the other factors loadings to their latest estimate, for any $i=1, \ldots, 7+|\mathscr{G}|$ and $j=0, \ldots, p-1$. Formally, at a generic $k+1$ iteration of the ECM algorithm, this derivative is equal to \begin{align*} &\frac{\hat{G}_{i+1,j+r,s}}{\varepsilon} - \frac{1}{\varepsilon} \sum_{t \in \mathscr{T}(s)} \Big(\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \ddot{A}_{i+1,l_1,t} + \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \underline{\Lambda}_{\,i,j+1} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t} \\ &\qquad + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j) \\ l_1 \leq 7}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1, l_2+1} \, \ddot{A}_{l_1+1, i+1, t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j) \\ l_1 > 7}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \CyEl_{l_1,i+1,t} \Big) \end{align*} when computed with respect to any factor loading associated to the macroeconomic aggregates, and \begin{align*} &\frac{\hat{\CyBe}_{i, j+r, s}}{\varepsilon} - \frac{1}{\varepsilon} \sum_{t \in \mathscr{T}(s)} \Big(\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \CyEl_{i,l_1,t} + \hat{F}_{r,j+r,t} \, \CyEl_{i,1,t} + \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t} \\ &\qquad + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j) \\ l_1 \leq 7}} \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \hat{F}_{l_2+r,j+r,t} \, \CyEl_{i,l_1+1,t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j) \\ l_1 > 7}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \CySha_{l_1,i,t} \Big) \end{align*} when computed with respect to any factor loading associated to the households data. These derivatives can be equivalently written in the compact forms \begin{align*} &\frac{\hat{G}_{i+1,j+r,s}}{\varepsilon} - \frac{1}{\varepsilon} \sum_{t \in \mathscr{T}(s)} \Big[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \ddot{A}_{i+1,l_1,t} + \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \underline{\Lambda}_{\,i,j+1} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t} \\ &\qquad + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1, l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \ddot{A}_{l_1+1, i+1, t} + \mathbb{I}_{l_1 > 7} \, \CyEl_{l_1,i+1,t} \big) \Big] \nonumber \end{align*} and \begin{align*} &\frac{\hat{\CyBe}_{i, j+r, s}}{\varepsilon} - \frac{1}{\varepsilon} \sum_{t \in \mathscr{T}(s)} \Big[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \CyEl_{i,l_1,t} + \hat{F}_{r,j+r,t} \, \CyEl_{i,1,t} + \underline{\Lambda}_{\,i,j+1} \, \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t} \\ &\qquad + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \CyEl_{i,l_1+1,t} + \mathbb{I}_{l_1 > 7} \, \CySha_{l_1,i,t} \big) \Big] \nonumber \end{align*} respectively. It follows that, when the penalty is not active and at a generic $k+1$ iteration of the ECM algorithm, \begin{align*} \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} &= \frac{1}{\sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t}} \, \Bigg\{ \hat{G}_{i+1,j+r,s} - \sum_{t \in \mathscr{T}(s)} \bigg[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \ddot{A}_{i+1,l_1,t} \\[0.5em] &\quad \qquad + \hat{F}_{j+r,r,t} \, \ddot{A}_{1,i+1,t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1, l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \ddot{A}_{l_1+1, i+1, t} + \mathbb{I}_{l_1 > 7} \, \CyEl_{l_1,i+1,t} \big) \bigg] \Bigg\} \end{align*} for $i=1, \ldots, 7$, and \begin{align*} \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} &= \frac{1}{\sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t}} \, \Bigg\{\hat{\CyBe}_{i, j+r, s} - \sum_{t \in \mathscr{T}(s)} \bigg[\sum_{l_1=1}^{NK} \sum_{l_2=1}^{r-1} \hat{F}_{j+r,l_2,t} \, B_{\,l_1,l_2} \, \CyEl_{i,l_1,t} \\ &\quad \qquad + \hat{F}_{r,j+r,t} \, \CyEl_{i,1,t} + \sum_{\substack{(l_1, l_2) \in \mathscr{U}_{\Lambda} \\ (l_1, l_2) \neq (i, j)}} \hat{F}_{j+r,l_2+r,t} \, \hat{\Lambda}^{\diamond}_{\,l_1,l_2+1} \, \big( \mathbb{I}_{l_1 \leq 7} \, \CyEl_{i,l_1+1,t} + \mathbb{I}_{l_1 > 7} \, \CySha_{l_1,i,t} \big) \bigg] \Bigg\} \end{align*} for $i=8, \ldots, 7+|\mathscr{G}|$. (iii) It follows directly from the results in step (i) and step (ii), and the proof of pellegrino2021factoraugmented that, at a generic $k+1$ iteration of the ECM algorithm, \begin{align*} \hat{\Lambda}^{k+1}_{\,i,j+1}(\boldsymbol{\mathbf{\gamma}}) = \begin{cases} \frac{\mathcal{S} \left[ \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t} \, ,\; \frac{\varepsilon \alpha}{2} \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \right]}{\varepsilon (1-\alpha) \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \;+\; \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \ddot{A}_{i+1, i+1, t}}, & \text{if } 1 \leq i \leq 7, \\[1em] \frac{\mathcal{S} \left[ \hat{\Lambda}_{\,i,j+1}^{QMLE, \, k+1} \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t} \, ,\; \frac{\varepsilon \alpha}{2} \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \right]}{\varepsilon (1-\alpha) \Gamma_{j+1, j+1} (\boldsymbol{\mathbf{\gamma}}, p) \;+\; \sum_{t \in \mathscr{T}(s)} \hat{F}_{j+r,j+r,t} \, \CySha_{i,i,t}}, & \text{if } 8 \leq i \leq 7+|\mathscr{G}|, \end{cases} \end{align*} for $j=0, \ldots, p-1$. \end{proofenv} \subsection{Initialisation} First, we compute group averages for the microeconomic data. Then we apply the procedure described in pellegrino2021factoraugmented on both the macroeconomic indices and group averages. For simplicity, we set $\lambda = 2.573$, $\alpha=0.667$ and $\beta=1.326$. These are the optimal values in pellegrino2021factoraugmented converted for quarterly frequency data. \subsection{Enforcing causality during the estimation} We enforce causality during the estimation following pellegrino2021factoraugmented. \section{Additional charts} \begin{figure}[!h] \caption{Core driver of RIPH computed as the sum of trend and business cycle. \\\textbf{Notes}: The model is estimated with quarterly data from October 1989 to December 2019.} \end{figure}