#--------------------------------# # # # week8a: wir reproduzieren # # die Bilder aus dem week7b # # # #--------------------------------# 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, 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