The Accuracy Smoothness Dilemma in Prediction: a Novel Multivariate M-SSA Forecast Approach

The Accuracy Smoothness Dilemma in Prediction: a Novel Multivariate M-SSA Forecast Approach

\\ \\ Keywords: Forecasting, signal extraction, smoothing, zero-crossings, sign accuracy, mean-squared error, business-cycles.

Introduction

Forecasting presents a complex estimation challenge, as it necessitates the consideration of various, often conflicting, priorities and objectives. In this, we focus on accuracy and smoothness (AS). Accuracy pertains to the estimation of the future level of a time series while smoothness serves to regulate `noisy' changes. Ideally, these two dimensions could be jointly optimized, resulting in a predictor that closely approximates the true, albeit unobserved, future observation, thereby minimizing the incidence of spurious changes. However, the relationship between accuracy and smoothness presents a predictive dilemma; when appropriately formalized, an enhancement in one dimension inevitably leads to a compromise in the other. In this study, we propose a novel multivariate framework that enables users to manage and balance both accuracy and smoothness effectively. Furthermore, we extend existing methodologies to incorporate the proposed tradeoff, allowing for a customization of traditional `benchmark' predictors in terms of AS performances.\\

Smoothness of a predictor can be characterized by various concepts, such as its curvature, which is a geometric measure based on the squared second-order differences of the predictor. In contrast, the SSA proposed in Wildi (2024,2026) introduces a novel concept of smoothness that focuses on the expected duration between sign changes or zero crossings of a (stationary, zero-mean) predictor. The analysis of zero-crossings in time series data was first introduced by Rice (1944), who established a correlation between the autocorrelation function (ACF) of a zero-mean stationary Gaussian process and the expected number of zero crossings within a specified time interval. In the context of non-stationary time series, sign changes in the growth rate signify transitions between phases of expansion and contraction. These alternating periods of growth and contraction can be linked to the business cycle, provided that the fluctuations exhibit sufficient magnitude and duration. To be classified as representative of a business cycle, consecutive zero-crossings of the growth rate must be separated by intervals potentially spanning several years. This requirement for duration necessitates a smooth trajectory of the corresponding economic indicator. The relationship between the smoothness of a time series and the mean duration between consecutive sign changes—termed the holding time (HT) by Wildi (2024)—is formalized in the work of Rice (1944). The new M-SSA integrates accuracy within its objective function and smoothness through its holding time constraint, resulting in the M-SSA predictor being derived as the outcome of a constrained optimization process. \\

The traditional mean-squared error (MSE) forecast paradigm predominantly emphasizes accuracy, often at the expense of smoothness. In certain applications, this focus on accuracy can lead to excessive noise leakage, thereby increasing the likelihood of false alarms, as demonstrated in Wildi (2024). To address this issue, McElroy and Wildi (2019) propose a forecasting approach based on a forecast trilemma; however, their solution does not take zero-crossings into account. Wildi (2024,2026) introduces the SSA optimization criterion, which explicitly regulates the HT of a predictor through a smoothing constraint derived from Rice's foundational work. We posit that regulating the rate of sign changes (in the growth rate), while simultaneously maintaining optimal tracking accuracy as measured by minimal MSE performance, provides a viable alternative to traditional filtering and smoothing methods, aligning with the decision-making and control logic pertinent during phases of alternating growth.Accordingly, we extend this approach to multivariate prediction problems, referred to as M-SSA, with a focus on a stationary framework to clarify the core concepts. An extension of the univariate SSA to non-stationary integrated processes has been explored in Wildi (2026), where non-stationary maximal monotone and minimal curvature predictors are derived. These findings can be further generalized within the M-SSA framework. \\

Our applications aim to highlight some of the distinctive features of M-SSA in forecasting, real-time signal extraction, and smoothing. All examples can be replicated using an open-source M-SSA package, which includes an R package along with comprehensive instructions, practical use cases, and theoretical results, available at https://github.com/wiaidp/R-package-SSA-Predictor.git. A notable feature of our approach is the ability to customize benchmark predictors in terms of accuracy and smoothness (AS) performance. The M-SSA package demonstrates applications of this customization to various filtering techniques, including the Hodrick-Prescott (HP) filter (Hodrick and Prescott, 1997), Hamilton's regression filter (Hamilton, 2018), the Baxter-King (BK) bandpass filter (Baxter and King, 1999) and a `refined' Beveridge-Nelson design (Kamber, Morley and Wong (2024)).\\

Section (ref) briefly reviews the (univariate) SSA approach and introduces to the M-SSA criterion; Section (ref) proposes solutions to the optimization problem; Section (ref) illustrates applications in forecasting, signal extraction and smoothing; finally, Section (ref) summarizes our main findings.

M-SSA Criterion

We provide a brief review of the univariate case and the various adaptations of the Smooth Sign Accuracy (SSA) optimization criterion introduced by Wildi (2024). Following this review, we derive the multivariate generalization M-SSA.

Univariate SSA

We consider a zero-mean stationary time series denoted as $x_t$ (processes with a non-vanishing mean are assumed to be centered) and a target variable $z_{t+\delta}$, where $ \delta \in \mathbb{Z}$, which depends on future values $x_{t-k}$ for $k < 0$. We then derive an optimal predictor $y_t$ for $z_{t+\delta}$ based on past values $x_{t-k}$ for $k \geq 0$ (a causal filter). This predictor is characterized by the property that the HT of $y_t$—defined as the mean duration between consecutive zero-crossings—can be specified by the user. For clarity and simplicity, we initially assume $x_t = \epsilon_t$ to be an independent and identically distributed white noise (WN) sequence. An extension of this approach to accommodate autocorrelated stationary processes will then be discussed.\\

Let

where $x_j=\epsilon_j, j \in \mathbb{Z}$, is WN. For the sake of clarity we may assume that $\epsilon_t$ is standardized. The sequence $\boldsymbol{\gamma} = (\gamma_k)$, where $k \in \mathbb{Z}$, is a real, square summable sequence, ensuring that $z_t$ constitutes a stationary process with zero mean and variance given by $\sum_{k=-\infty}^{\infty} \gamma_k^2$. We seek to construct a predictor defined as $y_t = \sum_{k=0}^{L-1} b_k \epsilon_{t-k}$ for the target $z_{t+\delta}$, where $\delta \in \mathbb{Z}$ and $b_k$ represents the coefficients of a one-sided causal filter of length $L$. This predictive framework is commonly categorized as fore-casting, now-casting, or back-casting, depending on whether $\delta > 0$, $\delta = 0$, or $ \delta < 0$, respectively. For illustration, consider the specific case where $\gamma_0 = 1$, $\gamma_1 = 0.5$, and $\gamma_k = 0$ for $k \notin {0, 1}$. In this scenario, $z_t = \epsilon_t + 0.5 \epsilon_{t-1}$ represents a moving-average process of order one. The task of one-step ahead forecasting for $z_t$ can be achieved by setting $\delta = 1$. The classic mean squared error (MSE) estimate for $z_{t+1}$ yields $y_t = 0.5 \epsilon_t$, resulting in the coefficients $b_0 = 0.5$ and $b_k = 0$ for $1 \leq k \leq L-1$\footnote{Furthermore, our target specification can also facilitate signal extraction problems, wherein the weights $\gamma_k$ are coefficients of a two-sided (potentially bi-infinite) filter.}.\\

The MSE predictor captures the level of future observations (or targets) optimally; however, it may produce forecasts that are 'noisy'. This phenomenon is evident in the previous moving average (MA(1)) example, where $y_t = 0.5\epsilon_t$ is white noise while the target $z_{t+1}$ is more persistent. Similarly, in real-time signal extraction, the classic one-sided (nowcast) concurrent filter tends to introduce markedly more `noise' relative to the acausal, two-sided target, as demonstrated in our subsequent examples. To mitigate this issue, Wildi (2024) proposes the following optimization problem:

Criterion (ref) is referred to as the SSA criterion, and its solution is denoted as SSA($\rho_1,\delta$). The constraints $\mathbf{b}'\mathbf{Mb}=l\rho_1$ and $\mathbf{b}'\mathbf{b}=l$ are termed the HT constraint and the length constraint, respectively. Smoothness of the predictor is obtained by specification of the hyperparameter $\rho_1$, which controls the first- order autocorrelation (ACF(1)) of the predictor\footnote{A formal link between ACF(1) and HT is established below.}. Here, $\mathbf{b}=(b_0,...,b_{L-1})'$ and $\boldsymbol{\gamma}_{\delta}=(\gamma_{\delta},...,\gamma_{\delta+L-1})'$ are column vectors of dimension $L$. The variable $l$ represents a scaling parameter, the specific choice of which will be discussed subsequently. This parameter is omitted from the criterion notation SSA($\rho_1,\delta$) due to its relatively minor importance in this context, see below for further reference. The autocovariance generating matrix $\mathbf{M}$ is defined as: \[ \mathbf{M}=\left(

\right), \] which has dimension $L\times L$. The matrix $\mathbf{M}$ is structured such that $\mathbf{b}'\mathbf{Mb}=\sum_{k=1}^{L-1}b_{k-1}b_k$, representing the first-order autocovariance of $y_t$ under the specified assumptions. This formulation aims to optimize the predictor while controlling for the noise present in the forecasts. Under the assumption of WN, the classic MSE predictor is defined as $y_{t,MSE}:=\boldsymbol{\gamma}_{\delta}'\mathbf{x}_t$, where $\mathbf{x}_t:=(x_t,...,x_{t-(L-1)})'$. The weights $\boldsymbol{\gamma}_{\delta}$ can be derived as a solution to the SSA criterion by setting $l:=\boldsymbol{\gamma}_{\delta}'\boldsymbol{\gamma}_{\delta}$ and $\rho_1:=\boldsymbol{\gamma}_{\delta}'\mathbf{M}\boldsymbol{\gamma}_{\delta}/\boldsymbol{\gamma}_{\delta}'\boldsymbol{\gamma}_{\delta}=:\rho_{MSE}$. However, the objective is to ensure that the SSA predictor $y_t := \mathbf{b}' \boldsymbol{\epsilon}_t$ exhibits reduced noise compared to $ y_{t,MSE}$ for which we can use the hyperparameter $\rho_1$ in the SSA criterion (ref). \\

SSA Interpretations

To streamline the terminology, we will refer to both $y_{t,MSE}$ and $\boldsymbol{\gamma}_{\delta}$ as the MSE predictor. Similarly, we will merge $y_t$ and $\mathbf{b}$ under the (M-)SSA designation, clarifying our intent when necessary. Under the assumption of white noise, $\boldsymbol{\gamma}_{\delta}$ represents the appropriate target $\boldsymbol{\gamma}$ in Criterion (ref), given that $\gamma_k$ is irrelevant for $k < \delta$ or $k > \delta + L - 1$. We will also assume $\boldsymbol{\gamma}_{\delta} \neq \mathbf{0}$. In this framework, the solution $\mathbf{b}_0$ to the SSA criterion can be interpreted as a constrained predictor for $z_{t+\delta}$. Alternatively, $\mathbf{b}_0$ can be viewed as a `smoother' for $y_{t,MSE}$, as discussed by Wildi (2024,2026). Thus, the SSA criterion simultaneously addresses prediction and smoothing. Moreover, under the imposed length constraint, the objective function $\mathbf{b}' \boldsymbol{\gamma}_{\delta}$ is proportional to: \[ \rho(y,z,\delta):=\mathbf{b}'\boldsymbol{\gamma}_{\delta}/\sqrt{l\boldsymbol{\gamma}'\boldsymbol{\gamma}} \] which represents the target correlation of $y_t$ with $z_{t+\delta}$, or to \[ \mathbf{b}'\boldsymbol{\gamma}_{\delta}/\sqrt{l\boldsymbol{\gamma}_{\delta}'\boldsymbol{\gamma}_{\delta}} \] the correlation of $y_t$ with $y_{t,MSE}$. Maximizing either of these functions consequently maximizes the other, allowing Criterion (ref) to be reformulated as:

Here, $\rho(y):=\mathbf{b}'\mathbf{Mb}/l=\mathbf{b}'\mathbf{Mb}/\mathbf{b}'\mathbf{b}$ denotes the ACF(1) of $y_t$. An increase in $\rho_1$ results in a stronger first-order ACF for the predictor, thus yielding a `smoother' trajectory for $y_t$, characterized by fewer zero-crossings. \\

Specifically, let $ht(y|\mathbf{b},i)$ be defined as $E[t_i-t_{i-1}]$, where $t_i$, $i\geq 1$ represents consecutive zero-crossings of the process $y_t$\footnote{These crossings satisfy $t_{i-1}<t_i$ with the condition $t_1\geq L$ and $y_{t_{i}-1}y_{t_{i}}<0$ for all $i$. Additionally, it holds that $y_{t-1}y_{t}>0$ if $t_{i-1}<t<t_i$.}. Under the assumptions of stationarity, $ht(y|\mathbf{b},i)=ht(y|\mathbf{b})$ is invariant with respect to $i$. Furthermore, assuming Gaussianity of $x_t$, we have:

establishing a connection between the HT and the first-order ACF, see Wildi (2024). It is noteworthy that this relationship remains relatively robust even when deviations from Gaussian assumptions occur, as discussed in Wildi (2024). Given that the trigonometric function arccos is monotonic in the interval $]-1,1[$, we can reformulate the SSA criterion as follows:

Since correlations, signs, and zero-crossings are invariant to the scaling of $y_t$, the parameter $l$ in the length constraint can be treated as a nuisance parameter, which primarily serves to ensure uniqueness. We typically assume $l = 1$ unless otherwise stated. If necessary, `static' adjustments for level and scale can be applied post hoc, once a solution $ \mathbf{b} = \mathbf{b}(l)$ has been determined for an arbitrary $l$. This adjustment can be performed, for instance, by regressing the predictor on the target variable.\\

However, our primary focus lies in the `dynamic' aspects of the prediction problem, particularly concerning the target correlation and sign accuracy (SA). Specifically, let $SA(y_t) = P(z_{t+\delta}y_t>0)$ denote the probability of the target and predictor sharing the same sign. Under the Gaussian assumption, it follows that: \[ SA(y_t)=0.5+\frac{\arcsin(\rho(y,z,\delta))}{\pi} \] thereby establishing a strictly monotonic (bijective) relationship between SA and target correlation. Consequently, the SSA criterion can be reformulated as:

In conclusion, the SSA framework effectively reconciles MSE, SA and HT, i.e. smoothing, in a flexible and interpretable manner.

Multivariate M-SSA

Let $\mathbf{x}_t = (x_{1t}, \ldots, x_{nt})'$ represent a multivariate stationary process of dimension $n$, characterized by a (purely non-deterministic) Wold decomposition of the form

where $\boldsymbol{\Xi}_0 = \mathbf{I}_{n \times n}$ denotes the identity matrix, and $\boldsymbol{\epsilon}_t = (\epsilon_{1t}, \ldots, \epsilon_{nt})'$ is a sequence of multivariate WN with a variance-covariance matrix denoted as $\boldsymbol{\Sigma}$. The entries of $\boldsymbol{\Sigma}$ are represented by $\sigma_{ij}$, and we assume that its eigenvalues $\tilde{\sigma}_j$ for $j = 1, \ldots, n$ are strictly positive, indicating full rank\footnote{If $\boldsymbol{\Sigma}$ is originally rank deficient, we assume that redundant noise components have been eliminated, resulting in a lower-dimensional system that retains full rank.}. The eigenvectors of $\boldsymbol{\Sigma}$ are denoted by $\mathbf{v}_{\sigma k}$ for $k = 1, \ldots, n$. Next, we define a multivariate stationary process $\mathbf{z}_t$ of dimension $n$ as follows: \[ \mathbf{z}_t=\sum_{|k|<\infty}\boldsymbol{\Gamma}_k\mathbf{x}_{t-k}, \] where $\boldsymbol{\Gamma}_k$ is a sequence of square-summable $n\times n$ dimensional filter matrices, with entries $(\gamma_{ijk})$ for $i, j = 1, \ldots, n$. The condition for square summability is expressed mathematically as: \[ \sum_{|k| < \infty} \text{tr}(\boldsymbol{\Gamma}_k' \boldsymbol{\Gamma}_k) < \infty, \] where $\text{tr}(\cdot)$ denotes the trace operator. We focus on the estimation of $\mathbf{z}_{t+\delta}$, utilizing the $n-$dimensional predictor defined by \[ \mathbf{y}_t = \sum_{k=0}^{L-1} \mathbf{B}_k \mathbf{x}_{t-k}, \] with $n\times n$-dimensional filter matrices $\mathbf{B}_k$, $k=0,...,L-1$. To facilitate the theoretical analysis, we simplify our model by assuming $\mathbf{x}_t = \boldsymbol{\epsilon}_t$, with further extensions to autocorrelated processes addressed in the subsequent section. Let $b_{ijk}$, $1\leq i,j\leq n$, denote the entries of the matrix $\mathbf{B}_k$ for $k = 0, \ldots, L-1$. We define several vector notations as follows:

Here, $\boldsymbol{\gamma}_{\delta}$ and $\mathbf{b}$ are matrices of dimension $(L\cdot n)\times n$. The $i$-th column, $i=1,...,n$, contains the filter weights associated with the $i$-th target (or $i$-th predictor), and the filter weights corresponding to the $j=1,...,n$ series are stacked into this column vector, of total length $n\cdot L$. The $i$-th column vector of $\boldsymbol{\gamma}_{\delta}$ (or $\mathbf{b}$) can be applied to the $L\cdot n$-dimensional $\boldsymbol{\epsilon}_{\cdot t}$, which stacks all $n$ series into a single long data vector. Specifically, we define \[ y_{ijt} := \mathbf{b}_{ij}'\boldsymbol{\epsilon}_{jt}, \] which allows us to express the $i$-th predictor as: \[ y_{it} = \mathbf{b}_i'\boldsymbol{\epsilon}_{\cdot t} = \sum_{j=1}^n y_{ijt}. \] Similarly, the multivariate predictor vector can be written as: \[ \mathbf{y}_t = \mathbf{b}'\boldsymbol{\epsilon}_{\cdot t} \] and the same applies to the target variables. Through orthogonal projection, the MSE predictor is derived as $\hat{\mathbf{y}}_{t,MSE} = \boldsymbol{\gamma}_{\delta}'\boldsymbol{\epsilon}_{\cdot t}$, with components $\hat{y}_{it,MSE}$, $i=1,...,n$. Consider the Kronecker products defined as $\tilde{\mathbf{I}}:=\boldsymbol{\Sigma}\otimes\mathbf{I}_{L\times L}$ and $\tilde{\mathbf{M}}:=\boldsymbol{\Sigma}\otimes\mathbf{M}$, of the $n\times n$ dimensional matrix $\boldsymbol{\Sigma}$, the $L\times L$ dimensional identity $\mathbf{I}_{L\times L}$ and the autocovariance generating matrix $\mathbf{M}$. We can express the following expectations:

The variance of $z_{it}$ is given by \[ E[z_{it}^2]=\sum_{|k|<\infty}\boldsymbol{\gamma}_{i k}'\boldsymbol{\Sigma}\boldsymbol{\gamma}_{i k}, \] where $\boldsymbol{\gamma}_{i k}:=(\gamma_{i1k},...,\gamma_{ink})'$. Consequently, we propose the following multivariate M-SSA criterion for $i=1,...,n$:

for $i=1,...,n$. Under the imposed length constraint, the objective function is proportional to the correlation between the predictor and the $i$-th target (or the $i$-th MSE predictor). Additionally, the HT constraint regulates the first-order ACF or, equivalently, the HT of the predictor. The relationships established between sign accuracy and target correlation, as well as between first-order ACF and HT, as discussed in the previous section, remain valid.

Dependence

We now relax the WN hypothesis and assume a stationary process given by Equation (ref). Then, the target and predictor can be expressed formally as:

Here, $(\boldsymbol{\Gamma}\cdot\boldsymbol{\Xi})_k =\sum_{m\leq k}\boldsymbol{\Gamma}_m\boldsymbol{\Xi}_{k-m}$ and $(\mathbf{B}\cdot\boldsymbol{\Xi})_j=\sum_{m=0}^{\min(L-1,j)}\mathbf{B}_m \boldsymbol{\Xi}_{j-m} $ represent the convolutions of the sequences $\boldsymbol{\Gamma}_k$ and $\mathbf{B}_j$ with the Wold-decomposition $\boldsymbol{\Xi}_m$ of $\mathbf{x}_t$. For given $(\mathbf{B}\cdot\boldsymbol{\Xi})_j$, $j=0,...,L-1$, the original coefficients $\mathbf{B}_{k}$ can be derived through deconvolution:

In our parametrization, $\boldsymbol{\Xi}_0^{-1} = \mathbf{I}$ simplifies this expression. Denote the elements of the matrices as follows: \[ (\boldsymbol{\Gamma}\cdot\boldsymbol{\Xi})_{lmk} \quad \text{and} \quad (\mathbf{B}\cdot\boldsymbol{\Xi})_{lmj}. \] We now assume $\delta\geq 0$\footnote{Similar but slightly more complex expressions would arise for $\delta < 0$ (backcasting), which are omitted here.} and $L$ to be sufficiently large to render $(\boldsymbol{\Gamma} \cdot \boldsymbol{\Xi})_{\delta+k} \approx 0$ and $(\mathbf{B} \cdot \boldsymbol{\Xi})_j \approx 0$ for $k,j \geq L$, so that their contributions to the variances of $\hat{\mathbf{y}}_{t,MSE}$ and $\mathbf{y}_t$ are negligible. We define:

It is important to note that if $\boldsymbol{\Xi}_k=\mathbf{0},~k>0$, resulting in $\mathbf{x}_t=\boldsymbol{\epsilon}_t$, then we have $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}=\boldsymbol{\gamma}_{i\cdot\delta}$ and $(\mathbf{b}\cdot\boldsymbol{\xi})_{i}=\mathbf{b}_i$, as defined in the preceding section. The generalized M-SSA criterion can then be expressed as:

for $i=1,...,n$. This formulation can be solved for $(\mathbf{b}\cdot\boldsymbol{\xi})_i$, as detailed in Theorem (ref) in the next section. The M-SSA coefficients $\mathbf{b}_i$ can subsequently be determined through deconvolution, as indicated in Equation (ref). Wildi (2026) posits an extension of SSA to non-stationary integrated processes, wherein the HT constraint pertains to stationary differences of the predictor. While a corresponding extension to M-SSA is feasible based on findings by McElroy and Wildi (2020), the present discussion is confined to stationary processes, omitting a lengthier treatment of rank-defficient designs (cointegration) in the non-stationary case.\\

Remark: The target filter $\boldsymbol{\Gamma}_k$, for $|k|<\infty$, and the MA-inversion $\boldsymbol{\Xi}_j$, for $j\geq 0$, are combined into the convolution $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}$, for $i=1,...,n$, in the M-SSA criterion (ref). Due to the commutative property of convolution, this combination is invariant to the ordering of the components. Consequently, M-SSA treats the target filter and the Wold decomposition as essentially equivalent `filters', with the distinction being largely semantic. Our notation emphasizes this difference to distinguish between forecasting, smoothing, and signal extraction, which pertain to specific target formulations $\boldsymbol{\Gamma}_k$ for a given process $\boldsymbol{\Xi}_j$. For further clarification, see Section (ref).

Embedding, Smoothing and Efficient Frontier

In the M-SSA criterion (ref), the terms $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}$, $(\mathbf{b}\cdot\boldsymbol{\xi})_i$, $(\mathbf{b}\cdot\boldsymbol{\xi})_i$, $\mathbf{\tilde{M}}$ and $\mathbf{\tilde{I}}$ depend on the data-generating process (DGP), characterized by $\boldsymbol{\Sigma}$ and $\boldsymbol{\Xi}_k$, and on the forecast target $\mathbf{z}_t$, specified by $\boldsymbol{\Gamma}_k$. Accordingly, for $i=1,...n$, the optimal M-SSA coefficient vector $\mathbf{b}_i$ is a function of the DGP, the target specification, and the constraint parameter $\rho_i$. When $\rho_i=\rho_{i,MSE}$, where $\rho_{i,MSE}$ denotes the first-order ACF of the MSE predictor $\hat{y}_{it,MSE}$, the solution reduces to $\mathbf{b}_i=\boldsymbol{\gamma}_{i\cdot\delta}$, thereby recovering the MSE predictor and rendering the HT constraint redundant (detailed in the subsequent section). Hence, the MSE predictor arises as a special case within the M-SSA framework (embedding).\\

The target $z_{it}$ may be causal or acausal. With a causal target, M-SSA is interpreted as a smoother; with an acausal target, it is interpreted as a predictor. Standard multivariate one-step or $h$-step-ahead forecasts are obtained by setting \[ \boldsymbol{\Gamma}_k=\left\{

\right. \] (see section (ref)). A nowcast ($h = 0$) or backcast ($h < 0$) is produced accordingly: the nowcast $y_{it}$ is smoother or noisier than $x_{it}$ depending on whether $\rho_i$ exceeds or falls short of the first-order ACF $\rho(x_i)$ of $x_{it}$ (see section (ref)). More general prediction and smoothing tasks can be accommodated by allowing an arbitrary (generally square summable) sequence $\boldsymbol{\Gamma}_k$ in the target specification $\mathbf{z}_t$ (see section (ref) for real-time trend extraction). In applications, the target specification can encode various research priorities by selecting the signals of interest. Finally, any acausal prediction problem can be expressed as a causal smoothing problem by replacing the target $z_{it}$ in the objective with the causal MSE predictor $\hat{y}_{it,MSE}$, as in (ref) and (ref). This reformulation is referred to as the equivalent causal M-SSA representation. \\

Under its primal specification (ref), or its extension (ref), M-SSA selects the predictor by maximizing correlation with the target, subject to a novel HT smoothing constraint. Conversely, the dual formulation detailed in the subsequent section characterizes the M-SSA predictor as the smoothest solution, with the largest HT, subject to a prescribed target correlation. Thus, for variable admissible $\rho_i$ in the HT constraint\footnote{It is assumed that the constraints under both primal and dual representations are feasible; technical details are provided in the subsequent section.}, M-SSA traces the efficient frontier linking accuracy (target correlation) and smoothness (zero-crossing rate) under the specified criterion.

Solution

For clarity and simplicity, we assume that the process $\mathbf{x}_t = \boldsymbol{\epsilon}_t$ is WN, acknowledging that autocorrelated processes may be accommodated through the extensions proposed in the preceding section. Consequently, M-SSA is adopted throughout in baseline (MSE) extension mode.\\

We consider the orthonormal Fourier eigenvectors defined as $\mathbf{v}_j:=\left(\sin(k\omega_j)/\sqrt{\sum_{k=1}^Lsin(k\omega_j)^2}\right)_{k=1,...,L}$ of the matrix $\mathbf{M}$, associated with the eigenvalues $\lambda_{j}=\cos(\omega_j)$ calculated at discrete Fourier frequencies $\omega_j=j\pi /(L+1)$, $j=1,...,L$, as detailed by Anderson (1975). We denote the solutions to the equation $\frac{\partial \rho(y_{i})}{\partial \mathbf{b}_i} = \mathbf{0}$ as stationary points of the first-order ACF $\rho(y_{i})$.

Proof: For the sake of clarity, we assume that the vector $\mathbf{b}_i$ is constrained to the unit-ellipse defined by $\mathbf{b}_{i}'\tilde{\mathbf{I}}\mathbf{b}_{i}=1$, which implies that $\rho(y_{i})=\mathbf{b}_i'\mathbf{\tilde{M}}\mathbf{b}_i$. To find a stationary point of $\rho(y_{i})$, we set the derivative of the Lagrangian $\mathfrak{L}=\mathbf{b}_i'\mathbf{\tilde{M}}\mathbf{b}_i-\lambda(\mathbf{b}_{i}'\tilde{\mathbf{I}}\mathbf{b}_{i}-1)$ to zero, leading to the equation $\left(\mathbf{\tilde{M}}-\lambda\mathbf{\tilde{I}}\right)\mathbf{b}_i=\mathbf{0}$. Given that $\boldsymbol{\Sigma}$ is assumed to be of full rank, we can manipulate the previous expression by multiplying it with $\boldsymbol{\Sigma}^{-1} \otimes \mathbf{I}_{L \times L} = \mathbf{\tilde{I}}^{-1}$, yielding \[ \left(\mathbf{I}_{n \times n} \otimes \mathbf{M} - \lambda \mathbf{I}_{nL \times nL}\right)\mathbf{b}_i = \mathbf{0}. \] This expression can be rewritten as $n$ distinct sub-space equations $\mathbf{M}\mathbf{b}_{ij} - \lambda \mathbf{b}_{ij} = \mathbf{0}$, for $j = 1, ..., n$. Therefore, $\lambda = \lambda_k$ must be an eigenvalue of $\mathbf{M}$ for some $k \in {1, ..., L}$. Furthermore, since the eigenvalues of $\mathbf{M}$ are distinct, the corresponding eigenvector $\mathbf{v}_k$ is uniquely determined, leading to the conclusion that $\mathbf{b}_{ij}\propto\mathbf{v}_k$, $j=1,...,n$. Additionally, we have \[ \rho(y_{ij})=\frac{\mathbf{b}_{ij}'\mathbf{Mb}_{ij}}{\mathbf{b}_{ij}'\mathbf{b}_{ij}}=\lambda\frac{\mathbf{b}_{ij}'\mathbf{b}_{ij}}{\mathbf{b}_{ij}'\mathbf{b}_{ij}}=\lambda_k. \] Since $\mathbf{b}_{ij}\propto\mathbf{v}_{k}$, we can express $y_{it}$ as follows \[ y_{it}=\mathbf{b}_i'\boldsymbol{\epsilon}_{\cdot t}=\sum_{j=1}^n\mathbf{b}_{ij}'\boldsymbol{\epsilon}_{jt}=\mathbf{v}_{k}'\boldsymbol{\tilde{\epsilon}}_{t}, \] where $\boldsymbol{\tilde{\epsilon}}_{t}=\sum_{j=1}^n a_j\boldsymbol{{\epsilon}}_{jt}$ is a linear combination of the components $\boldsymbol{{\epsilon}}_{jt}$. Given that $\boldsymbol{\tilde{\epsilon}}_{t}$ is also WN, we conclude that $\rho(y_{i})=\rho(y_{ij})=\lambda_k$, as asserted. Lastly, since the unit-ellipse is devoid of boundary-points, the extremal values $\rho_{min}(L)$ and $\rho_{max}(L)$ of the first-order ACF of $y_{it}$ must correspond to stationary points. Thus we have $\rho_{min}(L)=\min_k\lambda_k$ and $\rho_{max}(L)=\max_k\lambda_k$, which concludes the proof of the proposition. \qedsymbol\\

We will now derive the spectral decomposition of the MSE predictor $\boldsymbol{\gamma}_{i\cdot\delta}$. Let $\tilde{\lambda}_m$, $m=1,...,nL$, denote the eigenvalues of the matrix $\mathbf{\tilde{M}}$, with corresponding eigenvectors $\mathbf{\tilde{v}}_m$. Specifically, we have $\tilde{\lambda}_{(k-1)n+ j}=\lambda_k\tilde{\sigma}_j$ and $\mathbf{\tilde{v}}_{(k-1)n+ j}=\mathbf{v}_{\sigma j}\otimes\mathbf{v}_k$, where $\lambda_k$, $\tilde{\sigma}_j$ are the eigenvalues and $\mathbf{v}_k$, $\mathbf{v}_{\sigma j}$ are the eigenvectors of $\mathbf{M}$ and $\boldsymbol{\Sigma}$, respectively, as established by Horn and Johnson (1991). We will index the eigenvalues and eigenvectors of $\tilde{\mathbf{M}}$ using the notation $\tilde{\lambda}_{kj}$ and $\mathbf{\tilde{v}}_{kj}$, where the two-dimensional index $kj$ refers to the product $\lambda_k\tilde{\sigma}_j$. It is assumed that the eigenvalues $\lambda_k$ are ordered in decreasing magnitude, from largest to smallest. Furthermore, we assume that the eigenvectors $\mathbf{\tilde{v}}_{kj}$ are normalized to form an orthonormal basis of $\mathbb{R}^{nL}$. The eigenvalues of $\mathbf{\tilde{I}}$ are $\tilde{\sigma}_j$, each with multiplicity $L$, for $j=1,...,n$ (if $\tilde{\sigma}_j$ are not pairwise different then the multiplicities will be accordingly larger), and the eigenvectors $\mathbf{\tilde{v}}_{kj}$ of $\mathbf{\tilde{M}}$ are also eigenvectors of $\mathbf{\tilde{I}}$. By organizing the orthonormal basis $\mathbf{\tilde{v}}_{kj}$ into the $nL\times nL$-dimensional matrix $\mathbf{\tilde{V}}$, we can express the MSE predictor $\boldsymbol{\gamma}_{i\cdot\delta}$, for $i=1,...,n$, in the form of the spectral decomposition:

where the column vector $\mathbf{w}_i=(w_{ikj})_{k=1,...,L, j=1,...,n}$ of length $nL$ represents the spectral weights of the decomposition of $\boldsymbol{\gamma}_{i\cdot\delta}$. The predictor $\boldsymbol{\gamma}_{i\cdot\delta}$ is said to possess complete spectral support if the condition $\sum_{j=1}^n|w_{ikj}|\neq 0$ holds for $k=1,...,L$. This definition generalizes the concept proposed by Wildi (2024) in the univariate context, where $n=1$ and it is required that $w_{ikj}=w_k\neq 0$ for $1\leq k\leq L$.

Proof: The proof of the first assertion is derived from Proposition (ref). For the second assertion, we assume $i=1$ and analyze the M-SSA Criterion (ref) defined as:

where we assume $l=1$ in the length constraint. Applying spectral decomposition to $\mathbf{b}_1$ yields:

where the two-dimensional index corresponds to the eigenvalues $\tilde{\lambda}_{kj}=\lambda_k\tilde{\sigma}_j$ of $\mathbf{\tilde{M}}$, with $\lambda_k$ ordered from largest to smallest. The variable $\alpha_{kj}=\alpha_{1kj}$ depends on $i=1$, and we omit the complex three dimensional indexing for clarity. Since $\tilde{\mathbf{v}}_{kj}$ serves as an (normalized) eigenvector of $\mathbf{\tilde{I}}$, with strictly positive eigenvalue $\tilde{\sigma}_j>0$, the length constraint can be reformulated as follows:

which we refer to as the elliptic (length) constraint. In a similar manner, the HT constraint (assuming $l=1$) yields the following expression:

which we designate as the hyperbolic HT constraint. We can solve this equation for any $\alpha_{kj}$, at least if $\lambda_k\neq 0$. Let $\lambda_2\neq 0$\footnote{Since $\lambda_k=\cos(k\pi/(L+1))$, it follows that $\lambda_2\neq 0$ when $L>3$, which is virtually always the case in practical applications. For $L=3$ a similar proof applies when selecting $\lambda_3\neq 0$ instead of $\lambda_2$.} so that

where the notation $\sum_{(k,j)\neq(2,1)}$ indicates summation over all combinations of $(k,j)$ excluding the single pairing $(k,j)=(2,1)$. Similarly, we can express $\alpha_{11}$ in the context of the elliptic constraint as follows:

which can be rewritten as

Here, we substitute Equation (ref). This formulation characterizes the intersection of elliptic and hyperbolic constraints: if this intersection is empty, then the M-SSA problem lacks a solution. Assuming $\rho_1=\lambda_1=\max_k\lambda_k$ and isolating $\tilde{\sigma}_1\alpha_{11}^2$ to the left side of the last equation, we derive:

where the notation $\sum_{(k,j)\notin\{(2,1),(1,1:n)\}}$ signifies that the sum extends over all combinations of $(k,j)$ excluding $k=1$ and the single pairing $(k=2,j=1)$. We can rewrite this expression as:

Since $\tilde{\sigma}_j>0$ (full rank) and $\lambda_k<\lambda_j$, if $k>j$, we deduce \[ \tilde{\sigma}_j\frac{\lambda_2-\lambda_k}{\lambda_2-\lambda_1}\left\{

\right.. \] It follows that $\alpha_{kj}=0$ for all pairs $(k,j)$ with $k>2$ since otherwise the condition $\sum_{j=1}^n \tilde{\sigma}_j\alpha_{1j}^2>1$ implied by Equation (ref) would contradict the elliptic (length) constraint specified in Equation (ref). Consequently, we deduce:

where all zero-constraints are attributable to the elliptic constraint when $\boldsymbol{\Sigma}$ is of full rank. It is important to note that if $\rho_1 > \lambda_1$, then the expression \[ \frac{\lambda_2-\rho_1}{\lambda_2-\lambda_1}>1 \] in Equation (ref) would conflict with the elliptic constraint, thereby providing an alternative proof of the first claim (exterior point). In summary, when $\rho_1=\lambda_1$ (the first boundary case), the intersection of the ellipse and hyperbola manifests as a `reduced-rank' $(n-1)$-dimensional ellipse as delineated by Equation (ref). The spectral decomposition expressed in Equation (ref) then simplifies to \[ \mathbf{b}_1=\sum_{j=1}^n \alpha_{1j}\tilde{\mathbf{v}}_{1j}. \] Substituting this expression into the objective function represented in equation (ref) results in \[ \boldsymbol{\gamma}_{1\cdot\delta}'\tilde{\mathbf{I}}\mathbf{b}_{1}=\sum_{j=1}^n w_{11j}\tilde{\sigma}_j\alpha_{1j}, \] leading to the corresponding boundary (`reduced-rank') M-SSA problem defined as:

We can now specify the Lagrangian function \[ \mathcal{L}:=\sum_{j=1}^n w_{11j}\tilde{\sigma}_j\alpha_{1j}-\lambda_{\mathcal{L}}\left(\sum_{j=1}^n \tilde{\sigma}_j\alpha_{1j}^2-1\right) \] where $\lambda_{\mathcal{L}}$ is the Lagrange multiplier. The corresponding Lagrangian equations yield \[ \tilde{\sigma}_jw_{11j}-2\lambda_{\mathcal{L}}\tilde{\sigma}_j\alpha_{1j}=0~,~j=1,...,n \] which leads to the expression \[\alpha_{1j}=w_{11j}/(2\lambda_{\mathcal{L}}). \] Thus we can express \[ \mathbf{b}_1=\frac{1}{2\lambda_{\mathcal{L}}}\sum_{j=1}^n w_{11j}\tilde{\mathbf{v}}_{1j}\neq \mathbf{0} \] with the proportionality term \[ \frac{1}{2\lambda_{\mathcal{L}}}:=\frac{1}{\sqrt{\sum_{j=1}^n\tilde{\sigma}_jw_{11j}^2}} \] ensuring compliance with the elliptic length constraint (assuming $l=1$). A proof of the last assertion (second boundary case) follows a similar reasoning, by symmetry. \qedsymbol\\

Remark: If $|\rho_i|<\rho_{max}(L)$ (an interior point), then the intersection of the ellipse and hyperbola results in a subspace of dimension $Ln-2$ in contrast to the boundary scenario described by Equation (ref), which has dimension $n-1$. We will subsequently derive the corresponding M-SSA predictor for the `interior point' case, which is of primary interest in practical applications. For the purposes of this discussion, we will assume $Ln > 2$ to ensure that the optimization problem remains non-trivial (since otherwise the solution would be determined by the two constraints).

Proof: Consider the Lagrangian defined as

where $\lambda_{i1},\lambda_{i2}$ denote the Lagrange multipliers and we assume $l=1$ in the length constraint. Since the intersection of the ellipse (representing the length constraint) and hyperbola (representing the HT constraint) is devoid of boundary points, the solution $\mathbf{b}_{i}$ to the M-SSA problem must satisfy the Lagrangian equations given by \[ \tilde{\mathbf{I}}\boldsymbol{\gamma}_{i\cdot\delta}=\lambda_{i1}(\mathbf{\tilde{I}}+\mathbf{\tilde{I}}')\mathbf{b}_{i}+\lambda_{i2} (\mathbf{\tilde{M}}+\mathbf{\tilde{M}}')\mathbf{b}_{i}=\lambda_{i1} 2\mathbf{\tilde{I}}\mathbf{b}_{i}+\lambda_{i2} 2\mathbf{\tilde{M}}\mathbf{b}_{i}. \] The second regularity assumption (non-degenerate case) implies that the HT constraint, as expressed in Equation (ref), is `active', leading to $\lambda_{i2}\neq 0$. Dividing through by $\lambda_{i2}$ yields

where $D_i=1/\lambda_{i2}$ and $\nu_i=-2\frac{\lambda_{i1}}{\lambda_{i2}}$. Since the M-SSA solution is derived from the $nL-2>0$ dimensional intersection of the ellipse defined in Equation (ref) and the HT hyperbola characterized in Equation (ref), we can infer that the objective function is not overruled by the constraint, implying that $|\lambda_{i2}|<\infty$. Consequently, $D_i\neq 0$ in Equation (ref) holds true. The eigenvalues of $\mathbf{N}_i$ are given by $\tilde{\sigma}_j(2\lambda_{k}-\nu_i)$, $1\leq k\leq L$, $1\leq j\leq n$. It is important to note that if $\mathbf{b}_i$ is a solution to the M-SSA problem, then $\nu_i/2$ cannot be an eigenvalue of $\mathbf{M}$. To illustrate this, assume $\nu_i/2=\lambda_{k_0}$ for some $k_0\in\{1,...,L\}$. Under this assumption, $\mathbf{N}_i$ would map the eigenvectors $\tilde{\mathbf{v}}_{k_0j}$, for $j=1,...,n$, to zero in Equation (ref). We then look at the left-hand side of the equation: \[ D_i\tilde{\mathbf{I}}\boldsymbol{\gamma}_{i\cdot\delta}=D_i\tilde{\mathbf{I}}\sum_{k=1}^L\sum_{j=1}^nw_{ikj}\tilde{\mathbf{v}}_{kj}=D_i\sum_{k=1}^L\sum_{j=1}^n\tilde{\sigma}_jw_{ikj}\tilde{\mathbf{v}}_{kj}. \] If $\tilde{\mathbf{v}}_{k_0j}$, $j=1,...,n$, is mapped to zero, then $D_i\tilde{\sigma}_jw_{ik_0j}=0$, $j=1,...,n$. Consequently, $w_{ik_0j}=0$, $j=1,...,n$, (because $D_i\tilde{\sigma}_j\neq 0$), thereby contradicting the last regularity assumption (completeness). Thus, we conclude that $\nu_i/2$ cannot be an eigenvalue of $\mathbf{M}$, $\nu_i\in \mathbb{R}-\{2\lambda_k|k=1,...,L\}$, $\mathbf{N}_i^{-1}$ exists and \[ \mathbf{N}_i^{-1}=\mathbf{\tilde{V}}\mathbf{D}_{\nu i}^{-1}\mathbf{\tilde{V}}'. \] The diagonal matrix $\mathbf{D}_{\nu i}^{-1}$ is characterized by well-defined entries of the form $1/(\tilde{\sigma}_j(2\lambda_{k}-\nu_i))$. We can proceed to solve Equation (ref) for $\mathbf{b}_i$, yielding the following expression:

The proof of the second assertion is derived from Equation (ref) and relies on the orthonormality of $\tilde{\mathbf{v}}_{kj}$:

It is important to note that the first-order autocorrelation $\rho(\nu_i)$ of $y_{it}$ is influenced by $\boldsymbol{\Sigma}$, since the eigenvalues $\tilde{\sigma}_j$ do not cancel in Equation (ref). Consequently, selecting $\nu_i = \nu_{0i}$ such that $\rho_i(\nu_{0i}) = \rho_i$ (in accordance with the HT constraint) indicates that the SSA solution $\mathbf{b}_i(\nu_{0i})$ generally depends on $\boldsymbol{\Sigma}$. This is despite the fact that, as specified by Equation (ref), the expression for $\mathbf{b}_i$ appears to be independent of $\boldsymbol{\Sigma}$, owing to the cancellation of the eigenvalues $\tilde{\sigma}_j$ in this expression. As $\nu_i$ approaches $2\lambda_k$, we find that $\rho(\nu_i)$ converges to $\lambda_k$, for $1\leq k\leq L$, given the assumption that $\sum_{j=1}^{n}\tilde{\sigma}_jw_{ikj}^2>0$ for $k=1,...,L$. Thus, $\rho(\nu_i)$ can come arbitrarily close to the boundary values $\pm\rho_{max}(L)$. Therefore, the continuity of $\rho(\nu_i)$ along with the intermediate value theorem implies that any value $\rho_i$ satisfying $|\rho_i|<\rho_{max}(L)$ is permissible in the HT constraint, as asserted. \\ We now proceed to prove Assertion (ref). For this purpose, we define $\tilde{w}_{ik}^2:=\sum_{j=1}^{n}\tilde{\sigma}_jw_{ikj}^2$ in Equation (ref), resulting in the expression

utilizing the fact that $\sum_{k=1}^{L} \frac{1}{(2\lambda_k - \nu_i)^2} \tilde{w}_{ik}^2 = \mathbf{b}_i' \mathbf{\tilde{I}} \mathbf{b}_i = 1$ to derive the last equality. For the case when $k=j$, the terms cancel in Equation (ref). Therefore, we assume $k\neq j$ and consider the expression:

where we interchanged $k$ and $j$ in the last equality, owing to symmetry. Assuming $k\neq j$ and considering the addition of the two symmetric terms corresponding to $(k,j)$ and $(j,k)$, we derive the following expression:

where the strict inequality is justified by the conditions $|\nu_i|>2\rho_{max}(L)\geq 2|\lambda_k|$, $\lambda_k\neq \lambda_j$, and the assumption $\tilde{w}_{ik}^2=\sum_{j=1}^{n}\tilde{\sigma}_jw_{ikj}^2>0$. Consequently, We infer that \[ \frac{\partial}{\partial\nu_i}\rho(\nu_i)<0. \] Furthermore, from Equation (ref) we can deduce the limit \[ \lim_{|\nu_i|\to \infty}\rho(\nu_i)= \displaystyle{\frac{\sum_{k=1}^{L}\lambda_k\left(\sum_{j=1}^{n}\tilde{\sigma}_jw_{ikj}^2\right)}{\sum_{k=1}^{L}\left(\sum_{j=1}^{n}\tilde{\sigma}_jw_{ikj}^2\right)}}=\displaystyle{\frac{\boldsymbol{\gamma}_{i\cdot\delta}'\mathbf{\tilde{M}}\boldsymbol{\gamma}_{i\cdot\delta}}{\boldsymbol{\gamma}_{i\cdot\delta}'\mathbf{\tilde{I}}\boldsymbol{\gamma}_{i\cdot\delta}}}=\rho_{i,MSE}, \] which represents the first-order ACF of the MSE predictor. In conjunction with the identified strictly negative derivative, we conclude that $\rho_i(\nu_{1i})>\rho_i(\nu_{2i})$ for $\nu_{1i} > 2\rho_{max}(L)$ and $\nu_{2i}< -2\rho_{max}(L)$, as asserted.\\ For a proof of the last Assertion (ref) we analyze the expression

where we used $\mathbf{N}_i^{-1}~'\mathbf{N}_i^{-1}=\mathbf{N}_i^{-1}\mathbf{N}_i^{-1}=:\mathbf{N}^{-2}$, owing to symmetry. For the case where $\nu<-2\rho_{max}(L)$, the quotient is strictly positive. Consequently, the positivity of the objective function for the M-SSA solution implies that $\textrm{sign}(D_i)=-\textrm{sign}(\nu_i)=1$. Conversely, for $\nu>2\rho_{max}(L)$ the quotient becomes strictly negative, leading to $\textrm{sign}(D_i)=-\textrm{sign}(\nu_i)=-1$. Assuming $\nu<-2\rho_{max}(L)$ we deduce:

where we inserted $-\textrm{sign}(\nu)=1$. The last expression is obtained by recognizing that

where $\mathbf{N}_i^{-k}:=(\mathbf{N}_i^{-1})^k$, $(\mathbf{N}_i^{-1})'=\mathbf{N}_i^{-1}$ (by symmetry). The commutativity of the matrix product utilized in deriving the third equality arises from the simultaneous diagonalization of both $\mathbf{\tilde{M}}$ and $\mathbf{N}_i^{-1}$, which share the same eigenvectors. Additionally, standard rules for matrix differentiation were employed in the derivation of Equation (ref)\footnote{$\frac{d(\mathbf{N}_i^{-1})}{d\nu_i}=\mathbf{N}_i^{-2}$ and $\frac{d(\mathbf{N}_i^{-2})}{d\nu_i}=2\mathbf{N}_i^{-3}$. The first equation follows from the general differentiation rule $\frac{d(\mathbf{N}_i^{-1})}{d\nu_i}=-\mathbf{N}_i^{-1}\frac{d\mathbf{N}_i}{d\nu_i}\mathbf{N}_i^{-1}$, noting that $\frac{d\mathbf{N}_i}{d\nu_i}=-\mathbf{\tilde{I}}$. The second equation is derived by substituting the first equation into $\frac{d(\mathbf{N}_i^{-2})}{d\nu_i}=\frac{d(\mathbf{N}_i^{-1})}{d\nu_i}\mathbf{N}_i^{-1}+\mathbf{N}_i^{-1}\frac{d(\mathbf{N}_i^{-1})}{d\nu_i}$.}. Finally, the expression $2\mathbf{\tilde{M}}\mathbf{N}_i^{-k}= \mathbf{N}_i^{-k+1}+\nu\mathbf{N}_i^{-k}$ was incorporated into the numerator of Equation (ref), leading to the final equation after simplification. The aforementioned proof is also applicable in the scenario where $\nu > 2\rho_{\max}(L)$, albeit with a change in sign, such that $\text{sign}(D_i) = -\text{sign}(\nu) = -1$, as was to be demonstrated. \qedsymbol\\

Remarks\\ As $|\nu_i|\to \infty$ and $D_i=-\nu_i$, it follows that $-\mathbf{N}_i/\nu_i\to \mathbf{\tilde{I}}$. Consequently, we have \[ \mathbf{b}_i\to-\frac{D_i}{\nu_i}\tilde{\mathbf{I}}^{-1} \tilde{\mathbf{I}}\boldsymbol{\gamma}_{i\cdot\delta}=\boldsymbol{\gamma}_{i\cdot\delta}. \] This indicates that the M-SSA method asymptotically replicates the MSE predictor in the limit as $|\nu_i| \to \infty$, reflecting an asymptotically degenerate case where the HT constraint may be omitted. By Assertion (ref), we conclude that M-SSA can accommodate HT constraints $\rho_i<\rho_{i,MSE}$ when $\nu_i\in]-\infty,-2\rho_{max}(L)]$. Conversely, for $\nu_i\in[2\rho_{max}(L),\infty[$ M-SSA is capable of accommodating HT constraints $\rho_i>\rho_{i,MSE}$ (smoothing). Moreover, Equation (ref) elucidates a trade-off between Accuracy (target correlation) and Smoothness (first-order ACF) in the context of solving the M-SSA optimization problem. For a more comprehensive discussion, refer to Wildi (2024). Finally, Wildi (2026) also addresses the singular case of incomplete spectral support, which may be generalized to M-SSA; however, this topic is omitted here for the sake of brevity (and lack of practical relevance). \\

Proof: The first claim can be substantiated by multiplying both sides of Equation (ref) by $\mathbf{\tilde{I}}^{-1}=\boldsymbol{\Sigma}^{-1}\otimes \mathbf{I}_{L\times L}$. This yields the following equations:

A proof follows from the block-diagonal structure of $\mathbf{I}_{n\times n}\otimes\mathbf{{M}}$. The validity of the second claim follows directly from Equation (ref), achieved by multiplying both sides by $\boldsymbol{{\nu}}_i$. Specifically, the boundary equations $b_{ij1}-\nu_i b_{ij0}=D_i\gamma_{ij\delta}$ and $-\nu_i b_{ij(L-1)}+b_{ij(L-2)}=D_i\gamma_{ij(L-1+\delta)}$ at $k=0$ and $k=L-1$ substantiate the implicit boundary conditions $b_{ij,-1}=0$ and $b_{ijL}=0$.\\

Theorem (ref) presents a one-parameter formulation of the SSA solution, and the subsequent corollary elucidates this solution by establishing a relationship between the unknown parameter $\nu_i$ and the HT constraint.

A proof follows directly from Theorem (ref), with the observation that the scaling factor $s_i = \sqrt{l / \mathbf{b}_i(\nu_{0i})' \tilde{\mathbf{I}} \mathbf{b}_i(\nu_{0i})}$ does not interfere with either the objective function or the HT constraint. This scaling can be applied subsequently, after obtaining a solution under the assumption of an arbitrary scaling such that $|D_i| = 1$. Notably, assertion (ref) of the theorem guarantees the uniqueness of solutions within the set $\{\nu | |\nu| > 2 \rho_{\text{max}}(L)\}$, owing to the bijectivity of the first order ACF. \qedsymbol\\

Remark: The objective function $\rho(y_i,{z}_{i\delta},\delta)$ and the first-order ACF $\rho(y_{i})$ are generally influenced by $\boldsymbol{\Sigma}$. Interestingly, Equation (ref) indicates that the one-parameter form of the solution $\mathbf{b}_i(\nu_i)$ is independent of $\boldsymbol{\Sigma}$; this is corroborated by the cancellation of $\tilde{\sigma}_j$ in Equation (ref). However, the optimal solution $\mathbf{b}_i(\nu_{0i})$, which is derived from the optimal $\nu_{0i}$, remains dependent on $\boldsymbol{\Sigma}$ because $\tilde{\sigma}_j$ (for $j = 1, \ldots, n$) does not cancel in the first-order ACF (ref), which is critical for determining $\nu_{i0}$ within the context of the nonlinear HT constraint.\\

The next result addresses the distribution of the M-SSA predictor.

The proof follows directly from Equation (ref), acknowledging the symmetry of $\mathbf{N}_i^{-1}$. For a comprehensive derivation of the mean, variance, and (asymptotic) distribution of the MSE estimate, one may consult standard textbooks, such as Brockwell and Davis (1993).\\

The final result in this section introduces a dual reformulation of the M-SSA optimization criterion, which characterizes its solution as the smoothest predictor among all predictors that achieve the same tracking accuracy or target correlation. To facilitate this, we introduce some notation:

is the output of the filter with weights $\boldsymbol{\gamma}_{i\cdot\delta}^{\rho_{max}}:=\sum_{j=1}^n w_{i1j}\mathbf{\tilde{v}}_{1j}$, corresponding to the projection of the MSE predictor onto the subspace spanned by $\mathbf{\tilde{v}}_{11},...,\mathbf{\tilde{v}}_{1n}$. Under the regularity assumptions of Theorem (ref), specifically a complete spectral support, $\boldsymbol{\gamma}_{i\cdot\delta}^{\rho_{max}}\neq\mathbf{0}$. According to the proof of the second assertion in Proposition (ref), $\mathbf{y}_{it,MSE}^{\rho_{max}}$ maximizes the target correlation under the (extremal) HT constraint $\rho(\mathbf{y}_{it,MSE}^{\rho_{max}})=\rho_{max}(L)$ (boundary point). Consequently, if $y_{it}$ is such that $\rho(y_i,z_i,\delta)>\rho({y}_{i,MSE}^{\rho_{max}},z_i,\delta)$ for its target correlation, then $\rho(y_i)<\rho_{max}(L)$ for its first-order ACF (interior point). Moreover, under the assumption of complete spectral support, $\boldsymbol{\gamma}_{i\cdot\delta}^{\rho_{max}}\neq \boldsymbol{\gamma}_{i\cdot\delta}$ (since the latter's spectral weights $w_{ikj}$ cannot all vanish for $k>1$) and therefore the strict inequality $\rho({y}_{i,MSE}^{\rho_{max}},z_i,\delta)<\rho(y_{i,MSE},z_i,\delta)$ holds for the target correlation of the MSE predictor, with filter weights $\boldsymbol{\gamma}_{i\cdot\delta}$, owing to uniqueness of the MSE predictor.\\

Proof: Given the similarity in the problem structure, we may leverage the arguments employed in the proof of Theorem (ref) to verify the above claims. Accordingly, we here focus the analysis on the relevant deviations. The Lagrangian $\tilde{\mathcal{L}}_i$ of the dual problem is expressed as:

where we inserted $\rho(y_i,z_i,\delta)=\frac{\boldsymbol{\gamma}_{i\cdot\delta}'\tilde{\mathbf{I}}\mathbf{b}_{i}}{\sqrt{l\boldsymbol{\gamma}_{i\cdot\delta}'\mathbf{\tilde{I}}\boldsymbol{\gamma}_{i\cdot\delta}}}$ in the target correlation constraint. The two new regularity conditions ensure that the solution to the dual criterion is neither an unconditional maximum of the new objective function (non-degeneration) nor the maximum in the target correlation of the new constraint (interior point), see the reasoning following Eq.(ref). As a result, Lagrange multipliers are finite and non-vanishing. After differentiation of the Lagrangian $\tilde{\mathcal{L}}_i$, Equation (ref) is obtained, after re-arranging terms\footnote{The dependence of $\tilde{D}_i,\tilde{\nu}_i$ on the Lagrangian multipliers $\tilde{\lambda}_{i1},\tilde{\lambda}_{i2}$ differs slightly from that in the primal problem, a discrepancy that can be rectified through straightforward re-scaling, given that the multipliers are finite and non-vanishing.}. The feasible solution space is now defined by the intersection of the length-ellipse—corresponding to the original length constraint—and the target correlation plane—representing the new target correlation constraint (in contrast to the intersection of the ellipse and hyperbola in the primal problem).\\ For the proof of the second assertion, we first note that the main distinction between original M-SSA and dual criteria concerns the selection of $\tilde{\nu}_i$, such that it satisfies the new target correlation constraint $\rho(y_i(\tilde{\nu_i}),z_i,\delta)=\rho(y_i(\nu_{0i}),z_i,\delta)$ (instead of the original HT constraint). This distinction is significant insofar as the optimal $\tilde{\nu}_{i0}$ of the dual problem is generally not uniquely specified by the hyperparameter $\rho_{iyz}$ associated with the target correlation constraint\footnote{Unlike the first-order ACF, the target correlation is not a strictly monotonic function for $\nu_i\in\{\nu||\nu|>2\rho_{max}(L)\}$, since the sign of its derivative depends on the sign of $\nu_i$, see Equation (ref).}, in contrast to the primal (original M-SSA) problem, where the parameter $\nu_{0i}$ is uniquely determined by $\rho_i$ in the HT constraint (at least under the posited assumptions, see Corollary (ref)). However, we infer from the similar formal structure of the problem, as expressed by Equation (ref), that the strict monotonicity of the target correlation $\rho(y_i(\tilde{\nu_i}),z_i,\delta)$, when $\tilde{\nu}_i\in\{\nu|\nu>2\rho_{max}(L)\}$ or when $\tilde{\nu}_i\in\{\nu|\nu<-2\rho_{max}(L)\}$ (Assertion (ref) of Theorem (ref)), extends to the dual problem. To proceed with the proof of the second claim, we then initially assume that the search space for $\tilde{\nu}_i$ in the dual problem can be restricted to the set $\{\nu|\nu>2\rho_{max}(L)\}$. In this case, the solutions to the primal and dual problems must coincide due to strict monotonicity of $\rho(y_i(\tilde{\nu_i}),z_i,\delta)$, which guarantees the uniqueness of the solution. The requested extension of this result to the set $\{\nu||\nu|>2\rho_{max}(L)\}$ is further supported by Assertion (ref) of Theorem (ref), which applies in analogous manner to the dual problem and which establishes that $\rho(\tilde{\nu}_i)<\rho(\nu_{0i})$, when $\tilde{\nu}_i<-2\rho_{max}(L)$ (since $\nu_{0i}>2\rho_{max}(L)$, by assumption). Consequently, the maximum of the objective function of the dual criterion can be confined to the region $\{\nu|\nu>2\rho_{max}(L)\}$. \\ A similar line of reasoning applies to the proof of the last claim, where $\nu_{0i}<-2\rho_{max}(L)$ in the primal problem, with the necessary adjustment that the maximization in the dual Criterion (ref) must be replaced with minimization, as now $\rho(\tilde{\nu}_i)>\rho(\nu_{i0})$ if $\tilde{\nu}_i>2\rho_{max}(L)$\footnote{If the maximization were retained in the objective function of the dual criterion, it would generally result in swapping the solution $\nu_{0i}<-2\rho_{max}(L)$ of the primal problem from the left branch $\{\nu|\nu<-2\rho_{max}(L)\}$ to the right branch $\{\nu|\nu>2\rho_{max}(L)\}$, due to maximization of the first-order ACF. Such a swap would lead to a different solution to the dual problem in this case.}. \qedsymbol\\

Remark: Wildi (2026) demonstrates that the condition $|\nu_i|>2\rho_{max}(L)$, ensuring strict monotonicity and uniqueness in the previous results, is not a limitation in practical applications. Consequently, we may assume this condition holds except in highly unusual problem specifications.\\

Under the posited assumptions, the above theorem characterizes the M-SSA predictor as the smoothest predictor—exhibiting the lowest zero-crossing rate—for a given level of tracking accuracy (target correlation). This property distinguishes the M-SSA approach as a compelling alternative to traditional smoothing methods, as demonstrated in the next section.

Applications

Our examples highlight three key aspects: smoothing, forecasting, and nowcasting (real-time signal extraction). Smoothing is primarily concerned with causal targets, whereas traditional forecasting typically focuses on acausal, forward-looking allpass filters. Signal extraction encompasses various acausal filtering techniques, including lowpass filters for trend extraction, bandpass filters for cycle analysis, highpass filters for separating cycles from noise, and seasonal adjustment.

Forecasting

We apply M-SSA to a stationary VAR(1) process represented by the equation: \[ \left(

\right)=\mathbf{A}_1\left(

\right)+\left(

\right) , \mathbf{A}_1=\left(

\right) , \boldsymbol{\Sigma}=\left(

\right). \] We compute one-step ahead predictors for $z_{1t+1}$ and $z_{2t+1}$ with $\delta=1$ and a predictor length of $L=100$. In the classic forecast application, $\mathbf{z}_t=\mathbf{x}_t$, $\boldsymbol{\Gamma}_k=\left\{

\right.$ and the target $\mathbf{z}_{t+\delta}$ relies on an anticipative all-pass filter, for $\delta>0$. The HTs within the constraint are defined as $ht_1=3$, $ht_2=8$. The VAR process allows for a convergent MA-inversion with weights $\boldsymbol{\Xi}_k=\mathbf{A}_1^k$. Accordingly, the HTs of the MSE predictor can be computed by substituting the first-order ACFs $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}'\mathbf{\tilde{M}}(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}/(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}'\mathbf{\tilde{I}}(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}$, $i=1,2$, $\delta=1$, into Equation \eqref{ht} to yield $ht_{1MSE}=5.6$ and $ht_{2MSE}=4.6$. Here, $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}$ represents the convolution of the target and the MA-inversion, as discussed in Section \eqref{ext_stat}. Given that the target is a forward-shifting allpass filter, we conclude that $(\boldsymbol{\gamma}\cdot\boldsymbol{\xi})_{i\delta}$ corresponds to the shifted original MA-inversion, with the leading term at $-\delta=-1$ omitted (the forecast of the future noise term $\epsilon_{t+1}$ vanishes). \\

With our specifications we observe $ht_1<ht_{1MSE}$ , indicating that the M-SSA predictor must introduce additional zero crossings to `unsmooth' the MSE benchmark for the first series. Conversely, for the second series, $ht_2>ht_{2MSE}$ suggests that the M-SSA predictor must `smooth' the MSE benchmark by producing fewer zero crossings. M-SSA and MSE predictors are displayed in Fig.(ref), with the M-SSA predictor presented in the top panel and the MSE predictor in the bottom panel. The weights associated with the MSE predictor are defined as $\mathbf{A}_1^{\delta}=\mathbf{A}_1$. The unsmoothing effect is characterized by an alternating high frequency pattern exhibited by the M-SSA filter in the first series (top left panel). In contrast, the smoothing effect results in slowly monotonically decaying patterns that are observed in the second series (top right panel).

For confirmation, we proceed to compute sample estimates of the relevant performances measures for a lengthy sample, $N=100000$, which ensures a near convergence of empirical to expected measures, as presented in Table (ref). The final column of the table corroborates the HTs of the MSE predictor, as derived previously. The optimal values of $\nu_i$ specified in Corollary (ref) are $\nu_1$=-2.034$\in]-\infty,-2\rho_{max}(L)]$ (indicating unsmoothing) and $\nu_2$=2.001$ \in [2\rho_{max}(L),\infty[$ (indicating smoothing), as alluded to by the remark after the proof of Theorem (ref).

A comparative analysis of the sample predictors and targets is presented in Figure (ref), which demonstrates the effects of `unsmoothing' in the first series and smoothing in the second series, as dictated by the HT constraints.

Smoothing

If the target correlation in the previous example is computed against the acausal target $\mathbf{z}_{t+1}$, then M-SSA functions as a forecasting method. However, as discussed in Section (ref), the MSE predictor can be utilized in place of the effective target without impacting the optimization results. In this context, M-SSA is then considered a smoother for the causal MSE predictor, with smoothing referring specifically to the application of M-SSA to causal targets. In this context, we can examine the well-established Whittaker-Henderson (WH) graduation, expressed as follows:

where $\Delta^D$ represents the common difference operator applied $D$-times. Notably, for $D=2$, this formulation yields the Hoddrick-Prescott (HP) filter, as detailed by Hodrick and Prescott (1997). The equation delineates a tradeoff between conflicting objectives of data-fitting and data-smoothing terms, with the regularization parameter $\lambda$ serving to balance this trade-off. An alternative formulation (dilemma) can be formalized within the (univariate) SSA framework as follows:

where, in contrast to the previous section, $\delta\leq 0$. In this setting, the target $z_{t+\delta}=x_{t+\delta}$ corresponds to the output of a causal back-shifting all-pass filter. If $L$ is defined as an odd integer and we select $\delta=-(L-1)/2$, this ensures symmetry in the backcast; the coefficients of the causal filter are centered at $x_{t-(L-1)/2}$, with the tails reflecting symmetrically around this central point. For the range $-(L-1)/2\leq \delta\leq 0$, the design transitions from symmetry to asymmetry as it approaches the nowcast at $\delta=0$. To illustrate, Fig.(ref) presents the (univariate) SSA smoothers $\mathbf{b}_{\delta}$ based on Criterion (ref), scaled arbitrarily to unit length and plotted as a function of $\delta$. This analysis assumes WN data, with $L=201$ and $\rho_1=0.9986$, which is inserted into the HT constraint of Criterion (ref). Here, $\rho_1$ represents the first-order ACF of the two-sided HP filter with parameter $\lambda=14400$, as derived from the WH criterion (ref) and utilized in the subsequent Section (ref). This selection of $\lambda$ aligns with applications to monthly time series, as discussed by Ravn and Uhlig (2002). \\

The SSA smoothers in this example replicate exactly the HT of HP, regardless of $\delta$. Consequently, by Theorem (ref), they are expected to produce superior target correlations or achieve lower MSEs compared to HP (assuming white noise data). Conversely, for similar target correlation, we expect SSA to yield superior HT when compared to HP, owing to Theorem (ref) (duality). For illustration, Figure (ref) presents a truncated symmetric causal HP(14400) filter alongside two distinct symmetric SSA smoothers, with the assumption that $\delta = -\frac{(L-1)}{2}$. Additionally, Table (ref) provides a comparative analysis of the performance of all three designs, assuming the data follows a WN sequence.

For identical HT, the SSA1 design (represented by the blue line in the figure) demonstrates superior performance in target correlation compared to HP, as claimed. Similarly, the SSA2 design (depicted by the violet line in the figure) exhibits a maximal HT, surpassing the HP filter, given an equivalent target correlation. However, when evaluating smoothness based on the curvature (mean-squared second order differences), HP outperforms both SSA smoothers. The observed variations across the reported performance metrics are sufficiently substantial to warrant a carefully informed decision between minimizing the rate of zero-crossings with SSA or reducing curvature via WH. In summary, the simplified SSA optimization problem (ref) can be regarded as a novel smoothing algorithm that sets itself apart from traditional methods such as the WH or HP approaches derived from Criterion (ref). This distinction stems from the fact that the regularization term is explicitly formulated as a HT constraint, providing an alternative perspective on controlling smoothness in time series analysis. We argue that regulating the rate of sign changes offers a viable alternative to classical filtering and smoothing techniques, particularly within economic contexts characterized by alternating growth phases, where it aligns more closely with the underlying decision-making and control logic, as discussed in Wildi (2024).\\

The univariate smoothers in the above example are determined by the assumption of WN. We now extend the framework to multivariate autocorrelated data, as discussed in Section (ref). For illustration, consider the following three-dimensional VAR(1), \[ \mathbf{x}_t=\left(

\right)\mathbf{x}_{t-1}+\boldsymbol{\epsilon}_t , \boldsymbol{\Sigma}=\left(

\right) \] with Wold-decomposition $\boldsymbol{\Xi}_k=\mathbf{A}^k$, assuming a causal allpass target $z_t=x_t$ with $\delta=0$, which corresponds to the proper identity. We select $L=51$ and set the HTs $ht_1=8$, $ht_2=6$ and $ht_3=10$ within the HT constraints for the three M-SSA smoothers $y_{1t},y_{2t}$ and $y_{3t}$, as illustrated in Fig.(ref). The plots indicate that the second series $x_{2t}$ (in red) serves as a significant explanatory variable for all three target series (after suitable adjustment of scales). Fig.(ref) provides insight into this phenomenon, as the peaks of the CCF shown in the left panel, along with the original time series in the right panel, indicate that $x_{2t}$ is left-shifted or leading. This finding can be interpreted as supporting the importance of a leading indicator within a multivariate nowcasting framework. Table (ref) presents expected and sample performance measures, the latter in parentheses, based on an extensive sample size $N=100000$, which corroborates the posited constraints. Additionally, the HTs derived from the original (unfiltered) data, detailed in the last column, identify the second series as the smoothest among the three. The double strike of a lead as well as an increased smoothness of $x_{2t}$ further reinforce its importance as an explanatory variable in the context of endpoint smoothing (or forecasting). \\

Fig.(ref) presents a comparison between $x_{it}$ (black) and M-SSA endpoint smoothers $y_{it}$ (cyan), for $i=1,2,3$, thereby confirming the previously discussed findings. Specifically, the smoothing task proves to be significantly less challenging for $x_{2t}$ relative to $x_{3t}$, as evidenced by the more substantial difference between the HTs of the original data and the smoothers in the latter case. This observation is further supported by the last two columns of Table (ref). It is worth noting that if we selected $\delta=-(L-1)/2$ in the above example, instead of $\delta=0$, the resulting smoothers $y_{it}$ would simplify to (nearly) symmetric and univariate designs (not shown).

Signal Extraction

The M-SSA criterion is capable of addressing both uni- and multivariate forecasting and smoothing tasks. Here, we consider signal extraction, allowing for both the target and the MA inversion to deviate from allpass (identity) designs, representing the highest level of complexity for M-SSA. To illustrate this, we follow the business-cycle application outlined in Wildi (2024), which applies the HP(14400) filter to the monthly U.S. industrial production index\footnote{Board of Governors of the Federal Reserve System (US), Industrial Production: Total Index [INDPRO], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/INDPRO, October 23, 2024.}. We extend the scope of this analysis by incorporating an additional leading indicator, the Composite Leading Indicator (CLI) provided by the OECD\footnote{Organization for Economic Co-operation and Development, Leading Indicators OECD: Leading Indicators: Composite Leading Indicator: Amplitude Adjusted for United States [USALOLITOAASTSAM], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/USALOLITOAASTSAM, October 23, 2024.}, as depicted in the top left panel of Fig.(ref).

Given that the CLI fluctuates around a fixed level, we now consider the first differences of the log-transformed data. The top right panel of the figure presents the corresponding standardized time series, while the bottom left panel displays the actual data utilized for fitting. In this latter panel, singular pandemic outliers have been removed (trimmed) to ensure that absolute values do not exceed $5\sigma$, a threshold that aligns with extreme observations from earlier historical data. The peak of the CCF in the lower right panel indicates that the CLI leads the INDPRO by one quarter. We examine the application of the HP trend filter to differenced series, with a focus on the growth cycle to track U.S. recession episodes as identified by the National Bureau of Economic Research (NBER)\footnote{In contrast, the conventional output gap is defined as the difference between the actual (original non-stationary) time series and its trend, which addresses deviations from `potential output'. A comparative analysis of both designs is presented in Wildi (2024). However, to conserve space, this discussion will utilize the single growth cycle concept to illustrate the M-SSA.}. We then compare univariate and bivariate designs, highlighting target correlations, holding times (HTs) and leads (left shifts/advancement) or lags (right-shifts/retardation) of the corresponding nowcasts. Specifically, besides the bivariate MSE nowcast we consider the univariate SSA and the classic concurrent HP, as discussed in McElroy (2006), denoted as HP-C, as additional benchmarks. \\

To ensure clarity and parsimony, we employ a simple bivariate model, excluding data revisions by relying on the final releases of both indicators. Additionally, the model incorporates all available observations.\footnote{Sample truncation may affect model identification and parameter estimation, although the impact on the target correlation and the HT of the resulting filters remains limited unless truncation is severe (e.g., fewer than half of the available observations).} While a more comprehensive analysis would likely yield improved forecast performance and potentially insightful out-of-sample conclusions, our primary objective is to illustrate the accuracy-smoothness trade-off based on the M-SSA baseline extension discussed in Section (ref), avoiding unnecessary complications from modeling issues pertinent to the baseline (pseudo-)MSE specification. Heinisch et al. (2026) present a more complex five-dimensional M-SSA predictor for German GDP, emphasizing explicit modeling considerations; the corresponding empirical framework is available in our open source M-SSA package for further investigation.

Univariate: HP-C and SSA

We begin by deriving a univariate SSA framework, with the filter length set to $L=201$, see Wildi (2024). This analysis yields the following ARMA(2,1) model for the transformed INDPRO series illustrated in the bottom left panel of Fig.(ref):

According to standard diagnostic tests, the residiuals are white noise. This model is employed to derive the MA inversion $\xi_k$ for $k=0,...,L-1=200$, as illustrated in the left panel of Fig. (ref). Our acausal objective $z_{t+\delta}$ utilizes the (truncated) two-sided Hodrick-Prescott filter with $\lambda = 14400$, as illustrated in the middle panel of the figure. The (univariate) MSE nowcast of $z_{t+\delta}$, where $\delta=0$, is derived following McElroy and Wildi (2016)\footnote{An open source R-package is provided by McElroy and Livsey (2022) with further discussion available in McElroy, T. (2022).}

where $B$ denotes the backshift operator, $\boldsymbol{\xi}(B)=\sum_{k\geq 0}\xi_k B^k$, $\boldsymbol{\xi}^{-1}(B)$ represents the AR-inversion of the data generating process and $\gamma_{k}$, $|k|<\infty$, are the weights of the two-sided HP filter. The notation $[\cdot]_{|k|}^{\infty}$ signifies the omission of the first $|k|-1$ lags\footnote{The derivation of the MSE benchmark follows three distinct steps: first, compute the convolution of the (finite length) two-sided target filter $\gamma_k$ and the MA inversion $\xi_k$ of the ARMA model; second, truncate this expression by omitting the acausal terms associated with all future innovations $\epsilon_{t-k}$ for $k < 0$; lastly, deconvolute the MSE filter from the truncated convolution, as specified in Equation (ref). The M-SSA package provides comprehensive guidance on the sequence of transformations and operations.}. The resulting filter is depicted in the right panel of Fig. (ref) (green line). \\

The MSE nowcast filter exhibits a first-order ACF of 0.963, which corresponds to a HT of 11.508. Similarly, the classic HP-C filter has first-order ACF 0.967, which corresponds to a HT of 12.267. The (univariate) SSA, as illustrated in Fig.(ref) (blue line, right panel), is derived from the MSE nowcast, as outlined in Wildi (2024) (the proceeding corresponds to Section (ref) for the univariate case). Within the HT constraint, we set $ht_1=17.26$, which surpasses the MSE benchmark (and incidentally HP-C) by 50$\%$. As a result, we anticipate that the SSA nowcast will yield a reduction of 33$\%$ in zero-crossings over the long term. While this choice of the HT constraint is somewhat arbitrary, our selection is mainly intended to illustrate the accuracy–smoothness trade-off.

A comparison of SSA and HP-C nowcasts in Fig.(ref) reveals that the SSA exhibits increased smoothness, albeit at the cost of a one-month rightward shift or retardation relative to HP-C (sample performances are reported in Tables (ref) and (ref)). Thus , we anticipate that the bivariate leading indicator design, which incorporates the additional CLI, will mitigate this lag while maintaining comparable smoothness.

M-SSA: Bivariate Leading Indicator Design

We examine a bivariate design that incorporates the CLI as an additional explanatory variable for the industrial production index. The following VARMA(3,1) model is specified as:

with the covariance matrix given by $\boldsymbol{\Sigma} = \left(

\right)$, where the first series corresponds to INDPRO. The bivariate residuals are confirmed to exhibit white noise characteristics according to standard diagnostic tests. Subsequently, we focus attention on a INDPRO nowcast, discarding a corresponding analysis for the CLI. For consistency, we apply the same HT constraint as for the SSA in the previous section.

A simulation experiment utilizing an extensive sample of one million observations from the aforementioned VARMA process corroborates the convergence of sample performance metrics to their expected values, as detailed in Table (ref). The M-SSA exhibits the largest HT, which is $50\%$ greater than that of the MSE predictor (as constrained), and the second-largest target correlation, being only marginally outperformed by the MSE predictor in this dimension. It also outperforms the univariate benchmark HP-C on both accounts. It is important to emphasize that the results presented in the table are contingent upon the bivariate VARMA model. Specifically, for the classic HP-C the values would differ if the univariate ARMA model from the previous section were utilized instead. The derivation of these `true' performance metrics, presented in the first row of the table, is discussed in the appedix. In any case, the M-SSA package provides comprehensive guidance on the sequence of operations and transformations necessary to replicate all results. Note also that we exclude the univariate MSE from the subsequent benchmark comparisons since it is outperformed by the bi-variate (MSE) design.\\

The MA-inversion (impulse response function) of the VARMA-model, along with the nowcast weights, is presented in Fig.(ref). They indicate that the CLI serves as an important explanatory variable for INDPRO. The rightward shift of the majority of the CLI MA-weights (top left) further substantiates the leading characteristic of the indicator.

The sample and expected performances of the classic univariate HP-C filter and the univariate SSA, alongside the bivariate MSE and M-SSA nowcasts are compared in Tables (ref) (target correlations) and (ref) (HTs), with standard errors indicated in parentheses. All expected values are based on the bivariate VARMA model, except for the SSA (in the second columns of the tables), which are derived from the univariate ARMA model (therefore, direct comparisons of expected numbers with the SSA should be approached with caution). The sample estimates are derived from data spanning from 1955-02-01 to 2024-09-01, resulting in a net number of observations of $N=605$\footnote{We utilize the full available common sample of INDPRO and CLI to mitigate the impact of excessive variance in the sample HT. The net sample length is calculated as $T-L-L^{left-tail}=836-201-30=605$, where $L$ represents the filter length (initialization for all designs) and $L^{left-tail}$ denotes the length of the left `acausal' side of the target, thereby restricting comparisons to time points $t \leq T-L^{left-tail}$ near the end of the sample.}. A comparison of expected and sample HTs indicates that the latter are consistently larger, with a statistically significant finding for the MSE nowcast, partially due to the occurrence of protracted recession episodes\footnote{Alternative models employing shorter subsamples partially alleviate this issue without impacting our primary conclusions; thus, they are not included in this discussion.}. The SSA and M-SSA designs exhibit identical expected HTs by design (noting, however, that the underlying models differ). Nonetheless, sample HTs highlight the M-SSA as the smoothest nowcast, characterized by the fewest zero-crossings. This finding is further supported by the corresponding filter outputs displayed in Fig. (ref). The M-SSA demonstrates superior sample performances compared to the univariate designs, both in terms of target correlation and HT. Additionally, the CCF depicted in Fig. (ref) (bottom right panel) shows that M-SSA and HP-C are now coincident, unlike the SSA discussed in the previous section, which displayed a lagging behavior. The figure further highlights a two-month lead of the M-SSA design over the univariate SSA (top right). Finally, the bivariate MSE nowcast outperforms in terms of expected and sample target correlations, noting that the expected value of the SSA in Table (ref) cannot be directly compared, due to differing underlying models. \\

To conclude, it is interesting to examine the filtered components of the INDPRO M-SSA nowcast, which are derived from the sub-filters applied to both INDPRO and CLI. Figure (ref) illustrates these components alongside the aggregate for both the financial crisis and the `great lockdown'. Prior to and following the crises, the dynamics of the CLI component are relatively weak, and the nowcast (in blue) aligns with the INDPRO component (in green). However, during the crises, the CLI (in red) exhibits increased dynamics, leading to a divergence of the aggregate M-SSA from the INDPRO component. Notably, the lead introduced by the CLI results in a leftward shift of the aggregate nowcast in relation to the INDPRO component.

Conclusion

We propose a multivariate extension of the SSA, termed M-SSA, as introduced in Wildi (2024). M-SSA represents an innovative approach to prediction and smoothing, highlighting target correlation and the sign accuracy of the predictor while adhering to a novel HT constraint. The conventional MSE approach corresponds to unconstrained M-SSA optimization, subject to a specific scaling factor. However, our proposed criterion captures more nuanced performance metrics related to accuracy and smoothness (AS) characteristics.\\

In its primal formulation, M-SSA seeks to optimally track the target while enforcing noise suppression; conversely, in its dual formulation, the predictor minimizes zero-crossings for a predetermined level of tracking accuracy. A single hyperparameter controls the AS trade-off, tracing a new AS efficient frontier. Whereas the MSE predictor occupies a single spot on this curve, our approach spans the entire frontier.\\

The M-SSA predictor is both interpretable and attractive due to its inherent simplicity, as it integrates fundamental concepts of prediction, such as sign accuracy, mean squared error, and smoothing requirements. The predictor operates under a simple time-reversible unstable difference equation. Importantly, the stability of the predictor is maintained through the presence of implicit `zero' boundary conditions.\\

The M-SSA framework can effectively address forecasting, smoothing, and signal extraction applications, depending on whether its target specification is causal or non-causal, as well as whether it is all-pass or not. This methodology enables the customization of benchmark predictors in terms of smoothness and accuracy performance. Our applications in forecasting, smoothing, and signal extraction demonstrate the customization of MSE and HP benchmarks. Furthermore, our M-SSA package expands this analysis to incorporate additional contemporary business-cycle tools. Notably, the bivariate signal extraction application highlights the significance of leading indicators within a nowcasting framework.