---------------------------------------------------
#                                                 #
#     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) )




Option Pricing
Time Series Models















Hochschule RheinMain,
Applied Mathematics:


<<previous chapter       back to chapter overview       next chapter>>