#---------------------------------#
#     Week8: Das ARCH Modell      #         
#---------------------------------#


In week7.txt hatten wir das Modell

S(t_k) = S(t_{k-1}) * [ 1 + vol(t_{k-1})*phi_k ]                    (1)

mit 

vol^2(t_{k-1}) = stddev_d^2(t_{k-1})

               = 1/d * [ret^2(t_{k-1})+...+ret^2(t_{k-d})]          (2)

simuliert und festgestellt, dass die Volatilitaet 
vol(t_k) dann relativ bald nach 0 geht, die simulierten 
Pfade blieben konstant nach einer gewissen Weile und 
haben ihren Wert nicht mehr geaendert. Offensichtlich 
ist  vol(t_{k-1}) = 0  und  S(t_k) = const  eine 
Loesung der Gleichungen (1,2) fuer beliebige Zufalls-
zahlen phi_k. Auf der anderen Seite haben wir gesehen, 
dass das Modell 

S(t_k) = S(t_{k-1}) * [ 1 + bsvol * phi_k ]                         (3)

mit einer konstanten Black-Scholes Volatilitaet bsvol 
dieses Problem nicht hat, die simulierten Pfade sahen 
ganz vernunftig aus. Allerdings hat unsere Datenanalyse 
in week6.txt deutlich gezeigt, dass das Volatilitaets-
niveau eben nicht konstant ist, und ein gutes Modell 
sollte das beruecksichtigen. Deshalb machen wir jetzt 
den Ansatz

vol^2(t_{k-1}) = w0*bsvol^2 + w1*stddev_d(t_{k-1})^2                (4)

               = w0 * bsvol^2 + 
                 w1 * 1/d * [ret^2(t_{k-1})+...+ret^2(t_{k-d})]     

mit positiven Gewichten w0+w1=1. Das Preis-Model (1) 
mit der Volatilitaetsspezifikation (4) heisst das ARCH- 
Modell oder genauer das ARCH(d)-Modell, weil eben die 
Historie der letzten d Tage mit beruecksichtigt wird. 



Bevor wir ein paar Pfade simulieren, noch 2 Bemerkungen:


Bemerkung 1) Die Buchstaben A, R, C und H stehen fuer 

AR = AutoRegressive, das meint, dass zukuenftige returns 
                     von vergangenen returns abhaengen, und:
CH = Conditional Heteroskedastic, das meint, dass man 
     zeitabhaengige instantane Volatilitaeten hat. 

Etwas genauer: Man unterscheidet zwischen bedingten oder  
conditional oder instantanen Varianzen oder Volatilitaeten 

vol^2(t_{k-1}) = E[ ret^2(t_k) | t_{k-1} ]        # bis Zeit t_{k-1} ist alles bekannt,
                                                  # Erwartungswert nur ueber phi_k
und unbedingten Varianzen oder Volatilitaeten

volbar^2 = E[ ret^2(t_k) ]                        # Erwartungswert ueber alle phi_k,...,phi_1
 
Heteroskedastic meint einfach: nicht konstant, sondern 
zeitabhaengig.

Bemerkung 2) In der Vol-Spezifikation (4) koennte man 
natuerlich die Kombination der Konstanten w0*bsvol^2 
als eine neue Konstante alpha0 und etwa w1/d als ein 
alpha1 definieren oder allgemeiner einen Ansatz 

vol^2(t_{k-1}) = alpha0 + 
  alpha1 * ret^2(t_{k-1}) + ... + alphad * ret^2(t_{k-d})           (5)

machen mit d+1 Parametern alpha0,...,alphad. Typischer-
weise findet man die Definition (5) in den Buechern als 
die Definition des ARCH(d)-Modells. Wenn man das Modell 
allerdings an konkreten Zeitreihendaten fitten will, wird 
man nicht wirklich d+1 unabhaengige Parameter zulassen 
wollen, sondern eher eine Parametrisierung (4) waehlen, 
moeglicherweise nicht mit einer Gleich-Gewichtung der 
letzten d Returns, sondern vielleicht mit einer linear 
oder exponentiell abfallenden Gewichtung, das passiert 
dann in den sogenannten GARCH-Modellen, das "G" steht 
dann da einfach nur fuer "Generalized". Weiterhin haben 
die Parameter in der Spezifikation (4) eine intuitivere 
Bedeutung als die Parameter in (5). 


Schauen wir uns zunaechst einfach ein paar Pfade an. 
Dazu modifizieren wir unsere pathNaivArch-Funktion 
aus dem week7.txt:



###    Start R-Session    ### 
 

pathArchd = function( N , bsvol , w0 , d , vol0 )
{
  S   = rep(0,N)
  ret = rep(0,N)
  phi = rnorm(N)
  S0  = 100
  vol = rep(vol0,N)       # ARCH-vol

  sumretd = 0

  # the first step, k=1, is special:
  
  S[1]    = S0*( 1 + vol0 * phi[1] )
  ret[1]  = (S[1]-S0)/S0
  vol[1]  = sqrt( w0*bsvol^2 + (1-w0)*ret[1]^2 )
  sumretd = sumretd + ret[1]^2

  for(k in 2:N)
  {
    if(k<=d)
    {
      S[k]    = S[k-1]*( 1 + vol[k-1] * phi[k] )
      ret[k]  = (S[k]-S[k-1])/S[k-1]
      sumretd = sumretd + ret[k]^2
      vol[k]  = sqrt( w0*bsvol^2 + (1-w0)*sumretd/k )
    }
    else
    {
      S[k]    = S[k-1]*( 1+ vol[k-1] * phi[k] )
      ret[k]  = (S[k]-S[k-1])/S[k-1]
      sumretd = sumretd + ret[k]^2
      sumretd = sumretd - ret[k-d]^2
      vol[k]  = sqrt( w0*bsvol^2 + (1-w0)*abs(sumretd)/d )
    }
  }#next 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)
}



# let's check:
# this should be Black-Scholes:

res = pathArchd(2500, bsvol=0.01, w0=1, d=15, 0.01)
res = pathArchd(2500, bsvol=0.01, w0=1, d=1,  0.01)


# this is with variable vol:

res = pathArchd(2500, bsvol=0.01, w0=0.5,  d=15, 0.01)    # vol looks still quite uniform
res = pathArchd(2500, bsvol=0.01, w0=0.1,  d=15, 0.01)    # this shows vol clustering
res = pathArchd(2500, bsvol=0.01, w0=0.1,  d=30, 0.01)
res = pathArchd(2500, bsvol=0.01, w0=0.1,  d=1,  0.01)
res = pathArchd(2500, bsvol=0.01, w0=0.01, d=1,  0.01)
res = pathArchd(2500, bsvol=0.01, w0=0,    d=1,  0.01)    # w0>0 stabilizes the model
res = pathArchd(2500, bsvol=0.01, w0=0,    d=15, 0.01)    # w0>0 stabilizes the model
res = pathArchd(2500, bsvol=0.01, w0=0.001,d=15, 0.01)    # w0>0 stabilizes the model



Das sieht soweit ok aus und wir koennten jetzt anfangen 
zu ueberlegen, wie man dieses Modell an eine konkrete 
Zeitreihe fitten kann. Also wir wollen die Modellparameter 
w0, bsvol und d so bestimmen, dass sie "moeglichst gut", 
das ist dann mathematisch genauer zu praezisieren, zu den 
Zeitreihendaten passen. Das machen wir naechste Woche 
mit Hilfe der Maximum-Likelihood-Methode.


Hier wollen wir uns jetzt noch ein bischen mit dem 
contour()-Befehl vertraut machen, mit dem man die 
Hoehenlinien einer Funktion plotten kann. Den Befehl 
werden wir dann naechste Woche benutzen, um ein gewisses 
Gefuehl fuer die Likelihood-Funktion zu bekommen, die 
wir also naechste Woche berechnen wollen. 




# Beispiel zu Contour-Plot:
#
# Wir definieren eine Funktion f = f(x,y) von 2 Variablen:

f = function(x,y)
{
  result = ( -x^2 + y^2 ) * exp( -0.25*x^2 - 0.5*y^2 )
  return(result)
}

# der contour()-Befehl hat die Syntax 
# contour(  vector ,  vector , matrix ) = 
# contour( x-Achse , y-Achse , f(x,y) )


# x- und y-Achse:

x = seq(from=-4,to=4,by=0.1)
y = x


# die f(x,y)-Werte muessen wir in eine 
# Matrix schreiben, das koennen wir mit 
# der outer()-Funktion machen:

z = outer( x , y , FUN = f )         # damit dieser Befehl funktioniert, muss die Funktion
                                     # f Vektoren als input akzeptieren koennen; wenn sie 
                                     # das nicht tut, kann man das z natuerlich auch einfach
                                     # mit einem Doppel-Loop anlegen, ist weiter unten gemacht

# Und jetzt die Hoehenlinien:

contour( x , y , z )
contour( x , y , z , nlevels=50 )    # bessere Aufloesung


# anstatt mit dem outer()-Befehl kann man das z 
# auch mit einer doppelten Schleife anlegen:

nx = length(x)
ny = length(y)
z = matrix( 0 , nrow=nx , ncol=ny )

for(i in 1:nx)
{
  for(j in 1:ny)
  {
    z[i,j] = f(x[i],y[j])            # jetzt Skalare als Argument, keine Vektoren
  }
}


# Und jetzt wieder die Hoehenlinien:

contour( x , y , z )
contour( x , y , z , nlevels=50 )    # more refined


# also useful:

image( x , y , z )
persp( x , y , z )


# with plot3D package:
require(plot3D)

persp3D( x , y , z )