---------------------------------------------------
# #
# Chapter 7: Exponentially Weighted #
# Volatilities #
# and Introduction to GARCH #
# #
---------------------------------------------------
Fuer einen Zeithorizont d waren die d-Tages Volatili-
taeten definiert durch
vol^2(t_{k-1}) = [ret^2(t_{k-1})+...+ret^2(t_{k-d})] / d (1)
und wir hatten bei der ARCH(d)-Kalibrierung gesehen, dass
sowohl fuer die SPX-Daten als auch fuer die DAX-Daten
ein 3-Wochen Zeithorizont mit d=14 am besten passt.
Nun ist es plausibel anzunehmen, dass more recent
returns einen groesseren Impact haben als returns, die
schon laenger in der Vergangenheit liegen. Das heisst,
es waere plausibel, anstatt der Gleichgewichtung der
d returns in (1) mit Gewicht 1/d etwa folgende Gewich-
tung der returns vorzunehmen mit einem w < 1:
vol^2(t_k) = 1 / normalization * (2)
[ ret^2(t_k) + w * ret^2(t_{k-1}) + w^2 * ret^2(t_{k-2}) + ...
... + w^{d-1} * ret^2(t_{k-d+1}) ]
mit der Normierungskonstanten
normalization = 1 + w + w^2 + ... + w^{d-1}
= ( 1 - w^d ) / ( 1 - w )
Da fuer w<1 die Potenzen w^d recht schnell nach 0 konver-
gieren, koennen wir die Information, wieviele der returns
wir denn als relevant ansehen, auch vollstaendig durch
das w codieren und in Gleichung (2) einfach d=k setzen:
vol^2(t_k) = 1 / normalization * (3)
[ ret^2(t_k) + w * ret^2(t_{k-1}) + w^2 * ret^2(t_{k-2}) + ...
... + w^{k-1} * ret^2(t_1) ]
mit der Normierungskonstanten
normalization_k = ( 1 - w^k ) / ( 1 - w )
-> 1 / ( 1 - w ) for k to infinity (4)
Wir wuerden also folgende Volatilitaet mit `exponentially
weighted' returns betrachten:
vol^2(t_k) = (5)
(1-w)*[ ret^2(t_k) + w*ret^2(t_{k-1}) + ... + w^{k-1}*ret^2(t_1) ]
Das w parametrisieren wir durch ein d folgendermassen (siehe
Gleichung (9) unten):
Waehlen wir w deutlich kleiner als 1, etwa w=1/2, dann werden
die Potenzen w^i sehr schnell sehr klein und wir wuerden
effektiv nur 4 oder 5 returns zur Berechnung der Vol mit be-
ruecksichtigen. Waehlen wir jedoch w nahe bei 1, tun wir
offensichtlich eine lange Historie von returns mit beruecksich-
tigen.
Anstatt des etwas unintuitiven w's (w=0.92 oder w=0.98 hat
jetzt nicht unmittelbar eine anschauliche Bedeutung), moechten
wir jetzt das w durch einen anschaulicheren Parameter d aus-
druecken, der im wesentlichen folgende Bedeutung haben soll:
Bei der Berechnung der exponentially weighted vol (5) tun wir
`effektiv' d returns aufsummieren.
Vergleichen wir also die beiden unterschiedlichen Gewichtungen:
Alle returns mit gleichem Gewicht,
1/d 1/d 1/d . . . 1/d (d Stueck)
(6)
versus exponentially weighted returns:
(1-w) (1-w)w (1-w)w^2 . . . (1-w)w^{k-2} (1-w)^{k-1} (k Stueck)
(7)
Wir moechten das w als Funktion von d so waehlen, das man die
Volatilitaeten (5) und (1) `am ehesten' miteinander vergleichen
kann. Wenn wir die ersten d der exponentiellen Gewichte (7) auf-
summieren, erhalten wir
(1-w) * [ 1 + w + ... w^(d-1) ] = 1 - w^d (8)
Waehlen wir jetzt etwa
w = 1 - kappa / d (9)
bekommen wir
Summe in (8) = 1 - w^d
= 1 - (1-kappa/d)^d
~ 1 - e^(-kappa)
und das Gewicht fuer den d-ten return ist von der Groessen-
ordnung
w^d = (1-kappa/d)^d ~ e^(-kappa) (10)
Das wesentliche ist hier, dass das von d unabhaengig ist und
nicht etwa fuer grosse d ~ 100 nach 0 geht. Fuer kappa waeren
verschiedene Werte denkbar: Wenn man etwa fordert, dass der
jeweils letzte return ret^2(t_k) in (1) und in (5) mit dem-
selben Gewicht 1/d gewichtet werden soll, haben wir offensicht-
lich 1/d = 1 - w oder
w = 1 - 1/d (11)
Neben dieser Parametrisierung von w schauen wir uns ebenfalls
die Wahl kappa = log(4) an, also
w = 1 - ln(4)/d (12)
Bei dieser Wahl ist w^d ~ e^{-kappa} = 1/4 und die Summe der
ersten d Gewichte ist 1 - e^{-kappa} = 3/4. Bei der Wahl (11)
ist w^d ~ e^{-kappa} = 1/e und die Summe der ersten d Gewichte
ist 1 - e^{-kappa} = 1 - 1/e ~ 0.632 .
Plotten wir einmal die Gewichte (6) und (7) mit den Parametri-
sierungen (11) und (12):
# Start R:
#
k = 300
d = 60
w = 1 - 1/d # Parametrisierung (11)
equalweights = c( rep(1/d,d) , rep(0,k-d) )
expweights = (1-w)*w^(0:(k-1))
plot(expweights,type="l",col=2)
lines(equalweights)
k = 300
d = 60
w = 1 - log(4)/d # Parametrisierung (12)
equalweights = c( rep(1/d,d) , rep(0,k-d) )
expweights = (1-w)*w^(0:(k-1))
plot(expweights,type="l",col=2)
lines(equalweights)
# redo above for other d's: d=15
#
# End R
Wir wollen jetzt fuer die DAX-Zeitreihe die realisierten Vola-
tilitaeten gemaess Formel (1) und Formel (5) mit den Parametri-
sierungen (11) und (12) berechnen. Bevor wir das machen, beweisen
wir folgendes Lemma, welches uns erlaubt, die Volatilitaeten (5)
auf einfache Weise induktiv zu berechnen:
Lemma: Fuer die Groessen v_k (die Quadrate der Volatilitaeten)
und r_k (die Quadrate der returns) gelte folgende rekursive
Beziehung (alpha = 1-w):
v_k = w * v_{k-1} + alpha * r_k (13)
fuer k = 1,2,3,... . Dann gilt folgende explizite Formel fuer
die v_k:
v_k = w^k * v_0 + alpha*[ r_k + w*r_{k-1} + ... + w^{k-1}*r_1 ] (14)
Beweis: Mit Induktion, kurze Rechnung an der Tafel.
#-----------------------------------------------------------------
# Start R-Session: Equally weighted vs. Exponentially weighted Vol
#
dax = read.table("C:/Users/detlef/HSRM/Vorlesungen/WS201516/StatLearning/week8/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")
EquallyWeightedVol = function( d , ret )
{
n = length(ret)
vol = rep(0,n)
# we use code from last week:
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
# equally weighted vol:
vol = sqrt( sum_dretsquared )
return(vol)
}
ExpWeightedVol1 = function( d , ret )
{
n = length(ret)
vol = rep(0,n)
volsq = rep(0,n) # vol^2
vol0sq = 0
w = 1 - 1/d
# w = 1 - log(4)/d
volsq[1] = w * vol0sq + (1-w)*ret[1]^2
for(k in 2:n)
{
volsq[k] = w * volsq[k-1] + (1-w)*ret[k]^2
}
vol = sqrt( volsq )
return(vol)
}
ExpWeightedVol2 = function( d , ret )
{
n = length(ret)
vol = rep(0,n)
volsq = rep(0,n) # vol^2
vol0sq = 0
# w = 1 - 1/d
w = 1 - log(4)/d
volsq[1] = w * vol0sq + (1-w)*ret[1]^2
for(k in 2:n)
{
volsq[k] = w * volsq[k-1] + (1-w)*ret[k]^2
}
vol = sqrt( volsq )
return(vol)
}
d = 15
vol0 = EquallyWeightedVol( d , ret )
vol1 = ExpWeightedVol1( d , ret )
vol2 = ExpWeightedVol2( d , ret )
plot(vol0,type="l")
lines(vol1,col="red")
lines(vol2,col="green")
d = 60
vol0 = EquallyWeightedVol( d , ret )
vol1 = ExpWeightedVol1( d , ret )
vol2 = ExpWeightedVol2( d , ret )
plot(vol0,type="l")
lines(vol1,col="red")
lines(vol2,col="green")
# Finally, we use an artificially generated return set
# of Black-Scholes returns with constant bsvol = 1%
# and we put in by hand an up-return of 8% at k = 300
# and a down-return of 12% at k = 600. We generate
# N = 1000 returns corresponding to a time horizon
# of 4 years:
ret = rnorm( 1000 , mean=0 , sd=1/100 )
ret[300] = 8/100
ret[600] = -12/100
plot(ret,type="l")
d = 100
vol0 = EquallyWeightedVol( d , ret )
vol1 = ExpWeightedVol1( d , ret )
vol2 = ExpWeightedVol2( d , ret )
plot(vol0,type="l")
lines(vol1,col="red")
lines(vol2,col="green")
#
# End R-Session: Equally weighted vs. Exponentially weighted Vol
#-----------------------------------------------------------------
Die Verwendung von exponentially weighted returns scheint
also ziemlich vernuenftig zu sein: ein einziger grosser
return erzeugt nicht mehr so einen rechteckigen Block der
Laenge d in der Vol-Zeitreihe, sondern der Vol-Level tut
sofort wieder sinken nach Verstreichen dieses returns,
wobei die Abkling-Geschwindigkeit durch den Parameter w =
1 - kappa/d festgelegt ist (mit kappa=1 oder kappa=log4).
Nachdem wir aus der ARCH(d)-Kalibrierung d=14 als optimalen
Wert bei gleichgewichteter Vol-Spezifikation gefunden haben,
ist es jetzt natuerlich sehr interessant zu sehen, welchen
Wert von w oder d wir bei exponentially weighted returns
bekommen, und ob die Likelihood-Funktion sich auf diese Weise
noch vergroessern laesst, das Modell also noch statistisch
signifikanter wird. Schauen wir uns das an:
Wir benutzen also folgende Vol-Spezifikation, die im wesent-
lichen aequivalent ist zur Vol-Spezifikation des GARCH(1,1)-
Modells, was wir weiter unten als Spezialfall des GARCH(p,q)-
Modells definieren:
vol^2(t_k) = w0*bsvol^2 + (1-w0) * exponentially_weighted_vol (15)
= w0*bsvol^2 + (1-w0) * (1-w) * [ ret^2(t_k) +
w * ret^2(t_{k-1}) + ... + w^{k-1} * ret^2(t_1) ]
Die Likelihood-Funktion ist nach wie vor gegeben durch
log(tilde L) = (16)
- 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. Anstelle des w's nehmen wir das intuitivere d aus
Gleichung (9),
w = 1 - kappa / d
<=> d = kappa / ( 1 - w ) (17)
als Input-Parameter fuer die log-Likelihood Funktion, mit
kappa = 1 oder kappa = log(4) wie oben diskutiert.
Bevor wir jetzt die Likelihood-Funktion maximieren, geben
wir noch die allgemeine Definition des GARCH(1,1)- und des
GARCH(p,q)-Modells an, wie man sie ueblicherweise in den
Buechern findet. Das `G' in GARCH steht dabei fuer `Genera-
lized' Arch.
Das GARCH(1,1)-Modell ist definiert durch die Vol-Spezifi-
kation
vol^2(t_k) = (18)
alpha0 + alpha1 * ret^2(t_k) + beta1 * vol^2(t_{k-1})
Wir koennen das Lemma von oben anwenden und bekommen (kurze
Rechnung an der Tafel):
vol^2(t_k) = beta1^k * vol^2(t_0) + alpha0 * (1 - beta1^k)
+ alpha1*(1-beta1)*[ ret^2(t_k) + beta1 * ret^2(t_{k-1}) + ...
... + beta1^{k-1} * ret^2(t_1) ]
Da beta1 < 1 vorauszusetzen ist, koennen wir die beta1^k in
der ersten Zeile 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) ] (19)
und das hat genau die Form (15) von oben mit
beta1 = w
alpha0 = (1-w) * w0 * bsvol^2 (20)
alpha1 = (1-w) * (1-w0)
Schliesslich ist das GARCH(p,q)-Modell definiert durch die
Vol-Spezifikation
vol^2(t_k) = alpha0 (21)
+ alpha1 * ret^2(t_k) + ... + alphap * ret^2(t_{k-p+1})
+ beta1 * vol^2(t_{k-1}) + ... + betaq * vol^2(t_{k-q})
aus der sich dann, wie ueblich, der Preis-Prozess S(t_k) be-
rechnet gemaess
S(t_k) = S(t_{k-1}) * [ 1 + vol(t_{k-1}) * phi_k ] (22)
mit einer standard-normalverteilten Zufallszahl phi_k.
Zum Abschluss dieses week9.txt schauen wir uns ein paar
GARCH(1,1)-Pfade an:
pathGarch11 = function( N , bsvol , w0 , d , kappa )
{
# we put start vol = bsvol
S = rep(0,N)
ret = rep(0,N)
phi = rnorm(N)
S0 = 100
vol = rep(vol0,N) # GARCH-vol
ewvolsq = rep(0,N) # exp-weighted-vol^2
w = 1 - kappa/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 , kappa=1 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=14 , kappa=log(4) )