Measurement of Risk Based on QR-GARCH-EVT Model

Assume that the value of a financial position a moment t is V_t, the loss in the next phase 1 is: (1)Lt(l)=Vt+1−Vt{L_t}(l) = {V_{t + 1}} - {V_t}

Assume that the accumulative distribution function of random variable is L_t(l) under confidence level p, this paper defines VaR_1−p as the maximum loos potentially faced with: (2)VaR1−p=inf{x|FL(x)≥1−p}{VaR}_{1 - p} = \inf \{x|{F_L}(x) \ge 1 - p\}

(2) can be transformed to: (3)∫−∞VaRf(x)dx=1−p\int_{- \infty}^{VaR} f(x)dx = 1 - p

For the univariate function, VaR is the loss distribution function q = 1 − p fractile, accordingly: (4)Pr[Lt(l)>VaR1−p]≤p{P_r}\left[ {{L_t}\left(l \right) > {VaR}_{1 - p}} \right] \le p

Equation (4) can be transformed to: (5)∫VaR∞f(x)dx=p\int_{VaR}^\infty f(x)dx = p

So, assume that the mean value and variance of financial asset return are μ_t and σ_t, the fractile corresponding to confidence level p is Q_1−p, and the general expression of VaR is: (6)VaR1−p=μt+Q1−pσt{VaR}_{1 - p} = {\mu _t} + Q_{1 - p}^{{\sigma _t}}

Despite of simple definition, it is a challenging statistical question to calculate VaR. Statistically, VaR is the fractile of the expected return of an investment portfolio. So, VaR is closely correlative to the distribution of return. However, it is not enough to only assume the distribution function of the return. The reason is that the return does not meet the hypothesis of independent different distribution. Xiao and Koenker [22] proposed a fractile regression based QR-GARCH model: (7){rt=σtεt,εt∼iid(0,1)σt=ω+αxt−1+βσt−12xt=λ+δσt2+η1εt+η2εt−12+ut,εt∼iid(0,1)\left\{{\matrix{{{r_t} = {\sigma _t}{\varepsilon _t},{\varepsilon _t} \sim iid(0,1)} \hfill \cr {{\sigma _t} = \omega + \alpha {x_{t - 1}} + \beta \sigma _{t - 1}^2} \hfill \cr {{x_t} = \lambda + \delta \sigma _t^2 + {\eta _1}{\varepsilon _t} + {\eta _2}\varepsilon _{t - 1}^2 + {u_t},{\varepsilon _t} \sim iid(0,1)} \hfill \cr}} \right.

Where ω is a constant term and greater than 0, α represents the speed of the response to market impact made by the rate of return, β represents the amplitude of vitality compared with earlier stage, both α and β are greater than 0 and α + β < 1. Using iterative method to give the estimated value of σ_t: (8)σt=ω1−β+α(xt−1+βxt−2+β2xt−3+β3xt−4+…){\sigma _t} = {\omega \over {1 - \beta}} + \alpha \left({{x_{t - 1}} + \beta {x_{t - 2}} + {\beta ^2}{x_{t - 3}} + {\beta ^3}{x_{t - 4}} + \ldots} \right)

If the ARCH (∞) process in equation (8) is stable, and this paper can obtain: (9){σt2=α0+∑j=1∞αjxt−jrt2=(α0+∑j=1∞αjxt−j)εt2\left\{{\matrix{{\sigma _t^2 = {\alpha _0} + \sum\nolimits_{j = 1}^\infty {\alpha _j}{x_{t - j}}} \hfill \cr {r_t^2 = \left({{\alpha _0} + \sum\nolimits_{j = 1}^\infty {\alpha _j}{x_{t - j}}} \right)\varepsilon _t^2} \hfill \cr}} \right.

In equation (9), α0=ω1−β{\alpha _0} = {\omega \over {1 - \beta}} , α_j = β^j−1.

Let F_t−1 be the information at moment t-1, then the fractile function of rt2r_t^2 based on past information can be expressed as: (10)Qrt2(τ|Ft−1)=α0(τ)+∑j=1∞αjxt−j{Q_{r_t^2}}\left({\tau |{F_{t - 1}}} \right) = {\alpha _0}\left(\tau \right) + \sum\nolimits_{j = 1}^\infty {\alpha _j}{x_{t - j}}

Where, a_j(τ) = a_jQ_ε² (τ), j = 0,1,2,···, Q_e² (τ) are the τ fractile of εt2\varepsilon _t^2 .

Due to unobservability of fluctuation process, when there is no enough data, the estimation of extreme fractile regression is not accurate, so Xiao and Koenker [22] method shall be used to estimate (ω⌢,α⌢1,β⌢1)'{(\mathop \omega^\frown,{{\mathop \alpha^\frown }_1},{{\mathop \beta ^\frown} _1})^{'}} and the specific process is as follows:

First, according to equations (9) and (10) use fractile regression to estimate γ^(τk)\hat \gamma ({\tau _k}) of, specific equation is: (11)γ^(τk)=argmin∑t=m+2T[ρτ(rt2−a0(τ)−∑t=j+1Tajxt−j]\hat \gamma ({\tau _k}) = \arg \min \sum\limits_{t = m + 2}^T [{\rho _\tau}(r_t^2 - {a_0}(\tau) - \sum\limits_{t = j + 1}^T {a_j}{x_{t - j}}]

Use equation (11) to obtain a˜(τk)\tilde a({\tau _k}) , k = 1,2,···, k, ρ is loss function. Then use minimum distance method to estimate [22] to obtain (a^0,a^1,⋯a^m)'{({\hat a_0},{\hat a_1}, \cdots {\hat a_m})^{'}} , at the same time, the volatility σ_t can be estimated equation (9). Further, in the paper fractile regression of r_t is conducted based on (1,ρ^t−12,xt)'{(1,\hat \rho _{t - 1}^2,{x_t})^{'}} : (12)θ^(τ)=argmin∑t=m+2T[ρτrt2−ω(τ)−a(τ)xt−j]+β(τ)ρ^t−12\hat \theta (\tau) = \arg \min \sum\limits_{t = m + 2}^T [{\rho _\tau}r_t^2 - \omega (\tau) - a(\tau){x_{t - j}}] + \beta (\tau)\hat \rho _{t - 1}^2

Under (τ1,τ2,⋯,τk)'{({\tau _1},{\tau _2}, \cdots,{\tau _k})^{'}} fractile, it can be estimated that, (13)θ^(τk)=(ω^τk,α^τk,β^τk)'\widehat \theta \left({{\tau _k}} \right) = (\widehat \omega {\tau _k},\widehat \alpha {\tau _k},\widehat \beta {\tau _k}{)^{'}}

Finally, the evaluated value of (ω⌢,α⌢,β⌢)'{(\mathop \omega \limits^\frown,\mathop \alpha \limits^\frown,\mathop \beta \limits^\frown)^{'}} can be obtained, and the residuum term ε_t, of rate of return can be obtained.

Use QR-GARCH model to obtain residuum series, it is needed to transform ε_t into standard residuum term Z_t, let μ_t,σ_t be the condition mean value and condition variance of return series residuum term, then: (14)(Zt−n+1,⋅⋅⋅,Zt)=(rt−n+1−μ^t−n+1σ^t−n+1,⋅⋅⋅,rt−μ^tσ^t)({Z_{t - n + 1,}} \cdot \cdot \cdot,{Z_t}) = ({{{r_{t - n + 1}} - {{\hat \mu}_{t - n + 1}}} \over {{{\hat \sigma}_{t - n + 1}}}}, \cdot \cdot \cdot,{{{r_t} - {{\hat \mu}_t}} \over {{{\hat \sigma}_t}}})

Through the above equation the VaR value of r_t can be estimated through the VaR value of Z_t and the dynamic VaR of asset return r_t under confidence level p, which can be expressed as VaR1−ptVaR_{1 - p}^t , the calculation formula is: (15)VaR1−pt=μ+σtVaR(z)1−pVaR_{1 - p}^t = \mu + {\sigma _t}VaR{(z)_{1 - p}}

Where, μ is expected return, σ_t represents the fluctuation estimation of the day, VaR1−ptVaR_{1 - p}^t represents the risk value of residuum term Z_t when the fractile is p, equation (15) represents the fractile level VaR1−ptVaR_{1 - p}^t of standard extreme residuum series Z_t under p. During the process of solving, it is also the process of POT model to build a model for Z_t. The key to POT modelling is to determine the threshold, and the typical methods to determine the threshold includes the excess expectation function graph and Du Mouchel 10% principle [23]. In this paper, both methods this paper used to determine the threshold, so the function distribution of the standard residuum series of the part beyond threshold is: (16)Fu(y)=P(r−u≤y|r>u)=F(u+y)−F(u)1−F(u){F_u}\left(y \right) = P\left({r - u \le y|r > u} \right) = {{F(u + y) - F(u)} \over {1 - F(u)}}

When the threshold is great enough, according to extreme theory [24, 25], there is F_u(y) ≈ G_{ξ, σ} (y): (17)Fu(y)≈Gξ,σ(y)=1−(1−yξσ)−1ξξ≠0{F_u}\left(y \right) \approx {G_{\xi,\sigma}}(y) = 1 - {\left({1 - {{y\xi} \over \sigma}} \right)^{- {1 \over \xi}}}\xi \ne 0

G_{ξ, σ} (y) function is GPD distribution, Where, σ and ξ represent scale parameter and shape parameter, the greater ξ is, it represents that the fatter the tail is, and the values of σ and ξ can be obtained by using maximum likelihood estimation. Under the condition where F_u(y) ≈ G_{ξ, σ} (y), by combining equation (16) and (17) this paper can obtain: (18)F(y)=F(u)+Gξ,σ(y−u)[1−F(u)]F\left(y \right) = F\left(u \right) + {G_{\xi,\sigma}}(y - u)\left[ {1 - F\left(u \right)} \right]

In order for research convenience, this paper uses N to represent the total number of samples, use n to represent the number of samples beyond the threshold, the threshold distribution function F (u) is represented by (N − n)/N, putting it in equation (18), and combining the evaluated values of σ and ξ, we can obtain: (19)F^(y)=1−nN(1+ξ^y−uσ^)−1/ξ^\widehat F\left(y \right) = 1 - {n \over N}{\left({1 + \widehat \xi {{y - u} \over {\widehat \sigma}}} \right)^{- 1/\widehat \xi}}

By combining equation (6), it can be concluded that the risk value under confidence level p is: (20)VaR1−pt=F^−1(y)=u+σ^ξ^{[Nn(1−p)]−ξ^−1}VaR_{1 - p}^t = {\widehat F^{- 1}}\left(y \right) = u + {{\widehat \sigma} \over {\widehat \xi}}\left\{{{{\left[ {{N \over n}\left({1 - p} \right)} \right]}^{- \widehat \xi}} - 1} \right\}

For given confidence level p, by using equation (6), this paper can obtain VaR(Z)_1−P, then put VaR(Z)_1−P in equation (15), dynamic VaR value of the model to be solved can be obtained. So it can obtain the dynamic VaR model based on QR-GARCH-EVT: (21){rt=σtεt,εt∼iid(0,1)σt=ω+αxt−1+βσt−12xt=λ+δσt2+η1εt+η2εt−12+ut,εt∼iid(0,1)VaR1−pt=F^−1(y)=u+σ^ξ^{[Nn(1−p)−ξ^−1]}\left\{{\matrix{{{r_t} = {\sigma _t}{\varepsilon _t},{\varepsilon _t} \sim iid(0,1)} \hfill \cr {{\sigma _t} = \omega + \alpha {x_{t - 1}} + \beta \sigma _{t - 1}^2} \hfill \cr {{x_t} = \lambda + \delta \sigma _t^2 + {\eta _1}{\varepsilon _t} + {\eta _2}\varepsilon _{t - 1}^2 + {u_t},{\varepsilon _t} \sim iid(0,1)} \hfill \cr {VaR_{1 - p}^t = {{\widehat F}^{- 1}}(y) = u + {{\widehat \sigma} \over {\widehat \xi}}\left\{{\left[ {{N \over n}{{\left({1 - p} \right)}^{- \widehat \xi}} - 1} \right]} \right\}} \hfill \cr}} \right.

In order to test the model constructed previously, this paper selected HS300 index (HS300). Considering that HS300 index was launched beginning from April 4^th, 2005, for analysis convenience, this paper selected HS300 Index, the time span of Hang Seng Index sample is: from April 8^th, 2005 through March 19^th, 2017, totally nearly 12 years, for each index there are about 2,767 sample points. The source of the data is Tonghuashun software, the return per day is defined as: X_t = 100(ln p_t − ln p_t−1).

From Table 1 it can be found that all sample series are left-skewed form (skewed degree <0), in addition, they have leptokurtosis characteristics (kurtosis >3), and in J-B test, under 1% significance level every sample series is significantly different from normal distribution. In order to further test the distribution characteristics of sample series, in the paper, the frequency chart and QQ chart (Fig. 2) of the sample series were drawn, the upper tail and lower tail of HS300 series deviate from normal distribution apparently, showing thick tail characteristics, consequently it can be concluded that the sample series show typical “leptokurtosis” characteristics. In general, the sample series show asymmetrical distribution, with typical financial return “leptokurtosis” characteristics, and a certain burst and gathering characteristics. Ljung-Box statistic Q(10) shows that under 1% significance level, HS300 has no series correlation, ARCH-LM test indicates that HS300 return has ARCH effect, and ADF test shows that HS300 return is steady, So QR-GARCH-EVT model can be used to build model.

Fig. 1

Fig. 2

Considering excellent properties of QR-GARCH model, in the paper it used above QR-GARCH model to fit HS300 return series, for specific result please see Table 2.

From Table 2 it can be found that α + β < 1 indicates the stability of QR-GARCH model, the overall risk of HS300 Index return is apparently positively correlative to its volatility in the past, and the return series has a continuous fluctuation. A and η₁ are opposite in sign, indicating there is leverage effect, and negative impact has a greater influence that positive impact, β is 0.637, indicating that the volatility in QR-GARCH model is not only influenced by the fluctuation in previous stage.

Based on QR-GARCH model estimation, this paper obtained the residuum series of HS300 Index return. Standardizing the residuum series, the test showed that the standard residuum was steady, without self-correlation, So EVT model could be used to build a model for it.

A model was built for standard residuum series in three phases, including model exploration phase, model fitting and model diagnosis phase. In the first model exploration phase, the key to determine the threshold, currently popular threshold selection methods include excess expectation function graph and Du Mouchel 10% principle, etc., but there has been on a unified standard yet. Considering the importance of threshold selection attached to the model, this paper will combine both methods to determine the magnitude of the threshold.

Firstly, according to Du Mouchel 10% principle, through calculation, this paper has obtained that the threshold of HS300 series was 1.165, and made the mean life graph of HS300 Index return series to aid in selection of threshold. By observing its mean residual life, it can be found that the series began to have positive ramp near above 1, by combining the threshold obtained by Du Mouchel 10% principle, so this paper has concluded that the threshold was 1.569 respectively.

Fig. 3

Secondly, in model fitting phase, based on the obtained threshold, this paper used EVT model to estimate standard residuum series to obtain the scale parameter σ and shape parameter ξ of EVT mode, the result is shown in Table 2. From the table it can be found that the shape parameter ξ is greater than 0, indicating that extreme standard residuum series has apparent fat tail.

Finally, in model testing phase. In order to test the effectiveness of EVT model estimation, this paper made the fitting diagnosis chart of standard residuum series, as shown in Fig. 4, most points directly fell in the threshold distribution chart and tail distribution chart or nearby, only individual points deviate, without influencing fitting effect. Over fitting effect of EVT is good.

Fig. 4

Based on estimation of QR-GARCH-EVT model, this paper estimated risk value according to the risk measure model constructed previously. In Fig. 5 this paper lists the VaR values when the fractile is between 0.01 and 0.05.

Fig. 5

In order to test the accuracy of the model, this paper used Kupiec test (posterior test) proposed by Kupiec [26]. Kupiec test is a unconditional coverage test, which statistically check if the frequency of excess in the sample is enough statistically close to the selected confidence level. For example, for a confidence level of 99%, it is estimated that one excess occurs every 100 day on average. For the confidence level of 95%, at the same time, interval five excesses occur. For the confidence of other levels, the same situation occurs. Firstly, a binary variable is defined, it is 1 when VaR prediction fails; otherwise it is 0. During the back test the typical value of the variable is recorded, and the significance of unconditional coverage was tested using the statistical of likelihood ratio test, then the likelihood ratio test statistical is as follows:

Where, p is probability level, n is the number of tested samples and m is time numbers of failures. In theory, under given confidence level, the less the failure is, the higher the accuracy of the model is, the better the effect is. At the same time, for comparison, this paper calculated the results of Kupiec test of GARCH model, GARCH-EVT model as well as QR-GRACH model, for specific result please see Table 4.

From Table 4, through comparison of the model, it found that under 95% confidence level, the VaR estimated only by GARCH type model is different slightly, with the risk of underestimation, while combination of GARCH type and EVT has an apparently good effect. But compared with QR-GARCH-EVT model, in terms of failure rate, the latter has a higher accuracy.

Under 99% confidence level, the number of failures of every model is apparently greater than the theoretical level of 28. Underestimating ample risk to some extent, the model has an apparent inadaptability, with the low effectiveness of risk prediction. In terms of model comparison, QR-GARCH-EVT model still performs better than other models.

So, under 95% confidence level, the QR-GARCH-EVT model constructed in this paper has higher model accuracy, stronger effectiveness, but under 99% confidence level, the model has inadaptability. In general, QR-GARCH-EVT model performs better than other models.