#---------------------------#
#   Loesungen UeBlatt 12    #
#---------------------------#


---------------
#  Aufgabe 1  #
---------------

# Wir modifizieren die ARCH(1) log-Likelihood-Funktion 
# von week9.txt: dort hatten wir folgenden code:
logL = function( bsvol , w0 )
{
  w1 = 1 - w0
  vol = sqrt(w0*bsvol^2+w1*ret^2)     
                                      
  # shift vol by 1 element:
  n = length(ret)
  vol = vol[-n]                       # remove last one
  vol = c(sqrt(w0*bsvol^2),vol)       # add a new first one
  terms = log(vol) + 1/2*ret^2/vol^2
  result = -sum(terms)
  return(result)
}

# wir machen eine Modifikation gemaess der 
# neuen Modell-Spezifikation:


logLmod = function( bsvol , w0 , alpha )
{
  w1 = 1 - w0
  # der ifelse-Befehl ist sehr nuetzlich:

  wretsq = ifelse( ret>=0 , max(w1-alpha,0)*ret^2 , max(w1+alpha,0)*ret^2 )  

  vol = sqrt( w0*bsvol^2 + wretsq )     

  # shift vol by 1 element:
  n = length(ret)
  vol = vol[-n]                       # remove last one
  vol = c(sqrt(w0*bsvol^2),vol)       # add a new first one
  terms = log(vol) + 1/2*ret^2/vol^2
  result = -sum(terms)
  return(result)
}



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",main="DAX-returns")



alpha = seq( from=-0.98 , to=0.98 , by=0.02 )
w0 = seq( from=0.02 , to=0.98 , by=0.02 )
bsvol = sd(ret)
nalph = length(alpha)
nw = length(w0)
z = matrix( 0 , nrow=nalph , ncol=nw )
for(i in 1:nalph)
{
  for(j in 1:nw )
  {
    z[i,j] = logLmod( bsvol , w0[j] , alpha[i] )         
  }
}

# let's look at the result:
info = "logL_mod"
contour( alpha , w0 , z , nlevels = 50 , main=info )
contour( alpha , w0 , z , nlevels = 50 , zlim=c(9000,10000) , main=info )
contour( alpha , w0 , z , nlevels = 100 , zlim=c(9400,9600) , main=info )

# more refined
alpha = seq( from=0.01 , to=0.4 , by=0.01 )
w0 = seq( from=0.6 , to=0.99 , by=0.01 )
nalph = length(alpha)
nw = length(w0)
z = matrix( 0 , nrow=nalph , ncol=nw )
for(i in 1:nalph)
{
  for(j in 1:nw )
  {
    z[i,j] = logLmod( bsvol , w0[j] , alpha[i] )         
  }
}
contour( alpha , w0 , z , nlevels = 100 , zlim=c(9400,9600) , main=info )
max(z)                                       # 9548.758
logLmod( bsvol , w0=0.752 , alpha=0.17 )     # 9548.16


# check bsvol:
bsvol = seq( from=1/100 , to=2/100 , by=0.01/100 )
w0 = seq( from=0.5 , to=0.83 , by=0.01 )
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] = logLmod( bsvol[i] , w0[j] , alpha=0.17 )           
  }
}

contour( bsvol , w0 , z , nlevels = 100 , zlim=c(9000,10000) , main=info )  
contour( bsvol , w0 , z , nlevels = 100 , zlim=c(9400,9600) , main=info )  

max(z)                                          # 9549.034      
logLmod( sd(ret) , 0.7 , alpha=0.17 )           # 9548.108
sd(ret) # 0.01388
logLmod(  0.0141 , 0.7 , alpha=0.17 )           # 9548.922

# bsvol slightly larger..






# direktes Bestimmen des Maximums, 
# ohne die Bilder anzuschauen:

w0 = seq( from=0.5 , to=0.85 , by=0.01 )
alpha = seq( from=0.01 , to=0.4 , by=0.01 )
bsvol = seq(from=0.1/100,to=4/100,by=0.05/100)
nw = length(w0)
na = length(alpha)
nbs = length(bsvol)


# matrix z = logL for contour-plot:

z = array( 0 , dim=c(nbs,nw,na) )
for(i in 1:nbs)
{
  for(j in 1:nw )
  {
    for(k in 1:na)
    {
      z[i,j,k] = logLmod( bsvol[i] , w0[j] , alpha[k] )
    }
  }
}# ok, dauert jetzt ein bischen..

max(z)

# wo wird das Maximum angenommen?

which(z==max(z), arr.ind=TRUE )

ijk = which(z==max(z), arr.ind=TRUE )

bsvolmax = bsvol[ ijk[1] ]
w0max    =    w0[ ijk[2] ]
alphamax = alpha[ ijk[3] ]

bsvolmax          # 0.014
w0max             # 0.72
alphamax          # 0.17


logLmod( sd(ret) , 0.7 , alpha=0.17 )
logLmod( bsvolmax, w0max, alphamax )       # ok, noch ein kleines bischen mehr..