#----------------------------------------------# # # # Week11: Motivation des GARCH-Modells # # # #----------------------------------------------# 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 2-3 Wochen Zeithorizont mit d=14 oder d=11 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 exponentiell gewichteten 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 # entspricht Laenge der Zeitreihe d = 60 w = 1 - 1/d # Parametrisierung (11) equalweights = c( rep(1/d,d) , rep(0,k-d) ) # k Stueck expweights = (1-w)*w^(0:(k-1)) # k Stueck 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, schreiben wir noch das folgende Lemma hin, welches uns erlaubt, die Vola- tilitaeten (5) auf einfache Weise induktiv zu berechnen, mit alpha = 1-w : Lemma: Fuer die Groessen v_k (die Quadrate der Volatilitaeten) und r_k (die Quadrate der returns) gelte die folgende rekursive Beziehung: 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, Uebungsblatt11. #----------------------------------------------------------------- # Start R-Session: Equally weighted vs. Exponentially weighted Vol # dax = read.table("C:/Users/detlef/OneDrive/hochschule/Vorlesungen/WS2223/Datenanalyse-mit-R/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 = 250 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") 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") plot(vol0[700:900],type="l") lines(vol1[700:900],col="red") lines(vol2[700:900],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 9% at k = 300 # and a down-return of 15% 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] = 9/100 ret[600] = -15/100 plot(ret,type="l") d = 30 vol0 = EquallyWeightedVol( d , ret ) vol1 = ExpWeightedVol1( d , ret ) vol2 = ExpWeightedVol2( d , ret ) plot(vol0, type="l", ylim=c(0,4/100) ) 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=11 oder 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}) Durch Iteration bekommt man die folgende explizite Formel, die man leicht durch Induktion beweisen kann (Uebungsblatt 11): 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) ] Da beta1 < 1 vorauszusetzen ist, koennen wir die beta1^k in der ersten Zeile vernachlaessigen und bekommen fuer groessere k 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 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) # soll GARCH-vol werden 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.01 , w0 = 0.15 , d=14 , kappa=1 ) res = pathGarch11( N=2500 , bsvol=0.01 , w0 = 0.15 , d=14 , kappa=log(4) )