#-------------------------------------------------#
#                                                 #
#    Week12: Maximum Likelihood Estimation des    #
#              GARCH(1,1)-Modells                 #
#                                                 #
#-------------------------------------------------#




Das GARCH(1,1)-Modell ist definiert durch 

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


mit einer Vol-Spezifikation

  vol^2(t_k) =                                                       (2)

      alpha0 + alpha1 * ret^2(t_k) + beta1 * vol^2(t_{k-1})      


Mit Hilfe des Teil (b) der Aufgabe 2 vom UeBlatt11 koennen wir zeigen, 
dass diese rekursive oder implizite Definition der Volatilitaet,  
die vol taucht auf der linken und der rechten Seite von (2) auf, 
aequivalent ist zu der expliziten Darstellung

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) ]        (3)


Da beta1 < 1 vorauszusetzen ist, koennen wir die beta1^k in 
der ersten Zeile von (3) 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) ]      (4)

Um den Parametern eine intuitivere Bedeutung zu geben, repara-
metrisieren wir gemaess
 
  beta1  = w
  alpha0 = (1-w) * w0 * bsvol^2                                      (5)
  alpha1 = (1-w) * (1-w0) 
  
so dass (4) in 

vol^2(t_k) = w0*bsvol^2 + (1-w0) * exponentially_weighted_vol        (5)

           = w0*bsvol^2 + (1-w0) * (1-w) * [ ret^2(t_k) + 

                w * ret^2(t_{k-1}) + ... + w^{k-1} * ret^2(t_1) ]

uebergeht. Die uebliche GARCH(1,1)-Bedingung 

  alpha1 + beta1 < 1                                                 (6)

liest sich in den neuen Parametern

  1 - w0 * (1-w) < 1                                                 (7)

Im week11.txt hatten wir uns klar gemacht, dass das w im wesent-
lichen den Zeithorizont fuer die exponentially weighted vol fest-
legt: Parametrisieren wir das w gemaess 

      w = 1 - 1/d   <=>   d = 1/(1-w)                                (8)

dann ist das d im wesentlichen vergleichbar mit dem d aus dem 
ARCH(d)-Modell. Das GARCH(1,1)-Modell ist also im wesentlichen 
ein ARCH(d)-Modell mit exponentially weighted vol. Mit (8) liest 
sich die GARCH(1,1)-Bedingung (6):

  1 - w0/d < 1                                                       (9)

was offensichtlich automatisch erfuellt ist wegen w0 in (0,1) 
und d>=1. 

Schauen wir uns zunaechst einige GARCH(1,1)-Pfade an:

    

###  Start R-Session:  ### 

pathGarch11 = function( N , bsvol , w0 , d )  
{
  # we put start vol = bsvol
  S = rep(0,N)
  ret = rep(0,N)
  phi = rnorm(N)
  S0 = 100
  vol = rep(0,N)          # GARCH-vol
  ewvolsq = rep(0,N)      # exp-weighted-vol^2
  w = 1 - 1/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 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=2 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=1 )
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=50 )


# Wir wollen jetzt die Parameter fuer ein GARCH(1,1)-Modell 
# mit Hilfe der Maximum-Likelihood-Methode berechnen. Die 
# log-Likelihood-Funktion ist nach wie vor gegeben durch 
#
# log(tilde L) =                                                      (10)
#
# - 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. Da das Berechnen der exponentially weighted vol ein 
# bischen laenger dauert, tun wir fuer alle 
#
#    dd in { 1 , 2 , ... , d }
# 
# diese vol im voraus berechnen, so dass sie, wenn wir beim 
# Maximieren der Likelihood-Funktion etwa das w0 von 0 nach 1 
# variieren und das w oder das d festhalten, sie nicht jedesmal 
# neu berechnet werden muss:


GenerateEWVol = function( d )
{
  n = length(ret)
  bsvol = sd(ret)
  result = matrix( 0 , nrow=n , ncol=d )
  ewvolsq = rep(0,n)

  for(dd in 1:d)
  {
    w = 1 - 1/dd
    ewvolsq = rep(0,n)
    
    # k=1 is special
    ewvolsq[1] = w*bsvol^2 + (1-w)*ret[1]^2
    for(k in 2:n)
    {
      ewvolsq[k] = w*ewvolsq[k-1] + (1-w)*ret[k]^2
    }

    result[,dd] = sqrt(ewvolsq)
  }
  return(result)
}

# first we have to generate the return vector ret:
res = pathGarch11( N=2500 , bsvol=0.014 , w0 = 0.15 , d=14 )
ret = res$ret

EWVol = GenerateEWVol( 40 )
plot(EWVol[,14], type="l", ylim=c(0,0.05) )
plot(EWVol[,2] , type="l", ylim=c(0,0.05) )
plot(EWVol[,40], type="l", ylim=c(0,0.05) )


# jetzt berechnen wir die log-Likelihood-Funktion, 
# wir koennen die EWVol[,d]'s von oben benutzen:

logL = function( bsvol , w0 , d )
{
  w1 = 1 - w0
  n = length(ret)

  # exponentially weighted ARCH(d)-Vol:
  vol = sqrt( w0*bsvol^2 + w1*EWVol[,d]^2 )

  # shift vol-vector by 1 element ( there is a vol(t_{k-1}) in 
  # the likelihood-function, not a vol(t_k) )
  vol = vol[-n]                           # remove last one
  vol = c(bsvol,vol)                      # add a new first one

  # the actual calculation:
  terms = log(vol) + 1/2*ret^2/vol^2
  result = -sum(terms)

  return(result)
}


# before we calibrate to DAX-data, we check on an artificial 
# GARCH(1,1) return set where we know the actual model para-
# meters (bsvol,w0,d):

N = 2500
bsvol = 0.015
w0 = 0.3
d = 15
res = pathGarch11( N , bsvol , w0 , d )
ret = res$ret
summary(ret)
# before logL-calculation: precalc ew-vol! 
# the result has to have the name EWVol:
EWVol = GenerateEWVol( 40 )

# let's check:
logL( bsvol=0.015, w0=0.3, d=15 )          # alright, some number

# let's look for the maximum:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol

# we make a contour-plot in the (w0,d)-plane:

w0 = seq(from=0.02,to=0.98,by=0.02)
dd = 1:40
nw = length(w0)
nd = length(dd)

# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
  for(j in 1:nw )
  {
    z[i,j] = logL( bsvol , w0[j] , dd[i] )
  }
}# still very fast evaluation of logL since vol is precalculated!

 
# let's look at the result:
contour( dd , w0 , z )
contour( dd , w0 , z , nlevels = 50 )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )



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")
# before analysis: precalc ew-vol!
EWVol = GenerateEWVol( 40 )

# let's look for the maximum:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol

# we make a contour-plot in the (w0,d)-plane:

w0 = seq(from=0.02,to=0.98,by=0.02)
dd = 2:40
nw = length(w0)
nd = length(dd)

# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
  for(j in 1:nw )
  {
    z[i,j] = logL( bsvol , w0[j] , dd[i] )
  }
}# still very fast evaluation of logL since vol is precalculated!

 
# let's look at the result:
contour( dd , w0 , z )
contour( dd , w0 , z , nlevels = 50 )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
image( dd , w0 , z )
persp( dd , w0 , z )
require(plot3D)
persp3D( dd , w0 , z )


w0 = seq(from=0.01,to=0.3,by=0.01)
dd = 5:20
nw = length(w0)
nd = length(dd)

# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
  for(j in 1:nw )
  {
    z[i,j] = logL( bsvol , w0[j] , dd[i] )
  }
}

contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z)                        # 9883.602
logL( bsvol , 0.125 , 10 )    # 9883.629
image( dd , w0 , z )
persp( dd , w0 , z )
persp3D( dd , w0 , z )



# let's check for bsvol: we fix d=10 and vary bsvol and w0:
w0 = seq(from=0.01,to=0.3,by=0.01)
bsvol = seq(from=0.5/100,to=2.5/100,by=0.05/100)

nw = length(w0)
nbs = length(bsvol)

# matrix z = logL for contour-plot:
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] , d=10 )
  }
}
contour( bsvol , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z)                          # 9883.756
logL( sd(ret) , 0.125 , 10 )    # 9883.629 -> bsvol=sd(ret) is ok!



# let's recall the ARCH(d) logL-value:

logL_Archd = function( bsvol , w0 , d )
{
  w1 = 1 - w0
  n = length(ret)

  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
  
  # now we can calculate ARCH(d)-vol:
  vol = sqrt( w0*bsvol^2 + w1*sum_dretsquared )

  # shift vol-vector by 1 element ( there is a vol(t_{k-1}) in 
  # the likelihood-function, not a vol(t_k) )
  vol = vol[-n]                           # remove last one
  vol = c(sqrt(w0*bsvol^2),vol)           # add a new first one

  # the actual calculation:
  terms = log(vol) + 1/2*ret^2/vol^2
  result = -sum(terms)

  return(result)
}

bsvol = sd(ret)
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
  for(j in 1:nw )
  {
    z[i,j] = logL_Archd( bsvol , w0[j] , dd[i] )
  }
}

contour( dd , w0 , z , nlevels = 150 , zlim=c(9000,10000) )
contour( dd , w0 , z , nlevels = 150 , zlim=c(9700,10000) )
max(z)                             # 9873.058
logL_Archd( bsvol , 0.18 , 11 )    # 9873.058


# Also: 
#
# GARCH(1,1) = ARCH(d) mit exponentially weighted vol:
# Max[ logL_GARCH(1,1) ] = 9883.629
# at d=10, w0=0.125
# 
# besser als ARCH(d) mit gleichgewichteter vol:
# Max[ logL_ARCH(d) ] = 9873.058
# at d=11, w0=0.18