Introduzione al metodo monte carlo

Giovanni Della Lunga

Dipartimento di Chimica Centro Interdipartimentale per lo Studio dei Sistemi Complessi

Università degli Studi di Siena

Introduzione al metodo Monte Carlo

Concetti di Base e Applicazioni Finanziarie

INDICE* _______________________________________________ 1 Richiami di Teoria delle Probabilità.................................................................................................1

1.1 Definizioni....................................................................................................................................1 1.2 Variabili Aleatorie e Funzioni di Distribuzione ...........................................................................2 1.3 Momenti .......................................................................................................................................4 1.4 Distribuzioni Discrete ............................................................................................................8 1.5 Distribuzioni Continue.........................................................................................................12 1.6 Covarianza e Correlazione..........................................................................................................23 1.7 Che cos’è un Processo Stocastico...............................................................................................23

2 La Simulazione Monte Carlo...........................................................................................................25 2.1 Simulazione Monte Carlo e Integrazione Numerica...................................................................29

3 Generatori di Numeri Casuali .........................................................................................................33 3.1 Tre tipi di generatori ...................................................................................................................33 3.2 Il generatore Lineare Congruente..........................................................................................34 3.3 Le funzioni Rnd() e Randomize()...............................................................................................35 3.4 Test statistici sui generatori di numeri casuali ......................................................................37

4 Generazione di Distribuzioni non Uniformi...................................................................................41 4.1 Il Transformation Method ..........................................................................................................41

4.1.1 Utilizzo delle funzioni del foglio di lavoro di Microsoft Excel in Visual Basic.................42 4.1.2 Come Simulare una variabile distribuita Normalmente ...............................................43

4.2 La simulazione di un processo di Wiener ............................................................................45 4.3 Il Rejection Method ..............................................................................................................48 4.4 Generazione di variabili con correlazione assegnata ............................................................49 4.5 Come adattare una distribuzione non normale univariata ai dati empirici..................................51 4.6 Appendice A - Il condizionamento della matrice di covarianza.................................................53

5 La Valutazione dei Titoli Derivati...................................................................................................55 5.1 La stima del prezzo ...............................................................................................................55 5.2 La stima delle sensitività.......................................................................................................57

5.2.1 Approssimazione basata sulle differenze finite ..................................................................57 5.2.2 Stime dirette........................................................................................................................60

5.3 Un cenno alle problematiche dei derivati di tipo americano ......................................................62 5.4 Appendice A - Funzioni VB per la simulazione del Delta .........................................................64

6 La stima del Value-at-Risk ..............................................................................................................67 6.1 Le problematiche ........................................................................................................................67 6.2 Il rischio di mercato....................................................................................................................68

6.2.1 La rivalutazione parametrica del portafoglio......................................................................69

7 Metodi di riduzione della varianza .................................................................................................71 7.1 Introduzione................................................................................................................................71 7.2 Il metodo delle variabili antitetiche.......................................................................................72 7.3 Il metodo della variabile di controllo ....................................................................................74

* I paragrafi preceduti dal simbolo contengono argomenti sviluppati anche nelle cartelle Excel che accompagnano il testo.

Indice ii

7.3.1 Un esempio pratico: il pricing di un’opzione asiatica.........................................................75

7.4 Moment Matching ......................................................................................................................78 7.5 Importance Sampling..................................................................................................................79

7.5.1 Importance Sampling e stima del VaR per portafogli contenenti opzioni ..........................80 7.6 Campionatura uniforme..............................................................................................................83

7.6.1 Il Campionamento Stratificato ......................................................................................83 7.6.2 Il metodo dell’Ipercubo Latino .....................................................................................85 7.6.3 Le Sequenze a Bassa Discrepanza ................................................................................86

La sequenza di Halton ....................................................................................................................89 La sequenza di Faure ......................................................................................................................90

7.6.4 Metodi misti e decomposizione spettrale per problemi ad alta dimensionalità ..................91 7.7 Appendice A - Deduzione della formula (5.9) ...........................................................................92

BIBLIOGRAFIA........................................................................................................................... 97

1 Richiami di Teoria delle Probabilità

Possibile, ma non interessante – rispose Lönnrot -. Lei dirà che la realtà non ha il minimo obbligo di essere interessante. Io replicherò che se la realtà può sottrarsi a quest’obbligo, non possono sottrarvisi le ipotesi. Jorge Luis Borges, “La morte e la bussola”

1.1 Definizioni In questo paragrafo passeremo in rassegna i concetti del Calcolo delle Probabilità e della Teoria delle Distribuzioni. Naturalmente questa presentazione è ben lungi dall’essere completa e vuol soltanto essere di sussidio al lettore raccogliendo in maniera sintetica i principali strumenti indispensabili per la comprensione dei modelli di pricing di prodotti derivati e, più in generale, dei modelli di simulazione stocastica utilizzati in tale settore. Ogni tentativo di dare una definizione rigorosa dei concetti probabilistici più elementari si trova di fronte ad un problema; infatti, non solo esistono differenti formalizzazioni e assiomatizzazioni della probabilità ma a queste corrispondono, in generale, molteplici nozioni intuitive di probabilità spesso assai diverse fra loro. Con queste parole di avviso iniziamo a definire i principali concetti della teoria delle probabilità rinviando il lettore in cerca di formalizzazioni rigorose a testi specifici. Al di la delle differenze di carattere formale un elemento comune posseduto da tutte le forme di probabilità riguarda il suo significato intuitivo di valutazione della possibilità che un dato evento possa accadere o meno. Sia nelle scienze naturali sia in quelle economiche si è soliti assumere che un certo evento sia il risultato di un ipotetico esperimento intendendo con questo termine l’insieme di tutte “le azioni e le condizioni ambientali che conducono al determinarsi di un fatto1”. E’ un esperimento la misura di una grandezza fisica, il lancio di un dado o di una moneta, il verificarsi o meno di un particolare stato di natura (es. l’indice MIB30 supera il livello 50.000). Indicheremo con ω un particolare stato di natura esito di un dato esperimento e con Ω l’insieme di tutti gli stati possibili (spazio campione). Il concetto di evento é associato al verificarsi di uno o più stati di natura, esso verrà pertanto rappresentato come sottoinsieme di Ω. Lo spazio degli eventi, A, è quindi una famiglia di sottoinsiemi di Ω caratterizzata dalle seguenti proprietà:

1. Ω ∈ A; 2. se l’evento ω ∈ F allora anche il suo complemento Ω - ω ∈ A; 3. se ωn ∈ A, allora ωU

n

n ∈ A

1 O. Lessi “Corso di Probabilità”, pag. 21, Metria, 1993.

2 Introduzione al metodo Monte Carlo

Una famiglia di sottoinsiemi che soddisfa a queste proprietà viene chiamata σ-algebra. Diremo quindi che A è una σ-algebra di Ω. Definiamo funzione di probabilità una funzione P a valori reali che soddisfa le seguenti proprietà:

1. ∈∀≥ ωω ,0)(P A ; 2. ; 1)( =ΩP

3. ( )∑∞

=

∞

=

=

11 nn

nn PP ωωU

se gli ωn sono a due a due disgiunti cioè se ∅=∩ ji ωω per ogni ji . ≠Osserviamo che una funzione di probabilità così definita è anche una misura. La terna (Ω, A , P) viene detta spazio di probabilità. Esempio 1.1 – Consideriamo il caso di un indice azionario, ad esempio il MIB30, il cui valore in un dato intervallo di tempo possa solo aumentare o diminuire di una quantità assegnata (naturalmente si tratta di una schematizzazione). Indicando con u l’evento corrispondente all’aumento dell’indice e con d l’evento contrario, possiamo scrivere:

du, =Ω A = du , ,, ∅Ω Esempio 1.2 – Anche se questi sono solo dei brevi richiami non possiamo non citare l’esperimento aleatorio per antonomasia: il lancio di un dado. In questo caso lo spazio campione è formato dall’insieme dei sei numeri che possono risultare dal lancio stesso

6 5, 4, 3, 2, ,1 =Ω Vediamo il significato di alcuni elementi di A: l’elemento 2,1 1 =ω corrisponde all’evento “il numero risultante dal lancio è minore o uguale a 2”, 5,3,1 2 =ω corrisponde a “il numero risultante è dispari”, 6,4,2 3 =ω rappresenta l’evento “il numero uscente è pari”. 1.2 Variabili Aleatorie e Funzioni di Distribuzione Dato uno spazio di probabilità, una variabile aleatoria o casuale viene definita come una funzione tale che per ogni numero reale r si abbia ℜ→Ω:X

( ) ∈≤Ω∈ X rωω : A. La funzione ( ) ( ) ( ) ( )xXPxXPxFX ≤Ω∈=≤= ωω : , definita sull’insieme dei numeri reali, viene detta funzione di distribuzione cumulata o, più semplicemente, funzione di distribuzione. Esempio 1.3 – Con riferimento all’Esempio 1.1 sia X una funzione così definita X(u) = 0

Introduzione al metodo Monte Carlo 3

X(d) = 1 La funzione X è una variabile aleatoria in quanto

∅=≤< rXr )(: 0 ωω arXr )(: 10 =≤<≤ ωω Ω=≤≤ rXr )(: 1 ωω

Esempio 1.4 – Supponiamo che la variabile definita nell’esempio precedente assuma il valore 0 con probabilità p1 e 1 con probabilità p2 = 1 – p1. La funzione di distribuzione di questa variabile è data da

≤<≤

<=

xxp

xxFX

1:110:

0:0)( 1

La funzione indicatrice di un evento A, IA, definita come

∈∉

=AA

I A ωω

ω:1:0

)(

è un tipico esempio di funzione a due valori con p1 = P(A). Una variabile aleatoria è detta discreta se l’insieme dei valori che può assumere è numerabile. Sia (Ω, A , P) uno spazio di probabilità e X una variabile aleatoria discreta che può assumere i valori Definiamo la funzione di probabilità come KK ,,,, 21 kxxx

===

=altrimenti:0

,....2,1 qualcheper , se:)()(

ixxxXPxf i

X (1.1)

La funzione di probabilità e la funzione di distribuzione sono legate dalla relazione:

∑≤

=xix

iXX xfxF )()( (1.2)

Una variabile aleatoria X è detta continua se esiste una funzione reale fX tale che per ogni x reale sia soddisfatta la relazione

∫∞−

=x

XX dyyfxF )()( (1.3)

nei punti in cui la funzione di distribuzione è derivabile vale anche la relazione inversa


dxxdFxf X

X)()( = (1.4)

La funzione f(x) in questo caso viene detta funzione densità di probabilità (o semplicemente funzione densità). 1.3 Momenti Il valor medio o valore di aspettazione di X, che indicheremo con Xµ , è definito come

[ ] ∑==i

iXiX xfxXE )(µ (1.5)

In generale si definisce momento dall’origine (o momento grezzo) di ordine r, e si indica , la media della variabile aleatoria Xr

Xµ r. La definizione è naturalmente applicabile solo nel caso in cui tale media sia finita. La varianza di X, indicata con , è la media degli scarti quadratici rispetto alla media

2Xσ

Xµ e rappresenta una misura di dispersione di X. La sua radice quadrata è detta deviazione standard. La varianza è definita da

( )∑ −=i

iXXiX xfx )(22 µσ (1.6)

da cui è immediato ricavare

[ ] [ ]( 222 XEXEX −=σ ) (1.7) Un ruolo molto importante è rivestito dalla cosiddetta funzione generatrice dei momenti:

[ ]

∑=+++=

=

++++==

i

iiXXX

tXX

ti

tt

XtXtXtEeEtm

µµµ!1

211

)(!3

1)(211)(

22

32

K

K

(1.8)

dove abbiamo utilizzato lo sviluppo in serie della funzione esponenziale. E’ estremamente semplice verificare che

0

2

=

=t

Xr

r

X mdtdµ (1.9)

Si possono facilmente generalizzare al caso continuo i risultati per le distribuzioni discrete, il valore di aspettazione sarà pertanto definito come


∫+∞

∞−

= dxxxfXE X )(][ (1.10)

in maniera analoga si generalizzano le definizioni di varianza e degli altri momenti. La funzione generatrice dei momenti gode delle stesse proprietà del caso discreto. Un problema che citiamo solo per completezza ma che non può essere affrontato in questo testo è quello relativo alla stima dei parametri statistici. Nel seguito vedremo alcune formule (presenteremo quelle impiegate da Excel) che permettono di stimare il valore di uno o più momenti a partire da un campione di dati. Occorre sempre ricordare che campioni diversi forniscono stime diverse del parametro di riferimento. La stima di un parametro statistico non restituisce il valore esatto del parametro del campione, quanto una variabile casuale, dotata di una propria distribuzione di probabilità: la distribuzione di probabilità delle stime dei parametri è nota in statistica come distribuzione campionaria. Ci limitiamo a ricordare che per avere stime non distorte è necessario che il valore di aspettazione dello stimatore sia uguale al parametro stimato. Excel dispone di numerose funzioni per l’analisi statistica dei dati (vedi Tabella 1.1). Fra queste vi sono una serie di funzioni per il calcolo dei momenti che vale la pena descrivere in maggior dettaglio.

Tabella 1.1 – Le funzioni statistiche di Excel MEDIA.DEV Restituisce la media delle deviazioni assolute dei valori rispetto alla loro media MEDIA Restituisce la media degli argomenti MEDIA.VALORI Restituisce la media degli argomenti, inclusi i numeri, il testo e i valori logici DISTRIB.BETA Calcola la funzione densità di probabilità cumulativa beta INV.BETA Calcola l'inverso della funzione densità di probabilità cumulativa beta DISTRIB.BINOM Restituisce la distribuzione binomiale per il termine individuale DISTRIB.CHI Restituisce la probabilità ad una coda per la distribuzione del chi quadrato INV.CHI Restituisce l'inverso della probabilità ad una coda per la distribuzione del chi quadrato TEST.CHI Restituisce il test per l'indipendenza CONFIDENZA Restituisce l'intervallo di confidenza per la popolazione CORRELAZIONE Restituisce il coefficiente di correlazione tra due serie di dati CONTA.NUMERI Conta la quantità di numeri nell'elenco di argomenti CONTA.VALORI Conta il numero di valori nell'elenco di argomenti COVARIANZA Calcola la covarianza, la media dei prodotti delle deviazioni accoppiate CRIT.BINOM Restituisce il più piccolo valore per il quale la distribuzione cumulativa binomiale risulta maggiore o uguale

ad un valore di criterio DEV.Q Restituisce la somma dei quadrati delle deviazioni DISTRIB.EXP Restituisce la distribuzione esponenziale DISTRIB.F Restituisce la distribuzione di probabilità F INV.F Restituisce l'inverso della distribuzione della probabilità F FISHER Restituisce la trasformazione di Fisher INV.FISHER Restituisce l'inverso della trasformazione di Fisher PREVISIONE Restituisce i valori lungo una tendenza lineare FREQUENZA Restituisce la distribuzione di frequenza come matrice verticale TEST.F Restituisce il risultato di un test F DISTRIB.GAMMA Restituisce la distribuzione gamma GAMMAINV Restituisce l'inverso della distribuzione cumulativa gamma LN.GAMMA Restituisce il logaritmo naturale della funzione gamma, Γ(x) MEDIA.GEOMETRICA Restituisce la media geometrica CRESCITA Restituisce i valori lungo una linea di tendenza esponenziale MEDIA.ARMONICA Restituisce la media armonica DISTRIB.IPERGEOM Restituisce la distribuzione ipergeometrica INTERCETTA Restituisce l'intercetta della retta di regressione lineare CURTOSI Restituisce la curtosi di un insieme di dati


GRANDE Restituisce il k-esimo valore più grande in un insieme di dati REGR.LIN Restituisce i parametri di una tendenza lineare REGR.LOG Restituisce i parametri di una linea di tendenza esponenziale INV.LOGNORM Restituisce l'inverso di una distribuzione lognormale DISTRIB.LOGNORM Restituisce la distribuzione lognormale cumulativa MAX Restituisce il valore massimo in un elenco di argomenti MAX.VALORI Restituisce il valore massimo in un elenco di argomenti, inclusi i numeri, il testo e i valori logici MEDIANA Restituisce la mediana dei numeri specificati MIN Restituisce il valore minimo in un elenco di argomenti MIN.VALORI Restituisce il più piccolo valore in un elenco di argomenti, inclusi i numeri, il testo e i valori logici MODA Restituisce il valore più comune in un insieme di dati DISTRIB.BINOM.NEG Restituisce la distribuzione binomiale negativa DISTRIB.NORM Restituisce la distribuzione cumulativa normale INV.NOR Restituisce l'inverso della distribuzione cumulativa normale DISTRIB.NORM.ST Restituisce la distribuzione cumulativa normale standard INV.NORM.ST Restituisce l'inverso della distribuzione cumulativa normale standard PEARSON Restituisce il coefficiente del momento di correlazione di Pearson PERCENTILE Restituisce il k-esimo dato percentile di valori in un intervallo PERCENT.RANGO Restituisce il rango di un valore in un insieme di dati come percentuale PERMUTAZIONE Restituisce il numero delle permutazioni per un determinato numero di oggetti POISSON Restituisce la distribuzione di Poisson PROBABILITÀ Calcola la probabilità che dei valori in un intervallo siano compresi tra due limiti QUARTILE Restituisce il quartile di un insieme di dati RANGO Restituisce il rango di un numero in un elenco di numeri RQ Restituisce la radice quadrata del coefficiente di momento di correlazione di Pearson ASIMMETRIA Restituisce il grado di asimmetria di una distribuzione PENDENZA Restituisce la pendenza di una retta di regressione lineare PICCOLO Restituisce il k-esimo valore più piccolo in un insieme di dati NORMALIZZA Restituisce un valore normalizzato DEV.ST Restituisce una stima della deviazione standard sulla base di un campione DEV.ST.VALORI Restituisce una stima della deviazione standard sulla base di un campione, inclusi i numeri, il testo e i valori

logici DEV.ST.POP Calcola la deviazione standard sulla base di un'intera popolazione DEV.ST.POP.VALORI Calcola la deviazione standard sulla base sull'intera popolazione, inclusi i numeri, il testo e i valori logici ERR.STD.YX Restituisce l'errore standard del valore previsto per y per ciascun valore x nella regressione DISTRIB.T Restituisce la distribuzione t di Student INV.T Restituisce l'inversa della distribuzione t di Student TENDENZA Restituisce i valori lungo una linea di tendenza lineare MEDIA.TRONCATA Restituisce la media della parte interna di un insieme di dati TEST.T Restituisce la probabilità associata ad un test t di Student VAR Stima la varianza sulla base di un campione VAR.VALORI Stima la varianza sulla base di un campione, inclusi i numeri, il testo e i valori logici VAR.POP Calcola la varianza sulla base dell'intera popolazione VAR.POP.VALORI Calcola la varianza sulla base dell'intera popolazione, inclusi i numeri, il testo e i valori logici WEIBULL Restituisce la distribuzione di Weibull TEST.Z Restituisce il livello di significatività a due code per un test z In pratica vengono comunemente utilizzati i primi quattro momenti:

• media • varianza • skewness (o asimmetria) • curtosi

La funzione MEDIA()viene utilizzata per calcolare la media aritmetica di un insieme di valori. L’argomento di questa funzione deve indicare il blocco di celle che contiene i valori di cui calcolare il valor medio. Notate che le celle contenenti il valore 0 (zero) saranno incluse nel calcolo della media mentre le celle vuote (senza alcun dato) restano escluse da questo calcolo.


La varianza può essere calcolata tramite la funzione VAR(). Anche questa funzione prende come argomento il blocco di celle contenente i dati di cui si desidera stimare la varianza. Lo stimatore utilizzato da Excel per il calcolo è il seguente

)1(1

2

1

2

−

−

=∑ ∑

= =

nn

xxns

n

i

n

iii

(1.11)

Il momento centrale di ordine 3 ci dà informazioni sul grado di asimmetria di una distribuzione attorno alla sua media ed è comunemente indicato col termine skewness. L'asimmetria positiva indica una distribuzione con una coda asimmetrica che si estende verso i valori più positivi. L'asimmetria negativa indica una distribuzione con una coda asimmetrica che si estende verso i valori più negativi. Excel mette a disposizione la funzione ASIMMETRIA() che utilizza il seguente stimatore

∑=

−

−−

n

i

i

sxx

nnn

1

3

)2)(1( (1.12)

in cui s è lo stimatore della standard deviation e x è il valor medio. Se esistono meno di tre valori o se la deviazione standard sulla base di un campione è uguale a zero, ASIMMETRIA() restituirà il valore di errore #DIV/0!. La relazione tra momento del

terzo ordine e coefficiente di asimmetria, solitamente indicato con 1β , è data da

33

1 σµβ = (1.13)

Valori positivi dell’asimmetria indicano che la distribuzione è asimmetrica per valori crescenti della variabile x (a destra) mentre un’asimmetria negativa sta ad indicare una distribuzione asimmetrica a sinistra.

Infine CURTOSI() restituisce la curtosi, indicata con 2β , di un insieme di dati. Essa è legata al momento centrale di ordine 4 dalla relazione

44

2 σµβ = (1.14)

ed è caratteristica delle cosiddette “code grasse”. Lo stimatore utilizzato da Excel riporta in realtà la cosiddetta “curtosi in eccesso” ovvero la differenza fra 2β e 3. Questo è dovuto al fatto che la distribuzione normale o gaussiana ha curtosi pari a 3 e questo indicatore viene spesso utilizzato come indice per comprendere quando la distribuzione di un insieme di dati si allontani dalla normalità.


La formula utilizzata è

( ))3)(2(

13)3)(2)(1(

)1( 2

1

4

−−−

−

−

−−−+ ∑

= nnn

sxx

nnnnn n

i

i (1.15)

in cui s è lo stimatore della standard deviation e x è il valor medio. Se sono presenti meno di quattro dati o se la deviazione standard del campione è uguale a zero, CURTOSI() restituirà il valore di errore #DIV/0!. 1.4 Distribuzioni Discrete Vediamo ora alcune fra le più comuni funzioni di distribuzioni discrete. Il lettore può seguire la nostra analisi utilizzando il file “Richiami di Statistica 1.xls”. Distribuzione Uniforme Sia X una variabile aleatoria che assume valori nel dominio dei numeri naturali 1, 2, ... , n. Diremo che tale variabile ha una distribuzione uniforme se risulta

=

==altrimenti:0

,,2,1:1),()(

nxnnxfxf XX

K (1.16)

Valor medio e varianza sono dati da:

21

2)1(1][

1

+=

+== ∑

=

nn

nnin

XEn

i

(1.17)

( )12

14

)1(6

)12)(1(4

)1(1][][22

1

22222 −

=+

−++

=+

−=−= ∑=

nnnnnnin

XEXEn

iXσ (1.18)

Distribuzione Binomiale Dati n eventi indipendenti, tutti con uguale probabilità p, sia X la variabile casuale che conta il numero totale di eventi che si verificano fra quelli possibili. X ha una distribuzione binomiale con parametri n e p. La funzione di probabilità è

iniX pp

ininif −−−

= )1()!(!

!)( per i = 0, 1, 2, ..., n (1.19)


valor medio e varianza sono dati da

[ ] npXE = (1.20)

)1(2 pnpX −=σ (1.21)

Figura 1.1 – Funzione di Probabilità e Distribuzione Cumulata Binomiale per il caso n = 6 e p = 0.5 (le due figure differiscono solo nel formato grafico, il primo è quello predefinito prodotto da Excel tramite gli Strumenti di analisi, il secondo è ottenuto modificando la distanza fra le barre come spiegato nel testo) In Excel è possibile lavorare con la distribuzione binomiale tramite la funzione DISTRIB.BINOM() la cui sintassi è la seguente: DISTRIB.BINOM(num_successi; num_prove;probabilita_s;cumulativo) dove num_successi è il numero di successi in num_prove indipendenti, probabilita_s è la probabilità di successo per ciascuna prova e cumulativo è un valore logico che determina la forma assunta dalla funzione. Se il valore cumulativo è VERO, DISTRIB.BINOM restituirà la funzione distribuzione cumulata, ovvero la probabilità che venga restituito un numero massimo di successi pari al valore di num_successi. Se il valore cumulativo è FALSO, verrà restituita la funzione di probabilità. Esempio 1.5 - Il lancio di una moneta può produrre come risultato solo testa o croce. La probabilità che il primo lancio sia testa è uguale a 0,5, mentre la probabilità che esattamente sei lanci su dieci siano testa è espressa dalla formula: DISTRIB.BINOM(6;10;0;5;FALSO) che risulta pari a 0.205078


Distribuzione di Poisson Esistono numerosi eventi che accadono nel tempo con cadenza del tutto irregolare. Il numero di telefonate in arrivo ad un centralino, il numero di clienti che si presentano allo sportello di un ufficio, il numero di auto che giungono ad un casello autostradale sono tutti chiari esempi di questo tipo di processi. Indichiamo con µ il numero medio di occorrenze nell’unità di tempo e supponiamo che siano soddisfatte le seguenti proprietà:

1. La probabilità di avere esattamente un’occorrenza in un intervallo di tempo dt di ampiezza trascurabile è µ dt a meno di infinitesimi di ordine superiore mentre la probabilità di avere più di un’occorrenza è trascurabile;

2. I numeri di occorrenze in intervalli temporali disgiunti sono indipendenti.

Figura 1.2 – Funzione di Probabilità e Distribuzione Cumulata di Poisson per λ = 9 Consideriamo la variabile aleatoria X che rappresenta il numero di occorrenze in un dato intervallo . Dividiamo l’intervallo in n sotto-intervalli di ampiezza t / n. La probabilità di avere esattamente una occorrenza all’interno di uno di questi sotto-intervalli è per le ipotesi fatte pari a µ t / n; per la proprietà dell’indipendenza, e ricordando la definizione della distribuzione binomiale, otteniamo che la probabilità di k occorrenze è data da (a meno di infinitesimi di ordine superiore)

( t,0 )


kn

kk

knk

nt

ntt

nkknnn

nt

nt

kn

kXP−−

−

−

+−−=

−

==

µµµµµ 11)(!

)1()1(1)( L (1.22)

Supponiamo ora che n tenda all’infinito, per le ipotesi fatte la probabilità di occorrenza all’interno di un intervallo µ dt tende a zero ma il prodotto nµ dt è pari ad una costante λ = µ t, otteniamo così la cosiddetta distribuzione di Poisson

!);()(

xexfxf

xx

XXλλ

λ−

== con x = 0, 1, 2, ... (1.23)

Per il calcolo del valor medio e della varianza ricorreremo alla funzione generatrice dei momenti che, in questo caso, è data da

∑∞

=

−− ==1

)1(

!)()(

k

tekt

X ek

eetm λλ λ (1.24)

da cui si ricava facilmente

tteX edt

tdm +−= )1()( λλ (1.25)

)1()( )1(2

2

+= +− ttteX eedt

tmd λλ λ (1.26)

e quindi

[ ] λ==dt

dmXE X )0( (1.27)

[ ] )1()0(2

22 +== λλ

dtmdXE X (1.28)

La varianza è pertanto pari a

[ ] [ ]( ) λλλλ =−+=− 222 )1(XEXE (1.29) La media e la varianza di una distribuzione di Poisson coincidono e sono entrambe pari al parametro λ. In Excel è disponibile la funzione POISSON la cui sintassi è POISSON(x;media;cumulativo)


dove x è il numero degli eventi, media è il valore numerico previsto e, come al solito, cumulativo è un valore logico che determina la forma per la distribuzione di probabilità. Se cumulativo è VERO, POISSON restituirà la probabilità cumula di Poisson indicante la probabilità che il numero degli eventi casuali sia compreso tra zero e x inclusi. Se cumulativo è FALSO, verrà restituita la funzione di probabilità di Poisson indicante la probabilità che il numero di eventi sia esattamente uguale a x. 1.5 Distribuzioni Continue Distribuzione Uniforme Diremo che una variabile aleatoria X è uniformemente distribuita nell’intervallo reale [a, b] se la sua funzione di distribuzione cumulata è data da

>

≤≤−−

<

=

bxse

bxaseabax

axsexFX

1

0

)( (1.30)

a cui corrisponde una funzione densità di probabilità data da

>

≤≤−

<

=

bxse

bxaseab

axsexf X

0

1 0

)( (1.31)

La funzione generatrice dei momenti, il valore di aspettazione e la varianza sono semplicemente calcolabili

tabeedx

abetm

atbtb

a

txX )(

1)(−−

=−

= ∫ (1.32)

∫+

=−

=b

a

abdxab

xXE2

][ (1.33)

( ) ( )122

][][222

222 ababdxab

xXEXEb

aX

−=

+

−−

=−= ∫σ (1.34)

Come vedremo la distribuzione uniforme gioca un ruolo particolarmente importante nei metodi di simulazione in quando per generare le diverse distribuzioni si parte usualmente da generatori di variabili casuali uniformi. Excel mette a disposizione la funzione CASUALE() che restituisce un numero aleatorio distribuito in maniera uniforme nell’intervallo [0, 1) (si noti che i numeri prodotti sono strettamente minori di


1). Ogni volta che si effettua un calcolo nel foglio di lavoro viene restituito un nuovo numero casuale. Se si desidera utilizzare CASUALE per creare un numero casuale senza che i numeri vengano modificati ogni volta che si esegue il calcolo della cella, è possibile immettere = CASUALE() nella barra della formula e premere F9 per convertire la formula in un numero casuale. Per creare un numero reale casuale compreso tra a e b, occorre inserire la formula seguente: = CASUALE()*(b - a) + a Sebbene le informazioni derivanti dai momenti forniscano essenziali dettagli su un insieme di dati, spesso è più comodo rappresentare questi ultimi in modo da illustrare come i valori sono distribuiti all’interno del loro dominio. Un tipo di rappresentazione che si presta particolarmente bene a questo tipo di indagine è il cosiddetto istogramma o diagramma delle frequenze relative. Il primo passo per la creazione di un istogramma è la suddivisione dell’insieme dei dati in una serie di intervalli adiacenti della stessa ampiezza (questi intervalli sono anche chiamati Intervalli delle classi). Il primo intervallo deve iniziare da un valore minore o uguale al minimo valore presente nell’insieme dei dati e l’ultimo intervallo deve, ovviamente, estendersi oltre il limite superiore dei valori dell’insieme stesso. In Excel gli istogrammi possono essere costruiti molto facilmente utilizzando la funzione Istogramma del pacchetto degli Strumenti di Analisi. Per installare questi strumenti, selezionate Componenti Aggiuntivi dal menu Strumenti; nella finestra di dialogo che viene aperta selezionate le caselle Strumenti di Analisi e Strumenti di Analisi VBA (vedi Figura 1.3). Utilizzeremo un istogramma per analizzare il generatore di numeri casuali di Excel. Se la serie dei numeri casuali è realmente distribuita uniformemente si avrà un numero pressoché uguale di numeri “casuali” in ognuno degli intervalli. Aprite una cartella di lavoro e nella cella A1 inserite la formula = CASUALE(), premete INVIO e vedrete apparire nella cella stessa un numero reale compreso fra 0 e 1. In un’altra riga è necessario introdurre gli estremi superiori degli intervalli delle classi, l’ultimo intervallo viene considerato aperto nel senso che può contenere tutti i valori che superano il limite massimo. Copiate la formula nelle prime 1000 righe del foglio di lavoro e definite gli intervalli delle classi in modo da suddividere l’intervallo [0,1] in 25 classi equispaziate e adiacenti (vedi Figura 1.4). Notate che nella definizione degli intervalli delle classi è necessario introdurre soltanto l’estremo superiore dell’intervallo, per cui nel primo cadranno tutte le occorrenze comprese fra 0 e 0.04, nel secondo quelle comprese fra 0.04 e 0.08 e così via. A questo punto selezionate Analisi Dati dal menu Strumenti (Figura 1.5). Scegliete la voce Istogramma ed inserite le informazioni richieste nella finestra dialogo risultante (Figura 1.6). In particolare è necessario specificare i blocchi di celle che contengono i dati (casella Intervallo di input) e gli estremi degli intervalli (casella Intervallo della classe). Le due caselle di input sono degli oggetti di tipo RefEdit. Il controllo RefEdit è uno strumento estremamente interessante per i programmatori in Visual Basic per Excel per cui vale la pena spendere due parole di descrizione.


Figura 1.3 – La casella dei componenti aggiuntivi

Figura 1.4 – Valori della variabile aleatoria (colonna A) e classi dell’istogramma (colonna B) Questo controllo permette all’utente di inserire un riferimento ad una o più celle utilizzando il mouse. Quando l’utente preme il pulsante alla destra del componente RefEdit, la finestra viene nascosta e ne appare un’altra che permette di accedere al foglio e selezionare una cella o un insieme di celle: i riferimenti di tale selezione vengono copiati nella nuova finestra. Al termine di questa operazione l’utente deve premere INVIO o premere di nuovo il pulsante della casella, la finestra temporanea di


input si chiude ripristinando quella originaria. All’interno del controllo RefEdit che era stato attivato compariranno così i riferimenti scelti dall’utente.

Figura 1.5 – La finestra di dialogo Analisi Dati

Figura 1.6 – La finestra di dialogo per la costruzione di un istogramma Se l’istogramma deve essere inserito nel foglio corrente attivate l’opzione Intervallo di output e digitate l’indirizzo della cella dell’angolo superiore sinistro del blocco che dovrà contenere i dati di output. Dopo aver specificato le informazioni richieste premete il pulsante OK ed Excel creerà automaticamente la tabella delle frequenze. I valori presenti in questa colonna sono numeri interi, possiamo creare una colonna delle frequenze relative dividendo questi valori per il numero totale di dati. Qualora si voglia ottenere il grafico occorre anche selezionare la casella Grafico in output (vedi FIgura 1.7). Prestate attenzione al fatto che le etichette numeriche rappresentate sotto ogni intervallo nel grafico corrispondono ai limiti superiori degli intervalli e non agli intervalli stessi. Inoltre poiché le barre verticali dell’istogramma corrispondono ad intervalli contigui può essere meglio scegliere una rappresentazione grafica in cui le barre siano disposte una accanto all’altra senza spazi intermedi. Per apportare questa modifica è sufficiente fare clic con il tasto destro del mouse su una barra verticale qualsiasi del grafico e scegliere Serie dei dati selezionati/Opzioni dal menu Formato; impostate a 0 il campo Distanza fra le barre. L’utilizzo di uno strumento di analisi come l’istogramma è


estremamente utile in vari contesti per cui può essere comodo disporre del codice VBA per gestirlo da programma. Per questo risulta di grande utilità il registratore di macro che abbiamo presentato nel primo capitolo. Nel seguente riquadro trovate il codice prodotto automaticamente dal registratore all’interno di una macro realizzata durante la creazione di un istogramma. Application.Run "ATPVBAEN.XLA!Histogram", ActiveSheet.Range("$A$1:$B$1000") _ , ActiveSheet.Range("$D$1"), ActiveSheet.Range("$B$1:$B$25"), False, False _ , True, False

Figura 1.7 – Istogramma ottenuto generando 1000 numeri casuali col generatore di Excel I dati sono quelli dell’esempio riportato in Figura 1.7. Prima di tutto osserviamo che questo comando è in realtà una chiamata ad una funzione, precisamente al metodo Run dell’oggetto Application che corrisponde ad Excel. Il metodo Run esegue una macro o richiama una funzione e può essere utilizzato per eseguire una macro scritta in linguaggio Visual Basic, segue una lista di argomenti di cui tipicamente il primo rappresenta il nome della funzione da eseguire mentre gli altri rappresentano i parametri da passare alla funzione. Selezionando gli strumenti di analisi vengono caricate due aggiunte di Excel che contengono i riferimenti necessari per poter utilizzare da VB le rispettive funzionalità. In uno di questi, ATPVBAEN.XLA, è presente la subroutine Histogram che permette la realizzazione di un istogramma; si noti che per lanciare una qualunque macro presente in un aggiunta occorre far precedere il nome della macro dal nome del file in cui essa si trova separato dal simbolo “!”. I parametri successivi sono relativi a tre intervalli di celle e a quattro variabili logiche che caratterizzano altrettante opzioni. I parametri di tipo Range hanno il seguente significato:


1. intevallo di input dei dati; 2. intervallo di output della tabella delle frequenze (notate che è sufficiente in

questo caso indicare la cella corrispondente all’angolo superiore sinistro del range);

3. intervalli delle classi. Le quattro variabili logiche corrispondono alle opzioni che è possibile selezionare dalla finestra di dialogo dell’istogramma (se le variabili vengono impostate al valore VERO la corrispondente opzione risulta selezionata):

1. opzione Pareto (istogramma ordinato); 2. opzione percentuale cumulativa; 3. opzione grafico in output; 4. indica se le prime righe dei blocchi di dati in input contengono delle etichette.

Se provate a registrare alcune macro selezionando diverse modalità di presentazione dei risultati in output troverete che quando l’intervallo di output viene posto uguale a “” (stringa vuota) la tabella delle frequenze (ed eventualmente il grafico) vengono riportati su un nuovo foglio di lavoro della stessa cartella. Se scegliete di riportare i risultati dell’istogramma su una nuova cartella di lavoro il parametro relativo all’intervallo di output viene invece omesso del tutto. Nel riquadro seguente abbiamo riportato il codice di una semplice subroutine che, a partire dal codice prodotto dal registratore di macro, permette di gestire completamente da VBA la realizzazione di un istogramma. Sub Istogramma(ClassRange As String, _ InputRange As String, _ Optional OutputRange As Variant, _ Optional flgPareto As Boolean = False, _ Optional flgCumulata As Boolean = False, _ Optional flgGrafico As Boolean = False, _ Optional flgEtichette As Boolean = False) If IsMissing(OutputRange) Then Application.Run "ATPVBAEN.XLA!Histogram", ActiveSheet.Range(InputRange) _ , , ActiveSheet.Range(ClassRange), _ fl ta, gPareto, flgCumula flgGrafico, flgEtichette ElseIf Len t ge) = 0 Then (Outpu Ran Application Run V XLA!Histogram", ActiveSheet.Range(InputRange) _ . "ATP BAEN. , "", ActiveSheet.Range(ClassRange), _ flgPareto, flgCumulata, flgGrafico, flgEtichette Else Application.Run "ATPVBAEN.XLA!Histogram", ActiveSheet.Range(InputRange) _ , ActiveSheet.Range(OutputRange), ActiveSheet.Range(ClassRange) _ , flgPareto, flgCumulata, flgGrafico, flgEtichette End If End Sub

Si noti che tutti i parametri tranne i primi due sono opzionali; inoltre, poiché dobbiamo essere in grado di controllare se il parametro OutputRange viene passato o meno alla procedura tramite la funzione IsMissing(), è necessario definirlo di tipo Variant. Per evitare che ogni operazione condotta sul foglio produca il ricalcolo dell’intera


colonna di numeri casuali, abbiate cura di deselezionare l’opzione calcolo automatico tramite la voce Opzioni del menu Strumenti (vedi Figura 1.8).

Figura 1.8 – La finestra Opzioni attivabile dal menu Strumenti Nel file “Richiami di Statistica 2.xls” l’opzione calcolo manuale viene attivata automaticamente non appena si carica il file stesso utilizzando il metodo Open dell’oggetto Workbook. Il codice relativo è riportato nel riquadro seguente e non richiede particolari commenti. Private Sub Workbook_Open() Application.Calculation = xlManual End Sub

Nel file citato il lettore troverà un’applicazione del codice sopra discusso sia per la distribuzione uniforme, sia per la distribuzione normale che adesso andremo a discutere. Distribuzione Normale Una delle funzioni più importanti, sia nella teoria sia nella pratica, è la distribuzione normale o gaussiana la cui funzione densità è data da:

( )22

2

21),;()( σ

µ

πσσµ

−−

==x

XX exfxf (1.35)

dove i parametri µ e σ sono rispettivamente la media e la deviazione standard. L’espressione della funzione generatrice dei momenti è


[ ] ( )[ ]

∫

∫

∫

∞+

∞−

−−−+

∞+

∞−

−−−−

∞+

∞−

−−

−−

=

====

dxee

dxee

dxeeeeEeeEtm

txtt

xtxt

xxttXtttX

X

22

2)2(

2

22

22

)(222)(

22

2)()(

21

21

21)(

σ

σµσµ

σ

µσµµ

σ

µµµµµ

πσ

πσ

πσ

(1.36)

D’altra parte l’ultimo integrale è uguale ad 1 (si tratta infatti dell’integrale su tutto il dominio di una funzione gaussiana) per cui possiamo scrivere:

2

22

)(tt

X etmσµ+

= (1.37) Lasciamo al lettore verificare, utilizzando la funzione generatrice dei momenti, che effettivamente i parametri µ e σ 2 sono rispettivamente la media e la varianza della distribuzione. Una variabile aleatoria viene detta distribuita secondo una normale standard se la media è 0 e la standard deviation è 1. Nel corso del testo utilizzeremo anche una notazione abbastanza diffusa tramite la quale si indica che una generica variabile aleatoria X è distribuita come una normale con media µ e varianza σ 2: X ~ ),( σµN .

Figura 1.9 – Densità di Probabilità e Distribuzione Cumulata Normale (µ = 4, σ = 0.5)


Per l’analisi di distribuzioni gaussiane Excel mette a disposizione due funzioni: una per la distribuzione generica e l’altra specializzata al caso della normale standard. DISTRIB.NORM(). Restituisce la distribuzione normale con media e distribuzione standard specificate. La sintassi è la seguente DISTRIB.NORM(x;media;dev_standard;cumulativo) dove x è il valore per il quale si desidera la distribuzione, media è la media della distribuzione, dev_standard è la deviazione standard e cumulativo è un valore logico che determina la forma assunta dalla funzione. Se cumulativo è VERO, DISTRIB.NORM restituirà la funzione di distribuzione cumulata, se è FALSO restituirà la funzione densità di probabilità. DISTRIB.NORM.ST(). Restituisce la funzione di distribuzione normale standard. La distribuzione ha una media uguale a 0 e una deviazione standard uguale a 1 Sintassi: DISTRIB.NORM.ST(z) essendo z il valore per il quale si desidera la distribuzione. Distribuzione LogNormale Sia X una variabile aleatoria con distribuzione normale, allora la variabile Y = eX definisce una variabile aleatoria con distribuita in maniera log-normale. Se la variabile X ha media µ e standard deviation σ, allora la funzione densità di probabilità di z è data da

( )2ln22

1

21)(

µσ

πσ

−−

=z

Z ez

zf (1.38)

con . La media e la varianza della variabile Z possono essere espresse in funzione dei corrispondenti momenti di X tramite le relazioni

0≥z

2

21

][σµ+

= eZE (1.39)

( ) )1(1 22222 −=−= + ωωσ µσσµ eeeZ (1.40) avendo posto . I fattori di asimmetria e curtosi sono dati rispettivamente da

1)exp( 2 ≥= σω


( ) ( 21 2/1

1 +−= ωωβ ) (1.41)

332 2342 −++= ωωωβ (1.42)

Notate che per valori di σ non nulli, sia l’asimmetria è sempre maggiore di zero e la curtosi è sempre maggiore di 3. Questo vuol dire che la distribuzione log-normale è sempre asimmetrica a destra e leptocurtica.

Figura 1.10 – Densità di Probabilità e Distribuzione Cumulata Lognormale Nel grafico riportato in Figura 1.10 sono disegnati gli andamenti delle funzioni di densità e cumulata. Per quest’ultima è stata utilizzata la funzione Excel DISTRIB.LOGNORM() la cui sintassi è la seguente DISTRIB.LOGNORM(x;media;dev_standard) dove la media e la standard deviation sono relative al logaritmo di x. La funzione Excel in questo caso restituisce solo la distribuzione cumulata. Distribuzione Gamma Infine un breve cenno alla distribuzione Gamma che è utilizzata nella teoria della struttura a termine dei tassi di interesse. Questa funzione viene utilizzata in generale nello studio di variabili aleatorie con distribuzione asimmetrica e nell’analisi delle code. L’espressione della funzione densità di probabilità è:


βαα αβ

βα /1

)(1),( x

X exf −−

Γ= (1.43)

essendo

( ) ∫∞

−−=Γ0

1 dtetz tz (1.44)

In Excel è disponibile la funzione DISTRIB.GAMMA()la cui sintassi è la seguente DISTRIB.GAMMA(x;alfa;beta;cumulativo) in cui alfa e beta sono i parametri della distribuzione e cumulativo è, come al solito, un parametro logico che indica se si desidera calcolare la densità di probabilità o la distribuzione cumulata. La distribuzione Gamma è definita solo per valori positivi della x e strettamente positivi dei parametri alfa e beta.

Figura 1.11 – Densità di Probabilità e Distribuzione Cumulata Gamma


1.6 Covarianza e Correlazione Date due variabili aleatorie X ed Y con varianza finita, si definisce covarianza la quantità definita da

][][][),( YEXEXYEYXCovXY −==σ (1.45) Se la covarianza è nulla le due variabili si dicono non correlate. Solitamente viene introdotto un coefficiente di correlazione definito come

YX

XYXY σσ

σρ = (1.46)

I cui valori massimi e minimi dipendono dal tipo di distribuzione considerata 2. Due variabili si dicono indipendenti se la funzione di distribuzione congiunta FXY(x, y)3 è fattorizzabile nel prodotto delle marginali FX(x)FY(Y). Due variabili indipendendi con varianza finita sono anche non correlate ma non è vero il viceversa. La funzione COVARIANZA() di Excel ha la seguente sintassi COVARIANZA(matrice1; matrice2) dove matrice1 e matrice2 rappresentano due intervalli di celle contenenti le serie di cui si vuole valutare la covarianza. Date due variabili aleatorie X ed Y, questa funzione utilizza lo stimatore

∑=

−−n

iYiXi yx

n 1

))((1 µµ (1.47)

essendo le medie rispettivamente di X e Y. E’ possibile inoltre calcolare direttamente il coefficiente di correlazione tramite la funzione corrispondente di Excel. 1.7 Che cos’è un Processo Stocastico Consideriamo una successione discreta di istanti di tempo nttt <<< K21 . In generale possiamo descrivere il comportamento di un sistema che evolve nel tempo in maniera 2 E’ consuetudine (errata!) ritenere che l’indice di correlazione sopra definito abbia come dominio di variazione l’intervallo chiuso [-1,1]. In realtà ciò non è vero! Il dominio dell’indice di correlazione lineare è sicuramente contenuto nel precedente intervallo ma i valori massimo e minimo che ρ può assumere dipendono dalle distribuzioni di probabilità delle due variabili aleatorie X ed Y. Una descrizione esauriente di questo risultato è al di là degli scopi di questo libro, il lettore interessato è rinviato all’articolo di Embrechts et. al. citato in bibliografica. 3 Il significato da attribuire alla funzione di distribuzione congiunta di due o più variabili è analogo a quello visto per il caso univariato . Rinviamo ai testi di statistica citati in bibliografia per la definizione delle distribuzioni marginali.

) (),( yYexXPyxFXY ≤≤=


imprevedibile tramite una corrispondente sequenza di variabili aleatorie X1, X2, ..., Xn. Parleremo in questo caso di processo stocastico discreto. Naturalmente possiamo anche definire processi stocastici nel tempo continuo sia su un dominio finito, come ad esempio [0, 1], sia su un dominio infinito, ad esempio [0, ∞). Da un punto di vista formale consideriamo uno spazio di probabilità (Ω, A , P) e un insieme non vuoto, T, i cui elementi sono gli istanti che vengono presi in considerazione. Definiamo processo stocastico una funzione di due variabili TttXX ∈= :)( RTX →Ω×: tale che

è una variabile aleatoria per ogni .) ,t()( XtX = Tt ∈ . ( ) RT →:X , . ω viene chiamata realizzazione o traiettoria del processo stocastico considerato. Ogni realizzazione in pratica non è altro che un’osservazione dell’evoluzione temporale della quantità descritta dal processo.

Figura 1.12 – Due diverse realizzazioni di un processo stocastico La simulazione di un processo stocastico è dunque legata alla possibilità di generare successioni di numeri casuali. Una domanda sorge tuttavia spontanea: come può un elaboratore, macchina totalmente deterministica, generare numeri casuali e quindi per loro natura non deterministici? La risposta è molto semplice: non può! Qualunque sequenza di numeri prodotta da un computer viene sempre generata per mezzo di qualche algoritmo per cui non si può parlare di casualità essendo la sequenza predeterminata. Come vedremo nei prossimi capitoli, in compenso, con un computer si possono generare sequenze di numeri che sembrano casuali.

2 La Simulazione Monte Carlo

La ricerca nel campo della fisica ha dimostrato in modo lampante che, almeno nella stragrande maggioranza dei fenomeni naturali la cui regolarità e invariabilità ha condotto alla formulazione del postulato della causalità generale, è il caso la radice comune della stretta rispondenza alle leggi riscontrata. Erwin Schroedinger “Che cos’è una legge fisica”, 1922

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. Per il momento assumiamo che il lettore comprenda, per lo meno a livello intuitivo, che cosa sia un numero casuale. Più avanti avremo modo di ritornare sulla definizione di numero casuale e, come vedremo, questo porterà a delle questioni tutt’altro che banali. Storicamente la prima applicazione su larga scala di metodi computazionali basati sull'uso di numeri casuali furono gli studi di diffusione ed assorbimento di neutroni (una particella elementare che compone il nucleo atomico) in vari tipi di materiale nell'ambito del progetto per la realizzazione della prima bomba atomica. Prima di passare ad alcuni esempi pratici, cerchiamo di dare una definizione formalmente più precisa. 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. E’ opportuno soffermarsi sulla natura dei problemi che è possibile affrontare con il metodo Monte Carlo. Semplificando al massimo la questione, è possibile dividere tali problemi in due grandi categorie

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.

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.

Nel caso di un problema del primo tipo è possibile stabilire una corrispondenza diretta fra popolazione simulata e popolazione studiata. Tali applicazioni sono dette


simulazioni dirette. Nel secondo caso la formulazione del problema può essere tale da rendere l’applicazione del metodo Monte Carlo poco evidente soprattutto se il problema in esame è intrinsecamente analitico. Vediamo subito un esempio pratico. Con riferimento alla fig. 2.1 consideriamo dapprima dei problemi di natura probabilistica (la valutazione del premio per un’opzione di tipo europeo scritta su un titolo azionario che non paga dividendi e la determinazione del VaR) ed un problema di natura prettamente deterministica come la determinazione dell’area racchiusa da una figura piana, ad esempio un cerchio. Cominciamo dal primo problema. 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

[ ]TrT efEf −= ˆ (2.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 (2.1) diventa

[ TTr fEef ˆ−= ] (2.2)

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; l’assunzione più 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

SdzSdtdS σµ += (2.3) dove µ è il tasso di rendimento atteso e σ è la volatilità del prezzo dell’azione. 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 (2.3) è data da

tStSS ∆+∆=∆ εσµ (2.4) dove ∆S è la variazione di prezzo nell’intervallo ∆t ed ε è un numero casuale estratto da una distribuzione normale (supponiamo per il momento di avere una qualche procedura che ci permetta di generare numeri casuali con distribuzione di probabilità assegnata). I valori di µ e σ possono essere stimati da un’analisi delle serie storiche dei prezzi per il titolo sottostante l’opzione. 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

)0,max( KSf TT −= (2.5)

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 che compare nella (2.2). Lo stesso problema può essere risolto con tecniche analitiche sia tipo numerico, come l’approccio basato sull’uso di alberi binomiali, sia ricercando, laddove le ipotesi non siano troppo restrittive per un uso pratico, una soluzione esatta come avviene nel caso dell’analisi svolta da Black e Scholes che in un lavoro ormai storico hanno ricavato una formula chiusa per il pricing di opzioni europee. Restando all’interno del primo tipo di problematiche vediamo come formalizzare il problema del VaR di un portafoglio in maniera adatta all’utilizzo di tecniche Monte Carlo. Supponiamo che il valore corrente del portafoglio sia V(t), sia ∆t il periodo di tempo su cui intendiamo stimare il VaR e infine V )( tt ∆+ indichi il valore del portafoglio al termine di questo periodo. La perdita di valore del portafoglio nel periodo considerato è pari a e il VaR, x)()( tVtVVL −∆+=∆= p, associato ad un certo livello di probabilità p è definito dalla relazione

pxLP p => )( (2.6)

PROBLEMA FORMULAZIONE SOLUZIONE

VALUTAZIONE DELPREZZO DI UN’OPZIONE

(PROBABILISTICO)

CALCOLO DELL’AREA DIUN CERCHIO

(DETERMINISTICO)

TIPO A(STATISTICO)

TIPO B(ANALITICO)

ESTRATTA DALCAMPIONE

(APPROSSIMATA)

NUMERICA

ESATTA

Simulazionediretta

geometria

Numericasuali

Metodi divalutazione

“Crude”Monte Carlo

Ad esempio:Formula di Black & Scholes

S = πR2

Ad esempio:Alberi Binomiali

Integrazione numerica

Figura 2.1


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.

Figura 2.2 – Il calcolo di Pi greco.

Vediamo infine il problema di tipo geometrico che, al di là dell’ovvia banalità della soluzione analitica, non sembra presentare alcuna relazione con l’utilizzo di sequenze di numeri casuali. Consideriamo per semplicità un cerchio di raggio unitario di area pari a π. Poniamoci ora sotto una prospettiva leggermente diversa da quella usuale ed osserviamo che l’area di questo cerchio può essere misurata stimando la probabilità che un punto scelto a caso all’interno di un quadrato di lato pari a 2 cada all’interno del cerchio inscritto in esso. Con riferimento alla figura 2.2 (dove per semplicità ci siamo limitati alla considerazione del primo quadrante in un sistema di assi coordinati con origine nel centro del cerchio) la procedura per la determinazione del valore di π può essere descritta come segue:

1. si estraggono due numeri casuali, indipendenti, uniformemente distribuiti nell’intervallo [0,1]. Questa coppia di numeri rappresenta le coordinate di un punto del piano appartenente al quadrato di lato unitario;

2. si determina la distanza del punto così generato dall’origine degli assi, se tale

distanza è minore o uguale a 1 il punto cade all’interno della superficie del cerchio altrimenti no;

3. si ripete la procedura sopra descritta un numero N di volte e si misura il rapporto

fra il numero dei casi in cui il punto cade all’interno del cerchio, Nacc, e N stesso, questo rapporto rappresenta una stima della probabilità che un punto con coordinate scelte a caso cada all’interno del cerchio, tale probabilità è chiaramente uguale al rapporto fra l’area del quadrante del cerchio e quella del quadrato unitario quindi si ha

NNacc4=π


Non c’è bisogno di dire che nessuno di noi stimerebbe l’area di un cerchio con questa procedura, tuttavia problemi di natura analitica meno banali di quello presentato possono trarre grande giovamento dalla riformulazione in termini di processi aleatori al punto che, contrariamente a quello che si potrebbe pensare, tecniche di simulazione per la soluzione di problemi deterministici di tipo fisico e matematico erano già state sviluppate molto prima degli anni 40. Storicamente il primo esempio di simulazione Monte Carlo ante litteram risale addirittura al 1777 anno in cui il matematico e naturalista francese Buffon nel suo Essai d’arithmetique morale descrisse un esperimento di natura aleatoria che gli consentì di stimare, con una precisione per l’epoca accettabile, il valore di π. 2.1 Simulazione Monte Carlo e Integrazione Numerica Esiste una connessione formale fra l’utilizzo del metodo Monte Carlo e il concetto di integrazione di una funzione. Prima di tutto osserviamo come i problemi discussi nel paragrafo precedente possono essere ricondotti entrambi al calcolo di integrali. Il caso relativo all’area del cerchio è evidente

∫ −=1

0

214 dxxcerchio del Area (2.7)

Il prezzo di un’opzione come abbiamo visto non è altro che il valore scontato del valore di aspettazione del prezzo alla scadenza, il fattore di rischio sottostante (il prezzo dell’azione) risulta distribuito in accordo ad una distribuzione log-normale perciò abbiamo:

SdeSSf

TetSf TTrSS

Tr′

′′

= ∫∞

−+′−−

0

2/)])2/1([)/(log( 222)(2

),( σσ

πσ (2.8)

Infine, per quanto riguarda il calcolo del Value-at-Risk, osserviamo che indicata con I(L>x) la funzione indicatrice definita come

>

=>altrimenti

x L sexLI

01

)( (2.9)

e indicando con f(z) la densità di probabilità congiunta dei fattori di rischio, possiamo scrivere

])([)()()( ∫ >=>=> xLIEdzzfxLIxLP (2.10) Più in generale possiamo affermare che ogni estrazione di un campione di numeri casuali può essere utilizzata come stimatore di un integrale. Come esempio concreto


consideriamo il caso relativo all’integrazione di una funzione di una variabile

reale ; tramite un opportuno cambio di variabile possiamo sempre ricondurci al

caso più semplice in cui l’intervallo di integrazione è compreso fra 0 ed 1:

∫=b

adyygI )(

∫=1

0)( dxxfI (2.11)

Il punto fondamentale del nostro ragionamento consiste nel riconoscere che l’espressione scritta sopra rappresenta anche 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~ (2.12)

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)(1var)~var( σ

=−=

=

= ∫∑∑

== (2.13)

come si vede 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 n/1 . 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. E’ comunque conveniente a questo punto formalizzare il problema in modo più rigoroso. Consideriamo il seguente integrale

∫=D

xdxfI rr)( (2.14)

dove abbiamo indicato col simbolo x

r il vettore multidimensionale e con l’elemento di volume infinitesimo in n dimensioni , f(x) è una funzione

integrabile secondo Lebesgue definita sul dominio D. Ogniqualvolta l’integrale (2.14) esiste, il suo valore è uguale a quello dell’integrale di Lebesgue-Stieltjes

),,,( 21 nxxx K

xd r

ndxdxdx K21

∫=Z

LS zdFzI )()( rrκ (2.15)


dove ),( Dxxf ∈rr e D determinano la forma assunta dalla nuova regione di

integrazione Z, la funzione ),( Zzz ∈rrκ è misurabile e la funzione non decrescente

Zz ∈),( zF rr per ogni sottoinsieme misurabile di Z, con A e B tali che , soddisfa la relazione

ZBA ⊆⊆

1)()()(0 =≤≤≤ ZFBFAF

F(z) può quindi essere considerata una funzione di distribuzione e questo permette di interpretare l’integrale (2.15) come il valore di aspettazione di )(zrκ dove z è una variabile aleatoria n-dimensionale la cui funzione di distribuzione è proprio F(z). Questo ci permette di affermare in modo del tutto generale che la quantità

)(1

1∑=

=n

ii

nLS z

nI

rκ (2.16)

dove è un campionamento di n valori indipendenti di Z estratti dalla distribuzione F(Z), è uno stimatore privo di errori sistematici dell’integrale (2.15). Vale la pena notare che l’integrale di Lebesgue-Stiltjes associato all’integrale di Lebesgue non è unico. Infatti, se consideriamo una nuova funzione F* definita in modo che dF(z)/dF*(z) esista per ogni z appartenente a D, e poniamo

izr

)(*/)()()(* zdFzdFzz κκ = otteniamo

∫∫ =ZZ

zdFzzdFz )(*)(*)()( rrrr κκ (2.17)

Come vedremo più avanti questa opportunità risulta di grande aiuto nella ricerca di metodologie per la riduzione della varianza.

3 Generatori di Numeri Casuali

Random numbers should not be generated with a method chosen at random. Some theory should be used. Donald Knuth, “L’arte della Programmazione”, vol. 2

3.1 Tre tipi di generatori 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. In generale il loro uso è estremamente difficoltoso, le serie di numeri veramente casuali sono quasi sempre composte da un numero di dati insufficiente per effettuare simulazioni attendibili pertanto non le discuteremo ulteriormente.

• 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. Questi sono i numeri più utilizzati per le simulazioni Monte Carlo e verranno discussi più approfonditamente nel seguito del capitolo.

• 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. Tali sequenze non risultano appropriate per qualsiasi applicazione del metodo Monte Carlo, tuttavia laddove il loro uso è consigliato possono portare a dei miglioramenti impressionati nelle performance dei metodi numerici che ne fanno uso.

Le sequenze di numeri, sia pseudo che quasi casuali, sono caratterizzate dall’essere periodiche. Dopo un certo intervallo di numeri prodotti, detto periodo, i generatori iniziano a riproporre gli stessi numeri casuali della sequenza appena conclusa. Chiaramente una delle proprietà più importanti di qualunque tipo di generatore è quella di avere un periodo che risulti il più grande possibile. Discuteremo questa ed altre caratteristiche facendo riferimento ad un particolare metodo di generazione di numeri


casuali il cosiddetto generatore lineare congruente che ha il pregio di essere estremamente semplice e computazionalmente efficiente. 3.2 Il generatore Lineare Congruente Virtualmente tutti i generatori di numeri pseudo casuali impiegati in pratica sono basati sul generatore lineare congruente, da un valore iniziale, J0, si genera una sequenza di numeri basati sulla formula iterativa

m mod )caJ(J 1ii += − (3.1)

I parametri a, c e m determinano la qualità del generatore; a viene detto moltiplicatore, c incremento e m è il cosiddetto modulo. La formula appena descritta genera numeri interi compresi fra 0 e m. Usualmente si utilizzano generatori di numeri casuali uniformemente distribuiti fra 0 ed 1, per questo è sufficiente scegliere la sequenza

. Dalla formula risulta evidente che la sequenza di numeri casuali si ripeterà dopo un ciclo che, al più, potrà essere di lunghezza m. In generale se indichiamo con L la lunghezza di parola (word length) di un computer espressa in numero di bit, il massimo intero rappresentabile è pari a

mJU ii /=

12 −L , poiché per motivi pratici un bit è spesso perso (in quanto riservato per indicare il segno del numero), la scelta usuale è quella di porre . 11 −−L2=m Come abbiamo già accennato, il principale vantaggio del generatore lineare congruente va ricercato nella sua efficienza. Esso è molto veloce richiedendo pochissime operazioni per chiamata, questo è uno dei principali motivi che hanno finito per renderlo di uso universale. Fra i punti deboli di questo metodo il maggiore è rappresentato senz’altro dalla presenza di correlazione sequenziale fra i numeri generati; vediamo di capire di cosa si tratta. Come risulta chiaro dalla (3.1) ciascun numero prodotto dal generatore lineare congruente è completamente determinato dal numero generato al passo precedente. Questo fa sì che i valori dei numeri prodotti siano fortemente correlati fra loro. Ad esempio a numeri piccoli seguiranno preferibilmente numeri piccoli con una probabilità superiore a quella che dovremmo attenderci se i numeri fossero veramente casuali. La correlazione fra due numeri diminuisce all’aumentare della distanza che li separa espressa in termini di passo di generazione. Uno degli effetti più evidenti di questa correlazione consiste nel seguente fatto. Si generino n numeri casuali e si formino n/2 coppie di numeri scegliendo valori generati successivamente uno all’altro

. Consideriamo poi queste coppie come coordinate di punti su un piano. Se i numeri fossero veramente casuali dovremo aspettarci di veder riempire uniformemente il quadrato unitario. Invece quello che si osserva è che tali punti non si distribuiscono affatto in maniera uniforme ma, al contrario, tendono ad allinearsi lungo segmenti di retta.

),( 1+nn UU

Esempio 3.1 – Supponiamo di generare un migliaio di punti nel quadrato di lato unitario utilizzando come motore il generatore lineare congruente (3.1) con i seguenti


valori per i parametri: a = 1229, b = 1, c = 2048 e valore iniziale della sequenza pari a 0. Questi numeri sono tali da dar luogo ad un pessimo generatore di numeri casuali come si può facilmente verificare riportando in un grafico i punti così ottenuti (Fig.3 a). La correlazione sequenziale può essere facilmente rimossa con tecniche di mescolamento (shuffling), il numero prodotto al passo j-esimo non costituisce la risposta j-esima del generatore ma viene utilizzato per fornire ad un passo successivo un’uscita scelta in maniera casuale fra una serie di numeri precedentemente generati. Come si può vedere in Fig.3 b con questa semplice procedura il risultato migliora decisamente

a)

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

b)

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

Figura 3.1 – Risultati del generatore lineare congruente con a = 1229, b = 1, c = 2048 e valore iniziale nullo. Il generatore non utilizza la procedura di shuffling nel primo caso 3.3 Le funzioni Rnd() e Randomize() 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.


Esempio 3.2 - Il ciclo seguente For i = 1 To 5 Cells(i, 1).Value = Rnd(-1) Next

produce in output la serie costante 0.22400701 0.22400701 0.22400701 0.22400701 0.22400701 Se utilizzate Rnd() con un valore positivo For i = 1 To 5 Cells(i, 1).Value = Rnd(1) Next

otterrete invece 0.035845816 0.086352348 0.164263904 0.17973578 0.573200524 Se a questo punto richiamate Rnd(0) otterrete in output sempre l’ultimo numero generato 0.573200524 0.573200524 0.573200524 0.573200524 0.573200524 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. Per ripetere sequenze di numeri casuali, è sufficiente richiamare Rnd() con un argomento negativo immediatamente prima dell'utilizzo di Randomize con un argomento numerico. L'utilizzo di Randomize con lo stesso valore di numero non determinerà la ripetizione della sequenza precedente. Se Randomize non viene utilizzata, quando la funzione Rnd (senza argomenti) viene chiamata per la prima volta, utilizza come base lo stesso numero. Per le chiamate successive la funzione utilizzerà l'ultimo numero generato. 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) In questo contesto, limitesup indica il numero maggiore presente nell'intervallo, mentre limiteinf indica il numero minore. 3.4 Test statistici sui generatori di numeri casuali Nell’esempio 3.1 abbiamo visto come la presenza di correlazione sequenziale del generatore lineare congruente possa alterare in maniera drammatica i risultati di un esperimento che faccia uso di numeri casuali. Di fatto, come abbiamo già osservato, il generatore lineare congruente, non produce numeri veramente casuali. Diventa quindi essenziale disporre di test statistici che ci permettano di giudicare la qualità di una serie di numeri quasi-random anche al fine di poter paragonare le serie prodotte da diversi generatori per selezionare quello che meglio soddisfa le nostre esigenze. E’ opportuno ricordare che in tutte le considerazioni relative alla scelta del generatore di numeri casuali da utilizzare si dovrà sempre tenere conto anche dell’efficienza numerica del generatore stesso a causa dell’elevato numero di variabili aleatorie che è necessario generare nelle applicazioni pratiche del metodo Monte Carlo.

Esistono moltissimi test statistici che possono essere utilizzati per controllare la qualità di un generatore; tra i più usati vogliamo ricordare il test di χ2 e il test di Kolmogorov-Smirnov. Entrambi si basano su una comparazione fra la funzione di distribuzione cumulata estratta da un campione dei numeri casuali prodotti dal generatore con la funzione di distribuzione cumulata ricavata dalla funzione di densità di probabilità in accordo alla quale il generatore ha prodotto i numeri stessi. Vediamo un esempio concreto 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. 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. A questo punto calcoliamo il valore di χ2 e verifichiamo se il generatore è accettabile utilizzando i valori riportati in Tabella 1. La tabella riporta, per ciascun valore dei gradi di libertà k, il valore per cui la probabilità

2,1 kαχ −

( )2k,1

2αχχ −<P è 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. Una volta trovata la probabilità corrispondente al valore di χ2 misurato possiamo applicare il seguente criterio Knuth (1981):

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% ⇒ leggermente 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.

Come esempio pratico abbiamo applicato la procedura appena descritta al generatore lineare congruente caratterizzato dai parametri definiti poco sopra, i valori di χ2 misurati su tre serie indipendenti di 1000 numeri sono risultati pari a:

94.2 95.6 88.8

Utilizzando la formula asintotica con k = 99, troviamo i seguenti valori di χ2

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

61.66594 77.049554 81.38991 98.333333 108.13518 123.20325 134.73861


Come si vede i valori trovati cadono tutti nell’intervallo compreso fra il 25% e il 50% di probabilità, il nostro generatore ha quindi superato il test. La stessa prova eseguita scegliendo i parametri dell’esempio 1 ha dato invece i seguenti risultati

57.6 62.9 60.1

In questo caso tutti i valori sono tali che la probabilità corrispondente risulta inferiore all’1%, un’ulteriore conferma del fatto che il set di parametri che caratterizza questo generatore corrisponde a numeri molto correlati fra loro e quindi molto lontani dall’emulare un insieme di numeri veramente casuali.

Tabella 3.1 – 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 32

322 2 −++≈ pp xxkk

xp -2.33 -1.64 -0.675 0.00 0.675 1.64 2.33

4 Generazione di Distribuzioni non Uniformi

Non ho ancora trovato un problema che, per quanto fosse complicato, a considerarlo nel modo giusto non diventasse ancora più complicato. Poul Anderson

4.1 Il “Transformation Method” Nei precedenti paragrafi abbiamo visto come generare sequenze di numeri casuali uniformemente distribuite. Utilizzando un generatore di questo tipo la probabilità di ottenere un risultato con un valore compreso fra x e x+dx è data da

<<

= altrimenti

x per dxdxxp

010

)( (4.1)

Supponiamo ora di generare una serie di numeri casuali uniformemente distribuiti e di calcolarne una qualche funzione definita y=y(x). La distribuzione di probabilità di y, che indicheremo con g(y), può essere ricavata ricordando che

dxxpdyyg )()( = (4.2) da cui

dydxxpyg )()( = (4.3)

A scopo di esempio supponiamo che )ln()( xxy −= con x variabile distribuita uniformemente, si ha

dyedydydxdyyg y−==)( (4.4)

y risulta distribuita esponenzialmente. Il problema appena discusso può anche essere posto in maniera diversa. Supponiamo infatti di conoscere la distribuzione di probabilità della variabile y e cerchiamo di determinare quale deve essere la relazione funzionale che lega y a x. In altre parole, supponendo sempre di avere una distribuzione p(x) uniforme, il problema che dobbiamo affrontare è quello di risolvere l’equazione differenziale

dydxyg =)( (4.5)

la cui soluzione è

∫ == )()( yFdyygx (4.6) da cui


)(1 xFy −= (4.7)

Il metodo appena esposto va sotto il nome di Transformation Method. Seguendo questa procedura possiamo definire un semplice algoritmo per la generazione di variabili casuali da distribuzioni diverse da quella uniforme. L’algoritmo si riassume nei seguenti passi:

1. Si genera un numero casuale x distribuito uniformemente fra 0 ed 1; 2. Si determina, qualora risulti possibile, l’inversa della funzione di ripartizione

della distribuzione che si desidera replicare; 3. Si determina il valore di tale funzione in x; 4. Si torna al punto 1 e si ripete il ciclo n volte.

I numeri così ottenuti risulteranno distribuiti in accordo con la funzione desiderata. 4.1.1 Utilizzo delle funzioni del foglio di lavoro di Microsoft Excel in Visual

Basic Una caratteristica molto importante offerta dal Visual Basic è la possibilità di utilizzare direttamente nel codice funzioni del foglio di lavoro di Microsoft Excel. Per un elenco completo delle funzioni del foglio di lavoro disponibili si rinvia all’elenco presente nella guida in linea 4. In Visual Basic le funzioni del foglio di lavoro di Microsoft Excel sono disponibili tramite l'oggetto WorksheetFunction. Nella seguente routine Sub viene utilizzata la funzione del foglio di lavoro Min per determinare il valore minimo in un intervallo di celle. La variabile myRange viene innanzitutto dichiarata come oggetto Range, quindi viene associata all'intervallo A1:C10 di Sheet1. Ad un'altra variabile, answer, viene assegnato il risultato della funzione Min su myRange. Infine, il valore di answer viene visualizzato in una finestra di messaggio. Sub UseFunction() Dim myRange As Range Set myRange = Worksheets("Sheet1").Range("A1:C10") answer = Application.WorksheetFunction.Min(myRange) MsgBox answer End Sub

Se viene utilizzata una funzione del foglio di lavoro che richiede un riferimento di intervallo quale argomento, è necessario specificare un oggetto Range. Ad esempio, è possibile utilizzare la funzione del foglio di lavoro Match per eseguire ricerche in un 4 Alcune funzioni del foglio di lavoro sono superflue in Visual Basic. La funzione Concatenate, ad esempio, non è necessaria poiché in Visual Basic è possibile utilizzare l'operatore & per concatenare più valori di testo.


intervallo di celle. Occorre quindi immettere in una cella del foglio di lavoro una formula quale = CONFRONTA(9,A1:A10,0). Tuttavia, in una routine Visual Basic è possibile specificare un oggetto Range per ottenere lo stesso risultato. Sub FindFirst() myVar = Application.WorksheetFunction.Match(9, Worksheets(1).Range("A1:A10"), 0) MsgBox myVar End Sub

E’ importante notare che non è detto che una funzione Visual Basic con lo stesso nome di una funzione Microsoft Excel operi in modo identico. Ad esempio, Application.WorksheetFunction.Log e Log restituiranno valori diversi. E’ possibile anche immettere una funzione del foglio di lavoro in una cella. Per fare questo è necessario definire la funzione come valore della proprietà Formula del corrispondente oggetto Range. Nell'esempio seguente la funzione del foglio di lavoro CASUALE, che genera un numero casuale, viene assegnata alla proprietà Formula dell'intervallo A1:B3 di Sheet1 nella cartella di lavoro attiva. Sub InsertFormula() Worksheets("Sheet1").Range("A1:B3").Formula = "=RAND()" End Sub

4.1.2 Come Simulare una variabile distribuita Normalmente Come esempio del metodo descritto nel paragrafo precedente proviamo a generare una sequenza di numeri casuali con distribuzione normale. Per raggiungere il nostro scopo abbiamo bisogno della funzione inversa della distribuzione cumulata che, fortunatamente, Excel mette già a nostra disposizione. Le funzioni che potete inserire direttamente nei fogli di calcolo sono in realtà due: 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.


INV.NORM.ST(). Restituisce l'inversa della distribuzione normale standard cumulativa. La distribuzione ha una media uguale a 0 e una deviazione standard uguale a 1 Sintassi: INV.NORM.ST(probabilità) dove probabilità è la probabilità corrispondente alla distribuzione normale. Anche questa funzione utilizza un metodo di calcolo iterativo analogo a quello impiegato in INV.NORM. Le due funzioni sono richiamabili da VBA tramite l’oggetto WorksheetFunction con una sintassi leggermente diversa da quella utilizzata da Excel nei suoi fogli di calcolo. Se desiderate generare n variabili con distribuzione normale standard potete utilizzare il seguente codice Randomize For i = 1 To n Cells(i ,1).Value = Application.WorksheetFunction.NormSInv(Rnd()) Next

Si noti che in questo caso l’unico argomento richiesto dalla funzione è un numero casuale uniformemente distribuito fra 0 ed 1. Nel caso in cui desideriate generare variabili con distribuzione normale caratterizzate da una media M e standard deviation Sigma è sufficiente utilizzare la funzione NormInv() del tutto analoga a INV.NORM(); un esempio di codice in questo caso è Randomize For i = 1 To n Cells(i ,1).Value = Application.WorksheetFunction.NormInv(Rnd(), M,Sigma) Next

Nel file “Richiami di Statistica 2.xls” il lettore troverà un esempio di generazione di variabili aleatorie con distribuzione normale. In linea di principio il metodo appena visto si presta alla generazione di variabili aleatorie distribuite in accordo ad una qualunque funzione di densità di probabilità purché questa sia continua. In pratica tuttavia un’applicazione diretta di tale metodo può risultare molto onerosa sotto il profilo del calcolo. Questo accade in particolare per la distribuzione gaussiana per la quale non risulta possibile esprimere in forma analitica l’inversa della funzione di ripartizione. In casi come questo l’integrale della funzione deve essere calcolato numericamente con evidente perdita di tempo. Il metodo di Box-Muller per la generazione di variabili aleatorie distribuite normalmente evita questo problema. Questa procedura è basata sull’osservazione che se U1 e U2 sono due variabili aleatorie indipendenti distribuite uniformemente fra 0 ed 1 allora le variabili G1 e G2 definite dalle relazioni


)2()ln(2

)2()ln(2

212

211

UsinUG

UcosUG

π

π

−=

−= (4.8)

sono due variabili indipendenti distribuite normalmente.

Esempio 4.1 – Utilizziamo il generatore lineare congruente con i valori dei parametri definiti nell’esempio precedente per generare N coppie di punti nel piano con coordinate distribuite normalmente applicando la trasformazione di Box-Muller. Quello che si trova è abbastanza sorprendente e sicuramente diverso da quello che ci aspettavamo. Le coppie di numeri sono correlate clamorosamente una all’altra. Questo si verifica facilmente osservando che (dalla (3.1))

( )[ ] mm caUU ii /mod1 +=+ ricordando la definizione dell’operazione di “modulo” possiamo scrivere

( ) ( ) kUkmcU

macaU INT mcaU

mU iiiii −+=−+=+−+=+ γα1

1

essendo k un numero intero. A causa della periodicità delle funzioni trigonometriche, si ha:

( ) ( )[ ]γαππ +=+ ii UU 2cos2cos 1 ( ) ( )[ ]γαππ +=+ ii UsinUsin 22 1

Quindi in questo caso particolare le relazioni di Box-Muller diventano

( )( ))2()ln(2

)2()ln(2

112

111

γαπ

γαπ

+−=

+−=

UsinUG

UcosUG

Entrambe le nuove variabili dipendono solo dalla prima variabile U1 generata inoltre il luogo dei punti descritto dalle equazioni scritte sopra è una spirale logaritmica, Geweke (1995b). 4.2 La simulazione di un processo di Wiener Il processo di Wiener prende il nome dal fisico-matematico tedesco Norbert Wiener ed è una descrizione matematica del fenomeno fisico noto come moto browniano ovverosia del moto casuale a cui è sottoposta una piccola particella di materiale sospesa in una soluzione liquida. Il moto browniano trova la sua spiegazione nella teoria cinetica della materia in quanto il moto della particella altro non è che il risultato dei numerosissimi urti che questa riceve dalle molecole del liquido circostante. Nella letteratura i termini


processo di Wiener e moto browniano sono pertanto usati come sinonimi anche se, a rigore, il primo individua un modello matematico mentre il secondo andrebbe riservato per indicare il processo fisico corrispondente. Nella matematica finanziaria i processi di Wiener svolgono un ruolo estremamente importante per la descrizione della componente stocastica della dinamica di prezzi.

-3

-2

-1

0

1

2

3

-3 -2 -1 0 1 2 3

Figura 4.1 – Risultati dell’applicazione del metodo di Box-Muller alle coppie di numeri casuali prodotti dal generatore lineare congruente dell’esempio 1 senza shuffling.

Definiamo processo di Wiener standard 0),( ≥= ttWW un processo stocastico gaussiano nel tempo continuo con incrementi indipendenti tale che

• con probabilità 1; 0)0( =W• ; [ ] 0)( =tWE• per ogni [ ] stsWtW −=− )()(2σ ts ≤≤0 ;

Nella Figura 4.2 sono riportate alcune realizzazioni di un processo di Wiener realizzate utilizzando la procedura riportata nel riquadro seguente Sub Wiener() Dim i As Integer Dim nPunti As Integer Dim t0 As Single Dim tf As Single Dim t As Single Dim Delta As Single Dim x As Single


Dim wt As Single t0 = 0 tf = Range("t_finale").Value nPunti = Range("NrPunti").Value Delta = (tf - t0) / nPunti t = t0 - Delta wt = 0 Randomize (Timer) For i = 1 To nPunti t = t + Delta x = Application.WorksheetFunction.NormSInv(Rnd()) If x < 10 And x > -10 Then wt = wt + x * Sqr ta) (Del Cells(i + 1, 1).Value = t Cells(i + 1, 2).Value = wt End If Next Calculate Cells(k + 5, 16).Value = Cells(5, 15).Value Cells(5, 17).Value = k Cells(k + 5, 19).Value = Cells(5, 18).Value Cells(k + 5, 20).Value = Range("B202").Value k = k + 1 End Sub

Figura 4.2 – Alcune realizzazioni di un processo di Wiener


4.3 Il “Rejection Method” Un altro metodo estremamente diffuso per generare variabili aleatorie a distribuzione non uniforme è il cosiddetto Rejection Method. Esso si basa sul seguente teorema enunciato da von Neumann nel 1951, Fishman, (1996)

Teorema Sia f(z) con a bz ≤≤ una funzione di distribuzione di probabilità, supponiamo che sia possibile trovare una seconda funzione h(z) > 0 e con integrale finito tale che

f(z)=cg(z)h(z) essendo 1)(0 ≤≤ zg . La costante c può essere scelta in modo da normalizzare l’integrale della funzione f(z) per cui possiamo assumere senza perdere in generalità

∫=

b

adzzhzg

c

)()(

1

Sia infine Z una variabile aleatoria con funzione di distribuzione h(z) e U indichi un campione estratto da una distribuzione uniforme U(0,1), allora se noi selezioniamo le estrazioni di Z tali che U )(Zg≤ i valori così ottenuti saranno distribuiti secondo f(z).

Per dimostrare il teorema appena enunciato osserviamo prima di tutto che le variabili aleatorie U e Z hanno una funzione di distribuzione di probabilità congiunta data da

bza u zhzuf ZU ≤≤≤≤= 10)(),(,

Da questa possiamo ricavare la funzione di distribuzione di probabilità di Z condizionale a U )(Zg≤

( ))]([

)()()]([

)(

)]([

),()(|

)(

0

)(

0,

ZgUprzgzh

ZgUpr

duzh

ZgUpr

duzufZgUzh

zgzg

ZU

Z ≤=

≤=

≤=≤

∫∫

D’altra parte abbiamo

∫ ∫∫ ===≤)(

0,

1)()(),()]([zg b

aZU

b

a cdzzgzhzuf dudzZgUpr

da cui

( ) )()()()(| zfzgzh cZgUzhZ ==≤


come volevamo dimostrare. Da un punto di vista implementativo possiamo sfruttare questo teorema per definire un semplice algoritmo; per semplicità definiamo la funzione e(z) pari a c h(z). Cominciamo con l’osservare che dalla definizione di g(z) si ricava immediatamente che e(z) è una funzione maggiorante della distribuzione che intendiamo campionare f(z), ci riferiremo ad e(z) con il nome di funzione di paragone. In termini di e(z) la condizione di accettazione può essere riformulata come )(ZgU ≤ )()( ZgZe ≤ . A questo punto possiamo descrivere un algoritmo di generazione di numeri casuali con distribuzione assegnata secondo il rejection method:

1. Si genera un numero casuale y’ con distribuzione uniforme fra 0 e A essendo

∫=b

adzzeA )(

2. Si ricava il valore x = E-1(y’) essendo . I punti x saranno

distribuiti in accordo alla funzione e(z) (in pratica stiamo utilizzando il transformation method discusso prima);

∫=x

adzzexE )()(

3. Si genera un secondo numero casuale, y”, uniformemente distribuito fra 0 ed e(x);

4. Se risulta soddisfatta la condizione )(xfy ≤′′ accettiamo x; 5. Si ripete la procedura n volte

Un punto cruciale, da cui dipende il successo dell’applicazione di questo metodo, è rappresentato dalla scelta della funzione di paragone e(z). Tale funzione deve essere sempre maggiore o uguale alla funzione di distribuzione che intendiamo campionare e tale che sia facilmente computabile l’inversa della corrispondente funzione di ripartizione. 4.4 Generazione di variabili con correlazione assegnata Indichiamo con un vettore di variabili aleatorie indipendenti ciascuna delle quali distribuita secondo una normale standard, la matrice di varianza-covarianza di

),...,,( 21 nxxxx =r

xr

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

, non più indipendenti bensì dotato di matrice di varianza-covarianza assegnata Σ. Il nuovo insieme di variabili aleatorie

),...,,( 21 nyyyy =r

yr può essere ricercato come

combinazione lineare delle variabili indipendenti xr , cioè si pone

xAyrr

= (4.9)


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

Σ=TAA (4.10) La soluzione della (4.10) 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 la (4.10) 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

AAA

A

L

MOMM

L

K

21

2221

11000

(4.11)

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

∑−

=−=

1

1

2i

kikiiii aa σ (4.12)

∑−

=−=

1

1

1 i

kjkikij

iiij aa

aa σ (4.13)

Ad esempio per il caso semplice di due variabili con matrice di varianza-covarianza definita da

=Σ 2

221

2121

σρσσρσσσ

troviamo facilmente

−= 2

22

1

1

0

ρσρσ

σA

Esempio 4.2 – Consideriamo il modello di Fong e Vasicek della struttura per scadenza dei tassi di interesse. Si tratta di un modello a due fattori rappresentati rispettivamente dal tasso spot e dalla volatilità. La dinamica dei fattori è descritta dal sistema di equazioni differenziali stocastiche

1)( dZdtrkdr νµ +−=


2)( dZdtd νξναγν +−= Questo modello risulta particolarmente interessante in quanto la volatilità è stocastica e risulta correlata col livello del tasso di interesse (effetto GARCH), inoltre esso non ha soluzione in forma chiusa. I due processi di Wiener e non sono indipendenti ma presentano una correlazione istantanea pari a

1Z 2Zρ . Al fine di simulare il processo è

necessario prima di tutto trasformare la rappresentazione nel tempo continuo fornita dalle equazioni differenziali nella sua controparte nel tempo discreto. Applicando lo schema di Eulero Kloeden et. al. (1991) otteniamo:

ttrkrr tttt ∆+∆−=− −−− 1111 )( ενµ tt tttt ∆+∆−=− −−− 2111 )( ενξναγνν

Per generare 1ε ed 2ε occorre prima di tutto estrarre due campioni, ϕ1 e ϕ2 , da due distribuzioni normali univariate dopodiché possiamo applicare la decomposizione di Cholesky che, nel caso semplice di due variabili, equivale a porre (in questo caso 1σ e

2σ sono entrambe pari a 1):

11 ϕε = 2

212 1 ϕρρϕε −+=

Dai valori iniziali della volatilità e il tasso istantaneo di interesse si possono ricavare i valori di queste grandezze al valore t0 + ∆t. I valori così ottenuti vengono utilizzati in un processo iterativo che ci permette di ricavare i valori desiderati al tempo generico t0 + n∆ t.

4.5 Come adattare una distribuzione non normale univariata ai dati empirici

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. In questo paragrafo intendiamo riprendere tale questione ponendoci in un’ottica leggermente diversa. Infatti oltre al problema della scelta di una forma funzionale per la distribuzione e della successiva determinazione dei suoi parametri, adesso abbiamo anche il problema dell’utilizzo della distribuzione così trovata ai fini della simulazione Monte Carlo. Diventa pertanto fondamentale seguire un approccio che tenga in debito conto da un lato la semplicità implementativa degli algoritmi di calcolo e dall’altro l’efficienza della procedura in termini di tempo macchina. Una proposta molto interessante che si muove proprio lungo questa linea è quella sviluppata da Ramberg et. al. (1979). Com’è noto una distribuzione di probabilità continua può


essere descritta facendo riferimento alla sua funzione di densità o, alternativamente, alla funzione dei percentili (se tale funzione esiste). 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. Il punto di partenza del lavoro di Ramberg et. al. è la funzione lambda di Tukey definita da

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

λλ −−= (4.14)

con la quale consente di adattare ai dati una vasta categoria di distribuzioni simmetriche . Essi hanno generalizzato la distribuzione (4.14) in una distribuzione con quattro parametri definita dalla funzione percentile seguente (distribuzione lambda generalizzata)

1p0 ≤≤

λ

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

λλ

1

43 −−+= (4.15)

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. La funzione di densità corrispondente alla (4.15) è

( )[ ] 14

13

243 )1( −− −+

= λλ λλλ

pppRf (4.16)

Nel lavoro già citato gli autori hanno tabulato i valori dei quattro parametri λ in funzione dei parametri di asimmetria e curtosi per varie distribuzioni con media nulla e standard deviation unitaria. I valori di λ3 e λ4 dipendono esclusivamente e congiuntamente dai valori di curtosi e asimmetria; dai valori tabulati si possono calcolare i valori corretti di λ1 e λ2 validi nel caso in cui la media sia non nulla e la varianza diversa da 1:

µσλσµλ += )1,0(),( 1*1 (4.17)

σλσµλ )1,0(),( 2*

2 = (4.18)


I parametri λ*1, λ*2, λ3 e λ4 così ottenuti possono essere inseriti nella (4.17) o nella (4.18) per ottenere la distribuzione che meglio approssima i dati di partenza. 4.6 Appendice A - Il condizionamento della matrice di covarianza Può accadere, specialmente quando si lavora con problemi ad elevata dimensionalità, che la matrice di correlazione non risulti semi-definita positiva cioè che non sia soddisfatta la relazione

0≥tVCV 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

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

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

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

DC

DC 11

=′

dove D/1 è una matrice diagonale i cui elementi sono le radici quadrate degli elementi diagonali di C.

5 La Valutazione dei Titoli Derivati

He deals the cards to find the answer The sacred geometry of chance The hidden law of a probable outcome The numbers lead a dance Sting “The Shape of My Hart” Dall’album “Ten Summoner’s Tales”

5.1 La stima del prezzo La prima applicazione finanziaria del metodo Monte Carlo risale al 1977 quando Boyle presentò un metodo per la valutazione di opzioni europee su titoli che pagano dividendi. Già in questo primo lavoro vengono messi in luce gli aspetti peculiari delle applicazioni del metodo Monte Carlo in finanza, in particolare si sottolinea come la velocità del metodo stesso possa risultare talora insoddisfacente e come la conoscenza della soluzione esatta di un problema simile (in questo caso quello dell’opzione europea scritta su un titolo che non paga dividendi) possa essere usata con successo per ridurre la varianza del prezzo stimato numericamente. Nel 1990 Kemma e Vorst (Kemma e Vorst, 1990) affrontano con successo il problema di valutare un’opzione asiatica il cui valore dipende da una media aritmetica di prezzi. Il problema è interessante perché in questo caso è nota la soluzione analitica corretta per un problema simile cioè quello di un’opzione asiatica il cui valore dipende dalla media geometrica dei prezzi. Come vedremo in questo caso possiamo usare la soluzione esatta come indicatore per ridurre la varianza della stima ottenibile con il metodo Monte Carlo5. Nell’ipotesi di assenza di arbitraggio il prezzo di un derivato può essere espresso come il valore atteso dei flussi scontato al tasso risk-free. Come abbiamo già accennato nell’ introduzione, il metodo Monte Carlo può essere impiegato per stimare direttamente questo valore di aspettazione; da un punto di vista del tutto generale il calcolo può essere sviluppato come segue

1. Si genera un numero sufficientemente elevato di percorsi stocastici (path) che descrivono la dinamica delle variabili di stato del sottostante (ad esempio un prezzo o un tasso di interesse). La simulazione viene effettuata sull’orizzonte temporale di interesse che pertanto deve essere opportunamente discretizzato. I diversi path vengono generati in accordo con la misura risk-neutral6;

2. Per ogni realizzazione occorre valutare il valore attuale dei flussi scontato al tasso risk-free;

3. Si calcola la media dei valori scontati su tutti i percorsi.

5 Il problema della valutazione Monte Carlo del prezzo di un’opzione asiatica a media aritmetica e l’utilizzo di metodi di riduzione della varianza basati sul pricing dell’asiatica a media geometrica verranno discussi in dettaglio nella sezione dedicata al metodo della variabile di controllo. 6 Questo punto meriterebbe di essere discusso in maggior dettaglio, qui ci limitiamo ad osservare che tale assunzione equivale a ritenere che sia possibile (almeno in linea di principio) costruire una posizione coperta per l’opzione in esame. Quando questo non è possibile, la generazione dei path stocastici dovrebbe essere condotta secondo la misura di probabilità oggettiva e non secondo quella aggiustata per il rischio. Un esempio di questo caso è dato dal processo di pricing di una real option.


Il processo di generazione del prezzo sottostante avviene seguendo il modello di Black e Scholes a partire dalla relazione

( ) wtrSSST

T ∆+∆

−=

=∆

−

σσ 2

1 21lnln (5.1)

Per un processo di Wiener generalizzato si ha Ztw ∆=∆ essendo Z una variabile aleatoria distribuita secondo una normale standard. A partire dalla (5.1) si ottiene facilmente

Zttr

TT eSS∆+∆

−

−=σσ 2

21

1 (5.2a) Possiamo quindi generare iterativamente una serie di valori a partire da quello iniziale applicando la (5.2) ad ogni passo del calcolo. Il codice VB riportato sotto può essere utilizzato a tale scopo, la variabile t∆ è pari all’intervallo T complessivo diviso per il numero di passi.

Function BrownianPath(Prezzi() As Single, _ r As Single, _ dt As Single, _ sigma As Single) Dim i As Integer Dim NrStep As Integer Dim z As Double ' il Nr di passi viene posto uguale alla dimensione del vettore dei prezzi NrStep = UBound(Prezzi) For i = 2 To NrStep z = Application.NormSInv(Rnd(Timer)) Prezzi(i) = Prezzi(i - 1) * _ Exp((r - 0.5 * sigma * sigma) * dt + sigma * Sqr(dt) * z) Next End Function

Nella figura 5.1 e riportato un esempio di generazione di path ottenuto utilizzando la funzione sopra descritta. Nel caso di pricing di opzioni europee la (5.2) può essere impiegata per produrre direttamente il prezzo a scadenza

ZTTrT eSS σσ +−= )2/(

02

(5.2b)


Il metodo Monte Carlo può essere esteso con relativa facilità ai casi in cui vengano effettuati pagamenti durante la vita del derivato ma non può essere altrettanto facilmente esteso per valutare i derivati di tipo americano. Nell’ultimo paragrafo passeremo brevemente in rassegna alcuni approcci proposti per estendere l’applicazione del metodo anche a questo tipo di prodotti finanziari.

75

85

95

105

115

125

135

0 0,05 0,1 0,15 0,2 0,25

Figura 5.1 5.2 La stima delle sensitività In pratica la valutazione dei fattori di sensitività (greek letters) si dimostra di importanza almeno pari (e in taluni casi maggiore) della determinazione dei prezzi. Infatti mentre questi ultimi possono essere osservati sul mercato le sensitività ai vari fattori di rischio tipicamente non sono immediatamente disponibili ma devono essere calcolate. Nel resto di questa sezione discuteremo vari approcci alla stima delle sensitività e in particolare al delta. 5.2.1 Approssimazione basata sulle differenze finite Consideriamo il problema legato al calcolo del Delta di un’opzione europea

0SC

∂∂

=∆

essendo C il prezzo dell’opzione e S0 il livello corrente del prezzo del sottostante. Poniamoci per semplicità nell’ipotesi in cui valgano le condizioni del modello di Black e Scholes. E’ chiaro che, in questo caso, non ha senso stimare il Delta tramite una


simulazione poiché è disponibile una soluzione in forma chiusa al nostro problema ma l’esempio è utile a fini didattici. Un approccio diretto al problema può consistere nella generazione di due prezzi finali, il primo

ZTTrT eSS σσ +−= )2/(

02

(5.4) a partire dal valore S0 , e il secondo a partire dal valore perturbato S0 + ε

')2/(0

2)( ZTTr

T eSS σσε +−+= (5.5) dove Z e Z’ sono due estrazioni indipendenti da una normale standard. Per ogni valore del prezzo finale possiamo poi calcolare il valore dell’opzione corrispondente

( ) ),0max(0 KSeSC TrT −= −

( ) ))(,0max(0 KSeSC TrT −=+ − εε

una stima “brutale” del Delta può essere ottenuta approssimando la derivata con una differenza finita

( ) ( )ε

ε 00~ SCSC −+=∆ (5.6)

Generando n valori indipendenti per ST e ST(ε) possiamo calcolare la media campionaria di n valori indipendenti di ∆

~ . Nel limite in cui n tende all’infinito la media del campione tende al valore della derivata di C rispetto ad S. La discussione appena intrapresa sembra condurre alla conclusione che per avere una stima sufficientemente corretta del valore di Delta sia necessario considerare valori di ε sempre più piccoli. Tuttavia poiché i valori per ST e ST(ε) sono generati indipendentemente l’uno dall’altro abbiamo

( ) [ ] [ ] )()()( ~ 200

2 −− =++=∆ εεε OSCVarSCVarVar (5.7) Pertanto la varianza di ∆

~ 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. Un risultato del tutto generale dedotto da Glynn (1989) dimostra che il miglior rate di convergenza possibile con questo metodo è tipicamente dell’ordine di . Sostituendo la stima precedentemente discussa basata sulla (5.4) con quella basata sulla cosiddetta differenza centrata

4/1−n

( ) ( )

εεε

2~ 00 −−+

=∆SCSC (5.8)


migliora la velocità di convergenza ma il risultato rimane sempre dell’ordine di n in ogni caso inferiore al tipico andamento come caratteristico del metodo Monte Carlo.

3/1−

2/1−n

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 nella (5.4) che nella (5.5). 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(S

ˆ

∆0) e C(S0 + ε) non sono più indipendenti

( ) [ ] [ ] [ ] )(),(2)()(ˆ

00002 εεε +−++=∆ − SCSCCovSCVarSCVarVar (5.9)

Se C(S0) e C(S0 + ε) sono positivamente correlati allora ∆ ha una varianza minore di ∆~ 7. L’impatto della diminuzione della varianza è particolarmente significativo quando ε è piccolo. Un semplice calcolo mostra che, utilizzando lo stesso numero casuale Z sia nella (5.4) che nella (5.5), possiamo scrivere

ZTTrTT eSSSCSC σσεεε +−≤−≤−+ )2/(

00

2

)()()( (5.10) Poiché il limite superiore ha momento secondo finito possiamo concludere che

( ) )()()( 2200 εε OSCSCE =−+ (5.11)

pertanto

( ) ( ) )1(00 OSCSCVar =

−+

εε (5.12)

cioè la varianza di rimane limitata quando ε tende a zero. Pertanto tanto più precisamente cerchiamo di stimare il valore di Delta, prendendo ε piccoli, e tanto maggiore sarà il beneficio in termini di convergenza dovuto all’utilizzo del metodo appena discusso.

∆

Purtroppo i confortanti risultati appena visti non valgono in generale. In particolare non valgono quando i payoff delle opzioni coinvolte sono discontinui o presentano altre caratteristiche di irregolarità. Prendiamo ad esempio il caso di un’opzione digitale che paga un ammontare prefissato B se ST > K e 0 nel caso opposto. Il prezzo dell’opzione è dato da e e quindi uno stimatore per il metodo Monte Carlo è fornito da )( KSPB T >rT−

7 Nel caso in esame questo segue dalla monotonicità della funzione che lega Z a C.


( ) BeSC rTKST

−>= 1)( 0 (5.13)

poiché C(S0) e C(S0 + ε) differiscono solo quando )(εTT SKS ≤≤ in questo caso abbiamo

( ))())/1((

))(()()(

022

22200

εε

εε

OSSKSPeB

SKSPeBSCSCE

TTrT

TTrT

=+≤≤=

≤≤=−+−

−

(5.14)

di conseguenza la stima del Delta per un’opzione digitale risulta più difficile in quanto la stima della varianza non è più indipendente dal valore di ε. 5.2.2 Stime dirette Anche se l’utilizzo del metodo dei common random numbers può portare ad un miglioramento delle performance dei programmi di simulazione per la stima delle sensitività in ogni caso la stima basata sulle differenze finite soffre in ogni caso di una serie di controindicazioni. La più importante è senz’altro quella legata al fatto che per ogni parametro rispetto al quale vogliamo stimare la derivata dobbiamo eseguire due simulazioni, la prima con il valore del parametro al quale vogliamo effettuare il calcolo, la seconda assegnando a questo valore una perturbazione ε. Se questo non sembra di grande rilievo per la stima del Delta di un’opzione di tipo standard si pensi al caso in cui si vogliano stimare tutte le greche di un’attività finanziaria (come qualche strutturato) che può dipendere da oltre 10 fattori di rischio; è evidente che in questo caso riuscire a ridurre il numero di simulazioni richieste per eseguire una stima risulterebbe di grande aiuto. Per questo motivo, negli ultimi anni, sono stati sviluppati una serie di metodi che cercano di stimare una derivata tramite simulazione diretta ovvero senza ricorrere ad una duplice simulazione come nel caso della stima basata sulle differenze finite. Una delle tecniche più diffuse si basa sulla cosiddetta infinitesimal perturbation analysis (Glasserman, 1991). La stima del Delta avviene tramite il calcolo della derivata del prezzo campionato rispetto al valore del sottostante

εε

ε

)(ˆ)(ˆlim

ˆ00

00

SCSCdS

Cd −+=

→ (5.15)

se il limite esiste e se C(S0) e C(S0 + ε) sono calcolati con lo stesso Z, allora purché si abbia ST ≠ 0 possiamo scrivere

000

1ˆˆ

SSe

dSdS

dSCd

dSCd T

KSrTT

TT >

−== (5.16)


avendo usato la (5.2) per il calcolo della derivata di ST rispetto a S0 e tenendo presente che

<>

=−=−

−

KSkSe

KSdSde

dSCd

T

TrT

TT

rT

T , 0 ,

),0max(ˆ

(5.17)

Per ST > K, C non è derivabile ma poiché questo accade con probabilità zero la derivata è quasi sicuramente ben definita ovunque. La derivata così calcolata viene chiamata in letteratura pathwise derivative e può essere vista come caso limite dell’approccio basato sulla stima della differenza finita con common random numbers. Occorre naturalmente verificare se questo stimatore è corretto cioè se

[ ]CEdSd

dSdC

dSCdE ˆˆ

000

≡=

(5.18)

Pertanto la correttezza della derivata pathwise si riduce al fatto di valutare se è corretto scambiare fra loro gli operatori di derivazione e di valore di aspettazione. Nel mondo di Black e Scholes è estremamente facile dimostrare che tale relazione sussiste; per una dimostrazione più generale rinviamo a Broadie e Glasserman (1996). Sulla falsariga del calcolo del Delta possiamo stimare le altre greek letters8

)(1 KSrT

TKTe

drdC

≥−==ρ (5.19)

−−==Λ ≥

− TrSSSe

ddC TT

KSrT

T

2

0)( 2

1ln1 σσσ

(5.20)

−−−−==Θ ≥

−− TrSS

TSeKSre

dtdC TT

KSrT

TrT

T

2

0)( 2

1ln2

1)0,max( σ (5.21)

Ognuno di questi stimatori è corretto. Chiaramente nel contesto del modello di Black e Scholes non ha alcun senso utilizzare questi stimatori Monte Carlo in quanto esiste una soluzione in forma chiusa. Tuttavia, proprio per quest’ultimo motivo, è spesso utile avere un termine di riferimento. L’applicazione maggiore di questo tipo di tecniche si ha nel settore dei derivati esotici, nell’articolo già citato di Brodie e Glasserman vengono studiate in dettaglio (sia dal punto di vista analitico che da quello numerico) due applicazioni: opzioni asiatiche e opzioni con volatilità stocastica. Nel primo caso si dimostra che una stima corretta del Delta è data da

8 Sempre nel contesto del modello di Black e Scholes.


)~(0

1~

KSrT

SSe >

− (5.22)

essendo S~ la media dei prezzi del sottostante utilizzata per determinare il payoff dell’opzione. La valutazione di questa espressione richiede un tempo di calcolo trascurabile in confronto con quello necessario per una stima a partire da una simulazione con un valore del prezzo iniziale perturbato. L’applicabilità di questo metodo si basa essenzialmente sulla possibilità di scambiare l’operatore di derivazione col valore di aspettazione (vedi (5.18)), si può dimostrare che affinché questo requisito sia soddisfatto è necessario che il payoff del derivato sia una funzione continua nel parametro in questione. Un esempio in cui questo chiaramente non è vero è offerto dal caso dell’opzione digitale. A causa della discontinuità presente in ST = K il metodo appena visto non può essere applicato (almeno in questa forma semplificata). Un’altra situazione in cui ci imbattiamo in discontinuità è quella legata al calcolo del Gamma. Consideriamo infatti l’espressione per la stima del Delta nel caso di un’opzione europea (5.16) che qui riportiamo per convenienza

00

1ˆ

SSe

dSCd T

KSrT

T >−=

come si vede la presenza della funzione indicatrice rende d discontinua in S0/ˆ dSC T = K rendendo così inapplicabile una stima diretta del Gamma. Per risolvere questi problemi Broadie e Glasserman (1996) hanno costruito una serie di stimatori robusti che permettono il calcolo anche in presenza di un insieme finito di punti di discontinuità. Recentemente Fourniet et. al. (1999) e Benhamou (2000) hanno suggerito nuove metodologie di stima basate sul calcolo di Malliavin. 5.3 Un cenno alle problematiche dei derivati di tipo americano 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. Tuttavia ancora oggi non esistono soluzioni completamente soddisfacenti e il problema continua a rappresentare una sfida alle capacità degli studiosi del settore. 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). Il problema nell’applicazione del metodo Monte Carlo per la valorizzazione di derivati di questo tipo è rappresentato proprio dall’applicare una procedura intrinsecamente forward in time ad un problema che invece ne richiede una backward in time. A questo punto visto che esistono numerosi algoritmi di pricing basati su approssimazioni alle differenze finite ci si potrebbe chiedere qual è il vantaggio del metodo Monte Carlo e perché negli ultimi anni sforzi crescenti sono stati destinati alla soluzione del problema. Il punto fondamentale è da ricercarsi ancora una volta in una caratteristica peculiare del metodo Monte Carlo che lo rende estremamente attraente per chiunque si occupi di prodotti che dipendono da più sottostanti o, in generale, da più variabili aleatorie: la velocità di convergenza del metodo Monte Carlo è indipendente dalla dimensione del problema. Al contrario i metodi basati sulle differenze finite o sugli alberi binomiali richiedono risorse di calcolo e tempi crescenti (spesso in maniera proibitiva) all’aumentare delle dimensioni del problema. Indicando con h(Sτ) il payoff del derivato all’istante τ nello stato del mondo Sτ possiamo concludere che il problema è riconducibile a stimare la quantità

[ ])(max ττ

τSheEP r

T

−

≤= (5.23)

dove il massimo è calcolato su tutti gli stopping times τ inferiori o uguali a T9. Come abbiamo già accennato poco sopra il problema cruciale che distingue la simulazione di un prodotto di tipo americano da uno di tipo europeo è proprio la stima del tempo ottimale di esercizio. Se ammettiamo che i tempi di esercizio siano rappresentabili da un set discreto di valori ti con i = 1, 2, … , d, allora una soluzione ovvia potrebbe essere quella di generare un cammino casuale per ciascuno stato e da questi calcolare la quantità

itS

[ )(maxˆ

,,0 i

it

rt

diSheP −

==

K] (5.24)

tuttavia che questo stimatore non è privo di errore sistematico in quanto risulta

[ ] PPE ≥ˆ Questo risultato purtroppo sembra di natura generale, Broadie e Glasserman (1997) hanno analizzato una vasta classe di stimatori possibili e i loro risultati confermano che

9 La formula è facilmente generalizzabile al caso multi asset. Una certa attenzione invece deve essere posta alla discretizzazione temporale implicita in ogni simulazione. Questa discretizzazione riduce la possibilità di esercizio ad un set discreto di tempi tutti minori o uguali alla scadenza. Usualmente questo fatto non porta problemi di natura rilevante.


(almeno all’interno della classe da loro esaminata) non esistono stimatori corretti. Da questa constatazione gli autori hanno elaborato un metodo di simulazione che si basa su due stimatori entrambi non corretti il primo per eccesso il secondo per difetto. Utilizzando entrambi gli stimatori gli autori mettono a punto una procedura basata sulla programmazione dinamica che converge verso il corretto risultato. 5.4 Appendice A - Funzioni VB per la simulazione del Delta Le funzioni i cui listati sono riprodotti in questa appendice si trovano nel file “Simulazione delle Sensitività.xls”. La prima di queste funzioni, OPTDelta_MC_1 implementa il calcolo “brutale” con due simulazioni distinte per ogni stima della derivata, in essa è presente un flag che specifica se per la generazione di variabili uniformemente distribuite fra 0 ed 1 si desideri usare un generatore di numeri casuali standard oppure una sequenza deterministica basata sull’algoritmo di Halton discusso più avanti.

Public Sub OPTDelta_MC_1(Tipo As String, _ Flag As Integer, _ NSimulazioni As Long, _ PrezzoSottostante As Double, _ Strike As Double, _ r As Double, _ t As Double, _ sigma As Double, _ Opt_Delta As Double, _ Opt_StDev As Double) Dim s As Double ' Prezzo del sottostante Dim Prezzo_1 As Double ' Prezzo dell'opzione nella prima simulazione Dim Prezzo_2 As Double ' Prezzo dell'opzione nella seconda simulazione Dim Random_1 As Double Dim Random_2 As Double Dim epsilon As Double ' Perturbazione ReDim DeltaValues(NSimulazioni) As Double Opt_Delta = 0 Opt_StDev = 0 epsilon = 0.01 r = 0.01 * r sigma = 0.01 * sigma For i = 1 To NSimulazioni ' ' A seconda del valore del flag si utilizza un generatore uniforme ' oppure il generatore a bassa discrepanza basato sulla sequenza di ' Halton. In ogni caso si utilizza la funzione di Excel NormSInv ' per ottenere in output un numero distribuito secondo una normale ' standard. ' If Flag = 0 Then Random_1 = Application.NormSInv(Rnd(Timer)) Randomize (Timer) Random_2 = Application.NormSInv(Rnd(Timer)) Else Random_1 = Application.NormSInv(Halton(i, 3)) Random_2 = Application.NormSInv(Halton(i, 5)) End If ' ' Generazione del primo valore per il prezzo. Si noti che viene ' utilizzata la cosiddetta differenza centrata.


' s = (PrezzoSottostante + epsilon) * _ Exp((r - 0.5 * sigma * sigma) * t + sigma * Sqr(t) * Random_1) Prezzo_1 = Application.Max(0, s - Strike) ' ' Generazione del secondo valore per il prezzo. ' s = (PrezzoSottostante - epsilon) * _ Exp((r - 0.5 * sigma * sigma) * t + sigma * Sqr(t) * Random_1) Prezzo_2 = Application.Max(0, s - Strike) DeltaValues(i) = Exp(-0.01 * r * t) * _ ((Prezzo_1 - Prezzo_2) / (2 * epsilon)) Next ' Infine calcoliamo la media ... For i = 1 To NSimulazioni Opt_Delta = Opt_Delta + DeltaValues(i) Next Opt_Delta = Opt_Delta / NSimulazioni ' ... e la varianza For i = 1 To NSimulazioni Opt_StDev = Opt_StDev + (DeltaValues(i) - Opt_Delta) ^ 2 Next Opt_StDev = Sqr(Opt_StDev / (NSimulazioni - 1)) / Sqr(NSimulazioni) End Sub

La stima del Delta con il metodo dei numeri casuali comuni è effettuata con la funzione OPTDelta_MC_2 di cui qui riportiamo solo il loop di generazione (si noti l’utilizzo di un solo numero casuale per ogni simulazione del Delta).

For i = 1 To NSimulazioni If Flag = 0 Then random = Application.NormSInv(Rnd(Timer)) Else random = Application.NormSInv(Halton(i, 3)) End If s = (PrezzoSottostante + epsilon) * _ Exp((r - 0.5 * sigma * sigma) * t + sigma * Sqr(t) * random) Prezzo_1 = Application.Max(0, s - Strike) s = (PrezzoSottostante - epsilon) * _ Exp((r - 0.5 * sigma * sigma) * t + sigma * Sqr(t) * random) Prezzo_2 = Application.Max(0, s - Strike) DeltaValues(i) = Exp(-0.01 * r * t) * _ ((Prezzo_2 - Prezzo_1) / (2 * epsilon)) Next

La funzione OPTDelta_MC_3, infine, calcola il Delta applicando direttamente la (5.16) ad ogni simulazione (anche in questo caso riportiamo solo il codice essenziale rinviando il lettore al file Excel per lo studio del codice completo).

For i = 1 To NSimulazioni If Flag = 0 Then random = Application.NormSInv(Rnd(Timer)) Else random = Application.NormSInv(Halton(i, 3)) End If s = PrezzoSottostante * _ Exp((r - 0.5 * sigma * sigma) * t + sigma * Sqr(t) * random)


If s > Strike Then DeltaValues(i) = Exp(-0.01 * r * t) * s / PrezzoSottostante Else DeltaValues(i) = 0 End If Next

6 La stima del Value-at-Risk

Non c’è errore più comune di quello di presumere che, siccome sono stati effettuati calcoli matematici lunghi e accurati, l’applicazione del risultato a qualche realtà naturale sia assolutamente sicura. Alfred North Whitehead

6.1 Le problematiche 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. Con questo termine si indica una combinazione di due opzioni, una call e una put, con lo stesso prezzo di esercizio (strike price). Il payoff a scadenza del portafoglio composto da questi due titoli non è descrivibile con una funzione monotona. Questo strumento consente di scommettere sulla volatilità del mercato infatti la massima perdita per l’investitore si verifica se alla data di scadenza l’attività sottostante ha un valore pari al prezzo di esercizio, in questo caso minore è la volatilità del titolo sottostante e più alta sarà la perdita di valore del portafoglio. In questi casi una metodologia più precisa può essere rappresentata dalla scelta di opportuni metodi numerici tramite i quali simulare un certo numero (usualmente elevato) di scenari dai quali estrarre le stime di rischio del portafoglio. Come abbiamo visto nelle precedenti sezioni, le procedure basate sul metodo Monte Carlo permettono di trattare problemi che coinvolgono strumenti finanziari con payoff non lineari e possono essere applicate a derivati il cui valore dipende da una funzione dell’intero sentiero temporale seguito dalla variabile e non solo dal valore finale. Per contro tali tecniche possono essere molto onerose sotto il profilo computazionale (è indispensabile individuare caso per caso la tecnica ottimale per rendere il più efficiente possibile la simulazione) pertanto qualunque decisione operativa circa il loro impiego risulta quasi sempre l’esito di un compromesso fra precisione richiesta e tempo a disposizione. Sebbene 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. 6.2 Il rischio di mercato Da un punto di vista operativo nei modelli di valutazione basati su tecniche Monte Carlo il Value-at-Risk viene misurato direttamente dalla distribuzione di probabilità delle variazioni di valore del portafoglio, ottenute dalla simulazione, tagliando tale distribuzione al percentile desiderato. A scopo illustrativo consideriamo una posizione lunga in un’opzione il cui valore dipenda da due sottostanti: un indice azionario e un tasso di interesse. La procedura di calcolo del VaR può essere descritta come segue

1. Si genera un numero molto elevato di numeri casuali, ad esempio 10 000, per il valore dell’indice e del tasso di interesse alla data di scadenza dell’opzione. La generazione dei valori deve avvenire in maniera coerente con le proprietà statistiche delle quantità che intendiamo simulare in particolare volatilità e correlazione. Tali grandezze possono essere dedotte dall’analisi di serie storiche dei due sottostanti.

2. Dai valori precedentemente generati si valuta, per ciascuno scenario, il valore dell’opzione. Occorre notare che in questa fase è implicita la scelta di un modello finanziario di valutazione come discuteremo meglio nel paragrafo §.

3. Si genera l’istogramma corrispondente alla distribuzione dei risultati del valore dell’opzione.

4. Supponendo di voler calcolare il VaR al 95% di livello di confidenza si prende come risultato il valore della 500-esima peggiore osservazione (nell’ipotesi di aver generato 10 000 campioni).

Si noti che tale approccio rende irrilevante la presenza di titoli con payoff non monotono, infatti le variazioni di valore del nostro portafoglio vengono ordinate dalla peggiore alla migliore indipendentemente dal movimento del fattore di mercato che le ha generate. Il metodo Monte Carlo consente di superare il vincolo dell’ipotesi di normalità assunto dai modelli basati sull’approccio varianza-covarianza. Vale la pena di spendere alcune parole su questo punto oggetto, fra l’altro, di una certa confusione. Non è infrequente infatti leggere che il Metodo Montecarlo è intrinsecamente legato all’assunzione di


normalità e che, di conseguenza, sotto questo aspetto non è migliore degli altri metodi parametrici. In realtà finché ci limitiamo alla simulazione di portafogli sensibili ad un solo fattore di rischio o a più fattori privi di correlazione, tale limitazione è del tutto inesistente in quanto non solo è possibile generare campioni di numeri casuali estratti da una distribuzione descritta da una funzione qualsiasi ma esiste anche una scelta abbastanza ampia di tecniche per adattare una distribuzione non normale univariata ai dati storici. Diverso è il discorso qualora si debbano generare simulazioni per la variazione di valore di portafogli sensibili a più fattori di rischio fra loro correlati. In questo caso l’analisi numerica mette a disposizione tecniche numeriche relativamente semplici per la generazione di campioni di variabili aleatorie multivariate solo nel caso in cui la distribuzione sottostante sia normale. Il problema della generazione di variabili non normali multivariate presenta notevoli difficoltà analitiche legate sia alla generazione dei dati ma soprattutto al loro successivo adattamento ai dati misurati, per questo motivo il loro utilizzo pratico è estremamente limitato nella realtà operativa. 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) e 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. 6.2.1 La 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à. 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

( )∑=j

jj

mportfolio ttXPMPV ),()( (6.1)


dove Xj(t) sono i valori dei fattori di mercato da cui dipende il valore del titolo. La variazione di valore del portafoglio può essere quindi espressa come

L+∆∆∆∂∂∂

∂+

∆∆∂∂

∂+∆

∂∂

+

∆∆∂∂

∂+∆

∂∂

+∆∂∂

=∆

∑

∑∑

∑∑ ∑

≠≠

≠

≠

lkjlkj

lkj

kjkj

kjjj

j

kjkj kjj j

jj

jj

portfolio

XXXXXX

PV

XXXX

PVXXPV

XXXX

PVXXPVX

XPVPV

3

22

33

3

3

22

2

2

61

61

61

21

21

(6.2)

Il polinomio (6.2) è una forma di rivalutazione parametrica 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. 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.

7 Metodi di riduzione della varianza

Tutto ciò che esiste nell’universo è frutto del caso e della necessità. Democrito

7.1 Introduzione Come abbiamo già osservato esiste una corrispondenza fra il metodo Monte Carlo e il calcolo dell’integrale di Lebesgue-Stiltjes

∫=Z

LS zdFzI )()( rrκ (7.1)

La rappresentazione di questo integrale in termini di )(zk r e dF )(zr non è unica. Una generica copia di funzioni ZzzFzk ∈

rrr ),(),( che soddisfa la (7.1) rappresenta un piano di campionamento. Da una qualunque di queste coppie è possibile generare dei campioni che permettono di stimare il valore dell’integrale senza errori sistematici. In questo paragrafo mostreremo come sfruttare la possibilità di scegliere piani di campionamento alternativi al fine di generare serie di dati con varianza contenuta in maniera computazionalmente efficiente. Il problema dell’efficienza computazionale è molto importante, infatti, nella scelta di un piano di campionamento occorre tenere presenti non solo le esigenze di precisione desiderate ma anche il costo, in termini di tempo macchina, che tale scelta richiede. Vale pertanto la pena spendere due parole sull’argomento. 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

e . Supponiamo poi che entrambi gli stimatori siano corretti,

cioè valga

niPi ,...,1,1 = niP i ,...,1,2 =

[ ] P=1PE e [ ] P=PE 2 ma con 21 σσ < . 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 potrebbe richiedere molto più tempo che generare n replicazioni di .

1P

2P Per comparare fra loro stimatori diversi sotto il profilo dell’errore e delle risorse computazionali richieste è necessario quindi procedere ad un’analisi più attenta. 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 possa essere espresso da una costante che indicheremo con b

jPj10, avendo a disposizione un tempo totale di calcolo pari a t il numero

di replicazioni di che possiamo generare sarà pari a jP jbt / . I due stimatori possono pertanto essere riscritti introducendo esplicitamente il tempo di calcolo nelle formule

∑=

1/

1

11 ˆbt

iiP

tb e ∑

=

2/

1

22 ˆbt

iiP

tb

Per valori sufficientemente alti di t queste due quantità sono normalmente distribuite con media P e standard deviation

tb1

1σ e t

b22σ (7.2)

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

2221

21 bb σσ < (7.3)

L’equazione (7.3) ci permette di avere un’idea più precisa dei requisiti che devono guidare le nostre scelte in termini di performance di calcolo. Il 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 cioè se vale la (7.3). 7.2 Il metodo delle 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. La generazione di prezzi avviene come ormai è noto a partire dalla equazione (5.2b) tramite la quale possiamo generare un numero n di repliche del prezzo finale del sottostante

nieSS iZTTriT ,,1 )2/(

0)( 2

K== +− σσ (7.4)

10 In molti problemi anche questo tempo in realtà è una variabile aleatoria. Per semplicità noi considereremo solo il caso costante.


Se si adotta la tecnica della variabile antitetica, in ogni simulazione si devono determinare due valori. Il primo valore C è quello calcolato nel modo consueto

utilizzando le n repliche descritte dalla (7.4), il secondo valore Ci

i~ , viene calcolato

cambiando segno a tutti i campioni estratti casualmente dalle distribuzioni normali standardizzate. 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

( ) ([ iiiii CCCovCVarCCVar ~,

21

2)]

~+=

+ (7.5)

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 ] (7.6) Dalla (7.5) vediamo che questa condizione è equivalente a richiedere che

( ) 0~, ≤ii CCCov (7.7) Che questa condizione sia sempre verificata nel caso delle variabili antitetiche può essere dimostrato come segue. Indichiamo con ϕ la funzione definita dalla relazione

)( ii ZC ϕ= ; ϕ è la composizione di due funzioni monotone, la prima è quella che lega il valore del sottostante alla variabile aleatoria Z (Eq. (7.4)), 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. Pertanto anche ϕ è monotona. Utilizzando una disuguaglianza standard possiamo allora verificare che

[ ] [ ] [ ])()()()( iiii ZEZEZZE −≤− ϕϕϕϕ (7.8) da cui segue immediatamente

( ) [ ] [ ] [ ] 0)()()()(~, ≤−−−= iiiiii ZEZEZZECCCov ϕϕϕϕ Quindi il metodo delle variabili antitetiche è in ogni caso più efficiente del metodo standard.


7.3 Il metodo della variabile di controllo Come abbiamo accennato 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 Zzzk ∈

rr),( correlata a k(z) di cui è noto

esattamente l’integrale

∫=ℑZ

C zdFzv )()( (7.9)

Definiamo W(β, z) come

[ ]CzvzkzW ℑ−−= )()(),( ββ (7.10) E’ facile verificare che

[ ] [ ] ℑ== )(),( zkEzWE β (7.11)

( )[ ] 222 2, v(z)k(z),v(z)k(z) σββσσzWVar +−=β (7.12)

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

2)(

)(),(*zv

zvzk

σ

σβ = (7.13)

nel qual caso

( )[ ] 222 1, k(z)k(z),v(z)k(z) σ)ρ(σzWVar ≤−=∗β (7.14) come si vede l’effetto di riduzione della varianza è tanto più elevato quanto maggiore è la correlazione fra la funzione iniziale e quella ausiliaria. Vale la pena spendere alcune parole sulla determinazione del coefficiente ottimale β*. Infatti se utilizzassimo lo stesso campione di variabili Z per stimare contemporaneamente E[k(z)], E[v(z)] e β* si creerebbe una correlazione non desiderata tra il parametro β e la stima dei valori di aspettazione delle due funzioni che produrrebbe delle distorsioni nella stima del valore di aspettazione di W. Per questo solitamente si suddivide la simulazione in due parti, la prima destinata alla stima del fattore β* e la seconda destinata alla stima vera e propria del valore di aspettazione di W.

Il calcolo del VaR per portafogli contenenti opzioni offre un interessante esempio di applicazione di questa metodologia. Infatti se pensiamo di voler valutare P(L > x) possiamo utilizzare la funzione indicatrice I(Q > x) come variabile di controllo. Questo significa sostituire lo stimatore I(L > x) con


][ )()()( xQPxQIxLI >−>−> β (7.15)

dove L e Q sono valutate con lo stesso scenario di prezzi e P(Q > x) può essere calcolata numericamente. 7.3.1 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 e quello di una put scritta sul prezzo medio (average price put) è pari à 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. Se si assume che il prezzo dell’attività sottostante, S, sia distribuito in modo log-normale e che 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

( )XSave −,0

aveSX −,0max

aveS

( )aveS

m

m

jt j

SG/1

1

= ∏

=

(7.16)

è pari a

−

+−= )()(

21exp)exp( 21

2 dKNdNrTC GGG σµ (7.17)

dove

( ) mThhTqrSG / 22

1ln 20 =

+

−−+= σµ (7.18)

mmmhG 6

)1)(12(22 ++= σσ (7.19)

e

G

GG Kdσ

σµ 2

1)ln( +−

= Gdd σ−= 12

L’implementazione in Visual Basic dei calcoli riportati sopra si trova nel file “Riduzione Varianza 2.xls” nella funzione BS_Geometric il cui listato è riportato sotto


Function BS_geometric(sigma As Double, _ K As Double, _ m_ave As Integer, _ q As Double, _ r As Double, _ S_0 As Double, _ t As Double) As Double Dim h As Double ' intervallo fra due tempi di calcolo della media Dim mg As Double ' drift modificato Dim sg As Double ' standard deviation modificata Dim d1 As Double Dim d2 As Double h = t / m_ave mg = Log(S_0) + (r - q - sigma ^ 2 / 2) * (t + h) / 2 sg = Sqr(sigma ^ 2 * h * (2 * m_ave + 1) * (m_ave + 1) / (6 * m_ave)) d1 = (mg - Log(K) + sg ^ 2) / sg d2 = d1 - sg BS_geometric = Exp(-r * t) * _ (Exp(mg + sg ^ 2 / 2) * Application.NormSDist(d1) - _ K * Application.NormSDist(d2)) End Function

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. Nella figura 3 è riportato l’andamento del prezzo per tre opzioni (asiatica a media aritmetica, asiatica a media geometrica e europea standard) al variare della volatilità del sottostante e del tempo a scadenza (per i valori dei parametri impiegati si rinvia al file Excel). Nella simulazione Monte Carlo standard il prezzo dell’opzione viene calcolato come

[ ] niKAEeC irTi ,,1 ,)0,max( )()( K=−= − dove A(i) è la media aritmetica discreta campionata

∑=

=m

jt

iijS

mA

1

)()( 1

calcolata su un insieme discreto di punti

mjtmThhtt jj ,,2,1 ,0 , , 01 K===+= −

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


∑=

=n

i

iCn

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 AritmeticaEuropea Black & ScholesAsiatica 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 AritmeticaEuropea Black & ScholesAsiatica Geometrica

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


mm

j

it

ij

SG/1

1

)()(

= ∏

=

e il valore campionato dell’opzione asiatica a media geometrica

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

Questa volta però utilizzeremo lo stimatore (vedi Eq. (7.10))

∑=

+−=n

iG

iG

i CCCn

C1

)()( )(1ˆ

Nella figura 4 sono riportate le varianze campionarie di questo stimatore ottenute applicando il metodo Monte Carlo standard, il metodo delle variabili antitetiche e il metodo della variabile di controllo. Come si vede in ogni caso l’ultimo metodo è quello che produce risultati di gran lunga migliori. 7.4 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

smZZ i

i−

='

Correttamente la media dei campioni aggiustati è nulla e la deviazione standard è pari a 1. Questi campioni aggiustati vengono usati per effettuare le simulazioni. Il metodo dei momenti consente di ridurre i tempi di calcolo ma può creare problemi di memoria dato che ogni estrazione deve essere mantenuta in memoria fino alla fine della simulazione.


7.5 Importance Sampling Consideriamo due diversi piani di campionamento ZzzFzk ∈

rrr ),(),( e ZzzFzk ∈

rrr ),(),( ** soddisfacenti la (7.1). Supponiamo inoltre che )(zF r e )(z*F r siano entrambe differenziabili in Z e indichiamo con f(z) e f*(z) le rispettive densità di probabilità. Definiamo

Zz ∈== )()()(*

)()()(* zlzkzf

zfzkzk (7.20)

l(z) viene chiamata importance function. Dalla (7.20) segue immediatamente che

( )[ ] [ ] ℑ== )(** zkEzkE FFrr (7.21)

∫ ℑ−=Z

F dzzfzkzk 22 )()()(var r

r

(7.22)

∫ ℑ−=Z

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

da cui

∫ ℑ−−=−Z

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

Da queste formule risulta che, almeno in linea di principio, è possibile ridurre a zero la varianza campionata da )(* zF

r ; infatti è sufficiente scegliere

ℑ=

)()()(* zfzkzf (7.25)

per ottenere dalla (7.24) 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 (7.26)

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.


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,50Volatilità Sottostante

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,0Tempo a Scadenza

MC Standard MC ControlloMC Antithetic

Figura 7.2 7.5.1 Importance Sampling e stima del VaR per portafogli contenenti opzioni In questa sezione presentiamo un’interessante applicazione delle tecniche di Importance Sampling al caso particolare del calcolo del Value-at-Risk per portafogli contenenti titoli derivati con forte connotazione non lineare. L’approccio dovuto a P. Glasserman et. al. (1999a; 1999b; 1999c, 2000a, 2000b, 2000c) si basa sull’utilizzo dell’approssimazione Delta-Gamma non per la rivalutazione del portafoglio (per il cui scopo com’è noto tale approssimazione può essere estremamente insoddisfacente) bensì come guida per la determinazione di piani di campionamento efficienti sotto il profilo


)

della riduzione della varianza. Il punto essenziale è legato al fatto che l’utilizzo del metodo Monte Carlo per la stima del Value-at-Risk richiede un campionamento accurato delle code della distribuzione dei rendimenti. Tuttavia proprio perché nelle code cadono gli eventi meno probabili è chiaro che l’efficienza del metodo Monte Carlo è messa duramente alla prova. Lo scopo del lavoro che discutiamo in questo paragrafo è proprio quello di individuare, se esistono, dei cambiamenti di misura che rendano più probabili gli eventi delle code in modo da poter campionare con maggior efficienza proprio le regioni della distribuzione che risultato determinanti per una stima accurata del Value-at-Risk.

Riteniamo istruttivo seguire nei dettagli il calcolo che porta all’individuazione di una misura sotto la quale gli eventi che cadono nelle code non risultano più così rari come nel caso della distribuzione di partenza rinviando ai lavori già citati di Glasserman et. al. per ulteriori dettagli e una discussione approfondita di molte applicazioni pratiche. Consideriamo quindi il caso generale in cui si abbiano n fattori di rischio e sia

il loro valore al tempo t, con ∆S indicheremo la variazione nei fattori di rischio nell’intervallo [t, t +∆ t]. Qualora si consideri valida l’approssimazione Delta-Gamma, la variazione di valore di un portafoglio può essere sempre espressa nella forma (vedi Appendice)

))(,),(()( 1 tStStS nL=

(∑=

+=Λ+=n

iiiii

tt ZZbZZZbQ1

2λrrrr

(7.27)

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

Σ=tCC

CC tΓ21 , Cb δ

rr−= ,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 (7.28) dove a0 è una costante. Il VaR al livello di confidenza p è quel valore xp tale che

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

∫ >=> dzzfxLIxLP )()()( (7.30) 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 rendere più efficiente la simulazione (ad esempio riducendo la


varianza, o come nel caso attuale, rendendo più probabili gli eventi corrispondenti alle code della distribuzione originaria); con quest’obiettivo in mente poniamo dunque

∫ >=>=> )]()([~)()()(

)()( zlxLIEdzzgzgzf

xLIxLP (7.31)

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

21exp)( (7.32)

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. Questo approccio, diffuso nell’area di studio degli eventi estremi, è noto come exponential twisting. Poniamo quindi

( ) ( )[ ] ( ) ( )[

−−−= −− ϑµϑϑµϑ ZBZBzg t 12/1

21exp)( ] (7.33)

( ) ( )( ) ( ) 121 −Λ−=

=

ϑϑ

ϑϑϑµ

B

bB (7.34)

possiamo allora scrivere l’espressione della important function

( ) ( )[ ] ( ) ( )[ ]

−−−

−

=−− ϑµϑϑµϑ ZBZB

ZZzl

t

t

12/1

21exp

21exp

)( (7.35)

che, una volta sostituite le espressioni (7.34), si semplifica in

( ) ( )[ ] ( )[ ]ϑψϑϑψϑ +−=+Λ+−= QZZZbzl tt expexp)( (7.36) dove ( )ϑψ è data da

( ) ( ) (∑=

−−−

=n

ii

i

ib

1

2

21ln212

1 ϑλϑλ

ϑϑψ ) (7.37)

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 −−≤>+=>=> − ϑϑψϑϑψ (7.38) non è difficile ricavare anche un limite superiore per il momento secondo

[ ] ( )[ ]02

2 2)(2exp)()(~),( axzlxLIExm −−≤>= ϑϑψϑ (7.39) 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 −=′ ϑψ (7.40) Vediamo ora quanto vale il valore di aspettazione di Q sotto la nuova misura g(z)

[ ] ( ) ( ) ∫∫ ∫ −− === dzzfeQedzzfeQdzzgQQE QQ )()()(~ ϑϑψϑψϑ (7.41) Per definizione la funzione generatrice dei momenti di Q è data da

( ) [ ] ∫==Φ dzzfeeE QQ )(ϑϑϑ (7.42) da cui

( ) ∫=Φ dzzfeQdd Q )(ϑϑϑ

(7.43)

quindi, ricordando che ( ) ( )ϑϑψ Φ= ln ,

[ ] ( ) ( ) ( )ϑψϑϑ

ϑψ ′=Φ= −

ddeQE~ (7.44)

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. 7.6 Campionatura uniforme 7.6.1 Il 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 (vedi Fig. 7.3a).

senza campionamento stratificato

0

2

4

6

8

10

12

14

16

18

-3

-2,5 -2

-1,5 -1

-0,5 0

0,5 1

1,5 2

2,5 3

Classe

Freq

uenz

a

con campionamento stratificato

0

2

4

6

8

10

12

-3

-2,5 -2

-1,5 -1

-0,5 0

0,5 1

1,5 2

2,5 3

Classe

Freq

uenz

a

a) b) Figura 7.3

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 calcolare11

−+

= −

1001~ 1 i

iUiNU

Come si vede dalla figura 7.3b questo semplice accorgimento migliora molto i risultati della simulazione. Da un punto di vista più formale, se indichiamo con un insieme di

intervalli disgiunti di Z tali che , possiamo scrivere rZZZ ,...,, 21

ZZr

ii =

=U

1

11 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.


∫ ∑ ∫ ∑ ∫ ∑ ℑ====ℑZ

iZ

iZ

i iiii

i i

pzdFzkpzdFzkzdFzk )()()()()()( (7.45)

dove

ii p

zdFzdF )()( = (7.46)

∫=ℑiZ

ii zdFzk )()( (7.47)

e

∫=iZ

i zdFp )( (7.48)

Ci riferiremo all’intervallo Zi come allo strato i-esimo. Il metodo consiste nella generazione di ni campioni indipendenti Z per ogni strato con i che va da 1 a r. Si costruisce quindi lo stimatore

iij ,...., n j 1 ),( =

∑ ∑= =

=ℑr

i

n

j

ji

i

ir

i

zknpnn

1 11 )(),,( K (7.49)

Questo stimatore è non distorto ed ha varianza pari a

∑=

=ℑr

i i

iir n

pnn1

2

1 )],,(var[ σK (7.50)

dove

[ ]∫ ℑ−=ℑ−=iZ

ii

ii

i zdFzkp

zkE )()(1])([ 222σ (7.51)

Il grado di riduzione della varianza dipende dal comportamento della funzione integranda all’interno di ciascun intervallo; due problemi principali emergono in ogni applicazione del metodo che stiamo esaminando. Il primo è legato alla scelta della partizione , il secondo, non meno importante, riguarda il numero di campioni n

rZZZ ,...,, 21

i da generare per ciascuno strato. Entrambi i problemi sono discussi in termini generali in letteratura Sironi (1997) e Fishman (1996) per cui non verranno trattati in dettaglio in questa sede. 7.6.2 Il metodo dell’Ipercubo Latino In linea di principio non vi sono difficoltà ad estendere l’idea del campionamento stratificato al caso di più dimensioni. Per generare un campione stratificato nell’ipercubo a d dimensioni con n strati in ciascuna dimensione è sufficiente produrre una sequenza di vettori U e poi prendere le variabili definite da ),,( )()1( d

jjj UU K=


nUii

V jdj

1),( 1 −+=

K (7.52)

Il problema è che quando il numero di dimensioni aumenta diventa proibitivo generare anche un singolo strato. Per comprendere la natura del problema pensate al caso di una sola dimensione, se prendete il segmento [0, 1] e pensate di volerlo dividere in n parti uguali otterrete chiaramente n segmenti all’interno dei quali cadranno i vostri punti. Dov’è il problema si chiederà lettore. Bene, immaginiamo di prendere un quadrato di lato unitario e di suddividere ciascun lato in n segmenti, otteniamo evidentemente n2 quadrati, se ripetiamo la stessa cosa con un cubo otterremo n3 cubi e così via. In generale un ipercubo a d dimensioni contiene un numero di ipercubi di dimensione d – 1 che tende ad aumentare in maniera esponenziale all’aumentare del numero di dimensioni. Pertanto se il numero di dimensioni è molto elevato diventa molto probabile che una sequenza finita di punti sia tale che nessuno di essi cade all’interno di qualche ipercubo. Questo problema può essere in parte ovviato con un metodo che pur basandosi su un campionamento casuale del dominio delle variabili in esame cerca di mantenere una uniformità di campionamento. La tecnica nota come metodo dell’Ipercubo latino funziona nel modo seguente. Consideriamo d permutazioni indipendenti

dππ ,,1 K dell’insieme n,,1K distribuite uniformemente nell’insieme di tutte le n! permutazioni. Poniamo

nj dk n

jUV k

kjk

j ,,1,,,1,1)()(

)( KK ==−+

=π

(7.53)

L’aleatorietà ci assicura che ogni vettore Vj è uniformemente distribuito sull’ipercubo d dimensionale e, nello stesso tempo, le coordinate sono perfettamente stratificate nel senso che esattamente uno dei valori V cade fra (j-1)/n e j/n per ogni dimensione k.

)()(1 ,, k

nk VK

7.6.3 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à (il termine per ora va preso ancora nel suo senso intuitivo). 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 riatteremo 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 (il normale metodo Monte Carlo, ricordiamolo ancora una volta, ha un andamento legato all’inverso della radice quadrata del numero di punti e quindi risulta molto più lento in questo caso). 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. 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. Prima di andare avanti è indispensabile cercare di rendere più quantitativo il concetto di uniformità nella distribuzione di una serie di punti. Il grado di uniformità con cui una sequenza di punti riempie un intervallo assegnato è misurato dalla sua discrepanza: maggiore è l’uniformità minore la discrepanza. Questa grandezza è definita come segue. Supponiamo di avere una sequenza di n punti nxxx ,,, 21 K , nell’intervallo semiaperto d-dimensionale [ )ddI 1,0= , consideriamo poi un sottoinsieme J di Id, definiamo prima di tutto la quantità

)();();( JVn

nJAnJD −= (7.54)

dove A(J;n) è il numero di elementi della sequenza appartenenti a J e V(J) è il volume di J. La discrepanza della sequenza è definita come l’estremo superiore di D(J;n) per tutti i J appartenenti a Id

);(sup nJDDdIJ

n∈

= (7.55)

E’ utile definire anche la quantità (detta star discrepancy) come l’estremo superiore di D(J;n) calcolato su insiemi J della forma

*nD


[ )∏=

n

iiu

1,0 (7.56)

Poiché in realtà è questa la quantità che ricorre più sovente nello studio delle proprietà delle sequenze di punti “quasi-casuali” ci riferiremo d’ora in avanti a indicandola semplicemente come 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

*nD

10 21 ≤≤≤≤≤ nxxx K

nkx

nD k

nkn 2

12max21

,1* −

−+==

(7.57)

come si vede la discrepanza non è mai inferiore a 1/2n e il valore minimo è ottenuto per

nk n

kxk ≤≤−

= 12

12 (7.58)

nel caso di dimensioni maggiori non esiste in generale una formula semplice per il calcolo della discrepanza12. Tuttavia si può dimostrare che in generale (almeno per le sequenze a bassa discrepanza più utilizzate) l’andamento asintotico della discrepanza in uno spazio d-dimensionale è del tipo

≈

ndOD

n

nlog*

(7.59)

Per applicare la nozione di discrepanza all’analisi dei problemi legati all’integrazione di una funzione è necessario definire una misura della variabilità di una funzione assegnata all’interno di un certo dominio. Come al solito il caso più semplice da trattare è quello monodimensionale. DEFINIZIONE 7.1. La variazione totale di f, che indicheremo con V(f), sull’intervallo [0, 1] è definita da

∑=

−≤≤≤≤≤

−=n

jjj

xxxnxfxffV

n 11

10)()(supsup)(

10 K (7.60)

TEOREMA 7.1. Se f è a variazione totale finita, cioè V(f) < ∞, su I e la sequenza xj ∈ I, j = 1, … n ha discrepanza allora si ha *

nD

12 Una notevole eccezione è data da (FAURE)


)()()(1 *

1

1

0fVDdxxfxf

n nn

jj ≤−∑ ∫

= (7.61)

E’ possibile generalizzare questo risultato al caso di più dimensioni attraverso un’appropriata definizione della variazione totale di una funzione tuttavia non ci tratterremo oltre su questo risultato anche perché il suo valore è esclusivamente teorico risultando in pratica di pochissimo aiuto. Il limite teorico appena visto normalmente sovrastima l’errore che effettivamente si commette ed inoltre la funzione V(f) nel caso multidimensionale (che in pratica è quello che interessa) risulta difficilmente stimabile. Fra le sequenze più note ricordiamo la sequenza di Halton, quella di Faure e quella di Sobol. Delle prime due diamo di seguito una breve spiegazione mentre rinviamo a Press et. al. (1992) per una sintetica spiegazione del funzionamento della sequenza di Sobol. La sequenza di Halton Un semplice esempio è fornito dalla successione di Halton ottenibile con la seguente procedura

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

nnbdbdbdj +++= K1

10

0 • Si pone

1

02

11 −−−

−− +++= n

nnj bdbdbdH K 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, … )13. L’algoritmo per la generazione della sequenza di Halton è estremamente semplice e compatto, il listato VB è riportato nel riquadro seguente.

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

13 Si lascia come esercizio al lettore il compito di trovare la legge che determina la successione dei numeri primi; ci limitiamo a suggerire di dimostrare dapprima che ogni numero pari è esprimibile come somma di due numeri primi.


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

Figura 7.4 – 1000 punti di una sequenza di Halton bidimensionale La sequenza di Faure Descriviamo ora come generare una sequenza di Faure in s dimensioni. Sia r il più piccolo numero primo tale che r ≥ s e r ≥ 2, come nel caso precedente cominciamo con lo scrivere il generico numero n in base r

∑=

=m

j

jj rnan

0

1 )( (7.62)

il primo numero viene generato tramite riflessione rispetto al punto decimale come nel caso della sequenza di Halton


∑=

−−=m

j

jjr rnan

0

111 )()(φ (7.63)

i rimanenti elementi della sequenza vengono generati tramite una procedura ricorsiva. Noti i coefficienti allo step k – 1 , l’algoritmo calcola i coefficienti dello step attuale

r naCnam

ji

kij

ikj mod)()( 1∑

≥

−= (7.64)

a questo punto siamo pronti per calcolare i successivi punti della sequenza

sk rnanm

j

jkj

kr ,,2,)()(

0

1 K== ∑=

−−φ (7.65)

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

Figura 7.5 – 1000 punti di una sequenza di Faure bidimensionale

7.6.4 Metodi misti e decomposizione spettrale per problemi ad alta

dimensionalità La figura 7.6 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

Figura 7.6 – Coordinate 49 e 50 di una sequenza di Faure di 1000 punti in 50 dimensioni Per questo in molte applicazioni che richiedono la simulazione di un elevato numero di variabili l’approccio più diffuso consiste nel determinare un numero minimale di variabili critiche per le quali utilizzare una simulazione basata sulle sequenze a bassa discrepanza mentre il resto delle variabili viene campionato in maniera standard utilizzando sequenze aleatorie non deterministiche. Un caso particolarmente interessante di tecnica mista è quello messo a punto da Dupire e Savine (Dupire, 1998) per la generazione di path stocastici. Gli autori applicano la principal component analysis per determinare prima di tutto quali sono le dimensioni del problema che contribuiscono maggiormente alla varianza complessiva. Una volta individuati gli autovettori che spiegano la maggior parte della varianza questi vengono campionati utilizzando una delle sequenze deterministiche descritte sopra mentre i restanti autovettori vengono campionati con variabili casuali. 7.7 Appendice A - Deduzione della formula (7.9) Per semplicità assumiamo di avere un derivato sensibile a un solo fattore di rischio che identificheremo col prezzo del sottostante S. Utilizzando uno sviluppo in serie di Taylor arrestato al secondo ordine possiamo esprimere la variazione di prezzo dell’opzione in funzione della variazione di prezzo del sottostante come


( ) ( )22

2

2

21

21 SSS

SVS

SVV ∆Γ+∆=∆

∂∂

+∆∂∂

≈∆ δ (A.1)

I coefficienti di questa espansione δ e Γ rappresentano la sensitività al primo e secondo ordine del valore dell’opzione rispetto al valore del sottostante. Ai fini del nostro ragionamento è conveniente formulare le nostre espressioni in funzione dei rendimenti piuttosto che delle variazioni di prezzo. Dalla (A.1) si ottiene subito

22

21

∆

Γ+∆

=∆

SS

VS

SS

VS

VV δ (A.2)

da cui, ponendo

VS

VVr

SSr =

∆=

∆= η e ˆ ,

(A.3) si ha

22 ~21~

21ˆ rrrSrr Γ+=Γ+= δηδη (A.4)

essendo δηδ =

~e ηSΓ=Γ

~.

E’ fondamentale ricordare che questa approssimazione ha valore locale e non descrive il reale andamento del rendimento dell’opzione in funzione del rendimento del sottostante. Sebbene essa permetta, come vedremo tra breve, di ottenere una forma chiusa per la funzione di distribuzione di r , non bisogna mai dimenticare che un uso acritico di tale approssimazione può condurre a risultati errati. A questo punto, il passaggio successivo nella nostra analisi, consiste nel trovare un modo per calcolare la distribuzione di r o, quanto meno, per trovarne una buona approssimazione. Assumiamo che il rendimento del sottostante, r, sia distribuito secondo una normale con media 0 e standard deviation pari a σ

),0(~ σNr

e introduciamo una variabile normalizzata

σrr =∗

naturalmente avremo che r*~N(0,1). Possiamo riscrivere la (A.4) nella forma


222 ~

21~~

21~ˆ ∗∗∗∗∗∗ Γ+=Γ+= rrrrr δσσδ (A.5)

avendo posto

2~ ~ e ~~ σσδδ Γ=Γ= ∗∗ da cui

( )∗

∗

∗

∗∗∗

Γ−

Γ+Γ= ~

~

21

~~

~21ˆ

22

δδrr (A.6)

Nel caso in cui siano presenti più titoli possiamo facilmente generalizzare il nostro risultato riscrivendo la (A.5) in forma vettoriale

r rrr Γ′+⋅∆′= ~21~ˆ (A.21)

dove r è un vettore n-dimensionale essendo n il numero di fattori di rischio per il portafoglio in esame. Come al solito ipotizziamo che i rendimenti ri costituiscano un insieme di n variabili distribuite secondo una normale multivariata con matrice di correlazione Σ. Per la nostra analisi è conveniente ricondursi ad un set di variabili distribuite normalmente ma non correlate. A tal fine definiamo la seguente trasformazione

LrrAOr =′= −∗ 1 (A.22) dove A è una matrice tale che AA ′=Σ (A è la matrice risultante dalla decomposizione di Cholesky di Σ) mentre O è la matrice formata dagli autovettori della matrice simmetrica

AAA Γ′=Γ~~ (A.23)

La matrice di varianza e covarianza del nuovo set di variabili può essere facilmente calcolata dalla (A.22)

IOAAO =′Σ′=Σ −− 11* (A.24)


essendo I la matrice identità. La trasformazione (A.22) produce quindi un set di n variabili non correlate e, poiché sono distribuite secondo una normale multivariata, anche indipendenti14. Definiamo le nuove matrici di sensitività

∆′′=∆~~* AO e AOAO Γ′′=Γ

~~* Si noti che la nuova matrice *~

Γ risulta diagonale per costruzione. Con queste posizioni la (A.21) può essere riscritta come

***** ~21~ˆ rrr r Γ′+⋅∆′= (A.25)

14 Com’è noto n variabili non correlate non sono in generale indipendenti . Tuttavia nel caso nel caso di una distribuzione multivariata normale la funzione densità di probabilità è proporzionale a

( )[ ] [ ] ( ) ( )[ ] ∏∑==

−=

==Σ

n

ii

n

ii rrIrrrr

1

2*

1

2****1** expexp 'exp 'exp

come conseguenza le variabili ri risultano anche indipendenti.

Bibliografia* Benhamou E. (2000): “Determining Malliavin Weights in the Simulation of the

Greeks”, London School of Economics, Working Paper, March. Boyle P., Broadie M. e Glasserman P. (1997): “Monte Carlo Methods for Security

Pricing”, Journal of Economic Dynamics and Control, 21, 1267-1321. Broadie M., Glasserman P. (1996): “Estimating Security Price Derivatives by

Simulation”, Management Science, 42, 269-285. Broadie M., Glasserman P. (1997): “Pricing American Style Securities Using

Simulation”, Journal of Economic Dynamics and Control, 21, 1323-1352. Cherubini U., Della Lunga G. (2001): “Il Rischio Finanziario”, McGraw-Hill, Milano. Dupire B. (1998): (a cura di) Monte Carlo Methodologies and Applications for Pricing

and Risk Management, Risk Books, London. Fishman G. S. (1996): Monte Carlo Concepts, Algorithms and Applications, Springer-

Verlag, New York. Fournie E., Laszry J. M., Touzi N. (1999) : “Application of Malliavin Calculus to

Monte Carlo methods in finance”, Finance and Stochastics, 3, 391-412. Geweke J. (1995a): “Posterior Simulators in Econometrics”, Federal Reserve Bank of

Minneapolis, Research Department, Working Paper 555. Geweke J. (1995b): “Monte Carlo Simulation and Numerical Integration”, Federal

Reserve Bank of Minneapolis, Research Department Staff Report, 192. Glasserman P. (1991): “Gradient Estimation via Perturbation Analysis”, Kluwer

Academic Publishers, Norwell, MA. Glasserman P., Heidelberger P. e Shahabuddin P. (1999a): “Importance Sampling and

Stratification for Value-at-Risk”, IBM Research Report 21450, April 15. Glasserman P., Heidelberger P. e Shahabuddin P. (1999b): “Stratification Issues in

Estimating Value-at-Risk”, IBM Research Report 21548, August 18. Glasserman P., Heidelberger P. e Shahabuddin P. (1999c): “Variance Reduction

Techniques for Estimating Value-at-Risk”, IBM Research Report 21577, October 22.

Glasserman P., Heidelberger P. e Shahabuddin P. (2000a): “Efficient Monte Carlo Methods for Value-at-Risk”, IBM Research Report 21723, April 12.

Glasserman P., Heidelberger P. e Shahabuddin P. (2000b): “Variance Reduction Techniques for Value-at-Risk with Heavy-Tailed Risk Factors”, IBM Research Report 21816, August 3.

Glasserman P., Heidelberger P. e Shahabuddin P. (2000c): “Portfolio Value-at-Risk with Heavy-Tailed Risk Factors”, IBM Research Report 21817, August 3.

Hammersley J. M., Handscomb D. C. (1964): “Monte Carlo Methods”, Chapman & Hall, London.

Hull J. C. (1993): “Option, Futures and Other Derivatives”, 2nd Edition, Prentice Hall. Hull J. C. (2000): “Option, Futures and Other Derivatives”, 4th Edition, Prentice Hall. James F. (1968): Monte Carlo Phase Space, Lectures given in the Academic Training

Programme of CERN, Geneva.

* I lavori di Glasserman et. al. sono scaricabili dal sito: (www.research.ibm.com/people/b/berger/papers.html)


James F. (1988): “A Review of Pseudorandom Number Generators”, CERN Report

Data Handling Division DD/88/22, Geneva, 1988. Kloeden P. E., Platen E. e Schurs H. (1991): Numerical Solution of SDE through

Computer Experiment, Springer-Verlag, New York. Knuth D. (1981): The Art of Computer Programming, Vol. 2 Seminumerical

Algorithms, Addison-Wesley Publishing Company, Reading. Nelsen R. B. (1999): “An Introduction to Copulas”, Springer, New York. Press W. H., Flannery B. P., Teukolsky S. A. e Vetterling W. T. (1990) Numerical

Recipes in C, Cambridge University Press, Cambridge. Sironi A., Marsella M. (1997): (a cura di) La misurazione e la Gestione dei Rischi di

Mercato, Società Editrice Il Mulino, Bologna.

Introduzione al metodo monte carlo

Education

Transcript of Introduzione al metodo monte carlo