# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 13.1.2014 ## Thema: Reskalierte Markovketten und Diffusionsprozesse # # Stochastische "Differentialgleichungen" und ihre # Simulation via Euler(-Maruyama)-Verfahren # # Beispiele: # Irrfahrt -- Brownsche Bewegung (Reprise vom 9.12.2013) # Ehrenfestsches Urnenmodell -- Ornstein-Uhlenbeck-Prozess # Galton-Watson-Prozess -- Fellers verzweigende Diffusion # (s.a. Aufgabe 5 auf Blatt 10) # (2-Typ-)Wright-Fisher-Modell -- Wright-Fisher-Diffusion # (s.a. Aufgabe 5 auf Blatt 9) ################################################ ## ## Markovketten (mit Wertenbereich in den reellen Zahlen), Reskalierung ## ## Reskalierte Markovketten, Beispiel 1: Irrfahrt ## (vgl. auch Sitzung_9.12.2013.R, Abschnitt ## "Von der Irrfahrt zur Brownschen Bewegung ## -- zentraler Grenzwertsatz für zufällige Pfade") # Skalenparameter: N <- 20 # variieren, z.B. N <- 100; N <- 1000 pfadsim <- function(n) (cumsum(2*rbinom(n, size=1, prob=0.5)-1))/sqrt(n) zeiten <- (1:N)/N pfad <- pfadsim(N) plot(zeiten, pfad, ylim=c(-2,2), type="l", xlab="Zeit",ylab="", main=paste("Reskalierte Irrfahrtspfade:",N,"Schritte")) # ggf. mehrere Pfade einzeichnen: pfad2 <- pfadsim(N) points(zeiten, pfad2, type="l") ## ################## ## ## Reskalierte Markovketten, Beispiel 2: Ehrenfest-Modell # (vgl. Bsp. 6.14, Vorlesung vom 6.1.2014) # Skalenparameter: N <- 20 # variieren, z.B. N <- 100; N <- 1000 startwert <- round(N/2) pfadsimEM <- function(n) { pf <- numeric(n+1) pf[1]<-startwert for (i in 1:n) { # wenn aktuell pf[i] Kugeln links, so wechselt mit W'keit # pf[i]/N eine Kugel nach rechts: pf[i+1] <- pf[i]+ifelse(runif(1) <= pf[i]/N, -1, 1) } pf } pfad<-pfadsimEM(N) plot(0:N, pfad, xlab="Schritte", ylab="Anzahl links", type="l", main=paste("Ehrenfest-Prozess mit",N,"Kugeln"), ylim=c(0,N)) # ggf. mehrere Pfade einzeichnen: pfad2 <- pfadsimEM(N) points(0:N, pfad2, type="l") # Betrachte nun reskalierte Zeit- und Anzahl-Achsen: # Schritt entspricht Zeit 1/N, # eine Kugel entspricht "Masse" 1/sqrt(N), wir zentrieren ... x11() plot((0:N)/N, (pfad-N/2)/sqrt(N), xlab="Zeit", ylab="sqrt(N)*(Anz. links - N/2)", ylim=c(-2,2), type="l", main=paste("Renormierte Ehrenfest-Kette, N=",N,sep="")) points((0:N)/N, (pfad2-N/2)/sqrt(N), type="l") pfad2 <- pfadsimEM(N) # ggf. mehrere Pfade einzeichnen ## ############### ## ## Reskalierte Markovketten, Beispiel 3 (vgl. Aufgabe 5 auf Blatt 10): ## kritischer (Bienayme-)Galton-Watson-Prozess ## (Francis Galton, 1822-1911; Henry William Watson, 1827-1903; ## Irénée-Jules Bienaymé, 1796-1878) ## Aus Jane Austen (1775-1817), "Persuasion": ## "Sir Walter Elliot, of Kellynch Hall, in Somersetshire, was a man who, ## for his own amusement, never took up any book but the Baronetage; ... ## there his faculties were roused into admiration and ## respect, by contemplating the limited remnant of the earliest patents ..." # Skalenparameter: N <- 20 # variieren, z.B. N <- 100; N <- 1000 startwert <- N #1 #N pfadsimGW <- function(n) { pf <- numeric(n+1) pf[1]<-startwert for (i in 1:n) { neue.gen <- 0 if (pf[i]>0) { for (j in 1:pf[i]) { # jedes Ind. hat 0 oder 2 Kinder, je mit W'keit 0.5 neue.gen <- neue.gen+2*rbinom(1,size=1,prob=0.5) } } pf[i+1]<-neue.gen } pf } pfad<-pfadsimGW(N) plot(0:N, pfad, xlab="Generation", ylab="Anzahl", type="l", main=paste("Galton-Watson-Prozess: Starte mit",N,"Ind."), ylim=c(0,2*N)) # ggf. mehrere Pfade einzeichnen: pfad2 <- pfadsimGW(N) points(0:N, pfad2, type="l") # Betrachte nun reskalierte Zeit- und Anzahl-Achsen: # eine Generation entspricht Zeit 1/N, # ein Individuum entspricht "Masse" 1/N x11() plot((0:N)/N, pfad/N, xlab="Generation", ylab="Anzahl", type="l", main=paste("Galton-Watson-Prozess: Starte mit",N,"Ind."), ylim=c(0,2)) points((0:N)/N, pfad2/N, type="l") ######################################### ## ## Bericht: ## ## Die betrachteten reskalierten Markov-Ketten verbindet folgende ## Struktur: ## Aufgefasst als zufällige Funktion X(t), wo t ein ## reeller "Zeitparameter" ist, ist für kleine Zeitzuwächse h ## (wörtlich jedenfalls, wenn t und t+h im betrachteten diskreten ## Zeitgitter liegen), ## E[X(t+h)-X(t)| X(t)=x] = h*mu(x) (ggf. + "Korrekturterm", der mit höherer ## Ordnung als h^1 nach 0 geht, also ## viel kleiner als const*h ist), ## Var[X(t+h)-X(t)| X(t)=x] = h*(sigma(x))^2 (ggf. + "Korrekturterm", der mit ## höherer Ordnung als h^1 nach 0 geht), ## für gewisse Funktionen mu(x) und sigma(x), und die reskalierten ## Pfade sind approximativ (dem Augenschein nach zumindestens) stetig. ## ## Diese Eigenschaften charakterisieren die sog. Diffusionsprozesse, ## die "lokal" wie eine Brownsche Bewegung aussehen: ## Wenn mu(x)=mu, sigma(x)=sigma nicht vom Ort x abhängen, ## so hat X(t) = mu*t + sigma*B(t) mit B(t) (Standard-)Brownsche Bewegung die ## gewünschten Eigenschaften, wenn mu(x), sigma(x) stetige Funktionen ## sind, so legen obige Eigenschaften zusammen mit dem zentralen ## Grenzwertsatz nahe, dass für kleine Zeitinkremente h die ## Inkremente eines solchen zufälligen Pfades X(t+h)-X(t), gegeben X(t)=x, ## (approximativ) normalverteilt sind mit Erwartungwert h*mu(x) ## und Varianz h*sigma(x)^2. ## ## Diese Tatsache liegt der Theorie der ## stochastischen Differentialgleichungen zugrunde, die solche ## Prozesse als Lösungen einer "Gleichung" ## (*) dX(t) = mu(X(t)) dt + sigma(X(t)) dB(t) ## gewinnt und studiert. ## Die mathematisch rigorose Form dieser Theorie ist nicht Thema des ## Stochastik-Praktikums (sie übersteigt unsere technischen Möglichkeiten ## an dieser Stelle), wir betrachten hier eine zeitdiskretisierte ## Approximation an die Gleichung (*), das sogenannte ## Euler(-Maruyama)-Schema: ## Für eine gegebene Zeitschrittlänge h>0 definiert man ## längs des Zeitgiters 0,h,2*h,3*h,... rekursiv ## (**) X((j+1)*h) := X(j*h)+h*mu(X(j*h))+sigma(X(j*h))*Z(j,h), j=0,1,2,... ## wo Z(j,h) u.a. zentriert normalverteilt mit Varianz h. ## ## (Man kann mathematisch rigoros zeigen, dass die Lösung von (**) in ## geeignetem Sinne gegen die von (*) konvergiert; wir werden unten ## zumindest an Simulationen illustrieren, dass sich (**) für kleiner ## werdendes h "stabilisiert".) ################################################ ## ## Erinnerung/Exkurs (zu Analysis und Numerik): ## Euler-Verfahren für Anfangswertaufgaben gewöhnlicher Differentialgleichungen ## ## Gesucht: numerische Approximation der Lösung x(t) für t aus [0,T.max] ## der Differentialgleichung ## ## dx(t)/dt = rechte.Seite(x(t)) ## mit vorgegebenem Startwert x(0) ### (Bem.: dies entspricht (*) oben mit sigma(x)=0, rechte.Seite(x)=mu(x) rechte.Seite <- function(x) { 2*x } startwert <- 1.0 T.max <- 1.5 ## Euler.Verfahren.det <- function(f=function(x) {2*x}, startwert=1.0, T.max=1.5, schrittweite=0.1, nur.Endwert=FALSE) { anz.schritte <- ceiling(T.max/schrittweite) approx <- numeric(anz.schritte+1) approx[1] <- startwert for (j in 1:anz.schritte) { approx[j+1]<-approx[j]+schrittweite*f(approx[j]) } if (nur.Endwert) { approx[anz.schritte+1] } else { approx } } ######################################################## ## Illustration: ## Wir betrachten als Beispiel die Dgl. dx(t)/dt = 2*x(t) ## (mit expliziter Lösung x(t)=x(0)*exp(2*t) explizite.Loesung <- function(t) startwert*exp(2*t) h.1 <- 0.1 pfad.1 <- Euler.Verfahren.det(T.max=T.max, schrittweite=h.1) zeitgitter.1 <- (0:(ceiling(T.max/h.1)))*h.1 h.2 <- 0.05 pfad.2 <- Euler.Verfahren.det(T.max=T.max, schrittweite=h.2) zeitgitter.2 <- (0:(ceiling(T.max/h.2)))*h.2 h.3 <- 0.01 pfad.3 <- Euler.Verfahren.det(T.max=T.max, schrittweite=h.3) zeitgitter.3 <- (0:(ceiling(T.max/h.3)))*h.3 farben <- c("red", "darkgreen", "blue") ## Betrachte exakte Lösung und Euler-Approximation mit ## verschiedenen Schrittweiten: curve(explizite.Loesung(x), xlim=c(0,T.max), xlab="t", ylab="x", n=300, main="Lösung von x\'(t)=2x(t) und Euler-Approximationen") points(zeitgitter.1, pfad.1, pch=18, cex=0.8, col=farben[1]) points(zeitgitter.2, pfad.2, pch=18, cex=0.8, col=farben[2]) points(zeitgitter.3, pfad.3, pch=18, cex=0.8, col=farben[3]) legend("topleft", pch=18, col=farben, legend=c(paste("Schrittw.",h.1), paste("Schrittw.",h.2), paste("Schrittw.",h.3))) ## ## Betrachte Güte der Approximation ## (gemessen an der Abweichung zum Endzeitpunkt) ## als Funktion der Schrittweite h.werte <- c(0.1,0.03,0.01,0.003,0.001, 3e-4, 1e-4, 3e-5, 1e-5) endwerte <- numeric(length(h.werte)) abweichungen <- numeric(length(h.werte)) for (j in 1:length(h.werte)) { endwerte[j] <- Euler.Verfahren.det(T.max=T.max, schrittweite=h.werte[j], nur.Endwert=TRUE) abweichungen[j] <- abs( explizite.Loesung(T.max) - endwerte[j]) } numExperiment <- data.frame(Schrittw=h.werte, Endwert=endwerte, Abw=abweichungen) # Resultat als data frame anschauen: numExperiment # und als Grafik: attach(numExperiment) plot(-log10(Schrittw), log10(Abw), xlab="-log10(Schrittweite)", ylab="log10(Fehler)", main=c("Euler-Verfahren: Fehler in Abhängigkeit der Schrittweite", "(log-log-Skala)")) anpassung <- lm(log10(Abw)~ I(-log10(Schrittw))) # Bem.: Wir benutzen in der Definition obigen linearen Modells # die "Identitätsfunktion" I(...) (s.a. ?I ), # da R sonst die Formel log10(Abw) ~ -log10(Schrittw) in lm(...) als # kleinste-Quadrate-Anpassung _ohne_ Absolutglied interpretiert. anpassung # die Steigung ist (nahezu) -1, # das Euler-Verfahren hat Ordnung 1: # absoluter Fehler <= Const./(Schrittweite)^1 abline(anpassung) detach(numExperiment) #################################################### ## ## Euler-(Maruyama-)Verfahren zur Simulation der Lösung von ## stochastischen Differentialgleichungen (siehe (**) oben) ## (Bem.: Gisiro Maruyama, 1916-1986, japanischer Mathemematiker) Euler.verf <- function(startwert=0, T.max=1.0, zeitschritt, mu=function(x) { 0 }, sigma=function(x) { 1 }, rand=c(NA, NA), randtyp="absorb", randeps=1e-4, nur.Endwert=FALSE) { anz.schritte <- ceiling(T.max/zeitschritt) waz <- sqrt(zeitschritt) X <- numeric(anz.schritte+1) X[1] <- startwert for (j in 1:anz.schritte) { X[j+1] <- X[j]+zeitschritt*mu(X[j])+sigma(X[j])*rnorm(1, sd=waz) # Randbehandlung: unterer Rand if (!is.na(rand[1]) && X[j+1] <=rand[1]) { if (randtyp=="reflekt") { X[j+1] <- rand[1]+randeps } else { # für absorbierenden Rand: Rest das Pfads=Randwert X[(j+1):(anz.schritte+1)] <- rand[1] break } } if (!is.na(rand[2]) && X[j+1] >=rand[2]) { if (randtyp=="reflekt") { X[j+1] <- rand[2]-randeps } else { # für absorbierenden Rand: Rest das Pfads=Randwert X[(j+1):(anz.schritte+1)] <- rand[2] break } } } if (nur.Endwert) { X[anz.schritte+1] } else { X } } ######### ## ## SDGl mit Euler-Verfahren simulieren: Beispiele delta <- 0.01 ## SDGl-Beispiel: Brownsche Bewegung name <- "Brownsche Bewegung" x0 <- 0 # Startwert T.max <- 2.0; zeitgitter <- delta*(0:(ceiling(T.max/delta))) mu.x <- function(x) { 0 } sigma.x <- function(x) { 1 } pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x) plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(-3,3)) # ggf. weitere Pfade hinzu"malen": pfad2 <- Euler.verf(x0, T.max, delta, mu.x, sigma.x) points(zeitgitter, pfad2, type="l") # Betrachten wir die Verteilung zum Endzeitpunkt: wdh <- 100 endwerte <- replicate(wdh, Euler.verf(x0, T.max, delta, mu.x, sigma.x, nur.Endwert=TRUE)) hist(endwerte, prob=T) # in diesem Beispiel sollte Normalvert. mit Varianz T.max herauskommen: # (hier ist das Euler-Schema sogar "exakt") qqnorm(endwerte) ## ## SDGl-Beispiel: geometrische Brownbewegung name <- "geometrische Brownbewegung" x0 <- 1.0 # Startwert T.max <- 10.0; zeitgitter <- delta*(0:(ceiling(T.max/delta))) rendite <- 0.04; volatilitaet <- 0.14 # ggf. variieren mu.x <- function(x) { rendite*x } sigma.x <- function(x) { volatilitaet*x } pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA),randtyp="reflekt") plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(0,3)) # ggf. weitere Pfade hinzu"malen": pfad2 <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA),randtyp="reflekt") points(zeitgitter, pfad2, type="l") # Betrachten wir die Verteilung zum Endzeitpunkt: wdh <- 500 endwerte <- replicate(wdh, Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA), randtyp="reflekt", nur.Endwert=TRUE)) hist(endwerte, prob=T) hist(log(endwerte), prob=T) # in diesem Beispiel sollte log-Normalvert. herauskommen: # (exakt im Limes Schrittweite -> 0) qqnorm(log(endwerte)) qqnorm(endwerte) # zum Vergleich (endwerte selbst ist nicht normalvert.) ## ## SDGl-Beispiel: Ornstein-Uhlenbeck-Prozess name <- "Ornstein-Uhlenbeck-Prozess" x0 <- 1.5 # Startwert T.max <- 7.0; zeitgitter <- delta*(0:(ceiling(T.max/delta))) ruecktrpar <- 1.2 # ggf. variieren mu.x <- function(x) { -ruecktrpar*x } sigma.x <- function(x) { 1 } pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x) plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(-3,3)) # ggf. weitere Pfade hinzu"malen": pfad2 <- Euler.verf(x0, T.max, delta, mu.x, sigma.x) points(zeitgitter, pfad2, type="l") # Betrachten wir die Verteilung zum Endzeitpunkt: wdh <- 100 endwerte <- replicate(wdh, Euler.verf(x0, T.max, delta, mu.x, sigma.x, nur.Endwert=TRUE)) hist(endwerte, prob=T) # in diesem Beispiel sollte eine Normalvert. herauskommen: # (mit EW = x0*exp(-ruecktrpar*T.max), # Var = (1-exp(-ruecktrpar*T.max))/(2*ruecktrpar) ) # (exakt im Limes Schrittweite -> 0) qqnorm(endwerte) # Wo war das Teilchen: Betrachte Okkupationsmaß # eines Pfads über längere Zeitstücke # (bei positivem ruecktrpar) x0 <- 0 ruecktrpar <- 1.2 T.max <- 100; zeitgitter <- delta*(0:(ceiling(T.max/delta))) pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x) plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(-3,3)) x11() hist(pfad, prob=T, main="Okkupationsmaß eines OU-Prozesses") ## ## SDGl-Beispiel: Wright-Fisher-Diffusion name <- "Wright-Fisher-Diffusion" x0 <- 0.6 # Startwert, ggf. variieren T.max <- 2.0; zeitgitter <- delta*(0:(ceiling(T.max/delta))) mu.x <- function(x) { 0 } sigma.x <- function(x) { sqrt(x*(1-x)) } pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,1)) plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(0,1)) # ggf. weitere Pfade hinzu"malen": pfad2 <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,1)) points(zeitgitter, pfad2, type="l") # Betrachten wir die Verteilung zum Endzeitpunkt: wdh <- 100 endwerte <- replicate(wdh, Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,1), nur.Endwert=TRUE)) hist(endwerte, prob=T) # Bem.: Vert. ist (relativ) kompliziert, # hat Atome in 0 und in 1 # (für T.max nach unendlich approximiert sie # Atom in 0 mit Masse (1-x0) plus Atom in 1 mit Masse x0 ## ## SDGl-Beispiel: Fellers verzweigende Diffusion name <- "Fellers verzweigende Diffusion" x0 <- 1.0 # Startwert T.max <- 4.0; zeitgitter <- delta*(0:(ceiling(T.max/delta))) varianzpar <- 1.0 # ggf. variieren mu.x <- function(x) { 0 } sigma.x <- function(x) { sqrt(varianzpar*x) } pfad <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA)) plot(zeitgitter, pfad, type="l", xlab="t", ylab="", main=c("Simulation mit Euler-Verfahren", paste("Beispiel:", name)), ylim=c(0,3)) # ggf. weitere Pfade hinzu"malen": pfad2 <- Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA)) points(zeitgitter, pfad2, type="l") # Betrachten wir die Verteilung zum Endzeitpunkt: wdh <- 100 endwerte <- replicate(wdh, Euler.verf(x0, T.max, delta, mu.x, sigma.x, rand=c(0,NA), nur.Endwert=TRUE)) hist(endwerte, prob=T) # Bem.: Vert. ist (relativ) kompliziert, hat Atom in 0 # (das Gewicht des Atoms in 0 ist (exakt im Limes Schrittweite -> 0) # exp(-x0/(varianzpar*T.max/2)) ) # Wir können z.B. die "Aussterbewahrscheinlichkeit bis T.max" # (d.h. die Masse des Atoms in 0) schätzen: sum(endwerte==0.0)/wdh ################################################################## ## ## Exkurs, falls Zeit ist: ## ## Illustrationen zur Konsistenz des ## stochastischen Euler-Schemas: ## betrachte eine Schar ineinandergeschachtelter Gitter ## und benutze dieselbe Brownsche Bewegung für alle Gitter gleichzeitig driftfkt <- function(x) { 0 } difffkt <- function(x) { 1 } startwert <- 0.0 T.max <- 2.0 # Zeithorizont ## ggf. variieren, z.B. # driftfkt <- function(x) { 0.04*x }; difffkt <- function(x) { 0.1*x }; startwert <- 1.0; T.max <- 3.0 # Verfahren und Konsistenz auf Verfeinerung von Gittern betrachten: skalenexp.fein <- 15 # 4 #12 # ggf. variieren skalenexp.grob <- 3 # 1 # 3 # delta.t.fein <- T.max/2^skalenexp.fein anz.schritte.fein <- ceiling(2^skalenexp.fein) BBinkremente <- rnorm(anz.schritte.fein, 0, sd=sqrt(delta.t.fein)) delta.ts <- sapply(skalenexp.grob:skalenexp.fein, function(i) T.max/2^i) anz.schritte <- sapply(skalenexp.grob:skalenexp.fein, function(i) 2^i) Euler.schema <- vector("list", skalenexp.fein-skalenexp.grob+1) ## Euler-Schemata für verschiedene Skalen gleichzeitig betreiben: for (i in 1:(skalenexp.fein-skalenexp.grob+1)) { # stoch. Euler-Verfahren mit Schrittlänge delta.ts[i] Euler.schema[[i]]<-numeric(anz.schritte[i]+1) Euler.schema[[i]][1] <-startwert gru <- 2^(skalenexp.fein-skalenexp.grob-(i-1)) for (j in 1:(anz.schritte[i])) { BBinkr <- sum(BBinkremente[(gru*(j-1)+1):(gru*j)]) Euler.schema[[i]][j+1] <- Euler.schema[[i]][j]+ driftfkt(Euler.schema[[i]][j])*delta.ts[i]+ difffkt(Euler.schema[[i]][j])*BBinkr } } # Resultat der Schemata gemeinsam anschauen: farben <- rainbow(skalenexp.fein-skalenexp.grob+1) max.werte <- max(unlist(lapply(Euler.schema, max))) min.werte <- min(unlist(lapply(Euler.schema, min))) plot(0, xlim=c(0,T.max), ylim=c(min.werte, max.werte), xlab="t", ylab="", type="n") for (i in 1:(skalenexp.fein-skalenexp.grob+1)) { points((0:(anz.schritte[i]))*delta.ts[i], Euler.schema[[i]], type="l", col=farben[i]) } # "Farblegende" malen: x11() plot(0, type="n") legend("bottomleft", lty=1, col=farben, legend=skalenexp.grob:skalenexp.fein, cex=2) # Das Diagramm wird bei vielen Skalen gleichzeitig leicht # unübersichtlich: nur Auswahl zeichnen anz.skalen <- skalenexp.fein-skalenexp.grob+1 ausw <- c(1, round(anz.skalen/3), round(2*anz.skalen/3), anz.skalen) farben <- rainbow(skalenexp.fein-skalenexp.grob+1) #rainbow(length(ausw)) max.werte <- max(unlist(lapply(Euler.schema, max))) min.werte <- min(unlist(lapply(Euler.schema, min))) plot(0, xlim=c(0,T.max), ylim=c(min.werte, max.werte), xlab="t", ylab="", type="n") for (i in ausw) { points((0:(anz.schritte[i]))*delta.ts[i], Euler.schema[[i]], type="l", col=farben[i]) } # "Farblegende" malen: x11() plot(0, type="n") legend("bottomleft", lty=1, col=farben[ausw], legend=(skalenexp.grob:skalenexp.fein)[ausw], cex=2) ## folgendes ist nur sinnvoll, wenn mu und sigma nicht ## konstant sind (sonst ist das Euler-Verfahren nämlich exakt). ## Zur Beobachtung der "Konvergenz" betrachten wir die Endwerte: ## (Euler.schema[[i]][length(Euler.schema[[i]])] ist der ## Wert zum Zeitpunkt T.max im Schema mit Schrittweite delta.ts[i]) anz.skalen <- skalenexp.fein-skalenexp.grob+1 plot(1:anz.skalen, sapply(1:anz.skalen, function (i) Euler.schema[[i]][length(Euler.schema[[i]])]), xlab="Skala", ylab="Wert zum Endzeitpunkt") ## ######################################################## ## # Bem.: Es gibt auch ein R contributed package "sde" # das u.a. vorgefertigte Simulationsroutinen für stochastische # Differentialgleichungen enthält. # (siehe CRAN, das comprehensive R archive network, # Link von http://r-project.org)