Metoda Monte Carlo

Pojištění

Pojďme si spočítat nějaký příklad, který se vyskytuje ve finanční matematice:
Sαn=11n1+α.
Vyjádříme si:
Sα=n=11n(n+1)n+1nα.
Z analýzy víme, že n=11n(n+1)=1, takže si můžeme zavést náhodnou proměnnou ξα s diskrétním rozdělením
P[ξα=n+1nα]1n(n+1).
Potom bude platit Sα=Eξα. Můžeme tedy použít aproximaci
ξ¯α1Ni=1Nξi,
σ2=1N(N1)i=1N|ξ¯αξi|2.

Obecně

Mějme náhodnou veličinu ξ se střední hodnotou Eξ a rozptylem Dξ=σ2<. Tu aproximujeme jako
ξ¯1Ni=1Nξi,
Eξ¯=μ,
Dξ¯=σ2N2.
Platí odhad:
P(|ξ¯μ|ε)σ2N2ε2,
který si můžeme upravit jako
P(|ξ¯μ|<σNη)1η.

Převod rovnoměrného rozdělení na diskrétní

Mějme diskrétní rozdělení, kde P[x=xi]=pi. Máme-li k dispozici názhodnou proměnnou ru(0,1) a chceme pomocí ní náhodně vybrat x, stačí najít i takové, že
j=1i1pjr<j=1ipj.
Speciálně pro Poissonovo rozdělení Po(λ) existuje lepší algoritmus. Budeme postupně generovat náhodné proměnné γiR(0,1) a skončíme, jakmile i=0nγiexp(λ). Potom nPo(λ).

Více exponenciálně rozdělených veličin

Chceme vygenerovat exponenciálně rozdělenou veličinu:
FX(x)=1exp(λx),x0
Můžeme ji vygenerovat obecným způsobem, tedy vezmeme γR(0,1) a spočteme si x=ln(1γ)λ. K tomu potřebujeme jeden výpočet logaritmu. Ovšem existuje způsob, jak vygenerovat několik nezávislých exponenciálně rozdělených veličin s jen jedním logaritmem. Mějme γ1,,γ2n1R(0,1). Vezměme γn+1,,γ2n1 a seřaďme je: α0=0α1=γněcoαn1=γněcoαn=1. Potom pro každé in^ platí (bez důkazu)
(αiαi1)ln(γ1γn)exp(x).

Gamma rozdělení

Chceme mít veličinu s rozdělením
ξαfα(x)=xα1exΓ(α).
Toto rozdělení má zajímavou vlastnost, že ξα+ξβξα+β.
Důkaz Hustota součtu náhodných veličin je konvoluce hustot. Zajímá nás tedy
(fα*fβ)(x)=1Γ(α)Γ(β)0xyα1ey(xy)β1e(xy)dy=exΓ(α)Γ(β)0xyα1(xy)β1dy=|y=txdy=xdt|=exxα+β2Γ(α)Γ(β)01tα1(1t)β1.

Normální rozdělení

Z centrální limitní věty víme, že pokud budeme pro veličiny XiR(0,1) sčítat XEXDX, dostaneme normální rozdělení N(0,1). Stačí si tedy spočíst EXi=12 a DXi=112. Ale tím samozřejmě dostaneme normální rozdělení jenom přibližně. K získání přesného normálního rozdělení můžeme použít legrační tvrzení, které jsme si odvodili na cvičení. Mějme náhodné veličiny γ1,γ2R(0,1), potom z nich můžeme získat dvě nezávislé normálně rozdělené veličiny
η12lnγ1cos2πγ2,
η22lnγ1sin2πγ2.

χ-rozdělení

fχ2(x)=2n2Γ(n2)xn21exp(x2)
fχ(x)=2n2+1Γ(n2)xn1exp(x22)

Vícerozměrné náhodné veličiny

Jak generovat náhodnou veličinu s vícerozměrným rozdělením F(x1,,xn)? Pokud dokážeme spočítat všechny potřebné integrály, můžeme prostě generovat složky postupně s tím, že je vždy podmíníme již vygenerovanými složkami. Ale to většinou moc nefunguje. Další možnost je zamítací metoda, ale ta bude většinou strašně pomalá, protože je malá šance, že se trefíme. Například u dvacetirozměrné koule je to 2.5×108. Ale v určitých případech to jde lépe.

Kovarianční matice

Máme danou kovarianční matici 𝐂 a chceme vygenerovat vektor s rozdělením ηN(μ,𝐂), to znamená
𝐂=E((ηEη)(ηEη)𝖳).
Pro jednoduchost uvažejme μ=0. Vynásobíme-li to nějakou maticí A, máme
𝐂(Aη)=E(Aη(Aη)𝖳)=E(Aηη𝖳A𝖳)=AE(ηη𝖳)A𝖳=AA𝖳=𝐂.

Výpočet integrálu metodou Monte Carlo

Mějme nezápornou omezenou integrabilní funkci f:a,b0,c. Chceme najít Iabf. První způsob by byl prostě střílet do obdélníku a,b×0,c a počítat, kolikrát se trefíme. Formálně bychom měli náhodnou proměnnou
χ[γ2γ1],(γ1,γ2)R(a,b×0,c).
Porom Eχ=Ic(ba). Náš odhad bude vypadat jako
θ1c(ba)Ni=1Nχi.
Dθ1=(c(ba))2NDχ=(c(ba))2N(Ic(ba)I2c2(ba)2)=IN(c(ba)I).
Můžeme na to ovšem jít míň hloupě. Stačí si říct, že kdybychom měli náhodnou proměnnou γR(a,b), potom Ef(γ)=I. Náš odhad bude
θ2baNi=1Nf(γi),γiR(a,b).
Potom Eθ2=(ba)Ef(γi)=I,
Dθ2==baNabf2,
Nba(Dθ1Dθ2)=cIabf2=ab(cf(x))f(x)dx0.
Takže za našich předpokladů je lepší druhá metoda. Navíc potřebujeme dvakrát méně náhodných veličin. Mohli bychom to ještě nějak vylepšit? Máme-li odhad g(x)f(x)g(x)+d, kde funkci g umíme přesně zintegrovat, můžeme použít metodu vyloučení hlavní části: metodou Monte Carlo si spočteme ab(fg) a přičteme abg. Také pokud víme, že nějaká část funkce bude výrazně chaotičtější než jiná, můžeme si vzít nějakou distribuční funkci p, náhodnou proměnnou ξp a počítat
Ef(ξi)=abp(x)f(x)dx,
θ3=1Ni=1nf(ξi)p(ξi)ξi,
Eθ3==I,
Dθ3=1NDf(ξ)p(ξ)=1N(E(f(ξ)p(ξ))2I2)=1N(abf2(x)p(x)dxI2)
Virius zkoušel různými metodami spočítat 01exp=e1. Nejlepší možnost byla použít odhad 1+xexpxe1+x a hustotu p(x)=1+x.

Vícerozměrné integrály

Chceme najít integrál IMf, kde Md. Je-li M d-rozměrný kvádr, můžeme to dělat stejně jako u jednorozměrného, ale pokud to bude nějaká hnusná množina, máme problém. Pro malé d můžeme použít zamítací metodu k vytvoření rovnoměrného rozdělení na M, ale jak už víme, pro velké d to funguje blbě, protože může být malá šance, že se trefíme. Často se dá využít nějaká symetrie.

Neomezená funkce

Co když máme funkci, která sice je integrabilní, ale má singularitu? Kdybychom to zkoušeli spočíst vzorkováním funkce, mohli bychom se náhodou trefit hodně blízko, což by nám naprosto rozhodilo výpočet. Místo toho můžeme zkusit na nějakém malém okolí singularity integrál odhadnout a počítat ho jenom na zbytku. Další možnost je vhodně zvolit p nebo provést nějakou substituci, aby singularita zmizela.

Integrály na neomezené množině

Co kdybychom chtěli spočíst například Iex2dx? Použijeme hustotu ξp(x)1π(1+x2). Potom budeme mít
θ=πNi=1Neξi2(1+ξi2),ξ1π(1+x2),
Dθ=πNe2x2(1+x2)dx<.
Takže tato metoda je použitelná, přestože integrujeme na nekonečné množině. Kdybychom ale tímhle způsobem zkusili spočíst I11+x2 s náhodnou proměnnou ξN(0,1), vyšlo by nám Dθ2=.

Integrál přes kouli

Dá se ukázet, že pokud chceme mít rovnoměrné rozdělení v kouli, tak si stačí v nějakých upravených kulových souřadnicích vzít φR(0,2π),μR(1,1),r???. Pokud bychom ale i r volili rovnoměrně, dostaneme integrál z 1r2, což se může někdy hodit. Značí-li U1 jednotkovou kouli, zkusme například zintegrovat
U1U2f(|PQ|)dPdQ.
Označíme-li ρ|PQ|, pro ρ0f(ρ)ρ2 konečnou hustotu. Vzhledem k symetrii můžeme pro kulové souřadnice volit φP=ϑP=0,rP3r2. (Toto rP můžeme vygenerovat tak, že zvolíme tři rovnoměrně náhodná čísla a vezmeme největší z nich.) Stejně tak můžeme volit φQ=0. Zbylé dvě proměnné zvolíme jako ϑqR(1,1),rQR(0,l), kde lμQrP+1rP(1μQ2). Potom
1=U12π121l1ρ2ρ2dφdrdμ.

Zpět k pojištění

Na začátku jsme počítali sumáž
Sαi=11n1+α.
To můžeme také počítat pomocí integrálu:
011+x1+αdx=1+111+x1+αdx.
To shora odhadneme 1x1+α a spočteme pomocí vzorkování plochy nebo vzorkování funkce.

Markovské procesy

Mějme proces s konečně mnoha stavy, kde stav, do kterého se dostaneme v dalším kroku, závisí jen na stavu, ve kterém jsem teď. Takový proces můžeme vyjádřit pomocí matice přechodu:
vk+1=(p1,1p1,npn,1pn,n)vk=Pvk,
kde (vk)i značí pravděpodobnost, že v k-tém kroku jsme v i-tém stavu. Stav i je absorpční, pokud pi,j=δi,j, tedy se z něj nejde dostat do jiného stavu. Stavy jsou sdružené, pokud existuje cesta z jednoho do druhého. Řetězec je zanikající, pokud každý stav je sdružený s nějakým absorpčním stavem. Řetězec je ergodický, pokud všechny stavy jsou sdružené.

Lineární rovnice

Mějme lineární rovnici Ax=b. Tu můžeme přepsat do tvaru x=Bx+b, kde BIA. Je-li ρ(B)<1, můžeme ji řešit metodou postupných aproximací:
x=i=0Bif,
xmbm+i=1nBm,ibi+i=1nj=1nBm,iBi,jbj+
Je-li matice stochastická (všechny řádky mají součet 1), můžeme použít metodu Monte Carlo. Chceme-li aproximovat například druhý člen, vygenerujeme náhodnou veličinu podle m-tého řádku, čímž dostaneme i. Poté vygenerujeme náhodnou veličinu podle i-tého řádku, čímž dostaneme j. Do výsledného součtu přihodíme bj. Co když ale matice není stochastická? A jak se vypořádat s tím, že máme nekonečně mnoho sčítanců? Rozepíšeme si BPF,b=pf, kde i:j=1nPi,j+pj=1. Potom máme
xm=pmfm+i=1nPm,ipiFm,ifi+i=1nj=1nPm,iPi,jpjFm,iFi,jfj.
Budeme to modelovat tak, že si přidáme jeden absorpční stav a vezmeme přechodovou matici
(P1,1P1,np1Pn,1Pn,npn001).
Budeme simulovat proces s touto přechodovou maticí. Začneme s jednotkovým vektorem em. Pokud se po k+1 krocích dostaneme do absorpčního stavu a trajektorie předtím je i1,,ik, za hodnotu náhodné veličiny vezmeme η(j=1kFij1,ij)fik. Střední hodnota potom bude hledané x. Ale jaká bude disperze?
Dη=Eη2(Eη)2,
Eη2=pmfm+i=1nPm,ipiFm,i2fi2+.
To znamená, že aby disperze byla konečná, musí konvergovat metoda postupných aproximací pro rovnici x=Cx+d, kde CFFP,dffp. Takže musíme P,p zvolit šikovně, na což neexistuje žádná obecná metoda.

Laplaceova rovnice

Mějme rovnici 2u=0 na oblasti G neboli x2u+y2u=0,u|G=f. Při použití diferenční náhrady dostaneme
ui,j=14(ui+1,j+ui1,j+ui,j+1+ui,j1).
Když to interpretujeme jako střední hodnotu, tak to můžeme vzít jako náhodný proces, kde začneme v bodě (i,j) a budeme se náhodně pohybovat, až narazíme na hranici, a vezmeme hodnotu f v bodě, kam narazíme. Snadno vidíme, že k získání konečné disperze stačí, aby f byla omezená. Mimochodem, funkce u splňující vlastnost 2u=0 se nazývá harmonická. Pro ty existuje věta o průměrné hodnotě, která říká, že hodnotu v nějakém bodě můžeme spočíst jako průměr něčeho kolem něj. To je vlastně spojitá verze toho, co jsme si teď odvodili.

Co kdybychom měli rovnici s pravou stranou, 2u=g? Její diskrétní verze je

ui,j=14(ui+1,j+ui1,j+ui,j+1+ui,j1)+h24gi,j.
No to je jednoduché – stačí v každém kroku ještě připočíst h24gi,j. Jaká je střední hodnota počtu kroků, které provedeme? Dostáváme rekurentní rovnici
Ki,j=14(Ki+1,j+Ki1,j+Ki,j+1+Ki,j1)+1.
To jde vyřešit pomocí nějaké analýzy vlastních čísel, ale do toho nebudeme zabíhat.

Něčí rovnice

φ(x)=nk(y,x)φ(y)dy+f(x).
Postupné aproximace:
φ(x)=nk(y,x)f(y)dy+f(x),
φ2(x)=nnk(y,x1)k(x1,x)f(y)dydx1+φ1(x).
Nějaká magie s Lebesgueovými prostory, jádrem, sdruženými prostory a dalšími funkcionálními šílenostmi.
φ=𝔎φ+f
𝔎=supf+𝔎ffsupxnn|k(y,x)|dy