---------------------------------------------------
# #
# Chapter 8: Maximum Likelihood Estimation #
# of the GARCH(1,1)-Model #
# #
---------------------------------------------------
Das GARCH(1,1)-Modell ist definiert durch
S(t_k) = S(t_{k-1}) * [ 1 + vol(t_{k-1}) * phi_k ] (1)
mit einer Vol-Spezifikation
vol^2(t_k) = (2)
alpha0 + alpha1 * ret^2(t_k) + beta1 * vol^2(t_{k-1})
Mit Hilfe des Lemmas von week9.txt koennen wir zeigen, dass
diese rekursive oder implizite Definition der Volatilitaet,
die Vol taucht auf der linken und der rechten Seite von (2)
auf, aequivalent ist zu der expliziten Darstellung
vol^2(t_k) = beta1^k * vol^2(t_0) + alpha0 * (1-beta1^k)/(1-beta1)
+ alpha1*[ ret^2(t_k) + beta1 * ret^2(t_{k-1}) + ...
... + beta1^{k-1} * ret^2(t_1) ] (3)
-> kurze Rechnung an der Tafel.
Da beta1 < 1 vorauszusetzen ist, koennen wir die beta1^k in
der ersten Zeile von (3) vernachlaessigen und bekommen
vol^2(t_k) = alpha0/(1-beta1) + alpha1*[ ret^2(t_k) +
beta1 * ret^2(t_{k-1}) + ... + beta1^{k-1} * ret^2(t_1) ] (4)
Um den Parametern eine intuitivere Bedeutung zu geben, repara-
metrisieren wir gemaess
beta1 = w
alpha0 = (1-w) * w0 * bsvol^2 (5)
alpha1 = (1-w) * (1-w0)
so dass (4) in
vol^2(t_k) = w0*bsvol^2 + (1-w0) * exponentially_weighted_vol (5)
= w0*bsvol^2 + (1-w0) * (1-w) * [ ret^2(t_k) +
w * ret^2(t_{k-1}) + ... + w^{k-1} * ret^2(t_1) ]
uebergeht. Die uebliche GARCH(1,1)-Bedingung
alpha1 + beta1 < 1 (6)
liest sich in den neuen Parametern
1 - w0 * (1-w) < 1 (7)
In week9.txt hatten wir uns klar gemacht, dass das w im wesent-
lichen den Zeithorizont fuer die exponentially weighted vol fest-
legt: Parametrisieren wir das w gemaess
w = 1 - 1/d <=> d = 1/(1-w) (8)
dann ist das d im wesentlichen vergleichbar mit dem d aus dem
ARCH(d)-Modell. Das GARCH(1,1)-Modell ist also im wesentlichen
ein ARCH(d)-Modell mit exponentially weighted vol. Mit (8) liest
sich die GARCH(1,1)-Bedingung (6):
1 - w0/d < 1 (9)
was offensichtlich automatisch erfuellt ist wegen w0 in (0,1)
und d>=1.
Schauen wir uns zunaechst einige GARCH(1,1)-Pfade an:
#
# Start R-Session:
#
pathGarch11 = function( N , bsvol , w0 , d )
{
# we put start vol = bsvol
S = rep(0,N)
ret = rep(0,N)
phi = rnorm(N)
S0 = 100
vol = rep(0,N) # GARCH-vol
ewvolsq = rep(0,N) # exp-weighted-vol^2
w = 1 - 1/d
# the first step, k=1, is special:
S[1] = S0 * ( 1 + bsvol * phi[1] )
ret[1] = (S[1]-S0)/S0
ewvolsq[1] = w * bsvol^2 + (1-w)*ret[1]^2
vol[1] = sqrt( w0*bsvol^2 + (1-w0)*ewvolsq[1] )
for(k in 2:N)
{
S[k] = S[k-1] * ( 1 + vol[k-1] * phi[k] )
ret[k] = (S[k]-S[k-1])/S[k-1]
# use the Lemma to calc exp-weighted-vol:
ewvolsq[k] = w * ewvolsq[k-1] + (1-w)*ret[k]^2
# GARCH(1,1) vol-specification:
vol[k] = sqrt( w0*bsvol^2 + (1-w0)*ewvolsq[k] )
}
par(mfrow=c(1,2))
plot(S,type="l",ylim=c(0,250))
plot(ret,type="l",ylim=c(-0.08,0.08))
result = list(S,ret)
names(result) = c("S","ret")
return(result)
}
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=14 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=2 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=1 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=50 )
# Wir wollen jetzt die Parameter fuer ein GARCH(1,1)-Modell
# mit Hilfe der Maximum-Likelihood-Methode berechnen. Die
# log-Likelihood-Funktion ist nach wie vor gegeben durch
#
# log(tilde L) = (10)
#
# - sum_{k=1}^N { log[vol(t_{k-1})] + 1/2 * y_k^2 / vol^2(t_{k-1}) }
#
# wobei die y_k = ret(t_k) durch die realisierten returns gege-
# ben sind. Da das Berechnen der exponentially weighted vol ein
# bischen laenger dauert, tun wir fuer alle
#
# dd in { 1 , 2 , ... , d }
#
# diese vol im voraus berechnen, so dass sie, wenn wir beim
# Maximieren der Likelihood-Funktion etwa das w0 von 0 nach 1
# variieren und das w oder das d festhalten, sie nicht jedesmal
# neu berechnet werden muss:
GenerateEWVol = function( d )
{
n = length(ret)
bsvol = sd(ret)
result = matrix( 0 , nrow=n , ncol=d )
ewvolsq = rep(0,n)
for(dd in 1:d)
{
w = 1 - 1/dd
ewvolsq = rep(0,n)
# k=1 is special
ewvolsq[1] = w*bsvol^2 + (1-w)*ret[1]^2
for(k in 2:n)
{
ewvolsq[k] = w*ewvolsq[k-1] + (1-w)*ret[k]^2
}
result[,dd] = sqrt(ewvolsq)
}
return(result)
}
# first we have to generate the return vector ret:
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=14 )
ret = res$ret
EWVol = GenerateEWVol( 40 )
plot(EWVol[,14],type="l")
plot(EWVol[,2],type="l")
plot(EWVol[,40],type="l")
# jetzt berechnen wir die log-Likelihood-Funktion,
# wir koennen die EWVol[,d]'s von oben benutzen:
logL = function( bsvol , w0 , d )
{
w1 = 1 - w0
n = length(ret)
# exponentially weighted ARCH(d)-Vol:
vol = sqrt( w0*bsvol^2 + w1*EWVol[,d]^2 )
# shift vol-vector by 1 element ( there is a vol(t_{k-1}) in
# the likelihood-function, not a vol(t_k) )
vol = vol[-n] # remove last one
vol = c(bsvol,vol) # add a new first one
# the actual calculation:
terms = log(vol) + 1/2*ret^2/vol^2
result = -sum(terms)
return(result)
}
# before we calibrate to DAX-data, we check on an artificial
# GARCH(1,1) return set where we know the actual model para-
# meters (bsvol,w0,d):
N = 2500
bsvol = 0.015
w0 = 0.3
d = 15
res = pathGarch11( N , bsvol , w0 , d )
ret = res$ret
summary(ret)
# before logL-calculation: precalc ew-vol!
# the result has to have the name EWVol:
EWVol = GenerateEWVol( 40 )
# let's check:
logL( bsvol=0.015, w0=0.3, d=15 ) # alright, some number
# let's look for the maximum:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol
# we make a contour-plot in the (w0,d)-plane:
w0 = seq(from=0.05,to=0.95,by=0.05)
dd = 1:40
nw = length(w0)
nd = length(dd)
# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}# still very fast evaluation of logL since vol is precalculated!
# let's look at the result:
contour( dd , w0 , z )
contour( dd , w0 , z , nlevels = 50 )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
contour( dd , w0 , z , nlevels = 100 , zlim=c(10200,10400) )
dax = read.table("C:/Users/detlef/HSRM/Vorlesungen/WS201516/StatLearning/week10/DAX.txt",header=TRUE,sep=";")
#dax = read.table("C:/Users/Admin/HSRM/Teaching/WS201516/StatLearning/week9/DAX.txt",header=TRUE,sep=";")
dax = na.omit(dax)
S = dax$index
n=length(S)
ret = rep(0,n)
for(i in 2:n)
{
ret[i] = (S[i]-S[i-1])/S[i-1]
}
plot(ret,type="l")
# before analysis: precalc ew-vol!
EWVol = GenerateEWVol( 40 )
# let's look for the maximum:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol
# we make a contour-plot in the (w0,d)-plane:
w0 = seq(from=0.05,to=0.95,by=0.05)
dd = 2:40
nw = length(w0)
nd = length(dd)
# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}# still very fast evaluation of logL since vol is precalculated!
# let's look at the result:
contour( dd , w0 , z )
contour( dd , w0 , z , nlevels = 50 )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
image( dd , w0 , z )
persp( dd , w0 , z )
require(plot3D)
persp3D( dd , w0 , z )
w0 = seq(from=0.01,to=0.3,by=0.01)
dd = 5:20
nw = length(w0)
nd = length(dd)
# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z) # 9883.602
logL( bsvol , 0.125 , 10 ) # 9883.629
image( dd , w0 , z )
persp( dd , w0 , z )
persp3D( dd , w0 , z )
# let's check for bsvol: we fix d=10 and vary bsvol and w0:
w0 = seq(from=0.01,to=0.3,by=0.01)
bsvol = seq(from=0.5/100,to=2.5/100,by=0.05/100)
nw = length(w0)
nbs = length(bsvol)
# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nbs , ncol=nw )
for(i in 1:nbs)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol[i] , w0[j] , d=10 )
}
}
contour( bsvol , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z) # 9883.756
logL( sd(ret) , 0.125 , 10 ) # 9883.629 -> bsvol=sd(ret) is ok!
# let's recall the ARCH(d) logL-value:
logL_Archd = function( bsvol , w0 , d )
{
w1 = 1 - w0
n = length(ret)
cumsum_retsquared = cumsum(ret^2)
# remove last d elements:
cumsum_retsquared_shifted = cumsum_retsquared[-((n-d+1):n)]
# fill in d zeros at the first d places:
cumsum_retsquared_shifted = c(rep(0,d),cumsum_retsquared_shifted)
sum_dretsquared = cumsum_retsquared - cumsum_retsquared_shifted
# now we have to divide this sum by d, but the first d elements
# are special:
divisor1 = 1:d
divisor2 = rep(d,n-d)
divisor = c(divisor1,divisor2)
sum_dretsquared = sum_dretsquared/divisor
# now we can calculate ARCH(d)-vol:
vol = sqrt( w0*bsvol^2 + w1*sum_dretsquared )
# shift vol-vector by 1 element ( there is a vol(t_{k-1}) in
# the likelihood-function, not a vol(t_k) )
vol = vol[-n] # remove last one
vol = c(sqrt(w0*bsvol^2),vol) # add a new first one
# the actual calculation:
terms = log(vol) + 1/2*ret^2/vol^2
result = -sum(terms)
return(result)
}
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL_Archd( bsvol , w0[j] , dd[i] )
}
}
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z) # 9873.058
logL_Archd( bsvol , 0.15 , 14 ) # 9872.649
# Also:
#
# GARCH(1,1) = ARCH(d) mit exponentially weighted vol:
# Max[ logL_GARCH(1,1) ] = 9883.629
# at d=10, w0=0.125
#
# besser als ARCH(d) mit gleichgewichteter vol:
# Max[ logL_ARCH(d) ] = 9873.058
# at d=14, w0=0.15