#-----------------------------------------#
#                                         #
#   week12b: Die Schaetzer aus week12a    #
#     vs Maximum-Likelihood-Schaetzer     #
#                                         #
#-----------------------------------------#

# Wir benutzen den code mit den ML-Schaetzern 
# aus dem week8a.txt als Vorlage:


N = 10000              # Anzahl Pfade
n = 1000               # Anzahl Zeitreihendaten pro Pfad = Anzahl Delta_t's in [0,T]
T = 10                  # 1 oder 10 in week7b.pdf

kappa = 1              # 1 oder 10 in week7b.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)

hkappa2 = rep(0,N)           # fuer die neuen Schaetzer 
hmu2 = rep(0,N)              # aus dem week12a, also mit 
hsigma2 = rep(0,N)           # einer 2 hinten dran


# Loop ueber alle Pfade:
for(i in 1:N)
{ 
  y = x[i,]                  # der i-te OU-Pfad
  
  # Maximum-Likelihood-Schaetzer:
  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 ) )
  
  # die neuen Schaetzer:
  hmu2[i] = a12/n
  dy = diff(y)
  ny = length(dy)
  sigmasquare = 1/(ny*dt)*sum( dy^2 )
  hsigma2[i] = sqrt( sigmasquare )
  hs2 = a11/n - (a12/n)^2
  hkappa2[i] = 1/T * Finv( 2*hs2/(T*sigmasquare) )
}


# mu:
par(mfrow=c(1,2))
info = "hmu: neu in gelb, alt in blau"

hist(hmu,breaks=50,prob=TRUE,col="lightblue",main=info)
hist(hmu2,breaks=50,prob=TRUE,add=TRUE,col="lightyellow")

hist(hmu2,breaks=50,prob=TRUE,col="lightyellow",main=info)
hist(hmu,breaks=50,prob=TRUE,add=TRUE,col="lightblue")

summary(hmu)
summary(hmu2)


# sigma:
par(mfrow=c(1,2))
info = "hsigma: neu in gelb, alt in blau"

hist(hsigma,breaks=50,prob=TRUE,col="lightblue",main=info)
hist(hsigma2,breaks=50,prob=TRUE,add=TRUE,col="lightyellow")

hist(hsigma2,breaks=50,prob=TRUE,col="lightyellow",main=info)
hist(hsigma,breaks=50,prob=TRUE,add=TRUE,col="lightblue")

summary(hsigma)
summary(hsigma2)


# kappa:
par(mfrow=c(1,2))
info = "hkappa: neu in gelb, alt in blau"

hist(hkappa,breaks=50,prob=TRUE,col="lightblue",main=info)
hist(hkappa2,breaks=50,prob=TRUE,add=TRUE,col="lightyellow")

hist(hkappa2,breaks=50,prob=TRUE,col="lightyellow",main=info)
hist(hkappa,breaks=50,prob=TRUE,add=TRUE,col="lightblue")

summary(hkappa)
summary(hkappa2)





# Die Funktionen F und Finv := F^(-1)

FF = function(y)                        # F ist schon verbraucht fuer FALSE
{
  res1 = (exp(-2*y)-1+2*y)/(2*y^2)
  res2 = ( 4*exp(-y)-exp(-2*y)-3+2*y )/(y^3)
  return(res1-res2)
}

Fprime = function(y)
{
  res1 = ( -(y+1)*exp(-2*y)-y+1 )/(y^3)
  res2 = ( (2*y+3)*exp(-2*y)-(4*y+12)*exp(-y)-4*y+9 )/(y^4) 
  return(res1-res2)
}

Finv = function( z )
{
  # x = F^(-1)(z)  <=>  F(x)-z = 0
  # Newton Iteration:
  x = 0.01
  for(i in 1:10)      # Newton konvergiert sehr schnell
  {
    x = x - (FF(x)-z)/Fprime(x) 
  }
  return(x)
}


y = seq(from=0.01,to=50,by=0.01)
plot(y,FF(y))

y = seq(from=0.01,to=5,by=0.01)
plot(y,FF(y))

y = seq(from=0.01,to=50,by=0.01)
plot(y,Fprime(y))

y = seq(from=0.01,to=50,by=0.01)
plot(y,Finv(FF(y)))

x = seq(from=0.01,to=0.33,by=0.001)
plot(x,FF(Finv(x)))