Metodoloski zvezki, Vol. 15, No. 1, 2018, 59-72 Gumbel GARCH Model with Stock Application Mehrnaz Mohammadpour1 Fatemeh Ziaeenejad2 Abstract The paper proposes a new GARCH model with Gumbel conditional distribution. Several statistical properties of the model are established, like autocorrelation function and stationarity. We consider two methods for estimating the unknown parameters of the model and investigate properties of the estimators. The performances of the estimators are checked by a simulation study. We investigate the application of the process using a real stock data. 1 Introduction The generalized autoregressive conditionally heteroscedastic (GARCH) model has been found to be useful in many economic and financial studies which captures the tendency for volatility clustering (Bollerslev, 1986). In the classical GARCH model, normal distribution has been considered as a conditional distribution which is not quite logical consideration in the most financial studies. Modeling of GARCH time series was first introduced by Bollerslev (1986) based on the normal conditional distribution. Among the GARCH models, we cite the standardized t-student model (Bollerslev, 1986), the normal poisson mixture model (Jorion, 1988), the power exponential model (Baillie and Bollerslev, 1989), the normal-log normal mixture model (Hsieh, 1988), the generalized exponential model (Nelson, 1990b), the normal model (Nelson, 1990a, 1992), the threshold GARCH model (Glosten et al., 1993) and the stable GARCH model (Liu and Brorsen, 1995 and Calzolari et al., 2014). In this paper, a new GARCH model with Gumbel conditional distribution is introduced. The Gumbel GARCH model is justified by the need to model extreme observations more realistically than would be possible using the standard normal GARCH. The Gumbel model is the traditional model in extreme value analysis which has the same status as the normal model in other applications. The major advantage of the Gumbel model is that the distribution can be specified by location and scale parameters as in the Gaussian case. In the following, we briefly investigate Gumbel distribution and its features. The Gumbel distribution is referred to the distribution corresponding to extremes. The Gumbel distribution with the location parameter a and the scale parameter y (denoted by 1 Faculty of Mathematical Sciences, University of Mazandaran, Bbolsar, Iran; m.mohammadpour@umz.ac.ir 2 Faculty of Mathematical Sciences, University of Mazandaran, Bbolsar, Iran; ziaeene-jadsh@yahoo.com 60 Mohammadpour and Ziaeenejad Gumbel(a, y)) has the probability density function /» / \ 1 { x a , x a, fx(x) = - exp--exp(-) Y V Y Y The mean and variance of the distribution are p = a + vy, and 2 1 2 2 a = 7 n y , 6 where v is Euler-Mascheroni constant. In this work, we consider a = 0. The rest of the paper is organized as follows. In Section 2, we construct the Gum-bel GARCH model and investigate the stationarity condition and some properties of the model. Section 3 deals with the estimation of the model parameters by the Yule-Walker and maximum likelihood methods. The performances of the estimators are checked by a small Monte Carlo simulation. An application of the model for the stock data is given in Section 4. 2 Gumbel GARCH Let {Xt}teZ , Z the set of integers, be a discrete time second-order process and Ft-1 is a a field generated by {Xs}s 0, ai > 0, ft > 0, i = 1,... ,p, j = 1,..., q, p > 1, q > 0. The conditional probability mass function of {Xt} has the following form f (xt|Ft_i) = — exp( — )exp(-exp(-—)), _ Yt Yt Yt where Yt = at2) 1. n2 The conditional mean is 6 2\1 E (Xt |Ft-i) = vYt = v (—at)2. n2 In Proposition 1, we establish a necessary condition on the parameters of the model to ensure that the process is a second-order stationary process. Gumbel GARCH Model. 61 Proposition 1. For a second-order stationary process {Xt}teZ to satisfy (2.1) and (2.2), it is necessary that Y^,P=i ai(1 + v2 ) + E q=1 < 1. Proof. Let a2 = E(of). Since P ^¿) + ^T jE(at2-j) ¿=i j=i a2 = E(o2) = ao + J] E(X^) + J] jE(at2_j) we obtain p 6 q = ao + ^ "¿(1 + v2-2 )o2 + jo2, ¿=i j=i o2 = ?-a0-,. (2.3) (1 - EP=i ai(1 + v2n2) - E?=i j In (2.3), the parameters must necessarily satisfy the condition p 6 9 1 - T + v2^) - T jj > 0' ¿=i j=i which completes the proof. □ The following theorem presents a necessary and sufficient condition for the second order stationarity of the model. For the simplicity of notation, we assume that p > q. Theorem 1. A necessary and sufficient condition for the process {Xt} to be second-order stationary is that all roots of q p 1 - Aao - T(Aai + - T Aaizi = 0, (2.4) ¿=i ¿=q+i lie inside the unit circle, where A = (1 + v26/n2). Proof. Let Yi,t = E(XtXt-i), i = 1, 2,... ,p. The conditional second moment is obtained as E (Xt2|Ft_i) = Var(Xt|Ft_i) + E 2(X |Ft_i) = + v 2(6/n2at2) = Aa2. Then A ao + £ "¿E(X^) + J] jE(at2_j) Yo,t = E (X2) = E (E (Xt2|Ft_i) = E (Aat2) p q :t2_i) + T & E (a2_j) i=i j=i pq Aao + ^ Aai7o,t_i + ^ j 7o,t_j i=i j=i qp Aao + J^(Aai + A)Yo,t_i + Aai7o,t_i. ¿=i ¿=q+i 62 Mohammadpour and Ziaeenejad The necessary and sufficient condition for a non-homogeneous difference equation to have a stable solution is that all roots z\,... ,zp of eqn (2.4) lie inside the unit circle, (Goldberg, 1958). □ The following theorem gives the autocorrelation function (ACF) of {Xt2} which is used for Yule-Walker estimation. Theorem 2. Suppose that {Xt} following the Gumbel GARCH(p, q) model is second-order stationary. Then YX(k) := Cov(Xt2, Xt2_k) and y22 (k) := Cov(af, a2-k) satisfy the following equations p Y2(k) = £(l + V26/n2)aiY2(|i - k|) i=1 min(k_1,q) q + J & y2 (k - j) + ¿(1 + v26/n2m y22 (| k - j | ),k > 1, (2.5) j=i j=i and min(k,p) Y22(k) = E (1 + v26/n2)«iy22(|i - k|) i=1 + E ai/(1 + v26/n2)Y2(k - j) + J(|k - j|),k > 0. (2.6) i=k+1 j=1 Proof. See Appendix A for details. □ Corollary 1. Suppose that {Xt} following the Gumbel ARCH(p) model is second-order stationary. Then the autocovariance function y^O satisfies the following equation Yx2(k) = JE(1 + v26/n2)aiY2(|k - i|), k > 1. (2.7) i=1 The equations of Corollary 1 are obviously nearly identical to the Yule Walker equations of the standard AR(p) model. As a consequence, the model of order p can be identified with the help of the partial autocorrelation function (PACF). 3 Estimation and Simulation Comparison In this section, we will investigate two methods for parameter estimation of the Gumbel GARCH(p, q) model based on a realization X1,..., Xn of the process. These estimators are compared via Monte Carlo simulations in terms of their means and standard deviations. Gumbel GARCH Model. 63 3.1 Yule-Walker Estimator Let Y2 (k) = n (—2 — X2 j ^Xt2+fc — X2 j , 0 < k < n, be the sample autocovariance __n function of {—t2}, where X2 = nE —2 is the sample mean. The Yule-Walker (YW) n t=1 estimators of the model are obtained by substituting the sample yX (k) in eqn (2.7) and solving them. Example 1. Consider the Gumbel ARCH (2) as Xt|Ft-1 ~ Gumbel(0,Yft), 2 Var(Xt |Ft_i) = (n2^2 = at2 = a + ^ alX2t_l. i=1 The explicit ^YW estimators of a1 and a2 are PX(1)(1 — Px(2)) Q!1 a2 (1 — pX(l))(l + v 26/n2) Px(2) — (pX(1))2 (1 — (PX(1))2)2(1 + v2 6/n2) where pX(k) = ^. Also note that a0 can be estimated from eqn (2.3). Yx v0/ 3.2 Maximum Likelihood Here we derive the maximum likelihood estimator (MLE) of the unknown parameter B where B = [a0, a1,..., ap, p1, @2,..., ]. The MLE of the parameter is obtained by maximization of the conditional log-likelihood function n n X X £(B) = £ 4(B) = logYt + - — exp(—-)], ^ ^ yt Yt t=1 t=1 where and t 6 2n 1 Yt = ^ 2 : n2 ^t = ao + ai-t_i + Y1 # i=1 j=1 Solving the system of equations ^^ = 0, the MLE of B is obtained. This can be done by using standard nonlinear maximization procedures which may be found in most of the statistical and data analysis packages. 64 Mohammadpour and Ziaeenejad 3.3 Simulation Here we have carried out two simulation studies. In the first study, we compare the two method of estimation. To examine the performances of the YW and ML estimators, a Monte Carlo simulation is conducted for different sample sizes n = 100, 300, 500 with m = 200 replications for models Gumbel ARCH(1), Gumbel ARCH(2) and Gumbel GARCH(1,1) for true parameter values: • Gumbel ARCH(1): (a0,a1) = A1(3,0.5); A2(4,0.3); • Gumbel ARCH(2): (a0,a1,a2) = B1(3, 0.4, 0.2); B2(4, 0.3, 0.3); • Gumbel GARCH(1,1): (ao,ai,Pi) = C 1(2, 0.4, 0.3); C2(4, 0.3, 0.3). For the maximization of the log-likelihood function, the YW estimates were used as the initial values. Table 1 provides the mean and mean absolute deviation error (MADE) of the estimators for different values of the parameters and different sample sizes. As for the stationarity discussed, YJi=i «¿(1 + v2) + Pj in the models are all less than one and the stationary condition holds for the fitted models. It can be seen that as the sample size increases, the estimates seem to converge to the true parameter values. Two estimation methods seem to perform reasonably well but the MLEs provide better performance, which was expected. We compare the two models with respect to their general properties and coherent forecasting ability. To achieve this, we have simulated 200 series, each of size 100 from Gumbel GARCH(1,1) process with three sets of parameter values viz.; (a) a0 = 2, ai = 0.2, Pi = 0.3, (b) ao = 3, ai = 0.3, Pi = 0.4 and (c) ao = 3, ai = 0.4, Pi = 0.6. From Table 2 it can be observed that the values of the performance measures mean square prediction error (MSPE) for the Gumbel GARCH(1,1) model are relatively lower, supporting the fact that if the actual process is Gumbel GARCH(1,1), then it gives a better fit (less prediction error) than classic GARCH(1,1). It can be also seen that, as k in k-step ahead forecasts increases, the values of the measures increases, indicating that the error in forecast increases as lag increases. 4 Real Example In this section, we discuss some possible applications of the Gumbel GARCH model for a real time series of weekly data of Tehran Price Index (TEPIX). The data consist of 152 observation from 11 Nov. 2012 to 9 Sep. 2015. The data are obtained from the website of http://www.tse.ir. The sample paths, difference of log data, autocorrelation functions (ACFs) and partial autocorrelation functions (PACFs) of the difference log series are displayed in Figure 1. For selecting the model for the data series, we compare the classical ARCH(1), ARCH(2), ARCH(3), GARCH(1,1), GARCH(1,2), GARCH(1,3), Gumbel ARCH(1), Gumbel ARCH(2), Gumbel ARCH(3), Gumbel GARCH(1,1), Gumbel GARCH(1,2), Gumbel GARCH(1,3) models. For each model, based on MLEs we provide some well-known measures of goodness-of-fit statistics to check the adequacy of a time series model among a finite set of models. These statistics are the Akaike information criterion (AIC) and Bayesian information criterion (BIC). The obtained results, for Gumbel GARCH Model. 65 the data series, are shown in Table 3. As it can be seen from this table, the values of the goodness-of-fit statistics are the smallest for the Gumbel GARCH(1,3) model. The Pearson residual for Gumbel GARCH(1,3) is defined by Xn — E {Xn\Tn—l) where y/Var(Xn |Fn— l) E(Xn|Fn-1) = VTn; Var(X,n\Tn—i) = -n Yn. 6 Figure 1: (a) Sample path of the time series, (b) Difference log time series, (c) ACF, (d) PACF The residual analysis is shown in Figure 2. Figure 2 shows density and normal Q-Q plot for Pearson residuals which appear to be almost normally distributed. The Ljung-Box 66 Mohammadpour and Ziaeenejad statistic is 15.0825 by 15 lags (x0.05(14) = 23.6848). The results show that the residuals are independent. Figure 2: (a) ACF, (b) Q-Q plot, (c) density of Residuals 5 Conclusion This article discusses the financial time series modeling with potential extreme observations. The Gumbel GARCH model, a generalization of the classic GARCH model is proposed to modeling. Stationarity conditions are given as well as the autocorrelation function. For estimation, we present two approaches with the focus on the maximum likelihood method. Simulation results show that two estimation methods are sufficiently accurate but the MLEs provide better performance. Results on the real stock data indicate that the proposed method performs better than the classic GARCH model. Table 1: Mean (± MADE) of the estimators for different values of the parameters Model n Method ao ai a2 or ft Ei=i(1 + V2 )a + ft A1 100 YW 3.2348 ± 0.3991 0.4431 ± 0.1284 ML 3.1525 ± 0.3580 0.4773 ± 0.1209 0.2324 300 YW 3.1866 ± 0.3272 0.4593 ± 0.1075 ML 3.1540 ± 0.2660 0.4834 ± 0.0858 0.2368 500 YW 3.1475 ± 0.2621 0.4741 ± 0.0863 ML 3.0663 ± 0.2171 0.5012 ± 0.0512 0.2455 A2 100 YW 5.1435 ± 1.1690 0.2524 ± 0.0799 ML 4.6023 ± 0.8525 0.2660 ± 0.0587 0.1303 300 YW 4.7894 ± 0.8179 0.2721 ± 0.0612 ML 4.3769 ± 0.6936 0.2811 ± 0.0551 0.1377 500 YW 4.4469 ± 0.4858 0.2920 ± 0.0468 ML 4.1923 ± 0.4759 0.2964 ± 0.0436 0.1452 B1 100 YW 2.5347 ± 0.6879 0.3324 ± 0.1212 0.1438 ± 0.1233 ML 2.8489 ± 0.5592 0.3756 ± 0.1225 0.1772 ± 0.1020 0.3612 300 YW 2.6893 ± 0.5626 0.3466 ± 0.1028 0.1621 ± 0.0961 ML 2.9035 ± 0.4659 0.3723 ± 0.0838 0.1893 ± 0.0768 0.3717 500 YW 2.8767 ± 0.4542 0.3595 ± 0.0865 0.1794 ± 0.0740 ML 3.0288 ± 0.3821 0.3902 ± 0.0753 0.1988 ± 0.0598 0.3899 B2 100 YW 5.6798 ± 2.1657 0.2731 ± 0.1643 0.2432 ± 0.1905 ML 4.8586 ± 1.6984 0.2762 ± 0.1593 0.2463 ± 0.1732 0.3816 300 YW 5.0924 ± 1.6661 0.2794 ± 0.1528 0.2653 ± 0.1716 ML 4.9401 ± 1.5831 0.2804 ± 0.1468 0.2695 ± 0.1690 0.4068 500 YW 4.5712 ± 1.1649 0.2812 ± 0.1417 0.2851 ± 0.1563 ML 4.4499 ± 1.1028 0.2827 ± 0.1377 0.2876 ± 0.1539 0.4261 C1 100 YW 3.5619 ± 1.6097 0.2761 ± 0.0634 0.1579 ± 0.2611 continued ... continued Model n Method ao ai a2 or ft Ei=i(1 + v2 )a + ft ML 2.4378 ± 0.3580 0.3673 ± 0.0437 0.2379 ± 0.1842 0.4178 300 YW 2.8406 ± 0.9385 0.3022 ± 0.0420 0.2029 ± 0.2138 ML 2.2283 ± 0.4981 0.3767 ± 0.0208 0.2652 ± 0.1860 0.4498 500 YW 2.2464 ± 0.3405 0.3243 ± 0.0299 0.2571 ± 0.1746 ML 2.1563 ± 0.1184 0.3854 ± 0.0184 0.2983 ± 0.1545 0.4887 C2 100 YW 5.4571 ± 2.7082 0.2785 ± 0.1701 0.1376 ± 0.2537 ML 4.4554 ± 2.2111 0.2898 ± 0.1711 0.2617 ± 0.1796 0.4037 300 YW 4.9820 ± 2.2210 0.2876 ± 0.1518 0.1907 ± 0.1911 ML 4.4697 ± 2.0628 0.2967 ± 0.1546 0.2778 ± 0.1859 0.4231 500 YW 4.5606 ± 1.7099 0.2898 ± 0.1346 0.2327 ± 0.1339 ML 4.0485 1.5657 0.2904 0.1365 0.2978 0.1298 0.4400 © sa a © K a s a. i a' 01 S j a a. Table 2: MSPE of the data generated from Gumbel GARCH(1,1) model for k-step ahead forecasts k GARCH(1,1) Gumbel GARCH(1,1) GARCH(1,1) Gumbel GARCH( 1,1) GARCH(1,1) Gumbel GARCH( 1,1) (ao,ai,^i) = (2, 0.2, 0.3) (ao,ai,A) = (3,0.3, 0.4) (ao,«i,^i) = (3, 0.4, 0.6) 1 0.84 0.69 0.76 0.67 0.72 0.67 2 0.93 0.82 0.90 0.83 0.81 0.77 3 1.04 0.94 1.01 0.90 0.93 0.86 4 1.10 0.97 1.12 1.03 1.07 1.01 5 1.13 1.02 1.18 1.08 1.20 1.10 Table 3: Estimated parameters, AIC and BIC for the TEPIX time series Model ao ai a2 a3 ^i & AIC BIC ARCH(1) 0.4362 0.3748 350.6 356.6 ARCH(2) 0.4213 0.3106 0.2219 339.4 345.5 ARCH(3) 0.4147 0.3216 0.3105 0.2221 335.6 340.7 GARCH(1,1) 0.3927 0.2520 0.1984 350.9 361.3 GARCH(1,2) 0.3524 0.2212 0.1874 0.1523 351.5 363.9 GARCH(1,3) 0.3021 0.1980 0.1722 0.1281 0.1125 347.7 352.2 Gumbel ARCH(1) 0.4993 0.3216 318.2 322.7 Gumbel ARCH(2) 0.4325 0.2992 0.2579 315.3 319.7 Gumbel ARCH(3) 0.4091 0.2537 0.2280 0.1852 312.4 317.1 Gumbel GARCH(1,1) 0.3423 0.2758 0.1980 306.4 314.9 Gumbel GARCH(1,2) 0.3219 0.2227 0.1825 0.1653 304.7 311.5 Gumbel GARCH(1,3) 0.3051 0.2113 0.1812 0.1521 0.1210 303.2 307.4 K m s R H 0 01 CT\ VO 70 Mohammadpour and Ziaeenejad References [1] Baillie, R. T. and Bollerslev, T. (1989): Common stochastic trends in a system of exchange rates. The Journal of Finance, 44, 167-181. [2] Bollerslev, T. (1986): Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31, 307-327. [3] Calzolari, G., Halbleib, R. and Parrini, A. (2014): Estimating GARCH-type models with symmetric stable innovations: Indirect inference versus maximum likelihood. Computational Statistics and Data Analysis, 76, 158-171. [4] Goldberg, S. (1958): Introduction to Difference Equations, with Illustrative Examples from Economics, Psychology, and Sociology. New York: Wiley. [5] Glosten, L. R., Jagannathan, R. and Runkle, D. E. (1993): On the relation between the expected value and the volatility of the nominal excess return of stocks. Journal of Finance, 48, 1779-1801. [6] Hsieh, D. A. (1988): The statistical properties of daily foreign exchange rates, 19741983. Journal of International Economics, 24, 129-145. [7] Jorion, P. (1988): On jump processes in the foreign exchange and stock markets. Review of Financial Studies, 1(4), 427-445. [8] Liu, S. and Brorsen, B. (1995): Maximum likelihood estimation of a GARCH-stable model. Journal of Applied Econometrics, 10(3), 273-285. [9] Nelson, D. B. (1990a): ARCH models as diffusion approximation, Journal of Econometrics, 45, 7-38. [10] Nelson, D. B. (1990b): Conditional heteroskedasticity in asset returns: A new approach. Econometrica, 59, 347-370. [11] Nelson, D. B. and Cao, C. Q. (1992): Inequality constraints in the univariate GARCH model. Journal of Business and Economic Statistics, 10(2), 229-235. A Appendix Let It be the a-field generated by {at2, a''_1,...}, then we have E(X' |Ft-i, It) = E (Xt2|Ft_i) = + v26/n2), and a2 := E(X2) = E(E(X'|Ft_i)) = E(a'(1 + v26/n2)). Gumbel GARCH Model. 71 For k > 0, since Cov(Xt2 - *2(1 + v26/n2),*2_fc(1 + v26/n2)) = E[(X2 - *2(1 + v26/n2))(*2_fc(1 + v26/n2) - a2)] = E[(at2_fc(1 + v26/n2) - a2)E(X2 - *2(1 + v26/n2)|It)] = E[(a2_fc(1 + v26/n2) - a2)E(E(Xt2|Ft_i,Zt)|I - at2(1 + v26/n2))] = E(at2_fc(1 + v26/n2) - a2)E(at2(1 + v26/n2)|I - *2(1 + v26/n2)) we obtain Cov(Xt2,at2_fc (1 + v26/n2)) = Cov(*2(1 + v 26/n2),at2_fc (1 + v26/n2)). Similarly for k < 0, since Cov(Xt2,Xt2_fc - *2_fc(1 + v26/n2)) = E[(X2 - a2)(Xt2_fc - *2_fc(1 + v26/n2)))] = E[(X2 - a2)E((Xt2_fc - OU(1 + v26/n2))|Ft_fc_i)] = E[(X2 - a2)(*2_fc(1 + v26/n2) - E(*2_fc(1 + v26/n2)|Ft_fc_i)] = 0, we have Cov(Xt2,at2_fc (1 + v 26/n2)) = Cov(Xt2,Xt2_fc). Therefore CM*?,*2 ,(1 + v26/n2)) J*1 + ^^»■ k > 0 Let us now derive y22 (k). For k > 0, since p 722 (k) = Cov(*t2, ) = J] Cov(Xt2_i, *t2_fc) + &Cov(*t2_j, ) i=1 j=1 J(1 + v26/n2)/(1 + v26/n 2)aiC 1, we obtain YX(k) as 7X(k) = Cov(Xt°,X0-fc) = E (Xto,X0-fc) = E (E (Xt2,Xt2_k )|Ft_i) = E (E (XO |Ft-i, Xt2-fc)) = E ((1 + v2 6/n2 K° ,Xt2_k) p q = (1 + v2 6/n2 )E (ao + £ a Xt2_i + £ ft at2_j ,Xt2_fc) i=i j=i = (1 + v2 6/n2 )(£ aiE (Xt2_i,Xt2_k) + £ ft E (at2_j ,Xt2-fc i=i j=i p min(k_i),q = ¿(1 + v 2 6/n2aiCov(Xt2_i, Xt2_j) + £ ft Cov(Xt2_j, Xt2_k) i=i j=i + £(1 + v26/n 2)2ftCov(at2_j, at2_k)). □ j=i