In finance and econometrics, important information about past market movements is modelled through the conditional distribution of return data series using the autoregressive moving average (ARMA) models. However, in these models, the conditional distribution is assumed to be homoskedastic which poses challenges in the modelling and analysis of returns with time-varying volatility and volatility clustering, a commonly observed phenomena. This led to the development of autoregressive conditional heteroskedastic (ARCH) and generalized ARCH (GARCH) models introduced by Engle  and Bollerslev  respectively. In empirical finance, the most important and widely used model is the combination of ARMA with GARCH, referred to as an ARMA-GARCH model (see , ).
In general, the noise term of ARMA, GARCH and ARMA-GARCH models is assumed to be normal or Student-
distribution, where the normal distribution is known for the desirable property of stability but fails to capture the heavy-tailedness of the data. On the other hand, the Student-
distribution allows for heavier tails, but, when compared with the normal, lacks the desirable property of stability. These flaws, therefore, motivate us to use the family of stable distributions as a model for unconditional, conditional homoskedastic and conditional heteroskedastic return series distribution. Some of the significant and attractive features of stable distributions, apart from stability are heavy-tailedness, leptokurtic shape, domains of attraction and skewness. For more details on stable distributions, see. Hence, there is a need to explore the behaviour of ARMA, GARCH and ARMA-GARCH models with stable noise for effective modelling of the return data series which also involves estimation of the parameters of these models, and only a handful of estimation techniques based on stable noise are available. For details on the different estimation techniques used for such models, see . Therefore, we feel that this article is an important contribution to the literature available on these models.
In this article, we develop a method for estimating the parameters of ARMA model with symmetric stable noise and symmetric stable GARCH noise by modifying the Hannan-Rissanen Method . For the estimation of the GARCH parameters, we propose the modified empirical characteristic function method which is based on the method discussed in . The efficiency and effectiveness of the two proposed methods is validated through Monte Carlo simulation. For comparative analysis, we compare the modified Hannan-Rissanen method with the first two of the three specifc M-estimators of the parameters of ARMA models with symmetric stable noise introduced by Cadler and Davis , namely, the least absolute deviation (LAD) estimator; the least squares (LS) estimator, and the maximum likelihood (ML) estimator. Amongst the three M-estimators, LAD is an excellent choice for estimating ARMA parameters with stable noise due to its robustness, computational simplicity and good asymptotic properties. Our method works at par with the LAD estimator both in terms of accuracy and computational simplicity.
The paper is organised as follows. Section 2 gives a brief introduction to the stable distributions along with the necessary definitions and notations. Section 3 discusses the two new methods proposed for the estimation of ARMA, GARCH and ARMA-GARCH parameters. Section 4 deals with simulations and comparative analysis of the proposed methods. Section 5 discusses application of our proposed methods to the financial data. Finally, Section 6 gives some concluding remarks on the proposed methods.
2 Preliminaries and Notations
Stable distributions. These distributions form a rich class of heavy-tailed distributions, introduced by Paul Lévy 
, in his study on the Generalized Central Limit Theorem. Each distribution, in this class, is characterized by four parameters, namely, , , and , which, respectively, denote the index of stability, skewness, scale and shift of the distribution. Their respective ranges are given by , , and . For more details, see . In this paper, we deal with symmetric stable distributions. We say that the distribution is symmetric around zero if and only if and denote by . The characteristic function (cf) representation is given by
For simulations, we assume for distribution throughout.
Next, we discuss some important definitions and notations required for parameter estimation of ARMA, GARCH and ARMA-GARCH models. Let and
denote the parameter vectors used in the ARMA model. For the GARCH model, the parameter vectors used are, , and .
Let the time series be an ARMA() process with noise, denoted by -ARMA(), if
where denotes the order of autoregression, denotes the order of moving average and
denotes the noise sequence of iid random variables withdistribution with (cf) given in (1) for and . Further, in the unit disc , and have no common zeros. For more details on ARMA() process with noise, (see , pp. 376-380).
Let the time series be a GARCH() process with noise denoted by -GARCH(), if
where , and are non-negative and denotes the noise sequence of iid random variables with distribution with (cf) given in (1). It is known that the -GARCH() process has a unique strictly stationary solution if while -GARCH() process with or has a unique strictly stationary solution if and . For further details and information on the approximate values of used for different values of , see .
We say that the time series is an ARMA() process with -GARCH() noise, denoted by ARMA()--GARCH(), if
where , and are non-negative and denotes the noise sequence of iid random variables with distribution with (cf) given in (1) for and .
Note now that due to undefined covariance when , the classical autocovariance function cannot be considered as a tool for developing methods of parameter estimation of the processes defined in (2) and (4). Thus, autocovariation (normalized) and autocodifference serves as a measure of dependence best suited for heavy-tailed distributions and processes. In our proposed method of estimating the parameters of the processes defined in (2) and (4), we make use of the normalized autocovariation. The normalized autocovariation (see ) for process with lag is given by
and its estimator for a sample being a realization of a stationary process is given by
where and .
3 Parameter Estimation
We propose two methods for the estimation of the parameters , , and for the processes defined in (2), (3) and (4). For estimation of and , we suitably modify the Hannen-Rissanen method  and for estimating , , and we modify the empirical characteristic function method discussed in .
3.1 Modified Hannan-Rissanen Method (MHR)
Let be an AR() process given by
where constitutes sample of an iid random variables with .
Multiply both sides of (7) by the vector where and take the expectation to obtain equations of the form
Divide the th equation of (8) by to obtain
Finally, in order to estimate the values of the parameter vector , the matrix should be nonsingular and make use of the sample normalized autocovariation . Thus,
Using the obtained estimated coefficients , compute the estimated residuals from
Next, we estimate the vector of parameters, by least absolute deviation regression of onto by minimizing with respect to where,
The proposed MHR method is both computationally efficient and nearly optimal in estimating the ARMA coefficients of (2) or (4). It is important to note that the MHR method does not require information (or estimation) of the parameters of noise distribution which is also true for LAD method. The effectiveness of the proposed method in comparison to other methods is shown in Section 4.
3.2 Modified Empirical Characteristic Function Method (MECF)
We employ this method to obtain the estimates of , and for processes defined in (3) and (4). For simplicity, we shall discuss the method for the process defined in (3), for the case and relevant in empirical finance. The method can be extended to higher orders of and in a similar spirit. For and we get
where , and are non-negative satisfying the stationary condition and denotes the noise sequence of iid random variables with distribution. The parameter estimates for (10) are obtained by minimizing the function over , and defined as
where , , where are iid random variables with distribution and is the estimate of obtained from the noise sequence of iid random variables with distribution using the hybrid method discussed in . We make use of the function “optim” available in “stats” package of R for the minimization of the function over , and . The method used is “L-BFGS-B” introduced by Byrd et al. .
4 Simulation and Comparative Analysis
For the comparative analysis, we investigate the performance of LAD, LS and MHR for the estimation of the parameters of the processes defined in (2) and (4). The parameters of the process defined in (3) are obtained using MECF discussed in Section 3.2. We consider four models to check the efficiency of our proposed method for (2) and (4) in comparison to LAD and LS through Monte Carlo simulations.
-ARMA(1,1): , and
MA(1)--GARCH(1,1): , , with such that .
ARMA(1,1)--GARCH(1,1): , , ,
with such that .
The noise sequence is considered to be iid , where and , and are fixed to 0, 1, 0 respectively and is generated using the function “rstable” in the “stabledist” package of R. For a selected set of values , and , simulations are run where 1000 realisations of each model of length 1000 are generated. Tables 4 and 4 show the efficiency and effectiveness of our proposed MHR method over LAD and LS for models (M1) and (M2) respectively for different values of , and in terms of mean and root mean squared error (RMSE). The estimate of for model (M1) and (M2) is obtained from the estimated residuals given in Step 2 of Subsection 3.1 using the hybrid method discussed in  .
Tables 4 and 4 show the accuracy and efficiency of the proposed MHR method and the MECF used in obtaining the estimate of and respectively for model (M3). The estimate of for this model is obtained from the noise sequence using the hybrid method.
Tables 5 and 6 show the accuracy and efficiency of the proposed MHR method and the MECF used in obtaining the estimate of and respectively for model (M4). The hybrid method is employed in obtaining the estimate of from the noise sequence .
We observe that for all the four models (M1), (M2), (M3) and (M4), the least squares method (LS) performs the worst while the least absolute deviation (LAD) and the proposed modified Hannan-Rissanen method (MHR) are at par with each other in terms of the robustness, accuracy and efficiency of the estimates. The estimates of the GARCH parameters obtained via MECF are also precise and accurate.
5 Application to Financial Data
The time series real dataset that we have considered is the International Business Machines Corporation (IBM) obtained from Yahoo Finance for the period January 19, 2000 - March 19, 2005, comprising of 1297 daily log returns values of the adjusted closing price. From Figure 1
(a) , we observe the presence of time varying volatility of the log returns of the time series with mean almost zero and the minimum and maximum approximately around zero. The kurtosis is 5.8164 while the skewness is -0.030 which implies the distribution of the log returns is fairly symmetrical with heavier tails than the normal. These characteristics thereby suggest the application of ARMA()--GARCH() model to the given dataset. To determine whether the considered time series data is stationary, we implement the Augmented Dickey-Fuller Test (ADF) “adf.test” available in the “tseries” package in R. The -value obtained is 0.01 which confirms the stationarity of the data. We first try to implement the ARMA() model and check whether it fits the dataset well. In order to select the best order for fitting an ARMA() model according to either AIC (Akaike Information Criterion) or BIC (Bayesian Information Criterion), we make use of “auto.arima” available in the package “forecast” in R. The best model chosen by “auto.arima” was MA(1) model and we assume that the residuals obtained from this model constitute a sample from distribution. Thus, we can employ the MHR method to obtain the estimate .
In order to check if the residuals can be considered as an independent sample, we make use of the empirical autocovariation function instead of the classical autocovariance function. From Figure 1 (c), one can observe that the dependence in the residual series is almost unidentifiable. Finally, we prove that the distribution of the residual series is using the Kolmogorov-Smirnov (KS) test and estimate the parameters (, , , ) using the method in . The estimates obtained for the residual series are (, , , )=
. The KS test statistics is defined as:
denote the empirical cumulative distribution function for Sample 1 and Sample 2 respectively, both of length. The test derived from the KS statistics is called the two sample KS test. In our case, Samples 1 and 2 are the residuals and the random sample simulated from distribution with estimated parameters corresponding to the residuals respectively. Finally, we obtain 100 of the test and create a boxplot as shown in Figure 1 (d). Assuming 0.05 confidence level, we fail to reject the hypothesis that the samples are from the same distribution.
However, from the residuals plot of MA(1) with noise, we observe that the volatility decays with time, thereby suggesting the replacement of noise with -GARCH(1,1) noise in the -MA(1) model. Therefore, we will build MA(1)--GARCH(1,1) model. For the given time series, the initial values of , and is considered to be 0.001, 0.1 and 0.8 respectively to obtain the estimates. The GARCH(1,1) estimated parameters obtained after employing MECF are given in Table 7.
In order to evaluate if the considered model is correctly specified, we analyse the estimated standardized residuals using graphical diagnostics. We plot the standardized residuals to check the absence of conditional hetroskedascity. From Figure 1
(e), we observe that the standardized residual plot seems quite stable, indicating a constant variance over time. Further,we make use of the QQ plot to check if the distribution of the standardized residuals matches the assumed noise distribution (used in the estimation. From Figure 1 (f), we observe that the distribution of of the standardized residuals is almost ( with parameters (, , , )=.
6 Concluding Remarks
To conclude, we make the following observations in relation to our proposed methods.
The M-estimators namely, the LAD and ML perform well when the noise is heavy-tailed (stable). However, the LAD estimator is preferred over the ML as it is computationally efficient, robust and does not require information (or estimation) of the parameters of the noise distribution. The proposed MHR method also inherits the same advantages of LAD especially for , and . For details on the limitations of LS and ML estimator, see .
For the estimation of the GARCH coefficients for GARCH() process with noise defined in (3) and ARMA() process with -GARCH() defined in (4), we introduce and discuss MECF for the case and for simplicity. The method is computationally efficient, robust and can further be extended to higher orders of and .
Finally, we give an application of our proposed methods, where the noise distribution of the dataset is considered to be and later -GARCH due to volatility clustering. From Table 7, we observe the estimates are statistically significant and Figure 1 shows efficient modelling of the financial data using our proposed methods through residual analysis.
|0.0100 (0.0152)||0.0096 (0.0283)||0.0098 (0.0153)||1.5011 (0.0529)|
|0.0496 (0.0175)||0.0490 (0.0258)||0.0492 (0.0177)||1.5553 (0.0532)|
|0.0997 (0.0216)||0.0991 (0.0269)||0.0990 (0.0218)||1.6552 (0.0524)|
|0.3996 (0.0310)||0.3992 (0.0295)||0.3965 (0.0312)||1.8523 (0.0438)|
|0.4987 (0.0362)||0.4990 (0.0306)||0.4911 (0.0372)||1.9500 (0.0310)|
|0.0508 (0.0247)||0.0511 (0.0604)||0.0464 (0.0225)|
|0.1975 (0.0349)||0.1961 (0.0653)||0.1925 (0.0354)|
|0.3003 (0.0384)||0.3014 (0.0577)||0.2989 (0.0382)|
|0.3952 (0.0424)||0.3935 (0.0449)||0.3893 (0.0405)|
|0.5018 (0.0373)||0.5074 (0.0376)||0.4942 (0.0367)|
|True Values (|
|0.0085 (0.0037)||0.0196 (0.0035)||0.7086 (0.0275)||1.2011 (0.0471)|
|0.0499 (0.0032)||0.0399 (0.0032)||0.9060 (0.0246)||1.4017 (0.0511)|
|0.0121 (0.0879)||0.0498 (0.0032)||0.8046 (0.0234)||1.6548 (0.0552)|
|0.4966 (0.0335)||0.0607 (0.0033)||0.8102 (0.0272)||1.7584 (0.0496)|
|1.0507 (0.3165)||0.0301 (0.0033)||0.9115 (0.0267)||1.8574 (0.0433)|
|0.1004 (0.0497)||0.4028 (0.0422)||0.1013 (0.0850)||0.4008 (0.0764)||0.0873 (0.0476)||0.4126 (0.0438)|
|0.0973 (0.0509)||0.5979 (0.0415)||0.0917 (0.0763)||0.5977 (0.0584)||0.0921 (0.0510)||0.6015 (0.0419)|
|0.2007 (0.0568)||0.4944 (0.0417)||0.1929 (0.0932)||0.5025 (0.0651)||0.1920 (0.0564)||0.5007 (0.0456)|
|0.1957 (0.0505)||0.8971 (0.0155)||0.1965 (0.0683)||0.8964 (0.0192)||0.1912 (0.0482)||0.8964 (0.0172)|
|0.2967 (0.0436)||0.8005 (0.0219)||0.2995 (0.0420)||0.7970 (0.0228)||0.2935 (0.0428)||0.7998 (0.0280)|
|True Values (|
|0.0115 (0.0042)||0.0196 (0.0040)||0.7096 (0.0305)||1.5572 (0.0554)|
|0.0573 (0.0035)||0.0438 (0.0085)||0.9141 (0.0302)||1.6504 (0.0573)|
|0.1090 (0.0879)||0.0500 (0.0032)||0.8132 (0.0234)||1.7036 (0.0471)|
|0.5016 (0.0308)||0.0605 (0.0030)||0.8133 (0.0263)||1.7542 (0.0495)|
|1.0622 (0.3311)||0.0301 (0.0032)||0.9142 (0.0282)||1.8507 (0.0485)|
-  Engle, R.F. (1982). Autoregressive conditional heteroscedemticity with estimates of the variance of U.K. inflation, Econometrica 50, 987-1008.
-  Bollerslev (1986). Generalized autoregressive conditional heterockedasticity. Journal of Econometrics 31, 307-327 .
-  Cadler M. and Davis R.A. (1998). Inference for linear processes with stable noise. A practical guide to heavy tails: Birkhauser Boston Inc., Cambridge, MA, USA, 159-176.
-  Mittnik S. and Paolella M. S. (2003). Prediction of Financial Downside-Risk. Handbook of Heavy Tailed Distributions in Finance (S. Rachev, eds.), Elsevier, Amsterdam.
Gallagher, C.M. (2001). A method for fitting stable autoregressive models using the autocovariation function.
Statistics and Probability Letters53, 381-390.
-  Omelchenko, V. (2009). Parameter Estimation of the Stable GARCH(1,1)-Model. WDS’09 Proceedings of Contributed Papers, Part I, 137-142.
-  Kruczek, P., Wylomanska, A., Teuerle, M., and Gajda, J. (2017). The modified Yule-Walker method for -stable time series models. Physica A: Statistical Mechanics and Its Applications 469, 588-603.
-  Sathe, Aastha M. and Upadhye, Neelesh S. (2019). Estimation of the Parameters of Multivariate Stable Distributions. stat.CO, arXiv:1902.09796 .
-  Byrd, R. H., Lu, P., Nocedal, J. and Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16, 1190–1208.
-  Lévy, P. (1924). Th‘eorie des erreurs la loi de Gauss et les lois exceptionelles, Bulletin de la Soci‘et‘e de France 52:49-85.
-  Hannan, E. J. and Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving order. Biometrika 69(1):81-94.
-  Nolan, J. P. (2014). Financial modeling with heavy-tailed stable distributions. WIREs Computational Statistics 6:45-55.
-  Samorodnitsky, G. and Taqqu, M. S. (1994). Stable Non-Gaussian Random Processes, Stochastic Models with Infinite Variance. Chapman and Hall, New York, 632 pages.
-  Mikosch, Thomas, et al. (1995). Parameter Estimation for ARMA Models with Infinite Variance Innovations. The Annals of Statistics, 23(1),305-326.
-  Zhang, M., Li, C. (2007). Application of GARCH model in the forecasting of day-ahead electricity prices. Third International Conference on Natural Computation, 99-103.
-  Zhou, B., He, D. and Sun,Z. (2006). Traffic predictability based on ARIMA/GARCH model. Next Generation Internet Design and Engineering,199-207.