ARMA Models with GARCH/APARCH Errors

41 Pages·2006·1.07 MB·English

by  Diethelm Würtz

Journal of Statistical Software
JSS
MMMMMM YYYY, Volume VV, Issue II. http://www.jstatsoft.org/
Parameter Estimation of ARMA Models with
GARCH/APARCH Errors
An R and SPlus Software Implementation
Diethelm Wu¨rtz1, Yohan Chalabi2, Ladislav Luksan3
1,2Institute for Theoretical Physics
Swiss Federal Institute of Technology, Zurich
3Academy of Sciences of the Czech Republic
Institute of Computer Science, Praha
Abstract
We report on concepts and methods to implement the family of ARMA models with
GARCH/APARCH errors introduced by Ding, Granger and Engle. The software imple-
mentation is written in S and optimization of the constrained log-likelihood function is
achieved with the help of a SQP solver. The implementation is tested with Bollerslev’s
GARCH(1,1) model applied to the DEMGBP foreign exchange rate data set given by
Bollerslev and Ghysels. The results are compared with the benchmark implementation
of Fiorentini, Calzolari and Panattoni. In addition the MA(1)-APARCH(1,1) model for
the SP500 stock market index analyzed by Ding, Granger and Engle is reestimated and
compared with results obtained from the Ox/G@RCH and SPlus/Finmetrics software
packages. The software is part of the Rmetrics open source project for computational
finance and financial engineering. Implementations are available for both software envi-
ronments, R and SPlus.
Keywords: TimeSeriesAnalysis,Heteroskedasticity,ARCH,GARCH,APARCH,constrained
maximum log-likelihood, SQP solver, R, Rmetrics, SPlus, Finmetrics, Ox, G@RCH.
1. Introduction
GARCH, Generalized Autoregressive Conditional Heteroskedastic, models have become im-
portant in the analysis of time series data, particularly in financial applications when the goal
is to analyze and forecast volatility. For this purpose, we describe functions for simulating,
estimating and forecasting various univariate GARCH-type time series models in the condi-
tional variance and an ARMA specification in the conditional mean. We present a numerical
2 An R and SPlus Software Implementation
implementation of the maximum log-likelihood approach under different assumptions, Nor-
mal, Student-t, GED errors or their skewed versions. The parameter estimates are checked
by several diagnostic analysis tools including graphical features and hypothesis tests. Func-
tions to compute n-step ahead forecasts of both the conditional mean and variance are also
available.
The number of GARCH models is immense, but the most influential models were the first.
Beside the standard ARCH model introduced by Engle [1982] and the GARCH model in-
troduced by Bollerslev [1986], we consider also the more general class of asymmetric power
ARCH models, named APARCH, introduced by Ding, Granger and Engle [1993]. APARCH
models include as special cases the TS-GARCH model of Taylor [1986] and Schwert [1989],
the GJR-GARCH model of Glosten, Jaganathan, and Runkle [1993], the T-ARCH model of
Zakoian [1993], the N-ARCH model of Higgins and Bera [1992], and the Log-ARCH model of
Geweke [1986] and Pentula [1986].
Coupled with these models was a sophisticated analysis of the stochastic process of data gen-
erated by the underlying process as well as estimators for the unknown model parameters.
Theorems for the autocorrelations, moments and stationarity and ergodicity of GARCH pro-
cesseshavebeendevelopedformanyoftheimportantcases. Thereexistacollectionofreview
articlesbyBollerslev, ChouandKroner[1992], BeraandHiggins[1993], Bollerslev, Engleand
Nelson [1994], Engle [2001], Engle and Patton [2001], and Li, Ling and McAleer [2002] that
give a good overview of the scope of the research.
2. Mean and Variance Equation
We describe the mean equation of an univariate time series x by the process
t
x = E(x |Ω )+ε , (1)
t t t−1 t
whereE(·|·)denotestheconditionalexpectationoperator,Ω theinformationsetattimet−
t−1
1, andε theinnovationsorresidualsofthetimeseries. ε describesuncorrelateddisturbances
t t
withzeromeanandplaystheroleoftheunpredictablepartofthetimeseries. Inthefollowing
we model the mean equation as an ARMA process, and the innovations are generated from a
GARCH or APARCH process.
ARMA Mean Equation: The ARMA(m,n) process of autoregressive order m and moving
average order n can be described as
m n
X X
x = µ + a x + b ε +ε ,
t i t−i j t−j t
i=1 j=1
= µ + a(B)x + b(B)ε , (2)
t t
with mean µ, autoregressive coefficients a and moving average coefficients b . Note, that
i i
the model can be expressed in a quite comprehensive form using the backshift operator B
defined by Bx = x . The functions a(B) and b(B) are polynomials of degree m and n
t t−1
respectively in the backward shift operator B. If n = 0 we have a pure autoregressive process
and if on the other hand m = 0 we have a pure moving average process. The ARMA time
Journal of Statistical Software 3
series is stationary when the series a(B), which is the generating function of the coefficients
a , converges for |B| < 1 that is, on or within the unit circle.
i
GARCH Variance Equation: The mean equation cannot take into account for heteroskedas-
tic effects of the time series process typically observed in form of fat tails, as clustering of
volatilities, and the leverage effect. In this context Engle [1982] introduced the Autoregres-
sive Conditional Heteroskedastic model, named ARCH, later generalized by Bollerslev [1986],
named GARCH. The ε terms in the ARMA mean equation (2) are the innovations of the
t
timeseriesprocess. Engle[1982]definedthemasanautoregressiveconditionalheteroscedastic
process where all ε are of the form
t
ε = z σ , (3)
t t t
where z is an iid process with zero mean and unit variance. Although ε is serially uncor-
t t
related by definition its conditional variance equals σ2 and, therefore, may change over time.
t
All the GARCH models we consider in the following differ only in their functional form for
the conditional variance.
The variance equation of the GARCH(p,q) model can be expressed as
ε = z σ ,
t t t
z ∼ D (0,1) ,
t ϑ
p q
X X
σ2 = ω + α ε2 + β σ2 ,
t i t−i j t−j
i=1 j=1
= ω + α(B)ε2 + β(B)σ2 , (4)
t−1 t−1
where D (0,1) is the probability density function of the innovations or residuals with zero
ϑ
mean and unit variance. Optionally, ϑ are additional distributional parameters to describe
the skew and the shape of the distribution. If all the coefficients β are zero the GARCH
model is reduced to the ARCH model. Like for ARMA models a GARCH specification often
leads to a more parsimonious representation of the temporal dependencies and thus provides
a similar added flexibility over the linear ARCH model when parameterizing the conditional
variance. Bolerslev [1986] has shown that the GARCH(p,q) process is wide-sense stationary
with E(ε ) = 0, var(ε ) = ω/(1−α(1)−β(1)) and cov(ε ,ε ) = 0 for t 6= s if and only if
t t t s
α(1)+β(1) < 1.
The variance equation of the APARCH(p,q) model can be written as
ε = z σ ,
t t t
z ∼ D (0,1) ,
t ϑ
p q
X X
σδ = ω + α (|ε |−γ ε )δ + β σδ , (5)
t i t−i i t−i j t−j
i=1 j=1
where δ > 0 and −1 < γ < 1. This model has been introduced by Ding, Granger, and Engle
i
[1993]. It adds the flexibility of a varying exponent with an asymmetry coefficient to take the
4 An R and SPlus Software Implementation
leverage effect into account. A stationary solution exists if ω > 0, and Σ α κ +Σ β < 1,
i i i j j
where κ = E(|z|+γ z)δ. Note, that if γ 6= 0 and/or δ 6= 2 the κ depend on the assumptions
i i i
made on the innovation process. The family of APARCH models includes the ARCH and
GARCH models, and five other ARCH extensions as special cases:
• ARCH Model of Engle when δ = 2, γ = 0, and β = 0.
i j
• GARCH Model of Bollerslev when δ = 2, and γ = 0.
i
• TS-GARCH Model of Taylor and Schwert when δ = 1, and γ = 0.
i
• GJR-GARCH Model of Glosten, Jagannathan, and Runkle when δ = 2.
• T-ARCH Model of Zakoian when δ = 1.
• N-ARCH Model of Higgens and Bera when γ = 0, and β = 0.
i j
• Log-ARCH Model of Geweke and Pentula when δ → 0.
3. The Standard GARCH(1,1) Model
Code Snippet 1 shows how to write a basic S function named garch11Fit() to estimate
the parameters for Bollerslev’s GARCH(1,1) model. As the benchmark data set we use the
daily DEMGBP foreign exchange rates as supplied by Bollerslev and Ghysels [1996] and the
resultsobtainedbyFiorentini, CalzolarandPanattoni[1996]basedontheoptimizationofthe
log-likelihood function using analytic expressions for the gradient and Hessian matrix. This
setting is well accepted as the benchmark for GARCH(1,1) models. The function estimates
the parameters µ,ω,α,β in a sequence of several major steps: (1) the initialization of the
time series, (2) the initialization of the model parameters, (3) the settings for the conditional
distribution function, (4) the composition of the log-likelihood function, (5) the estimation
of the model parameters, and 6) the summary of the optimization results.
Given the model for the conditional mean and variance and an observed univariate return
series, we can use the maximum log-likelihood estimation approach to fit the parameters for
the specified model of the return series. The procedure infers the process innovations or
residuals by inverse filtering. Note, that this filtering transforms the observed process ε into
t
an uncorrelated white noise process z . The log-likelihood function then uses the inferred
t
innovations z to infer the corresponding conditional variances σ 2 via recursive substitution
t t
into the model-dependent conditional variance equations. Finally, the procedure uses the
inferred innovations and conditional variances to evaluate the appropriate log-likelihood ob-
jective function. The MLE concept interprets the density as a function of the parameter set,
conditionalonasetofsampleoutcomes. TheNormaldistributionisthestandarddistribution
when estimating and forecasting GARCH models. Using ε = z σ the log-likelihood function
t t t
of the Normal distribution is given by
LN(θ) = lnY p 1 e−2εσ2tt2 = lnY p 1 e−z2t2 (6)
(2πσ2) (2πσ2)
t t t t
1 X
= − [log(2π)+log(σ2)+z2] ,
2 t t
t
Journal of Statistical Software 5
or in general
Y
L (θ) = ln D (x ,E(x |Ω ),σ ) , (7)
N ϑ t t t−1 t
t
where D is the conditional distribution function. The second argument of D denotes
ϑ ϑ
the mean, and the third argument the standard deviation. The full set of parameters
θ includes the parameters from the mean equation (µ,a ,b ), from the variance equa-
1:m 1:n
tion (ω,α ,γ ,β ,δ), and the distributional parameters (ϑ) in the case of a non-normal
1:p 1:p 1:q
distribution function. For Bollerslev’s GARCH(1,1) model the parameter set reduces to
θ = {µ,ω,α ,β }. In the following we will suppress the index on the parameters α and
1 1
β if we consider the GARCH(1,1) model.
The parameters θ which fit the model best, are obtained by minimizing the“negative”log-
likelihood function. Some of the values of the parameter set θ are constrained to a finite or
semi-finite range. Note, that ω > 0 has to be constrained on positive values, and that α and
β have to be constrained in a finite interval [0,1). This requires a solver for constrained
numerical optimization problems. R offers the solvers nlminb() and optim(method="L-
BFGS-B") for constrained optimization, the first is also part of SPlus. These are R interfaces
to underlying Fortran routines from the PORT Mathematical Subroutine Library, Lucent
Technologies [1997], and to TOMS Algorithm 778, ACM [1997], respectively. In addition
we have implemented a sequential quadratic programming algorithm "sqp", Lucsan [1976],
which is more efficient compared to the other two solvers. And additionally, everything
is implemented in Fortran, the solver, the objective function, gradient vector, and Hessian
matrix.
Theoptimizerrequireavectorofinitialparametersforthemeanµ, aswellasfortheGARCH
coefficients ω, α, and β. The initial value for the mean is estimated from the mean µ of the
time series observations x. For the GARCH(1,1) model we initialize α = 0.1 and β = 0.8
by typical values of financial time series, and ω by the variance of the series adjusted by the
persistence ω = Var(x)∗(1−α−β). Further arguments to the GARCH fitting function are
the upper and lower box bounds, and optional control parameters.
3.1. How to fit Bollerslev’s GARCH(1,1) Model
In what follows we explicitly demonstrate how the parameter estimation for the GARCH(1,1)
model with normal innovations can be implemented in S. The argument list of the fitting
function garch11Fit(x) requires only the time series, the rest will be done automatically
step by step:
• Step 1 - Series Initialization: We save the time series x globally to have the values
available in later function calls without parsing the series through the argument list.
• Step 2 - Parameter Initialization: In the second step we initialize the set of model
parameters{θ},params,andthecorrespondingupperandlowerbounds. Intheexample
we use bounds lowerBounds and upperBounds which are wide enough to serve almost
every economic and financial GARCH(1,1) model, and define model parameters which
take typical values.
• Step 3 - Conditional Distribution: For the conditional distribution we use the Normal
distribution dnorm().
6 An R and SPlus Software Implementation
garch11Fit = function(x)
{
# Step 1: Initialize Time Series Globally:
x <<- x
# Step 2: Initialize Model Parameters and Bounds:
Mean = mean(x); Var = var(x); S = 1e-6
params = c(mu = Mean, omega = 0.1*Var, alpha = 0.1, beta = 0.8)
lowerBounds = c(mu = -10*abs(Mean), omega = S^2, alpha = S, beta = S)
upperBounds = c(mu = 10*abs(Mean), omega = 100*Var, alpha = 1-S, beta = 1-S)
# Step 3: Set Conditional Distribution Function:
garchDist = function(z, hh) { dnorm(x = z/hh)/hh }
# Step 4: Compose log-Likelihood Function:
garchLLH = function(parm) {
mu = parm[1]; omega = parm[2]; alpha = parm[3]; beta = parm[4]
z = (x-mu); Mean = mean(z^2)
# Use Filter Representation:
e = omega + alpha * c(Mean, z[-length(x)]^2)
h = filter(e, beta, "r", init = Mean)
hh = sqrt(abs(h))
llh = -sum(log(garchDist(z, hh)))
llh }
print(garchLLH(params))
# Step 5: Estimate Parameters and Compute Numerically Hessian:
fit = nlminb(start = params, objective = garchLLH,
lower = lowerBounds, upper = upperBounds, control = list(trace=3))
epsilon = 0.0001 * fit$par
Hessian = matrix(0, ncol = 4, nrow = 4)
for (i in 1:4) {
for (j in 1:4) {
x1 = x2 = x3 = x4 = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (garchLLH(x1)-garchLLH(x2)-garchLLH(x3)+garchLLH(x4))/
(4*epsilon[i]*epsilon[j])
}
}
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
}
(cid:4) Code Snippet 1: Example script which shows the six major steps to estimate the parameters
of a GARCH(1,1) time series model. The same steps are implemented in the“full version”of
garchFit() which allows the parameter estimation of general ARMA(m,n)-APARCH(p,q) with
several types of conditional distribution functions.
Journal of Statistical Software 7
• Step4-Log-LikelihoodFunction: ThefunctiongarchLLH()computesthe“negative”log-
likelihood for the GARCH(1,1) model. We use a fast and efficient filter representation
for the variance equation, such that no execution time limiting for loop will become
necessary.
• Step 5: - Parameter Estimation: For the GARCH(1,1) model optimization of the log-
likelihood function we use the constrained solver nlminb() which is available in R
and SPlus. To compute standard errors and t-values we evaluate the Hessian matrix
numerically.
• Step 6: - Summary Report: The results for the estimated parameters together with
standard errors and t-values are summarized and printed.
3.2. Case Study: The DEMGBP Benchmark
Through the complexity of the GARCH models it is evident that different software imple-
mentations have different functionalities, drawbacks and features and may lead to differences
in the numerical results. McCullough and Renfro [1999], Brooks, Burke and Persand [2001],
and Laurent and Peters [2003] compared the results of several software packages. These in-
vestigations demonstrate that there could be major differences between estimated GARCH
parameters from different software packages. Here we use for comparisons the well accepted
benchmark suggested by Fiorentini, Calzolari, and Panattoni [1996] to test our implementa-
tion.
ForthetimeseriesdatawetaketheDEMGBPdailyexchangeratesaspublishedbyBollerslev
and Ghysels [1996]. The series contains a total of 1975 daily observations sampled during the
period from January 2, 1984, to December 31, 1991. This benchmark was also used by
McCullough and Renfro [1999], and by Brooks, Burke, and Persand [2001] in their analysis
of several GARCH software packages. Figure 1 shows some stylized facts of the log-returns
of the daily DEM/GBP FX rates. The distribution is leptokurtic and skewed to negative
values. The log-returns have a mean value of µ = −0.00016, i.e. almost zero, a skewness of
ς = −0.25, and an excess kurtosis of κ = 3.62. The histogram of the density and the quantile-
quantile plot graphically display this behavior. Furthermore, the autocorrelation function of
the absolute values of the returns, which measures volatility, decay slowly.
We use the GARCH(1,1) model to fit the parameters of the time series using the function
garch11Fit() and compare these results with the benchmark computations of Fiorentini,
Calzolari, and Panattoni and with results obtained from other statistical software packages.
# Code Snippet 2: GARCH(1,1) Benchmark for the DEM/GBP Exchange Rates
> data(dem2gbp)
> garch11Fit(x = dem2gbp[, 1])
Coefficient(s):
Estimate Std. Error t value Pr(>|t|)
mu -0.00619040 0.00846211 -0.73154 0.46444724
omega 0.01076140 0.00285270 3.77236 0.00016171 ***
alpha 0.15313411 0.02652273 5.77369 7.7552e-09 ***
beta 0.80597365 0.03355251 24.02126 < 2.22e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
8 An R and SPlus Software Implementation
DEM/GBP FX Rate Histogram of demgbp
0.03 120
0.02 100
80
log Return 0.000.01 Density 4060
−0.01 20
−0.02 0
0 500 1000 1500 2000 −0.02 −0.01 0.00 0.01 0.02
index demgbp
Normal Q−Q Plot Series abs(demgbp)
0.03 1.0
0.02 0.8
mpirical Quantiles 0.000.01 ACF 0.40.6
E −0.01 0.2
−0.02 0.0
−3 −2 −1 0 1 2 3 0 5 10 15 20 25 30
Rmetrics FCP FinmetricsG@ARCH Shazam TSP
Normal Quantiles Lag
μ -0.006190 -0.006190 -0.006194 -0.006183 -0.006194 -0.006190
(cid:4) Fωigure0.10:107R6etu0r.0n1s0,76their distrib0.u01ti0o7n6, t0h.0e10q7u6antile-qua0.n0t1i0le76plo0t.0,10a7n6d the autocorrelation function of
α 0.1531 0.1531 0.1540 0.1531 0.1531 0.1531
the volatility for the DEMGBP benchmark data set.
β 0.8060 0.8060 0.8048 0.8059 0.8060 0.8060
T h e e s t i m a tRemdetmricsodel p a r a m e t e FrCsP aBreencthomamrkore tha n f o u r Gd@igAeRsCtHs Oexxactly the same numbers as
Estimate StdErrors EstimateStdErrors EstimateStdErrorsdErrors
theFCPbenchmark. Note,weusednumericallycomputedderivativeswhereasthebenchmark
useμd a0.n00a6l1y9t0i4c0a.0l0ly846c2a11lculated-0.0d0e6r19iv04a1t.i0v0e84s6.21T2 he obs-0e.0r0v6a1t84ion0.0t08h4a62t it is not rea lly nece ssary to
comωpu0t.e010t7h6e14d.0e0r2i8v5a27t0ives ana0ly.0t1i0c7a61l4ly.00w28a5s271already0m.01a0d76e0b0y.00D28o51ornik [2000] a nd Lau rent and
α 0.1531340.0265227 0.1531340.0265228 0.153407 0.026569
Peters [2001].
β 0.8059740.0335525 0.8059740.0335527 0.805879 0.033542
Rmetrics: FCP Benchmark: Ox/Garch: Splus/Finmetrics:
Estimate StdError t Value Estimate StdErrort Value Estimate StdError t Value Estimate StdErrort Value
μ -0.00619040.0084621 -0.732 -0.00619040.0084621 -0.732 -0.006184 0.008462 -0.731 -0.006053 0.00847 -0.715
ω 0.0107610.0028527 3.77 0.0107610.0028527 3.77 0.0107610.0028506 3.77 0.0108960.0029103 3.74
α 0.15313 0.026523 5.77 0.15313 0.026523 5.77 0.15341 0.026569 5.77 0.15421 0.026830 5.75
β 0.80597 0.033553 24.0 0.80597 0.033553 24.0 0.80588 0.033542 24.0 0.80445 0.034037 23.6
(cid:4) Table 1: Comparison of the results of the GARCH(1,1) parameter estimation for the DEMGBP
benchmarkdatasetobtainedfromseveralsoftwarepackages. Thefirsttwocolumnsshowtheresultsfromthe
GARCH(1,1)exampleprogram,thesecondgroupfromtheFCPbenchmarkresults,thethirdgroupfromOx’
G@RCH,andthelastgroupfromSPlus’Finmetrics. Thelefthandnumbersaretheparameterestimatesand
the right hand number the standard errors computed from the Hessian matrix.
Journal of Statistical Software 9
4. Alternative Conditional Distributions
Different types of conditional distribution functions D are discussed in literature. These are
the Normal distribution which we used in the previous section, the standardized Student-t
distribution and the generalized error distribution and their skewed versions. Since only the
symmetric Normal distribution of these functions is implemented in R’s base environment
and also in most of the other statistical software packages, we discuss the implementation of
the other distributions in some more detail.
The default choice for the distribution D of the innovations z of a GARCH process is the
t
Standardized Normal Probability Function
f?(z) = √1 e−z22 . (8)
2π
The probability function or density is named standardized, marked by a star ?, because f?(z)
has zero mean and unit variance. This can easily be verified computing the moments
Z ∞
µ? = zrf?(z)dz . (9)
r
−∞
Note, that µ? ≡ 1 and σ? ≡ 1 are the normalization conditions, that µ? defines the mean
0 1 1
µ ≡ 0, and µ? the variance σ2 ≡ 1.
2
An arbitrary Normal distribution located around a mean value µ and scaled by the standard
deviation σ can be obtained by introducing a location and a scale parameter through the
transformation
f(z)dz → 1f?(cid:16)z−µ(cid:17)dz = √1 e−(z2−σµ2)2dz . (10)
σ σ σ 2π
The central moments µ of f(z) can simply be expressed in terms of the moments µ? of the
r r
standardized distribution f?(z). Odd central moments are zero and those of even order can
be computed from
Z ∞ 2r (cid:16) 1(cid:17)
µ = (z−µ)2rf(z)dz = σ2rµ? = σ2r √ Γ r+ , (11)
2r 2r π 2
−∞
yielding µ = 0, and µ = 3. The degree of asymmetry γ of a probability function, named
2 4 1
skewness, and the degree of peakedness γ , named excess kurtosis, can be measured by nor-
2
malized forms of the third and fourth central moments
µ µ
3 4
γ = = 0 , γ = −3 = 0 . (12)
1 µ32/2 2 µ22
However, if we like to model an asymmetric and/or leptokurtic shape of the innovations
we have to draw or to model z from a standardized probability function which depends
t
on additional shape parameters which modify the skewness and kurtosis. However, it is
important that the probability has still zero mean and unit variance. Otherwise, it would be
impossible, oratleastdifficult, toseparatethefluctuationsinthemeanandvariancefromthe
10 An R and SPlus Software Implementation
Student−t: Density Student−t Distribution
0
1.
8
0.
0.8 nu=10
nu=2.5
6
0.
6
0.
f(z) 0.4 nu=5 F(z) 4 nu=5
0.
0.2 nu=10 0.2 nu=2.5
0 0
0. 0.
−4 −2 0 2 4 −4 −2 0 2 4
z z
(cid:4) Figure 2: The densi1ty000a0n Rdanddisotmri bDuetviioantesfor the Standardized StudKeunrtt-otsisdistribution with unit vari-
ance for three different degrees of freedom, ν = 2.5,5,10. The graph for the largest value of ν is almost
0
undistinguishable0.4from a Normal distribution. 1
8
3
fluctuations in t0.he shape of the density. In a first step we consider still symmetric probability
fcuonncstidioenrstbhueDensitytgwen0.2itehraalinzeadddeirtrioornadlissthraibpuetipoanraamndetKurtosistehrew6Shtiuchdemnto-dtedlsistthriebkuutirotnoswis.ithAsunexitamvaprlieasncwee,
4
both relevant in modelling GARCH processes.
1
0. 2
4.1. Student-t Distribution
0.0 0
Bollerslev [1987], Hsieh [1989)], Baillie and Bollerslev [1989], Bollerslev, Chou, and Kroner
−4 −2 0 2 4 6 0 5 10 15 20
[1992], Palm [1996], Pagan [1996)], and Palm and Vlaar [1997] among others showed that
r nu
the Student-t distribution better captures the observed kurtosis in empirical log-return time
series. The density f?(z|ν) of the Standardized Student-t Distribution can be expressed as
Γ(ν+1) 1
f?(z|ν) = 2 (13)
pπ(ν −2)Γ(ν) (cid:16) (cid:17)ν+1
2 1+ z2 2
ν−2
1 1
= √ ,
ν −2B(cid:0)1, ν(cid:1)(cid:16) (cid:17)ν+1
2 2 1+ z2 2
ν−2
where ν > 2 is the shape parameter and B(a,b) =Γ(a)Γ(b)/Γ(a +b) the Beta function.
Note, when setting µ = 0 and σ2 = ν/(ν−2) formula (13) results in the usual one-parameter
expression for the Student-t distribution as implemented in the S function dt.
Again, arbitrary location and scale parameters µ and σ can be introduced via the transfor-
mation z → z−µ. Odd central moments of the standardized Student-t distribution are zero
σ