#--------------------------#
#                          #  
#   Loesungen Ueblatt 6    #
#                          #
#--------------------------#



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

dax = read.table("C:/Users/detlef/OneDrive/hochschule/Vorlesungen/WS2223/Datenanalyse-mit-R/DAX5y.txt",header=TRUE,sep=";")

head(dax)                     # sieht ganz gut aus..
tail(dax)
summary(dax)                  # mmh..
str(dax)                      # bei mir: sind alles Faktoren.. (kann abhaengen von der R-Version die Sie benutzen)

names(dax)                    # funktioniert

# extrahieren wir die Zahlen aus der "Adj.Close"-Spalte: 

S = dax$Adj.Close
S                             # das sind 1264 "Levels"   (kann abhaengen von der R-Version)
plot(S,type="l")              # ok, offensichtlich keine Zahlen.. 

S1 = as.numeric(S)            # Befehl wird ausgefuehrt,       
S1                            # aber die Zahlen sind unsinnig..

S2 = as.numeric(as.character(S))    # gibt eine warning message..
S2                                  # ..aber die Zahlen machen Sinn 
plot(S2,type="l")

n = length(S2)
n

# wir berechnen die returns:

ret = rep(0,n)       
for(i in 2:n)
{
  ret[i] = (S2[i]-S2[i-1])/S2[i-1]
}
plot(ret,type="l")




#---------------#
#   Aufgabe 2   #
#---------------#

# Wir benutzen die Funktionen
#
#  dDayMean()
#  dDayStdDev()
#  dDayNormRet()
#  
# aus dem week6.txt: 


dDayMean = function( d , ret )
{
  n = length(ret)
  result = rep(0,n)
  summe = 0

  for(i in 1:n)
  {
    summe = summe + ret[i]
    if(i > d) 
    {
      summe = summe - ret[i-d]
      result[i] = summe/d
    }
    else
    {
      result[i] = summe/i
    }
  }
  return(result)
}

dDayStdDev = function( d , ret )
{
  n = length(ret)
  result = rep(0,n)
  summe = 0

  for(i in 1:n)
  {
    summe = summe + ret[i]*ret[i]
    if(i > d) 
    {
      summe = summe - ret[i-d]*ret[i-d]
      result[i] = summe/d
    }
    else
    {
      result[i] = summe/i
    }
  }
  # a theoretical 0 could numerically become slightly negative:
  result = abs(result)     # we want to take a square root
  result = sqrt(result)    # square root element by element
  return(result)
}

dDayNormRet = function( d , ret )
{
  stddev = dDayStdDev(d,ret)

  # stddev could be 0 if data are constant:
  stddev = pmax(stddev,0.00000001)     # "parallel max": maximum element by element

  result = ret/stddev                  # division element by element
  return(result)
}



#-------------------#
#  Analysis S&P500  #
#-------------------#

# Jetzt koennen wir die Daten laden und analysieren: 

spx = read.table("C:/Users/detlef/OneDrive/hochschule/Vorlesungen/WS2223/Datenanalyse-mit-R/SPX.txt",header=TRUE,sep=";")
head(spx)
tail(spx)
summary(spx)
str(spx)                      # Adj.Close is numeric, no Factor

names(spx)
S = spx[,2]                   # ist aequivalent zu S = spx$Adj.Close
plot(S,type="l")
plot(log(S/S[1]),type="l")    # straight line would be constant exponential growth
n=length(S)
n

# wir berechnen die returns:

ret = rep(0,n)       
for(i in 2:n)
{
  ret[i] = (S[i]-S[i-1])/S[i-1]
}
plot(ret,type="l")


# let's take a look at Stylized Fact 1:

mean20 = dDayMean(20,ret)
mean60 = dDayMean(60,ret)
mean180 = dDayMean(180,ret)
plot(mean20)
points(mean60,col="red")
points(mean180,col="yellow")


# stddev's = daily volatilities:

stddev20 = dDayStdDev(20,ret)
stddev60 = dDayStdDev(60,ret)
stddev180 = dDayStdDev(180,ret)
stddev_n = dDayStdDev(n,ret)
plot(stddev20)
points(stddev60,col="red")
points(stddev180,col="yellow")
points(stddev_n,col="green")


# normalized returns:
normret15 = dDayNormRet(15,ret)
normret60 = dDayNormRet(60,ret)
normret250 = dDayNormRet(250,ret)
normret_n = ret/stddev_n[n]

par(mfrow=c(2,2))            # set up 2 time 2 plot array
                             # that is, 4 pictures at once

hist(normret_n,breaks=150,xlim=c(-5,5),ylim=c(0,0.8),prob=TRUE)
curve(dnorm(x),from=-5,to=5,add=TRUE,col="red")
hist(normret250,breaks=80,xlim=c(-5,5),ylim=c(0,0.8),prob=TRUE)
curve(dnorm(x),from=-5,to=5,add=TRUE,col="red")
hist(normret60,breaks=50,xlim=c(-5,5),ylim=c(0,0.8),prob=TRUE)
curve(dnorm(x),from=-5,to=5,add=TRUE,col="red")
hist(normret15,breaks=50,xlim=c(-5,5),ylim=c(0,0.8),prob=TRUE)
curve(dnorm(x),from=-5,to=5,add=TRUE,col="red")

# thus: if returns are normalized with more recent volatility 
# data, then normalized returns are more Gaussian 



#-----------------------------#
#  Analysis General Electric  #
#-----------------------------#

GE = read.table("C:/Users/detlef/OneDrive/hochschule/Vorlesungen/WS2223/Datenanalyse-mit-R/GE.txt",header=TRUE,sep=";")
head(GE)
tail(GE)
summary(GE)                 # no NA's
names(GE)
S = GE[,2]                  
plot(S,type="l")
plot(log(S/S[1]),type="l")  # straight line would be constant exponential growth
n=length(S)
n
# jetzt berechnen wir die returns: -> redo above S&P500 code, conclusion also holds.