Lezione 3 metodo monte carlo

Giovanni Della Lunga

Università degli Studi di Bologna

Introduzione al Metodo Monte Carlo

Concetti di Base e Applicazioni Finanziarie

Introduzione

Generatori di Numeri Casuali

Sequenze a Bassa Discrepanza (Quasi Monte Carlo)

Misure di co-dipendenza: correlazione e copule

Metodi di Riduzione della Varianza

Applicazioni al Pricing di Attività Finanziarie

Applicazioni alla stima del Value-at-Risk

Definizioni

Da un punto di vista del tutto generale con il termine Monte Carlo si intende una tecnica numerica che faccia uso di numeri casuali per risolvere un problema.

il metodo Monte Carlo consiste nel rappresentare la soluzione di un problema come parametro di un’ipotetica popolazione ed usare una sequenza di numeri casuali per costruire un campione della popolazione dal quale possano essere estratte stime statistiche del parametro.

Problemi risolubili col Metodo Monte Carlo A. Problemi di natura intrinsecamente probabilistica in cui sono

coinvolti fenomeni legati alla fluttuazione stocastica di variabili aleatorie. Due problemi tipici che rientrano in questa categoria sono il pricing di un’opzione e la stima del VaR di un portafoglio.

B. Problemi di natura essenzialmente deterministica, del tutto

privi cioè di componenti aleatorie, ma la cui strategia di soluzione può essere riformulata in modo tale da risultare equivalente alla determinazione del valore di aspettazione di una funzione di variabili stocastiche.

Simulazione Monte Carlo e Integrazione

Un’estrazione da un campione di numeri casuali può essere utilizzata come stimatore di un integrale

1

0

)( dxxfI 1

0

)( dxxfI

Questa espressione può essere interpretata come il valore di aspettazione della funzione f di una variabile aleatoria a valori uniformemente distribuiti nell’intervallo [0, 1]


Diventa così possibile stimare il valore del nostro integrale tramite una media aritmetica di n valori di f(xi) dove ciascun xi rappresenta un campione estratto da una distribuzione uniforme in [0, 1]. In altre parole possiamo affermare che la quantità

n

i

in xfn

I1

)(1~

n

i

in xfn

I1

)(1~

rappresenta uno stimatore non distorto di I. La varianza di questa stima risulta pari a:

n

dxIxfn

xfn

xfn

In

ii

n

iin

21

0

2

12

1

)(1

)(var1

)(1

var)~

var(

n

dxIxfn

xfn

xfn

In

ii

n

iin

21

0

2

12

1

)(1

)(var1

)(1

var)~

var(


l’errore quadratico medio dello stimatore, che può essere interpretato come l’errore quadratico medio della simulazione Monte Carlo, decresce all’aumentare di n come

Questo risultato risulta del tutto indipendente dalla dimensione del problema.

E’ proprio quest’ultima caratteristica che rende attraente il metodo Monte Carlo per la risoluzione di problemi con un numero elevato di dimensioni. In questo caso tipicamente il metodo Monte Carlo risulta convergere verso il valore finale più velocemente dei metodi numerici tradizionali in cui il numero di iterazioni per raggiungere un’approssimazione prefissata cresce con l’aumentare del numero di dimensioni.

n/1 n/1

Pricing di strumenti derivati

Il pricing di un’opzione è affrontato usualmente nel contesto della cosiddetta valutazione neutrale rispetto al rischio. Indicando con fT il valore dell’opzione stessa alla scadenza T, il valore ad oggi, f, sarà dato da

(1)

essendo Ê il valore di aspettazione risk-neutral ed <r> il valore medio fra t = 0 e t = T del tasso privo di rischio. Qualora si assuma di conoscere con certezza il valore di <r> il problema può essere semplificato e la (1) diventa

(2)

TrT efEf ˆ Tr

T efEf ˆ

TTr fEef ˆ T

Tr fEef ˆ

Pricing di strumenti derivati La formulazione del problema rende evidente la sua natura intrinsecamente

probabilistica. L’applicazione del metodo Monte Carlo nel caso in esame si riduce essenzialmente

alla generazione di un numero sufficientemente elevato di stime di fT da cui estrarre il valore medio.

A tal fine è necessario innanzi tutto introdurre un’ipotesi sul modo in cui il prezzo del titolo sottostante si evolve nel tempo;

Un’assunzione abbastanza comune è che il prezzo dell’azione segua un moto geometrico browniano. Secondo questa ipotesi il tasso di variazione del prezzo in un intervallo di tempo infinitesimo è descritto da

(3)

dove m è il tasso di rendimento atteso e s è la volatilità del prezzo dell’azione.

SdzSdtdS SdzSdtdS

Pricing di strumenti derivati

La simulazione viene condotta dividendo l’intervallo di vita del

derivato in N intervalli ciascuno di ampiezza t.

Si può dimostrare che la versione discreta della precednte

equazione è data da

(4)

dove S è la variazione di prezzo nell’intervallo t ed e è un

numero casuale estratto da una distribuzione normale

tStSS tStSS

Pricing di strumenti derivati A questo punto siamo in grado di calcolare i valori assunti da S (e quindi da S)

agli istanti 0, t, 2 t, ….fino alla scadenza T. Si noti che il processo di simulazione richiede la generazione di N numeri casuali

indipendenti normalmente distribuiti.

Una volta simulato il valore del titolo sottostante al tempo T siamo in grado di ricavare il valore dell’opzione alla stessa data. Supponendo che l’opzione sia di tipo call avremo semplicemente

essendo K lo strike price. Ripetendo la procedura appena descritta un numero molto elevato di volte siamo in grado di ottenere una distribuzione di valori per fT

dalla quale è possibile estrarre il valore di aspettazione.

)0,max( KSf TT )0,max( KSf TT

Pricing Strumenti Derivati

SdtSdtdS SdtSdtdS

dzdtSd

2)ln(

2

dzdtSd

2)ln(

2

Lemma di ItoLemma di Ito

tzt2S

SSSS

2

00

ln)ln()ln()ln( tzt

2S

SSSS

2

00

ln)ln()ln()ln(

Dall’ultima equazione

ricaviamo direttamente

questa è l’espressione che utilizzeremo nel programma VBA

Pricing Strumenti Derivati

tzt2S

S 2

0

ln tzt

2S

S 2

0

ln

tzt

2SS

2

0 exp

tzt

2SS

2

0 exp

Esempio ProgrammazioneVBA


“Crude” Monte Carlo“Crude” Monte Carlo

Introduzione







Tre tipologie di numeri “casuali” I generatori di numeri casuali che possono essere utilizzati nel metodo Monte

Carlo appartengono a tre categorie generali classificate in accordo alla natura dei numeri prodotti:

Numeri veramente casuali (truly random numbers)sono numeri la cui sequenza è assolutamente impredicibile (anche in via teorica) e devono essere prodotti sfruttando opportuni fenomeni fisici di cui è nota l’intrinseca aleatorietà (es. decadimento radioattivo di un nucleo atomico). Esistono serie di numeri di questo tipo registrate su nastri o pubblicate su appositi manuali.

Numeri pseudo casuali (pseudorandom numbers)sono numeri prodotti da un computer per mezzo di un algoritmo. Ovviamente essendo generati da un processo di calcolo deterministico essi non sono numeri veramente casuali in quanto la sequenza può essere predeterminata. Tuttavia è possibile simulare in maniera efficiente una sequenza di numeri che appaia del tutto casuale a chiunque non conosca l’algoritmo usato per generarla.

Numeri quasi casuali (quasirandom numbers)si tratta di numeri prodotti anch’essi da un algoritmo numerico, tuttavia in questo caso non si cerca di riprodurre una sequenza aleatoria bensì di disporre di una serie di numeri distribuiti nella maniera più uniforme possibile al fine di ridurre la varianza dei risultati delle integrazioni Monte Carlo.

Generatori di Numeri Pseudo-Causali

Esistono numeri casuali ? Come può un elaboratore, macchina totalmente

deterministica, generare numeri casuali e quindi per loro natura non deterministici?

La risposta è molto semplice: non può! I numeri sono generati per mezzo di qualche algoritmo

per cui non si può parlare di casualità essendo la sequenza predeterminata;

In compenso con un computer si possono generare sequenze di numeri che sembrino aleatorie

http://random.mat.sbg.ac.at/links/rando.html

Generatori di Numeri Pseudocasuali Virtualmente tutti i generatori di numeri pseudo casuali impiegati in

pratica sono basati sul generatore lineare congruente

I parametri a, c ed m determinano la qualità del generatore. a viene detto moltiplicatore, c incremento ed m è il cosiddetto modulo.

Il generatore appena visto genera numeri interi compresi fra 0 ed m. Usualmente si utilizzano generatori di numeri casuali uniformemente distribuiti fra 0 ed 1, per questo è sufficiente scegliere

mcaJJ ii mod 1 mcaJJ ii mod 1

mJU ii / mJU ii /

Vantaggi E’ molto veloce richiedendo pochissime operazioni per chiamata,

questo lo rende di uso universale;

Svantaggi Il più grosso svantaggio è rappresentato dalla presenza di

correlazione sequenziale; La sequenza di numeri casuali si ripeterà dopo un ciclo che, al

più, potrà essere di lunghezza m. Può produrre risultati inaspettati quando viene usato per la

generazione di distribuzioni non uniformi.

Generatore Lineare Congruente

Generatore Lineare Congruente Se si generano n coppie di numeri casuali e si

associano ad esse n punti in un piano, i punti non si distribuiscono uniformemente ma tendono ad allinearsi lungo segmenti di retta.

0

0 .1

0 .2

0 .3

0 .4

0 .5

0 .6

0 .7

0 .8

0 .9

1

0 0 .2 0 .4 0 .6 0 .8 1

0

0 .1

0 .2

0 .3

0 .4

0 .5

0 .6

0 .7

0 .8

0 .9

1

0 0 .2 0 .4 0 .6 0 .8 1

Generatore Lineare Congruente La correlazione sequenziale può essere

facilmente rimossa con tecniche di mescolamento (“shuffling”);

Il numero prodotto allo step j non costituisce l’output j-esimo ma viene utilizzato per l’output ad uno step successivo scelto in maniera casuale;

Generazione di distribuzioni Uniformi

Microsoft Excel La funzione Rnd() restituisce un valore numerico di tipo Single che contiene un

numero casuale. La sintassi è la seguente:

Rnd[(num)] L'argomento facoltativo num può essere un valore Single o una qualsiasi

espressione numerica valida. I valori restituiti dalla funzione dipendono dal valore passato come argomento.

Per ogni base iniziale specificata, viene generata la stessa sequenza di numeri, in quanto ogni successiva chiamata alla funzione Rnd() utilizza il numero casuale precedente come base per il numero successivo nella sequenza. In particolare se il parametro num è minore di zero Rnd() genera sempre lo stesso numero, utilizzando num come

base; se num è maggiore di zero viene restituito il successivo numero casuale nella sequenza; se num è uguale a zero viene restituito il numero generato per ultimo; infine se il parametro in input viene omesso, Rnd() restituirà il successivo numero casuale nella

sequenza.

Prima di richiamare Rnd(), è consigliabile utilizzare l'istruzione Randomize senza argomento per inizializzare il generatore di numeri casuali con una base connessa al timer del sistema con la seguente sintassi

Randomize[(numero)]

Randomize utilizza il parametro numero per inizializzare il generatore di numeri casuali della funzione Rnd() assegnandogli un nuovo valore base. Se numero viene omesso, il valore restituito dal timer di sistema verrà utilizzato come nuova base.

Ricordate che la funzione Rnd() restituisce un valore minore di 1 ma maggiore o uguale a

zero. Per generare interi casuali in un dato intervallo, utilizzare la seguente formula:

Int((limitesup - limiteinf + 1) * Rnd + limiteinf)

dove limitesup indica il numero maggiore presente nell'intervallo, mentre limiteinf indica il numero minore.

Generazione di distribuzioni Uniformi Microsoft Excel

L’implementazione della funzione Rnd() dipende dalla particolare versione di Excel In Excel 97 (e fino ad XL 2002) la funzione Rnd() è implementata come

random_number=fractional part of (9821 * r + 0.211327), dove r = .5

per il primo numero

random_number=fractional part of (9821 * r + 0.211327), dove r = precedente numero casuale

per i successivi

(fonte: Microsoft Knowledge Base Articolo 86523)

Microsoft affermava che tale algoritmo è in grado di generare oltre un milione di numeri casuali. Successivamente ha riconosciuto che tale algoritmo non soddisfa i test statistici più sofisticati.

Per questo nella versione XL 2003 tale algoritmo è stato sostituito.


Microsoft Knowledge Base #828795 The RAND function in earlier versions of Excel used a pseudo-random number generation

algorithm whose performance on standard tests of randomness was not sufficient (...) the pseudo-random number generation algorithm that is described here was implemented for Excel 2003. It passes the same battery of standard tests.

The battery of tests is named Diehard (see note 1). The algorithm that is implemented in Excel 2003 was developed by B.A. Wichman and I.D. Hill (see note 2 and note 3). (...) It has been shown by Rotz et al (see note 4) to pass the DIEHARD tests and additional tests developed by the National Institute of Standards and Technology (NIST, formerly National Bureau of Standards).

Notes The tests were developed by Professor George Marsaglia, Department of Statistics, Florida State

University and are available at the following Web site: http://www.csis.hku.hk/~diehard Wichman, B.A. and I.D. Hill, Algorithm AS 183: An Efficient and Portable Pseudo-Random Number

Generator, Applied Statistics, 31, 188-190, 1982. Wichman, B.A. and I.D. Hill, Building a Random-Number Generator, BYTE, pp. 127-128, March

1987. Rotz, W. and E. Falk, D. Wood, and J. Mulrow, A Comparison of Random Number Generators

Used in Business, presented at Joint Statistical Meetings, Atlanta, GA, 2001.


The basic idea is to generate three streams of random numbers (in columns headed "ix", "iy", and "iz") by a common technique and then to use the result that if you take three random numbers on [0,1] and sum them, the fractional part of the sum is itself a random number on [0,1]. The critical statements in the Fortran code listing from the original Wichman and Hill article are:

C IX, IY, IZ SHOULD BE SET TO INTEGER VALUES BETWEEN 1 AND 30000 C BEFORE FIRST ENTRYIX = MOD(171 * IX, 30269)IY = MOD(172 * IY, 30307)IZ = MOD(170 * IZ, 30323)RANDOM = AMOD(FLOAT(IX) / 30269.0 + FLOAT(IY) / 30307.0 + FLOAT(IZ) / 30323.0, 1.0)

Therefore IX, IY, IZ generate integers between 0 and 30268, 0 and 30306, and 0 and 30322 respectively.

These are combined in the last statement to implement the simple principle that was expressed earlier: if you take three random numbers on [0,1] and sum them, the fractional part of the sum is itself a random number on [0,1].

Because RAND produces pseudo-random numbers, if a long sequence of them is produced, eventually the sequence will repeat itself. Combining random numbers as in the Wichman-Hill procedure guarantees that more than 10^13 numbers will be generated before the repetition begins. Several of the Diehard tests produced unsatisfactory results with earlier versions of RAND because the cycle before numbers started repeating was unacceptably short.


Results in Earlier Versions of Excel The RAND function in earlier versions of Excel was fine in practice for users who did not require a lengthy sequence of

random numbers (such as a million). It failed several standard tests of randomness, making its performance an issue when a lengthy sequence of random numbers was needed.

Results in Excel 2003 A simple and effective algorithm has been implemented. The new generator passes all standard tests of randomness

(?????).

The RAND function returns negative numbers in Excel 2003

SYMPTOMS

When you use the RAND function in Microsoft Office Excel 2003, the RAND function may return negative numbers.

CAUSE

This problem may occur when you try to use a large number of random numbers, and you update the RAND function multiple times. For

example, this problem may occur when you update your Excel worksheet by pressing F9 ten times or more.

RESOLUTION

How to obtain the hotfix

This issue is fixed in the Excel 2003 Hotfix Package that is dated January 12, 2004. For additional information, click the following article

number to view the article in the Microsoft Knowledge Base:

833618 Excel 2003 hotfix package released January 12, 2004

STATUS

Microsoft has confirmed that this is a problem in the Microsoft products that are listed in the "Applies to" section of this article.


Generatori di Numeri Pseudocasuali

Ran0, Ran1, Ran2 e Ran3 Press, Teulolsky, Vetterling e Flannery

“Numerical Recipes” Cambridge University Press

Ran0 Generatore “minimale” basato sulla relazione

(Park, S.K. and Miller, K.W., 1988, Comunication of the ACM, Vol. 31, pp. 1192)

214748364712m

168077a

maII

31

5

j1j

mod

214748364712m

168077a

maII

31

5

j1j

mod

Generatori di Numeri Pseudocasuali Run0

vantaggi facilità di implementazione; l’intero codice ammonta ad una

decina di righe di codice C; implementa l’algoritmo di Schrage per il calcolo del modulo

quindi non produce overflow su macchine a 32 bit; periodo sufficientemente elevato per molte applicazioni

svantaggi presenza di correlazione seriale

931 101212 .931 101212 .

Generatori di Numeri Pseudocasuali Run1

Implementa una procedura di shuffling per rimuovere la correlazione seriale di Run0

Generatori di Numeri Pseudocasuali

Generatori di Numeri Pseudocasuali:

Mersenne Twister

Un eccellente generatore di numeri casuali per scopi di simulazione (Makoto Matsumoto e Takuji Nishimura ).

Il Mersenne twister è veloce, ha periodo estremamente lungo (219937-1) , ha eccellenti proprietà dimostrate matematicamente, passa tutti i test conosciuti di casualità (non per scopi

crittografici). L'implementazione C ha velocità paragonabile a quella del

terribile rand(), e qualità neanche lontanamente paragonabile Una versione gratuita per Excel

http://www.numtech.com/NtRand/

Generazione di Numeri Casuali Generazione di Numeri con Funzione di Distribuzione Assegnata

Generazione di distribuzioni non uniformi univariate

Da un generatore di numeri distribuiti uniformemente si possono ricavare numeri distribuiti secondo una densità di probabilità prefissata. I due metodi più importanti sono Transformation Method; Rejection Method.

Metodo della trasformazione inversa

SCOPO: generare un campione di numeri Z distribuiti in accordo ad una funzione di distribuzione assegnata F(z).

INPUT: deve essere possibile valutare la funzione inversa di F(z).

OUTPUT: Z. METODO: Generare un set di numeri casuali U uniformemente

distribuiti fra 0 ed 1 e per ciascuno di questi calcolare Z = F-1(U)

Sia x una variabile aleatoria distribuita uniformemente fra 0 e 1, supponiamo di voler generare una variabile aleatoria con densità di probabilità g(y) essendo y=y(x). Dovremo avere

La soluzione di questa equazione differenziale è

)(

)()(

1 xFy

dzzfyFxy

)(

)()(

1 xFy

dzzfyFxy

dy

dxygdxdyyg )()(

dy

dxygdxdyyg )()(

Metodo della trasformazione inversa

Variabili Normali Univariate

Microsoft Excel

INV.NORM(). Restituisce l'inversa della distribuzione normale cumulativa per la media e la deviazione standard specificate. La sintassi é

INV.NORM(probabilità;media;dev_standard)

dove probabilità è la probabilità corrispondente alla distribuzione normale, media è la

media aritmetica della distribuzione, dev_standard è la deviazione standard della distribuzione.

INV.NORM utilizza una tecnica iterativa per il calcolo della funzione. Dato un valore di probabilità, INV.NORM applica il metodo delle iterazioni fino a quando la precisione del risultato non rientra in ± 3x10^-7. Se il risultato di INV.NORM non converge dopo 100 iterazioni, la funzione restituirà il valore di errore #N/D.

Distribuzioni Parametriche

Uno dei problemi più complessi nell’analisi statistica dei dati è senz’altro quello del del fitting di una distribuzione teorica su un insieme di dati estratti da una distribuzione empirica.

E’ necessario disporre di funzioni flessibili in grado di adattarsi al maggior numero di dati possibile

Es. Distribuzione Lambda Generalizzata

λ

]p)(1[pλ R(p)

2

λλ

1

43

λ

]p)(1[pλ R(p)

2

λλ

1

43

dove 1 è un parametro di localizzazione, 2 un parametro di scala, 3 e 4 sono parametri di forma della distribuzione. Questi parametri sono legati ai parametri di media, deviazione standard, asimmetria e curtosi della distribuzione empirica dei dati di partenza. Al variare di questi parametri si riesce ad ottenere un’enorme varietà di distribuzioni univariate in grado di approssimare una distribuzione empirica di dati.

Distribuzioni Parametriche

La descrizione basata sui percentili risulta particolarmente interessante dal punto di vista della simulazione Monte Carlo perché i percentili di una distribuzione continua sono uniformemente distribuiti fra 0 ed 1.

Questo significa che se indichiamo con u l’argomento della funzione percentile, R, possiamo generare n valori da una distribuzione uniforme in [0, 1] effettuare la trasformazione x = R(u) ed ottenere così una serie di valori per la variabile aleatoria x che risulta essere distribuita in modo tale da avere R come funzione dei percentili.

Generazione di Numeri Casuali Un cenno ai Test Statistici per Generatori di Numeri Casuali

Test sui generatori di numeri casuali Esistono moltissimi test statistici che possono

essere utilizzati per controllare la qualità di un generatore;

http://stat.fsu.edu/~geo/diehard.html Vediamo un semplicissimo esempio di

applicazione del test di 2 sul generatore lineare congruente per due diversi insiemi di parametri rimandando alla letteratura per una descrizione completa delle varie tipologie di test.

2 - test Dato l’algoritmo lineare congruente sono noti diversi

insiemi di parametri che permettono di ottenere buoni generatori di numeri casuali. Ad esempio il generatore definito da

a = 1812433253

c dispari

m = 232

produce un buona distribuzione di punti su macchine a 32 bit. Per eseguire un test su questo generatore possiamo impiegarlo per

produrre un numero elevato di valori (almeno 10000) raccolti in 100 intervalli equispaziati lungo il segmento [0,1]. Il numero attesto di valori all’interno di ciascun intervallo è pari a 100.

2 - test A questo punto calcoliamo il valore di 2 e verifichiamo se il

generatore è accettabile utilizzando i valori riportati in Tabella

La tabella riporta, per ciascun valore dei gradi di libertà k, il valore per cui la probabilità è pari al valore della colonna corrispondente. Ad esempio nella riga corrispondente a 10 gradi di libertà, sotto la colonna p = 95% troviamo il valore 18.3; questo significa che con 10 gradi di libertà dovremmo trovare un valore di 2 superiore a 18.31 in media 5 volte su 100. Quando il numero di gradi di libertà supera 30 possiamo usare una formula asintotica in funzione dei percentili xp riportati per completezza nell’ultima riga della tabella.

N

1i i

2ii2

Np

NpY

N

1i i

2ii2

Np

NpY

2,1 k

2k,1

2P

2 - testValori di probabilità per la distribuzione 2

p = 1% p = 5% p = 25% p = 50% p = 75% p = 95% P = 99%

k = 1 0.00016 0.00393 0.1015 0.4549 1.323 3.841 6.635

k = 2 0.02010 0.1026 0.5753 1.386 2.773 5.991 210

k = 3 0.1148 0.3518 1.213 2.366 4.108 7.815 11.34

k = 4 0.2971 0.7107 1.923 3.357 5.385 488 13.28

k = 5 0.5543 1.1455 2.675 4.351 6.626 11.07 15.09

k = 6 0.8720 1.635 3.455 5.348 7.841 12.59 16.81

k = 7 1.239 2.167 4.255 6.346 037 14.07 18.48

k = 8 1.646 2.733 5.071 7.344 10.22 15.51 20.09

k = 9 2.088 3.325 5.899 8.343 11.39 16.92 21.67

k = 10 2.558 3.940 6.737 342 12.55 18.31 23.21

k = 20 8.260 10.85 15.45 134 23.83 31.41 37.57

k = 30 14.95 18.49 24.48 234 34.80 43.77 50.89

k > 30 3

2

3

22 2 pp xxkk

xp -2.33 -1.64 -0.675 0.00 0.675 1.64 2.33

Valori di probabilità per la distribuzione 2

p = 1% p = 5% p = 25% p = 50% p = 75% p = 95% P = 99%

k = 1 0.00016 0.00393 0.1015 0.4549 1.323 3.841 6.635

k = 2 0.02010 0.1026 0.5753 1.386 2.773 5.991 210

k = 3 0.1148 0.3518 1.213 2.366 4.108 7.815 11.34

k = 4 0.2971 0.7107 1.923 3.357 5.385 488 13.28

k = 5 0.5543 1.1455 2.675 4.351 6.626 11.07 15.09

k = 6 0.8720 1.635 3.455 5.348 7.841 12.59 16.81

k = 7 1.239 2.167 4.255 6.346 037 14.07 18.48

k = 8 1.646 2.733 5.071 7.344 10.22 15.51 20.09

k = 9 2.088 3.325 5.899 8.343 11.39 16.92 21.67

k = 10 2.558 3.940 6.737 342 12.55 18.31 23.21

k = 20 8.260 10.85 15.45 134 23.83 31.41 37.57

k = 30 14.95 18.49 24.48 234 34.80 43.77 50.89

k > 30 3

2

3

22 2 pp xxkk

xp -2.33 -1.64 -0.675 0.00 0.675 1.64 2.33

I criteri di accettazione sono P(2) < 1% o P(2) > 99% rigetto 1% < P(2) < 5% o 95% < P(2) < 99% sospetto 5% < P(2) < 10% o 90% < P(2) < 95% leggero sospetto

Il criterio viene applicato ad almeno tre serie prodotte dal generatore e il risultato è definitivamente negativo se due serie su tre danno risultato negativo.

2 - test



Costruiamo un Add-In di Generatori per Excel Costruiamo un Add-In di Generatori per Excel

Introduzione







Le sequenze a bassa discrepanza L’idea che sta dietro a questo approccio è abbastanza intuitiva e può

essere facilmente spiegata almeno nel caso unidimensionale.

Supponiamo di voler integrare una funzione f(x) nell’intervallo [0, 1]

utilizzando una sequenza di n punti. Invece di scegliere una serie di

punti disposti a caso sul segmento unitario possiamo selezionare una

sequenza di punti distribuiti in maniera tale da garantire una certa

uniformità.

Possiamo vedere che con questa scelta l’accuratezza nella stima

dell’integrale è maggiore di quella che otterremo scegliendo i punti a

caso secondo l’approccio classico del metodo Monte Carlo. In particolare se scegliessimo una serie di punti posti all’interno di una

griglia equispaziata riotteremo il classico metodo numerico di integrazione

tramite approssimazione trapezoidale il quale ha un errore che va come

l’inverso del quadrato del numero dei punti

Le sequenze a bassa discrepanza La domanda che sorge spontanea a questo punto è perché non usare

sempre sequenze di punti equispaziati. La risposta va ricercata nel comportamento dell’errore dovuto alla

discretizzazione quando affrontiamo l’integrazione di funzioni di più variabili. Come sappiamo la convergenza del metodo Monte Carlo non dipende dal

numero di dimensioni del dominio di integrazione. I metodi numerici di integrazione come il metodo trapezzoidale, al contrario,

presentano la caratteristica di avere un errore che va tipicamente come n-2/d per cui è sufficiente che la dimensione sia superiore a 4 per dare un risultato peggiore di quello del Monte Carlo classico.

Inoltre l’utilizzo di una griglia equispaziata di punti presenta l’inconveniente che il numero di punti deve essere deciso in anticipo e questo può essere in molti casi una inaccettabile restrizione in quanto è più comodo disporre di procedure numeriche in grado di arrestarsi quando la precisione di calcolo giunge nell’intorno del valore desiderato.

Le sequenze a bassa discrepanza La caratteristica più importante delle sequenze deterministiche

che andremo ad introdurre è legata al fatto che i punti che via

via vengono introdotti tendono a disporsi il più possibile

all’interno degli spazi lasciati vuoti dagli altri punti.

In questo modo successive generazioni di punti tenderanno a

riempire nella maniera più uniforme possibile il dominio di

spazio in esame.

Per definire in maniera precisa una misura del “riempimento

uniforme” di un volume di spazio occorre introdurre il concetto

di DISCREPANZA.

Le sequenze a bassa discrepanza Generiamo N vettori {ri} in uno spazio a d dimensioni;

Ciascuno di questi vettori rappresenta un punto in uno spazio a

dimensioni le cui coordinate sono le singole componenti del vettore,

quindi ogni vettore individua un punto nell’ipercubo [0,1]d

Selezioniamo un sotto-ipercubo S(y) scegliendo un punto y che

individua lo spigolo superiore destro di un dominio iper-cubico da 0 ad

y;

In altre parole possiamo scrivere

),[),[)( d1 y0y0yS ),[),[)( d1 y0y0yS

Le sequenze a bassa discrepanza Indichiamo con nS(y) il numero di punti che si trovano in S(y)

Nel limite in cui N tende all’infinito richiedere la perfetta omogenità del

generatore equivale a richiedere che

per tutti gli y in [0,1]d

N

1i

d

1kry

N

1iySryS ikki

n 11 )()(

N

1i

d

1kry

N

1iySryS ikki

n 11 )()(

d

1ii

yS

Ny

N

n )(lim

d

1ii

yS

Ny

N

n )(lim

Le sequenze a bassa discrepanza Infatti questo equivale ad affermare che per una distribuzione

perfettamente uniforme su un ipercubo, la probabilità di estrarre un

punto appartenente ad una sottoregione dell’ipercubo è pari al volume

di tale regione;

Per avere un’idea della discrepanza complessiva nella distribuzione

dei punti possiamo quindi comparare le due quantità

per ogni y

))(( e )( ySVN

n yS ))(( e )( ySV

N

n yS

Le sequenze a bassa discrepanza Solitamente si sceglie la norma L il che porta alla seguente definizione di

discrepanza

Nel caso unidimensionale è facile calcolare una forma esplicita per la

discrepanza (star). Etichettiamo i punti della sequenza in modo tale da avere

la discrepanza della sequenza è quindi data da

d

1kk

yS

10y

dN y

N

nD

d

)(

],[

)( sup

d

1kk

yS

10y

dN y

N

nD

d

)(

],[

)( sup

10 21 nxxx

n2

1k2x

n2

1D k

n1k

1n

,

)( maxn2

1k2x

n2

1D k

n1k

1n

,

)( max

Le sequenze a bassa discrepanza nel caso di dimensioni maggiori non esiste in generale una formula semplice per il calcolo

della discrepanza Possiamo a questo punto dare una definizione formale di sequenza a bassa discrepanza:

Diremo che una sequenza in [0,1]d è a bassa discrepanza se per tutti gli N > 1 i primi N punti della sequenza soddisfano

dove c(d) è una costante che dipende solo da d. N

NdcD

dd

N

ln)()(

N

NdcD

dd

N

ln)()(

La sequenza di Halton La successione di Halton è ottenibile con la seguente

procedura Al passo j scrivere j stesso in base b con b numero primo

Il numero casuale generato al passo j viene calcolato come

nnbdbdbdj 1

10

0

nnbdbdbdj 1

10

0

10

21

1

nnnj bdbdbdH

10

21

1

nnnj bdbdbdH

La sequenza di Halton In sostanza la procedura consiste nel capovolgere le cifre

dell’espansione in base b e porre un punto decimale di fronte a quanto ottenuto.

Se sono necessari numeri casuali multidimensionali si può procedere generando per ciascuna componente una sequenza di Halton con radice b diversa. Solitamente si utilizzano in successione i numeri primi (2, 3, 5, … )

Es. consideriamo b = 2

1250212020H2020214

7502121H21213

2502120H20212

5021H211

3214

012

213

01

212

01

11

0

.

.

.

.

1250212020H2020214

7502121H21213

2502120H20212

5021H211

3214

012

213

01

212

01

11

0

.

.

.

.

La sequenza di Halton I numeri via via generati tendono a riempire i buchi

lasciati dai precedenti numeri; Quando vogliamo generare vettori in spazi a più

dimensioni, per evitare sovrapposizioni nel riempimento dello spazio, è sufficiente scegliere come base numeri incommensurabili;

Per questo solitamente per ogni dimensione si sceglie un diverso numero primo.

La seguenza di Halton

Public Function Halton(n As Integer, _ x As Integer) As Double ' n è il numero da trasformare ' x è la base ' H è il numero generato dall'algoritmo Dim H As Double Dim z As Double Dim m As Integer Dim na As Integer Dim nb As Integer H = 0 na = n z = 1 / x While (na > 0) nb = Int(na / x) m = na - nb * x H = H + m * z na = nb z = z / x Wend Halton = H End Function

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.2 0.4 0.6 0.8 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.2 0.4 0.6 0.8 1

Comportamento ad Alte Dimensionalità

La figura illustra un problema tipico delle sequenze deterministiche.

Quando si eseguono simulazioni su un numero elevato di dimensioni si scopre che le dimensioni più alte tendono ad essere estremamente correlate.

Questo problema non affligge in maniera uguale tutti gli algoritmi, in particolare le sequenze generate con l’algoritmo di Sobol tendono ad essere più robuste, tuttavia in tutti i casi si nota un progressivo aumento della discrepanza all’aumentare delle dimensioni.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Esempio Excel Esempio Excel

Sequenze a Bassa DiscrepanzaSequenze a Bassa Discrepanza

Introduzione







Misure di co-dipendenza

Distribuzioni Marginali Data la distribuzione congiunta di due variabili x ed y la funzione di densità

marginale di x è definita come

E, analogamente,

)(

),()(yD

x dyyxx )(

),()(yD

x dyyxx

)(

),()(xD

y dxyxy )(

),()(xD

y dxyxy


Indipendenza Due variabili x ed y si dicono indipendenti se la loro

funzione densità congiunta si fattorizza nel prodotto delle densità marginali

)()(),( , yxyxtiindipendenyx yx )()(),( , yxyxtiindipendenyx yx

Correlazione Lineare Ricordiamo la definizione di correlazione lineare tra

due variabili x ed y


2222

,

)()()()(

)()(),(

)()(

),cov(

dyyydyyydxxxdxxx

dyyydxxxdxdyyxxy

yx

yx

yyxx

yx

yx

2222

,

)()()()(

)()(),(

)()(

),cov(

dyyydyyydxxxdxxx

dyyydxxxdxdyyxxy

yx

yx

yyxx

yx

yx

Correlazione Lineare Spesso si ritiene che la conoscenza del coefficiente di correlazione

lineare unitamente alla specificazione delle distribuzioni marginali, permetta di determinare la distribuzione di probabilità congiunta.

In realtà questo è vero solo per certe classi di distribuzioni tra cui la distribuzione normale.

In generale quindi l’inferenza

Non è valida


),(),(),( , yxyx yxyx ),(),(),( , yxyx yxyx

Oltre l’indice di correlazione lineare La correlazione lineare è un buona misura di co-dipendenza per

variabili normali. Per distribuzioni che non sia allontanano troppo dalla “normalità” continua a fornire indicazioni utili ma all’allontanarsi da queste condizioni (in molti casi soltanto ideali) l’indice di correlazione lineare fornisce risultati sempre più fuorvianti!

L’indice di correlazione lineare non è invariante rispetto a trasformazioni non lineari delle variabili.

A differenza degli stimatori non-parametrici, l’indice di correlazione lineare può non coprire l’intero range da – 1 a + 1, rendendo problematica l’interpretazione del grado di dipendenza


Rho di Sperman Questo indice è definito come coefficiente di correlazione lineare fra le funzioni

di distribuzione delle due variabili. In altre parole, date due variabili x ed y e le loro distribuzioni marginali, calcoliamo prima di tutto le distribuzioni marginali cumulate

Utilizzando due semplici risultati…


y

yy

x

xx dydx )(:)( )(:)(

y

yy

x

xx dydx )(:)( )(:)(

1

0

221

0 3

1)()( ,

2

1)()( duudxxxududxxx xxx x

1

0

221

0 3

1)()( ,

2

1)()( duudxxxududxxx xxx x

Rho di Sperman … possiamo scrivere

Tau di Kendal


3),()()(12 dxdyyxyx yxS 3),()()(12 dxdyyxyx yxS

1),(),(4 dxdyyxyxK 1),(),(4 dxdyyxyxK


Il tau di Kendal e il Rho di Sperman appartengono alla categoria delle cosiddette rank correlation ;

Sono invarianti per un’ampia tipologia di trasformazione delle variabili;

Questo tipo di indicatori, a differenza del coefficiente di correlazione lineare, ha la proprietà che date due distribuzioni marginali esiste sempre una distribuzione congiunta per ogni valore dell’indice compreso nell’intervallo [-1, 1].

Generazione di Numeri Casuali Variabili Normali Multivariate

Variabili Normali Multivariate Cholescky Decomposition

Indichiamo con X un vettore di variabili aleatorie indipendenti ciascuna delle quali distribuita secondo una normale standard, la matrice di varianza-covarianza di X sarà pertanto data dalla matrice unità di dimensione n n. Supponiamo di voler derivare da questo insieme di variabili un secondo set di variabili, che indicheremo con Y, non più indipendenti bensì dotato di matrice di varianza-covarianza assegnata .

Il nuovo insieme di variabili aleatorie può essere ricercato come combinazione lineare delle variabili indipendenti , cioè si pone

Il problema si riconduce così alla determinazione di una matrice A di dimensione n n tale che

AXY AXY

tAA tAA

Cholescky Decomposition

La soluzione della precedente equazione non è unica nel senso che esistono più matrici A che, moltiplicate per la loro trasposta, danno come risultato . Se la matrice è definita positiva il metodo più efficiente dal punto di vista computazionale per risolvere il problema consiste nell’applicazione della scomposizione di Cholescky.

Il punto chiave di tale metodologia consiste nel ricercare A nella forma di una matrice triangolare inferiore, ovvero una matrice in cui tutti gli elementi sopra la diagonale sono nulli,

nnnn AAA

AA

A

A

21

2221

11

0

00

nnnn AAA

AA

A

A

21

2221

11

0

00

Variabili Normali Multivariate

Cholescky Decomposition

Sviluppando il prodotto AAt in componenti è facile verificare che gli elementi di A sono ricavabili dalle seguenti formule iterative

Ad esempio per il caso semplice di due variabili troviamo

1

1

2i

kikiiii aa

1

1

2i

kikiiii aa

1

1

1 i

kjkikij

iiji aa

aa

1

1

1 i

kjkikij

iiji aa

aa

2

22

1

1

0

A

2

22

1

1

0

A

Variabili Normali Multivariate

Il Condizionamento della Matrice di Varianza e Covarianza

Può accadere, specialmente quando si lavora con problemi ad elevata dimensionalità, che la matrice di correlazione non risulti semi-definita positiva;

In tal caso almeno uno degli autovalori della matrice risulterà negativo. Spesso questo fatto è riconducibile alle procedure che hanno portato alla

costruzione della matrice stessa; ad esempio è sufficiente avere serie storiche di prezzi non perfettamente allineate o prezzi

registrati a tempi diversi per produrre quasi sicuramente una matrice di correlazione mal definita.

In questi casi non è sempre possibile ricostruire la matrice stessa (si pensi al caso limite in cui la matrice viene scaricata da un provider esterno come nel caso di RiskMetrics) per cui è necessario disporre di metodi che possano rimuovere gli autovalori negativi in modo che la matrice risulti definita positiva col minor impatto possibile sui valori della matrice stessa.

Una procedura possibile è la seguente

Calcolare autovalori e autovettori della matrice di correlazione Porre gli autovalori negativi uguali a zero; Ricostruire la nuova matrice di correlazione;

L’ultimo passaggio è facilmente realizzabile a partire dalla matrice degli autovettori E e dalla matrice (diagonale) degli autovalori

Il Condizionamento della Matrice di Varianza e Covarianza

tEEC tEEC E’ molto probabile che la nuova matrice così ottenuta abbia elementi lungo la diagonale diversi da 1. Per E’ molto probabile che la nuova matrice così ottenuta abbia elementi lungo la diagonale diversi da 1. Per questo è necessario procedere ad una normalizzazione ponendo questo è necessario procedere ad una normalizzazione ponendo

dove Ddove D1/21/2 è una matrice diagonale i cui elementi sono le radici quadrate degli elementi diagonali di è una matrice diagonale i cui elementi sono le radici quadrate degli elementi diagonali di C.C.

DC

DC

11 D

CD

C11

Modello a due fattori di Fong e Vasicek

Il modello risulta particolarmente interessante in quanto: La distribuzione di probabilità del tasso di interesse non è normale; La volatilità è stocastica; La volatilità è correlata col livello del tasso di interesse (effetto GARCH)

Il modello non ha soluzione in forma chiusa;

Z1 e Z2 sono due processi di Wiener che presentano una correlazione

istantanea pari a ;

1

1

)(

)(

dzdtd

dzdtrdr

1

1

)(

)(

dzdtd

dzdtrdr

Al fine di simulare il processo occorre prima di tutto trasformare la rappresentazione nel tempo continuo fornita dalle equazioni differenziali nella sua controparte nel tempo discreto.

Applicando lo schema di Eulero otteniamo:

tt

ttrrr

tttt

tttt

2111

1111

tt

ttrrr

tttt

tttt

2111

1111


Nella versione discretizzata del modello 1 ed 2 sono due numeri casuali generati a partire da una distribuzione normale bivariata.

La decomposizione di Cholesky nel caso semplice di due sole variabili si riduce alle seguenti formule

22

12

11

1

22

12

11

1


L’esempio numerico 5 periodi temporali ciascuno pari ad un anno; due ZCB che producono un cash flow di 1$ al

periodo 3 e al periodo 4; un CAP che produce un flusso pari a max(100x(r5-

0.06),0), dove r5 è il tasso a breve al periodo 5.




Utilizzo Variabili Normali Multivariate Utilizzo Variabili Normali Multivariate

Generazione di Numeri Casuali Utilizzo delle funzioni di Copula

Funzioni di copula: concetti base Una funzione z = C(u,v) è detta copula se e solo se

1. z, u e v sono in [0,1]

2. C(0,v) = C(u,0) = 0, C(1,v) = v, C(u,1) = u

3. C(u2, v2 ) – C(u1, v2 ) – C (u2, v1) – C (u1, v1) 0 per tutti i valori

u2 > u1 e v2 > v1

Teorema di Sklar Ogni distribuzione congiunta può essere scritta come una funzione di

copula che abbia le distribuzioni marginali come argomenti e qualsiasi funzione di copula che abbia distribuzioni come argomenti è una distribuzione congiunta

Possiamo quindi scrivere

L’utilizzo delle funzioni di copula consente di specificare separatamente le distribuzioni marginali delle variabili e la loro struttura di dipendenza

Le funzioni di copula sono legate alle statistiche non-parametriche di dipendenza, ad esempio il di Kendall o il di Spearman.

Funzioni di copula: concetti base

)(),(),( yxCyx yx )(),(),( yxCyx yx

Algoritmi per la modellazione della dipendenza

La copula gaussiana

Calcolo della pseudo-radice di RA: R=AAT

INPUT

n funzioni di densità marginali

matrice di correlazione R

)(,),(),()( 2211 nxnxxx xxxx

Generazione di n variabili normalistandard non correlate

z

Da z si ricava un vettore di variabili vuniformemente distribuite fra 0 ed 1

ponendo

vi = N(z’i)

Introduzione della correlazionez’ = Az

Calcolo del vettore x utilizzando ladistribuzione cumulativa inversa

marginale

)(11 ixi vx

END

START

L’algoritmo si riassume nei seguenti passi generare un vettore di variabili normali con

correlazione assegnata; Trasformare tali variabili in variabili con distribuzione

uniforme in [0, 1] utilizzando la funzione comulata della normale;

Utilizzare queste variabili come base per generare il vettore x, secondo il metodo della trasformazione inversa utilizzando come funzioni di trasformazione le inverse delle comulate marginali specificate in input.

E’ importante ricordare che i coefficienti di correlazione che controllano la copula gaussiana, e che vengono utilizzati nell’algoritmo descritto a lato, possono essere molto diversi dalla correlazione lineare

Algoritmi per la modellazione della dipendenza

La copula t

Due variabili connesse da una copula di tipo gaussiano o di tipo t e con coefficiente di correlazione hanno un tau di Kendal pari dato da

arcsin2

K

arcsin

2K

Calcolo della pseudo-radice di RA: R=AAT

INPUT

n funzioni di densità marginali matrice di correlazione R

Generazione di n variabili normalistandard non correlate

z

Introduzione della correlazionez’ = Az

START

Generare una variabile, s, distribuitasecondo una chi-quadro con n gradidi libertà. Per n intero questo puòessere ottenuto semplicementegenerando n variabili normaliindipendenti e sommando i loroquadrati.

Da y si ricava un vettore di variabili vuniformemente distribuite fra 0 ed 1utilizzando la funzione di probabilità

comulativa della t di Student

Porre

zs

ny

Calcolo del vettore x utilizzando ladistribuzione cumulativa inversa

marginale

)(11 ixi vx

END



Esempio Copule Esempio Copule

Introduzione







Un Problema di Efficienza Immaginiamo di voler calcolare un certo parametro P (ad esempio il

prezzo di un’opzione) e di poter scegliere fra due diverse stime ottenibili con il metodo Monte Carlo rappresentate dalle due serie di valori ottenuti con il processo di simulazione

Supponiamo poi che entrambi gli stimatori siano corretti, cioè valga

ma con

niP i ,...,1,1̂ niP i ,...,1,1̂ niP i ,...,1,2̂ niP i ,...,1,2̂

PPE 1̂ PPE 1̂ PPE 2̂

PPE 2̂

21 21

Un Problema di Efficienza Chiaramente sulla base di queste sole informazioni

saremmo portati a scegliere il primo stimatore in quanto, a parità di numero di simulazioni, l’errore di stima risulterà senz’altro minore.

Tuttavia, come accennavamo poco sopra, questa conclusione rischia in realtà di non essere corretta in quanto non tiene conto del fatto che i due stimatori possono richiedere risorse computazionali molto diverse fra loro;

in particolare generare n replicazioni di P1 potrebbe richiedere molto più tempo che generare n replicazioni di P2.

Un Problema di Efficienza Un primo approccio al problema potrebbe essere quello di

introdurre esplicitamente nelle nostre considerazioni il tempo di calcolo richiesto.

Supponiamo che il tempo richiesto per generare una singola replicazione di Pj possa essere espresso da una costante che indicheremo con bj, avendo a disposizione un tempo totale di calcolo pari a t il numero di replicazioni di Pj che possiamo generare sarà pari a t/bj.

I due stimatori possono pertanto essere riscritti introducendo esplicitamente il tempo di calcolo nelle formule

1/

1

11 ˆbt

iiP

t

b

1/

1

11 ˆbt

iiP

t

b

2/

1

22 ˆbt

iiP

t

b

2/

1

22 ˆbt

iiP

t

b

Un Problema di Efficienza Per valori sufficientemente alti di t queste due quantità sono normalmente

distribuite con media P e standard deviation

pertanto per grandi valori di t il primo stimatore sarà preferibile al secondo se

l’inverso del prodotto della varianza per il tempo necessario ad eseguire un singolo run viene indicato in letteratura col nome di efficienza (Hammersley e Handscomb, 1964).

Usando l’efficienza come base per il confronto fra diversi stimatori possiamo concludere che lo stimatore a bassa varianza è preferibile all’altro solo se il rapporto delle varianza è più piccolo del rapporto fra i tempi di singola

replicazione.

t

b

t

b 22

11 t

b

t

b 22

11

2221

21 bb 2

221

21 bb

Variabili Antitetiche Uno dei metodi più semplici e più utilizzati in campo finanziario per

la riduzione della varianza è senz’altro il metodo delle variabili antitetiche.

Consideriamo di nuovo la procedura classica di stima del prezzo di un’opzione, per semplicità espositiva ci limiteremo ancora al contesto del modello di Black e Scholes.

Se si adotta la tecnica della variabile antitetica, in ogni simulazione si devono determinare due valori.

Il primo valore è quello calcolato nel modo consueto...

...mentre il secondo valore viene calcolato cambiando segno a tutti i campioni estratti casualmente dalle distribuzioni normali standardizzate.

Variabili Antitetiche I due stimatori hanno chiaramente le stesse proprietà

statistiche essendo estratti dallo stesso campione (in particolare hanno la stessa varianza).

Il valore campionario del derivato calcolato in ogni simulazione è la media di questi due valori e la sua varianza è data da

iii

iiii

CCCovCVar2

1

CCVar4

1

2

CCVar

~,

~~

iii

iiii

CCCovCVar2

1

CCVar4

1

2

CCVar

~,

~~

Variabili Antitetiche

Pertanto se iii CVarCCCov ~, avremo imedio CVarCVar .

Comunque occorre tenere presente che questa procedura richiede un numero di

simulazioni doppio rispetto al caso standard per cui è necessario ragionare in termini di efficienza.

Se supponiamo che la generazione dei numeri casuali richieda un tempo trascurabile

rispetto al calcolo del prezzo allora possiamo affermare che il tempo impiegato utilizzando le variabili antitetiche sia all’incirca doppio di quello richiesto nel caso standard.

In questo caso il metodo delle variabili antitetiche è preferibile in termini di efficienza

rispetto al metodo standard se si verifica la condizione

imedio CVarCVar 2


Questa condizione è equivalente a richiedere che

0~

, ii CCCov

Verifichiamo se questa condizione è valida.

Indichiamo con la funzione definita dalla relazione )( ii ZC ; supponiamo che sia

la composizione di due funzioni monotone, la prima è quella che lega il valore del sottostante alla variabile aleatoria Z, la seconda è la funzione che calcola il payoff come valore massimo fra 0 e la differenza fra il prezzo del sottostante e lo strike price. In queste condizioni anche è monotona.


Utilizzando una disuguaglianza standard possiamo allora verificare che

)()()()( iiii ZEZEZZE

da cui segue immediatamente

0)()()()(~

, iiiiii ZEZEZZECCCov

Quindi il metodo delle variabili antitetiche è più efficiente del metodo standard a patto

che siano verificate le condizioni di monotonicità citate.

Nel caso di payoff non monotoni il metodo non necessariamente fornisce prestazioni migliori del Monte Carlo standard, anzi, in alcune condizioni i risultati sono sensibilmente peggiori.

Moment Matching Il metodo dei momenti (moment matching) comporta l’aggiustamento dei campioni estratti da una distribuzione normale standardizzata in modo da assicurare l’uguaglianza tra i momenti campionari (in genere il primo e il secondo) e i corrispondenti momenti della distribuzione probabilistica.

Indichiamo con Zi i campioni estratti da una distribuzione normale usati per calcolare la variazione di valore di una certa variabile in un certo periodo di tempo. Per assicurare l’uguaglianza dei primi due momenti calcoliamo la media campionaria m e la deviazione standard campionaria s. Quindi definiamo nel modo seguente i campioni aggiustati

s

mZZ i

i

'

s

mZZ i

i

'


Confronto fra alcuni metodi di Riduzione della VarianzaConfronto fra alcuni metodi di Riduzione della Varianza

Greek Letters: Differenze Finite

Consideriamo il problema legato al calcolo del Delta di un’opzione europea

Un approccio diretto al problema può consistere nella generazione di due prezzi finali, il primo

a partire dal valore S0 , e il secondo a partire dal valore perturbato S0 +

0S

C

0S

C

ZTTrT eSS )2/(

02 ZTTr

T eSS )2/(0

2

')2/(0

2

)( ZTTrT eSS

')2/(0

2

)( ZTTrT eSS

Z e Z’ sono due estrazioni indipendenti da una normale standard

Greek Letters

Per ogni valore del prezzo finale possiamo poi calcolare il valore dell’opzione corrispondente

),0max(0 KSeSC TrT ),0max(0 KSeSC T

rT

))(,0max(0 KSeSC TrT ))(,0max(0 KSeSC T

rT

00~ SCSC

00~ SCSC

Greek Letters

Problema

Poiché i valori per ST e ST() sono generati indipendentemente l’uno dall’altro la varianza di delta diverge al diminuire del valore di .

Per ottenere uno stimatore che converga verso il valore effettivo del Delta occorre diminuire in maniera graduale (e lenta) il valore di all’aumentare di n. Complessivamente questo rallenta la velocità di convergenza fino a livelli del tutto inaccettabili.

)()()(~ 2

002 OSCVarSCVarVar )()()(

~ 200

2 OSCVarSCVarVar

Greek Letters Una stima migliore può essere ottenuta utilizzando il metodo dei Numeri

Casuali Comuni (Common Random Numbers) che nella fattispecie si traduce nell’utilizzare lo stesso numero casuale Z sia nel calcolo di S0 che di S0 + .

Se denotiamo con ^ la stima del Delta così calcolata; per un valore di fissato, la media di un campione di repliche indipendenti di ^ converge al valore effettivo di Delta ma il calcolo della varianza ora fornisce un valore diverso in quanto C(S0) e C(S0 + ) non sono più indipendenti

)1(00 O

SCSCVar

)1(00 OSCSC

Var

Cov > 0

)(),(2)()(ˆ0000

2 SCSCCovSCVarSCVarVar )(),(2)()(ˆ0000

2 SCSCCovSCVarSCVarVar


Common Random Numbers Stima del DeltaCommon Random Numbers Stima del Delta

Variabile di controllo Il metodo della variabile di controllo mantiene inalterata la funzione

di distribuzione campionata, il miglioramento dell’efficienza si ottiene in questo caso ricorrendo ad una funzione ausiliaria v(z) correlata a k(z) di cui è noto esattamente l’integrale

Il metodo della variabile di controllo funziona bene quando quest’ultima ha un elevato grado di correlazione con la variabile che intendiamo stimare. Una situazione di questo tipo si presenta quando vogliamo conoscere il prezzo di un’opzione asiatica se utilizziamo come variabile di controllo il prezzo dell’opzione asiatica a media geometrica corrispondente.

Z

C zdFzv )()(Z

C zdFzv )()(

Variabile di controllo

Definiamo W(, z) come

CzvzkzW )()(),(

E’ facile verificare che )(),( zkEzWE

222 2, v(z)k(z),v(z)k(z) σββσσzWVar

Possiamo scegliere il parametro in modo da minimizzare la varianza della funzione W

2)(

)(),(*

zv

zvzk

nel qual caso

222 1, k(z)k(z),v(z)k(z) σ)ρ(σzWVar

Definiamo W(, z) come

CzvzkzW )()(),(

E’ facile verificare che )(),( zkEzWE

222 2, v(z)k(z),v(z)k(z) σββσσzWVar

Possiamo scegliere il parametro in modo da minimizzare la varianza della funzione W

2)(

)(),(*

zv

zvzk

nel qual caso

222 1, k(z)k(z),v(z)k(z) σ)ρ(σzWVar

Variabile di controllo Un esempio pratico: il pricing di un’opzione asiatica

Le opzioni asiatiche sono opzioni il cui valore finale dipende dal prezzo medio dell’attività sottostante osservato, almeno in parte, durante la vita dell’opzione.

Il valore finale di una call scritta sul prezzo medio (average price call) è max(0,Save-X) e quello di una put scritta sul prezzo medio

(average price put) è pari a max(0, X – Save) essendo il prezzo

medio calcolato in un periodo determinato (Hull, 2000). Le opzioni average price sono meno care delle opzioni ordinarie in

quanto il calcolo della media diminuisce di fatto la volatilità del sottostante.

Variabile di Controllo Se si assume che il prezzo dell’attività sottostante, S, sia distribuito

in modo log-normale e che Save sia una media geometrica degli S, possiamo utilizzare delle formule analitiche per valutare le opzioni asiatiche di tipo europeo.

Ciò dipende dal fatto che la media geometrica di un insieme di variabili distribuite in modo log-normale è anch’essa log-normale.

Si può dimostrare che in un mondo neutrale verso il rischio il prezzo di un’opzione asiatica scritta su una media geometrica calcolata su m periodi temporalmente equidistanziati...

mm

jt j

SG

/1

1

mm

jt j

SG

/1

1

)()(

2

1exp)exp( 21

2 dKNdNrTC GGG

)()(

2

1exp)exp( 21

2 dKNdNrTC GGG

è pari a ...

Variabile di Controllo dove

mThhT

qrSG / 22

1ln 2

0

m

mmhG 6

)1)(12(22

G

GG Kd

2

1

)ln(

Gdd 12

Variabile di Controllo Nella simulazione Monte Carlo standard il prezzo dell’opzione viene calcolato come

niKAEeC irTi ,,1 ,)0,max( )()( dove A(i) è la media aritmetica discreta campionata

m

j

tii

jSm

A1

)()( 1

calcolata su un insieme discreto di punti

mjtm

Thhtt jj ,,2,1 ,0 , , 01

ed n è il numero di simulazioni. Questo porta allo stimatore

n

i

iCn

C1

)(1ˆ

Nella simulazione Monte Carlo standard il prezzo dell’opzione viene calcolato come

niKAEeC irTi ,,1 ,)0,max( )()( dove A(i) è la media aritmetica discreta campionata

m

j

tii

jSm

A1

)()( 1

calcolata su un insieme discreto di punti

mjtm

Thhtt jj ,,2,1 ,0 , , 01

ed n è il numero di simulazioni. Questo porta allo stimatore

n

i

iCn

C1

)(1ˆ


Utilizzando il metodo della variabile di controllo oltre alle variabili descritte sopra dobbiamo calcolare la media geometrica per ogni simulazione

mm

j

it

i

jSG

/1

1

)()(

e il valore campionato dell’opzione asiatica a media geometrica

)0,max()exp( )()( KGrTC iiG

Questa volta però utilizzeremo lo stimatore

n

iG

iG

i CCCn

C1

)()( )(1ˆ

Utilizzando il metodo della variabile di controllo oltre alle variabili descritte sopra dobbiamo calcolare la media geometrica per ogni simulazione

mm

j

it

i

jSG

/1

1

)()(

e il valore campionato dell’opzione asiatica a media geometrica

)0,max()exp( )()( KGrTC iiG

Questa volta però utilizzeremo lo stimatore

n

iG

iG

i CCCn

C1

)()( )(1ˆ


0,00

5,00

10,00

15,00

20,00

25,00

0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 0,45 0,50

Volatilità Sottostante

Asiatica Aritmetica

Europea Black & Scholes

Asiatica Geometrica

0,00

5,00

10,00

15,00

20,00

25,00

0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 0,45 0,50


Asiatica Aritmetica


Asiatica Geometrica


0,00

5,00

10,00

15,00

20,00

25,00

30,00

35,00

0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0

Tempo a Scadenza

Asiatica Aritmetica


Asiatica Geometrica

0,00

5,00

10,00

15,00

20,00

25,00

30,00

35,00

0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0

Tempo a Scadenza

Asiatica Aritmetica


Asiatica Geometrica


0,0000

0,1000

0,2000

0,3000

0,4000

0,5000

0,6000

0,7000

0,8000

0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 0,45 0,50


MC Standard MC Controllo MC Antithetic

0,0000

0,1000

0,2000

0,3000

0,4000

0,5000

0,6000

0,7000

0,8000

0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 0,45 0,50


MC Standard MC Controllo MC Antithetic


0,0000

0,1000

0,2000

0,3000

0,4000

0,5000

0,6000

0,7000

0,8000

0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0

Tempo a Scadenza

MC Standard MC Controllo

MC Antithetic

0,0000

0,1000

0,2000

0,3000

0,4000

0,5000

0,6000

0,7000

0,8000

0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0

Tempo a Scadenza

MC Standard MC Controllo

MC Antithetic


Pricing di un’Opzione AsiaticaPricing di un’Opzione Asiatica

Campionamento per Importanza (Importance Sampling)

Consideriamo due diversi piani di campionamento

ZzzFzk

),(),( e ZzzFzk

),(),( **

Supponiamo inoltre che )(zF

e )(* zF

siano entrambe differenziabili in Z e

indichiamo con f(z) e f*(z) le rispettive densità di probabilità.

Definiamo

Zz )()()(*

)()()(* zlzk

zf

zfzkzk

l(z) viene chiamata importance function.

Possiamo calcolare valore atteso e varianza del nuovo stimatore


)(** zkEzkE FF

Z

F dzzfzkzk 22 )()()(var

Z

F dzzfzlzkzk 22* )(*)()()(*var

Z

FF dzzfzlzkzkzk 22* )()](1)[()(*var)(var

Z

FF dzzfzlzkzkzk 22* )()](1)[()(*var)(var


Da queste formule risulta che, almeno in linea di principio, è possibile ridurre a zero la

varianza campionata da )(* zF ; infatti è sufficiente scegliere

)()()(*

zfzkzf

per ottenere dalla 0)(*var * zkF .

La sostituzione effettuata è comunque di interesse puramente teorico in quanto richiede

la conoscenza di che è proprio la quantità che intendiamo stimare.

In generale occorre scegliere una funzione f*(z) tale che

0)()](1)[( 22 Z

dzzfzlzk

Questo metodo risulta particolarmente efficace ogniqualvolta k(z) è una funzione non

limitata nel dominio di integrazione. In questi casi scegliere una f*(z) tale che k*(z) sia limitato può portare ad una considerevole riduzione della varianza.

Campionamento stratificato La campionatura stratificata (stratified

sampling) comporta la suddivisione in strati, o intervalli, della distribuzione probabilistica sottostante e l’estrazione di campioni da ciascun intervallo in base alla probabilità che è ad esso associata.

Consideriamo, per esempio, la generazione di 100 variabili distribuite normalmente. Il metodo più ovvio di fare ciò consiste nel

generare 100 numeri uniformemente distribuiti fra 0 ed 1 e calcolare per ciascuno di questi la funzione inversa della distribuzione normale.

Questo metodo tuttavia produce un risultato che risulta alquanto scadente poiché le code della distribuzione saranno sicuramente sottocampionate

senza campionamento stratificato

0

2

4

6

8

10

12

14

-3-2

.25

-1.5

-0.7

5 00.

75 1.5

2.25 3

ClasseF

req

ue

nza

senza campionamento stratificato

0

2

4

6

8

10

12

14

-3-2

.25

-1.5

-0.7

5 00.

75 1.5

2.25 3

ClasseF

req

ue

nza

Campionamento stratificato Un metodo alternativo, più efficace,

consiste nel forzare ciascun numero generato al passo i a cadere esattamente fra l’ (i-1)-esimo percentile e l’i-esimo.

Un modo estremamente semplice per ottenere questo risultato consiste nel generare 100 numeri Ui uniformemente distribuiti in [0, 1 ] e calcolare*

* I numeri così generati non sono più indipendenti e questo complica la stima dell’errore e del livello di confidenza. Questo problema è comune a tutte le tecniche di riduzione della varianza e non viene discusso in questa introduzione. Il lettore comunque deve essere consapevole della sua esistenza.

100

1~ 1 ii

UiNU

100

1~ 1 ii

UiNU

con campionamento stratificato

0

2

4

6

8

10

12

-3-2

.25

-1.5

-0.7

5 00.

75 1.5

2.25 3

Classe

Fre

qu

en

za

con campionamento stratificato

0

2

4

6

8

10

12

-3-2

.25

-1.5

-0.7

5 00.

75 1.5

2.25 3

Classe

Fre

qu

en

za

Introduzione







American-Style Options Nella seconda edizione del suo celebre libro sui derivati Hull (1993)

scriveva: “Monte Carlo simulation can only be used for European-style options”.

Sette anni e due edizioni dopo (Hull, 2000) quest’affermazione perentoria ha lasciato spazio ad una tesi più possibilista: “Monte Carlo simulation cannot easily handle situations where there are early exercise opportunities”.

Nell’intervallo di tempo trascorso fra le due edizioni molto lavoro è stato svolto per estendere il metodo Monte Carlo al trattamento di titoli derivati american-style , titoli cioè che incorporano caratteristiche legate ad opportunità di esercizio anticipato.

American-Style Options Il problema fondamentale è legato all’individuazione del tempo di

esercizio ottimale.

Negli algoritmi di simulazione classici esiste una direzione implicita del tempo. La conoscenza del prezzo ad un istante t permette il calcolo del prezzo

(o di un’altra variabile di stato) all’istante t + t, la freccia del tempo punta in avanti da cui la definizione inglese forward in time assegnata a questo tipo di simulazioni.

Al contrario le procedure di pricing di opzioni con esercizio anticipato richiedono algoritmi backward in time. Lavorando all’indietro nel tempo a partire dalla data di scadenza la

strategia ottimale di esercizio e il prezzo delle opzioni sono facilmente stimabili (ad esempio con un semplice albero binomiale).

American-Style Options Un approccio estremamente semplice e potente è quello dovuto a

Longstaff e Schwartz che discuteremo nel seguito di questa sezione Longstaff F. A., E. S. Schwartz

“Valuing American Options by Simulation: A Simple Least-Squares Approach”

The Review of Financial Studies, Spring 2001, Vol. 14, No. 1, pp.113-147

E estremamente suggestivo porre questo modello nel contesto più ampio della programmazione dinamica, per questo, seguendo P. Brandimarte, forniremo una breve introduzione a questa disciplina. P. Brandimarte

“Numerical Methods in Finance”Wiley Series in Probability and Statistics

Molto interessante il sito dell’autore www.polito.it/~brandimarte in cui si trovano numerosi supplementi ed aggiunte al libro oltre a tutti i sorgenti degli esempi matlab contenuti nel testo.

Pricing di Opzioni Americane La Programmazione Dinamica (cenni)

Perché il Monte Carlo

I metodi basati su alberi e quelli alle differenze finite, si estendono in maniera semplice e naturale al caso di esercizio anticipato.

Tuttavia ci sono limiti piuttosto stringenti al numero di fattori aleatori che questi modelli possono prendere in considerazione.

Al contrario la forza del metodo Monte Carlo risiede proprio nella possibilità di trattare problemi ad alta dimensionalità con un numero elevato di variabili aleatorie.

Programmazione Dinamica e Pricing di Opzioni Americane

Prezzare un’opzione con esercizio anticipato equivale a dover risolvere un problema di ottimizzazione dinamica;

Tuttavia anche i metodi classici della programmazione dinamica sono soggetti alla “maledizione” della dimensionalità;

Ecco perché sono estremamente importanti strategie di soluzione, sullo stile di quella di Longstaff e Schwartz, che di fatto risultano applicazioni ibride in cui al problema di programmazione dinamica viene affiancata una procedura di simulazione.

Il problema del cammino più breve

0

1

2

3

4

5

6 7

7

64

2

1

3

1

7

2

10

8

5

1

Una rete (network) è composta da un insieme di nodi e un insieme di archi che congiungono coppie di nodi.


Lo scopo è quello di ricercare il cammino minimo nella rete che congiunge il nodo 0 al nodo 7;

Il punto di partenza consiste nel trovare un’opportuna formulazione della soluzione ottimale;

Indichiamo con F(j) la lunghezza del cammino più corto che congiunge il nodo 0 con il nodo j, supponiamo poi che il nodo i si trovi lungo tale cammino. Si può allora dimostrare che il cammino ottimale che congiunge 0 ad i è un sottoinsieme del cammino ottimale che congiunge 0 a j;

In altre parole la soluzione ottimale al problema può essere trovata in modo “ricorsivo” assemblando soluzioni ottimali di problemi via via più semplici.


In altre parole possiamo scrivere

Dove A è l’insieme degli archi e N è l’insieme dei nodi.

Per trovare la soluzione ottimale, dovremmo considerare il cammino ottimale dal

nodo 0 a tutti i nodi che sono predecessori di j, dopodichè calcoliamo per ciascuno

di questi la somma del percorso ottimale fino ad i più il singolo step da i a j.

Questo tipo di equazione ricorsiva, la cui forma esatta dipende dal problema in

esame, è il cuore della programmazione dinamica ed è un esempio di equazione

funzionale.

NA

jciFjF ijji

)(min)(),( N

A

jciFjF ij

ji )(min)(

),(


16

810

511

108

min

)6(

)5(

)4(

min)7(

10111

19min

)5(

)3(min)6(

1129

76min

)3(

)2(min)5(

839

17min

)3(

)1(min)4(

946

27min

)2(

)1(min)3(

6}06min{)2(

7}07min{)1(

67

57

47

56

36

35

25

34

14

23

13

cF

cF

cF

F

cF

cFF

cF

cFF

cF

cFF

cF

cFF

F

F

5

01

7

3

Il problema del cammino più breve L’equazione ricorrente vista prima è un’equazione di tipo forward, possiamo

avere situazioni in cui è più conveniente avere equazioni di tipo backward in cui cioè si consideri la lunghezza del percorso da i ad N invece che da 0 ad i

La funzione valore B(i), valutata in un certo punto i, esprime il costo futuro ottimale della strategia condizionale al fatto di essere giunti nel nodo i-esimo.

Lo stesso tipo di formulazione può essere sviluppato in presenza di incertezza nel qual caso la funzione valore non sarà una funzione deterministica ma un valore di aspettazione. Connessione con gli alberi binomiali e trinomiali.

NA

ijBciB ijji

)(min)(),( N

A

ijBciB ij

ji )(min)(

),(

Processi decisionali sequenziali

Consideriamo un sistema dinamico discreto modellato dalla seguente equazione di stato

dove xt è il vettore delle variabili di stato alla fine dell’intervallo di tempo t e ut è il vettore delle variabili di controllo applicate durante l’intervallo t.

Per il momento non introduciamo alcuna incertezza, dato il valore delle variabili di stato all’istante t-1, xt-1, le funzioni ht e le variabili di controllo ut,

possiamo determinare con precisione il futuro stato del sistema.

T21tuxhx t1ttt ,,, ),(

T21tuxhx t1ttt ,,, ),(

Un processo decisionale sequenziale consiste nella selezione di una sequenza di controlli ut su un orizzonte temporale definito al fine di minimizzare una data funzione obiettivo con

l’eventuale presenza di vincoli sia sulle variabili di stato che sulle variabili di controllo. La funzione da minimizzare dipende in generale sia dalla traiettoria nello spazio delle fasi

del sistema (x1, x2, … ,xn) sia dalla sequenza delle variabili di controllo applicate

(u1, u2, …, un)

In generale possiamo descrivere il nostro problema come

Processi decisionali sequenziali

T21t0uxg

T21tuxhxts

uxf

t1tt

t1ttt

T

1tt1tt

u

,,, ),(

,,, ),( ..

),(min

T21t0uxg

T21tuxhxts

uxf

t1tt

t1ttt

T

1tt1tt

u

,,, ),(

,,, ),( ..

),(min

Il principio di ottimalità

La funzione obiettivo appena introdotta gode di un’importante proprietà:

essa è separabile.

Per un dato numero r il contributo degli ultimi r step dipende solo dallo

stato corrente xT-r e dagli r controlli (uT-r+1,…, uT) (Markovian State Property).

Una simile proprietà vale anche per la traiettoria nel senso che lo stato

xt+1 raggiunto a partire da xt applicando il controllo ut+1 dipende soltanto da

xt e da ut+1 e non dalla storia passata.

Una conseguenza di tale separabilità è il cosiddetto principio di ottimalità.

Il principio di ottimalità

Una strategia ottimale *** ,,, T21 uuu è tale che, comunque sia scelto lo stato iniziale 0x e il primo controllo *

1u , i controlli successivi ** ,, T2 uu individuano una strategia ottimale per il problema a 1T stadi con stato iniziale 1x ottenuto applicando il primo controllo *

1u .

1011T101u

0T uxhFuxfxF1

,,min)( 1011T101u

0T uxhFuxfxF1

,,min)(

)( 1ttt xu )( 1ttt xu

Il principio di ottimalità Presenza di incertezza; Non è possibile determinare una sequenza di controlli ottimale; Occorre introdurre una sequenza di controlli che dipende anche

dalla variabile di stato:

1011T01u

00 uxxVEuxfxV1

,)(),(min)( 1011T01u

00 uxxVEuxfxV1

,)(),(min)(

Pricing di Opzioni Americane Il modello di Longstaff-Schwartz(Least Squares Monte Carlo)

Il pricing di opzioni americane In ogni punto del path simulato il valore del contratto deve essere valutato

come il massimo fra il valore di continuazione del contratto stesso ed il payoff che conseguirebbe ad un esercizio immediato;

Il valore di continuazione è definito come il valore atteso dei flussi futuri scontati al tasso risk free assumendo che una strategia di esercizio ottimale venga esercitata nel futuro;

Ad esempio per un’opzione scritta su un singolo asset abbiamo

i1i1itrQ

iiiii SSVeESISV )(),(max)( i1i1itrQ

iiiii SSVeESISV )(),(max)(

C’è tuttavia una “piccola” difficoltà; Per definire la strategia ottimale dovremmo conoscere il valore di aspettazione del

valore futuro dell’opzione condizionale al prezzo allo stato attuale ma questo dipende dalla decisione assunta allo step successivo;

In altre parole per conoscere la strategia ottimale adesso dovremmo conoscere la strategia ottimale futura;

Una possibile soluzione consiste nel trovare un’approssimazione del valore di aspettazione della funzione valore;

Una possibilità è quella di approssimare il valore di aspettazione condizionale con una regressione verso il valore corrente del titolo utilizzando una opportuna base di funzioni;

Es. sviluppo polinomiale

Il pricing di opzioni americane

2i3i21i1i1i

trQi SaSaaSSVeE

)( 2i3i21i1i1i

trQi SaSaaSSVeE

)(

Il pricing di opzioni americane

Per ottenere i valori dei coefficienti incogniti, l’algoritmo proposto da Longstaff e Schwartz utilizza una regressione in cui i valori della variabile y sono rappresentati dai payoff

futuri su ciascun path e i valori della x sono rappresentati dal valore del

sottostante all’istante immediatamente precedente sullo stesso path;


Il metodo di Longstaff e SchwartzIl metodo di Longstaff e Schwartz

Introduzione







Stima del Value-at-Risk

Per valutare tale probabilità col metodo Monte Carlo possiamo simulare una serie di valori per tutti i fattori di rischio cui è sensibile il portafoglio.

La rivalutazione del portafoglio per ogni serie di nuovi valori dei fattori di rischio viene ripetuta un numero sufficientemente elevato di volte al fine di ottenere una distribuzione significativa della perdita.

Dalla distribuzione ottenuta si risale così al VaR come percentile al livello di probabilità predefinito.

Value-at-Risk I metodi di calcolo del Value-at-Risk basati sull’uso della matrice di

varianza-covarianza permettono un calcolo estremamente efficiente dal punto di vista computazionale (Cherubini e Della Lunga, 2001).

Tali metodi, tuttavia, non permettono un calcolo accurato nel caso in cui il portafoglio in esame contenga titoli con payoff non lineare come, ad esempio, delle opzioni.

L’introduzione di opportune correzioni (nel caso delle opzioni l’inclusione del Gamma) porta a stime migliori a scapito tuttavia di una maggior complicazione nella definizione degli algoritmi di calcolo e comunque risulta sempre insufficiente nel caso di portafogli fortemente non lineari.

Inoltre tale metodologia risulta del tutto inapplicabile in quei casi in cui all’interno del portafoglio siano presenti prodotti per i quali la relazione fra valore e fattori di rischio sia non monotona. In questi casi non è affatto detto che le perdite più elevate si verifichino in corrispondenza di movimenti estremi dei fattori di mercato. Un esempio di un prodotto di questo tipo è lo straddle.

Value-at-RiskSebbene le problematiche di base siano le stesse, l’utilizzo delle tecniche Monte Carlo per la stima del VaR presenta alcune peculiarità che lo distinguono rispetto all’utilizzo in altri settori della matematica finanziaria.

Nelle applicazioni di pricing, ad esempio, il Monte Carlo viene utilizzato prevalentemente per stimare il valore di aspettazione di una funzione di una o più variabili aleatorie.

Nella stima del VaR, invece, quello che di solito interessa è rappresentare nel modo più preciso possibile la distribuzione di probabilità delle variazioni di valore del portafoglio al variare dei fattori di rischio.

Due punti meritano di essere sottolineati; il primo è che le variazioni dei fattori di rischio devono essere generate sulla

base di distribuzioni il più possibile rappresentative delle reali distribuzioni empiriche osservate nella realtà;

il secondo punto riguarda la necessità di postulare delle relazioni di pricing tramite le quali ricavare le variazioni di valore del portafoglio dalle variazioni, simulate, dei fattori di rischio. Il calcolo del VaR richiede quindi la scelta di un preciso modello finanziario.

Simulazione Montecarlo

1. Scelta orizzonte temporale di riferimento2. Scelta modelli di pricing (relazione fattori di

rischio/singoli strumenti)3. Scelta processi stocastici 4. Generazione di un elevato numero di scenari

(10.000?) dei fattori di rischio (tassi,cambi,volatilità,…)

5. Calcolo delle ipotetiche variazioni di valore del portafoglio per ogni scenario (es: Taylor)

6. Ordinamento dei risultati7. Individuazione k-esimo percentile

Simulazione MontecarloDati

Storici/Impliciti

Dati Storici/Impliciti

Simulazioni Simulazioni

VaRVaR

Parametri dei modelli

Parametri dei modelli

Modelli Stocastici

Modelli Stocastici

Valutazione “completa”

Valutazione “completa”

Distribuzione P&LDistribuzione P&L

Posizioni Portafoglio

Posizioni PortafoglioModelli attività fin.Modelli attività fin.

Value-at-Risk Individuare la procedura appropriata per la simulazione delle

variazioni dei fattori di rischio rappresenta solo uno dei problemi in gioco.

Un altro aspetto fondamentale è senz’altro quello relativo alla metodologia di rivalutazione del portafoglio.

Due sono le principali alternative: rivalutazione esatta (full evaluation) rivalutazione parametrica.

Nel primo caso è necessario valutare attentamente le risorse di calcolo necessarie per ricalcolare migliaia di volte il valore di ciascun contratto.

Questo metodo, sebbene garantisca la miglior precisione, è spesso talmente oneroso in termini di risorse di calcolo da spingere verso metodi di tipo parametrico meno precisi ma sicuramente più efficienti in termini di velocità di elaborazione.

Value-at-Risk:

Rivalutazione Parametrica del Portafoglio In questo approccio la relazione funzionale fra la variazione del

valore del portafoglio e la variazione dei fattori di mercato è rappresentata in termini di un numero finito di parametri.

Tali parametri vengono ricavati da un numero relativamente piccolo di rivalutazioni esatte di ogni asset in portafoglio.

I due metodi più diffusi di rivalutazione parametrica del portafoglio sono basati su uno sviluppo in serie di Taylor

o, alternativamente, sulla cosiddetta griglia delle sensitività.

Value-at-Risk:

Rivalutazione Parametrica del Portafoglio Tramite lo sviluppo in serie di Taylor possiamo esprimere la variazione nel valore di

una funzione di n variabili sotto forma di polinomio con infiniti termini i cui coefficienti

sono le derivate parziali della funzione rispetto a ciascuna variabile.

Nelle applicazioni pratiche lo sviluppo del polinomio viene poi arrestato ad un numero

finito di termini in funzione della precisione richiesta.

Il valore di mercato del nostro portafoglio è naturalmente dato dalla somma dei valori

di ciascuna componente, quindi indicando con PVm la funzione di rivalutazione esatta

di ciascun asset possiamo scrivere

dove Xj(t) sono i valori dei fattori di mercato da cui dipende il valore del titolo

j

jj

mportfolio ttXPMPV ),()( j

jj

mportfolio ttXPMPV ),()(

Value-at-Risk:

Rivalutazione Parametrica del Portafoglio La variazione di valore del portafoglio può essere quindi

espressa come

lkjlkj

lkj

kjkj

kjjj

j

kjkj kjj j

jj

jj

portfolio

XXXXXX

PV

XXXX

PVX

X

PV

XXXX

PVX

X

PVX

X

PVPV

3

22

33

3

3

22

2

2

6

1

6

1

6

1

2

1

2

1

lkjlkj

lkj

kjkj

kjjj

j

kjkj kjj j

jj

jj

portfolio

XXXXXX

PV

XXXX

PVX

X

PV

XXXX

PVX

X

PVX

X

PVPV

3

22

33

3

3

22

2

2

6

1

6

1

6

1

2

1

2

1

I coefficienti del polinomio sono le derivate parziali del valore di portafoglio rispetto ai fattori di mercato. Questi coefficienti vengono calcolati come somma delle derivate parziali del valore di ciascun titolo rispetto ai fattori, ciascuna di queste derivate è calcolata dalle formule esatte di rivalutazione dei titoli. Il numero di termini da introdurre nello sviluppo in serie di Taylor dipende dal grado di accuratezza con cui è necessario stimare il VaR e dal grado di non linearità del portafoglio.

I coefficienti del polinomio sono le derivate parziali del valore di portafoglio rispetto ai fattori di mercato. Questi coefficienti vengono calcolati come somma delle derivate parziali del valore di ciascun titolo rispetto ai fattori, ciascuna di queste derivate è calcolata dalle formule esatte di rivalutazione dei titoli. Il numero di termini da introdurre nello sviluppo in serie di Taylor dipende dal grado di accuratezza con cui è necessario stimare il VaR e dal grado di non linearità del portafoglio.

Value-at-Risk:

Rivalutazione Parametrica del Portafoglio Nell’approccio basato sulla griglia di sensitività il punto di partenza è

costituito dalla valutazione esatta della variazione del portafoglio in funzione dei fattori di rischio per un insieme finito e discreto di valori di variazione dei fattori di rischio stesso.

Questo porta alla realizzazione di una griglia multidimensionale in cui ogni punto dell’iperpiano specifica una certa variazione congiunta dei fattori di rischio e il valore della funzione di rivalutazione in quel punto è il valore esatto della variazione del portafoglio corrispondente.

Una volta che la griglia è stata costruita le variazioni del valore di portafoglio rispetto ad una variazione arbitraria dei fattori di rischio possono essere ottenute interpolando linearmente le variazioni calcolate nei punti della griglia stessa.

Funzione Caratteristica

la funzione caratteristica di una variabile aleatoria x è definita come

è la Trasformata di Fourier della funzione densità di probabilità

si definisce poi la funzione generatrice dei momenti

dxexf xi )()(

dxexf xi )()(

)()( ,)()( idxexfs sx

)()( ,)()( idxexfs sx

Funzione Caratteristica

Dalla definizione è evidente che

inoltre

il che giustifica il nome di funzione generatrice dei momenti.

xisx eEeEs )( ,)( xisx eEeEs )( ,)(

)()()( )()( nnnsxnn mxE0exEs )()()( )()( nnnsxnn mxE0exEs

Value-at-RiskSimulazione Monte Carlo per la stima dello Shortfall (Glassermann & C.)

Qualora si consideri valida l’approssimazione Delta-Gamma, la variazione di valore di un portafoglio può essere sempre espressa nella forma

n

iiiii

tt ZZbZZZbQ1

2

essendo C una qualunque matrice per cui valga tCC ( è la matrice di varianza e covarianza), è la matrice diagonale degli autovalori di

CC t2

1 , Cb ,e, infine, Z indica un vettore di variabili indipendenti

distribuite secondo una normale standard.

In generale possiamo quindi definire, sempre nel contesto dell’approssimazione Delta-Gamma, una funzione di perdita (loss function) definita come

QaL 0 dove a0 è una costante. Il VaR al livello di confidenza p è quel valore xp tale che

pxLP p )( Indicando con I(x) la funzione indicatrice possiamo scrivere

dzzfxLIxLP )()()(

essendo f(z) la funzione densità di probabilità congiunta dei fattori di rischio.



Come abbiamo già discusso la tecnica di Importance Sampling consiste nell’effettuare un cambio di misura ricercando una funzione densità di campionamento g(z) che presenti caratteristiche tali da ridurre la varianza come nel caso in esame rendendo più probabili gli eventi corrispondenti alle code della distribuzione originaria

)]()([~

)()(

)()()( zlxLIEdzzg

zg

zfxLIxLP

avendo indicato con E

~ il valore di aspettazione sotto la nuova misura g(z) e con l(z) la importance function.


Nel nostro caso la funzione di distribuzione iniziale è la normale standard multivariata

ZZzf t

2

1exp)(

La nuova densità di probabilità g(z) verrà ricercata nella forma di una distribuzione multivariata normale con media e matrice di varianza-covarianza B parametrizzate in funzione di un parametro > 0.

ZBZBzg t 12/1

2

1exp)(

121

B

bB


possiamo allora scrivere l’espressione della important function

ZBZB

ZZ

zlt

t

12/1

2

1exp

2

1exp

)(

QZZZbzl tt expexp)(

dove è data da

n

ii

i

ib

1

2

21ln212

1

Si può verificare che è il logaritmo della funzione generatrice di Q.


Supponiamo ora che l’approssimazione quadratica sia esatta, cioè che sia L = a0 + Q allora poiché L > x implica Q > x – a0 possiamo ricavare un limite superiore per la probabilità P(L > x)

00 exp)()()(~

)( axxQaIeEzlxLIExLP Q

non è difficile ricavare anche un limite superiore per il momento secondo

02

2 2)(2exp)()(~

),( axzlxLIExm


Al fine di ottimizzare la simulazione occorre trovare il valore di per cui m2(x, ) è minimo o, in alternativa, il minimo del limite superiore che si ottiene per x tale che

0)( axx Vediamo ora quanto vale il valore di aspettazione di Q sotto la nuova misura g(z)

dzzfeQedzzfeQdzzgQQE QQ )()()(~


Per definizione la funzione generatrice dei momenti di Q è data da

dzzfeeE QQ )(

da cui

dzzfeQd

d Q )(


quindi, ricordando che ln ,

d

deQE

~

Il valore di aspettazione di Q valutato con x è pertanto pari a x – a0 e quindi il valore di aspettazione di L è proprio pari a x. Pertanto sotto la nuova misura l’evento L > x non è più un evento raro.

Date post:	17-Jan-2015
Category:	Education
Upload:	giovanni-della-lunga
View:	1,437 times
Download:	5 times

Lezione 3 metodo monte carlo

Education