Opdracht 2. Deadline maandag 28 september 2015, 24:00 uur. Deze opdracht bestaat uit vier onderdelen; in elk onderdeel wordt gevraagd een Matlabprogramma te schrijven. De vier bijbehorende bestanden stuur je naar
[email protected]. Er komt eerst wat uitleg over random number generators, en schatters met bijbehorende Matlab tips; daarna komt de beschrijving van de opdracht.
1 Random Number Generator Een random number generator (RNG) is een deterministisch algoritme dat een rij getallen U1 , U2 , . . . produceert die beschouwd mogen worden als onafhankelijke trekkingen uit de U (0, 1) verdeling (de uniforme, ook wel homogene, kansverdeling op (0, 1)). Een goede RNG is de basis van simulatiestudies. Wat betekent goed? 1. U1 , U2 , . . . passeert succesvol statistische toetsen. 2. Het gaat razendsnel, bv 10M getallen in één seconde. 3. Je moet de rij getallen kunnen herhalen! 4. Je laptop en de VU-computer moeten dezelfde getallen produceren. 5. Omdat computers eindige machines zijn, zijn er maar eindig veel machinegetallen (zie college 2). Dus zal elke RNG die door een computerprogramma wordt uitgevoerd een rij U1 , U2 , . . . produceren die zich op den duur herhaalt. Bij een goede RNG duurt het heel lang voordat deze herhaling optreedt (de zogeheten herhalingsperiode). Matlab heeft een aantal goede RNG algoritmes beschikbaar die aan al deze eisen voldoen. De standaard RNG in Matlab is de Mersenne Twister met herhalingsperiode van 219937 − 1 (zie eis 5). Het eerste getal dat geproduceerd wordt door een RNG, is afhankelijk van hoe het RNG algoritme wordt gestart. Dat gebeurt door een zogeheten seed, wat niets anders is dan een getal (of soms een vector van getallen). Telkens als je Matlab opstart staat die seed op hetzelfde getal, en dan krijg je dezelfde rij U1 , U2 , . . . als bij je vorige Matlabsessie. Je kunt je simulatie-experimenten ook telkens laten beginnen met dezelfde seed, namelijk door die steeds op hetzelfde getal te zetten (zie eis 3).
1.1 Matlabtip 1 De aanroep rand genereert een trekking uit de U (0, 1) verdeling. Wil je veel random getallen, bv n = 106 , dan gaat het veel sneller om die gevectorizeerd te genereren dan een voor een. Experimenteer zelf:
1
n = 10^6; tic; u=rand(1,n); toc tic; for i=1:n u = rand; end toc De seed van de RNG in Matlab kun je zetten op een natuurlijk getal, bv 127528, door rng(127528). Je kunt de seed ophalen door de aanroep sd = rng, zie de volgende Matlabsessies in het commandoscherm uitgevoerd vanuit de default waarde van de seed. >> x = rand(1,5) x = 0.8147 0.9058
0.1270
0.9134
0.6324
>> x = rand(1,5) x = 0.0975 0.2785
0.5469
0.9575
0.9649
>> s = rng; x = rand(1,5) x = 0.8147 0.9058 0.1270
0.9134
0.6324
>> rng(s); x = rand(1,5) x = 0.8147 0.9058 0.1270
0.9134
0.6324
>> rng(127528); x = rand(1,5) x = 0.2328 0.2870 0.3777
0.0060
0.4808
2 Schatten van Momenten van Stochasten Laat X een stochastische variabele zijn met cumulatieve kansverdelingsfunctie F (x) = P(X ≤ x), met verwachting µ = E[X], en variantie σ 2 = Var(X). In de vakken Probability Theory en Statistics heb je geleerd hoe je µ en σ 2 kunt schatten:
2
Als X1 , X2 , . . . , Xn gelijkverdeelde, onafhankelijke trekkingen zijn uit de kansverdeling van X, dan is n 1X Xi (1) µ b= n i=1 een zuivere schatter van µ [dwz E µ b = µ]. En n 2 1 X S = Xi − µ b n−1 2
(2)
i=1
is een zuiver schatter van σ 2 . Een speciaal geval is de Bernoulli stochast: B = 1 met kans p en B = 0 met kansP 1 − p. Omdat
E[B] = p, kun je p schatten uit de 0/1 waarnemingen B1 , . . . , Bn door pb = n1 ni=1 Bi . Een toepassing hiervan is een stochastische variabele X op R waarvan je de kans P(X ≤ c) wilt schatten voor een gegeven c. Als X1 , . . . , Xn realisaties van X zijn, dan is een schatting van P(X ≤ c): n 1X 1{Xi ≤ c}. n i=1
2.1 Matlabtip 2 Computersimulatie (bv in Matlab of Ox) is een methode om onafhankelijke gelijkverdeelde waarnemingen X1 , . . . , Xn te genereren van een stochastische variabele X. In het tweede deel van Numerical Methods (periode 2) wordt nader uitgelegd hoe dat gaat. Voor deze opdracht is voldoende dit te doen voor een 0/1 stochast. Dat wil zeggen, laat A een gebeurtenis zijn die optreedt met kans p, of niet optreedt met kans q = 1 − p. De bijbehorende stochastische variabele is de Bernoulli (p), X = 1{A} met
P(X = 1) = p, P(X = 0) = 1 − p. Laat U = U (0, 1) de uniforme (0, 1) stochast zijn waarvan we trekkingen krijgen door Matlab’s rand (zie tip 1). Omdat p = P(U < p) krijgen we een waarneming X = 1 met kans p door (rand
2.2 Matlabtip 3 Zoals je weet kun je de variantieschatter (2) herschrijven tot n n 1 X 2 Xi − µ b2 . S2 = n−1 n i=1
Dat geeft dat er twee manieren zijn om de schatters µ b en S 2 te implementeren.
3
(3)
(i). Genereer eerst alle data X1 , . . . , Xn ; sla op in een vector; en bereken dan formules (1) en (2). Als je in Matlab deze vector x hebt genoemd, volstaan de Matlabfuncties mean(x) en var(x). (ii). Na elk volgend datapunt Xk , bereken je direct via recursie k X i=1
Xi =
k−1 X
Xi + Xk ,
en
i=1
k X i=1
Xi2 =
k−1 X
Xi2 + Xk2 .
i=1
Bereken na de laatste xn pas de formules (1) en (3). Voorbeeld van de muntworp om verwachting en variantie te schatten. function methode_i start; function start n = 10^6; p = 0.5; x = (rand(1,n)
function methode_ii start; function start n = 10^6; p = 0.5; xsom = 0; x2som = 0; for i=1:n x = (rand
4
Opdracht 2: How To Gamble If You Must. Veronderstel dat je deelneemt aan een investeringsplan. Het plan heeft een looptijd van n = 60 maanden; aan het begin van elke maand leg je een bedrag in, zeg xk in maand k. De beheerder van dit investeringsplan is een gewiekste belegger die het lukt om aan het eind van de maand (voordat je de volgende inleg doet) jou met kans p = 0.36 het ingelegde bedrag in drievoud terug uit te keren, maar met kans q = 1 − p = 0.64 ben je helaas je inleg kwijt. De p en q kansen zijn in elke maand hetzelfde, en daarnaast zijn de maandelijkse handelingen van de beheerder onafhankelijk van van eerdere uitkomsten. Dit plan is inderdaad gunstig want bij een inleg van één euro in een maand is het verwachte nettoresultaat −1 + 0.36 × 3 = 0.08. Dat is, je verwacht een rendement van 8% per maand. Veronderstel dat je instapt in dit plan met een vermogen van V0 = 1000 euro. De vraag is, hoeveel ga je per maand inleggen om aan het eind van het plan aan een bepaalde criterium te voldoen. Noteer Vk voor je vermogen aan het eind van maand k, vlak na een eventuele uitkering en nog voor je de volgende inleg doet. Dus bij inleg xk in maand k is (natuurlijk kun je hoogstens inleggen wat je vermogen is: xk ≤ Vk−1 ) ( Vk−1 + 2xk , met kans p; Vk = (4) Vk−1 − xk , met kans q. Definieer de groeifactors G1 , G2 , . . . door Gk =
1 ln Vk /V0 k
⇔
Vk = V0 ekGk .
Criteria kunnen zijn: (a.) max
E[Vn ]: maximaal verwacht vermogen aan het eind van de looptijd.
(b.) min
Var[Vk ]: minimale volatiliteit van het vermogen gedurende de looptijd.
(c.) max
P(Vn > γ): maximale kans om een streefbedrag te bereiken, bv γ = 5000 euro.
(d.) min P(Vk < δ): minimaal risico (te weinig vermogen) gedurende de looptijd, bv δ = 100 euro. (e.) max
E[Gn ]: maximale groeifactor.
In de volgende onderdelen ga je door middel van simulatie een aantal van deze criteria onderzoeken, onder de volgende aanname De inleg xk is een constante fractie α van het huidige vermogen Vk−1 .
Opdracht (a) Simuleer een realisatie van het vermogen Vk , k = 0, 1, . . . , n volgens de recursie (4). Doe dit achtereenvolgens voor α = 0.05, 0.2, 0.5, gebruik makend van common random numbers; dat wil zeggen, voor alle α’s gebruik je dezelfde seed van de RNG. Schrijf een Matlabprogramma
5
om in één figuur de grafieken van deze drie realisaties te plotten. Geef de drie realisaties verschillende kleuren en voeg legenda toe. Schrijf je Matlabprogramma in een scriptfile genaamd opdr2a_achternaam.m (vervang achternaam door één van jullie twee achternamen met alleen kleine letters).
Opdracht (b) Simuleer 100000 realisaties van het vermogen Vk , k = 0, 1, . . . , n volgens de recursie (4), en (i) sla alle 100000 eindvermogens Vn , i = 1, . . . , 100000 op in een kolomvector. Doe dit achtereenvolgens voor α = 0.05 en α = 0.2. Zet de twee kolomvectoren in een matrix V . Schrijf een Matlabprogramma om in één figuur de twee histogrammen van deze twee datavectoren weer te geven [gebruik weer histc en bar toegepast op de matrix V ]. Gebruik als domein van de histogrammen het interval [0, 4000] en neem de binbreedte 500. Schrijf je Matlabprogramma in een scriptfile genaamd opdr2b_achternaam.m (vervang achternaam door één van jullie twee achternamen met alleen kleine letters).
Onderdeel (c) p Voor α ∈ [0.01, 0.1], schat het verwachte eindvermogen E[Vn ], en de standaardafwijking Var(Vn ) gebaseerd op 10000 gesimuleerde eindvermogens. Noem de twee schattingen µ b(α) en σ b(α). Schrijf een Matlabprogramma om in één figuur de grafieken van deze µ b(α) en σ b(α) te plotten als functie van α. Geef de grafieken verschillende kleuren en voeg legenda toe. Schrijf je Matlabprogramma in een scriptfile genaamd opdr2c_achternaam.m (vervang achternaam door één van jullie twee achternamen met alleen kleine letters).
Onderdeel (d) Voor α ∈ [0.01, 0.1], schat de kans dat minstens één keer Vk < 100 wordt bereikt, ook weer gebaseerd op 10000 gesimuleerde realisaties. Tevens, schat de verwachte groeifactor E[Gn ]. Noem de schattingen π b(α) en γ b(α). Schrijf een Matlabprogramma om de grafieken van µ b(α) en γ b(α) te plotten als functie van α. Elke grafiek komt in een afzonderlijk figuur [zet figure voorafgaand aan elke plot opdracht]. Schrijf je Matlabprogramma in een scriptfile genaamd opdr2d_achternaam.m (vervang achternaam door één van jullie twee achternamen met alleen kleine letters).
Conclusie Probeer nu voor jezelf te beredeneren wat een goede fractie α van je vermogen is om elke maand te beleggen.
6
