Quantizzazione: applicazioni al pricing di derivati · Questo lavoro ha l'obiettivo di illustrare...

POLITECNICO DI MILANO

FACOLTÀ DI INGEGNERIA

Corso di Laurea in Ingegneria Matematica

TESI DI LAUREA MAGISTRALE

Quantizzazione:

applicazioni al pricing di derivati

Relatore:Prof. Daniele Marazzina

Tesi di:Alessandra LietoMatricola: 799417

Anno accademico 2016/17

iii

Abstract

Questo lavoro ha l'obiettivo di illustrare le applicazioni della quantizzazione al pricing dei derivati nan-ziari.A seguito di un'introduzione teorica alla quantizzazione, presenteremo alcune nozioni generali sul pricingdi derivati, nonché sui modelli che tratteremo: Black & Scholes, Kou, Variance Gamma, Heston.Procederemo poi ad illustrare due metodi: Quantizzazione Marginale Ricorsiva (RMQ) e Quantizzazionecon trasformata di Fourier.Andremo quindi a stimare il prezzo di opzioni call vanilla per diverse dinamiche del sottostante, partendodai processi di Lévy (con e senza salti) no ad arrivare a modelli a volatilità stocastica.Per tutti i casi sopracitati procederemo a un confronto, in termini di accuratezza, con il metodo basatosulla FFT. Concluderemo con un'analisi dei risultati ottenuti ed un confronto dei metodi numerici uti-lizzati.

Abstract

The aim of this paper is illustrating the applications of quantization in the eld of Option Pricing.After a theoretical introduction about quantization, we will introduce some general notions about optionpricing and the models we will use: Black & Scholes, Kou, Variance Gamma, Heston.We will introduce two dierent methods: Recursive Marginal Quantization and Quantization via Fou-rier's transform.We will then estimate the price of vanilla call options for dierent dynamics of the underlying stochasticprocess, starting with Lévy processes (with and without jumps) and ending with stochastic volatilitymodels.For all those cases we will compare our prices with the ones calculated by the method based on the FFT.We will conclude with an analysis of our results and a comparison between the two methods.

Indice

0 Introduzione 1

1 La Quantizzazione 3

1.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Quantizzatori ottimi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2.1 Partizione di Voronoi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.2 La griglia ottima ΓN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.3 Esistenza e unicità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 La Quantizzazione quadratica: p = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.4 Metodi Numerici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4.1 Il metodo di Lloyd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.4.2 Metodi di discesa del gradiente - Il metodo di Newton Raphson . . . . . . . . . . 11

2 Option pricing: i modelli 13

2.1 Introduzione al pricing di derivati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2 Il modello di Black and Scholes: MBG . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Processi di Lévy: introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4 Processi di Lévy ad attivita nita: Kou . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5 Processi di Lévy ad attivita innita: Variance Gamma . . . . . . . . . . . . . . . . . . . 202.6 Processi a volatilità stocastica: Heston . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3 Quantizzazione Marginale Ricorsiva 23

3.1 Il metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2 Applicazione al modello di Black and Scholes . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2.1 Implementazione numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3 Applicazione al modello di Heston . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3.1 Implementazione numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4 Quantizzazione e trasformata di Fourier 35

4.1 Il metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1.1 Quantizzatori ottimi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1.2 Newton Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.1.3 Stima dell'errore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1.4 La scelta di p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4.2 Implementazione numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.2.1 Applicazione al modello di Black and Scholes . . . . . . . . . . . . . . . . . . . . 414.2.2 Applicazione al modello Kou . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.2.3 Applicazione al modello Variance Gamma . . . . . . . . . . . . . . . . . . . . . . 444.2.4 Applicazione al modello di Heston . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5 Conclusioni 47

A Codici Quantizzazione Marginale Ricorsiva 49

v

B Codici Quantizzazione con Fourier 61

Bibliograa 70

Elenco delle gure

1.1 Diagramma di Voronoi distribuzione Gaussiana. d = 2, N = 50. . . . . . . . . . . . . . 5

4.1 Quantizzatori p-ottimi, al variare di p, per Z ∼ N(0, 1) con N = 10 ssato. . . . . . . . 39

vi

Introduzione

I derivati sono strumenti nanziari complessi e, per la loro enorme diusione sui mercati nanziari, han-no ormai acquisito un ruolo di assoluta centralità nell'intera economia globale. Il problema del pricingdi derivati è dunque sempre più attuale e fondamentale.

Un'opzione (o derivato) è un contratto che dà al possessore il diritto, ma non l'obbligo, di acquistare(call) o vendere (put) il titolo sul quale si basa l'opzione (sottostante) ad un prezzo di esercizio ssato(strike price) e/o ad una determinata data. Nel nostro lavoro tratteremo opzioni Europee, ovvero conpossibilità di esercizio unicamente a scadenza (maturity).Il titolo sottostante, ovvero quello da cui dipende il valore dell'opzione, è modellizzato tramite un pro-cesso stocastico che può seguire dinamiche di diverso tipo. Tratteremo in particolare i seguenti casi:il moto browniano geometrico (GBM), un processo Lévy ad attività nita (Kou), un processo Lévyad attività innita (Variance Gamma) ed un modello a volatilità stocastica (Heston). Nel capitolo 2approfondiremo le caratteristiche di questi modelli.

Esistono ad oggi diverse tecniche numeriche per la stima del prezzo di titoli derivati, tra queste le piùutilizzate sono il metodo MonteCarlo, il metodo delle dierenze nite ed il metodo basato sulla FFT(Fast Fourier Transformation).L'applicazione dei metodi di quantizzazione in questo campo è più recente e costituisce tuttora un te-ma di ricerca. Vengono utilizzati attualmente due metodi dierenti, che proporremo nel nostro lavoro:Quantizzazione Marginale Ricorsiva (RMQ, capitolo 3) e Quantizzazione con trasformata di Fourier(capitolo 4).

Dato un segnale X, a valori in uno spazio continuo (quale ad esempio Rd ), la Quantizzazione permettedi ottenere un segnale, a valori in uno spazio discreto (x1...xN ), che approssima l'originale.L'insieme discreto di valori così ottenuto, ovvero la griglia N -dimensionale, è detto N -quantizzatore ead esso viene associato un vettore N -dimensionale contenente le probabilità associate a ciascun valorexi.Tema cruciale della quantizzazione è dunque la scelta della griglia ottimale, ovvero quella che permettela miglior approssimazione del segnale. Approfondiremo la teoria al riguardo nel capitolo 1.

La Quantizzazione ha origine negli anni '50, quando fu utilizzata per la prima volta per l'ottimizzazionedella compressione e trasmissione di segnali digitali. Solo all'inizio degli anni '90 si è iniziata ad utilizzarequesta metodologia in campo statistico, per la stima dei valori attesi.Dato un vettore aleatorio X, a valori in Rd, ed una funzione boreliana 1 f : Rd −→ R, una volta ottenutiil quantizzatore e i relativi pesi (probabilità) tramite quantizzazione, è infatti possibile approssimare ilvalore atteso di f(X) tramite il valore atteso di f(X), dove X rappresenta la versione quantizzata di X.

E[f(X)] =

N∑k=1

f(xk) ∗ P (Xk = xk) (1)

1Una funzione di Borel è una funzione misurabile rispetto alla σ-algebra di Borel. Quest'ultima è la più piccola σ-algebradenita su di un insieme dotato di struttura topologica che sia compatibile con la topologia stessa, ovvero che contengatutti gli aperti della topologia.

1

In questo lavoro andremo ad applicare la quantizzazione proprio per la stima di valori attesi. Nel nostrocaso X rappresenterà il prezzo del sottostante ed f il payo dell'opzione che andremo a prezzare.

Una caratteristica importante dei metodi numerici basati sulla quantizzazione è che, ssata la dimensionedi quantizzazione N , i quantizzatori per una data distribuzione di probabilità non cambiano. È quindipossibile sfruttarli nelle diverse applicazioni senza dover eseguire l'intera procedura di quantizzazione.

Tramite l'implementazione di codici numerici sarà possibile confrontare i risultati ottenuti tramite quan-tizzazione con quelli ottenuti tramite metodi più classici quali ad esempio il metodo basato sulla FFT,e di valutare quindi l'eettiva competitività dei metodi illustrati.

Capitolo 1

La Quantizzazione

In questo capitolo introduciamo le principali nozioni teoriche riguardanti la Quantizzazione, soerman-doci in particolar modo sulla quantizzazione quadratica (sezione 1.3), per poi passare, nella sezione 1.4,ad un'analisi dei metodi numerici utilizzati per la sua applicazione.

Faremo riferimento ai seguenti testi: [4], [6] e [17].

1.1 Introduzione

Con il termine quantizzazione si intende l'approssimazione di un segnale che ammette un insieme continuodi possibili valori, con un segnale che assume valori in un insieme discreto di dimensione N <∞.Quando lo spazio in cui X è costruito è uno spazio nito dimensionale, come ad esempio Rd, si parla diquantizzazione vettoriale.Nel nostro studio ci riferiremo sempre alla quantizzazione come quantizzazione vettoriale.

Consideriamo uno spazio di probabilità (Ω, F, µ) dove:

- Ω è uno spazio campionario

- F è una σ-algebra 1 su Ω

- µ: F−→[0,1] è una funzione di probabilità tale che µ(Ω) = 1

Consideriamo un vettore aleatorio X: (Ω, F, µ) −→ Rd, con distribuzione di probabilità µ

Denizione 1.1 (griglia di quantizzazione). Una griglia di quantizzazione, di livello N ≥ 1 è unsottoinsieme (nito) di Rd, d ≥ 1 della forma x1, ..., xN , di componenti distinte.

ΓN = x1, ..., xN.

Denizione 1.2 (funzione di quantizzazione). Una funzione di quantizzazione è una funzione di Borelq : Rd −→ ΓN che a un vettore X ∈ Rd associa X = q(X) ∈ ΓN .

1Si denisce σ-algebra, sull'insieme Ω, una famiglia F di insiemi con le seguenti proprietà:

- Ω ∈ F ;

- A ∈ F ⇒ AC ∈ F ;

- Aii famiglia numerabile di insiemi aperti t.c. Ai ∈ F∀i ⇒ Ui(Ai) ∈ F .

3

4 CAPITOLO 1. LA QUANTIZZAZIONE

Approssimare il vettore X tramite quantizzazione, signica approssimarlo con un valore X = q(X) chepuò assumere N valori in ΓN ⊂ Rd.

Questa procedura, in quanto procedura di approssimazione, comporta un errore.

Denizione 1.3 (errore di quantizzazione). Si denisce errore di quantizzazione il valore atteso delladierenza (in norma) tra il segnale originale e la sua quantizzazione.

Err = ‖X − q(X)‖. (1.1)

La scelta ottima di ΓN e q rappresenta dunque il fulcro della quantizzazione.

1.2 Quantizzatori ottimi

Il quantizzatore ottimo è, come è logico intuire, quello che minimizza l'errore.

Per quanto riguarda la scelta di q, ssato N , sia ΓN la griglia ottima, vale:

qN (X) quantizzatore ottimo di X ⇔ ‖X − qN (X)‖ = dist(X,ΓN ) := infi∈1,...,n ‖X − xi‖,

ovvero

qN (X) = ProjΓN (X) :=

N∑i=1

xi1Ci(ΓN )(X), (1.2)

dove Ci(ΓN )1≤i≤N è la partizione di Borel di (Rd, B(Rd)), anche detta partizione di Voronoi, generatadalla griglia ΓN .

1.2.1 Partizione di Voronoi

Anchè X = q(X) possa approssimare nel miglior modo possibile X, e quindi venga minimizzato l'er-rore, è essenziale scegliere in maniera ottimale la griglia ΓN , cui consegue la partizione Ci(ΓN ).

Data ΓN , Ci(ΓN )1≤i≤N è la partizione di Borel di Rd che soddisfa, per ogni i ∈ i, ..., N

Ci(ΓN ) ⊂ x ∈ Rd : |x− xi| = min

j=1:N|x− xj |. (1.3)

A seconda della norma 2 scelta, la partizione ottenuta varierà.Nel caso di norme euclidee 3, quali quelle che utilizzeremo nel lavoro di tesi, le celle di Voronoi hanno laparticolarità di essere insiemi convessi.Risulta dunque:

Ci(ΓN ) = X ∈ Rd : |x− xi| < min

1≤j≤N,j 6=i|x− xj |,

Ci(ΓN ) = x ∈ Rd : |x− xi| = min

1≤j≤N|x− xj |, (1.4)

Ci(ΓN ) ⊂ Ci(ΓN ).

2Dato uno spazio vettoriale V , su un campo K, una norma è una qualsiasi funzione N : V −→ R che rispetta:

- Denita positività (0 è il vettore nullo di V): N(v) ≥ 0 ∀v ∈ V , N(v) = 0⇔ v = 0;

- Positiva omogeneità: N(λv) = |λ|N(v) ∀λ ∈ K;

- Disuguaglianza triangolare (o subadditività): N(v1 + v2) ≤ N(v1) +N(v2) ∀v1, v2 ∈ V .

3dato x ∈ RN , la norma euclidea di x è ‖x‖2 =√∑N

i=1 x2i

1.2. QUANTIZZATORI OTTIMI 5

Figura 1.1: Diagramma di Voronoi distribuzione Gaussiana. d = 2, N = 50.

La partizione di Voronoi divide il dominio in N celle Ci(Γ) , ciascuna con un suo centroide xi.

Denizione 1.4 (centroide). Si denisce centroide l'elemento xi della griglia Γ.I centroidi (quantizzatori) xi sono i centri delle relative celle di Voronoi Ci(Γ) rispetto alla distribuzionedi probabilità µ.

I punti di una singola cella sono accumunati dalla caratteristica di essere più vicini al centroide di quellacella rispetto a tutti gli altri punti del piano non appartenenti alla stessa cella.

Nella gura 1.1 possiamo vedere un esempio di diagramma di Voronoi applicato alla distribuzione Gaus-siana con norma euclidea 4.

1.2.2 La griglia ottima ΓN

Per prima cosa deniamo l'errore in norma Lp, con p generico appartenente a [1,+∞), associato allagriglia di quantizzazione ΓN .

ep,N (X,Γ) := ‖X − X‖Lp = min1≤i≤N

‖X − xi‖Lp =

(∫Rd

min1≤i≤N

|ξ − xi|pµ(dξ)

)1

p. (1.5)

Dove abbiamo indicato, per semplicità di notazione: Γ = ΓN , X = XΓ = XΓN .

Cerchiamo dunque la griglia Γ, di dimensione massima pari ad N , che permetta di approssimare ladistribuzione µ di X in Lp

4Fonte: http://www.quantize.maths-.com


Denizione 1.5 (Distorsione). si denisce distorsione la funzione che, data una griglia Γ, vi associal'errore di quantizzazione corrispondente:

Dp : (x1, ..., xN ) 7−→ E[ min1≤i≤N

|X − xi|p] = ep,N (X,Γ)p,

Dp(Γ) = ep,N (X,Γ)p. (1.6)

Minimizzare l'errore, in relazione alla griglia Γ, signica dunque trovare il minimo della distorsione Dp

ad essa associata.

Possiamo cioè riscrivere il problema come segue:

ep,N (X,Γ) := infΓ⊂Rd,card(Γ)≤N

ep,N (X,Γ) = inf(x1,...,xN )∈(Rd)N

Dp(x1, ..., xN )1p . (1.7)

Ci interessa a questo punto studiare la dierenziabilità della funzione di Distorsione.

Un risultato importante al riguardo è il seguente, tratto da [17]:

Proposizione 1.2.1. Data una norma euclidea ‖ · ‖, dato p ∈ (1,+∞), la funzione Dp è dierenziabileper ogni ΓN = x1, ..., xN con componenti a due a due distinte e il gradiente è dato da:

∇Dp = p

(∫Ci(Γ)

xi − ξ|xi − ξ|

|xi − ξ|p−1µ(dξ)

)1≤i≤N

(1.8)

= p

(E[1X∈Ci(Γ)

xi −X|xi −X|

|xi −X|p−1]

)1≤i≤N

.

Introduciamo a questo punto una particolare categoria di quantizzatori, i quantizzatori stazionari.

Denizione 1.6 (quantizzatore stazionario). Dato X vettore aleatorio in Lp−1 con distribuzione µ, unagriglia di quantizzazione (o quantizzatore) Γ è p-stazionaria se:

a) la frontiera delle celle del diagramma di Voronoi ha misura nulla rispetto alla distribuzione di X

µ(∪1≤i≤N∂Ci(Γ)) = 0; (1.9)

b) i centroidi sono auto-consistenti, ovvero, ∀i ∈ 1, ..., N

E

[1X∈Ci(Γ)

xi −X|xi −X|

|xi −X|p−1

]= 0. (1.10)

Osservazione 1. I quantizzatori stazionari sono i punti critici della funzione di distorsione. Vale cioè:

Γ quantizzatore p-stazionario ⇔ ∇Dp(Γ) = 0.

Osservazione 2. I quantizzatori ottimi sono stazionari, il viceversa non sempre è vero.

1.2. QUANTIZZATORI OTTIMI 7

1.2.3 Esistenza e unicità

Avendo riscritto la ricerca del quantizzatore ottimo come un problema di minimo, risulta fondamentalestudiare, oltre alla convergenza dell'errore, l'esistenza del minimo e la sua unicità.

Proposizione 1.2.2. Esistenza e ConvergenzaSia p ∈ [1,+∞), sia X ∈ Lp, allora la funzione di distorsione Dp,N è nita ovunque su (Rd)N e :

a) Dp,N ha (almeno) un minimo nell'insieme di valori Γ∗ = x∗1, ..., x∗N;

b) se card(supp(µ)) > N , allora la corrispondente griglia Γ∗ ha dimensione N e, per ogni cella dellapartizione di Voronoi Ci(Γ

∗)1≤i≤N di Rd indotta da Γ∗, si ottiene:

P (X ∈ Ci(Γ∗)) > 0;

c) se card(supp(µ)) < N allora la sequenza N 7−→ ep,N (X) è strettamente decrescente e

limN→∞

ep,N = 0.

Per la dimostrazione si rimanda a [6].

Questo risultato garantisce che, come auspicabile, al crescere della dimensione N della griglia Γ, l'erroredi quantizzazione converge a zero.

Per quanto riguarda l'ordine di convergenza dell'errore di quantizzazione, invece, ci sono due risultatiteorici che illustreremo: il teorema di Zador e il lemma di Pierce.Il primo stabilisce l'ordine di convergenza, il secondo denisce un limite superiore non asintotico per ilquantizzatore ottimo (e viene utilizzato nella dimostrazione del teorema di Zador).

Teorema 1.2.3. Zador ([6])Sia X un vettore aleatorio in Rd tale che ∃ ν > 0 : E[Xp+ν ] < ∞, sia µ = fλd + λs la decomposizionedi Lebesgue di µ rispetto alla misura di Lebesue, dove λs rappresenta la sua parte singolare. Allora

limN→∞

N1d ep,N (X) = Qp(µ), (1.11)

dove Qp(µ) è la costante di Zador e si ha

Qp(µ) = Jp,d

(∫ d

R

fdd+p dλd

) 1p+ 1

d

= Jp,d‖f‖

1

pdd+p

∈ [0,+∞),

Jp,d = infN≥1

N1d ep,N (U [0, 1]d) ∈ (0,+∞).

Si rimanda a [6] per la dimostrazione.

Osservazione 3. Nel caso di norme euclidee, in dimensione 1 il valore di Jp,d è noto ed è pari a

Jp,1 =1

2p(p+ 1),

mentre

limd→+∞

Jp,d v

√d

2πe.

Il tasso di convergenza dell'errore è quindi O(N−1d ).


Lemma 1.2.4. Lemma di PierceSia ν > 0, ∃Kp,d,ν costante universale tale che, per ogni vettore aleatorio X: (Ω, F, µ) −→ Rd vale:

infΓ⊂Rd,card(Γ)≤N

‖X − X‖p 6 Kp,d,νσp,v(X)N−1d , (1.12)

doveσp,v(X) = inf

ξ∈Rd‖X − ξ‖p+ν ≤ +∞

rappresenta l' Lp-deviazione standard diX.

Abbiamo analizzato la convergenza dell'errore nella quantizzazione di un vettore aleatorio. Vediamo oracome si comporta l'errore quando, dati una funzione f , il vettore X e il suo quantizzatore X, vogliamoapprossimare E[f(X)] tramite E[f(X)].

Al riguardo presentiamo le seguenti proposizioni, introdotte in [5].

Proposizione 1.2.5. Sia Γ un quantizzatore stazionario, sia f una funzione di Borel su Rd, se f èconvessa allora vale:

E[f(X)] ≤ E[f(X)]. (1.13)

Proposizione 1.2.6. Sia Γ un quantizzatore stazionario, sia f una funzione di Borel su Rd, se f èLipschitz 5 e continua, allora vale:

|E[f(X)]− E[f(X)]| ≤ ‖f‖lip‖X − X‖2. (1.14)

Per maggiori dettagli si rimanda a [14]

Concludiamo questa sezione con un'importante risultato nel caso unidimensionale (d=1):

Proposizione 1.2.7. Sia µ è una funzione di densità log-concava, se X ha distribuzione AssolutamenteContinua (AC 6 ) allora esiste esattamente un unico quantizzatore stazionario e ottimo per ogni Nssato N ≥ 1.


Osservazione 4. la densità gaussiana è una funzione log-concava.

5Una funzione f denita su I ⊂ Rd a valori in Rd, si dice Lipschitz (o lipschitziana) se esiste L tale che ‖f(x1)−f(x2)‖ <L‖x1 − x2‖. La più piccola costante L per cui vale ciò è denominata costante di Lipschitz e viene indicata con ‖f‖lip.

6Una funzione f : R −→ R si dice AC se per ogni ε > 0 ∃δ(ε) t.c. ∀ successione di intervalli [ak, bk] disgiunti si ha:∑k(bk − ak) < δ(ε)⇒

∑k |f(bk)− f(ak)| < ε.

1.3. LA QUANTIZZAZIONE QUADRATICA: P = 2 9

1.3 La Quantizzazione quadratica: p = 2

In questa sezione tratteremo più dettagliatamente il metodo della quantizzazione nel caso p=2.

Per prima cosa osserviamo che l'eq. (1.8) , nel caso p=2, si può riscrivere come segue:

∇Dp = 2

(∫Ci(Γ)

(xi − ξ)µ(dξ)

)1≤i≤N

= 2(E[1X∈Ci(Γ)(xi −X)]

)1≤i≤N . (1.15)

Quindi un quantizzatore ΓN è stazionario se, ∀i ∈ 1, ..., N, si ha (autoconsistenza dei centoridi):

E[1X∈Ci(Γ)(xi −X)] = 0, (1.16)

che possiamo riscrivere anche come:

E[X1X∈Ci(Γ)]− xiµ(X ∈ Ci(Γ)) = 0, (1.17)

∫Ci(Γ)

ξµ(dξ)− xi∫Ci(Γ)

µ(dξ) = 0,

ovvero

xi =

∫Ci(Γ)

ξµ(dξ)

µ(Ci(Γ))= E[X|X ∈ Ci(Γ)]. (1.18)

La distorsione quadratica ha la forma:

D2(Γ) = E[|X − X|2] =

∫Rddist(ξ,Γ)2dξ =

N∑i=1

∫Ci(Γ)

|ξ − xi|2µ(dξ). (1.19)

A questo punto il problema di minimo da risolvere è il seguente:

trovare Γ∗ = x∗1, ..., x∗N tale che D2(x∗1, ..., x∗N ) ≤ D2(x1, ..., xN ) ∀(x1, ..., xN ) ∈ (Rd)N ,

dove la proposizione (1.2.2) garantisce l'esistenza di (almeno) una soluzione.

Procediamo ad un esempio pratico:

Esempio 1. d=1, µ = U [0, 1], N ssato

Ci(Γ) := [x−i , x+i ] con:

x−1 = 0 , x−i =xi−1 + xi

2∀i ∈ 2, ..., N;

x+i =

xi + xi+1

2∀i ∈ 1, ..., N − 1 , x+

N = 1.

Distorsione e Gradiente hanno la seguente forma:

D2(Γ) =∑Ni=1

∫ x+i

x−i

|ξ − xi|2dξ −→ ∇Dp = 2(∫ x+

i

x−i

(xi − ξ)µ(dξ))

1≤i≤N.


Il gradiente si può calcolare analiticamente e si ottiene:

∂D2

∂x1(Γ) = 1

8 (x1 + x2)(3x1 − x2),

∂D2

∂xi(Γ) = 1

4 (2xi − (xi+1 + xi−1))(xi+1 + xi−1) ∀i ∈ 2, ..., N − 1,

∂D2

∂xN(Γ) = 1

4 (3xN − xN−1 − 2)(2− (xN − xN−1)).

Ponendo a questo punto ∇D(Γ∗) = 0 otteniamo Γ∗ così denita:

Γ∗ = 2i− 1

2N, i = 1, ..., N;

mentre i pesi corrispondenti ai punti della griglia saranno tutti uguali e pari a

ωi =1

N, i = 1, ..., N,

con un errore di quantizzazione dato da

ep,N (X) =1

2(1− p)1pN

.

In questo esempio siamo riusciti ad ottenere la soluzione senza ricorrere a metodi numerici, nella maggiorparte dei casi però questo non è possibile.Introduciamo quindi la sezione seguente per analizzare le possibili soluzioni numeriche.

1.4. METODI NUMERICI 11

1.4 Metodi Numerici

Per la risoluzione del problema di minimizzazione esistono essenzialmente due categorie di metodi nu-merici: i metodi di punto sso e quelli basati sul gradiente.

Per quanto riguarda la prima categoria vedremo il metodo di Llyod, mentre per la seconda illustreremoil metodo di Newton Raphson.Entrambi questi metodi sono funzionali solo al caso monodimensionale (d = 1) e verranno trattati nelcontesto della quantizzazione quadratica (p = 2), negli altri casi esistono altri metodi basati sul gradientequali ad esempio il Competitive Learning Vector Quantization (CLQV).

Vedremo comunque come sarà possibile ricondurci al caso monodimensionale anche con i modelli avolatilità stocastica, potendo utilizzare quindi il metodo di Newton Raphson.

1.4.1 Il metodo di Lloyd

Il metodo di Llyod fu il primo ad essere implementato per la risoluzione del problema di minimizzazione1.3, nel caso p=2, d=1.

Si tratta di un metodo iterativo che sfrutta la proprietà di autoconsistenza (eq. 1.18 ) dei quantizzatoristazionari.

Sia XΓ(0) un quantizzatore iniziale di dimensione N , si aggiorna la griglia ad ogni passo in questo modo:

XΓ(k+1) = E[X|XΓ(k)] =(E[X|X ∈ Ci(XΓ(k))]

)1≤i≤N

=

E[X1X∈∈Ci(XΓ(k))

]µ(X ∈ Ci(XΓ(k)))

1≤i≤N

. (1.20)

Valgono le seguenti proposizioni, per la cui dimostrazione si rimanda a [6].

Proposizione 1.4.1. Sia e2,k = ‖X − XΓ(k)‖ l'errore di quantizzazione al passo k, la successionee2,k è una successione decrescente.

Proposizione 1.4.2. Sia µ una funzione strettamente log-concava, allora la successione XΓ(k) convergeal quantizzatore ottimo.

1.4.2 Metodi di discesa del gradiente - Il metodo di Newton Raphson

I metodi di discesa del gradiente consentono di trovare il minimo di una funzione.La tecnica di discesa secondo gradiente si basa sul fatto che, per una data funzione f(x), la direzione dimassima discesa in un assegnato punto x corrisponde a quella determinata dall'opposto del suo gradientein quel punto (−∇f(x)).Questa scelta per la direzione di discesa garantisce che la soluzione tenda a un punto di minimo di f .Il metodo del gradiente prevede dunque di partire da una soluzione iniziale x0 e di procedere iterativa-mente aggiornandola come:

xk+1 = xk − αk∇f(xk). (1.21)

Nel nostro contesto, la funzione da trattare è la distorsione D2 e la variabile è Γ. Quindi:

Γk+1 = Γk − αk∇D2(Γk).

L'algoritmo di Newton Raphson fa parte di questa famiglia.Lo utilizzeremo per la ricerca degli zeri della distorsione, nonchè i quantizzatori stazionari e, sottoopportune condizioni, ottimi.


In questo contesto il coeciente αk corrisponde all'inverso della matrice Hessiana calcolata al passoprecedente.In formula:

Γk+1 = Γk − [HD2(Γk)]−1∇D2(Γk). (1.22)

Questa procedura iterativa converge a Γ∗ tale che ∇D2(Γ∗) = 0.

Esempio 2. d = 1, µ = N(0, 1)

Deniamo, come nel caso della distribuzione uniforme (Esempio 1), le celle Ci(Γ) := [x−i , x+i ] con:

x−1 = −∞ , x−i =xi−1 + xi

2∀i ∈ 2, ..., N;

x+i =

xi + xi+1

2∀i ∈ 1, ..., N − 1 , x+

N = +∞.

Deniamo inoltre erf(x) =2√π

∫ x0e−t

2

.

La distorsione D2(Γ) ha la seguente forma:

∑Ni=1

((1 + x2

i )[erf(x+i )− erf(x−i )]− 1√

2π

[x+i e−

(x+i

)2

2 − (x−i )e−(x

−i

)2

2

]+ 2√

2πxi

[e−

(x+i

)2

2 − e−(x

−i

)2

2

]).

Possiamo calcolare il gradiente e risulta:

∂D2

∂xi(Γ) = xi

[erf(x+

i )− erf(x−i )]

+ 1√2π

[e−

(x+i

)2

2 − e−(x

−i

)2

2

].

Inne, anche l'Hessiana è calcolabile analiticamente e si ha:

∂2D2

∂xi∂xi−1(Γ) = 1

4√

2π(xi − xi−1)e−

(x−i

)2

2 ,

∂2D2

∂x2i

(Γ) = erf(x+i )− erf(x−i )− 1

4√

2π(xi+1 − xi)e−

(x+i

)2

2 − 14√

2π(xi − xi−1)e−

(x−i

)2

2 ,

∂2D2

∂xi∂xi+1(Γ) = − 1

4√

2π(xi+1 − xi)e−

(x+i

)2

2 .

In questo caso abbiamo quindi tutti gli strumenti per applicare il metodo di Newton Raphson (eq. 1.22)per il calcolo della soluzione.

Osservazione 5. Notiamo come la scelta della griglia iniziale rappresenti un problema cruciale.Griglie ottime per dimensioni N = 1÷ 5999 sono scaricabili a: http://www.quantize.maths-fi.com.

Capitolo 2

Option pricing: i modelli

In questo capitolo introduciamo alcune nozioni generali riguardanti le opzioni e l'option pricing (sezione2.1), per poi analizzare nel dettaglio i singoli modelli che verranno trattati nel lavoro di tesi.

Nello specico, introdurremo i seguenti modelli:

- Black and Scholes;

- Kou;

- Variance Gamma;

- Heston.

Principali fonti per questo capitolo saranno [10] e [7]. In particolare si farà riferimento a [10] per quantoriguarda le sezioni 2.1 e 2.2 e a [7] per le altre.

2.1 Introduzione al pricing di derivati

Scopo primario di un modello di pricing è quello di catturare lo stato del mercato dei derivati ad un datoistante. A questo scopo sono stati introdotti nel tempo una serie di modelli, con dierenti caratteristiche,per rispondere alle osservazioni dei dati di mercato.

Per prima cosa riprendiamo la denizione di opzione:

Denizione 2.1 (Opzione). Un'opzione (o derivato) è un contratto che da al possessore il diritto, manon l'obbligo, di acquistare (call) o vendere (put) il titolo sul quale si basa l'opzione (sottostante) adun prezzo di esercizio ssato (strike price) e/o ad una determinata data. Nel caso di opzioni Europee lapossibilità di esercizio si ha unicamente a scadenza (maturity).

A seconda degli istanti di possibile esercizio e/o della caratterizzazione del payo dell'opzione stessa, sidistinguono le diverse tipologie di opzioni.In particolare:

- opzioni europee (vanilla): possibilità di esercizio unicamente a scadenza, non dipende dal valoredel sottostante prima di T (maturity);

- opzioni americane : possibilità di esercizio in ogni istante temporale tra 0 e T , non dipende dalvalore del sottostante prima di T (maturity);

- opzioni esotiche: sono caratterizzate dal fatto di essere path dependent, ovvero il payo dell'opzionedipende dal valore del sottostante St anche in istanti di tempo precedenti a T;

13

14 CAPITOLO 2. OPTION PRICING: I MODELLI

Come anticipato, nel nostro lavoro tratteremo opzioni di tipo europeo e, in particolare, call europee.Tali opzioni hanno valore, nell'istante t = 0, determinato dalla seguente espressione:

Call = e−rT (ST −K)+ = e−rTmax((ST −K), 0), (2.1)

dove:

- K rappresenta lo strike price;

- S rappresenta il sottostante;

- T rappresenta la maturity (scadenza);

- r rappresenta il tasso d'interesse risk free del mercato.

Osservazione 6. Il termine (e−rT ) è detto fattore di sconto e serve per attualizzare il prezzo (calcolatoin t = T ) all'istante di tempo corrente (t = 0), tenendo conto del tasso d'interesse risk free del mercatoe della maturity dell'opzione. Il payo (valore) denito dalla (2.1) è per questo detto payo scontato.È possibile stimare il valore di un'opzione anche ad un istante t generico, in tal caso il suo valore saràdato dalla seguente espressione: Call(t) = e−r(T−t)(ST −K)+.

Il sottostante S cui fa riferimento l'opzione può essere modellizzato tramite diversi processi, ciascunocon una propria dinamica. A seconda del processo in questione, vanno svolti conti appositi per ricavareil valore di S a scadenza e poter quindi successivamente determinare il valore (prezzo) della Call.

Il primo modello che illustreremo è quello di Black and Scholes.Si tratta di un modello semplice e molto utilizzato che però presenta, come vedremo nella sezione 2.2,ipotesi poco realistiche. Ad esempio in tale modello la volatilità è una costante, caratteristica pococonforme all'eettivo andamento dei mercati.

Nel tempo sono stati dunque introdotte altre tipologie di modelli, per rispondere alle osservazioni sul-l'evoluzione dei prezzi dei titoli.

In questo lavoro tratteremo, oltre al modello di Black and Scholes, le seguenti tipologie di processi:

- processi di Lévy con salti, ad attività nita: sezione 2.4;

- processi di Lévy con salti, ad attività innita: sezione 2.5;

- processi a volatilità stocastica: sezione 2.6.

2.2. IL MODELLO DI BLACK AND SCHOLES: MBG 15

2.2 Il modello di Black and Scholes: MBG

Nel 1973, Fisher Black e Myron Scholes proposero un modello di valutazione delle opzioni 1 che vennesuccessivamente perfezionato e ampliato nello stesso anno da Merton 2. Tale modello valse ai due autorioriginari il premio Nobel per l'economia nel 1997.Ancora oggi, nonostante i suoi limiti che illustreremo nel seguito di questo capitolo, questo modellorappresenta un pilastro fondamentale nell'economia nanziaria.

Il modello di Black and Scholes prevede un'economia semplicata, che consiste in due soli asset (titoli):un'obbligazione priva di rischio B(t) con tasso di interesse r e scadenza T e un titolo rischioso S(t),sottostante a un titolo nanziario (derivato) che si vuole prezzare.

I due titoli seguono il seguente andamento:

St = S0e(r−σ2

2 )t+σWt , (2.2)

S(0) = S0;

Bt = B0ert (2.3)

B(0) = 1;

La dinamica di (St) è quella del moto browniano geometrico (MBG), dove:

- r è il tasso risk-free dell'obbligazione e corrisponde al termine di drift dell'EDS in rispetto delprincipio di non arbitraggio 3;

- σ è il termine di diusione dell'EDS e rappresenta la volatilità del sottostante;

- Wt è il Moto Browniano standard 4 continuo e reale sullo spazio di probabilità (Ω, F, P );

(Bt), invece, evolve nel tempo in maniera deterministica.

In questo lavoro ci concentreremo sul titolo rischioso S, in quanto sottostante dell'opzione Call Europeache andremo a prezzare.

L'EDS associata all'eq (2.2) è la seguente:

dSt = rStdt+ σStdWt, (2.4)

S(0) = S0.

Dall'eq. (2.2) possiamo facilmente calcolare il valore atteso e la varianza del processo St e risulta:

E[S(t)] = S0ert, (2.5)

V ar[S(t)] = S20e

2rt(eσ2t − 1). (2.6)

1Black, F. e Scholes, M. (1973), The Pricing of Options and Corporate Liabilities, Journal of Political Economy 81 (3),637-654.

2Merton, R. (1973), Theory of Rational Option Pricing, Bell Journal of Economics and Management Science 4 (1),141-183.

3Un arbitraggio è un'operazione che consiste nell'acquistare un bene o un'attività nanziaria su un mercato rivendendolosu un altro mercato, sfruttando le dierenze di prezzo al ne di ottenere un protto. Il principio di non arbitraggio è unprincipio base della matematica nanziaria e il suo rispetto è vincolante nel pricing degli strumenti derivati.

4Un moto Browniano standard (o processo di Weiner) è un processo stocastico gaussiano (Wt)t≥0 ∼ N(0, t) conincrementi stazionari, indipendenti e con distribuzione gaussiana: (Wt −Ws ∼ N(0, t− s).


Le ulteriori ipotesi su cui si basa il modello di Black and Scholes sono le seguenti:

- perfetta liquidità del mercato 5;

- possibilità di vendita allo scoperto 6 dei titoli;

- assenza di possibilità di arbitraggio;

- assenza di costi di transazione e/o tassazioni di altro tipo;

- scambio sul mercato in tempo continuo di sottostante e strumento derivato;

- medesimo tasso risk-free r per tutte le scadenze.

Per questo modello, esiste una formula chiusa per il prezzo della Call europea ed è la seguente [10]:

CBS(S(0), 0) = S(0)N(d1)−Ke−rTN(d2), (2.7)

dove N(·) indica la funzione di ripartizione della Normale standard e gli argomenti d1 e d2 sono datidalle seguenti espressioni:

d1 =ln(S(0)

K ) + (r − σ2

2 )T

σ√T

,

d2 =ln(S(0)

K )− (r − σ2

2 )T

σ√T

.

Per la dimostrazione di questo risultato si rimanda a [10].

In Matlab esiste una funzione apposita che, ricevendo in input i parametri del modello, restituisce ilvalore della Call: blsprice. Nel lavoro di tesi utilizzeremo questa funzione per stimare il prezzo esatto,da utilizzare come riferimento per valutare l'accuratezza delle nostre stime.

Il modello di mercato che sta alla base dell'analisi svolta da Black and Scholes ha avuto un note-vole successo applicativo sia per la sua semplicità che per la possibilità di ottenere formule esplicite divalutazione di molti titoli derivati; alcune assunzioni alla base del modello sono tuttavia poco realistiche.

Tra le problematiche di questo modello troviamo le seguenti:

- Le densità empiriche dei rendimenti presentano code più spesse della distribuzione normale edun'asimmetria negativa;

- la volatilità dei rendimenti è autocorrelata nel tempo ed inne si è osservata una correlazionenegativa tra rendimenti e volatilità (legerave eect).

- Dall'analisi delle serie storiche di mercato si osserva che σ = σ(K,T ), dierentemente da quantoipotizzato da Black and Scholes. Invertendo la formula (2.7) per ricavare σ, inserendo i prezzi dimercato delle opzioni, e rappresentando gracamente i prezzi rispetto allo strike K, si ottiene ilcosì detto smile di volatilità. Questa caratteristica forma è dovuta al fatto che la volatilità tendead essere maggiore per le opzioni con moneyness S

K molto alta o molto bassa.

Queste regolarità del comportamento dei prezzi evidenziano i limiti di fondo del modello lognormale nelcatturare gli aspetti statistici fondamentali delle serie storiche dei prezzi.

5un mercato à perfettamente liquido quando si possono scambiare i titoli, in quantità arbitraria, in ogni istante6La vendita allo scoperto di un titolo è un'operazione nanziaria che consiste nella vendita, eettuata nei confronti di

uno o più soggetti terzi, di titoli non direttamente posseduti dal venditore.

2.3. PROCESSI DI LÉVY: INTRODUZIONE 17

2.3 Processi di Lévy: introduzione

In questa sezione introduciamo i modelli di salto; ovvero che presentano l'introduzione dei salti aleatorinella dinamica del prezzo del sottostante all'opzione.Questa tipologia di modelli è stata introdotta in risposta alle osservazioni dei dati di mercato che hannoevidenziato traiettorie discontinue per i prezzi, i quali presentano appunto, saltuariamente, dei salti.

Esponiamo in questa sezione alcune nozioni generali riguardanti i processi di Lévy con salti, limitandocia quanto di interesse per il lavoro di tesi.

Per quanto presente in questa sezione e per eventuali approfondimenti si fa riferimento a [7].

Denizione 2.2 (Processo di Lévy). Un processo stocastico - cadlag 7 - (Xt)t>0 denito su (Ω, F, P ),a valori in Rd con X(0) = X0, si denisce processo di Lévy se rispetta le seguenti proprietà:

a) incrementi indipendenti: per ogni sequenza crescente t0, ..., tn le variabili aleatorie (Xt0), (Xt1 −Xt0), ..., (Xtn −Xtn−1

) sono indipendenti;

b) incrementi stazionari: la legge di (Xt+h −Xt) è indipendente da t ;

c) continuità stocastica: ∀ε > 0, limh→0 P (|Xt+h −Xt| > ε) = 0.

Denizione 2.3 (misura di Lévy). Dato un processo di Lévy (Xt)t>0 denito su (Ω, F, P ), a valoriin Rd, si denisce misura di Lévy del processo, la funzione ν : B(Rd) −→ R così denita sugli insiemiA ∈ B(Rd):

ν(A) = E[#t ∈ [0, 1] : ∆Xt 6= 0 ∧∆Xt ∈ A]; (2.8)

ν(A) rappresenta cioè il numero atteso di salti, per unità di tempo, la cui misura appartiene ad A.

Denizione 2.4 (misura di salto). Jx : Rd × [0,+∞) −→ R, Jx(B) = #(t,∆t) ∈ B è detta misuradi salto ed è la misura di un processo compound poisson con intensità di misura: ν(dx)dt.

Esempio 3. Un esempio di processo di Lévy è il processo Compound Poisson.Un compound Poisson è un processo Xt =

∑Nti=1 Yi =

∫ t0

∫RxJx(ds × dx), dove Nt è un processo di

Poisson 8, Yi sono variabili aleatorie e Jx è la misura di salto.

Introduciamo a questo punto alcune proprietà fondamentali di questa tipologia di processi.

Sia (Xt)t>0 un processo di Lévy denito su (Ω, F, P ), a valori in Rd, allora valgono le seguenti:

Teorema 2.3.1. (Lévy Ito decomposition [7])Xt si può scomporre come segue:

Xt = γt+Bt +X`t + lim

ε→0Xεt , (2.9)

dove:

7Una funzione f : [0, T ] −→ Rd si dice Cadlag se valgono le seguenti:

- f(t−) = lims→t,s<t esiste;

- f(t+) = lims→t,s>t = f(t).

Ogni processo continuo è cadlag, il viceversa non è necessariament vero.8Data una sequenza di variabili aleatorie (τi)i≥1, indipendenti e identicamente distribuite con distribuzione esponenziale

di parametro λ, data Tn =∑n

i=1 τi, si denisce processo di Poisson, di parametro λ, il processo Nt =∑

n≥1 1t>Tn Per

magigori dettagli si rimanda a [7].


- X`t =

∫ t0

∫|x|≥1

xJx(ds× dx) è la componente dei grandi salti (|x| ≥ 1);

- Xεt =

∫ t0

∫ε≤|x|≥1

xJx(ds× dx)− ν(dx)ds è la componente dei piccoli salti;

- γ ∈ Rd ;

- Bt ha matrice di varianza - covarianza A ∈ Rd x Rd.

Denizione 2.5. (tripletta caratteristica) La tripletta (γ,A, ν) à detta tripletta caratteristica delprocesso ed identica univocamente un processo di Lévy.

Proposizione 2.3.2. La funzione caratteristica 9 di (Xt)t>0 ΦXt(z), è moltiplicativa: ΦXt+s(z) =ΦXt(z)ΦXs(z).

Dimostrazione. Per la proprietà b dei processi di Lévy (denizione 2.2) vale:ΦXt+s(z) = E[eizXt+s ] = E[eiz(Xt+s)−(Xt)]E[eizXt ] = E[eizXs ]E[eizXt ] = ΦXs(z)ΦXt(z).

Proposizione 2.3.3. Esiste ΨXt : Rd −→ R continua tale che ΦXt(z) = etΨXt (z).

Denizione 2.6 (Esponente Caratteristico). ΨXt così denito è detto Esponente Caratteristico delprocesso Xt.

Un processo stocastico è univocamente identicato dalla sua funzione caratteristica, vale infatti ilseguente teorema:

Teorema 2.3.4. ([16])Date due variabili aleatorie X e Y vale che:

ΦX(z) = ΦY (z) per ogni z ⇔ X e Y hanno la distribuzione di probabilità.

Si può trovare una dimostrazione di questo risultato in [16].

Vedremo nel capitolo 4 come la funzione caratteristica di un processo può essere utilizzata per la stimadel prezzo del titolo modellizzato dallo stesso.

Esempio 4. Il Moto Browniano geometrico è un processo di Lévy ed è caratterizzato dal seguenteesponente caratteristico:

Ψ(u) = −σ2

2iu− σ2

2u2. (2.10)

I processi di Lévy si dividono essenzialmente in due categorie:

- processi di tipo jump diusion:sono caratterizzati dalla presenza di una componente Browniana e hanno distribuzione dei saltinota; fa parte di questa categoria il modello Kou (attività nita);

- processi di puro salto:si muovono essenzialmente per salti e risultano d'interesse solo quando ad attività innita; a questacategoria appartiene il processo Variance Gamma.

Approfondiremo meglio queste due categorie nelle prossime sezioni del presente capitolo.

9Dato un processo (Xt)t>0 denito su (Ω, F, P ), a valori in Rd, si denisce funzione caratteristica del processo lafunzione Φ : ΦXt (z) = E[eizXt ].

2.4. PROCESSI DI LÉVY AD ATTIVITA FINITA: KOU 19

2.4 Processi di Lévy ad attivita nita: Kou

In questa sezione andiamo ad approfondire le caratteristiche dei processi di Lévy ad attivita nita ditipo Jump Diusion, andando a descrivere in maniera più dettagliata il processo Kou.Anche per questa sezione, il testo di riferimento sarà [7].

Un processo di tipo Jump Diusion ha la seguente forma:

Xt = bt+ σWt +

Nt∑i=1

Yi, (2.11)

dove:

- Wt è il moto Browniano standard;

- b rappresenta il drift della parte continua del processo XCt = bt+ σWt;

- σ rappresenta la deviazione standard di XCt ;

- (Yi)i=1:Nt sono variabili aleatorie con distribuzione f nota (la scelta di f determina il modello);

- Nt è un processo di Poisson: Nt ∼ Poisson(λt), dove λ =∫ +∞−∞ ν(dx) è la misura di Lévy del

processo X.

Nel processo Kou le variabili (Yi)i=1:Nt seguono una distribuzione del tipo esponenziale asimmetrica,ovvero, per ogni i = 1 : Nt:

Yi ∼ f(y), (2.12)

con:f(dy) = [pλ+e−λ

+y1y>0 + (1− p)λ−eλ−y1x<0]dy, (2.13)

dove:

- λ+ governa il decadimento delle code per i salti positivi;

- λ− governa il decadimento delle code per i salti negativi;

- p ∈ [0, 1] è la probabilità di avere un salto positivo.

La densità del processo Xt è invece:

ν(x) = [pλλ+e−λ+x1x>0 + (1− p)λλ−e−λ

−x1x<0], (2.14)

dove λ è l'intensità di salto.

Un processo di tipo Kou è quindi identicato da 5 parametri: σ, p, λ, λ+ e λ−.

Allo stesso modo è possibile identicare il processo tramite il suo esponente caratteristico, il quale ha laseguente espressione:

Ψ(u) = −σ2u2

2+ ibu+ iuλ

[p

λ+ − iu+

1− pλ− + iu

]. (2.15)


2.5 Processi di Lévy ad attivita innita: Variance Gamma

Il processo Variance Gamma fa farte dei processi ottenuti tramite subordinazione browniana.In particolare viene utiizzato un processo tempered stable come subordinatore.

Per meglio chiarire tale espressione, introduciamo alcuni concetti teorici aggiuntivi riguardanti i processidi Lévy, per i cui approfondimenti si rimanda a [7].

Denizione 2.7. (Subordinatore) Si denisce subordinatore un processo di Lévy Xt per cui valgono leseguenti (equivalenti) proprietà:

a) Xt ≥ 0 ∀t ≥ 0 quasi certamente;

b) Xt > 0 per qualche t > 0 quasi certamente;

c) Xt non decrescente quasi certamente;

d) (γ,A, ν) tripletta caratteristica di Xt con: A = 0,∫|x|<1

xν(dx) < +∞, ν(−∞, 0) = 0, γ ≥ 0.

Osservazione 7. La proprietá d) indica in particolare che il processo é a variazione nita, con probabilitánulla di avere salti negativi e con drift positivo. Generalmente si indica la tripletta caratteristica di unsubordinatore con (b, 0, ρ).

Proposizione 2.5.1. (Subordinazione Browniana)Sia St un subordinatore con esponente di Laplace 10 `(u), siaWt un processo di Wiener tale cheWt ⊥ St;dato Yt = µt + σWt processo di Lévy, Xt = YSt = µSt + σWSt è anch'esso un processo di Lévy. Inparticolare vale:

ΨX(u) = `(ΨY (u)). (2.16)

Denizione 2.8. (Subordinatore tempered stable) Si denisce Subordinatore tempered stable un su-bordinatore caratterizzato dal seguente esponente di Laplace:

`(u) = cΓ(−α) [(λ− u)α − λα] Γ(z) =

∫ +∞

0

tz−1e−tdt, (2.17)

dove:

- c modula i salti;

- λ controlla i grandi salti;

- α determina il peso dei piccoli salti.

Denizione 2.9. (processo Variance Gamma) Un processo Variance Gamma é un processo della formaXt = YSt dove St é un subordinatore tempered stable con α = 0.

Proposizione 2.5.2. L'esponente caratteristico di un processo di tipo Variance Gamma è cosí denito:

Ψ(u) = −1

klog

(1 +

u2σ2k

2− iθku

). (2.18)

10dato un subordinatore con tripetta caratteristica (b, 0, ρ), il suo esponente di Laplace è così denito: `(u) = bu +∫+∞0 (eux − 1)ρ(dx).

2.6. PROCESSI A VOLATILITÀ STOCASTICA: HESTON 21

2.6 Processi a volatilità stocastica: Heston

I modelli a volatilità stocastica sono caratterizzati dal fatto che la volatilità non è più una costante o,più in generale, una funzione deterministica, bensì è essa stessa un processo stocastico.Questa scelta permette di riprodurre lo smile di volatilità tipicamente osservato nei mercati 11.

Il problema del pricing di derivati nel contesto dei modelli a volatilità stocastica è stato studiato a fondoin letteratura, a partire da Hull and White (1987), Chensey e Scott (1989), Stein and Stein (1991) edHeston (1993). Sono stati nel tempo introdotti numerosi modelli a volatilità stocastica, con o senza lapresenza di salti.

In questo lavoro di tesi andremo ad analizzare il modello di Heston [18].Tale modello rientra nella categoria dei così detti modelli ani, ovvero esiste una formulazione esplicitadella funzione caratteristica per il log-price. Questo ci permetterà, da un lato, di applicare il metodobasato sulla FFT per ssare il nostro benchmark e, dall'altro, di applicare la tecnica che combina quan-tizzazione e trasformata di Fourier, nel capitlo 4.

Il modello di Heston è caratterizzato dalla seguente coppia di EDS:per il (log)titolo Xt = St

S0:

dX(t) = X(t)rdt+X(t)√Vt(√

1− ρ2dWXt + ρdWV

t ), (2.19)

X(0) = X0 = log(S0) > 0;

per la volatilità V :

dV (t) = k(θ − V (t))dt+ ε√V tdWV

t , (2.20)

V (0) = V0 > 0.

dove:

- r rappresenta il tasso d'interesse privo di rischio;

- θ > 0 rappresenta la media di V ;

- ε > 0 è il parametro che determina la volatilità del processo V ;

- k > 0 rappresenta il tasso di reversione a θ ;

- WXt e WV

t sono due processi di Wiener scalari e indipendenti;

- ρ ∈ [−1, 1].

Il modello di Heston fu il primo in letteratura a introdurre una correlazione ρ tra il processo Xt e la suavolatilità Vt, rispondendo alle osservazioni di mercato che suggerivano una dipendenza tra l'andamentodella volatilità del prezzo di un titolo e quello del prezzo stesso.

Il processo caratterizzante V (t) è noto in letteratura come CIR 12 ed è caratterizzato da una dinamica di"mean reversion". Infatti il termine di drift k(θ−V (t)) fa tendere V (t) verso il valore θ, che rappresentaappunto la varianza media media a lungo termine del processo Xt, con una velocità di mean reversionpari a k.

11Dall'analisi delle serie storiche di mercato si osserva che σ = σ(K,T ), dierentemente da quanto ipotizzato da Blackand Scholes. Invertendo la (2.7) per ricavare σ, inserendo i prezzi di mercato delle opzioni, e rappresentando gracamentei prezzi rispetto allo strike K, si ottiene il così detto smile di volatilità. Questa caratteristica forma è dovuta al fatto chela volatilità tende ad essere maggiore per le opzioni con moneyness (S/K) molto alta o molto bassa.

12introdotto da Cox, Ingersoll e Ross, 1985


Osservazione 8. il processo che guida l'andamento di V non è lineare e non garantisce, senza ulteriorivincoli, la positività del processo, ovvero della varianza di X, la cui negatività non avrebbe senso inquesto contesto.

Per ovviare a questo problema introduciamo la seguente condizione, detta anche "Feller condition" [13]:

2kθ

ε2< 1, (2.21)

la quale garantisce, congiuntamente alla positività dei parametri k, θ, ε, la positività (stretta) del pro-cesso V .

Lo schema di Eulero associato al modello di Heston è il seguente:

Xi+1 = Xi + rXi∆t+Xi

√Vi(√

1− ρ2ZXi + ρZYi )√

∆t, (2.22)

V (i+ 1) = Vi + k(θ − Vi)∆t+ ε√ViZ

Vi

√∆t, (2.23)

dove (ZXi )i ∼ N(0, 1) e (ZVi )i ∼ N(0, 1) sono variabili normali standard indipendenti.

La funzione caratteristica del modello di Heston vale invece:

ΦX(u) = eC+DV0+iulog(X0), (2.24)

dove:

C =kθ

ε2

((k − iuρε− d(u))T − 2log

(c(u)− e−d(u)T

c(u)− 1

)),

D =k − iuρε+ d(u)

ε2

(1− e−d(u)T

c(u)− e−d(u)T

),

d(u) =√

(iuρε− k)2 + iuε2 + u2ε2,

c(u) =k − iuρε+ d(u)

k − iuρε− d(u).

Capitolo 3

Quantizzazione Marginale Ricorsiva

In questo capitolo introduciamo la quantizzazione marginale ricorsiva come metodo di pricing.

Introdurremo i concetti teorici fondamentali (sezione 3.1), per poi passare alle applicazioni nei modellidi Black and Scholes (sezione 3.2) e di Heston (sezione 3.3).

Per questo capitolo si è fatto riferimento a [4] e [5].

3.1 Il metodo

La quantizzazione marginale ricorsiva è un metodo "model free" che richiede unicamente la conoscen-za dello schema di Eulero (o Milstein) per l'equazione stocastica che governa l'evoluzione del sottostante.

Consideriamo lo spazio di probabilità (Ω, F, µ).Consideriamo un processo diusivo di Markov V continuo nel tempo, la cui evoluzione sia governatadalla seguente equazione dierenziale stocastica (EDS):

dVt = b(t, Vt)dt+ a(t, Vt)dWt, (3.1)

V (0) = V0,

dove:

- b(t, Vt) rappresenta la componente di diusione del processo;

- a(t, Vt) rappresenta la volatilità;

- Wt è un moto browniano standard;

- V0 > 0 ∈ R è noto e rappresenta il valore iniziale di V .

Fissiamo l'orizzonte temporale [0, T ] e la griglia di discretizzazione temporale 0 = t0, t1, ..., tM = T ,dove dt = T

M = (tk+1 − tk) è lo step di discretizzazione e ∀k = 0, ...,M si ha tk = kdt.

Lo schema di Eulero per il processo V risulta essere:

Vtk+1= Vtk + b(tk, Vtk)dt+ a(tk, Vtk)

√dtZk, (3.2)

Vt0 = V0 = V0,

dove (Zk)k ∼ N(0, 1) sono variabili aleatorie indipendenti con distribuzione normale standard.

23

24 CAPITOLO 3. QUANTIZZAZIONE MARGINALE RICORSIVA

Osservazione 9. Segue:√dtWk ∼ N(0, dt) ∀k = 0, ...,M .

Lemma 3.1.1. Dall'eq. (3.2) si evince che ∀k = 0, ...,M − 1 condizionatamente all'evento Vk = x,abbiamo:

(Vtk+1|Vtk = x) ∼ N(mk(x), σ2

k(x)), (3.3)

dove:

mk = x+ b(tk, Vtk)dt,

σ2k = a2(tk, Vtk)dt.

L'idea alla base della Quantizzazione Marginale Ricorsiva, introdotta da Pagès and Sagna nel 2014 [5],è di quantizzare il processo stocastico V dell'equazione (3.1) tramite approssimazione, lavorando sullevariabili aleatorie (marginali) Vtk dell'equazione (3.2).Applicheremo infatti la quantizzazione vettoriale alle singole variabili mono dimensionali Vtk , ottenendoi quantizzatori stazionari corrispondenti ∀k = 1, ...,M .

Osservazione 10. La quantizzazione sullo schema di Eulero porta con se un errore di approssimazione,rispetto alla reale distribuzione di probabilità delle variabili, ad ogni passo temporale.

Si parla di quantizzazione marginale ricorsiva perchè, per la stima di V tk non utilizzeremo la distribu-zione di Vtk+1

, generalmente sconosciuta, ma quella di Vtk+1|Vtk , nota dal lemma (3.1.1).

Notazione: per semplicità indicheremo con Vtk = V tk il quantizzatore di Vtk .

Osservazione 11. Vt0 = V0 , ovvero la griglia iniziale Γ0 consiste in un unico punto.

Deniamo a questo punto la distorsione quadratica (eq. 1.19) al passo k+ 1, relativa alla variabile Vtk+1

:

Dk+1(ΓNk+1) =

Nk+1∑i=1

∫Ci(Γ

Nk+1 )

(νk+1 − xNk+1

i )2P (Vtk+1∈ dνk+1), (3.4)

dove ΓNk+1 = xNk+1

1 , xNk+1

2 , ..., xNk+1

Nk+1 è l'Nk+1-quantizzatore al tempo tk+1 relativo a Vtk+1

.

Come già osservato la distribuzione di Vtk+1non è generalmente nota ma, grazie al lemma (3.1.1)

possiamo riscrivere la (3.4) in termini della legge condizionata di Vtk+1|Vtk .

Infatti, ssato k ≥ 1, la distribuzione di Vtk+1si può riscrivere come:

P (Vtk+1∈ dνk+1) = dνk+1

∫R

φmk(νk),σk(νk)(νk+1)P (Vtk ∈ dνk) = dνk+1E[φmk(Vtk ),σk(Vtk )(νk+1)

],

(3.5)dove φm,σ rappresenta la funzione di densità associata alla distribuzione normale N(m,σ2).

Inserendo (3.5) in (3.4) otteniamo:

Dk+1(ΓNk+1) =

Nk+1∑i=1

∫Ci(Γ

Nk+1 )

(νk+1 − xNk+1

i )2E[φmk(Vtk ),σk(Vtk )(νk+1)

]dνk+1. (3.6)

3.1. IL METODO 25

Sostituendo inne V con V otteniamo la seguente versione discreta della distorsione:

Dk+1(ΓNk+1) =

Nk+1∑i=1

∫Ci(Γ

Nk+1 )

(νk+1 − xNk+1

i )2E[φmk(Vtk ),σk(Vtk )(νk+1)

]dνk+1

=

Nk+1∑i=1

∫Ci(Γ

Nk+1 )

(νk+1 − xNk+1

i )2Nk∑j=1

φmk(x

Nki ),σk(x

Nki )

(νk+1)P (Vtk ∈ Cj(ΓNk))dνk+1, (3.7)

dove ΓNk = xNk1 , xNk2 , ..., xNkNk è l'Nk-quantizzatore al tempo tk relativo a Vtk .

Abbiamo quindi ottenuto una versione discreta e ricorsiva della distorsione, il che permette di calcolareil gradiente e la matrice Hessiana e quindi di applicare il metodo di Newton Raphson.


3.2 Applicazione al modello di Black and Scholes

Consideriamo la dinamica del modello di Black and Scholes:

dSt = rStdt+ σStdWt (3.8)

Lo schema di Eulero associato è:

Sk+1 = Sk + rSk∆k + σSk√

∆kZk, (3.9)

dove (Zk)k sono variabili aleatorie con distribuzione normale standard, indipendenti.

Possiamo calcolare il gradiente della distorsione come:

∂Dk+1

∂yk+1j

= 2

N∑i=1

(yk+1j −mk(yki ))

[Φ0,1(yk+1

j+ (yki ))− Φ0,1(yk+1j− (yki ))

]P (Ytk = yki ) (3.10)

+2

N∑i=1

σk(yki )[φ0,1(yk+1

j+ (yki ))− φ0,1(yk+1j− (yki ))

]P (Ytk = yki ),

dove

yk+1j+ (yki ) =

yk+1j + yk+1

j+1 − 2mk(yki )

2σk(yki )e yk+1

j− (yki ) =yk+1j−1 + yk+1

j − 2mk(yki )

2σk(yki )

Per l'Hessiana vale invece:

∂2Dk+1

∂2yk+1j

(yk+1) = 2

N∑i=1

[Φ0,1(yk+1

j+ (yki ))− Φ0,1(yk+1j− (yki ))

]P (Ytk = yki )

−N∑i=1

[1

2σk(yki )φ0,1(yk+1

j+ (yki ))(yk+1j+1 − y

k+1j )

]P (Ytk = yki )

−N∑i=1

[1

2σk(yki )φ0,1(yk+1

j− (yki ))(yk+1j − yk+1

j−1 )

]P (Ytk = yki ), (3.11)

∂2Dk+1

∂yk+1j ∂yk+1

j−1

(yk+1) = −1

2

N∑i=1

[φ0,1(yk+1

j− (yki ))

(yk+1j − yk+1

j−1

σk(yki )

)]P (Ytk = yki ),

∂2Dk+1

∂yk+1j ∂yk+1

j+1

(yk+1) = −1

2

N∑i=1

[φ0,1(yk+1

j+ (yki ))

(yk+1j+1 − y

k+1j

σk(yki )

)]P (Ytk = yki ).

Applicando il lemma (3.1.1) al modello di Black and Scholes otteniamo:

(Stk |Stk−1= x) ∼ N(mk−1(x), σ2k+1

(x)), (3.12)

con:mk−1(x) = x+ rx∆k,

σk−1(x) = σx√

∆k.

In ultimo, i pesi sono dati dalla seguente approssimazione:

P (Ytk+1∈ Ci(yk+1)) =

N∑j=1

[Φ0,1(yk+1

i+ (ykj ))− Φ0,1(yk+1i− (ykj ))

]P (Ytk = ykj ). (3.13)

3.2. APPLICAZIONE AL MODELLO DI BLACK AND SCHOLES 27

3.2.1 Implementazione numerica

Scriviamo in primo luogo una funzione che implementi il metodo di Newton Raphson.Tale funzione è così strutturata:

FUNZIONE: NEWTHON - RAPHSON

INPUT OUTPUT

parametro descrizione parametro descrizione

k n. di iterazione prob P aggiornatovect S al passo corrente vect S aggiornato

vectpast S al passo precedenteprobpast P al passo precedente

x0 spot pricer tasso d'interesse privo di rischioσ volatilitàdt passo temporaletoll tolleranza del metodo di N-R

PROCEDURA:

1. calcolo di media e varianza, tramite funzioni apposite che implementano le formule del lemma(3.12);

2. ciclo, no ad ottenere l'accuratezza desiderata, in cui dopo aver calcolato gradiente ed Hessiana (confunzioni apposite che implementino, rispettivamente, le formule (3.10) e (3.11) ) viene applicata laformula (1.22) per l'aggiornamento di vect;

3. aggiornamento del vettore dei pesi secondo la formula (3.13).


L'algoritmo implementato che applica la Quantizzazione Marginale Ricorsiva al il calcolo del prezzo diuna call europea, con sottostante che evolve secondo il modello di Black and Scholes è, in sintesi, cosìstrutturato:

ALGORITMO: RMQ - Call Europea - Black and Scholes

INPUT OUTPUT


K strike price (vettore) prezzo prezzo CallS0 spot pricer tasso d'interesse privo di rischioσ volatilitàT maturityN dimensioni quantizzatoreM n. di time stepstoll tolleranza del metodo di N-R

PROCEDURA:

1. inizializzazione di prob e vect (che saranno gli input per la funzione sopra descritta);

2. ciclo temporale: iterazione dell'algoritmo di Newton Raphson no ad ottenere il valore di vect (S)a scadenza;

3. calcolo del prezzo della Call per ogni valore di K.

Fissiamo i seguenti parametri:

Parametro Valore

S0 100Per il modello: r 0.04

σ 0.20

Per l'opzione: T 1

N 80Per le scelte numeriche 1: M 16

toll 10−2

3.2. APPLICAZIONE AL MODELLO DI BLACK AND SCHOLES 29

Di seguito un confronto dei prezzi ottenuti con questa procedura con quelli esatti calcolati tramite lafunzione di Matlab blsprice, per diversi valori dello strike-price K:

Strike Prezzo esatto Prezzo RMQ Errore relativo (%)

90 16.0604 16.0679 0.046392.5 14.3516 14.3544 0.019395 12.7575 12.7552 0.0174

97.5 11.2815 11.2757 0.0508100 9.9251 9.9136 0.1154

Il tempo di calcolo complessivo, per i 5 valori di K, è di 0.7819 secondi.

Analisi parametri numerici

Fissando ora K = 100 (opzione at the money), fermo restando i parametri di cui sopra, vediamo comecambiano l'errore e il costo computazionale al variare di N , toll ed M :

N Errore relativo (%) Tempo di calcolo

40 0.4490 0.247950 0.2888 0.407060 0.1981 0.543870 0.1480 0.631180 0.1154 0.760890 0.0976 1.0130100 0.0864 1.0527

Troviamo conferma della proposizione 1.2.2, ovvero della convergenza dell'errore all'aumentare di N . Dicontro, osserviamo che il costo computazionale aumenta considerevolmente all'aumentare delle dimen-sioni del quantizzatore.

toll Errore relativo (%) Tempo di calcolo

10−2 0.1154 0.716010−3 0.1220 0.854110−4 0.1223 0.9680

Osserviamo che al diminuire della tolleranza richiesta il tempo di calcolo aumenta (a causa del maggiornumero di interazioni all'interno del metodo di Newthon Raphson) e che, inoltre, non si ha un vantaggioin termini di errore chiedendo una tolleranza inferiore a 10−2.


M Errore relativo (%) Tempo di calcolo

4 0.2290 0.28956 0.1525 0.35978 0.1445 0.435110 0.1269 0.614812 0.1178 0.675914 0.1161 0.649816 0.1154 0.805518 0.1211 0.919520 0.1285 0.6965

Dalla teoria sappiamo che aumentare M porta, da un lato, a un minor errore dell'approssimazione delprocesso tramite schema di Eulero, dall'altro, a un maggior errore di quantizzazione, che viene propagatoad ogni passo temporale.Osserviamo che il punto di minimo per l'errore, in funzione di M , si trova per M = 16. Notiamo inoltreche, nonostante l'aumento del numero di passi temporali dovrebbe comportare un aumento del costocomputazionale, il metodo non risulta particolarmente sensibile a questo parametro; non si ha infattiuna variazione consistente del tempo di calcolo in funzione di M , se non per valori inferiori a M = 10.

Tutti i calcoli sono stati eettuati utilizzando Matlab su un computer con CPU 2.4 GHz e 8 GB di RAM.

3.3. APPLICAZIONE AL MODELLO DI HESTON 31

3.3 Applicazione al modello di Heston

In primo luogo riprendiamo la dinamica del modello di Heston vista nel capitolo (2.6), eq. (2.19) e (2.20):

dX(t) = X(t)rdt+X(t)√Vt(√

1− ρ2dWXt + ρdWV

t ),

X(0) = X0 > 0;

dV (t) = k(θ − V (t))dt+ ε√V (t)dWV

t ,

V (0) = V0 > 0

e lo schema di Eulero associato, eq. (2.22) e (2.23):

log(Xi+1) = log(Xi) + (r − 1

2V (t))∆t+

√Vi(√

1− ρ2ZXi + ρZYi )√

∆t,

V (i+ 1) = Vi + k(θ − Vi)∆t+ ε√ViZ

Vi

√∆t,

dove (Zi)i sono variabili aleatorie con distribuzione normale standard, indipendenti.

Per questo modello abbiamo, ad ogni passo ti, due vettori da quantizzare: Xti e Vti .Una strada percorribile sarebbe quello di estendere la teoria vista al caso d = 2 ma questo aumenterebbenotevolmente la complessità del problema, cerchiamo quindi una strada alternativa.

In questo contesto facciamo riferimento, per quanto esporremo, a [3] e [4].

Introduciamo il seguente lemma:

Lemma 3.3.1. Per ogni i = 0, ...,M − 1, condizionatamente a Sti = si, Vti+1= vi+1 e Vti = vi, la

variabile aleatoria Sti+1= vi+1 ha distribuzione gaussiana. In particolare vale:

(Sti+1|Sti , Vti+1

, Vti) ∼ N(mi(si, vi+1, vi), σ2i (si, vi)), (3.14)

dove:

mi(si, vi+1, vi) = si + rsi∆t+ ρ√visi

(vi+1 − vi − k(θ − Vi)∆t

ε√Vi

); (3.15)

σ2k(sk, vk) = (1− ρ2)vis

2i∆t. (3.16)

Sfruttando quest'ultimo risultato, possiamo dunque procedere come segue:

- Quantizzazione della variabile Vtk , data Vtk , con la procedura esposta per la quantizzazionemonodimensionale;

- Quantizzazione della variabile Stk , dati Stk−1, Vtk e Vtk−1

usando la conoscenza delle distribuzione,procedendo quindi anche qui con la quantizzazione monodimensionale.


3.3.1 Implementazione numerica

L'algoritmo implementato che applica la Quantizzazione Marginale Ricorsiva al il calcolo del prezzo diuna call europea, con sottostante che evolve secondo il modello di Heston è, in sintesi, così strutturato:

ALGORITMO: RMQ - Call Europea - Modello di Heston

INPUT OUTPUT


K strike price (vettore) prezzo prezzo CallS0 spot pricer tasso d'interesse privo di rischio

param ξ, k, ρ, θ, V0

T maturityNS dimensione quantizzatore per SNV dimensione quantizzatore per VM n. di time stepstoll tolleranza del metodo di N-R

PROCEDURA:

1. inizializzazione vettori;

2. quantizzazione processo V tramite RMQ, con la stessa procedura vista nella sezione 3.2.1 per laquantizzazione di S;

3. ciclo temporale: iterazione dell'algoritmo di Newton Raphson no ad ottenere il valore di S ascadenza;


ALGORITMO: RMQ - processo di volatilità - Heston

INPUT OUTPUT


param ξ, k, ρ, θ, V0 V V quantizzataT maturityN dimensioni quantizzatoreM n. di time stepstoll tolleranza del metodo di N-R

3.3. APPLICAZIONE AL MODELLO DI HESTON 33

Nella nostra implementazione abbiamo scelto un modello con i seguenti parametri per il sottostante eper il processo di volatilità:

Parametro Valore

S0 100V0 0.0456r 0.0400

Per il modello: ε 0.5087ρ -0.7495k 1.7029θ 0.0600

Per l'opzione: T 1

NS 20Per il metodo numerico 2: NV 10

M 16

Di seguito, un confronto dei prezzi ottenuti con questa procedura con quelli calcolati tramite il metodobasato sulla FFT, presi come Benchmark:

Strike Benchmark Prezzo RMQ Errore relativo (%)

90 17.3018 17.2267 0.434392.5 15.5079 15.4732 0.224195 13.7912 13.7444 0.3392

97.5 12.1584 12.1521 0.0513100 10.6161 10.6210 0.0463

Il tempo di calcolo complessivo. per i 5 valori di K, è di 12.9885 secondi.

I calcoli sono stati eettuati utilizzando Matlab su un computer con CPU 2.4 GHz e 8 GB di RAM.

Capitolo 4

Quantizzazione e trasformata di

Fourier

In questo capitolo introduciamo una nuova metodologia di pricing, che sfrutta sia la quantizzazione chela trasformata di Fourier.Per la teoria riguardante questa tecnica, si farà riferimento principalmente a [1].

Introdurremo i concetti teorici fondamentali (sezione 4.1), per poi passare alle applicazioni nei modellidi Black and Scholes (sezione 4.2.1), Kou (sezione 4.2.2), Variance Gamma (sezione 4.2.3) e di Heston(sezione 4.2.4).

4.1 Il metodo

L'idea alla base di questo metodo consiste nello sfruttare la tecnica di Fourier per il calcolo delle proba-bilità che il processo appartenga a una deteerminata cella di Voronoi.

Consideriamo lo spazio di probabilità (Ω, (Ft)t∈[0,T ], µ), dove (Ft)t∈[0,T ] è una ltrazione 1.Consideriamo un processo stocastico (St)t∈[0,T ].

Deniamo la funzione caratteristica condizionata di log(St) (log-prezzo):

Φt,T (u) = E[eiulog(S0)|Ft]. (4.1)

È noto dalla letteratura (si veda ad esempio [19]) che la densità e la distribuzione del processo (St)t∈[0,T ]

possono essere ricavate dalla funzione caratteristica, tramite la trasformata di Fourier.

In particolare vale:

P (ST ∈ dz) =

(1

π

1

z

∫ +∞

0

Re(e−iulog(z)Φ0,T (u))du

)dz =: f(z)dz (4.2)

e

P (ST ≤ dz) =1

2− 1

π

∫ +∞

0

Re

((e−iulog(z)Φ0,T (u)

iu

)du. (4.3)

1Si denisce ltrazione su (Ω, F, µ) una famiglia σ−algebre (Ft)t denite sullo spazio (Ω, F ) tali che Ft ⊂ F per ogni te Fs < Ft per ogni s e t tali che s < t.

35

36 CAPITOLO 4. QUANTIZZAZIONE E TRASFORMATA DI FOURIER

In questo lavoro sfrutteremo proprio queste espressioni per esprimere la densità di ST presente nellafunzione di distorsione (1.19).

Essenzialmente procederemo in questo modo:

- espressione della distorsione tramite la rappresentazione di Fourier della densità dei prezzi;

- dierenziazione della funzione di distorsione e calcolo dell'Hessiana;

- delineamento dell'equazione da risolvere tramite il metodo di Newthon Raphson.

4.1.1 Quantizzatori ottimi

Consideriamo il caso d = 1, in cui le celle di Voronoi sono intervalli di R+: Ci(Γ) := [x−i , x+i ],

dove: x−1 = 0 , x−i =xi−1 + xi

2∀i ∈ 2, ..., N; x+

i =xi + xi+1

2∀i ∈ 1, ..., N − 1 , x+

N = 1.

La distorsione ha la seguente forma:

Dp(Γ) =

N∑i=1

∫ x+i

x−i

|ξ − xi|2dPST (z), (4.4)

dove con PST (z) abbiamo indicato la distribuzione di ST .

Considerando che la densità di ST è continua e concentrata nell'intervallo (0,+∞), dovremo discretizzarela funzione di distorsione all'interno dell'integrale della (4.4).

Introduciamo il seguente teorema per l'identicazione dei quantizzatori stazionari:

Teorema 4.1.1. (Quantizzatori stazionari) ([1])Sia ΓN = x1, ..., xN una grigia di quantizzazione di dimensione N e sia dato p ∈ (1,+∞). ΓN èp-stazionaria per (4.4) se per ogni i = 1, ..., N vale:∫ +∞

0

Re

[e−iulog(xi)Φ0,T (u)

(β

(x−ixi,−iu, p

)− β

(xi

x+i

, 1− p+ iu, p

))]du = 0, (4.5)

dove la funzione β = β(x, a, b) è denita per ogni a ∈ C, bt.c.Re(b) > 0, x ∈ (0, 1), nel seguente modo:

β(x, a, b) :=

∫ 1

x

ta−1(1− t)b−1dt. (4.6)


Abbiamo ricondotto il nostro problema alla risoluzione di un sistema non lineare di equazioni (4.5).Anche in questo caso, come con la Quantizzazione Marginale Ricorsiva, procederemo dunque all'imple-mentazione di un algoritmo di Newton Raphson, come proposto in [5], [4] e [2].

4.1. IL METODO 37

4.1.2 Newton Raphson

Denotiamo con Li(Γ) la derivata della funzione di distorsione (4.4) rispetto a xi, ovvero:

Li(ΓN ) =

∂Dp(ΓN )

∂xi.

Vale:

Li(ΓN ) :=

p

πxp−1i

∫ +∞

0

Re

[e−iulog(xi)Φ0,T (u)

(β

(x−ixi,−iu, p

)− β

(xi

x+i

, 1− p+ iu, p

))]du. (4.7)

Trovare i quantizzatori stazionari consiste dunque nel risolvere il sistema L(ΓN ) = 0, che è equivalentealla (4.5).

Per applicare il metodo di Newton Raphson occorre calcolare ∇(L), ovvero la matrice Jacobiana 2 NxNdella funzione L.Osserviamo che Li dipende unicamente da xi−1, xi e xi+1, per cui ∇(L) risulta essere tridiagonale.

Teorema 4.1.2. (algoritmo di quantizzazione tramite trasformata di Fourier) ([1])Sia Γ0 una griglia di quantizzazione di dimensione N , sia p ∈ (i,+∞). L'algoritmo di Newthon Raphsonper la determinazione dei quantizzatori p-stazionari ha la seguente forma ricorsiva:

Γk+1 = Γk − [∇L(Γk)]−1L(Γk), (4.8)

dove L è denito dalla (4.7) e le componenti di ∇L sono così denite:

∇Li,i+1(Γ) = − p

2πx+i

(xi+1−xi

2

)p−1 ∫ +∞0

Re[e−iulog(x

+i )]du ;

∇Li,i−1(Γ) = − p

2πx−i

(xi−xi−1

2

)p−1 ∫ +∞0

Re[e−iulog(x

−i )]du ;

∇Li,i(Γ) = −p(p−1)π xp−2

i

∫ +∞0

Re[e−iulog(xi)

(β(x−i

xi,−iu, p− 1

)− β

(xix+i

, 2− p+ iu, p− 1))]

du

+∇Li,i+1(Γ) +∇Li,i−1(Γ).


Osservazione 12. Nel caso p = 2, ∇L risulta:

∇Li,i+1(Γ) = −(xi+1 − xi

2πx+i

)∫ +∞

0

Re[e−iulog(x

+i )]du;

∇Li,i−1(Γ) = −(xi − xi−1

2πx−i

)∫ +∞

0

Re[e−iulog(x

−i )]du; (4.9)

∇Li,i(Γ) = − 2

π

∫ +∞

0

Re

[e−iulog(xi)

(β

(x−ixi,−iu, 1

)− β

(xi

x+i

, iu, 1

))]du+∇Li,i+1(Γ)+∇Li,i−1(Γ).

2Data una funzione vettoriale f : RN −→ RM , la sua matrice di Jacobi (o matrice jacobiana) è la matrice i cui elementisono le derivate parziali prime della funzione. Il nome è dovuto a Carl Gustav Jacob Jacobi


4.1.3 Stima dell'errore

Come già osservato, i quantizzatori ottimi (quelli che minimizzano l'Lp-errore di quantizzazione) sonostazionari mentre i quantizzatori stazionari (i punti stazionari della funzione di distorsione) non sononecessariamente ottimi. È quindi importante fare un'analisi dell'errore di quantizzazione.

Ricordiamo innanzitutto la denizione dell'errore (1.5):

ep,N (X,Γ) := ‖X − X‖Lp = min1≤i≤N

‖X − xi‖Lp =

(∫Rd

min1≤i≤N

|ξ − xi|pµ(dξ)

)1

p.

Presentiamo in questa sezione due dierenti risultati riguardo il comportamento asintotico di ep,N (ST ,Γ).

Lemma 4.1.3. (Lemma di Pierce, 1970)Sia E[|ST |ρ+δ] < +∞ per qualche δ > 0 , allora l'Lp- errore di quantizzazione soddisfa la seguenteespressione:

Npep,N (ST ,Γ) ≤(C1E

(|ST |ρ+δ

)+ C2

), (4.10)

dove:

- N ≥ C3;

- le costanti C1, C2eC3 dipendono da p e da δ ma non da N o dalla distribuzione di ST .

Questo risultato fornisce un criterio di convergenza dell'errore dove le costanti sono legate alla regolaritàdel sottostante. In alcuni contesti questo può essere un problema e le costanti Ci come limite dell'errorediventano quindi inutili. Introduciamo allora un altro risultato in cui le costanti non dipendono dallaregolarità.

Teorema 4.1.4. ([1])Assumiamo che ST abbia p-esimo momento nito. Assumiamo che la sua densità (come da eq. (4.2))abbia un'espressione polinomiale in 0 e in +∞. Allora ep,N (ST ,Γ) ha il seguente comportamentoasintotico:

limN→+∞

Npep,N (ST ,Γ) ≤‖f‖ 1

p+1

2p(p+ 1). (4.11)


Osservazione 13. L'ipotesi di questo teorema è leggermente più forte di quella del lemma di Pierce.Infatti, chiedere che il p-esimo momento sia nito e che, allo stesso tempo, il comportamento della densitàsia nito in +∞, equivale chiedere l'esistenza del (p + δ)-esimo momento nito, per δ sucientementepiccolo.

Osservazione 14. Quest'ultima stima dell'errore non dipende più da δ ma in compenso, dipende dalladistribuzione del processo al tempo T (4.2).

4.1.4 La scelta di p

Per analizzare come variano i quantizzatori al variare di p, consideriamo la distribuzione di una variabilealeatoria normale standard N(0, 1), e un valore di N ssato pari a 10.

Nella Figura 4.1 possiamo osservare come, al crescere di p, i punti della griglia ottima tendano ad esseresempre più distanziati rispetto all'origine 3. Segue che la scelta di p è legata all'interesse che si ha nellostudio delle code della distribuzione.

Il valore più comunemente utilizzato di p è comunque, anche in questo caso, p = 2.In questo lavoro ci atterremo alla quantizzazione quadratica.

3Fonte: [1].

4.1. IL METODO 39

Figura 4.1: Quantizzatori p-ottimi, al variare di p, per Z ∼ N(0, 1) con N = 10 ssato.


4.2 Implementazione numerica

Per quanto riguarda l'implementazione numerica, questa si dierenzia da modello a modello solo nelladenizione della funzione caratteristica ad esso associata.

In questa sezione andiamo quindi ad illustrare la procedura generale per poi vedere, per ogni modello, irisultati ottenuti.

L'algoritmo che utilizzeremo è così strutturato:

ALGORITMO: QUANTIZZAZIONE + FOURIER -

INPUT OUTPUT


K strike price (vettore) prezzo prezzo CallS0 spot pricer tasso d'interesse privo di rischio

param ξ, k, ρ, θ, V0

T maturityN numero di punti della griglia

max_int limite per il calcolo degli integralitoll tolleranza del metodo di N-R

PROCEDURA:

1. inizializzazione griglie;

2. iterazione dell'algoritmo di Newton Raphson per ottenere il valore di S a scadenza con l'accuratezzadesiderata; in questa fase vengono calcolate Jacobiana ed Hessiana tramite apposita funzione cheimplementa le formule denite dalla 4.9;

3. calcolo dei pesi secondo la formula 4.2, tramite apposita funzione;


4.2. IMPLEMENTAZIONE NUMERICA 41

4.2.1 Applicazione al modello di Black and Scholes

Per questa applicazione si è scelto di utilizzare, per il modello, gli stessi parametri della sezione 3.2, inmodo da poter fare un confronto dei risultati ottenuti.

Riportiamo qui i parametri utilizzati:

Parametro Valore

S0 100Per il modello: r 0.04

σ 0.20

Per l'opzione: T 1

N 150Per le scelte numeriche: max_int 20

toll 10−5

Di seguito un confronto dei prezzi ottenuti con questa procedura con quelli esatti calcolati tramite lafunzione di Matlab blsprice:

Strike Prezzo esatto Prezzo quant + Fourier Errore relativo (%)

90 16.0604 16.0596 0.005192.5 14.3516 14.3510 0.004495 12.7575 12.7572 0.0022

97.5 11.2815 11.2814 0.0004100 9.9251 9.9250 0.0001

Il tempo di calcolo complessivo per i 5 valori di K, è di 1.2408 secondi.


Analisi parametri numerici

Fissando ora K = 100 (opzione at the money), fermo restando i parametri di cui sopra, vediamo comecambiano l'errore e il costo computazionale al variare di N , toll ed max_int:

N Errore relativo (%) Tempo di calcolo

110 0.0143 1.7488120 0.0063 1.7410130 0.0022 1.7138140 0.0004 1.7430150 0.0001 2.0134160 0.0001 2.2311

Osserviamo, come ci aspettavamo, che all'aumentare delle dimensioni del quantizzatore si ha una mag-giore accuratezza e un aumento del costo computazionale. Il tempo di calcolo rimane molto basso e convariazioni comunque poco rilevanti. La scelta ottimale risulta essere N = 150.

max_int Errore relativo (%) Tempo di calcolo

10 0.4825 0.695715 0.3791 0.817320 0.0001 1.137225 0.0001 1.481630 0.0001 1.5721

Anche qui si ha un riscontro coerente alla teoria, ovvero un aumento della precisione nel calcolo degliintegrali comporta un maggior costo computazionale ed una maggior precisione. Non si evidenzianovantaggi aumentando il valore oltre max_int = 20.

toll Errore relativo (%) Tempo di calcolo

10−2 0.00039 0.869910−3 0.00019 1.066710−4 0.00010 1.175910−5 0.00005 1.213310−6 0.00005 1.3798

Si evidenzia come anche solo una tolleranza di 10−2 comporti un errore molto contenuto, ad ogni modo,dato il non eccessivo aumento del costo computazionale, si valuta migliore utilizzare una tolleranza di10−5. Si evidenzia come diminuire ulteriormente la tolleranza non comporta vantaggi in termine di errore.

I calcoli sono stati eettuati utilizzando Matlab su un computer con CPU 2.4 GHz e 8 GB di RAM.


4.2.2 Applicazione al modello Kou

Nella nostra implementazione abbiamo i seguenti parametri:

Parametro Valore

S0 100r 0.01σ 0.2

Per il modello: λ 0.5λ+ 50λ− 10p 0.6

Per l'opzione: T 1


toll 10−5


Strike Benchmark Prezzo quant + Fourier Errore relativo (%)

90 14.5008 14.5006 0.001192.5 12.8810 12.8810 0.000395 11.3821 11.3823 0.0009

97.5 10.0055 10.0058 0.0032100 8.7495 8.7509 0.0159

Il tempo di calcolo complessivo per i 5 valori di K, è di 1.3715 secondi.I calcoli sono stati eettuati utilizzando Matlab su un computer con CPU 2.4 GHz e 8 GB di RAM.


4.2.3 Applicazione al modello Variance Gamma

Nella nostra implementazione abbiamo i seguenti parametri:

Parametro Valore

S0 100r 0.01

Per il modello: σV G 0.1243kV G 0.1686θV G 0.0436

Per l'opzione: T 1


toll 10−4



90 12.6872 12.6863 0.007392.5 10.7155 10.7154 0.000895 8.9155 8.9156 0.0011

97.5 7.3099 7.3102 0.0045100 5.9110 5.9119 0.0149



4.2.4 Applicazione al modello di Heston

Per questa applicazione si è scelto di utilizzare, per il modello, gli stessi parametri della sezione 3.3, inmodo da poter fare un confronto dei risultati ottenuti.Riportiamo qui i parametri utilizzati:

Parametro Valore

S0 100V0 0.0456r 0.0400

Per il modello: ε 0.5087ρ -0.7495k 1.7029θ 0.0600

Per l'opzione: T 1


toll 10−4



90 17.3018 17.2996 0.0129592.5 15.5079 15.5054 0.016895 13.7912 13.7887 0.0179

97.5 12.1584 12.1560 0.0190100 10.6161 10.6150 0.0103


Capitolo 5

Conclusioni

L'obiettivo del lavoro di tesi era analizzare le applicazioni della quantizzazione nel contesto del pricingdi derivati, per diverse dinamiche del sottostante, nel caso di opzioni plain vanilla (europee).

Le due metodologie esposte presentano delle dierenze sostanziali in termini di applicabilità.In particolare:

- La Quantizzazione Marginale Ricorsiva si basa su uno schema di Eulero per il sottostante, ed èdunque applicabile ai modelli per cui tale schema è disponibile. Segue che non risulta funzionaleai processi di salto.

- La Quantizzazione con trasformata di Forier sfrutta la funzione caratteristica, ed è dunque appli-cabile ai modelli ove essa sia nota in forma analitica.

Nell'insieme si è mostrato che la quantizzazione è eettivamente applicabile con successo al pricing diopzioni nanziarie, per tutte le diverse dinamiche del sottostante prese in considerazione.

Abbiamo riscontrato nello specico che, nell'ambito delle opzioni europee:

- Il metodo della Quantizzazione con trasformata di Fourier risulta accurato e non dispendioso intermini di costo computazionale, garantendo un errore di un ordine di grandezza di 10−3 (quasisempre, 10−2 nei casi peggiori).

- Il metodo della Quantizzazione Marginale ricorsiva risulta poco competitivo sotto entrambi i puntidi vista. Infatti, oltre a richiedere tempistiche più elevate, i risultati ottenuti sono poco precisi giàsulle seconde cifre decimali, talvolta anche sulle prime, con un ordine di grandezza dell'errore di10−1 (nei casi migliori 10−2).

Possibili approfondimenti

In questo lavoro, per mantenere il focus sulle diverse tipologie di sottostante, ci si è limitati al caso delleopzioni europee.Una possibile estensione dello studio può essere fatta analizzando le diverse tipologie di derivati, perentrambe le metodologie esposte.

Opzioni americane

Entrambi i metodi illustrati nel lavoro di tesi sono applicabili alle opzioni americane.Per la teoria e alcuni esempi di applicazione nel caso di Quantizzazione con trasformata di Fourier, sirimanda a [1], che mostra come sia possibile sfruttare la conoscenza delle probabilità di transizione per

47

48 CAPITOLO 5. CONCLUSIONI

una programmazione dinamica di tipo backward.Nel caso della quantizzazione marginale ricorsiva si rimanda invece a [8] e [15].

Opzioni path dependent

La quantizzazione marginale ricorsiva è applicabile alle opzioni esotiche. Per un primo studio in questosenso, relativo in particolare alle opzioni barriera, si rimanda a [9].Non è possibile invece applicare la quantizzazione con trasformata di Fourier in questo contesto, inquanto il valore del sottostante viene calcolato solo a scadenza.

Calibrazione

Un ulteriore campo di interesse può essere quello della calibrazione dei parametri dei modelli sui dati dimercato.La metodologia che sfrutta contemporaneamente quantizzazione e trasformata di Fourier risulta ecien-te ed accurata, dunque adatta ad essere utilizzata per la calibrazione. Si rimanda a [1] per una primaanalisi in questo senso.

Si conclude, in sintesi, che ci sono diverse possibili applicazioni della quantizzazione nel campo del pricingdi derivati e che, oltre ai risultati positivi già riscontrati, vi è ancora margine per un miglioramento.

Appendice A

Codici Quantizzazione Marginale

Ricorsiva

Modello di Black and Scholes

PRICING CALL EUROPEA

1 % datiS0 =100; T=1; r=0.04; sigma =0.2;

3 K=90:2.5:100;

5 % prezzo esattoexactprice=blsprice(S0 ,K,r,T,sigma);

7% prezzo RMQ

9 N=80; M=16; toll=1e-2;tic

11 call_marginal=bs_rmq(S0 ,K,r,T,sigma ,N,M,toll);toc

13disp('RMQ')

15 [exactprice ' call_marginal (exactprice '-call_marginal)./( exactprice ')]

BS RMQ

1 function call=bs_rmq(x0 ,K,r,T,sig ,N,M, toll)

3 %% inizializzazionedelta_t=T/M;

5 vect=x0 -(N/20 -0.1) :0.1: x0+N/20;x=( repmat(vect ,M+2,1))';

7 p=zeros(N,M+1);p(N/2,1)=1;

9%% NEWTHON RAPHSON

11 for i=1:M[p(:,i+1),x(:,i+1)]= nr_fast(i,x(:,i+1),x(:,i),p(:,i),x0 ,r,sig ,

delta_t ,toll);13 x(:,i+2)=x(:,i+1);

end

49

50 APPENDICE A. CODICI QUANTIZZAZIONE MARGINALE RICORSIVA

15%% pricing call

17 size=length(K);call=zeros(size ,1);

19 for i=1: sizepayoff=max(x(:,M+1)-K(i) ,0) '*p(:,M+1);

21 call(i,1) = exp(-r*T).*mean(payoff);end

23 end

FUNZIONI AUSILIARIE

1 function y=vcapplus(~,vect ,med ,sigm)N=length(vect);

3 vect(N+1)=inf;V=repmat ((vect (1:N)+vect (2:N+1))',N,1);

5 M=repmat (2*med ,1,N);S=repmat (2*sigm ,1,N);

7 y=(V-M)./S;end

9function y=vcapminus(~,vect ,med ,sigm)

11 N=length(vect);vect_new=zeros(N+1,1);

13 vect_new (1,1)=-inf;vect_new (2:N+1,1)=vect;

15 S=repmat (2*sigm ,1,N);V=repmat (( vect_new (1:N)+vect_new (2:N+1))',N,1);

17 M=repmat (2*med ,1,N);y=(V-M)./S;

19 end

1 % calcolo di media e varianzafunction y=media(~, x,r,deltat)

3 y=x.*(1+r*deltat);end

5function y=sigma(~,x,x0 ,sig ,deltat)

7 y=sig.*x.*sqrt(deltat);end

%% calcolo del gradiente2 function y_vect=g_dist(vect ,cxp ,cxm ,pxp ,pxm ,sigm ,med ,prob_past)

N=length(vect);4 V=repmat(vect ',N,1);

M=repmat(med ,1,N);6 S=repmat(sigm ,1,N);

NCDF=(cxp -cxm);8 NPDF=(pxp -pxm);

y_vect =(2*((V-M).*NCDF+S.*NPDF))'*prob_past;10 end

51

%% calcolo della 'hessiana2 function M=h_dist(vect ,xp ,xm ,cxp ,cxm ,pxp ,pxm ,sigm ,prob_past)

N=length(vect);4 %calcolo della diagonale principale

d=zeros(N,1);6 xjp1=xp(:,1);

xjm1=xm(:,1);8 xjpN=xp(:,N);

xjmN=xm(:,N);10 d(1,1) =(2.*(( normcdf(xjp1)-normcdf(xjm1))...

-1./(4.* sigm).*( normpdf(xjp1)*(vect (2)-vect (1))))) '*prob_past;12 d(N,1) =(2.*(( normcdf(xjpN)-normcdf(xjmN))...

-1./(4.* sigm).*( normpdf(xjmN)*(vect(N)-vect(N-1))))) '*prob_past;14 V1=repmat ((vect (3:N)-vect (2:N-1))',N,1);

V2=repmat ((vect (2:N-1)-vect (1:N-2))',N,1);16 S=repmat (1./(4* sigm),1,N-2);

d(2:N-1,1) =(2*(( cxp(:,2:N-1)-cxm(:,2:N-1))...18 -V1.*S.*pxp(:,2:N-1)-V2.*S.*pxm(:,2:N-1)))'*prob_past;

%calcolo della superdiagonale20 V=repmat ((vect (2:N)-vect (1:N-1))',N,1);

S=repmat (1./sigm ,1,N-1);22 dsup =( -1/2*(V.*S.*pxp(:,1:N-1))) '*prob_past;

%calcolo dell 'hessiana24 M=diag(d)+diag(dsup ,1)+diag(dsup ,-1);

end

1 %% algoritmo di newton raphsonfunction [prob ,vect]= nr_fast(k,vect ,vect_past ,prob_past ,x0 ,r,sig ,deltat

,toll)3 N=length(vect);

med=media(k,vect_past ,r,deltat);5 sigm=sigma(k,vect_past ,x0,sig ,deltat);

temp=zeros(N,1);7 iter =0;

while norm(vect -temp)>toll9 iter=iter +1;

xp=vcapplus(k,vect ,med ,sigm);11 xm=vcapminus(k,vect ,med ,sigm);

cxp=normcdf(xp); cxm=normcdf(xm);13 pxp=normpdf(xp); pxm=normpdf(xm);

temp=vect;15 gdist=g_dist(vect ,cxp ,cxm ,pxp ,pxm ,sigm ,med ,prob_past);

A=h_dist(vect ,xp,xm ,cxp ,cxm ,pxp ,pxm ,sigm ,prob_past);17 if rcond(A)<1e-08

display('error ')19 A=pinv(A); vect=vect - A*gdist;

else21 [L,U]=lu(A); vect=vect - U\(L\(gdist));

end23 vect=sort(vect);

end25 prob=(cxp -cxm)'*prob_past;

end


Modello di Heston


%% DATI2 T=1; K=90:2.5:100;

S0 =100; v0 =0.0456; r = 0.04;4 csi =0.5087; rho = -0.7495; kappa =1.7029; theta =0.0600;

param =[csi kappa rho theta v0];6

%% RMQ8 tic

NS=20; NV=10; M=16; toll=1e-4;10 Call_RMQ=heston_rmq(K,S0 ,param ,r,T,NS,NV,M,toll);

time_RMQ=toc12

%% FFT14 tic

Call_fft=heston_fft(param ,S0 ,K,r,T);16 time_fft=toc;

18 disp(' BENCHMARK - RMQ')[Call_fft ', Call_RMQ , (Call_fft ' - Call_RMQ)./( Call_fft ')]

20 [time_fft , time_RMQ]

VOLATILITY RMQ

function V=heston_volatility_rmq(N,M,toll ,T,param)2 deltat=T/M;

x0=param (5);4 step =0.001;

vect=x0 -(N/2*step -step):step:x0+N/2* step;6 x=( repmat(vect ,M+1,1))';

p=zeros(N,M+1);8 p(N/2,1)=1;

for i=1:M-110 [p(:,i+1),x(:,i+1)]= nr_fast(i,x(:,i+1),x(:,i),p(:,i),param ,deltat ,

toll);x(:,i+2)=x(:,i+1);

12 end[p(:,M+1),x(:,M+1)]= nr_fast(M,x(:,M+1),x(:,M),p(:,M),param ,deltat ,toll)

;14 x=x(:,2:end);

p=p(:,2:end);16 V=x;

return18 end

53

FUNZIONI AUSILIARIE (VOLATILITY RMQ)

%% calcolo di vcap+ e vcap -2 function y=vcapplus(~,vect ,med ,sigm)

N=length(vect);4 vect(N+1)=inf;

V=repmat ((vect (1:N)+vect (2:N+1))',N,1);6 M=repmat (2*med ,1,N);

S=repmat (2*sigm ,1,N);8 y=(V-M)./S;

end10 function y=vcapminus(~,vect ,med ,sigm)

N=length(vect);12 vect_new=zeros(N+1,1);

vect_new (1,1)=-inf;14 vect_new (2:N+1,1)=vect;

S=repmat (2*sigm ,1,N);16 V=repmat (( vect_new (1:N)+vect_new (2:N+1))',N,1);

M=repmat (2*med ,1,N);18 y=(V-M)./S;

end

1 %% calcolo di media e varianzafunction y=media(~,x,deltat , kappa , theta)

3 y=x+kappa *(theta -x)*deltat;end

5 function y=sigma(~,x,deltat , csi)y=csi*sqrt(abs(x))*sqrt(deltat);

7 end

1 %calcolo del gradiente.function y_vect=g_dist(vect ,cxp ,cxm ,pxp ,pxm ,sigm ,med ,prob_past)

3 N=length(vect);V=repmat(vect ',N,1);

5 M=repmat(med ,1,N);S=repmat(sigm ,1,N);

7 NCDF=(cxp -cxm);NPDF=(pxp -pxm);

9 y_vect =(2*((V-M).*NCDF+S.*NPDF))'*prob_past;return

11 end

1 %%Calcolo dell 'Hessianafunction M=h_dist(vect ,xp ,xm ,cxp ,cxm ,pxp ,pxm ,sigm ,prob_past)

3 N=length(vect);

5 %diagonale principaled=zeros(N,1);

7 xjp1=xp(:,1);xjm1=xm(:,1);

9 xjpN=xp(:,N);xjmN=xm(:,N);

11 d(1,1) =(2.*(( normcdf(xjp1)-normcdf(xjm1))...-1./(4.* sigm).*( normpdf(xjp1)*(vect (2)-vect (1))))) '*prob_past;


13 d(N,1) =(2.*(( normcdf(xjpN)-normcdf(xjmN))...-1./(4.* sigm).*( normpdf(xjmN)*(vect(N)-vect(N-1))))) '*prob_past

;15 V1=repmat ((vect (3:N)-vect (2:N-1))',N,1);

V2=repmat ((vect (2:N-1)-vect (1:N-2))',N,1);17 S=repmat (1./(4* sigm),1,N-2);

d(2:N-1,1) =(2*(( cxp(:,2:N-1)-cxm(:,2:N-1))...19 -V1.*S.*pxp(:,2:N-1)-V2.*S.*pxm(:,2:N-1)))'*prob_past;

21 %superdiagonaleV=repmat ((vect (2:N)-vect (1:N-1))',N,1);

23 S=repmat (1./sigm ,1,N-1);dsup =( -1/2*(V.*S.*pxp(:,1:N-1))) '*prob_past;

25%hessiana

27 M=diag(d)+diag(dsup ,1)+diag(dsup ,-1);end

%Algoritmo di Newton Raphson2 function [prob ,vect]= nr_fast(k,vect ,vectopt_past ,prob_past ,param ,

deltat ,toll)csi=param (1);

4 kappa=param (2);rho=param (3);

6 theta=param (4);v0=param (5);

8 N=length(vect);med=media(k,vectopt_past ,deltat , kappa , theta);

10 sigm=sigma(k,vectopt_past ,deltat , csi);vect_cp=zeros(N,1);

12 iter =0;while norm(vect -vect_cp)>toll

14 iter=iter +1;xp=vcapplus(k,vect ,med ,sigm);

16 xm=vcapminus(k,vect ,med ,sigm);cxp=normcdf(xp); cxm=normcdf(xm);

18 pxp=normpdf(xp); pxm=normpdf(xm);vect_cp=vect;

20 gdist=g_dist(vect ,cxp ,cxm ,pxp ,pxm ,sigm ,med ,prob_past);A=h_dist(vect ,xp,xm ,cxp ,cxm ,pxp ,pxm ,sigm ,prob_past);

22 if rcond(A)<1e-08display('error ')

24 A=pinv(A); vect=vect - A*gdist;else

26 [L,U]=lu(A);vect=vect - U\(L\( gdist));

28 endvect=sort(vect);

30 endprob=(cxp -cxm)'*prob_past;

32 end

55

HESTON RMQ

function Hcall=heston_rmq(K,x0 ,param ,r,T,N,N2,M,toll)2 format long

4 %% inizializzazionedeltat=T/M;

6 step =1;vect_test=x0 -(N/2*step -step):step:x0+N/2* step;

8 x=zeros(N,M);x(:,1)=vect_test ';

10 p=zeros(N,M);pdouble=zeros(N,N2,M);

12 p_trans=zeros(N,N2,N,N2,M-1);

14 %% quantizzazione del processo di volatilitÃ v=heston_volatility_rmq(N2 ,M,toll ,T,param);

16%% NEWTHON RAPHSON

18 [x(:,1),p(:,1),pdouble (:,:,1)]= nr_fastinit (1,x(:,1),v(:,1),x0 ,r,param ,deltat ,1e-2);

x(:,2)=x(:,1);20 for i=2:M-1

[x(:,i),p(:,i),pdouble (:,:,i),p_trans (:,:,:,:,i-1)]=...22 nr_fast(i,x(:,i),v(:,i),x(:,i-1),v(:,i-1),pdouble(:,:,i-1),r,

param ,deltat ,1e-4);x(:,i+1)=x(:,i);

24 end[x(:,M),p(:,M),pdouble (:,:,M),p_trans (:,:,:,:,M-1)]=...

26 nr_fast(i,x(:,M),v(:,M),x(:,M-1),v(:,M-1),pdouble(:,:,M-1),r,param ,deltat ,1e-4);

Hcall=zeros(length(K) ,1);28

%% pricing call30 for i=1: length(K)

payoff=max(x(:,M)-K(i) ,0) '*p(:,M);32 Hcall(i,1) = exp(-r*T).*mean(payoff);

end34 end


FUNZIONI AUSILIARIE (HESTON RMQ)

%% calcolo dei pesi per S e V e delle probabilitÃ di transizione2 function [M,M_]= prob_double(s_now ,v_now ,mx,my,sx,sy ,prob_past ,rho)

snew=[-Inf; s_now; Inf];4 vnew=[-Inf; v_now; Inf];

sp =0.5*( snew (3:end)+snew (2:end -1));6 sm =0.5*( snew (2:end -1)+snew (1:end -2));

vp =0.5*( vnew (3:end)+vnew (2:end -1));8 vm =0.5*( vnew (2:end -1)+vnew (1:end -2));

N1=length(sp);10 N2=length(vp);

[AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sp ,vp ,mx ,my);12 BB1=AA1 -AA3;

BB2=AA2 -AA4;14 CC1=reshape ([ reshape(BB1 ,1,[]),reshape(BB2 ,1 ,[])],[],2);

S__1=reshape(kron(sx,ones(N1*N2 ,1)) ,[],1);16 S__2=reshape(kron(sy,ones(N1*N2*N1 ,1)) ,[],1);

CC2=CC1./[S__1 ,S__2];18 DD1=internal.stats.bvncdf(CC2 ,rho ,1e-8);

[AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sm ,vm ,mx ,my);20 BB1=AA1 -AA3;




[AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sp ,vm ,mx ,my);28 BB1=AA1 -AA3;




[AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sm ,vp ,mx ,my);36 BB1=AA1 -AA3;




DDtot=DD1+DD2 -DD3 -DD4;44 M_=reshape(DDtot ,N1,N2,N1,N2);

prob=reshape(prob_past ,1,1,N1,N2);46 M=sum(sum(bsxfun (@times ,M_ ,prob) ,4) ,3);

end

1 % calcolo del gradientefunction y=g_dist(vect ,Fdiff ,Gdiff ,vectprob)

3 y=sum(sum(bsxfun (@times ,Fdiff -bsxfun (@times ,vect ,Gdiff),vectprob) ,3) ,2);

end

57

%calcolo della 'hessiana2 function M=h_dist(vect ,Gpplus ,Gpminus ,Gdiff ,prob)

N=length(vect);4 newvect=[-Inf; vect; Inf];

A1=bsxfun (@times ,Gpplus ,0.25*( newvect (3: end)-newvect (2:end -1)));6 A1(isnan(A1)==1) =0;

A2=bsxfun (@times ,Gpminus ,0.25*( newvect (2:end -1)-newvect (1:end -2)));8 A2(isnan(A2)==1) =0;

A=A1+A2-Gdiff;10 d=sum(sum(bsxfun (@times ,A,prob) ,3) ,2);

dsup=sum(sum(bsxfun (@times ,A1(1:N-1,:,:),prob) ,3) ,2);12 M=diag(d)+diag(dsup ,1)+diag(dsup ,-1);

return14 end

%% funzioni ausiliare per il calcolo dell 'Hessiana2 function y=G(scapplus ,scapminus)

y=0.5*( erf(scapplus)-erf(scapminus));4 return

end6 function y=F(scapplus ,scapminus ,mx_ ,sx_)

mx=reshape(mx_ ,1 ,[]);8 sx=reshape(sx_ ,1,size(sx_ ,1),size(sx_ ,2));

A1=1/( sqrt (2*pi))*bsxfun (@times ,(exp(-scapminus .^2)-exp(-scapplus .^2)),sx);

10 A2 =0.5* bsxfun (@times ,mx ,(erf(scapplus)-erf(scapminus)));y=A1+A2;

12 returnend

14 function y=Gp(x,sx_)sx=reshape(sx_ ,1,size(sx_ ,1),size(sx_ ,2));

16 y=1/( sqrt (2*pi))*bsxfun (@rdivide ,exp(-x.^2),sx);return

18 endfunction y=scap(x,mx_ ,sx_)

20 mx=reshape(mx_ ,1 ,[]);sx=reshape(sx_ ,1,size(sx_ ,1),size(sx_ ,2));

22 y=1/( sqrt (2))*bsxfun (@rdivide ,bsxfun (@minus ,x,mx),sx);return

24 end

%% medie e varianze per S (mux , six) e V (muy , siy)2

function y=mux(s_past ,deltat ,r)4 y=s_past .*(1+r*deltat);

return6 end

8 function y=muy(v_past ,deltat ,kappa ,theta)y=v_past+kappa *(theta -v_past)*deltat;

10 returnend


12function y=six(s_past ,v_past ,deltat)

14 vpast=reshape(v_past ,1,[]);y=sqrt(deltat)*bsxfun (@times ,s_past ,sqrt(abs(vpast)));

16 returnend

18function y=siy(v_past ,deltat ,csi)

20 y=csi*sqrt(abs(v_past))*sqrt(deltat);return

22 end

%% algoritmo di N-R per il primo passo2 function [vect , popt , popt_double ]=...

nr_fastinit (~,vect ,v_now ,s0,r,param ,deltat ,toll)4

csi=param (1);6 kappa=param (2);

rho=param (3);8 theta=param (4);

v0=param (5);10 N=length(vect);

12 mx=mux(s0 ,deltat ,r);my=muy(v0 ,deltat ,kappa ,theta);

14 mu=[mx , my];sx=six(s0 ,v0 ,deltat);

16 sy=siy(v0 ,deltat ,csi);

18 sig=[sx.^2, rho.*sx.*sy; rho.*sx.*sy , sy.^2];

20 vect_cp=zeros(N,1);iter =0;

22while norm(vect -vect_cp)>toll

24 iter=iter +1;vnew=[-Inf; vect; Inf];

26 vp =0.5*( vnew (3:end)+vnew (2:end -1));vm =0.5*( vnew (2:end -1)+vnew (1:end -2));

28 scapplus =(vp -mx)./( sqrt (2)*sx);scapminus =(vm-mx)./( sqrt (2)*sx);

30 Gpplus=exp(-scapplus .^2) ./( sqrt (2*pi)*sx);Gpminus=exp(-scapminus .^2) ./( sqrt (2*pi)*sx);

32 Gdiff =0.5*( erf(scapplus)-erf(scapminus));Fdiff=sx/(sqrt (2*pi))*(exp(-scapminus .^2)-exp(-scapplus .^2))+...

34 0.5*mx*(erf(scapplus)-erf(scapminus));vect_cp=vect;

36 gdist=Fdiff -vect.* Gdiff;dp=[ Gpplus (1:end -1) .*0.25.*( vnew (3:end -1)-vnew (2:end -2));0]+[0;

Gpminus (2:end).*0.25.*( vnew (3:end -1)-vnew (2:end -2))]-Gdiff;38 dsup=Gpplus (1:end -1) .*0.25.*( vnew (3:end -1)-vnew (2:end -2));

40 A=diag(dp)+diag(dsup ,1)+diag(dsup ,-1);

59


44 breakA=pinv(A);

46 vect=vect - A*gdist;else

48 [L,U]=lu(A);vect=vect - U\(L\( gdist));

50 endvect=sort(vect);

52 end

54 % calcolo la grigliavnew=[-Inf; vect; Inf];

56 sp =0.5*( vnew (3:end)+vnew (2:end -1));sm =0.5*( vnew (2:end -1)+vnew (1:end -2));

58 scapplus =(sp -mx)./( sqrt (2*pi)*sx);scapminus =(sm-mx)./( sqrt (2*pi)*sx);

60 Gdiff=sqrt(1-rho ^2) /(2*sy)*(erf(scapplus)-erf(scapminus));popt=Gdiff;

62 vnew=[-Inf; v_now; Inf];vp =0.5*( vnew (3:end)+vnew (2:end -1));

64 vm =0.5*( vnew (2:end -1)+vnew (1:end -2));N2=length(vp);

66 popt_double=zeros(N,N2);S__1=reshape(kron(sx,ones(N*N2 ,1)) ,[],1);

68 S__2=reshape(kron(sy,ones(N*N2 ,1)) ,[],1);S__3=[S__1 ,S__2];

70 [AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sp ,vp ,mx ,my);BB1=AA1 -AA3; BB2=AA2 -AA4;

72 CC1=reshape ([ reshape(BB1 ,1,[]),reshape(BB2 ,1 ,[])],[],2);DD1=internal.stats.bvncdf(CC1./S__3 ,rho ,1e-8);

74 [AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sm ,vm ,mx ,my);BB1=AA1 -AA3; BB2=AA2 -AA4;


78 [AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sp ,vm ,mx ,my);BB1=AA1 -AA3; BB2=AA2 -AA4;


82 [AA1 ,AA2 ,AA3 ,AA4]= ndgrid(sm ,vp ,mx ,my);BB1=AA1 -AA3; BB2=AA2 -AA4;


86 DDtot=DD1+DD2 -DD3 -DD4;popt_double=reshape(DDtot ,N,N2);

88 returnend

1 %% algoritmo di N-R per i passi successivi

3 function [vect , popt , popt_double ,ptrans ]=...nr_fast(~,vect ,v_now ,s_past ,v_past ,prob_past ,r, param ,deltat ,toll)

5


csi=param (1);7 kappa=param (2);

rho=param (3);9 theta=param (4);

v0=param (5);11

N=length(vect);13 prob=reshape(prob_past ,1,size(prob_past ,1),size(prob_past ,2));

15 mx=mux(s_past ,deltat ,r);my=muy(v_past ,deltat ,kappa ,theta);

17 sx=six(s_past ,v_past ,deltat);sy=siy(v_past ,deltat ,csi);

19vect_cp=zeros(N,1);

21 iter =0;

23 while norm(vect -vect_cp)>tolliter=iter +1;

25 snew=[-Inf; vect; Inf];vp =0.5*( snew (3:end)+snew (2:end -1));

27 vm =0.5*( snew (2:end -1)+snew (1:end -2));scapplus=scap(vp ,mx ,sx);

29 scapminus=scap(vm ,mx ,sx);Gpplus=Gp(scapplus ,sx);

31 Gpminus=Gp(scapminus ,sx);Gdiff=G(scapplus ,scapminus);

33 Fdiff=F(scapplus ,scapminus ,mx ,sx);vect_cp=vect;

35 gdist=g_dist(vect ,Fdiff ,Gdiff ,prob);A=h_dist(vect ,Gpplus ,Gpminus ,Gdiff ,prob);


39 breakA=pinv(A); vect=vect - A*gdist;

41 else[L,U]=lu(A);

43 vect=vect - U\(L\( gdist));end

45 vect=sort(vect);end

47%% calcolo griglia

49 snew=[-Inf; vect; Inf];sp =0.5*( snew (3:end)+snew (2:end -1));

51 sm =0.5*( snew (2:end -1)+snew (1:end -2));scapplus=scap(sp ,mx ,sx);

53 scapminus=scap(sm ,mx ,sx);Gdiff=G(scapplus ,scapminus);

55 %% calcolo pesipopt=sum(sum(bsxfun (@times ,Gdiff ,prob) ,3) ,2);

57 [popt_double ,ptrans ]= prob_double(vect ,v_now ,mx,my ,sx ,sy ,prob_past ,rho);return

59 end

Appendice B

Codici Quantizzazione con Fourier

Modello di Black and Scholes


1 %% parametri modelloS0 =100; r=0.04; sigma =0.2;

3%% parametri opzione

5 K=90:2.5:100; T=1;

7 %% parametri metodo numericoN=150; max_int =20; toll=1e-5;

9%% prezzo esatto

11 exactprice=blsprice(S0 ,K,r,T,sigma);

13 %% QUANTIZZAZIONE+FFTcall_quant_char=bs_quant_char(S0,K,r,T,sigma , N, max_int , toll);

15disp(' benchmark - quant + fourier - Errore relativo (%)')

17 [exactprice ' call_quant_char abs(( exactprice '-call_quant_char)./call_quant_char)*100]

61

62 APPENDICE B. CODICI QUANTIZZAZIONE CON FOURIER

BS QUANT CHAR

1 function call=bs_quant_char(S0 ,K,r,T,sigma , N, max_int , toll)

3 %% grigliaa=100; b=110;

5 x=(a:(b-a)/(N-1):b);x_ =[0.001 x 1000];

7 y=ones(1,N);

9 %% Newthon Raphsonwhile norm(y,2)>toll

11 xp =0.5*( x_(2:N+1)+x_(3:N+2));xm =0.5*( x_(1:N)+x_(2:N+1));

13 [y,A]= hess_der(xm ,xp ,x,max_int , T, r,sigma , S0);[L,U]=lu(A);

15 x=sort(x - (U\(L\(y')))');x_ =[0.001 x 1000];

17 endxp =0.5*( x_(2:N+1)+x_(3:N+2));

19 xm =0.5*( x_(1:N)+x_(2:N+1));

21 %% calcolo pesiG_=@(u) G(u,xm ,xp , T,r, sigma , S0);

23 prob_finale =1/pi*integral(G_ ,0,max_int ,'ArrayValued ',true);

25 %% pricing callsize=length(K);

27 call=zeros(size ,1);for i=1: size

29 call(i,1)=exp(-r*T).*(( max(x-K(i) ,0))*prob_finale ');end

31 end

63

FUNZIONI AUSILIARIE (BS QUANT CHAR)

1 % Calcolo matrice Hessiana e Jacobianafunction [y,M]= hess_der(xm ,xp ,x,max_int , T,r, sigma , S0)

3 funct =@(u) funz(u,xm,xp ,x, T, r,sigma , S0);integrals=integral(funct ,0,max_int ,'ArrayValued ',true);

5 part=integrals (1,:);d_up_p=integrals (2,:);

7 d_down_p=integrals (3,:);y=integrals (4,:);

9 d=d_up_p+d_down_p -part;d_up=d_up_p (1:end -1);

11 d_down=d_down_p (2:end);M=diag(d)+diag(d_up ,1)+diag(d_down ,-1);

13 returnend

% Funzione ausiliaria per i il calcolo di Jacobiana ed Hessiana2 function y=funz(u,xm ,xp ,x, T,r, sigma , S0)

PHI=PHI_T(u, T, r,sigma , S0);4 ep=exp(-1i*u*log(xp)).*PHI;

em=exp(-1i*u*log(xm)).*PHI;6 g= real((em-ep)./(1i*u));

dlp =0.5.* real(ep);8 dlm = -0.5.* real(em);

dgp=dlp./xp;10 dgm=dlm./xm;

l=real((xp.*ep-xm.*em)./(1 -1i.*u));12 y=[g; dlp -x.*dgp; dlm -x.*dgm; l-x.*g];

return14 end

% Funzione caratteristica risk - neutral di B&S2 function y=PHI_T(u, T, r,sigma , S0)

y=exp(T*(PSI_T(u, sigma)+1i*u*(r-PSI_T(-1i,sigma))) + 1i.*u*log(S0));4 return

end

1 % Esponente caratteristico di B&Sfunction y=PSI_T(u, sigma)

3 y= -0.5* sigma ^2*1i*u -0.5* sigma ^2*u.^2;return

5 end

1 % Funzione ausiliaria per il calcolo dei pesifunction y=G(u,xm,xp , T, r,sigma , S0)

3 y= real((exp(-1i*u*log(xm))-exp(-1i*u*log(xp)))*PHI_T(u, T, r,sigma , S0)./(1i*u));

return5 end


Modello KOU


1 %% PARAMETRIS0 =100; r=0.01;

3 sigma =0.2; lambda =0.5; lambdap =50; lambdam =10; p=0.6;param =[sigma lambda lambdap lambdam p];

5 K=90:2.5:100; T=1;N=150; max_int =20; toll=1e-5;

7%%Quantizzazione + Fourier

9 ticcall_quant_char=kou_quant_char(K,r,T,S0,param , N,max_int , toll);

11 toc

13 % FFT - BENCHMARKcall_fft=KOU_fft(K, S0 , T, r, param);

15disp(' benchmark - quant + fourier - Errore relativo (%)')

17 [call_fft call_quant_char abs(( call_fft - call_quant_char)./ call_fft)*100]

KOU QUANT CHAR

1 function call=kou_quant_char(K,r,T,S0,param , N,max_int , toll)

3 %% grigliaa=90; b=110;

5 x=(a:(b-a)/(N-1):b);x_ =[0.001 x 1000];

7 y=ones(1,N);


11 xp =0.5*( x_(2:N+1)+x_(3:N+2));xm =0.5*( x_(1:N)+x_(2:N+1));

13 [y,A]= hess_der(xm ,xp ,x,max_int ,S0 ,T,r,param);[L,U]=lu(A);

15 x=sort(x - (U\(L\(y')))');x_ =[0.001 x 1000];

17 endxp =0.5*( x_(2:N+1)+x_(3:N+2));

19 xm =0.5*( x_(1:N)+x_(2:N+1));

21 %% calcolo pesiG_=@(u) G(u,xm ,xp ,S0 ,T,r,param);

23 prob_finale =1/pi*integral(G_ ,0,max_int ,'ArrayValued ',true);size=length(K);


27 call(i)=exp(-r*T).*(( max(x-K(i) ,0))*prob_finale ');end

29 end

65

FUNZIONI AUSILIARIE (KOU QUANT CHAR)

1 % Calcolo matrice Hessiana e Jacobianafunction [y,M]= hess_der(xm ,xp ,x,max_int ,S0,T,r,param)

3 funct =@(u) funz(u,xm,xp ,x,S0 ,T,r,param);integrals=integral(funct ,0,max_int ,'ArrayValued ',true);





13 end

1 function y=funz(u,xm ,xp ,x,S0 ,T,r,param)PHI=PHI_T(u,S0,T,r,param);

3 ep=exp(-1i*u*log(xp)).*PHI;em=exp(-1i*u*log(xm)).*PHI;

5 g= real((em-ep)./(1i*u));dlp =0.5.* real(ep);

7 dlm = -0.5.* real(em);dgp=dlp./xp;

9 dgm=dlm./xm;l=real((xp.*ep-xm.*em)./(1 -1i.*u));

11 y=[g; dlp -x.*dgp; dlm -x.*dgm; l-x.*g];end

% Funzione caratteristica (risk neutral) KOU2 function y=PHI_T(u,S0 ,T,r,param)

sigma=param (1);4 k=param (2);

theta=param (3);6 y=exp(T*(PSI_T(u, param)+1i*u*(r-PSI_T(-1i,param))) + 1i.*u*log(S0));

end

1 % Esponente caratteristico KOUfunction y=PSI_T(u, param)

3 sigma=param (1);lambda=param (2);

5 lambdap=param (3);lambdam=param (4);

7 p=param (5);y= -0.5* sigma ^2*u.^2+1i*u*lambda .*(p./( lambdap -1i*u) -(1-p)./( lambdam +1i*

u));9 end

1 % Funzione ausiliaria per il calcolo dei pesifunction y=G(u,xm,xp ,S0 ,T,r, param)

3 y= real((exp(-1i*u*log(xm))-exp(-1i*u*log(xp)))*PHI_T(u,S0 ,T,r,param)./(1i*u));

end


Modello Variance Gamma


%% PARAMETRI2 sigmaVG =0.1243; kVG =0.1686; thetaVG =0.0436;

param =[ sigmaVG kVG thetaVG ];4 r=0.02; S0=100;

T=1; K=90:2.5:100;6 N=120; max_int =40; toll=1e-4;

8 %% FFT - BENCHMARKprice_fft=VG_fft(K, S0 , T, r, param);

10%% QUANTIZZAZIONE + FOURIER

12 ticprice_quant_char=VG_quant_char(K,r,T,S0 ,param , N , max_int , toll);

14 time=toc

16 disp(' fft - quant + fourier - Relative error (%)')[price_fft price_quant_char abs((price_fft -price_quant_char)./ price_fft

)*100]

VG QUANT CHAR

1 function call=VG_quant_char(K,r,T,S0,param , N,max_int , toll)%% griglia

3 a=90; b=110;x=(a:(b-a)/(N-1):b); x_ =[0.001 x 1000];

5 y=ones(1,N);


9 xp =0.5*( x_(2:N+1)+x_(3:N+2));xm =0.5*( x_(1:N)+x_(2:N+1));

11 [y,A]= hess_der(xm ,xp ,x,max_int ,S0 ,T,r,param);[L,U]=lu(A);

13 x=sort(x - (U\(L\(y')))');x_ =[0.001 x 1000];

15 endxp =0.5*( x_(2:N+1)+x_(3:N+2));

17 xm =0.5*( x_(1:N)+x_(2:N+1));

19 %% calcolo pesiG_=@(u) G(u,xm ,xp ,S0 ,T,r,param);





29 end

67

FUNZIONI AUSILIARIE (VG QUANT CHAR)

1 % Calcolo matrice Hessiana e Jacobianafunction [y,M]= hess_der(xm ,xp ,x,max_int ,S0,T,r,param)

3 funct =@(u) funz(u,xm,xp ,x,S0 ,T,r,param);integrals=integral(funct ,0,max_int ,'ArrayValued ',true);





13 end

1 function y=funz(u,xm ,xp ,x,S0 ,T,r,param)PHI=PHI_T(u,S0,T,r,param);

3 ep=exp(-1i*u*log(xp)).*PHI;em=exp(-1i*u*log(xm)).*PHI;

5 g= real((em-ep)./(1i*u));dlp =0.5.* real(ep);

7 dlm = -0.5.* real(em);dgp=dlp./xp;

9 dgm=dlm./xm;l=real((xp.*ep-xm.*em)./(1 -1i.*u));

11 y=[g; dlp -x.*dgp; dlm -x.*dgm; l-x.*g];end

% Funzione caratteristica (risk neutral) VG2 function y=PHI_T(u,S0 ,T,r,param)


theta=param (3);6 y=((1+0.5*k*sigma ^2*(u.^2)- 1i*theta*k*u).^(-T/k))...

.*exp(1i*u*T*(r-PSI_T(-1i,param))+ 1i*u*log(S0));8 end

% Esponente caratteristico VG2 function y=PSI_T(u, param)


theta=param (3);6 y=-1/k*log (1+0.5*k*u.^2* sigma^2-1i*theta*k*u);

end

1 % Funzione ausiliaria per il calcolo dei pesifunction y=G(u,xm,xp ,S0 ,T,r, param)

3 y= real((exp(-1i*u*log(xm))-exp(-1i*u*log(xp)))*PHI_T(u,S0 ,T,r,param)./(1i*u));

end


Modello di Heston


%% PARAMETRI2 T=1; K=90:2.5:100;

S0 =100; v0 =0.0456; r = 0.04;4 csi =0.5087; rho = -0.7495; kappa =1.7029; theta =0.0600;

param =[csi kappa rho theta v0];6 N=120; max_int =40; toll=1e-4;

8 %% quantizzazione + ffttic

10 Call_quant_char=heston_quant_char(K,r,T,S0, param , N,max_int ,toll);toc

12%% FFT - BENCHMARK

14 Call_fft=heston_fft(param ,S0 ,K,r,T);

16 disp(' BENCHMARK - QUANT+FOURIER - ERROR')[Call_fft ', Call_quant_char , abs((Call_fft ' - Call_quant_char)./

Call_fft ')*100]

HESTON QUANT CHAR

1 function call=heston_quant_char(K,r,T,S0, param , N,max_int , toll)%% griglia

3 a=90; b=110;x=(a:(b-a)/(N-1):b); x_ =[0.001 x 1000];

5 y=ones(1,N);


9 xp =0.5*( x_(2:N+1)+x_(3:N+2));xm =0.5*( x_(1:N)+x_(2:N+1));

11 [y,A]= hess_der(xm ,xp ,x,max_int , r,T,S0 , param);[L,U]=lu(A);

13 x=sort(x - (U\(L\(y')))');x_ =[0.001 x 1000];

15 endxp =0.5*( x_(2:N+1)+x_(3:N+2));

17 xm =0.5*( x_(1:N)+x_(2:N+1));

19 %% calcolo pesiG_=@(u) G(u,xm ,xp , r,T,S0 , param);





29 end

69

FUNZIONI AUSILIARIE (HESTON QUANT CHAR)

1 % Calcolo matrice Hessiana e Jacobianafunction [y,M]= hess_der(xm ,xp ,x,max_int ,r,T,S0,param)

3 funct =@(u) funz(u,xm,xp ,x,r,T,S0 ,param);integrals=integral(funct ,0,max_int ,'ArrayValued ',true);

5 part=integrals (1,:);%diagonale superiore

7 d_up_p=integrals (2,:);%diagonale inferiore

9 d_down_p=integrals (3,:);%diagonale principale


13 d_down=d_down_p (2:end);% Hessiana

15 M=diag(d)+diag(d_up ,1)+diag(d_down ,-1);%Jacobiana

17 y=integrals (4,:);end

% funzione ausiliaria per i il calcolo di Jacobiana ed Hessiana2 function y=funz(u,xm ,xp ,x,r,T,S0 ,param)

PHI=PHI_T(u,r,T,S0,param);4 ep=exp(-1i*u*log(xp)).*PHI;

em=exp(-1i*u*log(xm)).*PHI;6 g=real((em-ep)./(1i*u));

dlp =0.5.* real(ep);8 dlm = -0.5.* real(em);

dgp=dlp./xp;10 dgm=dlm./xm;

l=real((xp.*ep-xm.*em)./(1 -1i.*u));12 y=[g; dlp -x.*dgp; dlm -x.*dgm; l-x.*g];

end

1 % funzione ausiliaria per il calcolo dei pesifunction y=G(u,xm,xp , r,T,S0 , param)

3 y= real((exp(-1i*u*log(xm))-exp(-1i*u*log(xp)))*PHI_T(u,r,T,S0 ,param)./(1i*u));

end

%Funzione caratteristica del modello di Heston2 function y=PHI_T(u, r,T,S0 , param)

csi=param (1); kappa=param (2); rho=param (3); theta=param (4); V0=param(5);

4 M=sqrt((rho.*csi .*1i.*u - kappa).^2 + csi .^2.*(1i.*u + u.^2));N=(-rho.*csi .*1i.*u + kappa - M)./(-rho.*csi .*1i.*u + kappa + M);

6 A=r.*1i.*u.*T + kappa .*theta ./( csi .^2) .*((-rho.*csi.*1i.*u + kappa - M).*T...- 2.* log((1-N.*exp(-M.*T))./(1-N)));

8 C=(-exp(-M.*T)+1).*(-rho.*csi.*1i.*u + kappa - M)./(csi .^2.*(1 -N.*exp(-M.*T)));

y=exp(A + C.*V0 + 1i.*u.*log(S0));10 end

Bibliograa

[1] G. Callegaro, L. Fiorin, M. Grasselli. Quantization Meets Fourier: a New Technology for PricingOptions. Working Paper, available on SSRN, 2017.

[2] T.A. McWalter, R. Rudd, J. Kienitz and E. Platen. Recursive Marginal Quantization of Higher-Order Schemes. Working Paper, available on SSRN, 2017.

[3] T.A. McWalter, R. Rudd, J. Kienitz and E. Platen. Fast Quantization of Stochastic Volatility ModelsWorking Paper, available on SSRN, 2017.URL: https://arxiv.org/abs/1704.06388.

[4] G. Callegaro, L. Fiorin, M. Grasselli. Pricing via Recursive Quantization in Stochastic VolatilityModels. Quantitative Finance, Volume 17, Issue 6, 2017.

[5] G. Pagès and Sagna. Recursive marginal quantization of the Euler scheme of a diusion process.Applied Mathematical Finance, Volume 2, Issue 5, 2015URL: http://arxiv.org/abs/1304.2531.

[6] G. Pagès. Introduction to optimal vector quantization and its applications for numerics, ESAIM:Proceedings and Surveys, EDP Sciences, 48 (1), pp.29-79, 2015.URL: http://hal.archives-ouvertes.fr/INSMI/hal-0103419.

[7] R. Cont e P. Tankov. Financial Modelling with Jump Processes. Chapman e Hall/Crc FinancialMathematics Series, 2012.URL: http://books.google.it/books?id=-fZtKgAACAAJ.

[8] G. Pagès, B. Wilbertz. Optimal Delaunay and Voronoi quantization schemes for pricing Americanstyle options. Springer, 2012.

[9] A. Sagna. Pricing of barrier options by marginal functional quantization.Monte Carlo Methods andApplications, Volume 17, Issue 4, 2011.

[10] E. Barucci. Ingegneria nanziaria. Un'introduzione quantitativa. I Manuali. EGEA, 2009.URL: http://books.google.it/books?id=fWeEPgAACAAJ.

[11] G. Pagès e J. Printems. Optimal quantization for Finance: from random vectors to stochasticprocesses. Handbook of Numerical Analysis, volume 15, 2008.

[12] Leif Andersen. Ecient Simulation of the Heston Stochastic Volatility ModelAndersen, Leif B. G.,Ecient Simulation of the Heston Stochastic Volatility Model (January 23, 2007). Available atSSRN, 2007.URL: https://ssrn.com/abstract=946405

[13] H. Albrecher, P. Mayer, W. Schoutens e J. Tistaert. The little Heston Trap. Wilmott Magazine,2006.URL: https://perswww.kuleuven.be/u0009713/HestonTrap.pdf.

[14] G. Pagès and J. Printems. Functional quantization for numerics with an application to optionpricing. Monte Carlo Methods and Applications, 11(4), 407-446, 2005.

70

BIBLIOGRAFIA 71

[15] V. Bally, G. Pagès, J. Printems. A quantization tree method for pricing and hedgingmultidimensional American options. Mathematical Finance, 15(1):119 168, 2005.

[16] J. Jacod e P.E. Protter. Probability Essentials. Universitext(1979). Springer Berlin Heidelberg. isbn:9783540438717, 2003

[17] S. Graf, H. Luschgy. Foundations of quantization for probability distributions. Springer, 2000.

[18] S.L. Heston. A closed-form solution for options with stochastic volatility with applications to bondand currency options. The Review of Financial Studies, 6(2), pp. 327-343, 1993.URL: http://faculty.baruch.cuny.edu/lwu/890/Heston93.pdf.

[19] Shephard, N. G. From characteristic function to distribution function: A simple framework for thetheory. Econometric Theory, 1991

Quantizzazione: applicazioni al pricing di derivati · Questo lavoro ha l'obiettivo di illustrare...

Documents

Transcript of Quantizzazione: applicazioni al pricing di derivati · Questo lavoro ha l'obiettivo di illustrare...