Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
67,747 characters · 18 sections · 51 citation commands
Exponentially weighted estimands and the exponential family: Filtering, prediction and smoothing
Keywords: Exponential family; EWMA; Filtering; Likelihood; Time Series.
\baselineskip=20pt
We provide a simple way to carry out filtering, predicting and smoothing. It is based on a discounted version of a convex combination of a log-likelihood and the expected log-likelihood function. The discounting is carried out using an exponential function.
We focus on log-likelihoods based around the canonical exponential family. It has a similar structure to generalized linear regression models, but now for time series. Properties of the resulting procedures are established and illustrated by simulation and empirical work.
Section (ref) develops the new filter, predictor and smoother. Section (ref) looks at the canonical exponential family case. These are found by a convex numerical optimization.
Section (ref) looks at a flexible special case of the canonical exponential family which results in analytic expressions for the filters, predictors and smoothers. They do not need any form of iteration or approximation. This special case covers a vast number of familiar models, e.g. based on distributions such as, for example, the Gaussian, multinomial, Poisson, Pareto, exponential, von Mises and Dirichlet.
Section (ref) develops a quasi-likelihood based estimation strategy for the hyperparameters of this structure. This method is extended to allow for a two-step estimator, where an easy to use moment estimator reduces the dimension of the required numerical maximization.
Section (ref) provides an empirical illustration based on a seven category Dirichlet model. Section (ref) details some conclusions.
Think of time series data $y_{1:T} = \{y_1,...,y_T\}$ and a time series model for the corresponding random variables $Y_{1:T} = \{Y_1,...,Y_T\}$. Here $T$ is the length of the time series. The properties of the time series model $Y_{1:T}$ are determined by the abstract probability triple $(\Omega,\mathcal{F},P)$.
At no point do we regard this frame as true. It will index a class of filters, smoothers and predictors. The frame may depend on additional static parameters $\phi$, which we suppress notationally for clarity until Section (ref). The form of the frame can vary over time, e.g. it could be based around the random variable $Y_t \sim {\tt Multinomial}(n_t,p)$, where $p$ is determined by $\theta$ and $n_t$ is regarded as non-stochastic (e.g. we have formed the probability mass function for $Y_t$ by conditioning on $n_t$) or $Y_t \sim {\tt Pois}(e^{\theta^{\tt T} x_t})$ where $x_t$ is a non-stochastic vector (or conditioned on).
At various points in the paper some assumptions are made about the form of the frame.\footnote{We follow the statistical convention of taking the log-likelihood as ignoring terms which only involve data not also hyperparameters. One way of implementing this, for a model with a probability density function, for example, is to take $\log L_t(\theta;y_t) := \log f_t(y_t;\theta) - \log f_t(y_t;\theta_0)$ where $\theta_0 \in \Theta$ is some arbitrary fixed value as $\theta\in \Theta$ varies. As we need the expected log-likelihood to exist, removing these data terms lessens the requirements we need to prove the expected value exists.} Here we list some of them.
The expectation is under the probability measure for $Y_t$, not the probability model that would rationalize the frame.
The paper will be based around an exponential version of a weighting principle.\footnote{The extension to allow a prior density $\pi(\theta)$ for $\theta \in \Theta$, yielding an exponentially weighted MAP (i.e. mode of a type of posterior), means just adding $\log \pi(\theta)$ to our criterion function and then maximizing that with respect to $\theta$ for each $t$. This is particularly important if $\theta$ is high dimensional.}
The principle of weighting delivers three sequences of estimands through time: a filter, a predictor and a smoother. Notice that the estimands are random variables.
Note that Definition (ref), Assumptions (ref)--(ref) and $\alpha \in [0,1)$ imply that $Q_{\lambda,t}(\theta)$, $Q_{\lambda,t|t-s}(\theta)$ and $Q_{\lambda,t|T}(\theta)$ admit unique maximizers $\widetilde{\theta }_{t}$, $\widetilde{\theta }_{t|t-s}$ and $\widetilde{\theta }_{t|T}$, respectively.
As the likelihood function is invariant to reparameterization, like the maximum likelihood estimator (MLE), the exponentially weighted filter, predictor and smoother are invariant under one-to-one reparameterization, i.e. if $\gamma = \gamma(\theta)$, then $\widetilde{\gamma}_t = \gamma(\widetilde{\theta}_t).$ This will turn out to be a very powerful property in practice.
Intuitively, $\lambda \in [0,1]$ controls the rate at which older information is downweighted, while $\alpha \in [0,1]$ controls the weight placed on the anchor $\mathbb{E}[\log L_j(\theta; Y_j)]$ terms.
The rest of this paper will focus entirely on the case where the frame is based on a member of the canonical exponential family. This family includes most of the important probability models used in applied statistics. In our context, crucially, it allows either an analytic solution or an easy to compute solution for each $t$.
Start by recalling the definition of the canonical exponential family and three of its properties.
The cumulant function of the {\tt CEF} is $\log \mathbb{E}[e^{a^{\tt T}h(Y)};\theta] = \psi(\theta+a)-\psi(\theta)$, so
the finite expected value of $h(Y)$ under the ${\tt CEF}(\theta,h,\psi)$ model and
the corresponding finite variance. The minimal {\tt CEF} requires $\psi(\theta)$ to be strictly convex so ${\mathrm V}(h(Y);\theta)$ is positive definite under the model and Assumption (ref) is satisfied. It is helpful to go from $\mu = \mathbb{E}[h(Y)]$ to $\theta$ through the inverse function $ \psi'^{-1}(\mu) = \theta. $ This inverse has a unique solution (as $\psi$ is strictly convex), but in some cases it has to be solved numerically; a classic task in statistics as this is isomorphic to computing the MLE of $\theta$ in the {\tt CEF}. Famously, it can be carried out reliably and rapidly using a Newton-Raphson algorithm, which converges quadratically close to the solution. We state this formally as Algorithm (ref) for reference later.
There is a substantial literature on time series models based around the exponential family. Examples include ZegerQaqish(88), Li(94), BenjaminRigbyStasinopoulos(03) and the review by Davis03072021. The dynamic conditional score (DCS) filter and the generalized autoregressive score (GAS) filter of Harvey(13) and CrealKoopmanLucas(13) can be applied in this context to deliver recursive filters. Our predictive version is closest to the generalized {\tt ARMA} work of BenjaminRigbyStasinopoulos(03), which is built around the exponential family and generalized linear models.
In what follows throughout we take the frame $\log L_t(\theta;y) = \log f_t(y;\theta)$ as coming from ${\tt CEF}(\theta,h_t,\psi_t)$. Then, for each $t$, (ignoring the implied $b_t(y)$ term as it has no impact on $\theta$ and dropping it means we do not need to make an assumption that $\mathbb{E}[\log b_t(Y_t)]$ exists) $$ (1-\alpha) \mathbb{E}[\log L_t(\theta;Y_t)] + \alpha \log L_t(\theta;Y_t) = \theta^{\tt T} \{(1-\alpha)\mathbb{E}[h_t(Y_t)] + \alpha h_t(Y_t)\} - \psi_t(\theta). $$
In the case where the data is assumed to be strictly stationary, the frame is stable and $\mathbb{E}[h(Y_t)]$ exists, then the right hand side of this expression simplifies to $$\theta^{\tt T} \{(1-\alpha)\mathbb{E}[h(Y_1)] + \alpha h(Y_t)\} - \psi(\theta).$$
To compactly write the exponentially weighted filter it is helpful to denote the exponentially weighted moving sums $$ x_{\lambda,t} = \sum_{j=1}^t \lambda^{t-j}\mathbb{E}[h_j(Y_j)], \quad h_{\lambda,t} = \sum_{j=1}^t \lambda^{t-j}h_j(Y_j),\\ $$ as well as the double sided exponentially weighted moving sums $$ x_{\lambda,t|T} = \sum_{j=1}^T \lambda^{|t-j|}\mathbb{E}[h_j(Y_j)], \quad h_{\lambda,t|T} = \sum_{j=1}^T \lambda^{|t-j|}h_j(Y_j). $$ Each of these terms can be computed recursively, implying each entire series, e.g. $h_{\lambda,1|T},...,h_{\lambda,T|T}$, can be computed in $O(\sum_{t=1}^T \dim(h_t(y_t)))$ computations.
All the left hand side terms (e.g. $x_{\lambda,t}$ and $h_{\lambda,t}$) can be computed recursively. The right hand side, not so much, as they depend upon the value of $\theta$ that the $\psi_j$ function is evaluated at. Typically $\widetilde{\theta}_{t}$, $\widetilde{\theta}_{t|t-s}$ and $\widetilde{\theta}_{t|T}$ have to be found by numerical root solving. Due to the strict convexity of $\psi$ this is numerically straightforward, for each individual value of $t$, using Newton-Raphson, e.g. Algorithm (ref) computes $\widetilde{\theta}_t$, the exponentially weighted filter.
A downside is that the numerical procedure has to be run separately for each $t$, so computing, for example, the time series $\widetilde{\theta}_{1},...,\widetilde{\theta}_{T}$ costs $O(\{\sum_{t=1}^T \dim(h_t(y_t))^3\}^2)$ calculations.
We will use the following structure in our simulations later.
The use of the frame to derive a filter, which is then used as an input into a data generating process echos the DCS/GAS models of Harvey(13) and CrealKoopmanLucas(13).
Some examples of these types of simulations will be given shortly.
The following special case has an analytic solution. We focus on it in the rest of this paper.
This Example covers many interesting models.
We illustrate the exponentially weighted predictor on ten {\tt CEF} examples. Given $\widetilde{\mu}_{t|t-1}$, each yields a unique solution for $\widetilde{\theta}_{t|t-1} = (\psi')^{-1}(\widetilde{\mu}_{t|t-1})$. In cases 1--6 and 9, this can be solved analytically; in cases 7, 8, and 10, numerically via Algorithm (ref). Tables (ref)-- (ref) summarize the canonical parameterizations.
For case 2, the resulting predictor $\widetilde{\mu}_{t|t-1}$ is the stationary version of the {\tt ARMA}(1,1) process; we compare the filter with the Kalman filter in Section (ref). For case 5, note that $\mu = \mathbb{E}[Y^2] = \mathrm{V}(Y)$, and the resulting predictor $\widetilde{\mu}_{t|t-1}$ is the stationary version of the {\tt GARCH}(1,1) model of Engle(82) and Bollerslev(86). For case 9, $\widetilde{\sigma}^2_{t|t-1} = \widetilde{\mu}_{2,t|t-1} - \widetilde{\mu}^2_{1,t|t-1} \geq 0$ by invariance and Jensen's inequality. For the cases 7, 8, and 10, given $\mu$, the corresponding $\theta = (\psi')^{-1}(\mu)$, which is unique, can be found numerically via Algorithm (ref).
Following Assumption (ref), we simulate each process using $Y_t | Y_{1:t-1} \sim {\tt CEF}(\widetilde{\theta}_{t|t-1},h,\psi)$ with $$ \widetilde{\theta}_{t|t-1} = (\psi')^{-1}(\bar{m}_{t|t-1}) = (\psi')^{-1}\left( \frac{m_{\lambda,t|t-1}}{n_{\lambda,t|t-1}} \right), $$ where $m_{\lambda,t|t-1} = (1-\alpha)\mathbb{E}[h(Y_1)] + \alpha\lambda h(Y_{t-1}) + \lambda m_{\lambda,t-1|t-2}$. We impose stability by centering at $\mathbb{E}[h(Y_1)] = \mathbb{E}[Y_1] = 0.5$ (Bernoulli), $\mathbb{E}[h(Y_1)] = \mathbb{E}[Y_1] = 0$ (Gaussian with known $\sigma^2=1$), $\mathbb{E}[h(Y_1)] = \mathbb{E}[Y_1] = 1$ (Poisson, Exponential), $\mathbb{E}[h(Y_1)] = \mathbb{E}[Y_1^2] = 1$ (Gaussian with zero mean), and $\mathbb{E}[h(Y_1)] = \mathbb{E}[\log(Y_1)] = 1/3$ (Pareto with scale $m=1$, and shape $\theta_0 = -3$). For the bivariate cases, we take $\mathbb{E}[h(Y_1)] = \psi'(\theta_0)$ with initial values $\theta_0 = (2, 5)'$ for Beta (corresponding to ${\tt Beta}(2,5)$), $\theta_0 = (0, -1/2)'$ for Gaussian (corresponding to mean $0$, variance $1$), and $\theta_0 = (0, -2)'$ for von Mises (corresponding to mean direction $\pi$, concentration $2$).
Figures (ref) and (ref) show simulated time series $Y_{t}$ and conditional predictors for $t=5,...,T=2000$ with $\lambda=0.93$, $\alpha \in \{0.70, 0.95\}$, and $\widetilde{\theta}_{1|0} = (\psi')^{-1}(\mathbb{E}[h(Y_1)])$. Figure (ref) displays the conditional mean $\widetilde{\mathbb{E}[Y_t|Y_{1:t-1}]} = \psi'(\widetilde{\theta}_{t|t-1})$ for Bernoulli, Gaussian (known $\sigma^2=1$), and Poisson. Integer-valued observations are jittered with $\text{Unif}(-0.1, 0.1)$ noise for visualization. Figure (ref) shows the exponential distribution with conditional mean $\widetilde{\mathbb{E}[Y_t|Y_{1:t-1}]} = \psi'(\widetilde{\theta}_{t|t-1})$, the Gaussian (zero mean) with conditional standard deviation $\widetilde{\sigma}_{t|t-1} = \sqrt{\psi'(\widetilde{\theta}_{t|t-1})}$, and the Pareto with conditional mean $\widetilde{\mathbb{E}[Y_t|Y_{1:t-1}]} = \frac{\widetilde{\theta}_{t|t-1}}{\widetilde{\theta}_{t|t-1}+1}$ when $\widetilde{\theta}_{t|t-1} < -1$ and $\infty$ otherwise.
Figure (ref) shows simulation results for the three bivariate cases: Beta, Gaussian with time-varying mean and variance, and von Mises. Each column displays the conditional expectation $\widetilde{\mathbb{E}[Y_t|Y_{1:t-1}]}$ with simulated observations. For Beta, the conditional mean is $\widetilde{\theta}_{1,t|t-1}/(\widetilde{\theta}_{1,t|t-1}+\widetilde{\theta}_{2,t|t-1})$. For Gaussian, it is $-\widetilde{\theta}_{1,t|t-1}/(2\widetilde{\theta}_{2,t|t-1})$ and the conditional standard deviation $\widetilde{\sigma}_{t|t-1} = \sqrt{-1/(2\widetilde{\theta}_{2,t|t-1})}$ is shown as a red line. For von Mises, we present the mean direction $\widetilde{\mu}_{t|t-1} = \text{atan2}(\widetilde{\theta}_{1,t|t-1}, \widetilde{\theta}_{2,t|t-1})$, modulo $2\pi$ to ensure $\widetilde{\mu}_{t|t-1} \in [0,2\pi]$, where $\text{atan2}(y,x)$ is the four-quadrant arctangent of $(x,y)$. For $\alpha = 0.95$, the predictors oscillate more around the centering values than for $\alpha = 0.70$, due to weaker mean reversion.
The Gaussian example with known standard deviation, case 2, links strongly to the Kalman filter for the univariate Gaussian local level model
The Kalman filter implies $a_{t|t-1} = {\mathrm E}[\mu_t|Y_{1:t-1}] = {\mathrm E}[\mu_{t-1}|Y_{1:t-1}]$ is given by (e.g. Harvey(89) and DurbinKoopman(12))
where $K_t$ is called the Kalman gain and $\sigma^2_\varepsilon P_t = {\mathrm V}(\mu_t|Y_{1:t-1})$. If $q>0$, in steady state, $K_t$ and $P_t$ converge to $ K = \frac{P}{P + 1}$ and $ P=K+q. $ Then $K^2 + qK - q = 0$, implying $K = \frac{-q + \sqrt{q^2 + 4q}}{2}. $ In steady state, the one-step predictor is $$ a_{t+1|t} = (1-K)a_{t|t-1} + K Y_t= K \sum_{j=0}^{\infty} (1-K)^j Y_{t-j}, $$ which is an {\tt EWMA} with discount hyperparameter $$\lambda=1-K = \frac{2+q-\sqrt{(2+q)^2-4}}{2}.$$ The discount rate monotonically declines as $q$ increases, going from 1 to 0.
Our predictor in steady state, setting $\mathbb{E}[Y_1] = 0$ and $\alpha=1$, implies $$ \widetilde{\mu}_{t+1|t} = \lambda \widetilde{\mu}_{t|t-1} + (1-\lambda)Y_{t}. $$
A more subtle comparison is out of steady state. For the Kalman filter, take $P_{1}=\infty$, then $K_1 = 1$. In our predictor, $$ \widetilde{\mu}_{t+1|t} = \left(1-\frac{1}{n_{\lambda,t}}\right) \widetilde{\mu}_{t|t-1} + \frac{1}{n_{\lambda,t}}Y_{t}. $$ We note that if $\lambda>0$, then $ n_{\lambda,t} = \frac{1-\lambda^t}{1-\lambda}. $ Hence it is interesting to plot $K_t \times n_{\lambda,t}$ against $t$, when $P_{1|0}=\infty$, the so called diffuse initial conditions for a Kalman filter (e.g. Chapter 1 of DurbinKoopman(12)). If this product is less than one, then our recursion gives more weight on $Y_t$ than the Kalman filter, and subsequently, less on the past.
Then $K_t n_{\lambda,t}=1$ if $t=1$. It is also 1 for every value of $t$ when $q=0$ for then $n_{\lambda,t}=t$. The limit is also 1 as $t \rightarrow \infty$ for $q>0$ (steady state), but the product is not 1 for every $t$ and $q>0$.
The easiest way to see differences is to think of $q$ as small, then, for example, $ \lambda \approx 1-q^{1/2}$, $ n_{\lambda,2} \approx 2 - q^{1/2}$, $ K_2 \approx \frac{1}{2} -q,$ and $K_2 n_{\lambda,2} \approx 1-q^{1/2}/2.$ Figure (ref) plots the product $K_t \times n_{\lambda,t}$ for $q\in \{0.001,0.1,0.3,1,2,10\}$ against $t$. The worst case is when $q$ is tiny, as it takes quite a large $t$ for the effect to disappear (as we move to the steady state).
Appendix (ref) compares our predictor to the widely applied class of score-driven filters Harvey(13), CrealKoopmanLucas(13), lange2024robust.
The analysis in this Section will be based on a working model $$ Y_t|Y_{1:t-1} \sim {\tt CEF}(\widetilde{\theta}_{t|t-1};h_t,\psi_t),\quad t=1,...,T, $$ where $\widetilde{\theta}_{t|t-1}$ comes from the exponentially weighted estimands. Later we will explicitly note the dependence of $\widetilde{\theta}_{t|t-1}$ on some hyperparameters, denoted $\omega$.
Under the working model, write $$ \widetilde{\mu}_{t|t-1} := \psi_t'(\widetilde{\theta }_{t|t-1}),\quad \widetilde{\Sigma}_{t|t-1} := \psi_t''(\widetilde{\theta }_{t|t-1}), $$ the model's conditional mean and variance, respectively.
The working model yields a working log-likelihood via the prediction decomposition, ignoring constants,
To make sense of this working log-likelihood function we need some regularity assumptions.
As the working log-likelihood $$l_t,\quad t=1,...,T,$$ comes from the ${\tt CEF}(\widetilde{\theta}_{t|t-1},h_t,\psi_t)$, under weak regularity conditions, it will also be a quasi-likelihood, so long as the predictive mean is modeled correctly. This property does not depend on the particular features of $\widetilde{\theta}_{t|t-1}$. We demonstrate this for the canonical exponential family time series case. It applies when $\widetilde{\theta}_{t|t-1}$ is formed from the weighted predictor.
Now make some assumptions about the data.
Using strict convexity of $\psi_t$ and Assumption (ref)(a) there exists a unique $\theta_{t|t-1}$ such that $\mu_{t|t-1} =\psi'_t(\theta_{t|t-1})$.
Under Assumptions (ref)(a) and (ref)(a), define the two sequences
then $\{M_t\}_{t=1}^T$ is a supermartingale with respect to the data's natural filtration. The $\{C_t\}_{t=1}^T$ is a previsible drift with increments $$ c_t = \psi'(\theta _{t|t-1})^{\tt T} \{\widetilde{\theta }_{t|t-1}-\theta_{t|t-1}\} - [\psi\{\widetilde{\theta }_{t|t-1}\} -\psi\{\theta _{t|t-1}\}] \le 0, $$ by convexity of $\psi$ from Definition (ref) and then using p.69 of BoydVandenberghe(04). The equality is only obtained iff $ \widetilde{\theta }_{t|t-1} = \theta _{t|t-1}, $ that is only if the working model has the correct predictive mean. Thus under Assumptions (ref)(a) and (ref)(a) the log-likelihood process $\{l_t\}_{t=1}^T$ can be viewed as a quasi-likelihood process.
If the magnitude of $\widetilde{\theta }_{t|t-1}-\theta_{t|t-1}$ is small, then $$c_t \approx -\frac{1}{2}(\widetilde{\theta }_{t|t-1}-\theta_{t|t-1})^{\tt T} \psi_t''(\theta_{t|t-1}) (\widetilde{\theta }_{t|t-1}-\theta_{t|t-1}).$$
More broadly, under Assumptions (ref)(a) and (ref)(a), the $\{M_t-C_t\}_{t=1}^T$ is a martingale. Further, under Assumptions (ref)(b) and (ref)(b) the angle bracket of the log-likelihood process
the angle-bracket process of the $\{M_t-C_t\}_{t=1}^T$ process. For all cases where $\langle M-C\rangle_T\rightarrow \infty$ then $(M_T-C_T)/\langle M-C\rangle_T \rightarrow 0$ almost surely, by the martingale strong law of large numbers.
Now turn to the maximum likelihood estimator of the hyperparameters based on the above exponential family quasi-likelihood, where the predictor comes from exponential weighting and uses a {\tt CEF} frame.
The leading version of this is given in Example (ref).
The resulting maximum likelihood estimate formed from the quasi-likelihood is: $$ \hat{\omega}_{\tt MLE} = \underset{\omega \in \Omega}{\arg }\max \ l_T(\omega), $$ noting explicitly how the quasi-likelihood is impacted by the choice of $\omega$ through $\{\widetilde{\theta}_{t|t-1}(\omega)\}_{t=1}^T$. Under Assumption (ref), the corresponding $\omega^*$ is the pseudo-true value of $\omega$ for this quasi-likelihood.
Throughout we will posit that $\widetilde{\theta}_{t|t-1}(\omega)$ is infinitely differentiable for all $\omega \in \Omega^*$ and write $$ \widetilde{\mu}'_{t|t-1}(\omega) = \frac{\partial \widetilde{\mu}_{t|t-1}(\omega)^{\tt T}}{\partial \omega}. $$
The corresponding score up to time $1 \le t\le T$ is
Assume that for all $\omega \in \Omega$ that $\widetilde{\Sigma}_{t|t-1}(\omega)=\psi''(\widetilde{\theta}_{t|t-1}(\omega))$ is invertible. Then $$ s_t(\omega) = \widetilde{\theta}'_{t|t-1}(\omega)\{h(y_t) -\widetilde{\mu}_{t|t-1}(\omega) \}, \text{where } \widetilde{\theta}'_{t|t-1}(\omega):= \frac{\partial \widetilde{\theta}_{t|t-1}(\omega)^{\tt T}}{\partial \omega} = \widetilde{\mu}'_{t|t-1}(\omega)\widetilde{\Sigma}^{-1}_{t|t-1}(\omega), $$ recalling $\widetilde{\Sigma}_{t|t-1}(\omega):= \psi''(\widetilde{\theta}_{t|t-1}(\omega))$. The Hessian is $$ H_t(\omega) = -\sum_{j=1}^t l_j''(\omega),\quad \text{where} \quad l_t''(\omega)=\frac{\partial^2 l_t(\omega)}{\partial \omega \partial \omega^{\tt T}}. $$
An alternative to the MLE is a two step procedure (e.g. NeweyMcFadden(94), EngleMezrich(96), FrancqHorvathZakoian(13)), which has the following structure.
Here the numerical optimization of the likelihood is only $2+\dim(\phi)$ dimensional. This estimation strategy can be very attractive when $\mathbb{E}[h(Y_1)]$ is high dimensional.
The two step procedure can be viewed as a method of moments estimator Pearson(1894), based around a random function $$ S_t(\omega;y_{{1:t}}) := \sum_{j=1}^t s_j(\omega;y_{1:j}),\quad s_t(\omega;y_{1:t}) := \left(
\right) $$ for $t=1,...,T,$ then $\hat{\omega}_{\tt 2Step}$ is a method of moments estimator which solves $S_T(\hat{\omega}_{\tt 2Step};Y_{1:T})=0,$ while $\mathbb{E}[s_t(\omega^*;Y_{1:t})]=0$ for each $t=1,...,T$. The corresponding Hessian is $$ H_T(\omega) := -\sum_{t=1}^T \frac{\partial s_t(\omega)}{\partial \omega^{\tt T}} = \sum_{t=1}^T \left(
\right). $$
We apply the exponentially weighted predictor and smoother to monthly data on household financial situation expectations from the University of Michigan Survey of Consumers, covering January 1978 to September 2025 ($T = 573$ observations). The survey asks respondents two questions: first, whether they are better off or worse off financially than a year ago, and second, whether they expect to be better off or worse off a year from now. The combination of responses yields seven mutually exclusive categories:
For each month, only the sample proportions (percentages) falling into each category are reported, not the underlying individual responses or sample sizes.\footnote{While monthly sample sizes are available, our framework would require within-month i.i.d. responses, which is unlikely given the clustered sampling design. We leave them out of our analysis.}
We model the observed proportions $y_t = (y_{1,t}, \ldots, y_{7,t})^{\tt T}$ as realizations from a time-varying seven category Dirichlet distribution $ Y_t \mid Y_{1:t-1} \sim \text{Dirichlet}(\widetilde{\theta}_{t|t-1}),$ with $ \widetilde{\theta}_{t|t-1} = (\widetilde{\theta}_{1,t|t-1}, \ldots, \widetilde{\theta}_{7,t|t-1})^{\tt T}, $ with concentration parameters $\widetilde{\theta}_{j, t|t-1} > 0$ for all $j \in \{1, \ldots, 7\}$ and all $t$. Each observation is treated as a single draw of a proportion vector. The predictor is based on a Dirichlet frame.
The exponentially weighted predictor and smoother are computed via the recursions in Example (ref) using the Dirichlet distribution's CEF parameterization given in Section (ref). Following the two-step quasi-likelihood procedure from Definition (ref), we estimate the centering parameter $\mathbb{E}[h(Y_1)]$ by the sample mean $\widehat{\mathbb{E}[h(Y_1)]} = T^{-1} \sum_{t=1}^T (\log y_{1,t}, \ldots, \log y_{7,t})^{\tt T}.$ Next, we fix $\widehat{\mathbb{E}[h(Y_1)]}$ and estimate $(\alpha, \lambda)$ as maximizing $\sum_{t=1}^T \log f(y_t; \widetilde{\theta}_{t|t-1}(\alpha, \lambda)).$
For the household financial situation data with seven categories, we obtain the parameter estimates $(\widehat{\mathbb{E}[h(Y_1)]}^{\tt T}, \widehat{\alpha}, \widehat{\lambda}) \approx (-1.76, -1.41, -1.78, -1.77, -2.73, -2.23, -3.53, 0.95, 0.64)$. The high anchoring estimate $\widehat{\alpha} = 0.95$ indicates little long-run effect, suggesting that shifts in household financial expectations, once they occur, tend to persist a long time. The discount parameter $\widehat{\lambda} =0.64$ implies a half-life of approximately $1.56$ months for the exponential weights. Following Remark (ref)(e), the steady state process is a vector {\tt ARMA}(1,1)-MD process, with autoregressive root $\widehat{\lambda}/\{1-\widehat{\alpha} (1-\widehat{\lambda})\}\approx 0.97$, and moving average root $-\widehat{\lambda} = -0.64$. Hence, the process exhibits substantial memory, with substantial root cancellation.
Figure (ref) presents the estimation results. Each panel shows the evolution of one expectation category over the 47-year period. The gray circles represent observed monthly proportions $y_{k,t}$. The light colored lines show the exponentially weighted predictor, the dark colored lines the exponentially weighted smoother. Both track major shifts in household financial expectations. There are the sharp declines in both “Continuous Increase” and “Intermittent Increase” during the 2008 financial crisis, and, more strikingly, since around 2018. This deterioration in optimism about household finances is mirrored by corresponding increases in negative expectations: “Intermittent Decline” and “Continuous Decline” roughly track their inverse patterns. The categories “Remain Unchanged” and “Don't Know/No Answer” stay relatively stable at 15-20% and 3%, respectively, throughout the sample. Appendix (ref) uses the empirical results to provide simulation based evidence for the performance of the hyperparameter estimators.
The exponentially weighted estimands combined with the {\tt CEF} delivers simple exact filters, predictors and smoothers. They are in the same statistical spirit as generalized linear models (e.g. McCullaghNelder(89)) and are relatively easy to fit to data using a quasi-likelihood approach either all at once or using a two step approach.
The exponentially weighted estimands framework can be used on non-exponential family models and even for general loss functions. In such cases the optimization has to be carried out numerically or through a sequential approximation.
We are grateful for the comments and encouragement from Dick van Dijk, Siem Jan Koopman, Rutger-Jan Lange and Bram van Os.\ Donker van Heel conducted this research during a visiting research appointment at Harvard University's Department of Statistics, supported by the Fulbright Program. The authors report there are no competing interests to declare.