Implicit score-driven filters for time-varying parameter models

spacing{1.3} { Abstract} We propose an observation-driven modeling framework that allows model parameters to vary over time through an implicit score-driven (ISD) update. The ISD update maximizes the logarithmic observation density with respect to the parameter vector while penalizing the weighted \(\ell_2\) norm relative to a one-step-ahead predicted parameter. This yields an implicit stochastic-gradient update. We show that the popular class of explicit score-driven (ESD) models arises when the observation log density is linearly approximated around the prediction. By preserving the full density, the ISD update extends the favorable local properties of the ESD update to a global setting. For log-concave observation densities, whether correctly specified or not, the ISD filter is stable for all learning rates, and its updates are contractive in mean squared error toward the (pseudo-)true parameter at every time step. We demonstrate the usefulness of ISD filters in simulations and empirical applications in finance and macroeconomics.

bibunit[chicago] \section{Introduction} Empirical evidence shows that assuming constant model parameters over prolonged periods of time is often too restrictive. In economics and finance, parameters are often regime dependent or subject to structural breaks (e.g., stock1996evidence). They may also evolve gradually without an obvious pattern, making it unclear how to update them after observing new data. In some cases, ex-post estimators can be constructed; for example, in ARCH-type models (see terasvirta2009introduction), the squared shock is an unbiased proxy for the true variance. In general, however, such proxies may be difficult to derive, inefficient, or nonexistent. We propose a new framework that makes model parameters time-varying in an observation-driven setting via an implicit score-driven (ISD) filter. Analogous to kalman1960new's (kalman1960new) filter, the ISD filter alternates between prediction and update steps. The update solves an optimization problem that maximizes the current observation's log-likelihood contribution subject to a weighted $\ell_2$ penalty centered at the one-step-ahead prediction. The penalty weights are determined by a positive-definite matrix, whose inverse can be interpreted as a learning-rate matrix. Maximizing the log-likelihood contribution allows new information to be efficiently incorporated, while the penalty regularizes the parameter adjustment. The ISD setup also enables automatic coordination of multiple interacting parameters and the incorporation of parameter constraints without necessitating parameter transformations (e.g., link functions). The first-order condition of the ISD update can be written as an implicit stochastic-gradient step: implicit because the gradient is evaluated at the updated rather than predicted parameter, and stochastic because it uses noisy data. In the optimization literature, implicit updates are known as proximal-point methods and recognized as inherently more stable than their explicit counterparts, which arise as first-order approximations. For example, ISD updates guarantee log-likelihood improvements at each observation---unlike explicit score-driven (ESD) updates (Section (ref)). The ISD filter enjoys several attractive theoretical properties that are typically sought in observation-driven models, but rarely combined in a single framework. First, the ISD filter is invertible under easily verifiable conditions: concavity of the logarithmic postulated density is typically sufficient, even when misspecified (Theorem (ref)). As such, this invertibility result requires no knowledge of the data-generating process (DGP). Hence, any differences due to the initialization vanish almost surely and exponentially fast, a key property in the filtering literature (e.g., bougerol1993kalman; straumann2006quasi). Second, for concave logarithmic densities (again, even when misspecified), the ISD update is contractive in mean squared error toward a small region around the (pseudo-)true parameter (Theorem (ref)). On average, therefore, updates improve on predictions. In fact, the largest gains are expected for the worst predictions. This result requires the existence of a pseudo-true parameter, but is otherwise largely independent of the DGP. Only when the prediction is very close to the pseudo-true parameter may the update be less accurate, as is unavoidable for any stochastic-optimization method (e.g., patrascu2018nonasymptotic). Beyond being largely DGP agnostic, Theorems (ref)--(ref) are also robust in that they hold for an arbitrary (positive-definite) learning-rate matrix. In contrast, ESD filters require additional Lipschitz-type conditions and sufficiently small learning rates to avoid instability or divergence, especially in misspecified settings. Theorems (ref)--(ref) assume concavity of the researcher-postulated log density to obtain strong theoretical guarantees. While this reliance on concavity/convexity is standard in the optimization literature (e.g., boyd2004convex; nesterov2018lectures), throughout the article we highlight extensions indicating that the proposed method is broadly applicable (see Remark (ref)). We illustrate the properties of ISD filters through simulation experiments and empirical applications. In simulations, we show that under correct specification the static (hyper-)parameters are accurately recovered (Section (ref)). Under misspecification, ISD filters reliably track the pseudo-true parameter path even when the ESD filter diverges (Section (ref)), which can arise even in simple settings if the gradient is not Lipschitz continuous. Although this risk of divergence is well known in optimization, it appears to have been overlooked in the score-driven filtering literature (Section (ref)). In a high-dimensional network-flow design, we further demonstrate that, as dimensionality increases, ISD filtering outperforms particle filtering in both estimation accuracy and computational efficiency (Section (ref)). Finally, we present three empirical applications highlighting advantages of ISD over ESD filters. First, we estimate a linear regression of Microsoft returns on market returns with a time-varying beta. Second, we analyze growth-at-risk via lower quantiles of U.S. GDP growth, where the ISD filter yields an implicit version of engle2004caviar's (engle2004caviar) adaptive CAViaR model. Because ISD quantile updates cannot cross the observation on which they are based, simple restrictions on the prediction step ensure that jointly modeled quantiles remain properly ordered. Third, in a study of T-bill rate spreads, the ISD filter outperforms its explicit counterpart even when our concavity assumptions are violated. The paper is structured as follows. Section (ref) introduces ISD filters and compares them with their conventional ESD counterparts. Sections (ref) and (ref) develop our main theoretical results on filter stability and accuracy, respectively. Section (ref) discusses maximum-likelihood estimation of the static parameters, showing accurate parameter recovery in simulations. Section (ref) contains simulations investigating the filtering performance under misspecification. Section (ref) presents the empirical applications, while Section (ref) concludes. The appendix contains proofs of the main results (Appendix (ref)), additional theoretical results (Appendix (ref)), and further details for the simulations and empirical analyses (Appendix (ref)). \subsection{Positioning in the literature} This paper intersects with two strands of literature, summarized along two axes in Table (ref), differentiated by (a) the methodology (explicit vs.\ implicit gradient methods), and (b) the goal (learning vs.\ tracking). To our knowledge, this paper is unique in the fourth quadrant: it uses an implicit gradient method to track a dynamic parameter.\footnote{The overview in Table (ref) is not exhaustive; for instance, we omit simulation-based approaches such as particle filters (e.g., chopin2020introduction).} Along the methodology axis (implicit vs.\ explicit gradient methods), the ISD filter aligns with implicit methods for static optimization, notably rockafellar1976monotone's (rockafellar1976monotone) proximal-point algorithm, which couples a static target with a quadratic penalty around the previous iterate. Because our log likelihood uses random observations from the true density, the ISD filter at each step can be viewed as a stochastic proximal-point method (e.g., bauschke2003bregman; ryu2016stochastic; bianchi2016ergodic; patrascu2018nonasymptotic; asi2019stochastic). Proximal optimization is equivalent to an implicit stochastic-gradient step (e.g., toulis2015scalable; toulis2016towards; toulis2017asymptotic; toulis2021proximal). Our approach also ties in with online-learning methods that process data sequentially (e.g.,\ orabona2019modern; cesa2021online), including machine-learning applications (e.g.,\ kulis2010implicit). What sets our work apart from implicit gradient methods in the (stochastic) optimization literature is that the parameter to be estimated is dynamic. \begin{table}[t!] \caption{Overview of related methods. } \begin{threeparttable} \begin{tabular}{l|ll} & Explicit gradient method & Implicit gradient method \\ \hline Learning & \scriptsize SGD (e.g., robbins1951stochastic; & \scriptsize ISGD (e.g., patrascu2018nonasymptotic; \\ (static target) &\scriptsizeamari1993backpropagation; kushner2010stochastic; bottou2012stochastic)& \scriptsizeasi2019stochastic; toulis2021proximal)\\ \hline Tracking & \scriptsize Constant-gain tracking (e.g., benveniste2012adaptive) & \scriptsize {\textbf{ISD filter}} (this article) \\ (dynamic target) & \scriptsize ESD filter (e.g., creal2013generalized; harvey2013dynamic; \\ &\scriptsize for examples, see \href{https://www.gasmodel.com}{www.gasmodel.com}) \end{tabular} \begin{tablenotes} • Note: (I)SGD = (implicit) stochastic gradient descent. (I/E)SD = (implicit/explicit) score driven. \end{tablenotes} \end{threeparttable} \end{table} Along the goal-related axis (learning vs.\ tracking), our work aligns with engineering and optimization approaches that, following the seminal work of robbins1951stochastic, use explicit gradient methods with constant (i.e., non-vanishing) step sizes to track dynamic targets (e.g., ljung1977analysis; benveniste2003measure). The resulting procedures are usually referred to as constant-gain tracking (for an overview, see benveniste2012adaptive), although extensions with adaptive step sizes have also been considered (e.g., kushner2002analysis). In econometrics, dynamic conditional score (DCS; harvey2013dynamic) models and generalized autoregressive score (GAS; creal2013generalized) models likewise use the (explicit) gradient of the log-likelihood function, known as the \emph{score}, to update time-varying parameters. This framework encompasses many established models (e.g., GARCH) and is popular for its ease of use and strong forecasting performance (e.g., creal2014observation; harvey2014filtering; koopman2016predicting; harvey2017volatility; harvey2018modeling; opschoor2018new; gorgi2020beta). It has been used in ${\sim}400$ published articles; for a near-exhaustive list, see \href{https://www.gasmodel.com}{www.gasmodel.com}. Recent surveys (artemova2022score1,artemova2022score2; harvey2022score) have converged on the terminology of score-driven (SD) models. To align with this nomenclature while distinguishing our approach, we refer to this model class as using \emph{explicit} score-driven (ESD) filters. As this article demonstrates, ESD filters can be obtained within the ISD framework by locally linearizing the logarithmic observation density around the prediction at each time step, an apparently new insight in the econometrics literature. We will show that avoiding this local-linear approximation has both theoretical and practical benefits. \section{Implicit score-driven filters} We consider an $N\times 1$ variable of interest $y_t$, observed at times $t=1,\dots,T$, drawn from a data-generating process (DGP) characterized by a time-varying density, which we denote by $p^0_t(\cdot)$. We refer to $p^0_t(\cdot)$ as a density, though it could also be a mass function (i.e., we allow for discrete observations $y_t$). The dynamics and functional form of the true density $p_t^0(\cdot)$ are left, for the most part, unspecified. It could be parametric, in which case $p^0_t(\cdot)=p^0(\cdot |\theta^0_t)$, where $\theta^0_t$ is a $K_0 \times 1$ time-varying parameter vector taking values in some (non-empty) parameter space $\Theta^0$. Alternatively, $p^0_t(\cdot)$ could be a non-parametric time-varying density. The aim of this paper is to construct a modeling framework that tracks the true density $p^0_t(\cdot)$ using filters that alternate between prediction and update steps. Let $p(\cdot | \theta_{t}, x_{t}, \psi)$ denote the researcher-postulated observation density, which may be misspecified. Here, $\theta_t$ is a $K \times 1$ vector of time-varying parameters taking values in a (non-empty) convex parameter space $\Theta \subseteq \mathbb{R}^K$, $x_t$ is an $\mathcal{F}_{t-1}$-measurable variable, where $\mathcal{F}_{t-1}$ is the information set at time $t-1$, and $\psi$ is a vector of static parameters. The inclusion of $x_{t}$ allows the density to depend on exogenous variables and/or lags of $y_t$. For readability, we write $p(\cdot | \theta_t)$, suppressing the dependence on $x_t$ and $\psi$. We denote the researcher's predicted and updated parameter vectors by $\theta_{t\mid t-1} \in \Theta$ and $\theta_{t\mid t} \in \Theta$, which are based on $\mathcal{F}_{t-1}$ and $\mathcal{F}_t$, respectively. While the postulated observation density $p(\cdot|\theta_t)$ may be misspecified relative to the true density $p_t^0(\cdot)$, model selection remains relevant. This can be illustrated using a trivial case: if $\Theta$ is a singleton, filtering is useless because $\theta_{t|t}$ cannot vary, so the filter is severely misspecified whenever the true density $p^0_t(\cdot)$ is time varying. In practice, we aim to narrow the misspecification gap by postulating a sufficiently flexible density $p(\cdot|\theta_t)$ that can closely approximate $p_t^0(\cdot)$. For our theory development, however, we treat $p(\cdot|\theta_t)$ as generic but fixed. We impose conditions relating only to the postulated density $p(\cdot|\theta_t)$ (i.e., Assumptions (ref)--(ref)), thereby remaining agnostic about the true density $p^0_t(\cdot)$, or minimal conditions relating $p(\cdot|\theta_t)$ to $p^0_t(\cdot)$, such as the existence of a pseudo-true parameter (i.e., Assumptions (ref)--(ref)). Given some postulated density $p(\cdot|\theta_t)$, the main difficulty in designing filtering algorithms lies in specifying how the update $\theta_{t|t}$ should be constructed from the prediction $\theta_{t|t-1}$ and the observation $y_t$. Here we argue that a sound updating scheme should satisfy at least two criteria. First, the update should yield an improved fit of the observed data $y_t$ in terms of the postulated density; i.e., we want $p(y_t|\theta_{t|t}) \geq p(y_t|\theta_{t|t-1})$. As we shall see, explicit score-driven (ESD) filters, despite their popularity (Section (ref)), generally fail to meet this requirement. Second, as each observation $y_t$ is inherently noisy, it is desirable to regularize the extent to which the update $\theta_{t|t}$ deviates from the prediction $\theta_{t|t-1}$. Penalizing the magnitude of $\theta_{t|t}-\theta_{t|t-1}$ prohibits the filtered parameter path from becoming excessively volatile. To satisfy both criteria, we propose the new class of implicit score-driven (ISD) filters. These filters perform the parameter update at time $t$ by maximizing the researcher-postulated logarithmic observation density $\log p(y_t|\cdot)$ subject to a weighted $\ell_2$ penalty centered at the prediction $\theta_{t|t-1}$. That is, the parameter update is defined as \begin{equation} \theta_{t|t} \; :=\; \underset{\theta \in \Theta}{\mathrm{argmax}} \; f(\theta|y_t, \theta_{t|t-1}, P_t), \end{equation} where \begin{equation} f(\theta|y_t, \theta_{t|t-1}, P_t) \; :=\; \log p(y_t|\theta) \, -\, \frac{1}{2}\big\|\theta - \theta_{t|t-1}\big\|_{P_t}^2. \end{equation} Here, $f(\theta|y_t, \theta_{t|t-1}, P_t)$ denotes the “regularized” log-likelihood contribution and $\|z\|_{P_t}^2 = z' P_t z$ is the squared $\ell_2$ norm with respect to a $K\times K$ positive-definite penalty matrix $P_t$. Formulating the parameter update as a maximization problem yields several benefits: (a) it uses the full conditional density, not just moments, (b) the elements of $\theta_{t|t}$ are interdependent because they jointly solve the multivariate problem (ref), and (c) as an optimizer, $\theta_{t|t}$ automatically lies in the correct space $\Theta$ without necessitating link functions. For example, if $\theta_t$ contains positive shape parameters (as in Section (ref)), update (ref) keeps them positive; likewise, the ordering of multiple quantiles (as in Section (ref)) is preserved. More generally, the optimization framework easily accommodates additional constraints. The $\ell_2$ penalty yields tractable updates and can be interpreted as a second-order Taylor expansion around $\theta_{t|t-1}$ of a loss function, where $P_t$ acts as the Hessian. The ISD update (ref)–(ref) mirrors rockafellar1976monotone's (rockafellar1976monotone) classic proximal-point algorithm, which similarly considers the optimization of a target function---in our case, the log-likelihood contribution of the (\emph{a priori} random) observation $y_t$---subject to a quadratic penalty. Hence, for a fixed time step, the approach is a stochastic proximal-point algorithm (e.g., bauschke2003bregman; asi2019stochastic); the difference is that we consider a moving target (Section (ref)). Update (ref) can also be viewed as the posterior mode in a (possibly misspecified) Bayesian setup, with the quadratic penalty acting as a Gaussian prior. This resembles Laplace approximations in both the Bayesian (e.g., rue2009approximate) and frequentist approaches (e.g., koyama2010approximate). Update (ref) reduces to Kalman's level update when $p(\cdot|\theta)$ is Gaussian, the mean is linear in $\theta$, and the penalty matrix $P_t$ is the inverse predicted covariance (lange2024short). For a Gaussian linear regression with time-varying slopes, update (ref) specializes to the normalized least mean squares (NLMS) filter (nagumo1967learning); we return to this case in an empirical illustration (Section (ref)). Although the link between (least-squares) optimization and the Kalman filter has been known since bierman1977factorization and bertsekas1996incremental, it has recently attracted renewed interest in signal processing (e.g., akyildiz2019probabilistic), control (e.g., simonetto2024nonlinear), and econometrics (e.g., lange2024bellman,lange2024short). These links justify investigation of the proximal method (ref) in a more \emph{general} setting: the observation density may be non-Gaussian, $\theta$ need not relate to the mean, and the quadratic penalty need not represent a Gaussian prior. We aim to remain agnostic about the true density sequence $\{p_t^0(\cdot)\}$; i.e., we do not presume latent states with linear Gaussian dynamics. Accordingly, we postulate (ref) as part of a \emph{filter} or \emph{algorithm}, a conceptually distinct approach from imposing conditions on the DGP. We then investigate the algorithm's performance, especially when some or all of the classic assumptions fail. Despite its simplicity and close connection with existing methods, the proposed method is---at this level of generality---new. Whereas the DGP is generally unknown, the postulated density $p(\cdot|\theta)$ is under the researcher's control; typically, it is known in closed form. Assumptions (ref)--(ref) relate only to this postulated density and are, therefore, practically verifiable. Assumptions (ref) and (ref) are standard in optimization, ensuring existence and uniqueness of the maximizer in ((ref)). Assumptions (ref) and (ref) allow its characterization using standard first-order conditions; while not strictly necessary (e.g., subgradients could be used), this aids clarity and tractability. \begin{assumption} \textbf{\emph{(Existence)}} The solution set of $\smash{\underset{\theta \in \Theta}{\mathrm{argmax}} \: f(\theta|y_t,\theta_{t|t-1}, P_t)}$ is non-empty with probability one. \end{assumption} \begin{assumption} \textbf{\emph{(Strictly concave regularized log likelihood)}} $f(\theta|y_t,\theta_{t|t-1}, P_t)$ is proper concave and strictly concave in $\theta$, $\forall \theta \in \Theta$ with probability one. \end{assumption} \begin{assumption} \textbf{\emph{(Interior solution)}} $\theta_{t|t} \in \mathrm{Int}(\Theta)$ with probability one. \end{assumption} \begin{assumption} \textbf{\emph{(Differentiability)}} $\log p(y_t|\theta)$ is at least (a) once or (b) twice continuously differentiable in $\theta$, $\forall \theta \in \mathrm{Int}(\Theta)$ with probability one. When left unspecified, (b) holds. \end{assumption} Assumptions (ref) and (ref) can typically be satisfied by choosing a sufficiently large penalty $P_t$. In particular, even when the postulated log density is poorly behaved (e.g.,\ non-concave or multimodal), a sufficiently strong penalty ensures that update (ref) remains well-behaved. Under Assumptions (ref) through (ref)a, the first-order condition for the parameter update $\theta_{t|t}$ in the maximization problem ((ref)), i.e.,\ $0=\nabla( y_t | \theta_{t|t})-P_t(\theta_{t|t}-\theta_{t|t-1})$, can be rearranged as \begin{equation} \theta_{t|t}\; =\; \theta_{t|t-1}\, +\, H_t \, \nabla( y_t | \theta_{t|t}), \end{equation} where the inverse penalty $H_t := P_t^{-1}$ is referred to as the learning-rate matrix at time $t$ and $\nabla(y_t|\theta_{t|t}) := (\partial \log p(y_t|\theta) /\partial \theta)|_{\substack{\theta = \theta_{t|t}}}$ is the \emph{score} evaluated at the updated parameter $\theta_{t|t}$. Representation (ref) shows that the ISD update yields a gradient-type update, in which the learning-rate matrix $H_t$ controls the step size. Because the score is evaluated at the update $\theta_{t|t}$ instead of the prediction $\theta_{t|t-1}$, it is an \emph{implicit} method. Indeed, the update $\theta_{t|t}$ appears on both sides of equation (ref); hence, it is not immediately computable. As the update also depends on the (\emph{a priori} random) observation $y_t$, our method is also closely related to implicit stochastic-gradient methods (Section (ref)). While ISD update ((ref)) may not allow a closed-form solution, Assumptions (ref) and (ref)a ensure the global solution to (ref) can be found numerically using standard techniques. Specifically, for a scalar parameter ($K=1$), root-finding (e.g., bisection) applies. More generally, quasi-Newton methods with line search can be used, with projected or interior-point variants to handle constraints. For a toolbox of proximal schemes, see parikh2014proximal. In the optimization literature, $H_t$ is often set to decrease over time (e.g., $H_t=\mathcal{O}(t^{-1})$), so the estimated parameter converges to a constant vector. Here, we aim to track a moving target, so the parameter path must remain responsive, even asymptotically. We may therefore keep the learning rate constant, as is standard in tracking applications (e.g., kushner2010stochastic; benveniste2012adaptive). Indeed, we typically set $H_t=H$ for all $t$, where $H$ may contain static parameters to be estimated (Section (ref)). The ISD filter can thus be viewed as an implicit stochastic gradient method with a non-vanishing learning rate. It is more flexible in that $H_t$ can be a time-varying positive-definite matrix (rather than a scalar, as is typical in the optimization literature) and includes a separate prediction step, described below, which is crucial in time-series modeling. The prediction step generates one-step-ahead forecasts of the time-varying parameters. For simplicity, we consider a linear first-order specification as follows: \begin{equation} \theta_{t+1|t} \;=\; \omega\, + \, \Phi\, \theta_{t|t}, \end{equation} where $\omega$ is a $K\times 1$ vector of constants and $\Phi$ is a $K\times K$ autoregressive matrix. Conditions ensuring stable recursions are discussed in the next section. Prediction step (ref) could be generalized to allow for non-linear and/or higher order dynamics if these were found to be relevant for a particular application. In economics and statistics, however, mean reversion is often critical, while no additional information is available at the time of the prediction. For these reasons, a more complicated structure may not yield immediate benefits. Throughout, we require that the prediction step (ref) maps $\Theta$ to itself. For $\Theta=\mathbb{R}^K$, this requirement imposes no restrictions beyond $\omega$ and $\Phi$ having real-valued elements. For bounded parameter spaces, it can often be achieved by restricting $\omega$ and $\Phi$. For a scalar variable in $\Theta=\mathbb{R}_{>0}$ (e.g., a variance), it suffices to take $\omega,\Phi>0$. For a scalar variable in $\Theta=(-1,1)$ (e.g., a correlation), we can impose $\Phi\in[0,1)$ and $\tilde{\omega}:=\omega/(1-\Phi)\in(-1,1)$, so that $ \theta_{t+1|t}=(1-\Phi)\tilde{\omega}+\Phi\,\theta_{t|t}, $ is a convex combination of two elements in the convex set $\Theta=(-1,1)$. The same logic extends to higher dimensions; for example, if $\Theta$ is the positive orthant, then a diagonal matrix $\Phi$ with entries $\Phi_{ii} > 0$ and intercepts $\omega_i > 0$ for $i=1,\ldots,K$ suffices. Finally, when the prediction step is the identity mapping (i.e., $\omega=0_K$, $\Phi=I_K$), $\theta_{t+1|t}\in \Theta$ is immediate. As an alternative to imposing constraints on $\omega$ and $\Phi$, we can take $\Theta=\mathbb{R}^K$ and employ link functions mapping $\mathbb{R}^K$ into an appropriate domain. To illustrate, for $K=1$ we can use the exponential map from $\mathbb{R}$ to $\mathbb{R}_{>0}$ to enforce positivity. Although widely used, such non-linear transformations may complicate the analysis of the filter's theoretical properties. \subsection{Relationship with existing (explicit) score-driven filters} The implicit gradient update (ref) suggests a close link with existing (i.e., explicit) updating schemes. Here we show that linearizing the logarithmic observation density in the ISD optimization problem (ref) produces the familiar ESD update. To explain the relationship, we may approximate ((ref)) via a first-order Taylor expansion at $\theta_{t|t-1}$, i.e., using $\log p(y_t|\theta)\approx \log p(y_t|\theta_{t|t-1})+\langle\theta-\theta_{t|t-1},\nabla(y_t|\theta_{t|t-1})\rangle$, where $\langle x_1,x_2\rangle:=x_1'x_2$. To avoid boundary solutions, we may assume $\Theta=\mathbb{R}^K$. The regularized log-likelihood $f(\cdot|\cdot,\cdot,\cdot)$ in (ref) then becomes a linear objective with a quadratic penalty, so the optimization can be solved in closed form. Indeed, the resulting linearized version of ((ref)) and its associated closed-form solution read \begin{equation} \theta_{t|t}^\textnormal{ex} := \underset{\theta \in \mathbb{R}^K}{\mathrm{argmax}}\left\{ \underbrace{\log p(y_t|\theta_{t|t-1}) + \langle \theta - \theta_{t|t-1},\nabla( y_t | \theta_{t|t-1})\rangle }_{\textnormal{linear approximation of $\log p(y_t|\theta)$ at $\theta_{t|t-1}$}}- \frac{1}{2}\|\theta - \theta_{t|t-1}\|_{P_t}^2\right\}, \end{equation} \begin{equation} \theta_{t|t}^\textnormal{ex} \;=\; \theta_{t|t-1}\, +\, H_t \, \nabla( y_t | \theta_{t|t-1}), \end{equation} where the explicit update is denoted $\theta_{t|t}^\textnormal{ex}$ to differentiate it from the ISD update (ref). Combining the explicit gradient update (ref) with the linear prediction step (ref) reproduces the well-known class of dynamic conditional score (DCS; harvey2013dynamic) models or generalized autoregressive score (GAS; creal2013generalized) models. These equivalent model classes are collectively known as \emph{score driven} (Section (ref)). To align with this standard terminology, while emphasizing the difference with our approach, we refer to this existing model class as using \emph{explicit} score-driven (ESD) filters. Indeed, combining the explicit gradient step (ref) with the linear prediction step (ref) gives \begin{equation} \theta_{t+1|t}^\textnormal{ex} \;=\; \omega \,+\, \Phi\, H_t \,\nabla( y_t | \theta_{t|t-1}) \,+\, \Phi\, \theta_{t|t-1}. \end{equation} In this literature, it is standard to take $H_t = H S_t$, where $H$ is static, while $S_t$ is known as a \emph{scaling} matrix, which is often based on the Fisher information of the postulated density (e.g.,\ artemova2022score1). Subsequently writing $A := \Phi H$ yields the canonical prediction-to-prediction recursion in the ESD literature (e.g., creal2013generalized). While the relationship between explicit and implicit gradient steps is well-known in the optimization literature (e.g.,\ rockafellar1976monotone), it has not yet been recognized in the econometrics literature. The inherent locality of the first-order Taylor expansion suggests that ESD updates may perform satisfactorily only for small (or, technically, infinitesimal) parameter adjustments $\theta_{t|t}^\textnormal{ex}-\theta_{t|t-1}$. Whenever the parameter adjustment is sizable, no likelihood improvement $p(y_t|\theta_{t|t}^\textnormal{ex})>p(y_t|\theta_{t|t-1})$ is guaranteed. Indeed, the likelihood may regularly deteriorate. In contrast, for the ISD update with an arbitrary (positive-definite) learning-rate matrix $H_t$, the typical case $\theta_{t|t}\neq \theta_{t|t-1}$ automatically implies $p(y_t|\theta_{t|t})>p(y_t|\theta_{t|t-1})$. This property is due to the formulation of the ISD update as an optimization problem, where the value of the objective $f(\cdot|y_t,\theta_{t|t-1},P_t)$ at the optimizer $\theta_{t|t}$ must weakly exceed the value at the prediction $\theta_{t|t-1}$. Simple algebra yields the inequality $\log p(y_t|\theta_{t|t})-\log p(y_t|\theta_{t|t-1})\geq \sfrac{1}{2}\|\theta_{t|t}-\theta_{t|t-1}\|^2_{P_t}$. Hence, the ISD update cannot lead to a worse fit, and in fact produces a strictly better one whenever $\theta_{t|t} \neq \theta_{t|t-1}$. By the same inequality, $\|\theta_{t|t}-\theta_{t|t-1}\|^2_{P_t}$ is bounded if the likelihood improvement is bounded, which can be viewed as a robustness property of ISD updates. Generally, implicit and explicit gradients can point in opposite directions. When the maximization problem (ref) is concave, however, both recommend adjustments in roughly the same direction; i.e., the angle between $\theta_{t|t}-\theta_{t|t-1}$ and $\nabla(y_t|\theta_{t|t-1})$ cannot exceed $90^\circ$. \sloppy \begin{proposition} \textbf{\emph{(Relationship between ISD and ESD updates)}} Fix $t>0$ and let Assumptions (ref), (ref) and (ref)a hold. Consider a prediction $\theta_{t|t-1} \in \Theta$ and positive-definite penalty $P_t\in \mathbb{R}^{K \times K}$. Compute $\theta_{t|t}$ using the ISD update (ref). Then, with probability one, \begin{equation} \langle \theta_{t|t} - \theta_{t|t-1}, \nabla( y_t | \theta_{t|t-1}) \rangle \; \geq\; 0. \end{equation} If Assumptions (ref) and (ref)b also hold, we may write: \begin{equation} \theta_{t|t} = \theta_{t|t-1} + (P_t + \mathcal{I}_{t|t})^{-1}\nabla(y_t|\theta_{t|t-1}), \end{equation} where $\mathcal{I}_{t|t}$ denotes the negative average $K \times K$ Hessian between $\theta_{t|t-1}$ and $\theta_{t|t}$, \begin{equation} \mathcal{I}_{t|t} := -\int_{0}^{1} \left.\frac{\partial^2 \log p(y_t|\theta)}{\partial\theta \partial\theta'}\right|_{\substack{\theta \,= \, u \,\theta_{t|t-1} \,+\, (1-u)\,\theta_{t|t}}} \mathrm{d}u. \end{equation} \end{proposition} For a scalar time-varying parameter ($K=1$), equation (ref) implies that the implicit and explicit adjustments ($\theta_{t|t}-\theta_{t|t-1}$ and $\theta^\textnormal{ex}_{t|t}-\theta_{t|t-1}$) have the same sign; naturally, so do the gradients $\nabla(y_t|\theta_{t|t})$ and $\nabla(y_t|\theta_{t|t-1})$. For the ISD update, the derivative of $\log p(y_t|\theta)$ evaluated at the update thus has the same sign as the derivative at the prediction. Because the derivative cannot switch signs, the ISD update increases the value of $\log p(y_t|\theta)$ without surpassing the peak. For the ESD update, the derivative at the update and prediction can have opposite signs; hence, the explicit update can overshoot the peak. For this reason, ESD updates generally fail to ensure $ p(y_t|\theta^\text{ex}_{t|t})>p(y_t|\theta_{t|t-1})$. The second result of Proposition (ref) shows that the ISD update can be written as a “curvature-corrected” version of the ESD update. Here, $\mathcal{I}_{t|t}$ is the average negative $K\times K$ Hessian between $\theta_{t|t-1}$ and $\theta_{t|t}$; i.e., the average curvature of $\log p(y_t|\theta)$ between these points. The ISD update adjusts the step size by accounting for second-order effects. If the log-likelihood is linear in $\theta$, then $\mathcal{I}_{t|t}=O_K$ such that both updates are identical. If the log-likelihood is (multivariate) quadratic (e.g., as with a Gaussian distribution in terms of the mean), $\mathcal{I}_{t|t}$ is constant and both updates are equivalent, albeit for different penalty matrices. To further characterize the ISD update (ref), we focus on log-likelihoods $\log p(y_t|\theta)$ that are concave in $\theta$, yielding particularly strong stability and tracking results. \begin{assumption} \textbf{\emph{(Log-concave observation density)}} $\log p(y_t| \theta) + \frac{\alpha_t}{2}\,\|\theta\|^2$ is concave in $\theta$ for some $\alpha_t \geq 0$, $\forall \theta \in \Theta$, with probability one. \end{assumption} Assumption (ref) strengthens Assumption (ref) by imposing concavity on the log-likelihood function itself, rather than on its regularized version (ref). Concavity strength is indexed by $\alpha_t\ge0$, where $\alpha_t=0$ implies concavity, while $\alpha_t>0$ implies $\alpha_t$-strong concavity. Assumption (ref) allows us to show that, under a common penalty matrix $P_t$, the implicit gradient update is a “shrunken” version of the explicit gradient update. \begin{proposition} \textbf{\emph{(Step-size shrinkage)}} Fix $t>0$ and let Assumptions (ref) to (ref) hold. Take a prediction $\theta_{t|t-1} \in \Theta$ and positive-definite penalty $P_t\in \mathbb{R}^{K \times K}$ as given. Based on the observation $y_t$, compute $\theta_{t|t}$ using the ISD update (ref) and $\smash{\theta_{t|t}^\textnormal{ex}}$ using the ESD update (ref). Let $\lambda_{\max}(P_t)$ denote the largest eigenvalue of $P_t$. Then, with probability one, \begin{equation} \big\|\theta_{t|t}-\theta_{t|t-1}\big\|^2_{P_t} \; \leq \; \underbrace{\left(\frac{\lambda_{\max}(P_t)}{\lambda_{\max}(P_t)+ \alpha_t} \right)^2}_\textnormal{$\in [0,1]$, contraction coefficient}\big\|\theta^\textnormal{ex}_{t|t}-\theta_{t|t-1}\big\|^2_{P_t}. \end{equation} \end{proposition} The contraction coefficient in (ref), which returns multiple times in Sections (ref) and (ref), depends on the ratio between the measure of concavity and the penalization strength. Specifically, larger $\alpha_t$ or smaller $\lambda_{\max}(P_t)$ implies more shrinkage. To understand why the explicit step size is larger, note that for concave log-likelihoods, every tangent lies above the curve. Hence, the linear approximation underlying the explicit update \emph{over}states the achievable likelihood gain, implying that the larger step size of the explicit update is suboptimal. In practice, the shrinkage property in (ref) lets ISD filters use larger learning rates without sacrificing stability. As we will see, the update map is contractive for any positive-definite learning-rate matrix (Lemma (ref)). By contrast, ESD filters need additional (Lipschitz) conditions on the gradient and smaller learning rates to avoid overshooting (and possibly divergence), which can limit their responsiveness. Our results thus mirror well-known optimization results: implicit gradient methods typically can be more responsive without compromising stability (e.g., toulis2017asymptotic). \begin{remark}[Beyond log-concavity] Assumption (ref) is maintained mainly because it yields our strongest theoretical guarantees, in the same spirit that convexity and concavity assumptions are central in optimization theory (e.g., boyd2004convex; nesterov2018lectures). This should not be taken to mean that the proposed method is only useful in such settings: the optimization literature indicates that implicit gradient methods can remain effective for non-concave objectives (e.g., hare2009computing; grimmer2023landscape). To emphasize general applicability, we highlight extensions beyond the concave case throughout the paper. In Sections (ref) and (ref), we develop additional theoretical results on stability and mean squared error (MSE) improvements in non-concave settings. In Sections (ref) and (ref), we consider non-concave log densities in simulation studies and show accurate static \mbox{(hyper-)}parameter recovery and filtering performance. Finally, in an empirical illustration in Section (ref), we show that the implicit filter outperforms competing methods even when both Assumptions (ref) and (ref) fail. \end{remark} \section{Filter stability} We investigate the stability properties of the proposed ISD framework, providing sufficient conditions for filter invertibility, meaning that filtered paths based on identical data but with different initializations converge exponentially fast over time. We remain agnostic with regard to the DGP and use Assumptions (ref)--(ref), which relate to the postulated density only. Our results in this section are presented in three parts: (a) fixing $t$ and examining the update step (Lemma (ref)), (b) fixing $t$ and considering both the update and prediction steps (Lemma (ref)), and (c) proving invertibility by considering the composition of all prediction-to-prediction mappings (Theorem (ref)). We begin by evaluating update stability for a single, fixed time point $t$. Lemma (ref) shows that the ISD update (ref) is stable under Assumptions (ref) through (ref), while, absent further conditions, the same does not hold for the ESD update (ref). \begin{lemma} \textbf{\emph{(Prediction-to-update stability)}} Fix $t>0$ and let Assumptions (ref) to (ref) hold. Let $\theta_{t|t-1}$ and $ \tilde{\theta}_{t|t-1}$ denote two predictions in $\Theta$, which are combined with the observation $y_t$ in the ISD update step ((ref)) to yield corresponding parameter updates, $\theta_{t|t}$ and $\tilde{\theta}_{t|t}$. Then, with probability one, \begin{equation} \big\|\theta_{t|t} - \tilde{\theta}_{t|t}\big\|^2_{P_t} \; \leq\; \underbrace{\left(\frac{\lambda_{\max}(P_t)}{\lambda_{\max}(P_t)+ \alpha_t} \right)^2}_\textnormal{$\in [0,1]$, contraction coefficient}\big\|\theta_{t|t-1} - \tilde{\theta}_{t|t-1}\big\|^2_{P_t}, \end{equation} where $\lambda_{\max}(P_t)$ is the largest eigenvalue of $P_t$. For the ESD update ((ref)), under the additional assumptions that $\nabla(y_t|\theta)$ is $L_t$-Lipschitz continuous in $\theta$ with probability one and $\lambda_{\min}(P_t) \geq L_t/2$, where $\lambda_{\min}(P_t)$ is the smallest eigenvalue of $P_t$, with probability one, \begin{equation} \big\|\theta^\textnormal{ex}_{t|t} - \tilde{\theta}^\textnormal{ex}_{t|t}\big\|^2_{P_t} \leq \underbrace{\frac{\lambda_{\max}(P_t) - \alpha_t[2 - L_t/\lambda_{\min}(P_t)]}{\lambda_{\max}(P_t)}}_\textnormal{$\in [0,1]$, contraction coefficient} \big\|\theta_{t|t-1} - \tilde{\theta}_{t|t-1}\big\|^2_{P_t}. \end{equation} \end{lemma} The first part of Lemma (ref) shows that the ISD update is non-expansive in the squared $P_t$-norm: it never increases (and possibly shrinks) the distance between any two paths. The contraction coefficient matches that in Proposition (ref). With a strongly concave log-likelihood (i.e., $\alpha_t>0$), we obtain a strict contraction in the $P_t$-weighted norm whenever $\theta_{t|t-1}\neq\tilde{\theta}_{t|t-1}$. The second part shows an analogous result for the ESD update ((ref)) under two additional conditions: the score is $L_t$-Lipschitz and $\lambda_{\min}(P_t)\geq L_t/2$. Equivalently, eigenvalues of the learning-rate matrix $H_t$ do not exceed $2/L_t$. This general requirement for explicit gradient methods is well-known in optimization (see boyd2004convex, Eq. 9.17; nesterov2018lectures, Eqns. 1.2.18--22) and machine learning (e.g., wu2023implicit). To prove ESD filter stability under a non-Lipschitz gradient (i.e., $L_t=\infty$), we thus require a zero learning rate. If the DGP is unknown, this seems to be a near-necessary condition: in simulations (Section (ref)), ESD filters can be divergent when the learning rate is positive but the score non-Lipschitz. Below, we focus on proving ISD filter stability; parallel results for ESD filters under additional (Lipschitz) conditions are possible, but not pursued here. We now turn to the prediction-to-prediction mapping from time step $t$ to $t+1$. To obtain a strictly contractive prediction-to-prediction mapping for the ISD filter, it is sufficient for the update and prediction steps to be non-expansive in the ${P_t}$-weighted norm, provided at least one of them is strictly contractive. That is, when $\alpha_t = 0$, the prediction mapping from $\theta_{t|t}$ to $\theta_{t+1|t}$ must be strictly contractive. When $\alpha_t > 0$, on the other hand, it is sufficient for the prediction step to be non-expansive (e.g., allowing $\Phi=I_K$ and $\omega=0_K$). A sufficient condition for non-expansiveness (contractiveness) of the prediction step in the ${P_t}$-weighted norm is that $P_t \succeq \Phi' P_t \Phi$ ($P_t \succ \Phi' P_t \Phi$). Here, the notation $X \succeq Y$ ($X \succ Y$) indicates that $X - Y$ has non-negative (strictly positive) eigenvalues for two symmetric real-valued matrices $X$ and $Y$ of the same size. This requirement is equivalent to $\|\Phi\|_{P_t} \leq 1$ ($\|\Phi\|_{P_t} < 1$), where $\|X\|_{P_t}$ is the induced operator norm of a matrix $X \in \mathbb{R}^{K\times K}$, which is also closely related to the discrete Lyapunov equation (e.g., anderson2012optimal). \begin{lemma} \textbf{\emph{(Prediction-to-prediction stability)}} Fix $t>0$ and let Assumptions (ref) to (ref) hold. Let $P_t$ be given with $P_t \succeq \Phi' P_t \Phi$. Let $\theta_{t|t-1}$ and $ \tilde{\theta}_{t|t-1}$ denote two predictions in $\Theta$ that are used in the ISD update step ((ref)) to yield the corresponding parameter updates $\theta_{t|t}$ and $\tilde{\theta}_{t|t}$, which are subsequently passed on to the prediction step ((ref)) to yield predictions $\theta_{t+1|t}$ and $\tilde{\theta}_{t+1|t}$. With probability one, \begin{equation} \big \|\theta_{t+1|t} - \tilde{\theta}_{t+1|t}\big\|_{P_t}^2 \;\leq\; \kappa_t\, \big\|\theta_{t|t-1} - \tilde{\theta}_{t|t-1}\big\|_{P_t}^2, \end{equation} where the contraction coefficient $\kappa_t$ is \begin{equation} \kappa_t \;=\; \frac{\lambda_{\max}(P_t)[\lambda_{\max}(P_t) -\lambda_{\min}(P_t -\Phi'P_t\Phi)]}{(\lambda_{\max}(P_t) + \alpha_t)^2}\; \in\; [0,1]. \end{equation} Hence the prediction-to-prediction mapping is non-expansive. If either $\alpha_t > 0$ or $P_t \succ \Phi' P_t \Phi$, then, with probability one, $\kappa_t \in [0,1)$; as such, the mapping is contractive. \end{lemma} The contraction coefficient $\kappa_t\in[0,1]$ in the prediction-to-prediction mapping at time $t$ depends on the strength of concavity $\alpha_t$, the penalty matrix $P_t$, and the autoregressive matrix $\Phi$. It equals the contraction coefficients in Proposition (ref) and Lemma (ref) whenever $\lambda_{\min}(P_t-\Phi' P_t\Phi)=0$, as would be the case if $\Phi$ was the identity matrix. Whenever $\lambda_{\min}(P_t-\Phi' P_t\Phi)>0$, the contraction is further strengthened by the autoregressive matrix $\Phi$ in the prediction step, yielding $\kappa_t <1$ even if $\alpha_t=0$. That is, a strict contraction in the prediction-to-prediction mapping can be obtained through either $\alpha_t>0$ or $\lambda_{\min}(P_t-\Phi' P_t\Phi)>0$. If $\alpha_t=0$, it is sufficient that $\lambda_{\min}(P_t-\Phi'P_t\Phi)>0$. For a scalar time-varying parameter, this holds if $|\Phi|<1$. In the multi-parameter case, it holds under $\Phi'\Phi\prec I_K$ if (a) $\Phi$ and $P_t$ are both diagonal or (b) either $\Phi$ or $P_t$ is a scalar multiple of the identity matrix. If $\Phi$ is symmetric, the condition $\Phi'\Phi\prec I_K$ is in turn equivalent to $\varrho(\Phi)<1$, where $\varrho(\cdot)$ is the spectral radius. If $\alpha_t>0$, the above strict inequalities can be made weak. We can also allow for more general matrices $\Phi$ and $P_t$ by imposing the latter to solve the discrete Lyapunov equation $P_t-\Phi'P_t\Phi=\Delta_t\succ0$, which for $\varrho(\Phi)<1$ has a unique solution $P_t\succ0$ parameterized by $\Delta_t\succ0$ (e.g., bof2018lyapunov, Thm. 3.2). As we show in Appendix (ref), ISD filter stability can be established even without log-concavity of the observation density (i.e., without $\alpha_t\geq 0$ as in Assumption (ref)), provided the prediction mapping is sufficiently contractive to offset any expansion in the update step. We illustrate this for a Student's $t$ distribution with time-varying location, which is not log-concave; a more general treatment is developed in donkervanheel2025gradientbased. Finally, we study the composition of all prediction-to-prediction mappings, which must be contractive for the initialization effects to vanish exponentially fast. A sufficient condition is that all mappings are contractive in a single shared norm over time. Theorem (ref) ensures the existence of such a norm and provides an invertibility result that is crucial for maximum-likelihood (ML) estimation of static parameters (e.g., straumann2006quasi). \begin{theorem} \textbf{\emph{(Invertibility)}} For all $t>0$, let Assumptions (ref) to (ref) hold, with either (a) $P_t \succ \Phi' P_t \Phi$ or (b) $P_t \succeq \Phi' P_t \Phi$ and $\alpha_t > 0$. In addition, let there be some $\bar{P}, Q \in \mathbb{R}^{K\times K}$ with $\bar{P} \succ Q \succ O_{K}$ and a sequence $\{\rho_t > 0\}$ such that for all $t>0$, with probability one, \begin{equation} \kappa_t P_t \,+ \, \rho_t Q \;\preceq\; \rho_t \,\bar{P}\; \preceq \; P_t, \end{equation} where $\kappa_t$ is defined in ((ref)). Take two initial values $\theta_{0|0} \in \Theta$ and $\tilde \theta_{0|0} \in \Theta$, yielding two sequences $\{\theta_{t|t-1}\}$ and $\{\tilde \theta_{t|t-1}\}$, respectively. Then the ISD filter composed of ((ref)) and ((ref)) is \emph{invertible}, i.e., there exists a constant $c>1$ such that with probability one, \begin{equation} \lim_{t \to \infty} c^t\big\|\theta_{t|t-1} - \tilde \theta_{t|t-1}\big\|^2\; \rightarrow \; 0. \end{equation} \end{theorem} Theorem (ref) gives a sufficient condition for contraction of all prediction-to-prediction mappings in the common ${\bar P}$-weighted norm, where the (time-invariant) matrix $\bar P$ satisfies condition (ref). For a scalar time-varying parameter, (ref) holds for any sequence $\{P_t\}$ whenever $|\Phi|<1$. For the unit-root case $|\Phi|=1$, it suffices that $\{P_t\}$ is uniformly upper bounded, while ${\alpha_t}$ is uniformly bounded away from zero, thereby preventing $\kappa_t$ from approaching unity. In the multi-parameter case, ((ref)) prevents elements of ${P_t}$ from having drastically different dynamics. The condition holds automatically when their time variation is identical; for example, when $P_t=\zeta_t P$ with $\zeta_t>0$ and $P\succ O_{K}$ satisfying $P\succ \Phi'P\Phi$. In general, the presence of $Q$ in ((ref)) ensures the contraction coefficient remains strictly below unity. Interestingly, Theorem (ref) relies on the researcher-postulated density $p(\cdot|\theta_t)$, but not the true observation density $p^0_t(\cdot)$. Invertibility of the ISD filter can thus be guaranteed without restrictions on the DGP. Because the ISD-filtered path is asymptotically insensitive to initialization, it is likewise asymptotically insensitive to any individual data point. Hence ((ref)) implies that any two parameter paths $\{\theta_{t|t-1}\}$ and $\{\tilde\theta_{t|t-1}\}$ based on different initializations but (almost) the same data are exponentially almost surely (e.a.s.) convergent. As we show in Appendix (ref), Theorem (ref) also prevents numerical errors during implementation in practice from accumulating.\footnote{We thank two anonymous reviewers for suggesting the additional analyses in Appendices (ref) and (ref).} A longstanding concern among researchers is that numerical errors can induce instability, even for closed-form recursions like Kalman's filter (e.g., anderson2012optimal). For the ISD filter, Theorem (ref) allows us to demonstrate that differences stemming from (a) numerical (rounding) errors or (b) inexact ISD update steps remain uniformly bounded over time; i.e., paths cannot diverge. As we show in Appendix (ref), Theorem (ref) further implies the exponential \emph{Lyapunov} stability of the ISD filter: the limiting path $\{\theta_{t|t-1}\}$ equals the Lyapunov equilibrium. Moreover, if $\{y_t\}$ is stationary and ergodic and each $\theta_{t|t}$ is jointly measurable in $(y_t,\theta_{t\mid t-1})$, Theorem (ref) implies that the limiting path is also stationary and ergodic (krengel2011ergodic; brandt1986stochastic). This is critical to ensure the consistency of the ML estimator (Section (ref)). Finally, we can derive an invertibility result for ESD filters analogous to that in Theorem (ref) for ISD filters, but only under the additional Lipschitz and learning-rate conditions given in Lemma (ref). Doing so would yield a new multivariate, DGP-agnostic invertibility result for ESD filters, which contrasts with other approaches (e.g., blasques2022maximum) that rely on DGP-dependent expectations. Our simulations (Section (ref)) confirm the necessity of these additional conditions for ESD filters: if they fail, no positive learning rate may exist that guarantees the stability of the ESD-filtered path. In optimization, toulis2017asymptotic similarly show that explicit gradient methods can be divergent. Stability is even harder to achieve in a dynamic setting, as our filters must remain responsive rather than converge. Any filter with a positive probability of divergence will eventually diverge. Theorem (ref) guarantees that ISD filters based on log-concave densities are stable under DGP-agnostic conditions depending only on $\Phi$, $\{\alpha_t\}$, and $\{P_t\}$. \section{Filter accuracy} Ensuring filter stability, as in Theorem (ref), is necessary but insufficient for high accuracy. After all, even a trivial filter that yields the same output for all inputs is stable. Therefore, accuracy cannot be assessed solely from the internal characteristics of the filter; instead, the analysis must also take account of the true process. We consider a misspecified setting (e.g., the true observation density is unknown) and ask whether our updates still improve the quality of the filtered parameter path. In this case, we can only hope to track, as accurately as possible, the \emph{pseudo}-true parameter (e.g., beutner2023consistency). Assumptions (ref)--(ref) below thus link the postulated density $p(\cdot|\theta_t)$ to the true density $p^0_t(\cdot)$. Under these assumptions, we will show that the ISD update is contractive in mean squared error (MSE) toward a small, “noise-dominated” region of the pseudo-true parameter; moreover, under strong concavity, this contraction is geometric. This result distinguishes filters that are merely stable from those that are also accurate; for example, for the trivially stable filter mentioned above, such improvements are impossible. For the ESD update, an analogous MSE improvement guarantee can be shown to hold, but only if the gradient is Lipschitz continuous and the learning rate sufficiently small. If concavity of the logarithmic density fails, both updates may yield only a \emph{constant} (rather than proportional) MSE reduction. \begin{assumption} \textbf{\emph{(Existence and uniqueness of pseudo-true parameter)}} There exists a $\theta^\star_t \in \Theta$ such that $\smash{\underset{y_t}{\mathbb{E}}[\log p(y_t|\theta^\star_t)] > }$ $\smash{ \underset{y_t}{\mathbb{E}}[\log p(y_t|\theta)]}$ for all $\theta \in \Theta \setminus \{\theta^\star_t\}$. \end{assumption} \begin{assumption} \textbf{\emph{(Score moments at pseudo-true parameter)}} For each $t$, $\underset{y_t}{\mathbb{E}}[\nabla(y_t|\theta^\star_t)] = 0$ and $\underset{y_t}{\mathbb{E}}[\|\nabla(y_t|\theta^\star_t)\|^2] < \infty$. \end{assumption} Assumption (ref) posits the existence of a unique pseudo-true parameter $\theta_t^\star$ that maximizes the expected (postulated) log likelihood, $\mathbb{E}_{y_t}[\log p(y_t|\theta)]$. If the postulated density is differentiable and strongly concave with probability one (i.e., Assumptions (ref)a and (ref) hold with $\alpha_t>0$), Assumption (ref) is automatically satisfied. If $p(\cdot|\theta_t)$ coincides with $p^0_t(\cdot)$, then $\theta_t=\theta_t^\star$. Assumption (ref) thus underscores the need to choose a good candidate $p(\cdot|\theta_t)$, since no filter (including the ISD filter) can eliminate the misspecification gap. Assumption (ref) requires the first moment of the score at $\theta_t^\star$ to be zero and its second moment to be finite. If $\theta_{t|t-1}$ deviates substantially from $\theta_t^\star$, the update $\theta_{t| t}$ based on $y_t$ often represents an improvement. As $\theta_{t| t-1}$ approaches $\theta_t^\star$, however, further improvements become increasingly hard to achieve. If $\theta_{t\mid t-1}$ is already very accurate, the noisy observation $y_t$ may in fact pull the update $\theta_{t| t}$ away from the pseudo-truth. If $\theta_{t| t-1}=\theta_t^\star$, a deterioration is almost inevitable; this is inherent to stochastic optimization methods, not a limitation of our approach. The neighborhood of $\theta_t^\star$ is therefore called the \emph{noise-dominated region} (NDR; e.g.,\ ryu2016stochastic; patrascu2018nonasymptotic; lange2024bellman). Since improvements are not always guaranteed, Theorem (ref) explicitly characterizes the tug of war between contractive and expansive forces, which respectively decrease and increase the MSE relative to the pseudo-true parameter. Which of the two forces dominates largely depends on the accuracy of the prediction. While most authors (e.g., asi2019stochastic) establish upper bounds on the MSE after updating, our equations (ref)--(ref) below provide exact equalities (rather than inequalities). This approach enables us to identify the precise conditions under which updates lead to improvement (see further discussion below). \begin{theorem} \textbf{\emph{(Contraction to the NDR)}} Fix $t>0$ and let Assumptions (ref) to (ref) hold. Then, for the ISD update (ref), we have \begin{equation} \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\theta_{t|t} - \theta^\star_t\big\|_{P_t}^2 \right]}_\textnormal{MSE after update} = \underbrace{\big\|\theta_{t|t-1} - \theta^\star_t\big\|_{P_t}^2}_\textnormal{SE before update} -\underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\theta_{t|t} - \theta^\star_t\big\|_{2\mathcal{I}^{\star}_{t|t} + \mathcal{I}^{\star}_{t|t}P_t^{-1}\mathcal{I}^{\star}_{t|t}}^2 \right]}_\textnormal{$\geq 0$, contractive force} +\underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\nabla(y_t|\theta^\star_{t})\big\|_{P^{-1}_t}^2\right]}_\textnormal{$\geq 0$, expansive force}. \end{equation} For the ESD update (ref), we have \begin{equation} \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\theta_{t|t}^\textnormal{ex} - \theta^\star_t\big\|_{P_t}^2\right]}_\textnormal{MSE after update} = \underbrace{\big\|\theta_{t|t-1} - \theta^\star_t\big\|_{P_t}^2}_\textnormal{SE before update} - \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\theta_{t|t-1} - \theta^{\star}_{t}\big\|^2_{2 \mathcal{I}^{\star}_{t|t-1} }\right]}_\textnormal{$\geq 0$, contractive force} + \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\nabla(y_t|\theta_{t|t-1})\big\|^2_{P_t^{-1}}\right]}_\textnormal{$\geq 0$, expansive force}. \end{equation} Here $\mathcal{I}^{\star}_{t|t}$, $\mathcal{I}^{\star}_{t|t-1} \succeq \alpha_t I_K\succeq O_K$ denote the negative average $K \times K$ Hessians between $\theta_{t|t}$ or $\theta_{t|t-1}$ and $\theta_{t}^{\star}$, that is, \begin{equation} \mathcal{I}^{\star}_{t|t} := -\int_{0}^{1} \left.\frac{\partial^2 \log p(y_t|\theta)}{\partial\theta \partial\theta'}\right|_{\substack{\theta \,= \, u \,\theta_{t|t} \,+\, (1-u)\,\theta_t^{\star}}} \mathrm{d}u, \end{equation} \begin{equation} \mathcal{I}^{\star}_{t|t-1} := -\int_{0}^{1} \left.\frac{\partial^2 \log p(y_t|\theta)}{\partial\theta \partial\theta'}\right|_{\substack{\theta \,= \, u \,\theta_{t|t-1} \,+\, (1-u)\,\theta_t^{\star}}} \mathrm{d}u. \end{equation} \end{theorem} To the best of our knowledge, Theorem (ref) is new in the stochastic optimization literature: we did not find equations (ref)--(ref) in parikh2014proximal, polson2015proximal, ryu2016stochastic, bianchi2016ergodic, or asi2019stochastic. The comparison with this literature is relevant (see Section (ref)), because Theorem (ref) treats a single time step. If the second term on the right-hand side of (ref) is dropped and the equality replaced by an inequality, we obtain a result similar to Theorem 3.2 of asi2019stochastic. Since the expansive force in (ref) is bounded by Assumption (ref), asi2019stochastic conclude that implicit updates are “nondivergent,” whereas explicit updates lack this guarantee. Our contribution in Theorem (ref) is the inclusion of contractive forces, which are critical in improving updates over predictions; these terms allow us to write equalities rather than inequalities. Because updates should ideally be more accurate than predictions (at least outside the NDR), these contractive forces are key. We now analyze equations (ref)--(ref) in detail. On the left-hand side, we have the MSEs of the ISD and ESD updates in the squared $P_t$-weighted norm. Each right-hand side contains three terms. The first term, which is identical for ISD and ESD updates, is the $P_t$-weighted squared error (SE) before updating; i.e., the distance of the prediction from the pseudo-truth. The second term is the contractive force, which differs slightly between the two updates. For the ISD update, it is the \emph{post}-update MSE with weight matrix \(2\mathcal{I}^{\star}_{t|t} + \mathcal{I}^{\star}_{t|t}P_t^{-1}\mathcal{I}^{\star}_{t|t}\); for the ESD update, it is the \emph{pre}-update SE with weight matrix \(2\mathcal{I}^{\star}_{t|t-1}\). These forces are proportional to the negative average Hessians between \(\theta_{t|t}\) or \(\theta_{t|t-1}\), respectively, and \(\theta_t^{\star}\). Under strong concavity (\(\alpha_t>0\) in Assumption (ref)), we have \(\mathcal{I}^\star_{t|t},\mathcal{I}^\star_{t|t-1}\succeq \alpha_t I_K\), so the contractive forces grow quadratically with the distance of \(\theta_{t|t}\) (ISD) or \(\theta_{t|t-1}\) (ESD) from \(\theta_t^\star\). While Assumption (ref) guarantees \(\mathcal{I}^\star_{t|t},\mathcal{I}^\star_{t|t-1}\succeq O_{K}\) with probability one, it could be weakened by introducing a notion of “expected strong concavity,” as in optimization (e.g., toulis2021proximal). The third term on the right-hand sides of (ref)--(ref), the expansive force, reveals a key difference between the two updating methods. For the ISD update, it is the expectation of the $P_t^{-1}$-weighted squared norm of the score evaluated at the pseudo-true parameter $\theta^\star_t$, yielding a noise term that is uniformly bounded (Assumption (ref)) and can be made arbitrarily small by increasing $P_t$. For the ESD update, the same term is evaluated at the prediction $\theta_{t|t-1}$, so its magnitude typically increases with the distance of $\theta_{t|t-1}$ from $\theta_t^\star$. Under strong concavity (i.e., $\alpha_t>0$ in Assumption (ref)), the expansive force in the ESD update grows quadratically with this distance and, hence, may dominate the contractive force. Thus, explicit updates “can be unstable even for relatively simple problems” asi2019stochastic. A key takeaway of Theorem (ref) is that, absent further conditions, ESD updates are not universally beneficial. If the prediction is inaccurate and the expansive force dominates the contractive force, ESD filters may diverge (as we often observe in simulations). This is precluded for the ISD update by the fact that the expansive force is uniformly bounded. Moreover, under strong concavity (i.e., $\alpha_t>0$), the ISD update is geometrically contractive towards the NDR. \begin{corollary} \textbf{\emph{(Geometric contraction to the NDR)}} Fix $t>0$ and let Assumptions (ref) to (ref) hold, where Assumption (ref) holds for some $\alpha_t > 0$. Then the ISD update (ref) satisfies \begin{equation} \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\theta_{t|t} - \theta^\star_t\big\|_{P_t}^2\right]}_{\textnormal{MSE after update}} \leq \underbrace{\left(\frac{\lambda_{\max}(P_t)}{\lambda_{\max}(P_t)+ \alpha_t}\right)^2}_\textnormal{$\in [0,1)$, contraction coefficient} \left(\underbrace{\big\|\theta_{t|t-1} - \theta^\star_t\big\|_{P_t}^2}_{\textnormal{SE before update}} + \underbrace{\underset{y_t}{\mathbb{E}}\left[\big\|\nabla(y_t|\theta^\star_{t})\big\|_{P^{-1}_t}^2\right]}_{\textnormal{(bounded) irreducible noise}}\right). \end{equation} \end{corollary} The contraction in Corollary (ref) is geometric in the sense that, for large prediction errors, the ratio of the \emph{post}-update MSE to the \emph{pre}-update SE is equal to a constant, which is less than unity. This contraction coefficient, which also appears in Proposition (ref) and Lemma (ref), is stronger than comparable results in the literature (e.g.,\ lange2024bellman). To simplify, let $P_t=P I_K$ with $P\in \mathbb{R}_{>0}$, such that the $P_t$-weighted squared error becomes $P$ times the usual (i.e., Euclidean) squared error. Further, denote the (mean) squared error before and after updating as $\text{SE}_{t|t-1}:=\|\theta_{t|t-1}-\theta_t^\star\|^2$ and $\text{MSE}_{t|t}:=\mathbb{E}_{y_t}\|\theta_{t|t}-\theta_t^\star\|^2$, respectively. If $\alpha_t\geq \alpha>0$ for all $t$ in Assumption (ref), then inequality (ref) can be written as \begin{equation} \text{MSE}_{t|t}\;\leq\; \left(\frac{P}{P+\alpha} \right)^2 \,\left(\text{SE}_{t|t-1}+\frac{\sigma^2}{P^2}\right), \end{equation} where $\sigma^2\geq \sup_{t} \mathbb{E}_{y_t}[\|\nabla(y_t|\theta^\star_t)\|^2]$ bounds the gradient noise. It is easy to show that $\text{MSE}_{t|t}<\text{SE}_{t|t-1}$ whenever $\text{SE}_{t|t-1}>\sigma^2/(2P\alpha+\alpha^2)$. Thus, updates are helpful in expectation whenever the squared prediction error exceeds $\sigma^2/(2P\alpha+\alpha^2)$. For more accurate predictions, no MSE gain can be guaranteed; hence, $\sigma^2/(2P\alpha+\alpha^2)$ can be interpreted as the squared radius of the NDR, which increases with $\sigma^2$, but decreases with $\alpha$ and $P$. As such, the NDR can be made arbitrarily small by increasing $P$, but only to the detriment of the contraction rate. Figure (ref)(a) illustrates the geometric contraction to the NDR implied by inequality (ref). For strongly concave log densities, the contraction rate $P^2/(P+\alpha)^2<1$ means inaccurate predictions yield large MSE reductions; i.e., the improvement is roughly proportional to the squared prediction error. For non-concave logarithmic densities such as the Student's $t$ location model, as illustrated in Figure (ref)(b), updates may only yield constant rather than proportional MSE improvements (for details, see Appendix (ref)). \begin{figure}[bth] \begin{tabular}{cc} & \\ (a) MSE contraction (ref) with strong concavity & (b) MSE reduction with non-concavity \end{tabular} \caption{Illustration of mean squared error (MSE) improvements.} \end{figure} Returning to the case of log-concave densities as in Corollary (ref), the penalty $P_t$ affects both the contraction speed and the size of the NDR, thereby creating a trade-off. A larger penalty shrinks the NDR, but worsens the contraction rate. With a static penalty matrix $P_t=P \succ O_{K}$ for all $t$, we can have either a small NDR (large penalty) or rapid contraction (small penalty), but not both. In practice, the penalty matrix $P$ that optimizes this trade-off over longer time periods can be estimated via maximum likelihood (see next section). As we show in Appendix (ref), ESD updates are also contractive toward the NDR, but only under the additional conditions that (a) the score is Lipschitz in $\theta$ and (b) the learning rate is capped. These conditions are necessary to prevent the expansive force from dominating the contractive force if $\|\theta_{t|t-1}-\theta_t^\star\|$ is large. Finally, as in lange2024bellman, the geometric contraction of Corollary (ref) can be used to bound the long-run MSE of $\{\theta_{t|t}\}$ relative to $\{\theta_t^\star\}$, but this requires additional assumptions on the dynamics of the pseudo-true parameter. We leave this to future research; see donkervanheel2025gradientbased for related developments. \section{Estimation of static (hyper-)parameters} The ISD filter's static (hyper-)parameters are unknown and must be estimated. For simplicity, here we consider $P_t= P\succ O_{K}$ for all $t$.\footnote{Alternatively, $H_t=P_t^{-1}$ could be based on (powers of) the time-varying conditional information matrix as is standard in the ESD literature (e.g., creal2013generalized).} We must estimate the penalty matrix $P$ in the update step, $\omega$ and $\Phi$ in the prediction step (ref), and the static parameter $\psi\in\Psi\subseteq\mathbb{R}^M$ in the observation density. All static parameters can be jointly estimated by maximum likelihood (ML) using the standard prediction-error decomposition (e.g., creal2013generalized): \begin{equation} \hat{\xi} :=\; \underset{\xi \, \in \, \Xi}{\mathrm{argmax}} \; \sum_{t=1}^{T} \log p(y_t|\theta_{t|t-1}, \psi), \end{equation} where $\xi := [\mathrm{vech}(P)^\prime,\ \omega^\prime,\ \mathrm{vec}(\Phi)^\prime,\ \psi^\prime]^\prime$ stacks all static parameters, and $\mathrm{vec}(\cdot)$ and $\mathrm{vech}(\cdot)$ denote the (half-)vectorization operations. The optimization domain is $\Xi\subset\mathbb{R}^{\frac{3}{2}K(K+1)+M}$ with $P\succ O_K$ and $\psi\in\Psi$. The initialization $\theta_{0\mid 0}$ is treated as fixed but could be included in $\xi$. To ensure ISD filter stability by Theorem (ref), it is convenient to reparameterize $P$ in terms of $\Phi$ and $\Delta$ using the discrete Lyapunov equation $P-\Phi^\prime P\Phi=\Delta\succ O_K$ (as in Section (ref)). This has the advantage that the standard restrictions $\Delta\succ O_K$ and $\varrho(\Phi)<1$ ensure stability. Building on the existing theory for ESD filters (e.g., blasques2022maximum), we conjecture that, under standard regularity conditions, the ML estimator $\hat{\xi}$ is consistent for the pseudo-true parameter $\xi^\star$ and asymptotically normal. This requires that $\xi^\star$ is identified, $\{y_t\}$ is stationary and ergodic with finite moments and near-epoch dependence, while $p(y|\theta,\psi)$ has uniformly bounded derivatives (of a sufficiently low order) such that the law of large numbers and central limit theorem can be applied (blasques2022maximum). A key ingredient in standard proofs of consistency and asymptotic normality is invertibility in the sense of bougerol1993kalman and straumann2006quasi. For ESD filters, verifying the required contraction conditions is challenging because the maximal admissible learning rate depends on the (typically unknown) DGP; see the discussion following Theorem (ref). For ISD filters, by contrast, Theorem (ref) provides a strong invertibility result with no restriction on the learning rate and no required knowledge of the DGP. We therefore expect the asymptotic results in the ESD literature to carry over to our setting; in fact, they may be easier to verify given the ISD filter's DGP-agnostic stability properties. A potential caveat is that Theorem (ref) requires the postulated logarithmic density to be concave (Assumption (ref)). While sufficient, this condition is likely stronger than necessary. Below, we show that the static (hyper-)parameters can also be accurately recovered when Assumption (ref) fails (Section (ref)), while filtering accuracy remains high (Section (ref)). Although a full asymptotic theory is beyond the scope of the present paper, our simulations below support the conjectured consistency and asymptotic normality of the ML estimates. \subsection{Simulation study: Static (hyper-)parameter estimation} Here we show that the ML estimator (ref) can accurately recover the static parameters and present visual support for the conjecture of asymptotic normality around the true values. We use nine time-varying observation densities from koopman2016predicting; for full specifications, see Appendix (ref). We take $y_t\sim p^0(y_t|\theta_t^0)$ with scalar $\theta_t^0\in\mathbb{R}$; all distributions include link functions (e.g., mapping $\theta_t^0$ to $\mathbb{R}_{>0}$ for volatility). As in koopman2016predicting, the postulated density $p(\cdot|\theta_t)$ is correctly specified, meaning it matches the functional form of the true density (i.e., including the link function). However, some densities include additional shape parameters (see Appendix (ref)) that are treated as unknown and must be estimated: the degrees of freedom $\nu^0 = 6$ for both Student's $t$ distributions and the shape parameter in the negative binomial, gamma, and Weibull distributions (set to 4, 1.5, and 1.2, respectively). koopman2016predicting set $\theta_t^0=\theta_{t|t-1}^{\mathrm{ex}}$, making the ESD filter the DGP; here, we set $\theta_t^0=\theta_{t| t-1}$, making the ISD filter the DGP. In all nine DGPs, we set $\omega^0=0$, $\Phi^0=0.97$, and $H^0=0.10$. Assumption (ref) holds for the first seven densities, which are log-concave in the time-varying parameter. The remaining two are not, with Hessians bounded above by \(1/4\) rather than \(0\). To ensure that Assumption (ref) holds, we impose a lower bound of \(1/4\) on the penalty parameter, or equivalently an upper bound of \(4\) on the learning rate, although this restriction was never binding in practice. We solve (ref) using Newton steps with a simple line search. In fact, for five of the nine cases, the ISD update has a closed-form expression in terms of the Lambert \(W\) function; see Appendix (ref). In practice, however, evaluating the Lambert \(W\) function is no faster than our standard Newton-based routine. \begin{table}[t!] \caption{RMSEs of parameter estimates across $1{,}000$ replications.} \begin{threeparttable} \begin{tabular}{@llccccccl@} \toprule \textbf{DGP type} & \textbf{Distribution} & \textbf{Log} & \textbf{Sample size} & \multicolumn{4}{c}{\textbf{Static-parameter RMSEs}} & \\ && \textbf{concavity} & $T (\times 10^3)$ & $\hat{\omega}$ & $\hat{\Phi}$ & $\hat{H}$ & {$\hat{\psi}$}\\ \midrule Count & Poisson & \checkmark & 1 & .0053& .0266& .0249& \\ & && 4 &.0020& .0074& .0105\\ & && 16 &.0010& .0034 &.0051\\ \midrule Count& Negative bin. & \checkmark &1 &.0053& .0253 &.0284 &1.725\\ &&& 4 &.0020& .0081& .0130& .533\\ &&& 16 &.0010 &.0039 &.0064 &.238\\ \midrule Intensity & Exponential& \checkmark & 1& .0073& .0228& .0273& \\ &&& 4& .0025& .0070 &.0124\\ &&& 16& .0012 &.0033 &.0061\\ \midrule Duration &Gamma & \checkmark &1 &.0072 &.0168 &.0220 &.029 \\ &&& 4 &.0028& .0062& .0100 &.016\\ &&& 16 &.0014 &.0030 &.0053 &.007\\ \midrule Duration & Weibull & \checkmark &1& .0064& .0161& .0193 &.064\\ &&& 4 &.0028 &.0063 &.0103& .030\\ &&& 16 &.0013 &.0029 &.0046& .015\\ \midrule Volatility & Gaussian & \checkmark &1& .0114 &.0666 &.0704 \\ &&& 4 &.0024 & .0109 &.0211\\ &&& 16 & .0010 &.0046 &.0101\\ \midrule Volatility & Student's $t$ & \checkmark & 1 & .0132 &.0801& .0840 &10.233\\ &&& 4 &.0023& .0149 &.0250& .605\\ &&& 16 & .0010 & .0059 & .0117& .275 \\ \midrule Dependence & Gaussian & \ding{55} & 1 & .0143 & .1234 & .1412\\ &&& 4 &.0023& .0223 &.0367\\ &&& 16 &.0008 & .0074 &.0156\\ \midrule Dependence &Student's $t$ & \ding{55} & 1 & .0191 & .1497 &.1740 & .630\\ && & 4 & .0031 & .0370 & .0463 & .318\\ &&& 16 & .0009 & .0096 & .0181 & .146\\ \bottomrule \end{tabular} \end{threeparttable} \end{table} Table (ref) reports RMSEs of static (hyper-)parameter estimates for all DGPs and sample sizes $T\in\{10^3,4{\cdot}10^3,16{\cdot}10^3\}$ with $1{,}000$ replications. Static parameters are generally well-recovered; RMSEs roughly halve when the sample size $T$ quadruples from $4{,}000$ to $16{,}000$. Recovery remains adequate for the final two densities, but the asymptotics appear to set in later: RMSEs still fall by more than half even for the largest sample size. This may be due to the non-concavity and/or general difficulty of estimating dynamic-dependence models. Figure (ref) shows the distribution of static-parameter estimation errors to be approximately Gaussian for the dynamic Poisson model based on the largest sample size. The same is true for the other DGPs (see Appendix (ref)), providing visual support for our conjecture of asymptotic normality around the true values. \begin{figure}[tb] \caption{Kernel estimate (solid) of estimation errors of the ISD filter's static parameters for the dynamic Poisson distribution with $T=16{,}000$. } \end{figure} \section{Simulation studies: Filtering performance} We perform simulations comparing the proposed ISD filter with its ESD counterpart. We consider misspecified settings in which the filters are score driven (ISD or ESD), while the DGP is a state-space model. The various settings below are chosen to highlight specific advantages of ISD over ESD filters. Unless stated otherwise, we use $1{,}000$ replications of length $T$; the first $R$ observations are used to estimate the static parameters. Performance is evaluated through MSEs of the predicted parameters $\{\theta_{t\mid t-1}\}$ and $\{\theta_{t\mid t-1}^{\textnormal{ex}}\}$ relative to the true parameters $\{\theta_t^0\}$, computed in-sample for $t=1,\ldots,R$ and out-of-sample for $t=R+1,\ldots,T$. \subsection{Dynamic Poisson distribution: Non-Lipschitz gradient} For each time $t$, $y_t\in \mathbb{N}$ is drawn from a Poisson distribution with a time-varying intensity $\lambda^0_t:=\exp(\theta^0_t)$, i.e.,\ $p(y_t|\theta_t^0)=(\lambda^0_t)^{y_t}\exp(-\lambda^0_t)/y_t!$. The score with respect to the log-intensity parameter $\theta\in \mathbb{R}$ is $y_t-\exp(\theta)$. The negative Hessian and Fisher information are both $\exp(\theta)>0$; hence, the density is log-concave in $\theta$. We specify the state dynamics as $\theta_t^0=0.98\theta_{t-1}^0+\eta_t$, where $\eta_t \mathbin{\overset{\text{i.i.d.}}{\kern\z@\sim}}\mathrm{N}(0,\sigma_\eta^2)$, and we vary the value of $\sigma_\eta$. We consider ISD and ESD filters based on the (correctly specified) Poisson distribution with exponential link, with $R=2{,}000$ and $T=10{,}000$. We follow koopman2016predicting, who used a time-varying learning rate (for the ESD filter) that scales with the inverse square root of the (model-based) predicted Fisher information quantity; similarly, we set $H^{j}_t=H^{j} \exp(-\theta^{j}_{t|t-1}/2)$ with a static learning-rate coefficient $H^{j}>0$ for $j\in \{\text{im},\text{ex}\}$. For the ESD filter, the scaled score $\exp(-\theta/2)(y_t-\exp(\theta))$ is non-Lipschitz in $\theta\in \mathbb{R}$. \begin{table}[h!] \caption{MSEs of predicted states with dynamic Poisson distribution. } \begin{threeparttable} \begin{tabular}{llrrrrr} & & $\sigma_\eta=0.10$ & 0.15 & 0.20 & 0.25 & 0.30 \\ \hline \multicolumn{1}{l}{ISD filter} & In-sample & 0.088 & 0.149 & 0.221& 0.304 & 0.398 \\ & Out-of-sample & 0.091 & 0.154 & 0.227 & 0.311 & 0.407\\ \hline \multicolumn{1}{l}{ESD filter} & In-sample & 0.088 & 0.149 & 0.220 & 0.311& 0.439 \\ & Out-of-sample & 0.091& 0.153 & 0.227 & $\infty$ & $\infty$ \end{tabular} \end{threeparttable} \end{table} Table (ref) shows the in-sample MSEs to be comparable. For both filters, MSEs increase with the state variability $\sigma_{\eta}$. The out-of-sample performance is similar only for $\sigma_{\eta}\lesssim 0.20$. For $\sigma_\eta=0.25$ and $0.30$, the ISD filter's out-of-sample performance remains aligned with the in-sample results, but the ESD filter diverged in 2% and 26% of replications, respectively. Larger state innovations induce larger prediction errors, such that the (non-Lipschitz) ESD update may lead to divergence. In their simulation study, koopman2016predicting used $\sigma_\eta=0.15$, such that the potential instability of the ESD filter went unnoticed. \subsection{Dynamic GED: Non-concave and non-Lipschitz gradient} At each time $t$, $y_t\in\mathbb{R}$ follows a generalized error distribution (GED) with time-varying mean $\theta_t^0$. Hence $p^0(y_t|\theta^0_t)=\upsilon \exp(-|(y_t-\theta^0_t)/\sigma|^\upsilon)/(2\sigma\Gamma(\upsilon^{-1}))$, where $\Gamma(\cdot)$ is the Gamma function and $\sigma,\upsilon>0$ are static shape parameters. Here $\sigma$ is a scale parameter, but not the standard deviation. We set $\sigma^2=\Gamma(\upsilon^{-1})\Gamma(3\upsilon^{-1})$, such that, conditional on $\theta_t^0$, the variance of $y_t$ is unity. We vary $\upsilon$, with $\upsilon=1$ corresponding to Laplace and $\upsilon=2$ to Gaussian. For $\upsilon>1$, the log-density is $C^1$ and concave in $\theta\in\mathbb{R}$; for $\upsilon<1$ it is neither. At $\upsilon=1$, it is concave but not $C^1$. The gradient is Lipschitz only for $\upsilon\in(1,2]$. A comparative advantage of the ISD update $\theta_{t|t}$ is that it must lie between $\theta_{t|t-1}$ and $y_t$ at every step. In contrast, the ESD update may overshoot $y_t$, meaning $ |\theta^{\textnormal{ex}}_{t\mid t}-\theta^{\textnormal{ex}}_{t\mid t-1}|>|y_t-\theta^{\textnormal{ex}}_{t\mid t-1}| $ may regularly occur. The true parameter evolves as $\theta_t^0=0.98\,\theta_{t-1}^0+\eta_t$, with $\eta_t\mathbin{\overset{\text{i.i.d.}}{\kern\z@\sim}}\mathcal{N}(0,1)$. The signal-to-noise ratio is roughly one, since the state-innovation and observation-noise variances are both unity. We simulate series with $T=10{,}000$ and use the first $R=1{,}000$ observations to estimate the autoregressive parameter $\Phi\in(-1,1)$ in (ref) and the learning rate $H>0$. For simplicity, we set $\omega=\omega^0=0$ and $\upsilon=\upsilon^0$ (i.e., using the true values). \begin{table}[h!] \caption{ MSEs of predicted states with dynamic GED.} \begin{tabular}{rlrrrrrrrr} & & ${\upsilon}=0.5$ & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1.5} & \multicolumn{1}{c}{2} & \multicolumn{1}{c}{2.5} & \multicolumn{1}{c}{3} & \multicolumn{1}{c}{3.5} & \multicolumn{1}{c}{4} \\ \hline \multicolumn{1}{l}{ISD filter} & In-sample & 1.55 & 1.64 & 1.58& 1.57 & 1.57 & 1.58 & 1.58 & 1.59 \\ & Out-of-sample & 1.64 & 1.67 &1.59 &1.58 & 1.58 & 1.59& 1.59& 1.60 \\ \hline \multicolumn{1}{l}{ESD filter} & In-sample & 4.60 & 1.90 & 1.64 & 1.57 & 1.61 & 1.72 & 1.86 & 2.08 \\ & Out-of-sample & 22.68 & 2.11 & 1.65 & 1.58 & 1.62 & $\infty$ & $\infty$ & $\infty$ \\ \end{tabular} \end{table} Table (ref) reports in- and out-of-sample MSEs. For the Gaussian case (i.e., $\upsilon=2$), the performance is identical. Generally, the ISD filter attains lower MSEs (in the range ${\sim}1.55$–${\sim}1.67$ for all $\upsilon$), with out-of-sample MSEs closely tracking in-sample MSEs. For the ESD filter, in-sample MSEs are always higher, and the ESD filter diverged in $2\%$, $58\%$, and $80\%$ of out-of-sample replications for $\upsilon=3$, $3.5$, and $4$, respectively, yielding infinite MSEs. This occurs because the gradient is polynomial of degree $\upsilon-1$ in the prediction error---i.e., excessively large for inaccurate predictions---so poor predictions induce poor updates and, possibly, divergence. For $\upsilon<1$, the MSE of the ESD filter remains finite (as the gradient is unbounded only for near-perfect predictions), but MSEs can still be very large (e.g., ${\sim}23$ for $\upsilon=0.5$). \subsection{Dynamic Gamma distribution: Two time-varying parameters} For each time $t$, an observation $y_t\in \mathbb{R}_{>0}$ is drawn from a Gamma distribution with two dynamic parameters $a_t\in \mathbb{R}_{>0}$ and $b_t\in \mathbb{R}_{>0}$, which are collected in the (true) state vector $\theta^0_t=(a_t,b_t)'\in \mathbb{R}_{>0}^2$, such that $p(y_t|\theta^0_t)=(b_t)^{a_t} y_t^{a_t-1}\exp(-b_t y_t)/\Gamma(a_t)$. The same parametrization is used in fearnhead2004exact, albeit in a static context. Conditional on $\theta_t^0$, the mean and variance of $y_t$ are $a_t/b_t$ and $a_t/b_t^2$, respectively. The ISD filter can be applied directly to $\theta_t=(a_t,b_t)'$, because optimization (ref) guarantees that both elements remain positive: the optimization domain $\Theta$ is the positive quadrant in $K=2$ dimensions. For each $y_t$, the Gamma log density is (jointly) concave in $(a_t,b_t)\in \mathbb{R}_{>0}^2$; hence, Theorem (ref) guarantees ISD filter stability. For the ESD filter, positivity of the time-varying parameters must be enforced through exponential link functions. The resulting log density is neither concave nor $L$-smooth in the transformed (i.e., logarithmic) parameters, meaning no theoretical guarantees can be made for the ESD filter. \begin{figure}[h!] \center \begin{tabular}{cc} & \\ (a) True and predicted mean of $y_t$ & (b) True and predicted variance of $y_t$ \end{tabular} \caption{Illustration of filtering performance with dynamic Gamma distribution.} \end{figure} Figure (ref) shows filtering results for a single time series \(\{y_t\}\) of length \(T=5{,}000\) generated by the non-smooth DGP \(a_t=2+\text{sign}\{\sin(2\pi t/400)\}\) and \(b_t=8+3\,\text{sign}\{\cos(2\pi t/1000)\}\). Only the first \(R=500\) observations are used to estimate the static parameters \(\omega\), \(\Phi\) (diagonal), and \(P\) (positive definite). The figure plots the theoretical mean and variance, \(a_t/b_t\) and \(a_t/b_t^2\), together with their ISD and ESD predictions. The two filters perform similarly in sample, but out of sample ESD diverges, whereas ISD continues to track both moments relatively well. \subsection{Dynamic Dirichlet distribution: Fat-tailed state innovations} At each time $t$, the observation $y_t\in[0,1]^N$ is drawn from a homogeneous Dirichlet distribution with dynamic concentration parameter $\lambda_t^0=\exp(\theta_t^0)>0$ and density $p^0(y_t|\lambda^0_t)=\Gamma(\lambda_t^0 N) / \Gamma(\lambda_t^0)^N \prod_{i=1}^N y_{it}^{\lambda^0_t-1}$, where $\Gamma(\cdot)$ denotes the Gamma function. The $N$ elements of $y_t$ are nonnegative and sum to one. The log density is concave in $\theta_t^0$, but the gradient is not Lipschitz, meaning no stability guarantees exist for the ESD filter. The true process is $\theta_t^0=\omega^0+\Phi^0 \,\theta_{t-1}^0+\sigma_\eta\eta_t$, with i.i.d.\ increments $\{\eta_t\}$ of unit variance, distributed normally or as a Student's \emph{t} with $\nu^0=5$. We set $\omega^0=0.1$, $\Phi^0=0.95$, and $\sigma_\eta^2=0.195$; hence, $\mathbb{E}(\theta_t^0)=\mathbb{V}(\theta_t^0)=2$. We vary $N$ from $2$ to $100$. For each $N$, we simulate $100$ series of length $T=10{,}000$, using the first $R=1{,}000$ observations to estimate $\omega$, $\Phi$ and $H$. \begin{figure}[h!]\center \caption{MSE of out-of-sample predicted states for $N$-dimensional observations from dynamic Dirichlet distribution with Gaussian or Student's $t$ state innovations.} \end{figure} Figure (ref) shows the out-of-sample MSEs of predicted states for ISD and ESD filters under Gaussian and fat-tailed state innovations. Across both settings and for all $N$, the ISD filter outperforms its ESD counterpart. Under fat-tailed state innovations, the ISD filtering performance remains accurate, whereas the ESD filter exhibits unreliable behavior. Moreover, ESD filtering performance does not consistently improve as $N$ increases. In particular, the ESD filtering MSEs sometimes exceed the unconditional variance $\mathbb{V}(\theta_t^0)=2$. Abrupt state changes combined with a non-Lipschitz gradient appear to prevent the ESD filter from exploiting the greater information content of higher-dimensional draws. \subsection{Dynamic network flows: High-dimensional state space} At each time \(t\), we observe an \(N\times N\) matrix of bilateral trade or information flows, \(Y_t=(Y_{ij,t})\), on a directed network with \(N\in\{3,5,10,20,30,40,50\}\) nodes; self-flows are excluded. Node-specific sender and receiver effects are denoted by \(\theta_{1,t}^0,\theta_{2,t}^0\in\mathbb{R}^N\). The full state vector is \(\theta_t^0=\big((\theta_{1,t}^0)',(\theta_{2,t}^0)'\big)'\in\mathbb{R}^{2N}\), so the state dimension is \(K=2N\). Conditional on \(\theta_t^0\), the flows are independent Gamma random variables with common shape parameter \(\psi_0=5\) and mean matrix \(M_t^0=\exp\!\big(\bar{m}+\theta_{1,t}^0\iota' + \iota(\theta_{2,t}^0)'\big)\), where \(\bar m\in \mathbb{R}^{N\times N}\) is a baseline matrix of logarithmic means, \(\iota\) is a vector of ones, and the exponential is applied elementwise. This specification of \(M_t^0\) is closely related to standard gravity models of bilateral trade with exporter/importer effects; e.g., anderson2003gravitas, head2014gravity, and fally2015structural. The latent state evolves as \(\theta_t^0=\Phi_0\theta_{t-1}^0+\eta_t\), with \(\eta_t\mathbin{\overset{\text{i.i.d.}}{\kern\z@\sim}}\mathcal N(0,\Sigma_{\eta})\), where \(\Phi_0\) and \(\Sigma_{\eta}\) are diagonal. For each replication, the diagonal entries of \(\Phi_0\) are drawn independently from a uniform random variate \(\mathcal U[0.80,0.995]\), the standard deviations are drawn from \(\mathcal U[0.2,0.5]\), and the initial state is drawn from its unconditional distribution. The off-diagonal entries of \(\bar m\) are drawn from \(\bar m_{ij}\mathbin{\overset{\text{i.i.d.}}{\kern\z@\sim}} \mathcal N(0,0.40^2)\). As our focus here is on filtering performance, we treat \(\bar m\), \(\Phi_0\), and \(\psi_0\) as known, while the particle filter (PF) also uses \(\Sigma_{\eta}\). We compare the PF with \(10^3\) and \(10^4\) particles (denoted PF3 and PF4) to ISD and ESD filters with estimated scalar learning-rate parameters. The ISD update is implemented by Newton's method with backtracking (see Appendix (ref)). Since the Gamma log density is concave in \(\theta_t^0\), our theoretical guarantees apply to the ISD filter. However, since the gradient is not Lipschitz, no guarantees for the ESD filter are available. For each \(N\), we generate \(100\) replications of length \(T=1{,}000\). For the ISD and ESD filters, the first \(R=50\) observations are used to tune the scalar learning-rate parameters, and the remaining \(T-R=950\) observations are used for out-of-sample evaluation. Accuracy is measured by computing the MSE of the predicted log-intensity surface, \(\theta_{1,t|t-1}\iota' + \iota\theta_{2,t|t-1}'\), excluding diagonal elements, relative to its true simulated counterpart. \begin{figure}[t!] \center \begin{tabular}{cc} & \\ (a) MSE & (b) Time (seconds) \end{tabular} \caption{Out-of-sample prediction MSE and filtering time for the dynamic network model with $N$ nodes, \(N\times N\) observations and \(K=2N\) states.} \end{figure} Figure (ref)(a) shows that the ISD filter is competitive for all \(N\) while clearly outperforming the alternatives for \(N\geq 10\). For \(N=50\), its MSE is about three times lower than that of the ESD filter and about four to five times lower than that of PF3 and PF4. It is also the only method whose performance improves with \(N\). For \(N=10,20,30,40,\) and \(50\), the ESD filter exhibits instability: among the 100 replications, between 1 and 8 runs diverge, while no divergence occurs for \(N=3\) or \(5\). The reported MSE for the ESD filter is based only on the non-divergent runs. Figure (ref)(b) shows that the ESD filter is fastest, as expected, followed by the ISD filter. For \(N=50\), despite solving a 100-dimensional optimization problem over \(1{,}000\) time steps, the ISD filter takes about \(2.5\) seconds per replication, making it three times faster than PF3 and thirty times faster than PF4. To interpret these results, we distinguish statistical accuracy from computational cost. In terms of accuracy, the particle filter suffers from the curse of dimensionality: although the conditional mean is optimal under MSE loss, approximating it becomes prohibitively difficult as \(N\) grows (e.g., bengtsson2008curse). By contrast, the ISD update exploits the richer cross-sectional information contained in each observation, so its performance improves with \(N\). In terms of computation, the ISD filter solves at each time step a smooth, strongly concave optimization problem, initialized at the one-step-ahead prediction, which is typically close to the optimum. As a result, only a few Newton steps are usually needed, even in high dimensions. Hence, in this design, the ISD filter is both more accurate and faster than the PF. These findings are consistent with the arguments and simulation evidence in lange2024bellman, which suggest that, in high-dimensional settings, optimization-based methods such as the ISD filter can outperform integration-based methods such as the PF, even when the object of interest is the conditional mean. \section{Empirical illustrations} We consider three applications for which the ISD update (ref) allows a closed-form solution. \subsection{Linear regression with time-varying slope} The capital asset pricing model (CAPM), a standard benchmark in finance, links expected excess returns of individual assets linearly to those of the market. However, empirical evidence (e.g., jagannathan1996conditional) shows that assuming a constant market coefficient $\beta$ can be unrealistic, especially in equity markets. We study possible time variation in the CAPM market $\beta$ using the ISD filter. We model the excess asset return $y_t$ as \begin{equation} y_t \;= \;\alpha \,+\, \beta_{t}\,x_t \,+\, \varepsilon_t, \qquad \varepsilon_t \mathbin{\overset{\text{i.i.d.}}{\kern\z@\sim}}\textnormal{N}(0,\sigma^2), \end{equation} where $\alpha$ is a static intercept, $x_t$ is the excess market return, and $\varepsilon_t$ is an i.i.d.\ Gaussian disturbance with mean zero and variance $\sigma^2>0$. The ISD update (ref) applied to a prediction $\beta_{t|t-1}$ for $\beta_t$ in (ref) can be solved analytically (see Appendix (ref) for details), yielding \begin{equation} \beta_{t|t} \;= \;\beta_{t|t-1} \,+\, \frac{\sigma^2}{\sigma^2 + H x_t^2}\, H\, \nabla(y_t | \beta_{t|t-1}, x_t), \end{equation} where $H = P^{-1}>0$ is a scalar learning-rate parameter and $\nabla(y_t | \beta_{t|t-1}, x_t):=x_t( y_t - \beta_{t|t-1}\,x_t)/\sigma^2 $ is the explicit score (i.e., evaluated at $\beta_{t|t-1}$). For the prediction step, we use the linear first-order specification (ref). Equation (ref) illustrates Proposition (ref) for linear regression: its right-hand side includes the shrinkage factor $\sigma^2/(\sigma^2 + H x_t^2)\in(0,1]$, which is absent (i.e., equal to one) in the ESD update. In the ISD update, larger values of $x_t^2$ induce greater shrinkage. Indeed, the specific functional forms of the explicit score $\nabla(y_t | \beta_{t|t-1}, x_t)$ and the shrinkage factor as a function of $x_t$ imply that $\beta_{t|t}\to 0$ as $|x_t|\to \infty$. In fact, the ISD update (ref) is bounded in both $x_t^2$ and the learning rate $H>0$. The latter follows from $H$ appearing both in the shrinkage factor and in front of the score. This boundedness in $H$ is not limited to the present regression example: ISD filters are robust under a wide range of learning rates, whereas ESD filters require more careful fine-tuning. In stochastic-gradient methods, similar arguments arise when comparing the least-mean-square (LMS) filter to its normalized version (NLMS; see nagumo1967learning and diniz1997adaptive). We apply the ISD filter to the dynamic regression model (ref) with daily excess returns of Microsoft (MSFT) from 14 March 1986 to 29 April 2022, obtained from Yahoo Finance.\footnote{\url{https://finance.yahoo.com/quote/MSFT/history?p=MSFT}} For the market return and risk-free rate, we use the series from Kenneth French's database.\footnote{\url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}} Figure (ref)(a) shows the evolution of $\beta_{t|t-1}$ for the ISD and ESD filters. In Figure (ref)(b), we plot the slope adjustments $\beta_{t|t}-\beta_{t|t-1}$ as a function of $x_t$ for a fixed $y_t = 0$ and two predictions ($\beta_{t|t-1} = 1$ and $\beta_{t|t-1} = -0.5$); we refer to these as “impact curves.” \begin{figure}[t!] \center \begin{tabular}{cc} & \\ (a) Time evolution of $\beta_{t|t-1}$ & (b) Estimated impact curve \end{tabular} \caption{Dynamic CAPM using the ISD and ESD filters for MSFT from March 1986 until April 2022. Vertical dotted lines in panel (a) mark Black Monday on October 19, 1987.} \end{figure} In panel (a), the ISD and ESD filters display similar paths $\{\beta_{t\mid t-1}\}$, although the ESD filter recovers more slowly from large shocks such as Black Monday (1987). Panel (b) explains why: the ESD impact curve is quadratic (and hence unbounded) in the exogenous variable $x_t$, which necessitates a small learning rate to curb outliers. Accordingly, its estimated rate is only $\hat{H}=0.0092$, suppressing responsiveness in the rest of the sample (e.g., around 1994 and 2004, when the ESD filter lags behind its ISD counterpart). By contrast, the ISD filter is not overly sensitive to outliers like Black Monday and is more responsive to persistent movements (e.g., the downward trend around 2004). Because the ISD impact curve redescends to zero as $|x_t|\to\infty$, it admits a larger estimated learning rate, $\hat{H}=0.0169$; roughly twice that for the ESD filter. This allows the ISD filter to pick up gradual changes more quickly, without being thrown off by large values of $x_t^2$. \subsection{Time-varying growth-at-risk} For policymakers, monitoring macroeconomic downside risk is crucial. A popular approach is growth-at-risk (GaR), which focuses on the conditional lower quantiles of GDP growth. GaR and its link to financial/economic conditions are typically estimated via quantile regressions (QRs; koenker2001quantile) with exogenous predictors (e.g., adrian2019vulnerable). We update a time-varying conditional quantile by postulating an asymmetric Laplace distribution with a time-varying location. Maximizing this density is equivalent to minimizing koenker1978regression's (koenker1978regression) QR check function (see koenker1999goodness). The ESD update for the $\tau$-level quantile with $0<\tau<1$, denoted by $q^\textnormal{ex}_{t|t}(\tau)$, is \begin{equation} q^\textnormal{ex}_{t|t}(\tau) = q_{t|t-1}(\tau) - 1[y_t < q_{t|t-1}(\tau)] \frac{H(1-\tau)}{\sigma} + 1[y_t > q_{t|t-1}(\tau)] \frac{H\tau}{\sigma}, \end{equation} where $y_t$ is GDP growth, $1[\cdot]$ is the indicator, and $H>0$ and $\sigma>0$ are constant learning-rate and dispersion parameters. The ESD update (ref) adjusts the quantile by $-H(1-\tau)/\sigma$ when $y_t<q_{t|t-1}(\tau)$ and by $+H\tau/\sigma$ when $y_t>q_{t|t-1}(\tau)$. If $y_t=q_{t|t-1}(\tau)$, no adjustment is made, such that $q^\textnormal{ex}_{t|t}(\tau)=q_{t|t-1}(\tau)$; in this case, the ESD update allows a subgradient interpretation. Aside from this probability-zero event, the ESD quantile update $q^\textnormal{ex}_{t|t}(\tau)$ coincides with (the limiting version of) engle2004caviar's (engle2004caviar) adaptive CAViaR update. Due to their design, ESD/CAViaR updates can overshoot $y_t$. While the ISD update $q_{t|t}(\tau)$ is similar to the ESD/CAViaR update, it cannot surpass the observation $y_t$: \begin{equation} q_{t|t}(\tau) = \left\{ \begin{array}{l@l} \min\{y_t,q^\textnormal{ex}_{t|t}(\tau)\}, & y_t > q_{t|t-1}(\tau),\\ \max\{y_t,q^\textnormal{ex}_{t|t}(\tau) \}, & y_t\leq q_{t|t-1}(\tau). \end{array} \right. \end{equation} As before, the ISD update is a shrunken version of its ESD version: if $y_t$ is above (below) $q_{t|t-1}(\tau)$, ISD follows the ESD update upward (downward), but never beyond $y_t$. This capping of ISD updates is advantageous, as updating beyond $y_t$ would decrease model fit. When tracking multiple quantiles with QRs, quantile crossing poses a serious problem. For ESD/CAViaR updates, crossings cannot generally be avoided. For the ISD update (ref), as we show in Appendix (ref), quantiles remain properly ordered if all updates use the same $H$ and $\sigma$. To see why, note that for an observation lying between any two predicted quantiles, neither update can surpass the observation, thus preventing the quantiles from crossing. To ensure that quantiles also remain ordered in the prediction step, we specify \begin{equation} q_{t+1|t}(\tau) \;=\; c(\tau)\,(1-\Phi) \,+\, \Phi\, q_{t|t}(\tau)\, +\, \gamma\, x_t, \end{equation} with a common autoregressive parameter $\Phi\in[0,1)$, long-run quantile levels $c(\tau)$ that are strictly ordered in $\tau$, and a common slope $\gamma$ relating to the exogenous predictor $x_t$.\footnote{While it may be useful to allow different quantiles to have different sensitivities to the exogenous input $x_t$, this has the potential to introduce quantile crossings. Moreover, for our application, the likelihood improvement of quantile-specific slopes $\gamma(\tau)$ is too small to justify the additional model complexity.} We estimate the 5, 10, 25, and 50 percent GaR using ISD and ESD filters. For $\{y_t\}$, we use quarterly U.S.\ GDP growth from 1971Q1 to 2021Q4, while for $\{x_t\}$ we take the National Financial Conditions Index (NFCI), with quarterly values obtained by averaging weekly data (cf. adrian2019vulnerable). Both series are available in the FRED database.\footnote{See \url{https://fred.stlouisfed.org/series/GDP} and \url{https://fred.stlouisfed.org/series/NFCI}.} To reduce the number of parameters, we use targeting (e.g., engle2002dynamic) and set $c(\tau)$ to its full-sample empirical counterpart. The remaining static parameters are estimated using composite likelihood, comparable to zou2008composite. We fix $\sigma=1$, which can be treated as a nuisance parameter (geraci2007quantile). The postulated log-likelihood equals the sum of four Laplace log-densities, of which three are asymmetric and one is symmetric. \begin{figure}[t!] \hbox \caption{Growth-at-risk estimates for the ISD and ESD models for $\tau = 0.05$, $\tau = 0.10$, $\tau = 0.25$, and $\tau = 0.50$, 1971Q1 to 2021Q4.} \end{figure} Figure (ref) displays the 5, 10, 25, and 50 percent GaR estimates from the ESD/CAViaR and ISD filters, involving (ref) and (ref), respectively. The ISD-filtered quantiles are more responsive, showing larger downward adjustments at the onset of COVID-19 (2020Q2) as well as faster post-crisis recovery. The enhanced stability of the ISD update (ref) allows its estimated learning rate to exceed that of its explicit counterpart by a factor of five ($\hat{H}=4.002$ vs.\ $\hat{H}=0.804$). While the ESD-filtered quantiles exhibit regular crossings, the ISD-filtered versions remain properly ordered. For instance, the ESD median ($\tau=0.50$) often overshoots $y_t$ and occasionally crosses the first quartile ($\tau=0.25$). In line with adrian2019vulnerable (adrian2019vulnerable), we find a negative NFCI effect ($\hat{\gamma}=-0.052$ for ISD; $\hat{\gamma}=-0.019$ for ESD), implying that tighter financial conditions exert downward pressure on the GDP quantiles. \subsection{T-bill rate spreads} Monitoring term spreads---leading indicators of business-cycle recessions---is essential for companies and policymakers alike. We use harvey2014filtering's (harvey2014filtering) Student's $t$ location model, which has been widely used (e.g., caivano2016robust; blasques2018feasible; blasques2022maximum), including for T-bill rate spreads (artemova2022score1). However, the Student's $t$ density is not log-concave in the location parameter (i.e., Assumption (ref) fails). Moreover, if the penalty parameter $P$ is sufficiently small, optimization (ref) is not strongly concave (i.e., Assumption (ref) may also fail). Here we show that the ISD filter can perform well even if Assumptions (ref) and (ref) fail. The ISD update for the Student's $t$ location model admits a closed-form solution: it reduces to solving a cubic (for details, see Appendix (ref)). Standard software provides closed-form (albeit lengthy) expressions for up to three roots, of which at least one is real valued. The lack of concavity means that up to three stationary points (two of them local maxima) may exist; in this case, we select the global maximum by comparing objective values in (ref). Multiplicity of stationary points is precluded if the penalty term is strong enough for Assumption (ref) to hold, i.e., if $P>1/8$ or, equivalently, $H<8$ (for details, see Appendix (ref)). As the update is well-defined for all $P>0$, we do not impose this condition. The ISD update $\theta_{t|t}$ is a weighted average of the prediction $\theta_{t|t-1}$ and the observation $y_t$. The weight of $y_t$ lies in the interval $(0,H/(1+H)]\subseteq (0,1)$ for any $H>0$ (see Appendix (ref)), such that $\theta_{t|t}\in [\min\{y_t,\theta_{t|t-1}\},\max\{y_t,\theta_{t|t-1}\}]$. In contrast, the ESD update with $H>1$ can overshoot the data and lead to “zig-zagging,” as we will see below. Following artemova2022score1, we analyze the spread between three- and six-month T-bill rates from the FRED-QD database (series “TB6M3M”) scaled by a factor of ten. The sample comprises 249 quarterly observations from 1959:Q1 to 2021:Q1. For these data, artemova2022score1 found that the ESD filter with a Student's $t$ density performed best. \begin{table}[H] \caption{Parameter estimates for Student's $t$ location model with T-bill spread data.} \begin{threeparttable} \begin{tabular}{cc@ccccc@c} Distribution & Filter & $\hat{\omega}$ & $\hat{\Phi}$ & $\hat{H}$ & $\hat{\sigma}^2$ & $\hat{\nu}$ & LL \\ \hline Student's $t$ & ISD & 0.235 & 0.751 & 23.717 & 0.387 & 2.060& $-364.4$ \\ & ESD & 0.353 & 0.714 & 2.194 & 0.516 & 2.632 & $-370.6$ \\ \hline \end{tabular} \end{threeparttable} \end{table} Table (ref) reports ML estimates and log-likelihood values for ESD and ISD filters based on the Student's $t$ with scale $\sigma>0$ and $\nu>0$ degrees of freedom. The ISD filter achieves a higher likelihood and a larger learning rate ($\hat{H}\approx 24$ vs.\ $\hat{H}\approx 2.2$). Interestingly, it fits the data best when the regularized objective (ref) is non-concave (as $\hat{H}\approx 24>8$). \begin{figure}[H] \begin{tabular}{c@c} & \\ (a) T-bill rate spread and predicted paths & (b) Impact curves ($y_t-\theta_{t|t-1}$ vs.\ $\theta_{t|t}-\theta_{t|t-1}$) \end{tabular} \caption{T-bill spread data with filtered paths and impact curves.} \end{figure} Figure (ref) shows the associated predicted-parameter paths and impact curves. The ESD path is relatively unresponsive in the 1960s–1990s, but overly responsive after 2010, producing a zig-zagging pattern around the data, consistent with $\hat{H}\approx 2.2>1$. By contrast, the ISD filter's impact curve is more responsive to large shocks without overshooting the data, allowing it to track peaks and troughs in the volatile decades, without zig-zagging after 2010. Both filters largely ignore the 1982:Q3 outlier (10.6). For the ISD impact curve, the discontinuity around this observation reflects a bifurcation point, where the second local maximum starts to dominate. The non-concavity of the objective function (ref) due to $\hat{H}>8$ and the associated multiplicity of stationary points allow the ISD filter to largely disregard this outlier. Relative to the ESD filter, the ISD filter lowers the in-sample MSE by roughly 18%, which may be economically significant. \section{Conclusion} We have introduced an implicit score-driven (ISD) framework for updating time-varying parameters in observation-driven models. At each time step, the ISD filter maximizes the logarithmic observation density subject to a quadratic penalty centered at the one-step-ahead prediction; the associated first-order condition is an implicit stochastic-gradient update. For the class of (possibly misspecified) log-concave densities, we showed that the ISD filter is invertible, while its updates are contractive in mean squared error toward the \mbox{(pseudo-)}true parameter at every time step. The popular class of explicit score-driven (ESD) filters---also known as dynamic conditional score (DCS; harvey2013dynamic) models or generalized autoregressive score (GAS; creal2013generalized) models---arises within our framework by locally linearizing the logarithmic density around the prediction. However, to obtain guarantees comparable to those of the ISD filter, especially under model misspecification, ESD filters require the score to be Lipschitz continuous and the learning rate to be sufficiently small. We illustrated the theoretical benefits of our approach in simulations and its practical relevance in empirical applications to asset prices, growth-at-risk, and T-bill rates. { \putbib[ref] }