Olivier Darne

Amélie Lemna

Charles

The accuracy of asymmetric GARCH model estimation

Keywords: EGARCH, GJR-GARCH, TARCH, APARCH, accuracy, forecasting, software. JEL Classification: C22

This paper reviews eight software packages when estimating asymmetric GARCH models (from their default option). We consider the numerical con-sistency of GJR-GARCH, TGARCH, EGARCH and APARCH estimations with Normal and Student distributions as well as out-of-sample forecasting accuracy, using the model confidence set procedure. We show that results are clearly software-dependent for both asymmetric volatility models, especially for the t-ratios. The out-of-sample forecast results show that the differences in estimating symmetric and asymmetric GARCH models imply slight differences in terms of forecast accuracy, not statistically significant, except in few cases from the QLIKE loss function. Further, the results indicated that the different specifi-cations of the asymmetric GARCH-type models used by the different packages appear to have no significant effect on their forecast accuracy.

Introduction

Understanding the behavior of volatility in financial assets is important to risk management, derivative pricing and hedging, market making, market timing, portfolio selection, and many other financial activities. Autoregressive conditionally heteroscedastic (ARCH) models introduced by Engle (1982) and extended to GARCH models, independently, by Bollerslev (1986) and Taylor (1986), were developed to capture the two most important stylized facts of returns, which are heavy-tailed distributions and volatility clustering. GARCH estimation is widely available in a number of software packages (e.g., EViews, Gauss, Gretl, Matlab, Ox, R, Rats, Stata,S-Plus, TSP).

A number of reviews have investigated the accuracy, ease to use, availability of documentation and other attributes of the softwares available for the estimation of univariate (symmetric) GARCH models (see, e.g., Brooks, 1997; McCullough and Renfro, 1999; Brooks et al., 2001). McCullough and Renfro (1999) and Brooks et al. (2001, 2003) discussed numerical accuracy issues associated with maximizing the GARCH log-likelihood in several commercial softwares. They found that parameter starting values, initializations for the error and conditional variance-covariance series, optimization algorithm choice, use of analytic or numerical derivatives, and convergence criteria all influence the resulting numerical estimates of the GARCH parameters. As argued by Brooks et al. (2003) "a thorough examination of all of these issues is virtually impossible since the packages on the whole simply do not give sufficient detail on these points". [START_REF] Alexander | Market Risk Analysis[END_REF] There is extensive empirical evidence about the existence of an asymmetric response of volatility to positive and negative past returns. In particular, increases in volatility are larger when previous returns are negative than when they have the same magnitude but are positive. The asymmetric volatility property is explained in [START_REF] Alexander | Market Risk Analysis[END_REF] McCullough and Renfro (1999) compared seven (anonymous) packages by using the FCP GARCH benchmark of Fiorentini et al. (1996) on the daily percentage nominal returns for the Deutschemark/British pound exchange rate of Bollerslev and Ghysels (1996). Brooks et al. (2001) used the same dataset by comparing nine software packages (EViews 3.1, Gauss 3/Fanpac, LIMDEP 8.0, Matlab 11, Microfit 4.0, Rats 4.3, SAS 6.12, Shazam 8.0, and TSP 4.5.). Laurent and Peters (2002) also compared GARCH accuracy between their Ox-G@RCH 2.3 package and four econometric softwares (EViews 4, PcGive 10, S-Plus 6 and TSP 4.5).

2

the literature in terms of the leverage effect (Black, 1976; Christie, 1982; Schwert, 1989) and the volatility feedback effect (French et al., 1987; Campbell and Hentchel, 1992; Bollerslev et al., 1992). [START_REF] Black | Studies of stock market volatility changes[END_REF] There is a large number of non-linear GARCH models, including EGARCH (Nelson, 1991), GJR-GARCH (Glosten et al., 1993), TGARCH (Rabemananjara and Zakoian, 1993) and APARCH (Ding et al., 1993) which allows for asymmetry or leverage effect, so that the conditional variance can be affected differently by positive and negative shocks of the same magnitude. Asymmetric GARCH estimation is also widely available in a number of commercial software packages (e.g., EViews, Gauss, Matlab, Ox, Stata) and there are also a few free open source implementations (e.g., R, Gretl).

However, the asymmetric GARCH models can be expressed in various different forms according to the software packages. This is a severe problem that makes interpretation and diagnostic checking of the model specification extremely difficult. To the best of our knowledge, Brooks et al. (2001) and Alexander (2008) are the only studies which compare the numerical consistency of asymmetric GARCH models, especially EGARCH and/or GJR-GARCH models. Brooks et al. (2001) examine EGARCH estimation and forecasting for four econometric softwares (EViews, Microfit, Rats and SAS), and Alexander (2008) compares the estimation of GJR-GARCH and EGARCH models from EViews and Matlab using daily log returns on the FTSE 100.

The aim of this paper is to review a number of the most widely used software packages when estimating asymmetric GARCH models. In this sense, our contribution extends the work of Brooks et al. (2001) who have shown how different software packages can lead to quite different results in the standard GARCH model using updated versions of software packages and new free packages. We consider the numerical consistency of the four most widely estimated asymmetric univariate models of conditional volatility, namely the GJR-GARCH, TGARCH, EGARCH and APARCH models with Normal and Student distributions [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF] , using a number of [START_REF] Black | Studies of stock market volatility changes[END_REF] Hansen and Lunde (2005) compared 330 GARCH-type models and concluded that in the case of the Deutschmark-Dollar exchange rate, no model beats the GARCH(1,1) model, whereas the GARCH (1,1) model is clearly inferior to asymmetric GARCH models in the case of IBM returns. [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF] See Table 3 for a presentation of the GARCH-type models and distributions available in the different software packages.

different econometric software packages. [START_REF] Bollerslev | A conditionally heteroskedastic time series model for speculative prices and rates of return[END_REF] We thus evaluate eight software packages:

five commercial packages, including Financial Analysis Package 3.0 (FANPAC) for GAUSS 17, Econometrics Toolbox for Matlab 13a, G@RCH 8 package in Ox 8.1, Stata 14 and EViews 9, one semi-commercial package with the free Oxford MFE Toolbox for Matlab, and two free packages with the rugarch package for R 3.4 and Gretl 2018a.

To evaluate the software accuracy in estimating a GARCH model Brooks et al. (2001) used the GARCH model proposed by Fiorentini et al. (1996) as benchmark. Brooks et al. (2001) argue that "there is a need for a benchmark and also some consensus

as to what forms of the model are preferable, and the latter should be estimable by all packages" (p.54). Unfortunately, it is currently not the case as the asymmetric GARCH models can be expressed in various different forms according to the software packages. Brooks et al. (2001) argue that the establishment of benchmarks for other standard non-linear models is long overdue. Therefore, we compare the estimation of asymmetric GARCH models from financial time series by using the default option of the software packages. All the packages analyzed use conditional maximum likelihood to estimate the models. We also address the issue of whether the differences in parameter estimation between packages make a difference from a practical perspective, by evaluating the out-of-sample forecasting accuracy. For that, we extend the work The remainder of the paper is structured as follows. Section 2 introduces the GARCH, discusses the maximum likelihood estimation of these models. The comparison in estimating GARCH-type models from the default option of the packages is presented in Section 4. Section 5 proposes a robustness check with a Student distribution, small and medium sample sizes, and 'harmonized' options for all the packages. The forecasting accuracy is displayed in Section 6. Section 7 proposes a discussion and concludes.

GARCH-type models

The GARCH model was developed independently by Bollerslev (1986) and Taylor (1986a). The GARCH model extends the ARCH model by allowing the conditional variance to be dependent upon previous own lags.

Consider the returns series r t , defined by r t = log P t -log P t-1 , where P t is the observed price at time t, follows a Normal GARCH(1,1) model defined as

r t = c + ε t (1)
ε t = σ t z t , z t ∼ i.i.d.N(0, 1), (2)
σ 2 t = ω + αε 2 t-1 + βσ 2 t-1 (3)
where ε t denotes the 'market shock' or 'unexpected return', and is commonly taken as the mean deviation (r tr), with r is the sample mean. [START_REF] Bollerslev | Quasi-maximum likelihood estimation inference in dynamic models with time-varying covariance[END_REF] The parameters should satisfy ω > 0, α ≥ 0 and β ≥ 0 to guarantee the positivity of the conditional variance.

The stationarity of the process (second-order moment condition) is achieved when the restriction α + β < 1 is satisfied. Ling and McAleer (2002a, 2002b) have derived the regularity conditions of a GARCH(1,1) model, defined as follows:

E[ε 2 t ] = ω 1-α-β < ∞ if α + β < 1, and E[ε 4 t ] < ∞ if kα 2 + 2αβ + β 2 < 1,
where k is the conditional fourth moment of z t . [START_REF] Brooks | GARCH modeling in finance: A review of the software options[END_REF] Ng and McAleer (2004) show the importance to verify these conditions. The sum of α and β quantifies the persistence of shocks to conditional variance, meaning that the effect of a volatility shock vanishes over time at an exponential rate.

The GJR-GARCH model

The GJR-GARCH model developed by Glosten, Jagannnathan et Runkle (1993) is constructed to capture the asymmetric leverage volatility effect. Specification for the conditional variance of GJR-GARCH(1,1) model is

σ 2 t = ω + α ε 2 t-1 + γ I(ε t-1 < 0) ε 2 t-1 + β σ 2 t-1 = ω + (α + γ I(ε t-1 )) ε 2 t-1 + β σ 2 t-1 (4)
where

I(ε t-1 < 0) = 1 if ε t-1 < 0, and 0 otherwise. The volatility is positive if ω > 0, α > 0, γ ≥ 0, α + γ ≥ 0 and β ≥ 0.
The process is defined as stationary if the constraint Ling and McAleer (2002b) have derived the regularity conditions for a GJR-GARCH(1,1), defined as follows:

α + β + (γ/2) < 1 is satisfied.
E[ε 2 t ] < ∞ if α + β + δγ < 1, and E[ε 4 t ] < ∞ if kα 2 + 2αβ + β 2 + βγ + kαγ + kδγ 2 < 1. 8
Asymmetry exists if γ > 0, i.e. positive and negative shocks of equal magnitude have different effects on conditional volatility. The asymmetry is observed as the impulse (α + γ) of negative shocks, which is larger than the impulse (α) of positive shocks. In this model, good news and bad news have different effects on the conditional variance: good news has an impact of α while bad news has an impact of (α+γ).

The GJR-GARCH model nets the GARCH model when γ = 0.

We can find a different parametrization for the GJR-GARCH model in other software packages. In Stata the dummy variable is defined as I(ε t-1 > 0) = 1 if ε t-1 > 0, and 0 otherwise. Therefore, we have for Stata

σ 2 t = ω + (α Stata + γ Stata ) ε 2 t-1 + βσ 2 t-1 if ε t-1 ≥ 0 = ω -γ Stata ε 2 t-1 + βσ 2 t-1 if ε t-1 < 0 (5)
while we have from the equation ( 4) for Matlab-MFE, Matlab-Ek, Gauss-Fanpac, Rrugarch and Ox-G@rch

σ 2 t = ω + γ Matlab ε 2 t-1 + βσ 2 t-1 if ε t-1 ≥ 0 = ω -(α Matlab + γ Matlab ) ε 2 t-1 + βσ 2 t-1 if ε t-1 < 0 (6)
Finally, we obtain the following correspondence

α Matlab = α Stata + γ Stata γ Matlab = -γ Stata
The GJR-GARCH model is specified in Gretl without dummy variable as follows:

σ 2 t = ω + α (|ε t-1 | -γε t-1 )) 2 + β σ 2 t-1
We have the following correspondence9

α Matlab = α Gretl (1 -γ Gretl ) γ Matlab = 4α Gretl γ Gretl

The TGARCH model

A similar way to model asymmetric effects was introduced by Zakoian (1994), and developed further in Rabemananjara and Zakoian (1993), by defining the threshold GARCH (TGARCH) model:

σ t = ω + α + ε + t-1 + γ -ε - t-1 + βσ t-1 = ω + αI(ε t-1 ≥ 0)ε t-1 + γI(ε t-1 < 0)ε t-1 + βσ t-1 (7)
where

ε + t = ε t if ε t > 0, ε + t ≡ 0 otherwise, ε - t ≡ ε t -ε + t , and I(ε t-1 < 0) = 1 if ε t-1 <
0, and 0 otherwise. The TGARCH model is similar to the GJR-GARCH model but the volatility is specified in terms of σ t instead of σ 2 t . Therefore, the TGARCH model does not impose any positivity restrictions on the volatility coefficients. However, in such a case, σ t is no longer assumed to be positive which can be problematic in terms of probabilistic properties on the standardized innovations. The conditions for existence of moments can be found in He and Teräsvirta (1999a).

We can find different parametrizations for the TGARCH model among the software packages. In Gretl, Eviews 10 and R-rugarch the parametrization is given by

σ t = ω + α Gretl (|ε t-1 | -γ Gretl ε t-1 ) + βσ t-1 (8)
which can be re-written as

σ t = ω + α Gretl (1 -γ Gretl )ε t-1 + βσ t-1 if ε t-1 ≥ 0 ω -α Gretl (1 + γ Gretl )ε t-1 + βσ t-1 if ε t-1 < 0
Matlab-MFE use a dummy variable in the equation of the TGARCH(1,1) model, defined as

σ t = ω + α Matlab |ε t-1 | + γ Matlab I(ε t-1 < 0)|ε t-1 | + βσ t-1 (9)
where I(ε t-1 < 0) = 1 if ε t-1 < 0, and 0 otherwise. This equation can be re-written as

σ t = ω + α Matlab ε t-1 + βσ t-1 if ε t-1 ≥ 0 ω -(α Matlab + γ Matlab )ε t-1 + βσ t-1 if ε t-1 < 0
We have the following correspondence between Gretl/Eviews/R-rugarch and Matlab-MFE

α Matlab = α Gretl (1 -γ Gretl ) γ Matlab = 2α Gretl γ Gretl
In Stata the dummy variable is defined as I(ε t-1 > 0) = 1 if ε t-1 > 0, and 0 otherwise. Therefore, the parametrization is given by

σ t = ω + (α Stata + γ Stata ) ε t-1 + βσ t-1 if ε t-1 ≥ 0 = ω -γ Stata ε t-1 + βσ t-1 if ε t-1 < 0 ( 10
)
10 Note that EViews proposes the estimation of a TGARCH model but its TGARCH model is estimated on σ 2 t and not on σ t as defined in Zakoian (1994) and Rabemananjara and Zakoian (1993). Its TGARCH specification is rather similar to the GJR-GARCH model of Glosten et al. (1993). However, it is possible to estimate the TGARCH model in Eviews from the APARCH model by setting δ = 1.

The relationship between the parametrization of the TGARCH defined by Rabemananjara and Zakoian (1993) and that used in Gretl, Eviews and R-rugarch, which is derived from an APARCH model, is given by: α + = α Gretl (1γ Gretl ) and γ -= -α Gretl (1 + γ Gretl ).

We obtain the following correspondence between Matlab-MFE and Stata 11

α Matlab = α Stata + γ Stata γ Matlab = -γ Stata 2.3
The EGARCH model Nelson (1991) proposes the exponential GARCH model (EGARCH), which can also capture volatility leverage. The EGARCH(1,1) model can be expressed in various different forms according to the software packages. This is a severe problem that makes interpretation and diagnostic checking of the model specification extremely difficult.

• Ox-G@RCH:

ln(σ 2 t ) = ω + (1 -βL) -1 (1 + αL)g(z t-1 ) • Gauss-Fanpac / Matlab-MFE / Matlab-Econometrics / R-rugarch / Stata: ln(σ 2 t ) = ω + θ 1 (|z t-1 | -E(|z t-1 |)) + θ 2 z t-1 + β ln(σ 2 t-1 )
• SAS: ln(σ 2 t ) = ω + αg(z t-1 ) + β ln(σ 2 t-1 )

• EViews / Gretl: [START_REF] Carnero | Effects of outliers on the identification and estimation of the GARCH models[END_REF] The EGARCH(1,1) specification expressed in Ox-G@RCH is given by:

ln(σ 2 t ) = ω + θ 1 (|z t-1 |) + θ 2 z t-1 + β ln(σ 2 t-1 ) with g(z t ) = θ 1 (|z t | -E(|z t |)) + θ 2 z t ,
ln(σ 2 t ) = ω(1 -β) + (1 + αL)g(z t-1 ) + β ln(σ 2 t-1 )
implying that ln(σ 2 t ) is a function of two lags of standardized residuals:

ln(σ 2 t ) = ω(1 -β) + θ 1 (|z t-1 | -E(|z t-1 |)) + θ 2 z t-1 +αθ 1 (|z t-2 | -E(|z t-2 |)) + αθ 2 z t-2 + β ln(σ 2 t-1 )
To estimate the conditional variance as a function of one lag of standardized residuals the EGARCH(0,1) specification has to be chosen in Ox-G@RCH. Note there is a difference in the intercept term specification between Ox-G@RCH and Matlab-MFE, Matlab-Ek, Gauss-Fanpac, R-rugarch and Stata, with the following correspondence:

ω Matlab = ω G@RCH (1 -β G@RCH ).
The EGARCH model with the Normal distribution is slightly differently estimated in EViews and Gretl, with the term E(|z t |) = 2/π removed from the following specification: [START_REF] Carnero | Estimating GARCH volatility in the presence of outliers[END_REF] ln(

σ 2 t ) = ω + θ 1 |z t-1 |) + θ 2 z t-1 + β ln(σ 2 t-1 )
Therefore, the intercept term in the conditional variance equation will differ by a factor of θ 1 2/π compared with that of Matlab/Gauss/R/Stata, with the following correspondence: ω Matlab = ω EViews/Stata + θ 1 2/π.

12 Some authors make a distinction between asymmetry, referred to as the different impacts on conditional volatility of positive and negative shocks of equal magnitude, and leverage effect, regarded as the negative correlation between returns shocks and subsequent shocks to volatility (McAleer, 2014;

Martinet and McAleer, 2018); Chang and McAleer (2017) propose the regularity conditions that an EGARCH(1,1) model obtained from a random coefficient complex nonlinear moving average process, should fulfill to capture asymmetry and/or leverage effects. [START_REF] Carnero | Estimating GARCH volatility in the presence of outliers[END_REF] The software packages advance some computational advantages to remove the term 2/π out of the summation operator.

The APARCH model

Ding, Granger et Engle (1993) propose the asymmetric power ARCH (APARCH) model which nests several ARCH-type models. The APARCH(1,1) model is defined as follows

σ δ t = ω + α (|ε t-1 | -γε t-1 ) δ + βσ δ t-1
where ω > 0, α ≥ 0, and β ≥ 0. Parameter δ (δ > 0) plays the role of a Box-Cox transformation of the conditional standard deviation σ t . γ, with -1 < γ < 1, reflects the asymmetric effect. A positive (resp. negative) value of the γ means that past negative (resp. positive) shocks have a deeper impact on current conditional volatility than past positive shocks. The properties of the APARCH model have been studied by He and Teräsvirta (1999b) and He et al. (2008).

The APARCH model includes several ARCH extensions as special cases, including the GARCH(1,1) model when δ = 2 and γ = 0, the TGARCH(1,1) model when δ = 1, and the GJR-GARCH(1,1) model when δ = 2. [START_REF] Carnero | Leverage effect in energy futures revisited[END_REF] This equation is the parametrization used in Eviews, Ox-G@arch, Gretl and Rrugarch. In MFE-Matlab and Stata the parametrization is slightly different with a positive sign for the asymmetric parameter γ, given by

σ δ t = ω + α (|ε t-1 | + γε t-1 ) δ + βσ δ t-1
3 Estimating GARCH models

Maximum likelihood estimation

The most commonly used method in estimating the vector of unknown parameters θ of GARCH-type models is the method of maximum likelihood (ML) estimation. The function of the parameters set is called the log-likelihood function. [START_REF] Carnero | Identification of asymmetric conditional heteroscedasticity in the presence of outliers[END_REF] Given that ε t = σ t z t , the log-likelihood function L of the standard normal distribution is

L = T ∑ t=1 l t = - 1 2 T ∑ t=1 ln(2π) + ln(σ 2 t ) + z 2 t
where l t = -1 2 ln(2π) + ln(σ 2 t ) + z 2 t , and T is the number of observations. As discussed in McCullough and Renfro (1999) and Brooks et al. (2001), there are several practical issues to consider in the maximization of likelihood function. Starting values for the model parameters θ = (ω, α, β) need to be chosen and an initialization of ε 2 t and σ 2 t must be supplied. [START_REF] Chang | The correct regularity condition and interpretation of asymmetry in EGARCH[END_REF] Once the log-likelihood is initialized, it can be maximized using numerical optimization techniques. The most popular methods of optimization make use of exact or approximated Hessian matrices, and iterative procedure that updates the current values of the estimates θ i+1 at iteration i + 1. A standard approach based on the exact Hessian matrix is the Newton-Raphson algorithm. [START_REF] Clemens | The meaning of failed replications: A review and proposal[END_REF] The iterative procedure is given by

θ i+1 = θ i -λ i H( θ i ) -1 G( θ i ) (11)
where θ i denotes the vector of estimated model parameters, λ i is a scalar (called step) which may be fixed or variable, G(

θ i ) = ∂L( θ i ) ∂ θ i
is the gradient (or score) vector of the log-likelihood function (first derivative of the log-likelihood function with respect to the vector of unknown parameters), and

H( θ i ) = ∂ 2 L( θ i ) ∂ θ i ∂ θ i
is the Hessian matrix (or a suitable approximation of it) of the log-likelihood function (second derivative of the log-likelihood function with respect to the vector of unknown parameters) all evaluated at iteration i.

An alternative procedure based on approximated Hessian matrix is the Berndt-Hall-Hall-Hausman (BHHH) algorithm (Berndt et al., 1974) which only uses first [START_REF] Chang | The correct regularity condition and interpretation of asymmetry in EGARCH[END_REF] Almost all packages that estimate GARCH-type models use the sample variance as default for initializing the variance process (σ 2 t = T -1 ∑ T t=1 ε 2 t ), except EViews which uses backcast exponential smoothing, and Matlab-MFE which uses a local average. Pelagatti and Lisi (2009) analyzed the precision of QML estimates under different choices of initialization and sample dimension. [START_REF] Clemens | The meaning of failed replications: A review and proposal[END_REF] Analytic or numerical derivatives may be used for the Newton-Raphson algorithm. Fiorentini et al.

(1996) provided algorithms for computing analytic derivatives for GARCH models.

derivatives. The BHHH algorithm is similar to the Newton-Raphson algorithm but instead of using the Hessian matrix of second derivatives H( θ i ) it is based on an approximation formed by the sum of outer product of the gradient (OPG) vectors for the contribution of each observation to the objective function. [START_REF] Davidson | Moment and memory properties of linear conditional heteroscedasticity models, and a new model[END_REF] The Hessian matrix may be approximated by

-H( θ i ) ≈ B( θ i ) = T ∑ t=1 G t ( θ i )G t ( θ i ) = T ∑ t=1 ∂l t ∂θ i ∂l t ∂θ i
where G t ( θ) is the gradient of the log-likelihood function for each observation. The BHHH algorithm can be modified by employing the Marquardt correction which adds a correction matrix to the sum of the OPG vectors (Hessian approximation) as follows

B( θ i ) = T ∑ t=1 G t ( θ i )G t ( θ i ) -aI
where I is the identity matrix, and a a positive number chosen by the algorithm.

Another common approach based on approximated Hessian matrix is the Quasi-Newton methods which are to build up an estimate of the Hessian at each iteration, starting from an initial estimate. The inverse of the Hessian matrix is given by

H( θ i ) = H( θ i-1 ) + C i-1
where C is a correction matrix. The Quasi-Newton algorithms differ only in their choice of this matrix, and the most widely used are the DFP algorithm proposed by Davidson (1959) and Fletcher and Powell (1963) or the BFGS algorithm developed by Broyden (1970), Fletcher (1970), Goldfarb (19790) and Shanno (1970). [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF] 18 Computational speed is increased by not calculating the actual Hessian matrix at each iteration for each time step, but the approximation can be weak when the log-likelihood function is far from its maximum, thus requiring more iterations to reach the optimum. [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF] The Quasi-Newton BFGS and DFP methods are like Newton-Raphson method in that they use both first and second derivative information, however, the Hessian matrix is approximated, reducing the computational requirements. These methods take more iterations than the Newton-Raphson algorithm but the use of an approximation produces a gain because it can be expected to converge in less overall time.

The estimates that maximize the conditional log-likelihood L are called the maximum likelihood (ML) estimates. Under suitable regularity conditions, the ML estimates are consistent and asymptotically Normally distributed and an estimate of the asymptotic covariance matrix of the ML estimates is constructed from an estimate of the final Hessian matrix from the optimization algorithm used. However, the appropriate regularity conditions have only been done for some GARCH models. [START_REF] Efron | Assessing the accuracy of the maximum likelihood estimator: Observed versus expected Fisher information[END_REF] Typically, the GARCH models are represented as unconstrained optimization problems. However, if the likelihood maximization must be achieved subject to nonlinear constraints in the vector of unknown parameters, i.e. lower and upper bounds on the parameters, like the stationary constraint 0

< α + β < 1 in the GARCH(1,1) model or the constraint -1 < γ < 1 in the GJR-GARCH(1,1) model, constrained
nonlinear programming (NLP) methods can be used, such as the sequential quadratic programming (SQP) algorithms. [START_REF] Engle | Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation[END_REF]

Estimators of variance-covariance matrix

It is important to know how are computed the standard errors of the parameter estimators. The different methods for computing standard errors for GARCH coefficients can lead to different standard error estimates, and thus to important implications for test significance on the parameters.

The variance-covariance matrix of the parameter estimators when the variables are assumed to be Normally distributed is calculated as the inverse of the Fisher information matrix I( θ), defined as

Q( θ) = I( θ) -1
It is diagonal in the case of a Normal density but otherwise it can be a complicated non-diagonal matrix, and thus when MLE is applied to non-Normal distributions the computations of the covariance matrix of estimators can be difficult. In general two estimators of the Fisher information matrix can be calculated. The first estimator is the expected values of the first derivatives of the log-likelihood function. It is known as the BHHH estimator or OPG estimator:

I( θ) ≈ E ∂lnL ∂θ ∂lnL ∂θ
The second estimator, called Hessian-based estimator, is the expected values of the negative of the second derivatives of the log-likelihood function: 22

I( θ) = -E ∂ 2 lnL ∂θ∂θ = -E H( θ)
These two estimators are asymptotically equivalent, but they could give different results in finite samples. Available evidence suggests that in small or moderate sized samples, the Hessian is preferable (Greene, 2007). However, in most cases, the BHHH estimator will be the easiest to compute.

If the log-likelihood function is maximized using Hessian-based algorithms, the Fisher information matrix can be obtained as a by-product of the optimization procedure, such as Newton-Raphson and BHHH methods. On the other hand, if Quasi-Newton methods, such as BFGS method, are employed, the estimated Hessian matrix is not guaranteed to converge to its exact value when approaching the optimum. Using those methods to obtain the estimated Fisher information matrix may thus be misleading. [START_REF] Engle | Semiparametric ARCH models[END_REF] The assumption of Normally distributed standardized innovations is often violated by the data, especially in finance. This has motivated the use of alternative distributional assumptions (e.g., Student-t or Skewed-Student-t distributions). Alternatively, [START_REF] Engle | Modelling the persistence of conditional variances[END_REF] Often the expected value of the second derivatives of the log likelihood have to be approximated by taking their actual values, rather than the expected values at the maximum likelihood estimates (Efron and Hinkley, 1978).

the MLE based on the Normal density may be given by a quasi-maximum likelihood (QML) interpretation. Bollerslev and Wooldridge (1992) showed that the maximization of the Normal log-likelihood function results in QML estimates of the parameter vector θ that are consistent and asymptotically Normally distributed even when the distribution of z t is non-Normal. In addition, Bollerslev and Wooldridge (1992) derived an asymptotic covariance matrix for the QMLEs that is robust to conditional non-Normality. [START_REF] Engle | A permanent and transitory component model of stock return volatility[END_REF] This matrix is estimated using both matrices of second order derivatives and outer products of first derivatives:

Q( θ) = H( θ QML ) -1 B( θ QML )H( θ QML ) -1
where θ QML is the QML estimator of θ, which is also called the "sandwich" estimator, and the coefficient standard errors as the (robust) Bollerslev-Wooldridge standard errors.

Under correct specification the QML estimator is equivalent to the other covariance estimators, as it gives asymptotically the inverse of the Fisher information matrix (thus, block-diagonal).

We evaluate eight of the most widely used software packages when estimating asymmetric GARCH models: Financial Analysis Package 3.0 (Fanpac) for Gauss 17

(Gauss-Fanpac), Econometrics Toolbox for Matlab 13a (Matlab-EK), the free Oxford MFE Toolbox for Matlab (Matlab-MFE) [START_REF] Engle | The Spline-GARCH model for low-frequency volatility and its global macroeconomic causes[END_REF] (Sheppard, 2009), G@RCH 8 package in (Galanos, 2018). [START_REF] Fiorentini | Analytic derivatives and the computation of GARCH estimates[END_REF] Table 1 summarizes the MLE in the default option of the econometric packages in terms of numerical optimization techniques and standard errors of the estimated parameters. In terms of optimization routine, Gauss-Fanpac, Matlab-MFE, Gretl and Ox-G@RCH use a version of the BFGS algorithm whereas EViews uses the BHHH algorithm with the Marquart correction. Stata switches between the BHHH and BFGS algorithms. [START_REF] Fletcher | Practical Methods of Optimization[END_REF] Finally, Matlab-Econometrics and R-rugarch apply a SQP algorithm. [START_REF] Francq | Maximum likelihood estimation of pure GARCH and ARMA-GARCH processes[END_REF] In terms of standard errors for GARCH coefficients, Gauss-Fanpac, EViews, Stata It is important when comparing econometric softwares to have an appropriate dataset, especially to avoid the misleading effects of outliers on the estimation of the volatility. [START_REF] Ghalanos | Rsolnp: General non-linear optimization using augmented Lagrange multiplier method[END_REF] For that, we have checked that all observations are lesser than 5 σ t , where σ t is an estimate of the conditional standard deviation. [START_REF] Glosten | On the relation between the expected value and the volatility of the nominal excess return on stocks[END_REF] All packages should be able to es-Table 1: Maximum likelihood estimation in the default option of econometric packages. The GARCH models. The intercepts in the conditional variance equation (ω) are similar for all the packages for both returns. However, it is the GARCH parameters where the differences across packages appear. The parameters on the lagged squared errors (α) are close for Gauss-Fanpac, Matlab-Econometrics, Ox-G@RCH, Stata, Gretl and R-rugarch at around 0.71, and higher than for EViews and Matlab-MFE at around 0.61. We find similar results for the parameters on the lagged conditional variance, but the values are higher for EViews and Matlab-MFE than those for the other packages. Note that the levels of persistence (α + β) are broadly in agreement for all the packages, with a persistence around 0.988.

For the t-ratios associated with the coefficients (Tables 4 and5) it is evident that the differences across packages are more marked than they were for the parameter estimates. For example, the t-ratios for the GARCH parameters given by EViews, Stata and Matlab-Econometrics (based on OPG estimator) are around two times higher than those by Ox-G@rch, Gretl and R-rugarch (based on robust QMLE estimator). However, none of these differences are important for tests of significance since all of the parameters are statistically significant at the 1% level under all packages, except the intercept for Matlab-MFE and Ox-G@rch. The differences in t-ratios are arguably unsurprising since a similar result was found by Brooks et al. (2001) for GARCH(1,1) models, who argued that differences in default calculation methods for the standard errors can explain these differences. However, these results should be take with caution because the robust QML estimator is prone to give very different results in the case where the standardized innovations are not Normal. As a consequence, it is likely that GARCH models.

the standardized innovations issued from the series are not Normal because softwares that implement robust standard errors give very far away results. Inference, and thus t-ratio must be biased for the other software packages.

The GJR-GARCH models. As for the GARCH model the intercepts are broadly in agreement for all the packages. The parameters on the lagged conditional variance are close for Gauss-Fanpac, Matlab-Econometrics, Ox-G@RCH, Stata, Gretl and R-rugarch, and lower than for EViews and Matlab-MFE. The parameters on the lagged squared errors are also the same for Gauss-Fanpac, Matlab-Econometrics, Ox-G@RCH, Stata (with α * ), Gretl and R-rugarch, but its values are negative or null for EViews and Matlab-MFE, which do not satisfy the condition α > 0 to have the volatility positive. [START_REF] Hamadeh | Asymptotic properties of LS and QML estimators for a class of nonlinear GARCH processes[END_REF] The asymmetric parameters (γ) are close for Gauss-Fanpac, Matlab-Econometrics, Ox-G@RCH, Stata (with a negative sign), Gretl and R-rugarch, and for EViews and Matlab-MFE.

When one considers the t-ratios, the discrepancies between methods used by each of the packages become far more apparent. For example, the t-ratio on the asymmetric term in the conditional variance equation (γ) varies from 1.44 to 5.79. These strong differences are important for test significance since the asymmetric parameters are not significant at the 10% level under Ox-G@RCH and R-rugarch, and significant for the other packages. This is strong difference can be explained by the (possible) non Normality of standardized innovations.

The EGARCH models. Even if the constant terms in the conditional variance equation are different across the packages due to the different forms to express the EGARCH models, once we modify the values of the coefficients (ω * ) in a consistent manner across the packages as defined in Section 2.2 we find that the intercepts are broadly in agreement for all the packages and returns. [START_REF] Hansen | A test for superior predictive ability[END_REF] The parameters on the lagged [START_REF] Hamadeh | Asymptotic properties of LS and QML estimators for a class of nonlinear GARCH processes[END_REF] For Matlab-MFE the zero value for the lagged squared errors can be explained by the fact that its maximization package fminunc stopped the optimization with the following message "because the size of the current step is less than the selected value of the step size tolerance". For Eviews by changing its default option in initializing of the conditional variance by the sample variance we found the same values of the parameters that those obtained from the other packages. [START_REF] Hansen | A test for superior predictive ability[END_REF] Alexander (2008) compared the estimation of GARCH(1,1), GJR-GARCH(1,1) and EGARCH(1,1) models from EViews and Matlab using daily log returns on the FTSE 100 over the period 2003-2007.

20

squared errors are close for all the packages. However, the discrepancies become apparent for the two asymmetric parameters in the conditional variance equation. For

Gauss-Fanpac, Matlab-Econometrics, Gretl, Stata and R-rugarch the coefficient values of the magnitude effect (θ 1 ) are close to one another, while those of EViews, Matlab-MFE and Ox-G@RCH are typically up to 30% away from the three others. We also find differences for the parameters of the sign effect (θ 2 ) for all the packages, which vary between -0.070 and -0.037, except for Gauss-Fanpac, Gretl, Stata and R-rugarch where the coefficients are very close. The differences across the packages for the EGARCH models can be explained by both the various specifications of the model and the estimation methods used by each of the packages.

As found previously the results also show differences across packages for the tratios. For example, the t-ratios for the lagged conditional variance displayed by EViews are around 3 times higher than those by Ox-G@rch and R-rugarch. However, these differences are not important for tests of significance since all of the parameters are statistically significant at the 1% level under all packages, except Stata with the parameter θ 2 .

The TGARCH models. All the coefficients are different between all the packages, except for the parameter ω. When taking the corresponding coefficients for α and γ we find similar values for Gretl and Stata. The differences are more important in terms of t-ratios, and the TGARCH model is rejected for R-rugarch since the asymmetric parameter is not significant at the 10% level.

The APARCH models. As for the TGARCH model the coefficients are different across the packages. Nevertheless, we find similarities for the values of parameters for some packages, such as α, β and γ are close for Stata and Gretl, α for Matlab-MFE and R-rugarch, β and γ for Eviews and Matlab-Econometrics. Note that the power parameter δ is close between the packages, ranging between 1.4 and 1.8, except for R-rugarch which is almost two times higher. The results also show differences for the t-ratios,

Her results showed that the parameter estimates are similar, except that the EViews optimizer for EGARCH model converges to a solution that is not sensible because the estimated long term volatility,

given by ln(σ 2 ) = ω 1-β , is too low. However, if we use the intercept term modified ω * we find that the long term volatility estimated by both packages is close.

and Gretl (based on robust QMLE estimator) which reject the APARCH model. [START_REF] Hansen | A forecast comparison of volatility models: Does anything beat a GARCH(1,1)[END_REF] Overall, for the asymmetric GARCH-type models the GJR-GARCH and APARCH models are not rejected only for Gauss-Fanpac/Matlab-Econometrics and Eviews/Stata, respectively, whereas all EGARCH and TGARCH models are not rejected, whatever the packages, except Stata and R-rugarch, respectively. [START_REF] Hansen | Model confidence sets for forecasting models[END_REF] Further, the strong differences in terms of t-ratios between the packages using the OPG estimator and those using the robust QML estimator is likely due to the fact taht the standardized innovations are not Normal.

Robustness

Student distribution

As shown in Table 2 the FTSE returns are characterized by non-Normality, with evidence of significant negative excess skewness and kurtosis. We re-estimate the five GARCH-type models with a Student distribution, which is the non-Normal distribution available for all the software packages (see Table 3). The results are given in Tables 10 and11.

The GARCH models. We find similar results than those obtained with a Normal distribution, namely the parameter ω is again similar for all the packages and the parameter α and β are similar for Ox-G@rch, Matlab-Ek, R-rugarch, Stata and Gretl whereas they are lower and higher, respectively, for Eviews and Matlab-MFE. Opposed to the results obtained when the standardized innovations are assumed Normal, the results for the t-ratios are now slightly different when the standardized innovations follow a Student distribution since the t-ratios based on QML estimator are slightly lower than those based on OPG estimator. The parameter η of the Student distribution is close for all the packages.

The asymmetric models. The results are also identical to those with a Normal distribution in terms of comparison of parameters for the GJR-GARCH, EGARCH and TGARCH models. The differences in terms of t-ratios are lesser across the packages. However, some t-ratios based on the QML estimator are strongly affected.

For example, Ox-G@rch and Gretl do not reject the GJR-GARCH model as the asymmetric parameter γ is become significant at the 10% level or Matlab-MFE rejects the TGARCH model as the parameter α is become not significant at the 10% level.

Finally, all the packages reject the APARCH model with a Student distribution.

Sample size

We now evaluate the effect of the sample size on the estimations of the five GARCHtype models by using medium and small samples, namely 500 and 100 observations, respectively.

For the medium sample size the GARCH model is well estimated by all the packages with similar values for all the parameters, except for Eviews and Matlab-MFE (Tables 8 and9). For these both packages their values of the parameter on the lagged squared errors (α) are very different (0.013 and 0.004 for Eviews and Matlab-MFE, respectively) compared to those of the other packages (0.025) and their t-ratios are not significant. However, only few asymmetric models are not rejected with a Normal distribution (the GJR-GARCH model for Matlab-Econometrics and the EGARCH models for Eviews and Ox-G@rch) because some parameter estimations are not significant and positivity or stationarity constraints are not satisfied.

When analyzing Tables 6 and7 for the small sample size we observe that the results are worst since almost all the GARCH-type models are rejected, except the EGARCH (Matlab-MFE and R-rugarch) and TGARCH (R-rugarch) models, whatever the distribution. [START_REF] He | Properties of moments of a family of GARCH processes[END_REF] Note that in most models Gauss-Fanpac, Ox-G@rch and Stata do not succeed to estimate the parameters due to problem of convergence. This result confirms those obtained by Ng and Lam (2006) and Hwang and Pereira (2006) showing that the estimation of the GARCH-type models is dramatically affected by the sample size.

Identical algorithm and variance-covariance estimator

We now compare the accuracy of the packages by using the same optimization algorithm and variance-covariance estimator. This comparison can be interesting when a study estimate the volatility of asset returns using the GARCH-type models available in an econometric package from the previous options and another study compare these estimations with other GARCH-type models not available in the previous package (see Table 3), as a robustness study by reanalysis (Clements, 2017). For that we use the BFGS algorithm and the Hessian-based estimator. The choice of these options is guided by their availability in all the packages, except for Matlab-Econometrics.

We also use the sample variance as option for initializing the variance process. The other options, such as parameter starting values, are not harmonized because they are impossible to modify in some packages. The results are given in Tables 12 and13.

For the GARCH model all the packages estimate the same values for the parameters and similar values for the t-ratios. We obtain similar results for the GJR-GARCH model, with one difference in terms of t-ratios for the parameter α estimated by Stata likely due to the different parametrization used by this package. The differences for the t-ratios are more important for the EGARCH models, even if the values of the parameters are similar. For example, the parameter θ 2 is not significant for R-rugarch whereas it is significant for the other packages. For the TGARCH and APARCH models the estimations are close for Eviews, Stata and Gretl whereas they are slightly different for Matlab-MFE and R-rugarch, respectively. These results

show that the parametrization of the asymmetric models can be important, sometimes dramatically with the rejection of the models. [START_REF] He | Statistical properties of the asymmetric power ARCH process[END_REF] One way to produce dynamic out-of-sample volatility forecasts, is to output an outof-sample forecast for a one-day horizon, and then assess model forecast performance ).

The one-step ahead forecast for T + 1 is given for each model by but this issue is beyond the scope of this paper. Note that we have also estimated the GARCH-type models for the medium and small samples and obtained that the asymmetric models are rejected for both samples. The results are available in the online Appendix.

• GARCH model σ2 T +1 = ω + αε 2 T + β σ2 T • GJR-GARCH model σ2 T +1 = ω + α + γ Î(ε T ) ε2 T + β σ2 T • EGARCH model σ2 T +1 = exp( ω) exp ( ĝ(z T )) σ2 β T • TGARCH model σT+1 = ω + α + γ Î(ε T ) εT + β σT
Note that the one-step ahead forecast for the asymmetric GARCH-type models can be slightly different according to the parametrization of these models (see Section 2).

Most of packages produce out-of-sample conditional variance forecasts automatically, namely EViews (Static forecasting) [START_REF] He | Higher-order dependence in the general power ARCH process and a special case[END_REF] , Gauss-Fanpac (uforecast), Matlab-Ek(forecast), Ox-G@RCH (FORECAST), Stata (predict) and R-rugarch (ugarchforecast). Nevertheless, we also produce conditional variance forecast from Matlab-MFE and Gretl. [START_REF] He | Fourth moment structure of a family of first-order exponential GARCH models[END_REF] Patton (2011) shows that only two loss functions are robust to noise in the volatility proxy: the mean squared error (MSE) and the loss implied by a Gaussian likelihood (QLIKE). We thus assess volatility forecast performance by the MSE and QLIKE loss functions, defined by

MSE = 1 H H ∑ t=1 σ 2 t -σ 2 t 2 QLIKE = 1 H H ∑ t=1 σ 2 t σ 2 t -ln σ 2 t σ 2 t -1
where σ 2 t is 5-min calendar-time RV, and σ 2 t is the forecast volatility. to contain the best models with a given level of confidence. This MCS produces a number of models with the same forecasting performance, and therefore more robust [START_REF] He | Higher-order dependence in the general power ARCH process and a special case[END_REF] In EViews we have the possibility of Dynamic or Static forecasting. Dynamic forecasting performs a multi-step forecast whereas Static forecasting performs a series of one-step ahead forecasts. [START_REF] He | Fourth moment structure of a family of first-order exponential GARCH models[END_REF] The package GIG 2.21 for Gretl does not provide forecasts for the conditional variance for the EGARCH model.

forecasting than that using only one model. The t-statistic is defined as

T max M = max i∈M t i with t i = di var( di ) (12)
where var( di ) denotes the estimate of var( di ), di = m -1 ∑ j∈M di j , and di j = n -1 ∑ n t=1 d i j,t , with d i j,t = L i,t -L j,t for all i, j ∈ M 0 , and L i,t is a loss function (here MSE). [START_REF] Heber | Oxford-Mann Institute's realized library version 0.2[END_REF] The tstatistic is associated with the null hypothesis of equal predictive ability (EPA, Hansen, 2005) H 0,M : E( di ) = 0 for all i ∈ M , where M ⊂ M 0 . The MCS is a sequential testing procedure, eliminating at each step the worst model from M , until the null hypothesis of EPA is accepted for all the models. If the null of EPA is rejected for M = M 0 , the worst-performing model is excluded from the set M . [START_REF] Hentschel | All in the family nesting symmetric and asymmetric GARCH models[END_REF] The iterative procedure stops when the null hypothesis of EPA of the models still included in the set cannot be rejected. If H 0,M is accepted at level α then the MCS is the set M * 1-α . [START_REF] Higgins | A class of nonlinear ARCH models[END_REF] The MSE and QLIKE values as well the MCS associated with each loss function from the packages are presented in Tables 141516, with a Normal distribution, a Student distribution and using the same optimization algorithm and variance-covariance estimator for each package, respectively. The initial model space M 0 consists of 8 (GARCH) to 3 (GJR-GARCH) models estimated from the different packages. We have excluded the APARCH models because the conditions are not satisfied for most of the packages. [START_REF] Hwang | Small sample properties of GARCH estimates and persistence[END_REF] We set the confidence level for the MCS to α = 0.20 and 0.50.

Table 14 gives the results for the Normal distribution. For the GARCH model, the MSE values are very close for each package, which is confirmed by the MCS as all the packages are in M * 50% , suggesting that all the packages forecast similarly in this 40 di j measures the relative sample loss between the i-th and j-th models, while di is the sample loss of the i-th model relative to the average across the models in M . [START_REF] Hentschel | All in the family nesting symmetric and asymmetric GARCH models[END_REF] The choice of the worst model to be eliminated uses the following elimination rule:

e max M = arg max i∈M t i .
case. For the QLIKE loss function the values are slightly different, where Matlab-Econometrics displays the lowest value and is also the only package to be in M * 20% . For the asymmetric GARCH-type models the MSE and QLIKE values are close for each package. Further, most of the packages are in M * 20% for the MSE loss function. However, it is not the case for the QLIKE loss function where only one package is in M * 20% , such as Matlab-Econometrics for the GJR-GARCH model, Gauss-Fanpac for the EGARCH model, and Gretl for the TGARCH model.

The results obtained with a Student distribution are similar than those obtained with the Normal distribution, whatever the GARCH-type model (Table 15). The findings show that we do not obtain strong differences in terms of forecast accuracy between the packages, but only slight differences, not statistically different from the MSE loss function, and but sometimes statistically different from the QLIKE loss function.

These differences can be due to the various specifications of the GARCH-type models and/or the estimation methods used by each of the packages.

Finally, Table 16 displays the results by using the same optimization algorithm (BFGS), variance-covariance estimator (Hessian-based estimator) and initialization of the variance process to examine whether the different parametrizations of the asymmetric GARCH models can have an effect on the forecast performance. The results show that the MSE values are very close for all the packages and GARCH-type models, with similar forecast accuracy as they are in M * 50% . We obtain similar results from the QLIKE loss function, with most of the package in M * 20% , except for the GJR-GARCH model. Overall, the different parametrizations used in the packages seem to have no impact on the forecast performance.

Discussion and conclusion

The aim of this paper was to review a number of the most widely used software packages when estimating symmetric and asymmetric GARCH models (using their default option). We considered the numerical consistency of GARCH, GJR-GARCH, TGARCH, EGARCH and APARCH estimations, with Normal and Student distributions, as well as out-of-sample forecasting accuracy, using the model confidence set procedure.

As found by Brooks et al. (2001) results are clearly software-dependent for both asymmetric volatility models, especially for the t-ratios, which makes the results from two different papers, which probably used different packages, not comparable. Further, the different specifications of the asymmetric GARCH models can have an impact on the significance of the parameters, especially for the EGARCH model.

The out-of-sample forecast results show that the differences in estimating symmetric and asymmetric GARCH models imply slight differences in terms of forecast accuracy, not statistically significant, except in few cases with the QLIKE loss function.

The results indicated that the different specifications of the asymmetric GARCH-type models used by the different packages appear to have no significant effect on their forecast accuracy.

The aim of this study was not 'pass judgement' on the packages that are good or not, since defaults may differ between packages in terms of initializations, starting values for parameters, convergence tolerances, optimization routines . . . but we can do some recommendations for the users. Our results showed that these choices can be important in terms of parameter values, t-ratios and, sometimes, forecasting, and thus the results can be not comparable between the packages. However, in most of times, the users do not give information on these options. In few cases the optimization algorithm and the method for computing the standard errors are given, and sometimes the package employed but not its version. The package used for the estimation is also important to know what parametrization is used for the asymmetric GARCH-type models. In this study we showed that the various specifications can lead to different values of parameters for the same model, but we provided the correspondences between the packages, making these estimations comparable. The version of the package is also important because some options can be improved or modified. For example, before the version 9 of Eviews the optimization algorithm in its default option was the BHHH algorithm, and now its is the BFGS algorithm.

Overall, we think that, based on the results obtained in our study, software users should mention the name (and version) of the package used in estimation, what the default estimation methods are, such as method to optimize the likelihood function, how the standard errors are computed, and so on. This information will be also important for reproducibility, replicability and robustness which are fundamental characteristics of scientific studies (Burman et al., 2010; Clemens, 2017).

Our estimation results on the small, medium and large sample sizes showed that the estimation of GARCH-type models is sensitive to the sample size. As argued by Ng and Lam (2006) and Hwang and Pereira (2006) we recommend using at least 700 observations for estimating the (symmetric and asymmetric) GARCH-type models.

Larger samples would be preferable but increase the presence of outliers, jumps or breaks which can lead to misleading effects on the estimation of the volatility (see, e.g., Franses and Ghijsels, 1999; Carnero et al., 2007, 2012, 2016).

There are no consensus on the choices of initial values, optimization routines and standard errors computation methods. However, it would be possible to propose a 'benchmark' for these options in order to make the results obtained from different packages comparable. These choices depend on the options proposed by the packages as well as their 'flexibility' to change the options. As we do in Section 5.3, a standard approach could be to use the BFGS optimization algorithm, the Hessian-based estimator for the standard errors, and the sample variance as option for initializing the variance process. [START_REF] Kyriakopoulou | Asymptotic normality of the QML estimator of the EGARCH(1,1) model[END_REF] As suggested by Brooks et al. (2001) there is a need for some consensus as to what forms of the asymmetric GARCH-type models are preferable, and should be estimable by all packages, as for the GARCH model. Further, a benchmark for each asymmetric GARCH model (GJR-GARCH, TGARCH, APARCH and EGARCH models) is also necessary as the benchmark for the GARCH model proposed by Fiorentini et al. (1996), namely including a specification of the parameter starting values and the initializations for the errors and conditional variances.

  Brook et al. (2001) in several ways: (i) we use updated versions of the software packages which have been modified (improved) since their study; (ii) we analyze four asymmetric GARCH-type models with Normal and Student distributions; (iii) we compare the estimation by trying to "harmonize" the options of the packages; (iii) we estimate the models on small, medium and large samples; and (iv) we compare their forecast performance by using the model confidence set (MCS) procedure proposed by Hansen et al. (2011) 5 and the one-step ahead out-of-sample volatility forecasts produced automatically by the packages and by own codes.

Ox 8 . 1 (

 81 Ox-G@RCH)(Laurent and Peters, 2001; Laurent, 2009), EViews 7, Stata 14, and two free software packages with GIG (Garch In Gretl) 2.21 package for Gretl 2018a (Lucchetti and Balietti, 2011) and rugarch 1.4-0 package for R 3.4 (R-rugarch)

  and Matlab-Econometrics compute standard ML standard-error estimates (OPG or Hessian-based estimators) whereas Ox-G@RCH and Gretl displays the Bollerslev-Wooldridge robust QML standard errors. Matlab-MFE and R-rugarch provide both standard (Hessian-based estimator) ML and robust QML standard errors.4 Comparison in estimating GARCH-type modelsIn this section, the five GARCH models previously described are fitted to represent the evolution of the volatility of two series of daily returns. The series analyzed are daily returns of the FTSE100 index observed from April 1, 2003 toDecember 31, 2007 (1195 observations). This series is downloaded from the webpage of the Oxford-Man Institute of Quantitative Finance(Heber et al., 2009).[START_REF] Ghalanos | Introduction to the rugarch package (version 1.3-8[END_REF]

  Robust QMLE denotes the Bollerslev-Wooldridge robust covariance matrix of the QML estimated parameters (sandwich estimator). OPG denotes outer product of the gradients. SQP denotes sequential quadratic programming. * denotes that the variance estimator depends on optimization algorithm since Stata switches between the BHHH and BFGS algorithms. timate symmetric and asymmetric GARCH models, i.e. a constant with GARCH(1,1)type errors.Tables4 and 5give default estimation results for the GARCH, GJR-GARCH, TGARCH, EGARCH and APARCH models. As inBrooks et al. (2001) we do not intend to 'pass judgement' on the packages that are good or not, since defaults may differ between packages in terms of initializations, starting values for parameters, convergence tolerances . . . . Rather, we only display default results for comparison with one another.

  by out-of-sample one-step ahead prediction errors. The models are re-estimated every trading day, thus aiming at a closer to real world forecast accuracy. The models are estimated using the first 4-year period, from April 1st 2003 to December 31st 2006 (T = 942 observations), and the out-of-sample forecasts of the conditional variance are calculated for the last year, from January 2, 2007 to December 31, 2007 (H = 253

  where ε2t and σ2 t are the fitted values for the squared error and the conditional variance for observation t = 1, . . . , T , with T is the last observation in the sample. One-day out-of-sample volatility forecasts σ 2 T +1 are obtained for the forecast horizon H, and compared to the realized volatility (RV) as a proxy for actual volatility. Following Liu et al. (2015) we use the 5-min calendar-time RV as RV measure, obtained from the Oxford-Man Institute of Quantitative Finance.

  However, based only on the loss function criteria, it is difficult to conclude that the forecasting performance of one model estimated from a package dominates that of the same model estimated from another one. To draw such conclusions, we use statistical tests that can provide more reliable information, namely the model confidence set (MCS) procedure proposed by Hansen et al. (2011). The MCS procedure proposed by Hansen et al. (2011) determine the set, M * , that consists of a subset of equivalent models in terms of superior predictive accuracy (SPA, Hansen, 2005) over the other competing collections of models, M 0 . The MCS procedure yields a model confidence set, M * , which is a set of models constructed

Table 3 :

 3 GARCH-type models and the econometric softwares. IGARCH: Integrated GARCH (Engle and Bollerslev, 1986); NGARCH: Nonlinear GARCH (Higgins et Bera, 1992); AVGARCH: Absolute Value GARCH (Taylor, 1986b); CGARCH: Component GARCH (Engle and Lee, 1993); ACGARCH: Asymmetric Component GARCH (Engle and Lee, 1993); fGARCH: family GARCH (Hentchel, 1995); HYGARCH: HYperbolic GARCH (Davidson, 2004); FIEGARCH: Fractionally Integrated EGARCH (Bollerslev and Mikkelsen, 1996); FIAPARCH: Fractionally Integrated APARCH (Tse, 1998); Spline-GARCH (Engle and Rangel, 2008). * : the rugarch package propose four other conditional distributions (Generalized Hyperbolic, GH; Normal Inverse Gaussian; GH Skew-Student; Johnson's reparametrized SU).

  where θ 1 (|z t | -E(|z t |)) denotes the magnitude effect, and θ 2 z t the sign effect. Thus, if θ 2 < 0 then negative innovations induce higher

	positive values. As EGARCH model does not impose any positivity restrictions on
	the volatility coefficients a sufficient condition for the stationarity of the EGARCH
	model is |β| < 1. He et al. (2002) have considered the fourth moment structure of the
	EGARCH(1,1) model.
	volatility than positive innovations of the same magnitude (asymmetry). E(|z t |) de-
	pends on the assumption made on the unconditional density of ε t . For the Normal
	distribution E(|z t |) = 2/π. The specification of the volatility in terms of its loga-
	rithmic transformation implies that the parameters in this model are not restricted to
	11 Rodríguez and Ruiz (2012) use another parametrization for the TGARCH(1,1) model, given
	by σ t = ω + α|ε t-1 | + γε t-1 + βσ t-1 . Carnero and Perez (2018) give the relationship between the
	parameters of this parametrization and those used in Matlab-MFE, Stata and Gretl (see Appendix).

Table 2 :

 2 Descriptive statistics of FTSE returns (1195 observations).

	mean (%)	Std. dev. (%)	Skewness	Kurtosis	JB test
	0.020	0.643	-0.214 *	5.278 *	148.2 *

Notes: * denotes significance at the 5% level.

Table 4 :

 4 Default estimation for volatility models with Normal distribution on large sample (1195 observations).Notes: ω * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. b means that the positivity or stationarity constraint is not satisfied.

					Packages
	Parameters		Ox-G@RCH Eviews Matlab-MFE Matlab-Ek Gauss-Fanpac R-rugarch	Stata	Gretl
	GARCH
	ω	coef.	0.005	0.005	0.005	0.005	0.005	0.005	0.005	0.005
		std.	0.003	0.002	0.003	0.002	0.002	0.003	0.002	0.003
		t-stat	1.60 *	2.65	1.81	2.41	1.99	1.54 *	2.37	1.60 *
	α	coef.	0.071	0.063	0.061	0.072	0.071	0.071	0.071	0.071
		std.	0.024	0.011	0.020	0.011	0.016	0.024	0.012	0.024
		t-stat	3.02	6.00	3.01	6.26	4.37	2.92	6.21	3.01
	β	coef.	0.917	0.923	0.926	0.917	0.917	0.917	0.917	0.917
		std.	0.027	0.014	0.023	0.015	0.020	0.029	0.015	0.028
		t-stat	33.4	68.6	39.8	62.7	46.6	31.8	62.6	33.3
	LL		-1045.0	-1039.2	-1039.3	-1045.1	-1041.5	-1045.0 -1045.0 -1045.0
	GJR-GARCH
	ω	coef.	0.005	0.003	0.003	0.005	0.005	0.005	0.005	0.005
		std.	0.003	0.001	0.002	0.002	0.002	0.003	0.002	0.003
		t-stat	1.48 *	2.99	1.73	2.44	1.96	1.42 *	2.43	1.47 *
	α	coef.	0.046	-0.001 b	0.000	0.048	0.046	0.046	0.088	0.046
	α *	coef.	0.046	-0.001	0.000	0.048	0.046	0.046	0.047	0.046
		std.	0.021	0.013	0.000	0.014	0.016	0.025	0.014	0.021
		t-stat	2.23	-0.09 *	0.84 *	3.49	2.83	1.85 * *	6.20	2.24
	β	coef.	0.921	0.952	0.955	0.921	0.921	0.921	0.921	0.921
		std.	0.028	0.012	0.013	0.015	0.020	0.031	0.014	0.029
		t-stat	32.5	81.4	74.5	63.7	47.0	30.0	64.1	32.4
	γ	coef.	0.043	0.075	0.068	0.039	0.043	0.043	-0.041	0.043
		std.	0.028	0.013	0.020	0.015	0.020	0.030	0.015	0.028
		t-stat	1.52 *	5.79	3.46	2.64	2.12	1.44 *	-2.76	1.51 *
	LL		-1042.7	-1031.1	-1030.3	-1043.1	-1039.1	-1042.7 -1042.9 -1042.7
	EGARCH
	ω	coef.	-0.939	-0.091	-0.012	-0.013	-0.013	-0.013	-0.012	-0.135
	ω *	coef.	-0.009	-0.011	-0.012	-0.013	-0.013	-0.013	-0.012	-0.012
		std.	0.333	0.021	0.006	0.007	0.008	0.009	0.007	0.034
		t-stat	-2.83	-4.36	-2.10	-1.89	-1.66 * *	-1.47 *	-1.76 * *	-4.02
	β	coef.	0.990	0.988	0.990	0.987	0.986	0.986	0.987	0.987
		std.	0.006	0.004	0.005	0.006	0.006	0.008	0.006	0.007
		t-stat	176.4	260.7	191.7	174.5	155.7	130.0	174.2	134.0
	θ 2	coef.	-0.056	-0.070	-0.061	-0.037	-0.042	-0.042	-0.042	-0.042
		std.	0.019	0.012	0.019	0.011	0.016	0.025	0.012	0.021
		t-stat	-2.95	-5.76	-3.29	-3.29	-2.67	-1.78	-3.57	-1.97
	θ 1	coef.	0.115	0.100	0.101	0.154	0.154	0.154	0.154	0.154
		std.	0.029	0.024	0.027	0.022	0.028	0.040	0.022	0.037
		t-stat	3.92	4.24	3.79	6.87	5.42	3.82	6.89	4.15
	LL		-1034.8	-1030.8	-1030.7	-1043.6	-1039.5	-1043.1 -1043.1 -1043.1

Table 5 :

 5 Default estimation for volatility models with Normal distribution on large sample (1195 observations).Notes: α * and γ * denote the parameters modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. a means no convergence. b means that the positivity or stationarity constraint is not satisfied.

					Packages
	Parameters		Ox-G@RCH Eviews Matlab-MFE Matlab-Ek Gauss-Fanpac R-rugarch	Stata	Gretl
	TGARCH
	ω	coef.	-	0.007	0.007	-	-	0.010	0.008	0.005
		std.	-	0.002	0.004	-	-	0.006	0.003	0.003
		t-stat	-	3.16	1.60 *	-	-	1.71 * *	2.38	1.69 * *
	α	coef.	-	0.055	0.024	-	-	0.099	0.105	0.082
	α *	coef.	-	0.016	0.024	-	-	0.080	0.060	0.057
		std.	-	0.013	0.015	-	-	0.028	0.015	0.021
		t-stat	-	4.39	1.65 * *	-	-	3.48	7.22	3.87
	β	coef.	-	0.945	0.944	-	-	0.906	0.923	0.923
		std.	-	0.012	0.018	-	-	0.030	0.013	0.022
		t-stat	-	79.1	51.7	-	-	30.5	71.7	42.3
	γ	coef.	-	0.705	0.064	-	-	0.196	-0.045	0.303
	γ *	coef.	-	0.078	0.064	-	-	0.039	0.045	0.050
		std.	-	0.218	0.021	-	-	0.167	0.013	0.136
		t-stat	-	3.24	3.09	-	-	1.18 *	-3.57	2.22
	LL		-	-1031.0	-1032.0	-	-	-1051.3 -1043.5 -1042.9
	APARCH
	ω	coef.	-a	0.005	0.003	-	-	0.001	0.006	0.005
		std.	-	0.002	0.002	-	-	0.003	0.003	0.003
		t-stat	-	2.33	1.35 *	-	-	0.45 *	1.93 * *	1.55 *
	α	coef.	-	0.041	0.023	-	-	0.020	0.074	0.077
		std.	-	0.022	0.034	-	-	0.047	0.017	0.026
		t-stat	-	1.82 * *	0.82 *	-	-	0.43 *	4.50	2.95
	β	coef.	-	0.949	0.954	-	-	0.939	0.922	0.922
		std.	-	0.013	0.017	-	-	0.038	0.014	0.025
		t-stat	-	75.2	55.5	-	-	24.7	66.2	36.8
	γ	coef.	-	0.769	-0.740	-	-	0.152	-0.193	0.223
		std.	-	0.449	0.753	-	-	0.091	0.074	0.170
		t-stat	-	1.72 * *	-0.98 *	-	-	1.67 * *	-2.60	1.31 *
	δ	coef.	-	1.412	1.777	-	-	3.500	1.612	1.476
		std.	-	0.368	0.558	-	-	2.400	0.497	0.630
		t-stat	-	3.83	3.18	-	-	1.46 *	3.25	2.35
	LL		-	-1030.4	-1030.2	-	-	-1040.3 -1042.8 -1042.3

Table 8 :

 8 Default estimation for volatility models with Normal distribution on medium sample (500 observations).Notes: ω * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. means no convergence. b means that the positivity or stationarity constraint is not satisfied.

					Packages
	Parameters		Ox-G@RCH Eviews Matlab-MFE Matlab-Ek Gauss-Fanpac R-rugarch Stata	Gretl
	GARCH
	ω	coef.	0.001	0.002	0.002	0.001	0.001	0.001	0.001 0.001
		std.	0.002	0.001	0.001	0.001	0.001	0.001	0.001 0.002
		t-stat	0.34 *	2.64	2.40	0.50 *	0.41 *	0.53 *	0.50 *	0.32 *
	α	coef.	0.025	0.013	0.004	0.025	0.025	0.025	0.025 0.025
		std.	0.015	0.009	0.021	0.009	0.011	0.005	0.009 0.015
		t-stat	1.66 * *	1.45 *	0.17 *	2.89	2.21	5.28		1.64 * *
	β	coef.	0.971	0.974	0.985	0.971	0.971	0.971	0.971 0.971
		std.	0.017	0.009	0.021	0.010	0.013	0.005	0.010 0.018
		t-stat	56.5	114.6	48.1	97.0	74.7	186.2	97.0	55.1
	LL		-446.9	-438.6	-436.7	-447.0	-443.5	-446.9	-447.0 -447.0
	GJR-GARCH
	ω	coef.	-a	0.002	0.001	0.001	0.001	0.001	0.001 0.001
		std.	-	0.001	0.001	0.001	0.001	0.001	0.001 0.001
		t-stat	-	3.33	1.57 *	0.69 *	0.59 *	0.73 *	0.75 *	0.45 *
	α	coef.	-	-0.047 b	0.000	0.026	0.026	0.026	-0.000 0.026
	α *	coef.	-	-0.047	0.000	0.026	0.026	0.026	0.025 0.026
		std.	-	0.005	0.000	0.008	0.011	0.009	0.012 0.015
		t-stat	-	-9.90	1.17 *	3.14	2.40	2.96	-0.03 * 1.73 * *
	β	coef.	-	1.005 b	0.982	0.981	0.981	0.982	0.982 0.980
		std.	-	0.000	0.008	0.008	0.013	0.002	0.008 0.017
		t-stat	-	11680.0	128.6	119.1	77.8	636.9	123.6	57.1
	γ	coef.	-	0.053	0.017	-0.024	-0.024	-0.025	0.025 -0.023
		std.	-	0.014	0.015	0.014	0.016	0.018	0.013 0.015
		t-stat	-	3.79	1.13 *	-1.72 * *	-1.53 *	-1.36 *	1.96	-1.52 *
	LL		-	-433.5	-433.5	-446.1	-442.6	-446.0	-445.9 -446.1
	EGARCH
	ω	coef.	-a	0.030	-0.006	-0.007	-0.007	-0.007	-0.006 -0.041
	ω *	coef.	-	-0.006	-0.006	-0.007	-0.007	-0.007	-0.006 -0.006
		std.	-	0.000	0.001	0.007	0.004	0.002	0.004 0.048
		t-stat	-	248.9	-105.4	-1.63 *	-1.59 *	-4.24	-1.69 * * -0.87 *
	β	coef.	-	0.997	0.999	0.997	0.997	0.997	0.997 0.997
		std.	-	0.002	0.000	0.003	0.004	0.001	0.003 0.005
		t-stat	-	484.5	24131.4	307.0	269.8	5742.0	302.9 220.8
	θ 2	coef.	-	-0.062	-0.053	0.020	0.021	0.021	0.021 0.022
		std.	-	0.006	0.001	0.015	0.019	0.016	0.016 0.025
		t-stat	-	-10.1	-105.9	1.29 *	1.10 *	1.35 *	1.38 *	0.88 *
	θ 1	coef.	-	-0.046	-0.043	0.047	0.045	0.044	0.045 0.044
		std.	-	0.000	0.001	0.022	0.035	0.002	0.022 0.057
		t-stat	-	-1118.6	-35.8	2.17	1.28 *	23.5	2.05	0.77 *
	LL		-	-425.7	-423.7	-446.9	-443.4	-446.8	-446.8 -446.8

* means no convergence. a * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. a means no convergence. b means that the positivity or stationarity constraint is not satisfied.

Table 11 :

 11 Default estimation for volatility models with Student distribution on large sample (1195 observations).Notes: ω * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. a means no convergence. b means that the positivity or stationarity constraint is not satisfied.

					Packages
	Parameters		Ox-G@RCH Eviews Matlab-MFE Matlab-Ek Gauss-Fanpac R-rugarch	Stata	Gretl
	TGARCH
	ω	coef.	-	0.005	0.006	-	-	0.008	0.007	0.004
		std.	-	0.003	0.004	-	-	0.006	0.004	0.003
		t-stat	-	2.18	1.55 *	-	-	1.37 *	1.79 * *	1.52 *
	α	coef.	-	0.047	0.018	-	-	0.085	0.093	0.072
	α *	coef.	-	0.010	0.018	-	-	0.067	0.053	0.050
		std.	-	0.014	0.013	-	-	0.027	0.018	0.019
		t-stat	-	3.43	1.39 *	-	-	3.14	5.12	3.92
	β	coef.	-	0.954	0.952	-	-	0.920	0.932	0.933
		std.	-	0.013	0.015	-	-	0.029	0.015	0.020
		t-stat	-	72.3	61.8	-	-	31.3	60.2	47.8
	γ	coef.	-	0.783	0.061	-	-	0.211	-0.040	0.303
	γ *	coef.	-	0.074	0.061	-	-	0.036	0.040	0.044
		std.	-	0.320	0.018	-	-	0.177	0.017	0.132
		t-stat	-	2.45	3.43	-	-	1.20 *	-2.30	2.29
	ν	coef.	-	13.0	12.9	-	-	10.5	11.8	11.9
		std.	-	4.42	3.80	-	-	2.61	3.71	3.34
		t-stat	-	2.93	3.39	-	-	4.01	3.19	3.57
	LL		-	-1025.0	-1025.1	-	-	-1040.8 -1035.7 -1035.4
	APARCH
	ω	coef.	-a	0.003	0.002	-	-	0.001	0.003	0.003
		std.	-	0.019	0.001	-	-	0.003	0.002	0.003
		t-stat	-	1.61 *	1.55 *	-	-	0.25 *	1.24 *	1.18 *
	α	coef.	-	0.023	0.014	-	-	0.014	0.051	0.055
		std.	-	2.065	0.011	-	-	0.067	0.023	0.035
		t-stat	-	0.01 *	1.29 *	-	-	0.21 *	2.21	1.57 *
	β	coef.	-	0.956	0.960	-	-	0.954	0.937	0.935
		std.	-	0.014	0.010	-	-	0.034	0.016	0.024
		t-stat	-	70.6	94.8	-	-	27.8	57.7	38.3
	γ	coef.	-	0.999	-0.992	-	-	0.171	-0.162	0.182
		std.	-	105.8	0.051	-	-	0.109	0.109	0.158
		t-stat	-	0.01	-19.5	-	-	1.57 *	-1.48 *	1.15 *
	δ	coef.	-	1.730	2.098	-	-	3.500	2.090	1.932
		std.	-	0.561	0.567	-	-	4.555	0.805	1.032
		t-stat	-	3.08	3.70	-	-	0.77 *	2.60	1.87 * *
	ν	coef.	-	12.2	11.7	-	-	11.0	11.2	11.3
		std.	-	3.95	3.30	-	-	2.92	3.38	3.09
		t-stat	-	3.09	3.53	-	-	3.78	3.31	3.67
	LL		-	-1023.4	-1021.9	-	-	-1030.8 -1034.4 -1034.1

Table 12 :

 12 Estimation for volatility models with Normal distribution on large sample (1195 observations) from BFGS algorithm and Hessian-based estimator.Notes: ω * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. a means no convergence. b means that the positivity or stationarity constraint is not satisfied. c means failure to improve likelihood (non-zero gradients).

					Packages
	Parameters		Ox-G@RCH Eviews Matlab-MFE Matlab-Ek Gauss-Fanpac R-rugarch	Stata	Gretl
	GARCH
	ω	coef.	0.005	0.005	0.005	-	0.005	0.005	0.005	0.005
		std.	0.002	0.002	0.002	-	0.002	0.002	0.002	0.002
		t-stat	1.99	1.98	2.01	-	1.99	2.12	1.98	1.98
	α	coef.	0.071	0.071	0.072	-	0.071	0.072	0.071	0.071
		std.	0.016	0.016	0.016	-	0.016	0.016	0.016	0.016
		t-stat	4.37	4.36	4.37	-	4.37	4.50	4.36	4.36
	β	coef.	0.917	0.917	0.917	-	0.917	0.917	0.917	0.917
		std.	0.020	0.020	0.020	-	0.020	0.019	0.020	0.020
		t-stat	46.5	46.5	46.5	-	46.6	48.2	46.5	46.5
	LL		-1045.0	-1045.0	-1045.1	-	-1041.5	-1045.2 -1045.0 -1045.0
	GJR-GARCH
	ω	coef.	0.005	0.005	0.004	-	0.005	-a	0.005	0.005
		std.	0.002	0.002	0.002	-	0.002	-	0.002	0.002
		t-stat	1.95 * *	1.94 * *	1.89 * *	-	1.96	-	1.94 * *	1.94 * *
	α	coef.	0.046	0.047	0.043	-	0.046	-	0.088	0.046
	α *	coef.	0.046	0.047	0.043	-	0.046	-	0.047	0.046
		std.	0.016	0.016	0.016	-	0.016	-	0.021	0.016
		t-stat	2.82	2.87	2.62	-	2.83	-	4.10	2.83
	β	coef.	0.921	0.921	0.925	-	0.921	-	0.921	0.921
		std.	0.020	0.020	0.020	-	0.020	-	0.020	0.020
		t-stat	46.6	46.6	47.4	-	47.0	-	46.5	46.5
	γ	coef.	0.043	0.041	0.042	-	0.043	-	-0.041	0.043
		std.	0.020	0.020	0.020	-	0.020	-	0.020	0.020
		t-stat	2.12	2.04	2.15	-	2.12	-	-2.03	2.10
	LL		-1042.7	-1042.9	-10341.4	-	-1039.1	-	-1042.9 -1042.7
	EGARCH
	ω	coef.	-0.939	-0.135	-0.013	-	-0.013	-0.014	-0.012	-0.135
	ω *	coef.	-0.009	-0.012	-0.013	-	-0.013	-0.014	-0.012	-0.012
		std.	0.275	0.027	0.007	-	0.008	0.008	0.008	0.027
		t-stat	-3.42	-5.11	-1.72 * *	-	-1.66 * *	-1.77 * *	-1.64 * *	-5.11
	β	coef.	0.990	0.987	0.987	-	0.986	0.986	0.986	0.986
		std.	0.005	0.006	0.006	-	0.006	0.006	0.006	0.006
		t-stat	203.4	155.1	160.8	-	155.7	153.2	155.0	154.5
	θ 2	coef.	-0.056	-0.042	-0.038	-	-0.042	-0.037	-0.042	-0.042
		std.	0.015	0.016	0.015	-	0.016	0.025	0.016	0.016
		t-stat	-3.84	-2.66	-2.53	-	-2.67	-1.49 *	-2.65	-2.66
	θ 1	coef.	0.115	0.154	0.146	-	0.154	0.154	0.154	0.154
		std.	0.026	0.029	0.028	-	0.028	0.039	0.029	0.029
		t-stat	4.50	5.39	5.26	-	5.42	3.96	5.39	5.39
	LL		-1034.8	-1043.1	-1043.3	-	-1039.5	-1043.7 -1043.1 -1043.1

Table 14 :

 14 Out-of-sample results of GARCH-type models with a Normal distribution. Notes: MCS-MSE and MCS-MAD denote the p-value of the statistic T max M of Hansen et al. (2011) based on the MSE and QLIKE loss functions respectively. * and * * mean that the forecasts are in M * 20% and M * 50% , respectively. The lowest MSE and QLIKE values appear in bold.

	Packages

Table 15 :

 15 Out-of-sample results of GARCH-type models with a Student distribution. Notes: MCS-MSE and MCS-MAD denote the p-value of the statistic T max M of Hansen et al. (2011) based on the MSE and QLIKE loss functions respectively. * and * * mean that the forecasts are in M * 20% and M * 50% , respectively. The lowest MSE and QLIKE values appear in bold.

	Packages

Table 16 :

 16 Out-of-sample results of GARCH-type models with BFGS algorithm. Notes: MCS-MSE and MCS-MAD denote the p-value of the statistic T max M of Hansen et al. (2011) based on the MSE and QLIKE loss functions respectively. * and * * mean that the forecasts are in M * 20% and M * 50% , respectively. The lowest MSE and QLIKE values appear in bold.

	Packages

In this paper we do not deal with detailed descriptions of the packages, their user-friendliness, flexibility, speed . . . See Brooks (1997) for these issues.

[START_REF] Brooks | Benchmarks and the accuracy of GARCH model estimation[END_REF] is the only study which compares the out-of-sample forecasts of the conditional variance produced by a GARCH model estimated from different packages. They only compare the forecast values from the packages for one to eight step ahead forecasts.

Note that the OLS estimate of c is r.

Under the assumption of Normal distribution k = 3 and thus the condition becomes 3α 2 + 2αβ + β 2 < 1. SeeLing and McAleer (2002a, 2002b) for other distributions.

Under a Normal distribution and a Student-t(ν) distribution, with ν > 5, δ = 1 2 . SeeLing and McAleer (2002a, 2002b) for other distributions.

Note that Gretl provides estimation of GJR-GARCH parameters for both parametrizations.

Note that some packages specify the GJR-GARCH (Gretl) and TGARCH (Eviews, Gretl and Rrugarch) models from the parametrization of the APARCH model

.[START_REF] Carnero | Identification of asymmetric conditional heteroscedasticity in the presence of outliers[END_REF] See[START_REF] Zivot | Practical issues in the analysis of univariate GARCH models[END_REF] and Xekalaki and Degiannakis (2010) for a discussion on the ML estimations for the GARCH-type models.

The asymptotic properties, in particular the asymptotic normality, of the quasi-maximum likelihood estimator (QMLE) hold under mild conditions for GARCH models (e.g.,[START_REF] Lee | Asymptotic theory for the GARCH(1,1) quasi maximum likelihood estimator[END_REF] Berkes et al., 2003;[START_REF] Francq | Maximum likelihood estimation of pure GARCH and ARMA-GARCH processes[END_REF]) and for GJR-GARCH models(Hamadeh and Zakoian, 2011). However, the statistical properties for the QMLE of the EGARCH(1,1) parameters are not available under general conditions, but rather only for special cases under highly restrictive conditions[START_REF] Wintenberger | Continuous invertibility and stable QML estimation of the EGARCH(1,1) model[END_REF][START_REF] Kyriakopoulou | Asymptotic normality of the QML estimator of the EGARCH(1,1) model[END_REF][START_REF] Martinet | On the invertibility of EGARCH(p,q)[END_REF].

An overview of SQP methods can be found in[START_REF] Nocedal | Numerical Optimization[END_REF] and[START_REF] Fletcher | Practical Methods of Optimization[END_REF].

Note that, although the Fisher information matrix is known to be block diagonal, the efficient implementation of Newton or Newton-like algorithms requires the full Hessian matrix, which is not block-diagonal in small samples[START_REF] Calzolari | Alternative estimators of the covariance matrix in GARCH models[END_REF].

The QMLE is generally close to the exact MLE for symmetric departures from Normality.However, Engle and Gonzales-Rivera (1991) showed that the loss in efficiency may be quite high for non-symmetric distributions.

The Oxford MFE Toolbox is the follow on to the UCSD_GARCH toolbox.

Another package is available for GARCH modelling in R with fGarch, which is a part of the Rmetrics suite. However, this package only estimates GARCH and APARCH models.

By default, Stata performs 5 iterations with the BHHH algorithm followed by 15 iterations with the BFGS algorithm, and then switch back to BHHH for five iterations, and so on.Lombardi and Gallo (2002) show that the use of a mixed algorithm can provide some gains in terms of iterations and RMSE.

We use the following optimization packages: the maximization package maxlikmt2.0 for Gauss 17; the maximization package fminunc for Matlab-MFE; the optimization toolbox fmincon for Matlab-Econometrics, where the function solves a quadratic programming subproblem at each iteration, and updates an estimate of the Hessian of the Lagrangian at each iteration using the BFGS formula; the SQP solver solnp of[START_REF] Ye | Interior point algorithms: Theory and analysis[END_REF] based on an augmented Lagrange multiplier method with a SQP interior algorithm, implemented in R byGhalanos and Theussl (2011).

[START_REF] Ghalanos | Introduction to the rugarch package (version 1.3-8[END_REF] The data are available on the following website: http://realized.oxford-man.ox.ac.uk/

data.[START_REF] Ghalanos | Rsolnp: General non-linear optimization using augmented Lagrange multiplier method[END_REF] It is well known that the outliers may pose difficulties for the identification and estimation of GARCH models governing the conditional volatility of returns (e.g., Franses and Ghijsels, 1999;[START_REF] Carnero | Effects of outliers on the identification and estimation of the GARCH models[END_REF][START_REF] Carnero | Estimating GARCH volatility in the presence of outliers[END_REF]

, 2016).[START_REF] Glosten | On the relation between the expected value and the volatility of the nominal excess return on stocks[END_REF] [START_REF] Carnero | Estimating GARCH volatility in the presence of outliers[END_REF] and[START_REF] Rodriguez | Revisiting several popular GARCH models with leverage effect: Differences and similarities[END_REF] advocate this approach to detect outliers in

Note that Ox-G@rch does not estimate this model due to problem of convergence.

By analyzing different asymmetric models[START_REF] Rodriguez | Revisiting several popular GARCH models with leverage effect: Differences and similarities[END_REF] show that when the parameters satisfy the positivity, stationarity, and finite kurtosis conditions, the dynamics that the GJR-GARCH model can represent are heavily limited while those of the EGARCH and TGARCH models are less restricted.

Note that all the GARCH-type models with a Student distribution are rejected. The results are given in the online Appendix.

It would be interesting to determine a benchmark for the parametrization of the asymmetric models[START_REF] Engle | A permanent and transitory component model of stock return volatility[END_REF]

The MCS p-values are calculated using bootstrap implementation with 10,000 resamples[START_REF] Hansen | Model confidence sets for forecasting models[END_REF]. The MCS test is carried out using the Ox software package MULCOM ofHansen and Lunde (2007).

The APARCH models have been not rejected for only two packages with a Normal distribution and no package for the Student one.

A comparison from a Monte Carlo experiment by trying various choices of initial values, optimization routines and standard errors computation methods would be informative but it is beyond the scope of this study. Future research is encouraged to address these issues.

Notes: ω * denotes the intercept term modified for comparison. * and * * mean significant at the 5% and 10% level, respectively. a means no convergence. b means that the positivity or stationarity constraint is not satisfied.