#------------------------------#
#   Cramer-Rao Abschaetzung    #
#     fuer den OU-Prozess,     #
#        R-Simulation          #
#------------------------------#


###  Wir nehmen die Loesung8-Aufg2.txt als Vorlage  ###


N = 10000              # Anzahl Pfade
n = 10000              # Anzahl Zeitreihendaten pro Pfad = Anzahl Delta_t's in [0,T]

T = 80                 # hier jetzt 10, 20, 40 oder 80

kappa = 1              # 1 oder 10 in week9a.pdf
mu = 5                 # lassen wir fest
sigma = 1              # lassen wir fest



dt = T/n 

x = matrix(0,N,n+1)     # in der i-ten Zeile steht der i-te (von N) simulierte Pfad
                        # x_0 -> 1.Spalte, x_1 -> 2.Spalte, ... , x_n -> n+1 -te Spalte
x[ ,1] = rep(mu,N)


# Schleife ueber die Zeiten:
for( k in 2:(n+1) )
{
  phi = rnorm(N)
  x[,k] = x[,k-1] + kappa*(mu - x[,k-1])*dt + sigma*sqrt(dt)*phi
}



# Berechnung der Maximum-Likelihood-Schaetzer:

hkappa = rep(0,N)
hmu = rep(0,N)
hsigma = rep(0,N)

# Loop ueber alle Pfade:
for(i in 1:N)
{ 
  y = x[i,]
  a11 = sum( (y[1:n])^2 )
  a12 = sum( y[1:n] ) 
  A = rbind( c(a11,a12) , c(a12,n) )
  b1 =  sum( y[1:n]*y[2:(n+1)] )
  b2 =  sum( y[2:(n+1)] )
  b = c(b1,b2)
  res = solve(A,b)
  alpha = res[1]
  beta = res[2]
  hkappa[i] = (1-alpha)/dt
  hmu[i] = beta/(1-alpha)
  hsigma[i] = sqrt( 1/(n*dt) * sum( ( y[2:(n+1)] - (alpha*y[1:n] + beta) )^2 ) )
}


# wir merken uns das hkappa mit dem entsprechendem T, 
# genau eins auswaehlen:

hkappa10 = hkappa
hkappa20 = hkappa
hkappa40 = hkappa
hkappa80 = hkappa

# jetzt alle 4 Histogramme in einem Plotfenster:
par(mfrow=c(2,2))
hist(hkappa10,breaks=50,xlim=c(0,5),prob=TRUE)
hist(hkappa20,breaks=40,xlim=c(0,5),prob=TRUE)
hist(hkappa40,breaks=30,xlim=c(0,5),prob=TRUE)
hist(hkappa80,breaks=20,xlim=c(0,5),prob=TRUE)

# ..etwas andere Darstellung:
info = "Histogram  of  hkappa10,  hkappa20,  hkappa40,  hkappa80"
hist(hkappa80,breaks=50,xlim=c(0,5),prob=TRUE,col="lightyellow",main=info)
hist(hkappa40,breaks=80,xlim=c(0,5),prob=TRUE,add=TRUE,col="lightgreen")
hist(hkappa20,breaks=160,xlim=c(0,5),prob=TRUE,add=TRUE,col="turquoise")
hist(hkappa10,breaks=240,xlim=c(0,5),prob=TRUE,add=TRUE,col="lightblue")

# ..und mit umgekehrter Reihenfolge:
hist(hkappa10,breaks=240,xlim=c(0,5),ylim=c(0,2.5),prob=TRUE,col="lightblue",main=info)
hist(hkappa20,breaks=160,xlim=c(0,5),prob=TRUE,add=TRUE,col="turquoise")
hist(hkappa40,breaks=80,xlim=c(0,5),prob=TRUE,add=TRUE,col="lightgreen")
hist(hkappa80,breaks=50,xlim=c(0,5),prob=TRUE,col="lightyellow",add=TRUE)


# Und die numerischen Werte:

summary(hkappa10)
summary(hkappa20)
summary(hkappa40)
summary(hkappa80)