#-----------------------#
#    Simulation von     # 
#     N OU-Pfaden:      #
#-----------------------#

N = 20
T = 10
dt = 0.01

kappa = 1
mu = 5
sigma = 1
x0 = 2

t = seq(0,T,by=dt)
nt = length(t) 

X = matrix(0,nrow=N,ncol=nt)
X[,1] = x0

sqrdt = sqrt(dt)
EXmc = rep(0,nt)

for(k in 2:nt)
{
  phi = rnorm(N) 
  X[,k] = X[,k-1] + kappa*(mu-X[,k-1])*dt + sigma*sqrdt*phi 
}


# Monte Carlo Berechnung von E[Xt]: (N dann groesser waehlen)

for(k in 1:nt)
{
  EXmc[k] = sum(X[,k])/N
}

mycol = rainbow(N)   # wir machen die Pfade bunt..

plot(t, X[1,], type="l", ylim=c(0,10) )
for(i in 2:N)
{
  points(t, X[i,], type="l", col=mycol[i] )
}

EXtheo = mu + (x0-mu)*exp(-kappa*t)      # theoretischer Wert
points(t, EXtheo, col="red")
points(t, EXmc, col="green", cex=0.5)    # cex = character expansion, 
                                         # cex = 0.5 -> nur halb so gross