Abifunktsioonid keskväärtuse ja veahinnangu arvutamiseks (kliki, kui tahad näha!)

Abifunktsioonid keskväärtuse arvutamiseks:

MC <- function(gen,g,n,alfa){
  X <- gen(n)
  Y <- g(X)
  EY <-  mean(Y)
  viga <- -qnorm(alfa/2)*sd(Y)/sqrt(n)
  return(c(EY,viga))
}
#abifunktsioon etteantud täpsusega arvutamiseks
MC_viga <- function(gen,g,alfa, viga){
  #arvutab vastuse etteantud veaga
  n <- 1000
  vastus <- MC(gen,g,n,alfa)
  while(vastus[2]>viga){
    #arvutame uue n väärtuse väikese varuga
    n <- ceiling(1.1*n*(vastus[2]/viga)^2)
    vastus <- MC(gen,g,n,alfa)
  }
  return(c(vastus,n))
}

Teema 1: MC meetodi rakendused

Finantsmatemaatika: optsioonide hinnad. Finantsoptsioonid on lepingud, mis võimaldavad omanikul saada tulevikus tulu, mille suurus sõltub aktsiaturu käitumisest vaadeldaval tulevikuperioodil. Näiteks ostuoptsioon, mis võimaldab osta ühe Amazon.com aktsia 3 kuu pärast hinnaga 3200 dollarit, võimaldab teenida tulu, kui turuhind on 3 kuu pärast suurem kui 3200; tulu suurus on sel juhul \(S(T)-3200\), kus \(S(T)\) tähistab aktsiahinda kolme kuu pärast. Kui tuleviku aktsiahind on väiksem kui 3200, siis on leping väärtusetu ja sellelt teenitav tulu on null. Seega saab tulevikus saadavat tulu arvutada funktsiooni \(p(s)=max(s-3200,0),\ 0\leq s<\infty\) abil kujul \(p(S(T))\). Optsioonide teooria üks põhitulemusi on see, et optsiooni turuhind sõlmimismomendil peab olema kujul \[H=E(e^{-r\,T} p(S(T)))),\] kus \(r\) on riskivaba intressimäär ning \(S(T)\) on nn riskineutraalsele turumudelile vastav tuleviku aktsiahind. Samasugune valem kehtib ka juhul, kui tulevikus teenitav tulu sõltub mitme aktsia tulevikukäitumisest. Lihtsamate turumudelite korral on \(S(T)\) jaotus teada.

Harjutus 1

Leia MC meetodil eelnevalt kirjeldatud amazon.com ostuoptsiooni hind täpsusega 1 dollar (tõenäosusega 0.99), kui hinna esituses kasutatav tulevikuhind \(S(T)\) on lognormaalse jaotusega \(LN(\mu,\sigma)\), kus \[\mu=\ln(3185)+(r-\frac{0.24^2}{2})\cdot T,\ \sigma=0.12,\] riskivaba intressimäär on 2% ja \(T=\frac{1}{4}\). NB! tulevikutulu arvutamisel (ehk funktsioni \(p\) defineerimisel) on otstarbekas kasutada R funktsiooni pmax, mis töötab korrektselt vektoritega (st kui \(s\) asemel anda ette vektor, tahame saada tulemuseks vektorit, mis iga etteantud väärtuse korral on leidnud sellele hinnale vastava funktsiooni väärtuse). Lognormaalse jaotusega juhuslike suuruste genereerimiseks on käsk rlnorm.

Kliki vastuse nägemiseks!

Õige hind nelja komakoha täpsusega on 152.8419

Kliki näidislahenduse nägemiseks! 

    r=0.02
    T=1/4
    mu <- log(3185)+(r-0.24^2/2)*T
    sigma <- 0.12
    r=0.02
    T=1/4
    alfa <- 0.01
    gen1 <- function(n){
      return(rlnorm(n,mu,sigma))
    }
    g1 <- function(s){
      return(exp(-r*T)*pmax(s-3200,0))
    }
    n <- 10000
    vastus1 <- MC(gen=gen1,g=g1,n,alfa)
    n <- ceiling(n*(vastus1[2]/1)^2)
    MC(gen=gen1,g=g1,n,alfa)

Kindlustusfirma laostumistõenäosus Kindlustusfirmal on alguses teatud rahasumma \(u\) (et maksta kahjunõudeid), olemasolevad kliendid maksavad kindlustuspreemiaid ning firma maksab välja kindlustujuhtumite kahjusid. Lisaks on muidugi veel firma püsikulud (töötajate palgad, kontorikulud jms), aga need jätame vaatluse alt välja. Eeldame, et preemiad laekuvad ajas pidevalt kiirusega \(c\) rahaühikut ajaühikus ja et kahjud \(X_i,\ i=1,2,\ldots\) on sõltumatud sama jaotusega juhuslikud suurused ning et kahjude toimumise ajad on juhuslikud suurused \(0\leq T_1<T_2<T_3\) jne (st välistame kahe kahju korraga toimumise). Siis ajamomendil \(t\) on kindlustusfirmal järgi rahasumma \[V_t=u+c\,t -\sum_{i:T_i\leq t} X_i\] Kindlustusfirma laostub, kui mingil ajamomendil muutub suurus \(V_t\) negatiivseks.

Harjutus 2

Olgu kindlustusfirma algne reserv \(u=50000\), preemiate laekumise kiirus \(c=50000\) (eurot aastas), olgu kahjude suurused lognormaalse joatusega \(LN(5,2)\) ning olgu kahjude toimumisaegade vahed \(T_1, T_2-T_1,T_3-T_2,\ldots\) sõltumatud juhuslikud suurused eksponentjaotusega \(Exp(45)\). Leida \(n=10000\) simulatsiooni põhjal tõenäosus, et kindlustusfirma läheb järgmise aasta jooksul pankrotti ning leida ka tõenäosuse hinnangu tõenäoline viga.

Vaata vihjet!

Kõige töökindlam, kuigi ka suhteliselt aeglane, on teostada simulatsioone kahekordse tsükliga: iga kord genereeri järjest nii kahjude vahele jäävaid aegasid kui kahjusid ning uuenda reservi väärtust, kuni järgmine kahju toimub juba väljaspool vaadeldavat ajaperioodi või kuni firma laostub.

Kliki, et näha oodatavat vastust!

Võimalik vastus on laostumistn=0.0883 tõenäoline_viga=0.00191383. Saadavate vastuste erinevused võiksid enamasti olla väiksemad kui kolmekordne tõenäoline viga.

Kliki, et näha näitelahendust! 

     c <- 50000
     u <- 50000
     gen2 <- function(n){
       laostumine <- rep(NA,n)
       for(i in 1:n){
         V<-u
         T <- 0
         while(T<1 &V>=0){
           vahe <- rexp(1,45)
           T <- T+vahe
           if(T<=1){
             kahju <- rlnorm(1,5,2)
             V<-V+c*vahe-kahju
           }
         }
         laostumine[i]=(V<0)*1
       }
     return(laostumine)
     }
     g2 <- function(x){return(x)}
     MC(gen = gen2,g=g2,n=10000,alfa = 1/2)

Kliki, et näha veidi efektiivsemat koodi! Allpool on toodud ka lahendus, kus generaatori sees on ainult üks tsükkel. Siin lahenduses tuleb kasuks funktsioon cumsum, mis tekitab vektorist osasummade vektori.

     c <- 50000
     u <- 50000
     gen2 <- function(n){
       laostumine <- rep(NA,n)
       m <- 100 #genereeritavate vaheaegade arv
       #peab olema piisavalt suur, et summa oleks suurem kui 1 väga suure tõenäosusega
       #mõelge ooteaegade summa prognoosiintervalli peale ...
       for(i in 1:n){
         vahed <- rexp(m,45)
         ajad <- cumsum(vahed)
         kahjud <- rlnorm(m,5,2)
         kahjusummad <- cumsum(kahjud)
         V <- u+c*ajad-kahjusummad
       #jätame järele ainult 1 aasta omad
         V <- V[ajad<=1]
         laostumine[i]=(min(V)<0)*1
       }
     return(laostumine)
     }
     g2 <- function(x){return(x)}
     MC(gen = gen2,g=g2,n=1000000,alfa = 1/2)

Teema 2: Antiteetiliste muutujate kasutamine

Antiteetilisteks nimetame juhuslikke suurused \(X_1\) ja \(X_2\), mis on sama jaotusega ja negatiivselt korreleeritud. Võimalused saada antiteetlisi juhuslikke suurused:

Kui on soovi leida MC meetodiga \(EY\), kus \(Y=g(X)\), siis antiteetiliste muutujate meetod seisneb MC meetodi rakendamises antiteetiliste muutujate \(X\) ja \(X_2\) abil defineeritud juhusliku suuruse \(\tilde{Y}=\frac{g(X)+g(X_2)}{2}\) kesväärtuse leidmiseks. Meetod annab ajalist võitu siis, kui \(g(X)\) ja \(g(X_2)\) on samuti antiteetilised juhuslikud suurused.

Harjutus 3

Kasuta ühtlast jaotust ja antiteetilisi juhuslikke suuruseid selleks, et arvutada integraal \[\int_1^3 e^{x^2-4x}\cos(x)\,dx\] veaga, mis on väiksem kui 0.0001 tõenäosusega 0.9. Leia, kui suur oli võit genereeritud juhuslike suuruste väärtuste arvus võrreldes tavalise MC meetodiga.

Kliki, et näha tavalise MC tulemust! Lahendus Kõigepealt tavalise MC meetodiga

g <- function(x){return(2*exp(x^2-4*x)*cos(x)*(x>=1)*(x<3))}
gen <- function(n){return(runif(n,1,3))}
vastus1 <- MC_viga(g=g,gen=gen,alfa=0.1,viga=0.0001)
vastus1
## [1]     -0.01767243585      0.00009338846
## [3] 339921.00000000000

Kliki, et näha antiteetiliste muutujate abil saadavat tulemust! 

gen2 <- function(n){
  X<-runif(n,1,3)
  X2 <- 4-X
  return(cbind(X,X2))
}
g2 <- function(M){
  X <- M[,1]
  X2 <- M[,2]
  return((g(X)+g(X2))/2)
}
vastus2 <-MC_viga(g=g2,gen=gen2,alfa=0.1,viga=0.0001)  
vastus2
## [1]   -0.01784947719    0.00009376251
## [3] 1484.00000000000
round(vastus1[3]/vastus2[3],1)
## [1] 229.1

Nagu näha, oli võit genereerimiste arvus 229.1-kordne. Kuna ühe väärtuse genereerimiseks tehakse antiteetiliste muutujate korral ligikaudu 2 korda rohkem arvutusi, on võit tööajas umbes 114.55-kordne

Harjutus 4

Kasuta antiteetiliste juhuslike suuruste ideed selleks et leida ligikaudselt keskväärtus \(E e^{\sqrt{X}}\) täpsusega \(0.001\) (juhul \(\alpha=0.05\)), kus \(X\) on eksponentjaotusega \(Exp(2)\) juhuslik suurus. Kui suur on võit vajaminevas genereerimiste arvus?

Kliki, et näha vihjet!

Eksponentjaotus ei ole sümmeetriline oma keskväärtuse suhtes. Seega kasuta antiteetiliste suuruste saamiseks jaotusfunktsiooni pöördfunktsiooni kaudu genereerimise ideed.

Kliki, et näha lahendust! 

n<- 1000000
gen <- function(n){return(rexp(n,2))}
g <- function(x){
  return(exp(sqrt(x)))
}
vastus1 <- MC_viga(gen=gen,g=g,alfa=0.05,viga=0.001)
vastus1
## [1]       1.9811095403       0.0009519438
## [3] 2321423.0000000000
gen2 <- function(n){
  U <- runif(n)
  return(cbind(qexp(U,rate=2),qexp(1-U,rate=2)))
}
g2 <- function(x){
  return((g(x[,1])+g(x[,2]))/2)
}
vastus2 <- MC_viga(gen=gen2,g=g2,alfa=0.05,viga=0.001)
vastus2
## [1]      1.9823306141      0.0009625894
## [3] 300987.0000000000
round(vastus1[3]/vastus2[3],1)
## [1] 7.7
Jällegi, kui tahta hinnata võitu tööajas, siis tuleks viimane number jagada kahega.

Antiteetiliste muutujate ideed saab kasutada ka mitmest juhuslikust suurusest sõltuva keskväärtuse leidmisel, kuid siis peab otsustama, kuidas on kõige suurem lootus saada keskväärtuse all oleva funktsiooni vastandlikke väärtuseid. Näiteks \(E(\sqrt{X^2+Y})\) leidmisel juhul \(X\sim N(0,1),\ Y\sim U(0,1)\) ei ole mingit mõtet asendada suurust \(X\) antiteetilise suurusega \(-X\), sest see ei ei anna tulemuseks vastandlikult käituvaid väärtuseid. Küll aga võib loota mingit võitu, kui lisaks paarile \((X,Y)\) kasutada paari \((X,1-Y)\), sest suurem \(Y\) väärtus vastab suuremale keskväärtuse all olevale funktsioni väärtusele ja seega saame vastandlikult käituvad tulemused. Mõnikord on kasulik asendada mitme juhusliku suurure väärtused antiteetilistega ning võimalikke üldistusi mitmemõõtmelisele juhule on veelgi (aga neid me siin kursusel ei vaatle)

Harjutus 5

Kasuta sobivat antiteetilistel juhuslikel suurustel põhinevat ideed, et leida keskväärtus \[E|X+Y-1|,\ X\sim U(0,\,2),Y\sim U(0,\,3)\] nii, et 10000 genereeritud arvupaari korral oleks vastus võimalikult täpne. Leida vastuse 95% usaldusintervall.

Leia võrdlemiseks leiame tulemuse tavalis MC meetodiga

Kliki, et näha koodi ja oodatavat vastust! 

gen<-function(n){
  X<- runif(n,min=0,max=2)
  Y<- runif(n,min=0,max=3)
  return(cbind(X,Y))
}
f<-function(x,y){
  return(abs(x+y-1))
}
g <-function(M){
  return(f(M[,1],M[,2]))
}
vastus1<-MC(gen = gen,g = g,alfa = 0.05,n = 10000)
vastus1
## [1] 1.55398922 0.01869039

Kliki, et näha vihjet, kuidas antiteetilisi väärtuseid tekitada!

Kuna tulemus sõltub summast, siis võiks kõige parem olla vaadata paare (X,Y) ja (2-X,3-Y)

Kliki, et näha vastust! 

gen_anti<-function(n){
  X<- runif(n,min=0,max=2)
  Y<- runif(n,min=0,max=3)
  X1 <- 2-X
  Y1 <- 3-Y
  return(cbind(X,Y,X1,Y1))
}
g_anti1 <-function(M){
  return((f(M[,1],M[,2])+f(M[,3],M[,4]))/2)
}
#genereerime poole vähem nelikuid, et oleks võrreldav
vastus2<-MC(gen = gen_anti,g = g_anti1,alfa = 0.05,n = 10000/2)
vastus2
## [1] 1.550951720 0.004165144

Vastav usaldusintervall on

print(sprintf("(%.3f, %.3f)",vastus2[1]-vastus2[2],vastus2[1]+vastus2[2]))
## [1] "(1.547, 1.555)"