#------------------------------------------------#
#                                                #
#    Week9: Maximum Likelihood Estimation des    #
#                ARCH(1)-Modells                 #
#                                                #
#------------------------------------------------#


Letzte Woche hatten wir das ARCH(d)-Modell definiert 
und uns einige Pfade angeschaut und gesehen, dass man 
durch Variation der Parameter bereits eine recht 
grosse Vielfalt an Preisdynamiken generieren kann. 

Die Frage ist nun, wie man festellen kann, welche 
dieser Preisdynamiken, d.h. welcher Set von Para-
metern, denn am ehesten zu einer gegebenen Zeitreihe 
passt. Das wollen wir hier zunaechst fuer das ARCH(1)-
Modell diskutieren. Die Logik fuer das ARCH(d)-Modell 
ist dann voellig analog, und die entsprechenden 
numerischen Rechnungen dazu schauen wir uns dann 
naechste Woche an (man hat dann 3 Parameter, nicht 
nur 2).


Das ARCH(1)-Model war gegeben durch 

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

mit der Vol-Spezifikation 

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

mit Gewichten w0+w1=1 und einer konstanten, zeit-
unabhaengigen Volatilitaet bsvol. Fuer w0=1, w1=0 
erhaelt man das zeitdiskrete Black-Scholes Modell 
mit der Volatilitaet bsvol (und Drift-Parameter = 0).

Wir moechten jetzt dieses Modell an die DAX-Zeitreihe 
fitten. Also wir moechten die Parameter w0 (und damit 
w1=1-w0) und bsvol finden, die 'am besten' zu den 
Daten der DAX-Zeitreihe passen. Dazu schreiben wir 
zunaechst (1) um zu 

ret(t_k) = vol(t_{k-1}) * phi_k                                   (3)

und erinnern uns daran, dass phi_k eine normalver-
teilte Zufallszahl ist. Wenn man also bei Zeit t_{k-1} 
ist, dann ist die W'keit, dass sich ein return im 
Intervall [y_k,y_k+dy_k] realisiert, gegeben durch 

Prob[ ret(t_k) in [y_k,y_k+dy_k] | t_{k-1} ]                      (4)

 = 1/sqrt{2pi} 1/vol(t_{k-1}) * exp{ -1/2*y_k^2/vol^2(t_{k-1}) } dy_k

Die W'keit, dass sich dann alle returns der DAX-Zeit-
reihe in Intervallen [y_k,y_k+dy_k] befinden, k=1,2,...,N, 
ist dann einfach gegeben durch das Produkt der W'keiten 
aus (4), also:

Prob[ DAX-ret in [y_k,y_k+dy_k] fuer alle k ] 

  = prod_{k=1}^N 1/sqrt{2pi} 1/vol(t_{k-1}) *                     (5)

    prod_{k=1}^N exp{ -1/2*y_k^2/vol^2(t_{k-1}) } dy_k

Setzen wir jetzt die Formel (2), 

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

auf der rechten Seite von (5) ein, koennen wir die 
W'keit aus (5),

Prob[ DAX-ret in [y_k,y_k+dy_k] fuer alle k ] 

auch ausschliesslich als Funktion von bsvol und w0 auf-
fassen, denn die y_k sind ja die bekannten returns, die 
sich in den letzten 10 Jahren realisiert haben. Um das 
deutlicher zu machen, dass man das jetzt als Funktion der 
Modell-Parameter auffasst, benutzt man den Buchstaben 
L fuer Likelihood und schreibt 

Prob[ DAX-ret in [y_k,y_k+dy_k] fuer alle k ] 

  =  L( w0 , bsvol )                                              (6)

  =  L( w0 , bsvol | y_k=ret(t_k) )                                   

wobei die letzte Zeile mit dem senkrechten Strich, "gegeben 
die y_k=ret(t_k)", nochmal deutlich machen soll, dass die 
y_k gegebene Daten sind, naemlich die Returns, die sich 
realisiert haben. Dieses L heisst dann die Likelihood-Funktion. 

Ein Standard-Schaetz-Verfahren der Statistik besteht dann 
ganz einfach darin, diejenigen Modell-Parameter w0 und bsvol 
zu nehmen, die die Likelihood-Funktion L maximieren. 


Da wir fuer eine 10-Jahres Zeitreihe etwa N=2500 returns 
haben, ist L gegeben durch ein Produkt von 2500 Faktoren, 
L wird sich also verhalten wie

  L ~ exp( +N*irgendwas )  oder   L ~ exp( -N*irgendwas ), 

das ist entweder sehr, sehr gross oder sehr, sehr klein. 
Um numerische Probleme zu vermeiden, betrachtet man deshalb 
typischerweise den Logarithmus der Likelihood Funktion und 
benutzt die Tatsache, dass das Maximieren von L aequivalent 
ist zu dem Maximieren von log(L), da der Logarithmus ja 
streng monoton steigend ist. Also schaut man sich an 

log(L) = sum_{k=1}^N { log[1/vol(t_{k-1})] - 

                       1/2 * y_k^2 / vol^2(t_{k-1}) } 
         + const

wobei die Konstante gegeben ist durch 

const = N*log[1/sqrt{2pi}] + sum_{k=1}^N log[dy_k] 

und unabhaengig ist von den Modell-Parametern w0 und bsvol, 
also kann man sie beim Maximieren weglassen. Wir muessen 
also das Maximum von 


log(tilde L) = log(tilde L)( bsvol , w0 )                         (7)

= - sum_{k=1}^N { log[vol(t_{k-1})] + 1/2 * y_k^2 / vol^2(t_{k-1}) }        

mit y_k = ret(t_k) und

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


bestimmen. Wenn man das analytisch probiert, mit Ableitung 
gleich Null setzen, stellt man nach einer Weile fest, dass 
das nicht so ohne weiteres geht, es gibt keine analytische 
closed form solution. Das ist nun leider eine ziemlich ty-
pische Situation, das Maximum einer Likelihood-Funktion muss 
typischerweise mit einem numerischen Algorithmus berechnet 
werden.



###   Start R-Session:   ### 

# before we fit to DAX-data, we check on an artificially 
# generated ARCH(1) return set: let's use the pathArchd 
# function from last week with d=1 to generate these 
# artificial return data:

# from week8.txt:
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)
}
# END from week8.txt

N     = 2500
bsvol = 0.01
w0    = 0.5
res = pathArchd(N,bsvol,w0,d=1,vol0=bsvol)
ret = res$ret
summary(ret)

# fuer diese returns ret wollen wir jetzt also die log-
# Likelihodd Funktion berechnen und maximieren, und es 
# sollte dann also rauskommen, dass das Maximum bei 
# bsvol = 0.01 = 1% und w0 = 0.5 liegt (oder zumindest 
# in der Naehe davon; es sind ja Zufallszahlen und durch
# Zufall koennte es ja auch mal sein, dass etwa 1.4% 
# und w0 = 0.3 noch besser passen):


# we calculate the log-Likelihood function:

logL = function( bsvol , w0 )
{
  w1 = 1 - w0
  vol = sqrt(w0*bsvol^2+w1*ret^2)     # man kann den Vektor ret benutzen, obwohl 
                                      # der nicht als Argument uebergeben wurde
  # shift vol by 1 element:
  n = length(ret)
  vol = vol[-n]                       # remove last one
  vol = c(bsvol,vol)                  # add a new first one
  terms = log(vol) + 1/2*ret^2/vol^2
  result = -sum(terms)
  return(result)
}


# let's check:
logL( bsvol=0.01 , w0=0.5 )           # ok, we get some number


# bevor man einfach blind einen Optimierungs-Algorithmus 
# anwirft, ist es immer gut, ein gewisses Gefuehl fuer die 
# zu maximierende Funktion zu bekommen, sofern das denn 
# moeglich ist. Hier haben wir nur 2 Parameter und wir koen-
# nen uns die Hoehenlinien der Funktion anschauen. Machen 
# wir jetzt also einen Contour-Plot der log-likelihood 
# Funktion:

# bsvol von 0.1% bis 4%, Schrittweite 0.1%:
# (recall that the model is unreasonable for w0=0 or bsvol=0)

bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0    = seq( from=0.02    , to=0.98  , by=0.02 )


# the following does not work, since logL does not accept 
# vectors as input; the reason is that in logL there is a line
# vol = sqrt( w0*bsvol^2 + w1*ret^2 ) 
# which does not make sense if all w0, w1 and ret are vectors, 
# in particular, of different lengths:

z = outer( bsvol , w0 , FUN = logL )          # does not work


# let's do this with a loop:

nbs = length(bsvol)
nw  = length(w0)
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] )         # this is ok
  }
}

# let's look at the result:
contour( bsvol , w0 , z )
contour( bsvol , w0 , z , nlevels = 60 )
contour( bsvol , w0 , z , nlevels = 60 , zlim=c(10000,11000) )       # eventually adjust zlim-values

# okay, that works!!  -> we recover correct parameter values!!


# some more pictures (would require more fine-tuning):
image( bsvol , w0 , z )
persp( bsvol , w0 , z )

require(plot3D)
persp3D( bsvol , w0 , z )




##  now lets fit to DAX-data:  ##

dax = read.table("C:/Users/detlef/OneDrive/hochschule/Vorlesungen/WS2223/Datenanalyse-mit-R/DAX.txt",header=TRUE,sep=";")
dax = na.omit(dax)
summary(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")


# calculate DAX-logL and make setup for contour-plot:

bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0    = seq( from=0.02    , to=0.98  , by=0.02 )
nbs = length(bsvol)
nw  = length(w0)
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] )         
  }
}

# let's look at the result:
contour( bsvol , w0 , z )
contour( bsvol , w0 , z , nlevels = 60 )
contour( bsvol , w0 , z , nlevels = 60 , zlim=c(8000,10000) )
contour( bsvol , w0 , z , nlevels = 60 , zlim=c(9000,9600) )


# also: w0 etwa 0.7, bsvol schauen wir noch etwas genauer:

w0 = 0.7
bsvol = seq( from=0.01 , to=0.02 , by=0.001 )
logLL = rep(0,length(bsvol))
for(i in 1:length(bsvol))
{
  logLL[i] = logL(bsvol[i],w0)
  print(c(bsvol[i],logLL[i]))
}

plot(bsvol,logLL,type="l")

# also: bsvol = 1.4%
# das ist auch tatsaechlich identisch mit der Standardabweichung 
# der DAX-returns ueber den gesamten 10-Jahres-Zeitraum:

sd(ret)

# sieht also alles recht vernuenftig aus.