Pracovní text a úkoly ke cvičením MF002 1 Wienerův proces (Brownův pohyb) Základním stavebním kamenem simulací náhodných procesů popsaných pomocí stochastických diferenciálních rovnic je generování trajektorií Wienerova procesu Wt = W(co, t). Definice říká, že počáteční hodnotou je W0 = 0, ale s jakoukoliv jinou známou hodnotou se pracuje podobně. Přírůstky Wienerova procesu za časové intervaly délky Ař jsou náhodné veličiny AW ~ N(0, Ař), a pro disjunktní časové intervaly jsou tyto přírůstky stochasticky nezávislé. Vlastní simulace trajektorie Wienerova procesu spočívá v dělení požadovaného časového intervalu [0, T] na velmi krátké podintervaly délek Ař. V těchto jemně zvolených časech postupně (s rostoucím časem) počítáme hodnotu Wienerova procesu tak, že k hodnotě v bezprostředně předchozím známém čase připočítáváme realizaci náhodného přírůstku AW. Pro hodnotu W(co, ŕ + Ař) tedy platí: W(co, ř + Ař) = W(co, ř) + AW , přičemž AW = í/v/ÄT, kde Í/~N(0,1), tj. realizace standardizované náhodné náhodné veličiny, jsou pro různá ř vzájemně nezávislé. Projděte si následující skript, který Wf generuje pomocí for-cyklu: 1 dt <- 0.001 2 W0 <- 0 3 t <- seq (0, 1, by=dt) 4 s N <- length (t) 6 W <- rep (0 , N) 7 W [1] <- W0 s Wakt <- W0 9 for (k in 2:N) { ío dW <- sqrt (dt) * rnorm (1) n Wakt <- Wakt + dW 12 W[k] <- Wakt 13 } 14 15 plot (t, W, type="l", col="red", xlab="t", ylab="W") 16 abline (h = W0 , lty = 2) Výše uvedený postup se však v moderních matematických programech (R, Matlab) označovaných jako SVM (Support Vector Machine) nepoužívá. Místo toho se jen jedním průchodem generuje celý náhodný výběr (vektor) přírůstků (nezávislost z definice Wf), který se kumulativně přičítá (funkcí cumsum) k výchozí hodnotě W0: 1 dt <- 0.001 2 W0 <- 0 3 t <- seq (0, 1, by=dt) 4 s N <- length (t) 6 dW <- rnorm (N 1) * sqrt (dt) 7 W <- cumsum (c (W0 , dW)) 8 Ondřej Pokora, ÚMS PřF MU Brno, aktualizace 3. března 2014 1 9 plot (t, W, type="l", col="red", xlab="t", ylab="W") 10 abline (h = W0 , lty = 2) 0.0 0.2 0.4 0.6 0.8 1.0 t Řádky skriptu generující Wt je vhodné si uložit jako vlastní funkci. Následně totiž můžeme jen změnou parametrů této funkce snadno generovat nové trajektorie: i generuj.Wp <- function (t, dt, W0) { dW <- rnorm (length (t) 1) * sqrt (dt) 3 W <- cumsum (c (W0, dW)) 4 } 5 6 W <- generuj.Wp (t, dt, W0) 7 8 plot (t, W, type = "1", col = "red", xlab = "t", ylab = "W") 9 abline (h = W0 , lty = 2) Výhoda vytvořené funkce se především ukáže, když chceme vygenerovat více trajektorií: 1 M <- sapply (1:10 , function (k) { 2 generuj.Wp (t, dt, W0) 3 }) 4 s matplot(t,M,type="l",lty=l,xlab="t",ylab="W",main="trajektorie W") 6 abline (h = W0 , lty = 2) 2 trajektorie W 0.0 0.2 0.4 0.6 0.8 1.0 t Pouhou změnou rozsahu prvního parametru funkce sapply lze měnit počet trajektorií. Výsledné vektory hodnot pro jednotlivé trajektorie se přitom skládají vedle sebe, takže výsledná proměnná je matice. Čemu v této matici odpovídají jednotlivé řádky, resp. jednotlivé sloupce? Zkoušejte si vygenerovat jiné počty trajektorií (až 1000, 10 000) Wt. Vykreslete je do grafu, seznamte se s parametry příkazů plot a matplot pro rozsahy a popisy souřadnicových os, a pro volbu stylu, tloušťky, symbolu a barvy čar zobrazených dat. Měli byste dostat podobný graf: trajektorie W n-1-1-1-1-r 0.0 0.2 0.4 0.6 0.8 1.0 t Dále budeme zkoumat rozdělení pravděpodobnosti realizací Wt projeden konkrétní čas ř = 0.6. Nejprve potřebujeme zjistit, který řádek tomuto času odpovídá. 3 1 index <- which (t==0.6) 2 vyber <- M[index,] Vybrali jsme z matice trajektorií jen jeden řádek odpovídající požadovanému času t = 0.6, který je z pohledu statistiky náhodným výběrem. Následující příkazy shrnují nejčastěji používané nástroje pro ověřování rozdělení pravděpodobnosti: histogram, QQ-plot, empirickou distribuční funkci (ecdf) a Kolmogorovův-Smirnovův test. Pomocí nápovědy prozkoumejte použití příkazů a především interpretaci výsledků. 1 index <- which (t == 0.6) 2 vyber <- M[index,] 3 4 hist (vyber, breaks = 50, freq = FALŠE) 5 6 qqnorm (vyber, pch = "+") 7 qqline (vyber, col = "blue") 8 9 plot (ecdf (vyber)) 10 n x <- seq (-2 ,2 , by = 0.1) 12 y <- pnorm (x , mean = 0 , sd = sqrt (0.6)) 13 lineš (x, y, col = "red") 14 15 ks. test (vyber, pnorm, mean = 0 , sd = sqrt (0.6)) 16 ks.test (vyber, pnorm, mean = 1 , sd = sqrt (0.6)) 17 ks. test (vyber, pnorm, mean = 0 , sd = 1) Měli byste dostat podobné grafy: Histogram of vyber C/5 CD Q CD O o CM O o o -2 -1 0 H-rTU_n n -1 2 vyber 4 Normal Q-Q Plot 1.1 Úkoly Generujte postupně 10, 100, 1000, 10000 trajektorií Wt a vykreslujte je do grafů příkazem matplot. K matici realizací Wienerova procesu spočítejte střední hodnotu a směrodatnou odchylku pro každý čas (tedy pro každý řádek matice). K tomu se bude hodit funkce apply pro aplikaci na každý řádek či sloupec matice, a funkce mean pro střední hodnotu a sd pro směrodatnou odchylku. Do obrázku trajektorií pak přidejte křivku pro střední hodnoty a 95% interval spolehlivosti realizací. Měli byste dostat obrázky podobné následujícím, které jsou pro 100 a 1000 trajektorií: trajektorie W t 6 t Pro jednotlivé varianty (10, 100, 1000, 10000 trajektorií) si pro realizace Wienerova procesu v časech ŕ = 0.2 a t = 0.8 zobrazte histogram a QQ-plot. Vykreslete si také empirickou distribuční funkci a graficky ji porovnejte s teoretickou distribuční funkcí. Jaké rozdělení pravděpodobnosti a s jakými parametry má Wt pro t = 0.2 a t = 0.8? Definujme náhodné procesy Y(co, t) = min W(co,s) , s