Model Specification

The score-driven model is specified as in Creal et al. (2013) and is given by

where y_t is the time series observation at time t, for $t=1,\dots ,T$, $p_y\left(y_t|f_t,{\mathrm{ℱ}}_{t-1};\theta \right)$ denotes the predictive conditional density for observation y_t, f_t is the time-varying parameter for the conditional density, with${\mathrm{ℱ}}_{t-1}$denoting the information set based on the observations$y_1,\dots ,y_{t-1}$, and θ‎ is the parameter vector that contains fixed and unknown coefficients, including ω‎, α‎, and β‎, but also parameters that index the predictive density$p_y$. It is assumed that the observations become available sequentially over time. When the observation y_t arrives, the time-varying parameter f_t is updated using the updating equation $f_{t+1}=\omega +\alpha s_t+\beta f_t$. The innovation term for the updating equation is the scaled score $s_t=S_t\cdot {\nabla }_t$, where ∇_t denotes the score of the predictive conditional density $p_y\left(y_t|f_t,{\mathrm{ℱ}}_{t-1};\theta \right)$ with respect to f_t and $S_t\equiv S\left(f_t,{\mathrm{ℱ}}_{t-1};\theta \right)$ is a scaling function (discussed in the section “Possible Scaling Measures for the Score”). A similar formulation of the score-driven model was given by Harvey (2013).

The updating of f_t using the score of the conditional predictive density is intuitive. Namely, if a local ascent algorithm was used to find the local maximum of the predictive density at time t over f_t for a given parameter vector θ‎, then the score ∇_t indicates the direction in which the time-varying parameter must move, given the current position of f_t, according to the algorithm. For this reason, it is natural to use the score to determine the optimal updating step of f_t. This was underlined by Blasques et al. (2015), who proved the information-theoretic optimality of a simple class of score-driven models under certain regularity conditions. In particular, they showed that time-varying parameter updates based on the score always reduce the local Kullback-Leibler divergence in expectation and in every step. These results hold for misspecified models as well. Also, they argued that any parameter update that is not based on the score does not have this property. These are all limit results. However, Blasques et al. (2020) showed in a Monte Carlo study that the optimality also holds in finite samples for score-driven conditional volatility models.

Creal et al. (2008, 2013) discussed more options for generalizing the specification in Equation 1. For example, the model can be specified with a more extensive lag structure and with exogenous covariates $x_t$ in the updating equation for f_t. Also, the model can include other forms of nonlinearity in the updating equation, such as a regime-switching process and a long-memory process. Furthermore, a setting where the predictive density depends on past values of y_t and on some exogenous covariates$x_t$can be considered. Finally, y_t and f_t can be vectors, or even matrices, rather than scalars. If f_t is a time-varying parameter vector, the score$s_t$is a vector and the scaling factor$S_t$is a matrix. Having an observation vector y_t results in having a multivariate score-driven model. A range of extensions, including multivariate score-driven models, is discussed in Artemova et al. (2022).

Possible Scaling Measures for the Score

The scaling function$S_t$allows us to determine how the score ∇_t impacts the updating at time t, to obtain the next $f_{t+1}$. It must be emphasized that a different choice for$S_t$results in an inherently different model, with a different set of statistical properties. Hence, the most suitable scaling function for a specific setting must be determined with some care.

Typically, a natural choice for the scaling factor is a function of the variance of the score ∇_t. Then the scaling is based on the curvature of the conditional log density of the tth observation. More specifically, it is common to use the inverse asymptotic variance of the score, which is equal to the conditional information matrix under standard regularity conditions; so $S_t={\mathrm{ℐ}}_{t\mid t-1}^{-1}$ where ${\mathrm{ℐ}}_{t\mid t-1}={\mathbb{E}}_{t-1}\left[{\nabla }_t{\nabla }_t^{\prime }\right]$. This leads to the scaled score$s_t$having a variance of ${\mathrm{ℐ}}_{t\mid t-1}^{-1}$. It is also an intuitive scaling choice, because in this case the updating equation of f_t in Equation 1 can be interpreted as a Gauss-Newton algorithm for estimating f_t over time, as was pointed out by Creal et al. (2011). For this choice of $S_t$, well-known models become special cases of the score-driven model. For example, a score-driven volatility model $y_t=f_t{\epsilon }_t$, with error${\epsilon }_t$being standard normally distributed, reduces to the GARCH model of Bollerslev (1986) for this choice of scaling. Another example of a model that is encompassed by score-driven models for this choice of scaling is the multiplicative error (MEM) model of Engle and Gallo (2006), which in turn encompasses the autoregressive conditional duration and intensity (ACD and ACI) models of Engle and Russell (1998) and Russell (2001), respectively (see Creal et al., 2013, for a more extensive discussion).

Another possible scaling matrix is$S_t={\mathcal{J}}_{t\mid t-1}$, where ${\mathcal{J}}_{t\mid t-1}$is the square root of the

(pseudo-)inverse information matrix such that ${\mathcal{J}}_{t\mid t-1}^{\prime }\cdot {\mathcal{J}}_{t\mid t-1}={\mathrm{ℐ}}_{t\mid t-1}^{-1}$. This choice of $S_t$is particularly convenient, because it standardizes the score ∇_t such that$s_t$itself has a unit variance, which improves the tractability of the statistical properties of the model.

Finally, a typical last-resort option is to set $S_t=I$. For this choice of scaling, or essentially the lack of scaling, the statistical properties of the model will typically become more complicated. This is a clear disadvantage. However, when the information matrix ${\mathrm{ℐ}}_{t\mid t-1}$ is a constant (scalar), or it does not depend on the time-varying parameter, the scaling is less relevant, because the scaled score is scaled again through its multiplication by α‎ in the updating equation. Such cases include log-scale models of the form $y_t=exp\left(0.5f_t\right){\epsilon }_t$ and location models of the form $y_t=f_t+{\epsilon }_t$, where $S_t={\mathrm{ℐ}}_{t\mid t-1}^{-1}$ is proportional to $S_t=1$.

Estimation, Filter Invertibility, and Asymptotic Properties

As highlighted already, a convenient property of observation-driven models, and hence of score-driven models, is that an explicit expression of the log likelihood is available. Hence the estimation of the static parameter vector θ‎ via the method of maximum likelihood (ML) is straightforward. Consider a sample of T observations ${\left\{y_t\right\}}_{t=1}^T$ generated as described in Equation 1 under some true parameter${\theta }_0$. The true values of ${\left\{f_t\right\}}_{t=1}^T$ that are used to generate the observations are not observed. Therefore, a filtered sequence ${\left\{{\widehat{f}}_t\left(\theta \right)\right\}}_{t=1}^T$ is constructed recursively using the updating equation of f_t for some value of θ‎ and some starting value ${\widehat{f}}_1$. The application of the prediction error decomposition and the use of the filtered sequence ${\left\{{\widehat{f}}_t\left(\theta \right)\right\}}_{t=1}^T$ provide the maximization problem:

It is straightforward to evaluate the value of the log likelihood for some θ‎ because it requires only two steps: calculating the filtered time-varying parameter ${\left\{{\widehat{f}}_t\left(\theta \right)\right\}}_{t=1}^T$ via the score-driven updating in Equation 1, and calculating the log likelihood contribution of every $y_t$ given ${\widehat{f}}_t\left(\theta \right)$, for $t=1,\dots ,T$.

Before turning to the asymptotic properties of this ML estimator, an important property that must be examined when using observation-driven filters is filter invertibility. Invertibility essentially means that the filtered sequence ${\left\{{\widehat{f}}_t\left(\theta \right)\right\}}_{t=1}^T$ will “forget” its starting value in the limit. The starting value of the time-varying parameter f_t is unobserved. Hence, the true starting point $f_1$ is unknown and is not available for initializing the sequence ${\left\{{\widehat{f}}_t\left(\theta \right)\right\}}_{t=1}^T$. It needs to be replaced by some arbitrary value ${\widehat{f}}_1$. Sometimes it is possible to estimate the starting value alongside the other parameters, but typically this is not preferred, especially when f_t represents a high-dimensional vector. It is important to notice that when filter invertibility fails, the true path $f_t\left({\theta }_0\right)$ will not be retrieved in the limit, even if the true static parameter ${\theta }_0$ is known. Clearly, this will be problematic if the filtered values ${\widehat{f}}_t\left(\theta \right)$ are used to construct the log likelihood that is used for ML. This crucial point is sometimes overlooked by those who base their theoretical results on the implicit assumption that the value of $f_1$ is known exactly. Also, this assumption does not appear to be realistic because the rest of the sequence $f_2,f_3,\dots$ is assumed to be unobserved.

The concept of filter invertibility was discussed by Straumann and Mikosch (2006), who stressed that it is a condition needed for applicability of the quasi ML estimation of a general class of GARCH models. The importance of filter invertibility for the EGARCH model was highlighted by Wintenberger (2013). Blasques et al. (2018) provided a way to determine invertibility regions when no feasible analytical invertibility conditions on the parameters are available. Blasques et al. (2022) gave sufficient conditions for filter invertibility for score-driven models. As is usual in the stationary observation-driven literature, invertibility is established by showing that the filtered process converges almost surely to some unique stationary and ergodic limit process. They followed the approach of Straumann and Mikosch (2006), which adapted Bougerol (1993, Theorem 3.1) to obtain low-level conditions for filter invertibility for score filters and their first- and second-derivative processes. The latter results are useful for establishing asymptotic normality of the ML estimator.

After establishing the invertibility results for the score-driven model, Blasques et al. (2022) further provided sufficient conditions for consistency and asymptotic normality of the ML estimator ${\widehat{\theta }}_T$. Two conditions are crucial in this development. First, the updating equation of the “true” f_t must be contracting on average to ensure that the true time-varying parameter is stationary and ergodic. Second, the filtering equation of ${\widehat{f}}_t$ must be uniformly contracting to ensure filter invertibility. The asymptotic properties derived by Blasques et al. (2022) are global and they rely on low-level conditions in terms of “building blocks” of score-driven models, which are represented by the equations provided in Equation 1. For example, the derivatives of the score with respect to the parameters must have bounded moments up to some specific order. These results are particularly helpful for researchers who want to establish theoretical properties of the ML estimator for a specific score-driven model. Given certain conditions for the observations, the asymptotic properties remain applicable under potential model misspecification. These theoretical results do not trivially generalize to settings where the observations y_t and/or the time-varying parameter f_t are no longer univariate. However, specific asymptotic results can be established for particular multivariate score-driven models (see the discussions in Blasques et al., 2016, and Bazzi et al., 2017).

Location Model Specification

When the aim is to filter the conditional mean $f_t=\mathbb{E}\left(y_t|{\mathrm{ℱ}}_{t-1}\right)$ of a (univariate) sample of observed data ${\left\{y_t\right\}}_{t=1}^T$, consider the score-driven filtering model as given by

where ${\left\{{\epsilon }_t\right\}}_{t\in \mathrm{\mathbb{Z}}}$ is assumed to be an independent and identically distributed (i.i.d.) random sequence with mean 0 and probability density function$p_{\epsilon }$, and where f_t is updated according to Equation 1. Assume that ${\epsilon }_t$is Gaussian with variance ${\sigma }^2$. Then $S_t={\mathrm{ℐ}}_{t\mid t-1}^{-1}={\sigma }^2$, leads to a scaled score $s_t=y_t-f_t$ which yields the linear updating equation for f_t as given by

The updating for f_t reduces to an exponentially weighted moving average (EWMA) recursion when coefficients have values ω‎ = 0, β‎ = 1, and 0 < α‎ < 1. When α‎ > 0, the updating of the conditional expectation f_t is intuitive, since f_t is effectively the one-step-ahead forecast of y_t and $y_t-f_t$ is the corresponding forecast error. Hence, the updating equation, Equation 3, takes into account the forecast error to construct $f_{t+1}$, which is defined as the one-step-ahead forecast of the next observation $y_{t+1}$. Specifically, for α‎ > 0, if f_t has a lower (higher) value than y_t, then the scaled error $\alpha \left(y_t-f_t\right)$ will ensure that $f_{t+1}$ increases (decreases) compared to f_t.

Substituting Equation 2 into this equation shows that ${\left\{f_t\right\}}_{t\in \mathrm{\mathbb{Z}}}$ follows an autoregressive process of order 1, an AR(1) process, with autoregressive coefficient β‎. Hence, a necessary and sufficient condition for stationarity of this sequence is that $\mid \beta \mid <1$. The updating equation, Equation 3, implies that f_t is a weighted average of all past observations, where observation $y_{t-j}$ has weight $\alpha {\left(\beta -\alpha \right)}^{j-1}$, for $j=1,\dots ,t-1$. It follows immediately that the filtered sequence ${\widehat{f}}_t$ is stationary in the limit, if and only if $\mid \beta -\alpha \mid <1$. This is also the condition for filter invertibility.

Finally, by substituting Equation 3 into Equation 2, it follows that y_t is implicitly generated by an ARMA(1,1) process as given by

with autoregressive coefficient β‎ and moving average coefficient α‎ – β‎. This is an interesting result and it applies to any distribution for ${\epsilon }_t$. Notice that the process reduces to an AR(1) process when α‎ = β‎ , and to an MA(1) process when β‎ = 0. Higher-order ARMA processes are obtained when higher-order lags for f_t and $s_t=y_t-f_t$ are considered for the updating equation $f_{t+1}$ in Equation 1 or, more specifically, in Equation 3.

Robust Filtering

The properties of the score-driven location model are almost identical to those of a parameter-driven Gaussian (stationary) signal plus noise model, which is obtained by replacing $s_t$with a Gaussian random error sequence in the updating Equation 3. As shown by Harvey and Luati (2014), this parallel vanishes when considering a non-Gaussian predictive conditional density $p_y\left(y_t|f_t,{\mathrm{ℱ}}_{t-1};\theta \right)$in Equation 1. For instance, consider a fat-tailed conditional density such as the Student’s t density. The data-generation process for this model will lead to many observations that would be referred to as “outliers” in a Gaussian context. The updating of $f_{t+1}=\omega +\alpha \left(y_t-f_t\right)+\beta f_t$ will not work well, because large observations from previous times will be incorporated in the filtered level f_t: it will take time for the outlying observations to “work through the system.” Hence, this dynamic model is not robust to large observations that are induced by a heavy-tailed conditional distribution. Therefore, it is preferred to adopt a model that accounts for these fat tails. Score-driven models are designed to do this by considering fat-tailed densities for $p_y\left(y_t|f_t,{\mathrm{ℱ}}_{t-1};\theta \right)$in Equation 1. A more elaborate discussion was provided by Caivano et al. (2016, Section 2.3 and 2.4), who demonstrated that the conditional score reflects the tail shape of distributions and how this connects to robustness.

For instance, let ${\epsilon }_t$in Equation 2 be Student’s t distributed with v > 0 degrees of freedom and scale parameter $\sigma >0$. Then the conditional observation distribution becomes

implying that the score-driven updating equation will become

according to Equation 1, using $S_t={\left(1+{\nu }^{-1}\right)}^{-1}{\sigma }^2$, which is proportional to the inverse of the conditional information matrix ${\mathrm{ℐ}}_{t\mid t-1}=\left(\nu +1\right){\sigma }^{-2}/\left(\nu +3\right)$. Equations 2 and 4 together form the score-driven Student’s t model proposed by Harvey (2013, Section 3.1) and Harvey and Luati (2014). Notice that the degrees of freedom parameter v is only required to be positive, which is why the model is referred to as a location model and not as a mean model. If the degrees of freedom $\nu \to \infty$, the update becomes identical to the Gaussian update of Equation 3, which is not surprising because in that case the Student’s t distribution approaches a normal distribution. The score is plotted for various different choices of v in Figure 1. For finite values of v, it is clear that the update in Equation 4 downweights large values of $y_t-f_t$, and this downweighting is more severe for smaller values of v. In that way, the updating takes into account the heavy tails of the disturbances ${\epsilon }_t$, because it implies that large prediction errors might be due to the nature of the innovations. As $\mid y_t-f_t\mid$ goes to infinity, the score $s_t$ goes to 0. Because of the redescending nature of the score, Caivano and Harvey (2014) referred to the weighting induced by this score-driven filter as a parametric form of trimming. In the robustness literature, trimming is a common technique in which observations above some threshold receive a weight of 0. Hence, the Student’s t location filter essentially induces a soft form of trimming.

Open in new tab

Figure 1. Student’s t score plots for the location model. The score s_t of the Student’s t location model in Equation 4 for scale $\sigma =1$ and different degrees of freedom v.

Harvey and Luati (2014) provided asymptotic results for the ML estimator of this model, including an explicit asymptotic variance matrix in terms of the model’s parameters, but invertibility of the filter was not explicitly discussed. Blasques et al. (2022) have formulated conditions for consistency and asymptotic normality of the ML estimator of this model while taking account of the invertibility conditions. Blasques et al. (2018) developed a weaker version of the invertibility condition for this model. It is also possible to adopt a nonstationary version of the score-driven local level model as given by

where $s_t$is defined as in Equation 4; this model was suggested by Harvey and Luati (2014). It is a special case of the model in Equations 2 and 4 with ω‎ = 0 and β‎ = 1. This updating equation for f_t also implies an EWMA scheme but now with time-varying weights that account for the shape of the distribution (see Caivano et al., 2016). Furthermore, this model can be extended to a more general unobserved components model with score-driven time-varying trend and seasonal components. Such specifications can also be extended with explanatory variables and fixed effects (see Harvey & Luati, 2014).

Leptokurtosis and Asymmetry

Whereas the score-driven location filter based on Student’s t innovations induces a form of trimming, the filter of a score-driven location model based on the exponential generalized beta distribution of the second kind (EGB2), as proposed by Caivano and Harvey (2014) and Caivano et al. (2016), has a Winsorizing property. Winsorizing is a well-known approach for the treatment of large observations in the robustness literature. It entails an updating equation for $f_{t+1}$ where $s_t$is set to a constant value when the observation passes some given threshold (see Maronna et al., 2006, Chapter 2).

The EGB2 distribution allows for asymmetry and leptokurtosis, but unlike the Student’s t distribution, it has exponential tails instead of heavy tails. If ${\epsilon }_t$in Equation 2 is EGB2 distributed with mean 0, variance ${\sigma }^2>0$ and nonnegative shape parameters $\xi$ and $\varsigma$, then

where $\mathrm{\Delta }=\psi \left(\xi \right)-\psi \left(\varsigma \right)$, $h=\sqrt{{\psi }^{\prime }\left(\xi \right)+{\psi }^{\prime }\left(\varsigma \right)}$, $\psi \left(\cdot \right)$ is the digamma function, ${\psi }^{\prime }\left(\cdot \right)$ is the trigamma function, and$B\left(\cdot ,\cdot \right)$is the beta function (see, for instance, Wang et al., 2001, for further details). The distribution is symmetric if $\xi =\varsigma$, positively skewed if $\xi >\varsigma$, and negatively skewed if $\xi <\varsigma$. The kurtosis decreases as $\xi$ and $\varsigma$ increase and the distribution encompasses the normal distribution (if $\xi =\varsigma \to \infty$) and the Laplace distribution (if $\xi =\varsigma =0$). The scaled score $s_t$in Equation 1, corresponding to the predictive density in Equation 5, is given by

where the scale is set to $S_t={\sigma }^2$, which is proportional to the inverse of the conditional information matrix. It is clear that the fraction in the scaled score $s_t$is uniformly bounded between $-\mathit{h\sigma \xi }$ and $\mathit{h\sigma \varsigma }$, and $s_t$converges to these values as $y_t-f_t\to -\infty$ and $y_t-f_t\to \infty$, respectively. Hence, this updating mechanism is effectively subject to a gentle or soft form of Winsorizing. The contribution of large observations is bounded, but it is not redescending, as for Student’s t distributed errors; see Figure 2 for a visual representation of this.

Open in new tab

Figure 2. Plots of the score ∇_t corresponding to the location model for different distributions. The unscaled score ∇_t corresponding to the location model (Equation 2) for the Gaussian distribution (with ${\sigma }^2=1$), the Student’s t distribution (with ${\sigma }^2=6/8$ and v = 8) and the symmetric EGB2 distribution (with ${\sigma }^2=1$ and $\xi =\varsigma \approx 0.7689$), which all have variance 1 and excess kurtosis 1.5 (except for the Gaussian distribution).

When $\xi \ne \varsigma$, the distribution is skewed and the score update becomes asymmetric. For example, if the distribution is negatively skewed, so if $\xi <\varsigma$, then large positive values get a larger weight than large negative values of the same magnitude. This is intuitive because negative skewness implies that negative spikes are more likely to occur than positive spikes, and therefore the robust filter should be less sensitive to negative spikes.

Naturally, score-driven location models based on many other distributions can also be considered. This will lead to location models with different properties, because the score will take into account the shape of the distribution. Take, for example, a general error distribution (GED), also known as the exponential power distribution, with parameter v (see Harvey, 2013, Section 3.10). The GED encompasses the normal distribution (for v = 2) and the Laplace distribution (for v = 1). This GED is more peaked in comparison to the EGB2 distribution and has super-exponential tails for v > 1. This leads to the score being unbounded, unlike the score of the EGB2 distribution, but it will diverge at a lower rate than the score of the normal distribution if v < 2. Other distributions to consider for robust filtering are, for example, the generalized t distribution of McDonald and Newey (1988) and skewed versions of all aforementioned symmetric distributions, which can for instance be obtained by using the approach of Fernández and Steel (1998). The section “Score-Driven Scale Models” considers these distributions in the context of score-driven scale models.

Illustration: Treasury Bill Rate Spreads

An empirical illustration is presented to show how the score-driven location models can perform in practice.¹ The data set under consideration consists of quarterly observations of the difference between the 3-month and the 6-month treasury bill (T-bill) rates. The sample ranges from 1959:Q3 until 2021:Q1.² The T-bill rate for a certain maturity is the yield received by investors for T-bills of that particular maturity in the secondary market. The difference between rates of different maturities is referred to as the spread. The left panel of Figure 3 shows a plot of the resulting T-bill rate spreads. This time series is rather noisy, because it has some sudden and temporary spikes and drops, especially in the early 1980s. This seems to indicate the need for a robust filter if the goal is to filter some underlying location parameter.

Open in new tab

Figure 3. Data of T-bill rate spreads between 3 and 6 months maturity with the corresponding filtered location f_t of the normal, Student’s t, and EGB2 models (left), and the corresponding estimated response curves (right).

To accommodate these concerns, score-driven models with Student’s t and EGB2 innovations are considered, together with a score-driven model with Gaussian innovations for comparison. The corresponding updating functions are provided in Equations 4, 6, and 3, respectively. The static parameters are estimated by the method of ML. It is anticipated that the robust Student’s t and EGB2 models will be the most suitable here, because of the occasional erratic behavior of the data. The importance of robustness is confirmed by the estimation results reported in Table 1. Both the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) are the lowest for the model with Student’s t innovations and highest for the model with the Gaussian innovations. The EGB2 model has worse AIC and BIC values than the Student’s t model, but they are both better than those of the Gaussian model. These results appear to suggest that the Student’s t distribution, which has fat tails, fits the data better than the EGB2 distribution, which has exponential tails.

Flexible distributions like the EGB2 distribution can lead to optimization problems during parameter estimation, especially for this illustration, where the sample size is small. For this illustration, it appears that likelihood optimization for the EGB2 model enters into a saddlepoint, even when the process is restarted with new starting values. Hence, no standard errors are reported in this case. It appears that this numerical issue is caused by the joint estimation of $\xi$ and$\varsigma$. Therefore, Table 1 also reports the estimation results for the symmetric EGB2 model, which is obtained through the restriction $\xi =\varsigma$. For this model, we do not encounter these numerical issues during the estimation process. This fitted model has a slightly worse AIC, but a slightly better BIC than its unrestricted version.

The estimated degrees of freedom parameter of the Student’s t distribution is somewhat small. This estimation result implies that the error distribution is heavy-tailed, which in turn causes the score s_t to become more robust against outliers. This is clearly visible in the plotted response curve of the models presented in the right panel of Figure 3. For small to moderate values, the Student’s t filter reacts stronger than the Gaussian filter, but for larger values this is no longer the case. Indeed, the response of the former filter redescends back to 0 in the limit, whereas that of the latter diverges to infinity. For the EGB2 models, the score follows a path similar to that of the Student’s t, but it is monotonically increasing instead of redescending. The unrestricted EGB2 model has a score function that is asymmetric, since $\xi >\varsigma$ implies a mildly positively skewed distribution (the skewness is$0.33$). The left panel of Figure 3 presents the filtered location corresponding to the fitted models, where the filtered path of the symmetric EGB2 model is omitted because of its similarity to that of the asymmetric EGB2 model. It is clear that the Gaussian filter reacts more to large observations than the Student’s t and the EGB2 filters. The filtered paths of the Student’s t and EGB2 models are very similar, but the latter tends to have a slightly stronger reaction to large shocks than the former.

Univariate Scale Models

For a univariate (demeaned) observation y_t, the conditional univariate scale is denoted by f_t and the basic score-driven scale model is formulated as

where ${\left\{{\epsilon }_t\right\}}_{t\in \mathrm{\mathbb{Z}}}$ is an i.i.d. sequence with mean 0 such that $\mathbb{E}\left(y_t|{\mathrm{ℱ}}_{t-1}\right)=0$. When the variance of ${\epsilon }_t$is equal to 1, the model (Equation 7) reduces to a volatility model where f_t is equal to the conditional variance, that is,$\mathrm{Var}\left(y_t|{\mathrm{ℱ}}_{t-1}\right)=f_t$. However, this restriction is not necessary; the variance of ${\epsilon }_t$ could even be infinite. Instead of using demeaned observations, it is straightforward to extend the model to allow for a non-zero mean$\mu$and for autoregressive or ARMA dynamics in the observation equation. To keep this treatment simple, assume that observations$y_t$are generated by the model in Equation 7.

A basic illustration is obtained by taking ${\epsilon }_t\sim \mathcal{N}\left(0,1\right)$ for every t. If the conditional variance f_t is updated according to the score-driven framework in Equation 1, with scaling factor $S_t={\mathrm{ℐ}}_{t\mid t-1}^{-1}=2f_t^2$, the updating equation becomes

where ω‎ > 0, α‎ ≥ 0, and β‎ ≥ α‎ to ensure positivity of f_t. It is not hard to see that the resulting model is equivalent to a regular GARCH model (Bollerslev, 1986) with parameters α‎ and β‎ ˗ α‎. There are certain parameter restrictions that impose stationarity and filter invertibility (Straumann & Mikosch, 2006). If instead a Student’s t distribution with v > 0 degrees of freedom is considered, without changing the updating equation of f_t, the resulting model is the GARCH-$t$ model of Bollerslev (1987). However, if a score-driven updating equation is used, the resulting model departs from the regular GARCH framework, because then:

where ω‎ > 0, α‎ ≥ 0, and β‎ ≥ α‎, to impose positivity, and where the used scaling factor $S_t=2f_t^2$ is proportional to the inverse of conditional information matrix (see Table 2). Notice that this model is not based on a standardized Student’s t distribution and that v > 2 is not required, which is why this model is referred to as a scale model rather than a volatility model. Using a standardized Student’s t distribution with unit variance is also possible (see Creal et al., 2013, for an example). Similar to the Student’s t location model in Equations 2 and 4, it is clear that for finite values of v, large values of $y_t^2$ are downweighted, unlike in the Gaussian scale model. For $\nu \to \infty$, the updating function becomes identical to Equation 8, which also follows from the fact that the Student’s t distribution reduces to a Gaussian distribution in that case. Blasques et al. (2022) gave (parameter) restrictions for filter invertibility, consistency, and asymptotic normality of the MLE for this model in their main example.

Univariate Log-Scale Models

A point of concern for scale models is that filtered conditional volatilities should be strictly positive. In score models, the positivity of the conditional volatility can be typically ensured by imposing appropriate parameter restrictions. To avoid having to impose positivity restrictions on the parameters, it may be convenient to choose a different parametrization, where the logarithm of the scale is modeled instead of the scale itself. An advantage of this choice of parametrization is also that the stationarity condition is relatively straightforward, as it will simply be $\mid \beta \mid <1$. So instead of modeling the scale $f_t={\sigma }_t^2$ in $y_t={\sigma }_t{\epsilon }_t$, the log scale $f_t=log\left({\sigma }_t^2\right)$ is modeled, which gives the log scale model:

where${\epsilon }_t$is again i.i.d. with mean 0 and where f_t is updated according to Equation 1 based on the distribution of the innovations. So if${\epsilon }_t$has unit variance, then $exp\left(f_t\right)=\mathrm{Var}\left(y_t|{\mathrm{ℱ}}_{t-1}\right)$. Creal et al. (2013) stated that the score-driven log scale model is equivalent to the well-known EGARCH model of Nelson (1991) if ${\epsilon }_t$has an asymmetric Laplace distribution. For a selection of other distributions, Table 2 reports the score ∇_t and the conditional information matrix corresponding to the log scale parametrization. As was pointed out, the conditional information matrix does not depend on f_t in this case, so the choice of scaling is less critical here. For example, using the Student’s t distribution leads to:

for a scaling factor equal to$S_t=2$(see Creal et al., 2011, 2013; Harvey, 2013, Chapter 4). In the latter reference, the model is referred to as the Beta-$t$-EGARCH model since the fraction in the score function has a Beta distribution, if evaluated at the true parameter. This property makes it a theoretically appealing model because the asymptotic properties of the ML estimator of the parameter vector in this model can be derived in a convenient manner (see Harvey, 2013, Section 4.2). This model has been adopted in many empirical studies in the literature, and it has been shown to outperform the regular GARCH and GARCH-t models in most of the studies (see, for example, Harvey & Sucarrat, 2014, Blazsek et al., 2016, and Catania & Nonejad, 2020).

Table 2 also gives the score of the EGB2 log scale model. The corresponding score-driven model is discussed briefly in Caivano and Harvey (2014), for a slightly different parametrization and with the main focus on a model where both the location and scale are score-driven.

Figure 4 plots the news impact curves (NIC) of score-driven log scale models with a normal, Student’s t, and EGB2 distribution with the same variance. The NIC is calculated by evaluating the score ∇_t as a function of the standardized observations$y_t/exp\left(f_t/2\right)$. The bounded score of the Student’s t distribution reflects the fact that it is a fat-tailed distribution and therefore large observations do not imply a proportionally large increase in the underlying volatility process. On the contrary, the score of the EGB2 distribution does diverge, which reflects the exponential tails of the distribution. This divergence is linear instead of a quadratic, which reflects that it has excess kurtosis. It is interesting to compare the shape of these NICs to those of the location model that were plotted in Figure 2.

Open in new tab

Figure 4. News impact curves of different score-driven log scale models. News impact curves— i.e., the plots of the score ∇_t as a function of the standardized observations $y_t/exp\left(f_t/2\right)$— corresponding to the log scale model (Equation 10) for the Gaussian distribution (with ${\sigma }^2=7/5$), the Student’s t distribution (with v = 7), and the symmetric EGB2 distribution (with $\xi =\varsigma \approx 1.874$), all with variance 7/5.

A more flexible alternative to the Student’s t distribution is the generalized t distribution proposed by McDonald and Newey (1988). Harvey and Lange (2017) considered a score-driven scale model based on this distribution. The corresponding predictive density for two positive shape parameters v and h is given in Table 2. For h = 2, this is a regular Student’s t distribution with v degrees of freedom, while a GED distribution with parameter h is obtained if $\nu \to \infty$. The GED with v = 2 is the normal distribution, whereas v = 1 results in a Laplace distribution. Table 2 presents the predictive density, score, and information matrix for a log scale model under the GED. The updating of the log scale parameter f_t of the score-driven model based on the generalized t distribution is given by:

using$S_t=2$. So it is clear that this update has the same robustness properties as the Student’s t scale filter in Equation 11 as long as $\nu <\infty$, because in that case the score is bounded. Harvey and Lange (2017) provided explicit asymptotic results for the distribution of the ML estimator. They also proposed likelihood ratio (LR) and Lagrange multiplier (LM) tests for the null hypothesis of thin tails (that is, $\nu \to \infty$) against the alternative hypothesis of fat tails (that is, $\nu <\infty$).

Univariate Log-Scale Models with Skewness

Another way to introduce flexibility into the score-driven model is to use a skewed error distribution. For instance, Harvey and Sucarrat (2014) considered constructing score-driven scale models based on distributions skewed by the method of Fernández and Steel (1998). Any probability density function$p_{\epsilon }\left(\epsilon \right)$ that is unimodal and symmetric around 0 can be converted to a skewed density function as follows:

where $0<\gamma <\infty$ is the skewing parameter, and for $\gamma =1$, $\gamma <1$, and $\gamma >1$, the distribution is symmetric, left skewed, and right skewed, respectively. If this skewing method is applied to the Student’s t distribution, the score is given by

It follows that skewing the distribution directly induces asymmetry in the score and therefore in the update of the log scale parameter f_t. Figure 5 presents a plot of the NIC for some different choices of $\gamma$. It is clear that for $\gamma <1$, negative values are downweighted more, whereas positive values are downweighted less than for the symmetric distribution. Hence, under negative skewness, the score update takes into account that positive values are relatively less common than negative values of the same magnitude and should therefore receive a relatively higher weight. The asymptotic results of the symmetric score-driven Student’s t log scale model generalize straightforwardly to this skewed variant. The same skewing approach can be used for other distributions, such as the GED.

It is important to realize that because of the skewness, the expectation of the innovations is no longer equal to 0, which can be solved by reformulating the observation equation of the model as

where${\mu }_{\epsilon }=\mathbb{E}\left[{\epsilon }_t\right]$. The score and therefore the updating equation of f_t must be altered accordingly. A more detailed discussion is provided in the study by Harvey and Sucarrat (2014, Section 4); in their empirical study, they showed that the score-driven skewed t model generally outperforms competing models, including alternative skewed models, such as the GJR GARCH model (Glosten et al., 1993) with skewed innovations and the normal mixture GARCH model (Alexander & Lazar, 2006).

Open in new tab

Figure 5. News impact curves of score-driven log scale model for the skewed t distribution. News impact curves—i.e., the plots of the score ∇_t as a function of the standardized observations $y_t/exp\left(f_t/2\right)$—corresponding to the log scale model (Equation 10) for the skewed Student’s t distribution with v = 6 degrees of freedom and for different values of skew parameter $\gamma$.

There are also alternative ways to skew a distribution. For example, Harvey and Lange (2017) considered a score-driven scale model based on a skewed generalized t distribution, using the skewing method of Zhu and Galbraith (2010), which uses a slightly different parametrization than that of Fernández and Steel (1998). Furthermore, not only skewness is introduced, but also the shape parameters v and h are allowed to have a different value above and below 0. So apart from skewness, a more explicit form of asymmetry is introduced to the distribution as well. Harvey and Lange (2017) demonstrated that the resulting score-driven model works well in practice, although it is recommended to simplify it by testing different parameter restrictions using LR, Wald, or LM tests.

Volatility in Mean

Multiple extensions of score-driven scale models have been proposed in the literature. As an illustrative example, consider the volatility-in-mean model proposed in Harvey and Lange (2018). This model alters the log specification of the ARCH in mean (ARCH-M) model of Engle et al. (1987) to have score-driven updating of the standard deviation. So instead of the regular log scale specification in Equation 10, now the time-varying scale parameter$exp\left(0.5f_t\right)$ also occurs in the mean equation of the observations$y_t$:

where the innovation sequence ${\left\{{\epsilon }_t\right\}}_{t\in \mathrm{\mathbb{Z}}}$ is i.i.d. with mean 0 and fixed variance, and where the dynamics of f_t are again score-driven as in Equation 1. Here μ‎ and λ‎ are fixed but unknown coefficients that determine the level and the ARCH-M effect, respectively. The rationale for including the conditional scale parameter in the mean equation is that volatility of stock returns tends to be positively correlated with the level of the returns, because the equity risk premium increases if uncertainty rises.

When the innovations${\epsilon }_t$are Student’s t distributed with v degrees of freedom, Harvey and Lange (2018) showed that

If λ‎ = 0 (and μ‎ = 0), the regular log scale update of Equation 11 is recovered. In practice, the last term of the score has little impact on the filter, because λ‎ is typically very small. Harvey and Lange (2018) derived theoretical properties for the resulting model, including moments and asymptotic results of the ML estimator.

Illustration: Electricity Spot Price Returns

To demonstrate the robustness property of score-driven scale models with heavy-tailed distributions, this section compares the Gaussian and Student’s t volatility model for a time series of electricity spot price returns.³ The data under consideration is a sequence of daily returns of the PJM electricity market, which serves 13 states in the United States. The data⁴ span from April 6, 2008, to December 31, 2013 (2,096 observations). Electricity prices are known to show occasional extreme behavior in the form of sudden positive spikes, after which they quickly move back to the pre-spike level. The primary cause of this is the nonstorability of electricity

(see Escribano et al., 2011). In effect, the returns based on these prices also display these incidental “outliers.” Figure 6 clearly shows that these properties are observable for the prices and returns currently under consideration. Models that are robust to sporadic extreme observations, such as the Student’s t score-driven model, are thus most likely to be valuable in modeling conditional volatility dynamics of these data.

Open in new tab

Figure 6. Electricity prices and returns. The average daily electricity price in the PJM market from April 6, 2008, to December 31, 2013 (left panel) and the corresponding returns—i.e., $\left(p_t-p_{t-1}\right)/p_{t-1}$ where p_t denotes the price—scaled by a factor of ten (right panel).

Consider the log scale specification in Equation 10 with a non-zero mean μ‎: $y_t=\mu +exp\left(\frac{1}{2}{f}_t\right){\epsilon }_t.$ The scaling factor is chosen to be equal to the inverse of the conditional information matrix, $S_t={\mathrm{ℐ}}_{t\mid t-1}^{-1}$. Table 2 reports the score and information matrix of the Gaussian and Student’s t distributions that were used to construct the models.

The resulting ML estimates are reported in Table 3 alongside the corresponding likelihood-based criteria. According to both reported information criteria, AIC and BIC, the Student’s t model has a better in-sample log likelihood value than the Gaussian model, taking into account that the former model has one extra parameter. The estimates of α‎ and β‎ are similar for both models, but because the estimated degrees of freedom parameter v is small, the two resulting models are inherently different.

Open in new tab

Figure 7. News impact curves of estimated models log scale models. The news impact curves— i.e., the estimated response $\alpha s_t$ as a function of the standardized observations $y_t/exp\left(f_t/2\right)$—corresponding to the Gaussian and Student’s t score-driven volatility models fitted to PJM electricity spot price returns.

Figure 7 shows the NICs of the estimated models. The estimated degrees of freedom parameter v is low, so the estimated error density has heavy tails, and therefore the score update of the Student’s t model is robust to large values. This is visible in the figure, as the response value corresponding to this model is bounded by a constant value. For the Gaussian model, this is not the case, because the response increases at a quadratic rate as y_t grows. Figure 8 shows the filtered volatility corresponding to the two fitted models alongside the absolute returns. As the NICs indicate, the filtered volatilities of the Gaussian model are impacted more by large observations than those of the Student’s t model are. After a spike in the returns, it takes a while for the Gaussian filtered volatility to move back to the pre-spike level. The filtered volatilities of the Student’s t model also react to large observations, but in a less extreme manner. A spike in the returns does not indicate a large change in the underlying volatility process, given the fact that prices quickly move back to their original level. Therefore, the robust volatility filter of the score-driven Student’s t model seems to be more suitable than the nonrobust filter of the Gaussian model, since the former is less sensitive to these temporary increases in the returns.

Open in new tab

Figure 8. Filtered volatilities. The absolute PJM electricity price returns and the corresponding filtered volatilities of the Gaussian model—i.e., $exp\left(\frac{1}{2}{f}_t\right)$—and Student’s t model—i.e., $\sqrt{\nu /\nu -2}exp\left(\frac{1}{2}{f}_t\right)$.