Anticipating extreme losses using score-driven shape filters

2.1 Econometric models with score-driven shape parameter filters

The score-driven models of the daily log-return of the S&P 500 index y _t are formulated as:

where μ _t and exp(λ _t ) are the location and scale parameters, respectively. The error term ϵ _t is specified according to the following conditional distributions: (i) ϵ t | F t − 1 ∼ EGB 2 [ 0,1 , e x p ( ξ t ) , exp ( ζ t ) ] . If ξ _t = ζ _t then the probability distribution is symmetric, if ξ _t > ζ _t then the probability distribution is skewed to the right, and if ξ _t < ζ _t then the probability distribution is skewed to the left. Different values of ξ _t and ζ _t may imply tail-thickness on the left or right side of the probability distribution. For this conditional distribution all moments are finite, and the log of the conditional density of y _t is:

where Γ(x) is the gamma function. (ii) ϵ t | F t − 1 ∼ NIG [ 0,1 , e x p ( ν t ) , exp ( ν t ) tanh ( η t ) ] , where tanh(x) is the hyperbolic tangent function, and the absolute value of parameter exp(ν _t )tanh(η _t ) is less than parameter exp(ν _t ), as required for NIG. For this error term ν _t controls for tail-thickness. Moreover, if η _t = 0 then the probability distribution is symmetric, if η _t > 0 then the probability distribution is skewed to the right, and if η _t < 0 then the probability distribution is skewed to the left. For this conditional distribution all moments are finite, and the log of the conditional density of y _t is:

where K⁽¹⁾(x) is the modified Bessel function of the second kind of order 1. (iii) The third probability distribution is ϵ t | F t − 1 ∼ Skew – Gen – t [ 0,1 , t a n h ( τ t ) , exp ( ν t ) + 4 , exp ( η t ) ] , where the degrees of freedom parameter, exp(ν _t ) + 4, is higher than four. For this error term, ν _t controls for tail-thickness and η _t controls for the peakedness of the center of the probability distribution. Moreover, if τ _t = 0 then the probability distribution is symmetric, if τ _t > 0 then the probability distribution is skewed to the right, and if τ _t < 0 then the probability distribution is skewed to the left. For this conditional distribution the fourth moment is finite, and the log of the conditional density of y _t is:

where sgn(x) is the signum function. The density functions for the EGB2, NIG, and Skew-Gen-t distributions with alternative shape parameters are presented in Figure 1. Further technical details for each distribution are presented in Appendix A.

Figure 1:

Density functions for the EGB2, NIG, Skew-Gen-t (thick lines), and standard normal (thin lines) distributions.

For the shape parameters, we use the general notation ρ_k,t for k = 1, …, K, for each distribution. For example, EGB2 [0, 1, exp(ξ _t ), exp(ζ _t )] = EGB2 [0, 1, exp(ρ_1,t), exp(ρ_2,t)], where K = 2 is the number of shape parameters. In the specifications of the conditional distributions of ϵ _t , we use F t − 1 = [ μ 1 , λ 1 , ( ρ 1,1 , … , ρ K , 1 ) , ( y 1 , … , y t − 1 ) ] , i.e. we condition on the initial values of all filters.

In the following, the score-driven filters for μ _t , λ _t , and ρ_k,t are presented. First, μ _t is specified as:

where |ϕ| < 1 and u_μ,t is the scaled score function of the log-likelihood (LL) with respect to μ _t (Appendix A and Section 2.4). The score-driven filter for μ _t is named the first-order quasi-autoregressive QAR (1) model (Harvey 2013). Second, λ _t is specified as:

where |β| < 1, u_λ,t is the score function of the LL with respect to λ _t (Appendix A and Section 2.4). This specification measures leverage effects (Black 1976), by using parameter α* in the score-driven EGARCH model (Harvey 2013). The score-driven EGARCH models with constant shape parameters that use EGB2, NIG, and Skew-Gen-t distributions are named EGB2-EGARCH (Caivano and Harvey 2014), NIG-EGARCH (Blazsek, Ho, and Liu 2018), and Beta-Skew-Gen-t-EGARCH (Harvey and Lange 2017), respectively. Third, ρ_k,t is specified as:

where |γ _k | < 1, and u_ρ,k,t is the score function of the LL with respect to ρ_k,t (Appendix A and Section 2.4). For each distribution the constant shape parameter model, i.e. ρ_k,t = δ _k , is used as the benchmark, which we name the econometric tool of the existing method. Lags of exogenous explanatory variables may also be included in Eqs. (5)–(7), to extend the models of this paper.

Different ways of initialization are considered for the filters. For the results reported in the present paper, μ _t is initialized by using pre-sample data, λ _t by using parameter λ₀, and ρ_k,t by using its unconditional mean δ _k /(1 − γ _k ). Nevertheless, the results of this paper are robust to other ways of initialization. For example, the results for the case where parameters μ₀ and ρ_k,0 were used in this study for the initialization of μ _t and ρ_k,t are similar to the results reported in this paper.

Finally, we refer to the following two model specifications, which we studied as alternatives to the specifications of this section. First, we studied two-component λ _t filters, motivated by the work of Harvey and Lange (2018), in which two-component score-driven volatility models are suggested with time-invariant shape parameters. Second, we also studied multivariate filters involving the vector (μ _t , λ _t , ρ_1,t, …, ρ_K,t)′, motivated by the work of Blazsek, Escribano, and Licht (2022), in which multivariate score-driven filters are suggested for (μ _t , λ _t )′ with time-invariant shape parameters. For both alternatives, we find that the ML estimates are not reliable due to numerical problems in the maximization of the LL function. Hence, the results for those alternatives are not reported.

2.2 ML estimator

All models of this paper are estimated by using the ML method:

where Θ = (Θ₁, …, Θ _S )′ is the vector of time-invariant parameters, and ln ⁡ f ( y t | F t − 1 ; Θ ) for the EGB2, NIG, and Skew-Gen-t probability distributions is presented in Appendix A.

The T × S matrix of contributions to the gradient G (y₁, …, y _T ; Θ) is defined by its elements:

for period t = 1, …, T, and parameter i = 1, …, S (Wooldridge 1994). The tth row of G (y₁, …, y _T ; Θ) is denoted by using G _t (Θ), which is the score vector for the tth observation. The asymptotic covariance matrix of Θ ̂ is estimated by using the information matrix ( 1 / T ) ∑ t = 1 T G t ( Θ ̂ ) ′ G t ( Θ ̂ ) − 1 (Blasques et al. 2022; Creal, Koopman, and Lucas 2013; Harvey 2013). We use the following maintained assumptions:

1. Asymptotically, f ( y t | F t − 1 ; Θ 0 ) = p 0 ( y t | F t − 1 ; Θ 0 ) for Θ₀ from the parameter set Θ ̃ ⊂ I R S , where p₀ is the true conditional density, and Θ₀ represents the true values of Θ.

2. Asymptotically, f ( y t | F t − 1 ; Θ 0 ) is a dynamically complete density (Wooldridge 1994).

3. For all Θ ∈ Θ ̃ , |λ _t | < λ_max < ∞, and |ρ_k,t| < ρ_k,max < ∞

for all t and k, where λ_max ∈ IR⁺ and ρ_k, max ∈ IR⁺ do not depend on ϵ _t for all t. (S1) assumes an asymptotically correct model specification, which is studied empirically by using the specification tests of Section 2.3. (S2) assumes that the dynamics of the dependent variable are asymptotically correctly formulated for each period. One of the consequences of (S2) is that, asymptotically, G _t (Θ₀)′ is a MDS (Wooldridge 1994). (S3) assumes the uniform boundedness of the log-scale and all shape parameters, which implies finite moments for the score functions. (S2) and (S3) have consequences on the properties of the score functions of dynamic parameters (Section 2.4).

2.3 Specification tests

We empirically study assumption (S1) for the first four conditional moments for each score-driven model. The first four conditional moments of ϵ _t for the EGB2, NIG, and Skew-Gen-t distributions are reported in Appendix A. We define the following auxiliary error term as:

This transformation reduces the importance of those outliers that appear within ϵ _t . The robustness of model specification tests is increased when residuals are standardized according to Eq. (10) (see Li 2004, Chapter 4). The model specification test of the present paper uses the following properties:

The random variables within the expectations of Eqs. (11)–(14) are MDSs, which is tested by using the MDS test of Escanciano and Lobato (2009) for all model specifications of our paper. The rejection of the MDS null hypothesis of the Escanciano–Lobato test for any of the conditional moments of Eqs. (11)–(14) is evidence against the correct model specification assumption (S1).

2.4 Score functions

We present the asymptotic properties of the score functions u_μ,t, u_λ,t, and u_ρ,k,t at Θ = Θ₀.

First, we prove that all score functions, asymptotically at the true values of parameters, are MDSs. Score function u_μ,t, asymptotically at the true values of parameters, is a MDS due to the following arguments. Due to (S2), G _t  (Θ₀)′, asymptotically, is a MDS:

where index t − 1 indicates expectations that are conditional on F t − 1 . Since ∂ μ t / ∂ Θ ′ ≠ 0 ,

where k _t is defined in Appendix A for each distribution. Thus, E_t−1 (u_μ,t) = 0. Score function u_λ,t, asymptotically at the true values of parameters, is a MDS due to the following arguments. Due to (S2), asymptotically, G _t  (Θ₀)′ is a MDS:

Since ∂ λ t / ∂ Θ ′ ≠ 0 ,

Score functions u_ρ,k,t for k = 1, …, K, asymptotically at the true values of parameters, are MDSs because, due to (S2), G _t  (Θ₀)′, asymptotically, is a MDS:

Since ∂ ρ k , t / ∂ Θ ′ ≠ 0 ,

Hence, E(u_μ,t) = 0, E(u_λ,t) = 0, and E(u_ρ,k,t) = 0, asymptotically at the true values of parameters, due to the law of iterated expectations.

Second, we prove the finiteness of the second moments and covariances of the score functions and their derivatives, which is needed for the finiteness of the elements of the information matrix. The following results are true for all Θ ∈ Θ ̃ .

The score-driven models are robust to extreme observations, because the score functions in those models reduce the effects of outliers on the location, scale, and shape filters. Hence, outliers appear within the error term ϵ _t rather than within the score functions that update the dynamic equations. By using the S&P 500 dataset, the outlier transformation of the score functions is presented in Figure 2, which shows that extreme observations are never accentuated by the location, scale, and shape score functions. According to the figure, the transformations of ϵ _t asymptotically go to zero, are linear, or are in accordance with a slowly increasing function (e.g. the natural logarithm function ln(x)).

Figure 2:

Score functions for models with score-driven location, scale, and shape parameters. The score functions are presented as a function of ϵ _t ∈ [−250, 250], to show the asymptotic properties of the score functions for |ϵ _t | → ∞. In the score functions the dynamic parameters are replaced by their unconditional means.

In Appendix B, the derivatives of u_μ,t, u_λ,t, and u_ρ,k,t for k = 1, …, K, with respect to μ _t , λ _t , and ρ_k,t for k = 1, …, K are presented. Figures B1–B3 indicate that the derivatives of each score-function go to zero, are linear, or are slowly increasing functions as |ϵ _t | → ∞. In the illustrations of the score functions and their derivatives, the unconditional mean estimates of the score-driven variables are used, but the same functional forms of the score functions and their derivatives hold if the score-driven variables are replaced by their finite boundaries in accordance with (S3).

The score functions and their derivatives are nonlinear transformations of ϵ _t , for which the highest rate of increase is a linear increase, as |ϵ _t | → ∞ (Figures 2 and B1–B3).

This implies that, for the second moments and covariances of the score functions and their derivatives, the highest rate of increase of the squared score function or its derivative is a quadratic increase, as |ϵ _t | → ∞. For EGB2 and NIG all moments of ϵ _t are finite, and for Skew-Gen-t we assume that the fourth moment is finite. Hence, there exist positive numbers a and b that define the intervals − ∞ , − a and b , ∞ , respectively, for which the integrands of the fourth moment of ϵ _t bound the integrands of the second moments and covariances of all score functions and their derivatives. This implies that the variances and the covariances of all score functions and their second moments are finite for all score-driven models of this paper.

Third, scaled score function, asymptotically at the true values of parameters, u_μ,t is white noise, because u_μ,t, asymptotically, is a MDS and Var(u_μ,t) < ∞. Score function u_λ,t, asymptotically at the true values of parameters, is white noise, because u_λ,t, asymptotically, is a MDS and Var(u_λ,t) < ∞. For each k = 1, …, K, score function u_ρ,k,t, asymptotically at the true values of parameters, is white noise, because u_ρ,k,t, asymptotically, is a MDS and Var(u_ρ,k,t) < ∞. Hence, if |ϕ| < 1, |β| < 1, and |γ _k | < 1, then, asymptotically at the true values of parameters, μ _t , λ _t , and ρ_k,t are covariance stationary and are F -measurable functions of ϵ _s for all s < t.

2.5 MC simulation experiments

For all MC experiments, zero mean μ _t = 0, unit scale exp(λ _t ) = 1, and score-driven shape parameters are used for t = 1, …, T. Two sets of true parameter values are used: The first set assumes high persistence for the shape parameters (i.e. γ _k = 0.95 for all k). The second set assumes low persistence for the shape parameters (i.e. γ _k = 0.15 for all k). The true values of all parameters are presented in Table 1. By using those true values, 1,000 trajectories are simulated, and each trajectory includes T = 10,000 periods. By using the ML method, the parameters of the score-driven models with dynamic shape parameters are estimated for each trajectory.

We study the parameter estimates by using non-parametric methods. In Table 1, the 5, 50, and 95% quantiles of the 1,000 parameter estimates are reported. For the high-persistence case, the medians give a good approximation of the true values, and the 90% confidence intervals of the quantiles include all true values. For the low-persistence case, the medians give a good approximation of the majority of the true values; the only exceptions are some of the constant parameters, for which the true value is not within the 90% confidence interval (i.e. δ₂ for the EGB2 and NIG models). For the score-driven Skew-Gen-t model, for all parameters, the true value is within the 90% confidence interval for both the high- and low-persistence cases. The medians indicate that γ _k and κ _k are consistently estimated for all cases, supporting the use of the score-driven Skew-Gen-t model.

4.1 In-sample estimation and forecasting

The in-sample parameter estimates for all models for the period of February 14, 1990 to October 21, 2021 are presented in Tables C1–C4 of Appendix C. The in-sample parameter estimates for all models for the period of January 3, 2000 to October 21, 2021 are presented in Tables 3–6. For most of the specifications the ϕ and θ parameters which measure the expected return dynamics, are significantly different from zero. For all specifications, highly significant ω, α, α*, and β parameters are found for the scale. For most of the specifications, the dynamic parameters of shape (i.e. γ₁, γ₂, and γ₃), and the scaling parameter of the score function with respect to shape (i.e. κ₁, κ₂, and κ₃) are significant. All estimates of ϕ, β, γ₁, γ₂, and γ₃ are less than one in absolute value.

Moreover, we also report specification and statistical performance test results in Tables 3–6 and C1–C4. (i) For most of the specifications, the MDS null hypothesis of the Escanciano–Lobato test is not rejected up to the fourth moment, which supports the correct the model specification assumption (S1). (ii) In-sample statistical performances are compared by using the LR test. For many cases, we find that the performance of the score-driven specifications with dynamic shape parameters are superior to the performance of the score-driven model with constant shape parameters.

For the most general score-driven Skew-Gen-t model of this paper, the evolution of the shape parameters tanh(τ _t ), exp(ν _t ) + 4, and exp(η _t ), and the evolution of the scale parameter exp(λ _t ) are presented in Figure 4. The lowest values of exp(ν _t ) + 4 indicate the extreme events which significantly impacted the US stock market, and the circumstances of those events are described in Appendix D.

Figure 4:

Evolution of the shape and scale parameters for the score-driven Skew-Gen-t model with dynamic shape parameters τ _t , ν _t , and η _t . The lowest extreme values of ν _t indicate the extreme events (Appendix D).

In Tables 3–6, we also report in-sample volatility forecasting accuracy test results. As aforementioned, we use the 5 min realized volatility σ t * as a proxy of true volatility (Liu, Patton, and Sheppard 2015; Harvey and Lange 2018). The estimated conditional volatility for each model is denoted σ ( y t | F t − 1 ; Θ ) , for which closed-form formulas are reported in Appendix A, for all probability distributions of this paper. We use the loss functions SE i , t = σ t * − σ ( y t | F t − 1 ; Θ ) 2 and AE i , t = | σ t * − σ ( y t | F t − 1 ; Θ ) | , where i denotes the EGARCH specification for which the conditional volatility is σ ( y t | F t − 1 ; Θ ) . We note that as alternatives we also use SE and AE loss functions in which σ t * 2 and σ 2 ( y t | F t − 1 ; Θ ) are included; we obtain similar in-sample forecasting accuracy results to the ones reported in Tables 3–6.

Figure 5:

Log-returns on the S&P 500 (solid circles) and out-of-sample VaR (1 day, 99%) estimates for the score-driven Skew-Gen-t model are presented. All shape parameters are dynamic for the Skew-Gen-t dynamic-shape specification. After an extreme S&P 500 observation, for most of the cases, the VaR for dynamic shape (thick solid) predicts a more severe potential loss than the VaR for constant shape (thin solid). For all out-of-sample VaR estimates we use an expanding window, which starts at February 14, 1990 and ends on the trading day before each VaR estimate.

For each probability distribution, we compare the in-sample volatility forecasting performance of the score-driven model with constant shape parameters with that of the score-driven models with time-varying shape parameters. We define d _t (SE) = SE_c,t − SE_d,t, where i = c indicates the score-driven model with constant shape parameters and i = d indicates a score-driven model with time-varying shape parameters. We also define d _t (AE) = AE_c,t − AE_d,t, where i = c indicates the score-driven model with constant shape parameters and i = d indicates a score-driven model with time-varying shape parameters. We regress d _t (SE) and d _t (AE) on a constant parameter and we estimate its standard error by using the Newey–West heteroskedasticity and autocorrelation consistent (HAC) estimator (Newey and West 1987). A significantly positive estimate of the constant parameter indicates that the in-sample volatility forecasting performance of the score-driven model with dynamic shape parameters is superior to that of the score-driven model with constant shape parameters.

In Tables 3–6, we report the root mean SE (RMSE) and the mean AE (MAE) estimates for each specification, and for those statistics we also report in parentheses the p-values corresponding to the constant parameter in the linear regressions of the aforementioned loss function differences. We find two cases for the score-driven EGARCH models with dynamic shape parameters for which the in-sample volatility forecasting performance is significantly superior to that of the corresponding score-driven EGARCH models with time-invariant shape parameters: (i) score-driven Skew-Gen-t model with dynamic τ _t , dynamic ν _t , and dynamic η _t ; (ii) score-driven Skew-Gen-t model with constant τ _t , dynamic ν _t , and dynamic η _t . These specifications are also supported by the LR tests, hence we consider them as the best-performing score-driven EGARCH specifications of the present paper. We note that for the score-driven EGB2 and NIG specifications, the realized volatility-based model performance analysis does not indicate superior performances of the score-driven models with dynamic shape parameters.

4.2 Out-of-sample VaR backtesting

Financial institutions frequently evaluate and update their risk management systems. One of those evaluation methods is the backtesting of VaR models, by which the predictive performance of VaR models is tested by using out-of-sample forecasting and evaluation approaches.

In this section VaR backtesting applications are presented, for which the VaR measurements of the score-driven models with constant and dynamic shape parameters are compared. If the score-driven shape parameters predict tail-shape dynamics, then we will expect that the VaR predictive performance of the score-driven models with dynamic shape parameters is superior to the VaR predictive performance of the score-driven models with constant shape parameters. The results provide the following insight on the S&P 500 for practitioners about the quality of the VaR measurements for the new score-driven specifications. We find that the score-driven models with dynamic shape parameters anticipate better extreme losses than the score-driven models with constant shape parameters.

Extreme observations in the S&P 500 log-returns are concentrated during the period of the dot.com boom, 2008 US Financial Crisis, and COVID-19 pandemic, which are some of the most recent high-volatility periods presented in Appendix D. In Table 2, we present the dates of the periods of the dot.com boom, 2008 US Financial Crisis, and COVID-19 pandemic. Moreover, in Figure 3A and B, we present the concentration of the outliers for those subperiods of the sample. This motivates the consideration of those periods in the VaR backtesting applications. The design of the VaR backtesting procedure of this paper is in accordance with the framework of the Basel Committee (1996), in which a 1 day VaR is estimated out-of-sample at the 99% confidence level. In the present paper, a VaR (1 day, 99%) is estimated for each of the trading days of the backtesting period (Table 2).

An extending-window estimation approach is used for all score-driven models. For all VaR estimates, the first observation of the extending-window is February 14, 1990 and the last observation is for the last trading day before the day of the VaR estimate. VaR is approximated after the parameter estimation by using MC simulations, for which 10,000 possible log-returns are simulated for the trading day after the last observation of each rolling window. VaR (1 day, 99%) is defined by the 1% quantile of the log-return simulations. Motivated by the results presented in Section 4.1, the performance of VaR is compared for the score-driven Skew-Gen-t models with constant and dynamic shape parameters. All shape parameters are time-varying for the dynamic shape filters. See the VaR estimates in Figure 5.

To evaluate the VaR performance of different models, we use the Kupiec test to evaluate whether the proportion of VaR (1 day, 99%) failures is significantly higher than 1% during the backtesting period. The null hypothesis of the Kupiec test is that the observed VaR failure rate is equal to the failure rate suggested by the VaR confidence interval (i.e. 99% in this paper). The Kupiec test statistic is:

where p = 1%, T _b is the sample size for the backtesting period, and X is the number of VaR failures. Under the null hypothesis, the probability distribution of the Kupiec test statistic is the chi-squared distribution with 1 degrees of freedom.

The number of VaR failures and the Kupiec test results for the dot.com boom, 2008 US Financial Crisis, and COVID-19 pandemic backtesting periods are presented in Table 7. The results indicate a clear difference between the VaR forecasting accuracy of the score-driven EGARCH model with constant shape parameters and the score-driven EGARCH model with dynamic shape parameters: (i) For the score-driven EGARCH model with constant shape parameters, the null hypothesis of the Kupiec test is rejected at the 1% level of significance for the dot.com boom and it is rejected at the 10% level of significance for the 2008 US Financial Crisis and COVID-19 pandemic. (ii) For the score-driven EGARCH model with dynamic shape parameters, the same null hypothesis is never rejected. These results indicate for practitioners that financial institutions or investors can anticipate extreme losses better by using score-driven models with dynamic shape parameters than by using score-driven models with constant shape parameters.

We also performed the Christoffersen test. The null hypothesis of the Christoffersen test is that the arrival times of VaR failures are independent. If that null hypothesis is rejected, for example, due to consecutive VaR failures within the backtesting period, then the econometric model is not updated correctly after extreme observations. For all cases of Table 7 we find that the Christoffersen test does not provide evidence against the model specifications.

The VAR backtesting application indicates that, during periods of high market volatility, the VaR measurements of the score-driven models with dynamic shape parameters are improved, compared to the VaR measurements of the score-driven models with constant shape parameters.