#--------------------------------#
#                                #
#   Aufg2:  wir reproduzieren    #
#   die Bilder aus dem week9a    #
#                                #
#--------------------------------#



N = 10000              # Anzahl Pfade

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

T = 10                 # 1 oder 10 in week9a.pdf

kappa = 1              # 1 oder 10 in week9a.pdf, oder 0 fuer Brownsche Bewegung

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 ) )
}


hist(hkappa,breaks=50)
hist(hmu,breaks=50,xlim=c(0,10))
hist(hsigma,breaks=50)

summary(hkappa)
summary(hmu)
summary(hsigma)


# alle 3 Histogramme in einem 1x3 Plot-Array:
par(mfrow=c(1,3))
hist(hkappa,breaks=50)
hist(hmu,breaks=30,xlim=c(0,10))
hist(hsigma,breaks=50)



# Bild 20 Pfade:
plot(x[1,], type="l",ylim=c(2,8) )
for(i in 2:20)
{
  points(x[i,],type="l",col=i)
}






# Test Gleichungssystem loesen:
# x=2, y=3
#
# 4x - 3y = -1
#  x + 5y = 17

A = rbind( c(4,-3) , c(1,5) )
b = c(-1,17)
res = solve(A,b)
res                      # ok, das passt