1 Introduzione

1.1 Due software a confronto

In questi appunti mostrerò come utilizzare R e Python per alcune applicazioni finanziarie, scrivendo a sinistra il codice R e a destra il codice Python. I due linguaggi hanno molte caratteristiche in comune e si diversificano principalmente per la sintassi.

1.2 Operazioni matematiche

Le operazioni di base sono gestite nello steso modo dai due software R e Python, con la differenza di notazione sull’elevamento a potenza.

2+3
## [1] 5
2*3
## [1] 6
2^3
## [1] 8
2/3
## [1] 0.6666667
2^(1/2)
## [1] 1.414214
2+3
## 5
2*3
## 6
2**3
## 8
2/3
## 0.6666666666666666
2**(1/2)
## 1.4142135623730951

Per qualsiasi calcolo matematico avanzato su Python si deve caricare il pacchetto numpy (numerical python) al quale si può dare il nome che si preferisce.

Questo è un punto debole di un linguaggio di programmazione. Se ogni utente può cambiare il nome del pacchetto, pezzi di codice potrebbero essere incompatibili tra loro. Nella pratica si stanno sempre di più diffondendo nomi tandard per i pacchetti. Il pacchetto numpy è, quasi ovunque, chiamato “np”.

Se si usa il sfotware RStudio/Posit per scrivere in Python, allora bisogna usare, all’interno dell’ambiente R, i comandi seguenti: library(reticulate) e py_install('nome pacchetto')

log(2)
## [1] 0.6931472
exp(2)
## [1] 7.389056
import numpy as np
np.log(2)
## np.float64(0.6931471805599453)
np.exp(2)
## np.float64(7.38905609893065)

1.3 Sequenze di numeri

Può essere utile, per alcune applicazioni, creare una sequenza di valori numerici. Per creare tale sequenza bisogna sempre specificare il valore iniziale, il valore finale e il “passo” della serie (cioè di quanto si increment un elemento per ottenere il successivo).

seq(2,5,0.2)
##  [1] 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0
np.arange(2,5,0.2)
## array([2. , 2.2, 2.4, 2.6, 2.8, 3. , 3.2, 3.4, 3.6, 3.8, 4. , 4.2, 4.4,
##        4.6, 4.8])

Notiamo una differenza importante! R fa terminare la sequenza, se possibile, con il numero indicato come valore finale, mentre Python lascia questo numero fuori dalla lista. Questo è l’effetto di un’altra caratteristica di Python che vedremo meglio a breve.

Si possono effettuare operazioni sulle sequenze che, ovviamente, vengono aplicate a tutti i termini della sequenza.

log(seq(2,5,0.2))
##  [1] 0.6931472 0.7884574 0.8754687 0.9555114 1.0296194 1.0986123 1.1631508
##  [8] 1.2237754 1.2809338 1.3350011 1.3862944 1.4350845 1.4816045 1.5260563
## [15] 1.5686159 1.6094379
np.log(np.arange(2,5,0.2))
## array([0.69314718, 0.78845736, 0.87546874, 0.95551145, 1.02961942,
##        1.09861229, 1.16315081, 1.22377543, 1.28093385, 1.33500107,
##        1.38629436, 1.43508453, 1.48160454, 1.5260563 , 1.56861592])

1.4 Lavorare con le matrici

Supponiamo di voler inserire la seguente matrice:

\[A=\left[\begin{array}{cc} 1 & 5\\ 2 & 3\\ 3 & 7 \end{array}\right]\]

In R si inseriscono tutti i numeri colonna per colonna e si indicano, poi, le dimensioni della matrice. In Python si possono inserire i numeri già ordinati.

In R un insieme di qualsiasi valori è sempre indicata con il comando c() che è l’iniziale del verbo “concatenate”.

A=array(c(1,2,3,5,3,7),dim=c(3,2))
A
##      [,1] [,2]
## [1,]    1    5
## [2,]    2    3
## [3,]    3    7
A=np.array([[1,5],[2,3],[3,7]])
A
## array([[1, 5],
##        [2, 3],
##        [3, 7]])

I singoli elementi di una matrice si indicano mediante parentesi quadrate che ne contengono le coordinate (riga e colonna), eparate da una virgola.

Inoltre:

  • se si indicano, per una coordinata, due numeri separati da “due punti”, allora si considerano tutti i numeri compresi nell’intervallo delle coordinate
  • se non si indica una dimensione, si intende che devono essere presi tutti gli elemenit di quella dimenione
  • si posono indicare due elementi separati

Python ha una caratteristica importante, che non ha in comune con nessun altro software.

Le coordinate di un vettore o di una matrice partono con l’emento di ordine “zero”. Quindi A[1] rappresenta il secondo elemento del vettore A.

A[3,2]
## [1] 7
A[2:3,2]
## [1] 3 7
A[2:3,]
##      [,1] [,2]
## [1,]    2    3
## [2,]    3    7
A[c(1,3),]
##      [,1] [,2]
## [1,]    1    5
## [2,]    3    7
A[2,1]
## np.int64(7)
A[1:3,1]
## array([3, 7])
A[1:3,:]
## array([[2, 3],
##        [3, 7]])
A[[0,2],:]
## array([[1, 5],
##        [3, 7]])

1.5 Operazioni tra matrici

Sia in R sia in Python le operazioni tra matrici avvengono “elemento per elemento” e se si vuole effettuare un calcolo matriciale occorre indicarlo espressamente. Di seguito vediamo come fare la seguente operazione

\[\left[\begin{array}{cc} 1 & 5\\ 2 & 3 \end{array}\right]\left[\begin{array}{c} 3\\ 7 \end{array}\right]\]

X=A[1:2,1:2] #matrice
y=A[2:3,2] #vettore

# Prodotto standard
X*y
##      [,1] [,2]
## [1,]    3   15
## [2,]   14   21
# Prodotto matriciale
X%*%y
##      [,1]
## [1,]   38
## [2,]   27
X=A[0:2,0:2] #matrice
y=A[1:3,1] #vettore

# Prodotto standard
X*y
## array([[ 3, 35],
##        [ 6, 21]])
# Prodotto matriciale
np.matmul(X,y)
## array([38, 27])

L’inverione di una matrice su R avviene mediante il comando solve (questo proviene dal fatto che invetire una matrice è l’operazione con cui si risolve un sistema di equazioni lineari). Se su R si applica l’elevamento a potenza ^(-1) su una matrice, tutti gli elementi vengono invertiti (uno a uno).

In Python per invertire una matrice bisogna usare un sottopacchetto di numpy. Questo si chiama linalg (intendendo, ovviamente, “linear albegra”) e, al suo interno, si richiama la funzione “inv”.

solve(X)
##            [,1]       [,2]
## [1,] -0.4285714  0.7142857
## [2,]  0.2857143 -0.1428571
X%*%solve(X)
##      [,1]         [,2]
## [1,]    1 1.110223e-16
## [2,]    0 1.000000e+00
np.linalg.inv(X)
## array([[-0.42857143,  0.71428571],
##        [ 0.28571429, -0.14285714]])
np.matmul(X,np.linalg.inv(X))
## array([[1.00000000e+00, 5.55111512e-17],
##        [0.00000000e+00, 1.00000000e+00]])

1.6 Generazione di numeri casuali

In R e Python si possono generare numeri casuali da molte distribuzioni. Bisogna utilizzare un comando diverso per ogni distribuzione e, al suo interno, indicare i parametri della distribuzione e il numero di “estrazioni” da generare.

Nel caso di una distribuzione normale, per esempio, oltre al numero di estrazioni bisogna anche indicare la media e la varianza (i due parametri della distribuzione).

Nei due software, sfortunatamente, l’ordine degli input è invertito.

rnorm(n=10,mean=0,sd=1)
##  [1]  0.11822497  0.10227075 -0.04163615 -0.06913379 -0.07062279 -0.39402102
##  [7] -0.09787932  1.83095972 -0.19537799  1.18599659
np.random.normal(loc=0,scale=1,size=10)
## array([ 1.24798524,  0.29866558,  1.00488371, -0.0037097 , -0.6582682 ,
##         0.81954773,  0.49922077,  0.15522686,  0.61829483,  0.64124149])

Quando in una funzione si mantiene l’ordine dato dai programmatori, allora si può anche omettere il nome dell’input che si sta per inserire. Se l’ordine viene cambiato, invece, i nomi sono necessari.

In R, per esempio, si può scrivere rnorm(10,0,1) e, quindi, si genereranno 10 valori estratti da una normale con media 0 e varianza 1, oppure si può scrivere rnorm(mean=0,sd=1,n=10).

Esistono dei valori che sono attribuiti di default ai parametri. Questi valori sono specificati nelle descrizioni delle funzioni.

Nel caso di Python, per esempio, la funzione di creazioni di numeri casuali ha la forma random.normal(loc=0.0, scale=1.0, size=None). In questo caso, quindi, non indicando né il primo né il secondo elemento, Python dà per scontato che si vogliano dare, rispettivamente, i valori 0 e 1. Il valore di size, invece, va specificato perché per default non assume alcun valore.

2 Rappresentazione grafica dei dati

R è già equipaggiato per creare dei grafici, mentre in Python bisogna installare e richiamare uno specifico pacchetto.

Vediamo, quindi, come creare e rappresentare 10’000 valori da una normale standard.

N=rnorm(10000)
plot(N,type='l')

# il "tipo" di grafico è "l", cioè "a linee"
# per default R crea un grafico a punti
import matplotlib.pyplot as plt
N=np.random.normal(size=10000)

# Si apre la figura
plt.figure()
# Si danno i comandi
plt.plot(N)
# Si chiude la figura
plt.show()

Vediamo come aggiungere al grafico le linee orizontali della media e della media aumentata e diminuita dello scarto quadratico medio. Queste linee saranno tratteggiate, di colore diverso e di spessore diverso.

plot(N,type='l',xlab='estrazioni',main='10000 estrazioni da una normale')
# xlab=etichetta dell'asse x
# main=titolo del grafico
abline(h=mean(N),lty=2,lwd=2,col='blue')
# lty=line type (1 = continuo, 2 = tratteggiato, 3 = puntato)
# lwd=line width
abline(h=mean(N)+sd(N),lty=2,lwd=2,col='red')
abline(h=mean(N)-sd(N),lty=2,lwd=2,col='red')

plt.figure()
plt.plot(N)
plt.xlabel('estrazioni')
plt.title('10000 estrazioni da una normale')
plt.axhline(y = np.mean(N),color='blue',linestyle = '--',linewidth=3)
plt.axhline(y = np.mean(N)+np.std(N),color='red',linestyle = '--',linewidth=3)
plt.axhline(y = np.mean(N)-np.std(N),color='red',linestyle = '--',linewidth=3)
plt.show()

Se vogliamo rappresentare questi dati mediante un istrogramma i comandi sono i seguenti.

hist(N)

plt.figure()
plt.hist(N)
plt.show()

Possiamo cambiare i parametri dell’istogramma. Il più significativo riguarda l’unità di misura dell’asse delle ordinate. Essa può essere la frequenza (come nel caso di default) oppure la frequenza relativa (che deve sommare a 1 su tutto il dominio). Per ottenre la frequenza relativa il comando è il seguente.

hist(N,freq=F)

plt.figure()
plt.hist(N,density='true')
plt.show()

Si può rapperesentare sul grafico una densità della funzione normale di cui conosciamo la formula algebrica. Si crea una sequenza di valori (per esempio da -4 a 4 con un certo intervallo) e poi si genera la funzione normale.

x=seq(-4,4,0.1)
Gauss=1/sqrt(2*pi)*exp(-x^2/2)
hist(N,xlim=c(-4,4),freq = F,breaks = 100)
lines(x,Gauss,col='red')

x=np.arange(-4,4,0.1)
Gauss=1/np.sqrt(2*np.pi)*np.exp(-x**2/2)
plt.figure()
plt.hist(N,bins=x,density='true')
plt.plot(x,Gauss,color='red')
plt.show()

3 Caricare dati esterni

3.1 Dati da Yahoo

In R e Python esistono pacchetti dedicati a connettere il software con database esterni in modo che i dati possano essere scaricati direttamente.

Esistono pacchetti per i più comuni database gratuiti: Yahoo, FRED, OECD, World Bank, Eurostat.

In R il pacchetto quantmod consente di scaricare dati da Yahoo e da FRED. In Python i due database sono accessibili con due pacchetti diversi. Per Yahoo si usa yfinance.

Se il pacchetto viene installato per la prima volta, in R il comando è install.packages('...') mentre per Python useremo, all’interno di un ambiente R, il pacchetto reticulate e, poi, il comando py_install('pacchetto').

#install.packages('quantmod')
library(quantmod)
## Caricamento del pacchetto richiesto: xts
## Caricamento del pacchetto richiesto: zoo
## 
## Caricamento pacchetto: 'zoo'
## I seguenti oggetti sono mascherati da 'package:base':
## 
##     as.Date, as.Date.numeric
## Caricamento del pacchetto richiesto: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
import yfinance as yf

Questi pacchetti permettono di caricare dati direttamente da Yahoo. Vediamo, per esempio, come richiamare i valori dell’indice S&P500 e FTSE.MIB. Sul sito di Yahoo prendiamo i codici degli indici che sono ^GSPC e FTSEMIB.MI.

In R chiediamo che i dati vengano scaricati nella forma di una serie “zoo”. Si tratta di un pacchetto che associa dati e date e che consente di lavorare in modo molto efficiente con serie storiche.

getSymbols(c('^GSPC','FTSEMIB.MI'), # simboli da Yahoo
           src='yahoo', # database sorgente
           return.class='zoo', # tipo di dati in output
           from='1995-01-01', # data di inizio
           to='2019-12-31', # data di fine
           periodicity = 'daily') # per default
## Warning: FTSEMIB.MI contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.
## [1] "GSPC"       "FTSEMIB.MI"
S=yf.download(['^GSPC','FTSEMIB.MI'], # simboli da Yahoo
              start="1995-01-01", # data di inizio
              end='2019-12-31', # data di fine
              interval='1d') # per default
## <string>:1: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [                       0%                       ][*********************100%***********************]  2 of 2 completed

I dati scaricati con R sono archiviati in due variabili separate, mentre su Python appartengono già a un unico data frame.

Se ci interessano solo i valori delle chiusure aggiustate, possiamo selezionarle e rappresentarle graficamente con i comandi seguenti.

In R il comando merge permette di fondere tra loro due o più serie storiche (confrontadole ripetto alle date). Con l’opzione all=FALSE si eliminano dalla nuova serie tutte le date che non hanno osservazioni comuni.

In R si utilizza il simbolo del dollaro $ per identificare una particolare colonna di un data frame. In Python si può utilizzare un punto se la colonna ha un nome “semplice”, ma se nel nome ci sono spazi o caratteri particolare, allora bisogna utilizzare le parentesi quadre (come nel caso seguente).

S=merge(GSPC$GSPC.Adjusted,
        FTSEMIB.MI$FTSEMIB.MI.Adjusted,all=FALSE)
plot(S)

S=S['Close']
plt.figure()
plt.plot(S)
plt.show()

Vediamo, così, che Python pone le due variabili sullo stesso grafico, mentre R le mette in grafici sovrapposti.

Per fare un grafico sovrapposto anche in R occorre specificare che si vuole utilizzare un solo “schermo” per le serie storiche.

plot(S,screens=1,col=c('red','blue'),
     xlab='',ylab='')
legend('topright',legend=c('S&P500','FTSE MIB'),
       lty=1,col=c('red','blue'),bty='n')
grid()

plt.figure()
plt.plot(S['^GSPC'],color='red',label='S&P500')
plt.plot(S['FTSEMIB.MI'],color='blue',label='FTSE MIB')
plt.legend()
plt.grid()
plt.show()

I dati appartengono a due scale diverse. Per poterli confrontare in modo più “chiaro” possiamo trasformare le serie storiche in modo che alla data iniziale entrambe valgano 1. In questo modo possiamo vedere comi si sarebbe evoluta un’unità monetaria investita in entrambi gli indici.

Il primo passaggio è quello di eliminare tutte le date per la quali non ci sono valori in comune ai due indici. Si fa, in R, con il comando na.omit e in Python con il comando dropna.

# Togliamo i valori non disponibili
S=na.omit(S)

# Dividiamo le due colonne per il primo dato
S[,1]=S[,1]/S[[1,1]]
S[,2]=S[,2]/S[[1,2]]

plot(S,screens=1,col=c('red','blue'),xlab='',ylab='')
legend('topleft',legend=c('S&P500','FTSE MIB'),
       lty=1,col=c('red','blue'),bty='n')
grid()

# Togliamo i valori non disponibili
S=S.dropna()

# Dividiamo le due colonne per il primo dato
S['^GSPC']=S['^GSPC']/S['^GSPC'][0]
## <string>:3: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
S['FTSEMIB.MI']=S['FTSEMIB.MI']/S['FTSEMIB.MI'][0]
## <string>:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
# Si usa la funzione '.index' per usare sulle ascisse
# i valori delle date
plt.figure()
plt.plot(S['^GSPC'],color='red',label='S&P500')
plt.plot(S['FTSEMIB.MI'],color='blue',label='FTSE MIB')
plt.legend()
plt.grid()
plt.show()

Dai grafici vediamo che fino al 2010 i due indici si muovono con una fortissima correlazione. Dopo la crisi dei debiti sovrani in Europa, l’indice italiano non cresce più, mentre quello USA riprende la crescita.

3.2 Indicazione di date sui grafici

Per analizzare le serie storiche può essere utile inserire linee verticali che identifichino date significative.

Sui dati scaricati, per esempio, può essere opportuno individuare le segeunti date:

  • 11 settembre 2001: l’attentato alle Torri Gemelle di New York
  • 15 settembre 2008: il fallimento di Lehman Brothers
  • 27 aprile 2010: il declassamento di rating dei titoli di Stato emessi dalla Grecia
# Definizione delle date
date=as.Date(c('2001-09-11',
               '2008-09-15',
               '2010-04-27'))

# Grafico
plot(S,screens=1,col=c('blue','red'),
     xlab='',ylab='')
abline(v=date,lty=2)

plt.figure()
plt.plot(S)
plt.axvline(np.datetime64('2001-09-11'), color='black', linestyle='dotted')
plt.axvline(np.datetime64('2008-09-15'), color='black', linestyle='dotted')
plt.axvline(np.datetime64('2010-04-27'), color='black', linestyle='dotted')
plt.show()

A questo punto possiamo inserire delle scritte sui grafico così da indicare a quali eventi si riferiscan le date. Il comando prevede che si diano le coordinate del punto in cui inserire il testo.

Per mandare il testo a capo si usa il comando \n all’interno della stringa di caratteri.

testo=c('Twin \nTowers',
        'Lehman \nBrothers',
        'Crisi \ngreca')
plot(S,screens=1,col=c('blue','red'),
     xlab='',ylab='')
abline(v=date,lty=2)
text(date,max(S,na.rm=TRUE),labels=testo,pos=1)

plt.figure()
plt.plot(S)
plt.axvline(np.datetime64('2001-09-11'), color='black', linestyle='dotted')
plt.axvline(np.datetime64('2008-09-15'), color='black', linestyle='dotted')
plt.axvline(np.datetime64('2010-04-27'), color='black', linestyle='dotted')

plt.text(np.datetime64('2001-09-11'),np.max(np.max(S)), 'Twin \ntowers', ha='center',va='top')
plt.text(np.datetime64('2008-09-15'),np.max(np.max(S)), 'Lehman \nBrothers', ha='center',va='top')
plt.text(np.datetime64('2010-04-27'),np.max(np.max(S)), 'Crisi \ngreca', ha='center',va='top')
plt.show()

3.3 Dati da FRED

Per scaricare dati da FRED si può utilizzare, in R, ancora il pacchetto quantmod oppure si può ricorrere a un nuovo pacchetto. In Python, invece, bisogna ricorrere a un altro pacchetto.

I pacchetti dedicati a FRED funzionano dopo avere effettuato una registrazione sul sito FRED per ottenere il codice chiave API dalla FED (in queste lezioni inserisco il mio).

# install.packages('fredr')
library(fredr)
fredr_set_key('c5b150d7ef18ec656ac2f4a3541dd60f')
from fredapi import Fred
fred = Fred(api_key='c5b150d7ef18ec656ac2f4a3541dd60f')

Se, per esempio, vogliamo scaricare il tasso di inflazione per Italia e USA (definiti come la variazione mensile, rispetto allo stesso mese dell’anno precedente, dell’indice dei prezzi al consumo) possiamo dare i seguenti comandi, dopo aver individuato i codici delle variabili su FRED: per l’Italia ITACPIALLMINMEI e per gli USA USACPIALLMINMEI.

In R, una volta scaricati i dati, creiamo una serie storica con il comando zoo (all’interno del pacchetto quantmod) indicando come argomenti i dati e le date.

X=fredr('ITACPIALLMINMEI', # Codice FRED
        frequency='m', # Frequenza dei dati
        units='pc1') # Var. % sull'anno precedente
inf_ITA=zoo(X$value,X$date) # Creazione serie storica 

X=fredr('USACPIALLMINMEI',
        frequency='m',
        units='pc1')
inf_USA=zoo(X$value,X$date)

inf=merge(inf_ITA,inf_USA,all=FALSE)
plot(inf,screens=1,col=c('blue','red'),lty=1,lwd=2,
     xlab='',ylab='tasso di inflazione')
legend('topright',legend=c('Italia','USA'),
       col=c('blue','red'),lty=1,lwd=2,bty='n')
grid()

inf_ITA=fred.get_series('ITACPIALLMINMEI', # Codice FRED
                        frequency='m', # Frequenza dei dati
                        units='pc1') # Var. % sull'anno precedente
inf_USA=fred.get_series('USACPIALLMINMEI',
                        frequency='m',
                        units='pc1')

plt.figure()
plt.plot(inf_ITA,label='Italia')
plt.plot(inf_USA,label='USA')
plt.legend()
plt.show()

4 Simulazione di processi stocastici

4.1 Simulazione di un percorso

Nei modelli in tempo continuo il prezzo di un titolo \(S_t\) si evolve secondo un’equazione differenziale stocastica della forma

\[\frac{dS_t}{S_t}=\mu dt+\sigma dW_t\]

dove \(dW_t\) è un processo di Brown (o processo di Wiener) distribuito normalmente con media nulla e varianza \(dt\).

Abbiamo già visto come simulare una variabile normale. Adesso dobbiamo capire come simulare una equazione differenziale stocastica. Il primo passo è quello di discretizzare l’equazione. Ricordando che vale \(dS_t=S_{t+dt}-S_t\), si può scrivere

\[S_{t+dt} = S_t+S_t\mu dt + S_t\sigma dW_t\]

Partendo da un valore inizial \(S_0\), quindi, si può generare il valore del periodo successivo e poi, usando questo, il valore del periodo ancora successivo e così via.

Questo tipo di operazione ricorsiva si effettua mediante un ciclo for.

I dati di cui abbiamo bisogno sono:

  • i parametri del modello (in questo caso \(\mu\) e \(\sigma\));
  • il valore iniziale \(S_0\);
  • il numero di periodi che vogliamo simulare; se simuliamo \(T\) anni suddivii in \(dt\) sottoperiodi, allora bisogna simulare \(\frac{T}{dt}\) osservazioni (se prendiamo valori giornalieri con \(dt=\frac{1}{250}\) e vogliamo simulare due anni \(T=2\), allora le osservazioni sono \(500\)).
S0=10 #valore iniziale
T=5 #anni da simulare
dt=1/250 #frequenza dei dati (giornalieri)
mu=0.08 #tasso di crescita dell'8%
sigma=0.2 #volatilità del processo

S=array(NA,dim=c(T/dt,1)) #creazione di un vettore vuoto per i prezzi simulati
S[1]=S0 #definizione del primo valore del vettore dei prezzi

dW=rnorm(T/dt,0,sqrt(dt)) #definizione del processo di Wiener

# ciclo
for (i in 2:(T/dt)){
  S[i]=S[i-1]*(1+mu*dt+sigma*dW[i-1])
}

plot(S,type='l')

S0=10 #valore iniziale
T=5 #anni da simulare
dt=1/250 #frequenza dei dati (giornalieri)
mu=0.08 #tasso di crescita dell'8%
sigma=0.2 #volatilità del processo

dW=np.random.normal(0, np.sqrt(dt),int(T/dt))
S=np.zeros(int(T/dt))
S[0]=S0

for i in range(1,int(T/dt)):
  S[i]=S[i-1]*(1+mu*dt+sigma*dW[i-1])

# Mostro il grafico delle N simulazioni
plt.figure()
plt.plot(S)
plt.show()

4.2 Simulazione di più percorsi

Possiamo simulare più di un processo e questo ci permetterà di calcolare statitiche sui processi (per esempio la media, la varianza, i quantili).

Nel calcolo precedente si può aggiungere una variabile (\(N\)) che indica quanti percorsi simulare.

Le simulazioni, quindi, non occuperanno più un vettore, bensì una matrice che avrà \(\frac{T}{dt}\) righe e \(N\) colonne.

La procedura più efficace consiste nel generare la matrice una riga per volta (ogni riga indica un istante di tempo).

S0=10 #valore iniziale
T=5 #anni da simulare
N=100 #numero di simulazioni
dt=1/250 #frequenza dei dati (giornalieri)
mu=0.08 #tasso di crescita dell'8%
sigma=0.2 #volatilità del processo

S=array(NA,dim=c(T/dt,N)) #creazione di una matrice vuota per i prezzi simulati
S[1,]=S0 #definizione della prima riga dei prezzi

dW=array(rnorm(T/dt*N,0,sqrt(dt)),dim=c(T/dt,N)) #definizione della matrice dei processi di Wiener

# ciclo
for (i in 2:(T/dt)){
  S[i,]=S[i-1,]*(1+mu*dt+sigma*dW[i-1,])
}

matplot(S,type='l',col='lightgray')

S0=10 #valore iniziale
T=5 #anni da simulare
N=100 #numero di simulazioni
dt=1/250 #frequenza dei dati (giornalieri)
mu=0.08 #tasso di crescita dell'8%
sigma=0.2 #volatilità del processo

dW=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))
S=np.zeros((int(T/dt),N))
S[0,:]=S0

# ciclo
for i in range(1,int(T/dt)):
  S[i,:]=S[i-1,:]*(1+mu*dt+sigma*dW[i-1,:])

plt.figure()
plt.plot(S,color='lightgray')
plt.show()

4.3 Statistiche sulle simulazioni

Possiamo aggiungere al grafico alcune linee che, per esempio, indichino la media, la mediana o i quantili.

Con il comando quantile si possono generare i quantili relativi a ogni livello di probabilità.

Sulle simulazioni vogliamo applicare questo comando alle righe della matrice. Ogni quantile, infatti, va calcolato su ognuna delle righe (per le \(N\) simulazioni delle colonne).

Se, per esempio, vogliamo generare i quantili 0.1, 0.5 (cioè la mediana) e 0.9, il comando è il seguente.

In R e Python gli N quantili applicati alle M righe di una matrice generano una matrice NxM. Per i nostri fini grafici, invece, vogliamo una matrice MxN e, quindi, applichiamo l’operatore “trsposto” t().

Q=apply(S,1,quantile,probs=c(0.1,0.5,0.9))
dim(Q)
## [1]    3 1250
Q=t(Q)
dim(Q)
## [1] 1250    3
Q=np.quantile(S, [0.1,0.5,0.9], axis=1)
np.shape(Q)
## (3, 1250)
Q=np.transpose(Q)
np.shape(Q)
## (1250, 3)

Adesso possiamo aggiungere queste tre curve al grafico delle simulazioni.

matplot(S,type='l',col='lightgray',xlab='giorni')
matlines(Q,lty=1,col=c('red','black','darkgreen'))
legend('topleft',legend=c('Quantile 0.1','Mediana','Quantile 0.9'),
       lty=1,col=c('red','black','darkgreen'),bty='n')
grid()

plt.figure()
plt.plot(S,color='lightgray')
plt.plot(Q[:,0],color='red',label='Quantile 0.1')
plt.plot(Q[:,1],color='black',label='Mediana')
plt.plot(Q[:,2],color='darkgreen',label='Quantile 0.9')
plt.legend()
plt.show()

5 Stima dei parametri di un processo stocastico e previsione del futuro

5.1 Il moto browniano geometrico (MBG)

Dato un processo stocastico nella forma di un moto browniano geometrico (MBG)

\[\frac{dS_t}{S_t}=\mu dt+\sigma dW_t\]

i rendimenti logaritmici (mediante il lemma di Ito) sono dati da

\[d\ln S_t = \left(\mu-\frac{1}{2}\sigma^2\right)dt+\sigma dW_t\]

La media e la varianza di questo processo sono:

\[\mathbb{E}_t\left[d\ln S_t\right]=\left(\mu-\frac{1}{2}\sigma^2\right)dt\] \[\mathbb{V}_t\left[d\ln S_t\right]=\sigma^2 dt\]

Dall’equazione della varianza si ricava lo stimatore di \(\sigma\)

\[\hat\sigma=\sqrt{\frac{\mathbb{V}_t\left[d\ln S_t\right]}{dt}}\]

e una volta ottenuta questa stima, si può ottenere uno stimatore della media \(\mu\)

\[\hat\mu=\frac{\mathbb{E}_t\left[d\ln S_t\right]}{dt}+\frac{1}{2}\hat\sigma^2\]

Questo metodo di calibrazione, appena presentato, si definisce “Metodo dei momenti” e consiste nell’uguagliare i momenti teorici di una distribuzione ai suoi momenti empiric.

5.2 Scaricare e suddividere i dati

Scarichiami i dati dell’indice S&P500 per tutta la durata che abbiamo a disposizione.

In R dobbiamo richiamare il pacchetto quantmod e scaricare da Yahoo l’indice che ha codice “^GSPC”. All’interno dei dati scaricati, poi, ci interessa solo la chisurua aggistata.

getSymbols('^GSPC',src='yahoo',return.class='zoo',from='1920-01-01')
## [1] "GSPC"
SP500=GSPC$GSPC.Adjusted
S=yf.download('^GSPC',start="1920-01-01")
## <string>:1: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [*********************100%***********************]  1 of 1 completed
SP500=S['Close']

Vogliamo dividere i dati in due parti:

  • una parte (in sample) che utilizzeremo per calibrare il modello iniziale;
  • una parte (out of sample) che utilizzeremo per verificare che il modello calibrato sulla prima parte sia adeguato per descrivere anche qualcosa al “di fuori” del campione iniziale.

L’idea è di “far finta” di essere arrivati a una certa data e non conoscere il futuro rispetto a quella data.

Nel caso dell’indice S&P500 vediamo che dal primo gennaio 1995 l’andamento dei prezzi assume un tasso di crescita decisamente più elevato rispetto al passato.

plot(SP500,xlab='',ylab='')
abline(v=as.Date('1995-01-01'),col='red',lty=2)
text(as.Date('1995-01-01'),max(SP500),
     labels=c('in sample <-','-> out of sample'),pos=c(2,4))

plt.figure()
plt.plot(SP500)
plt.axvline(np.datetime64('1995-01-01'), color = 'red',linestyle='dotted')
plt.text(np.datetime64('1995-01-01'),np.max(SP500),'in sample <-',ha='right')
plt.text(np.datetime64('1995-01-01'),np.max(SP500),'-> out of sample',ha='left')
plt.show()

Per creare un sottoinsieme di valori relativi solo ad alcune date si possono effettuare due procedure diverse.

# La funzione "time(x)" restituisce le date relative alla varaibile temporale "x". Prendiamo i dati SP500 prima della fine del 1994:
SP=SP500[time(SP500)<='1994-12-31'] # dati in sample

# Si prendono poi i dati "futuri", successivi al 1994:
SPF=SP500[time(SP500)>'1994-12-31'] # dati out of sample
# La funzione x.index restituisce l'indice della variabile x che, in caso di serie storiche, è la "data"
SP=SP500[SP500.index<'1995-01-01'] # dati in sample
SPF=SP500[SP500.index>='1995-01-01'] # dati our of sample

5.3 Calibrazione “in sample” del modello

Come visto nelle lezioni precedenti, possiamo calibrario i parametri mu e sigma del modello stocastico calcolandi i rendimenti logaritmici dell’indice. Lo facciamo sui dati “in sample”.

dt=1/250 # dati giornalieri
dlnS=diff(log(SP)) # differenze logaritmiche
sigma=sqrt(var(dlnS)/dt)
mu=mean(dlnS)/dt+0.5*sigma^2
mu
## [1] 0.06638521
sigma
## [1] 0.1894119
dt=1/250 # dati giornalieri
dlnS=np.diff(np.log(SP['^GSPC'])) # differenze logaritmiche
sigma=np.sqrt(np.var(dlnS)/dt)
mu=np.mean(dlnS)/dt+0.5*sigma**2
mu
## np.float64(0.06638414723242678)
sigma
## np.float64(0.18940631332079627)

5.4 Simulazione del modello

Avendo calibrato il modello teorico sui dati empirici fino al 1995, possiamo ora simulare delle traiettorei (ne facciamo 1’000) per capire come il modello avrebbe previsto il futuro.

Vogliamo simulare tanti periodi quanta è la lunghezza dei dati out of sample e vogliamo che le simulazioni partano dall’ultimo valore conosciuto (cioè l’ultimo valore dell’insieme “in sample”).

N=1000 # simulazioni
T=length(SPF)*dt # anni da simulare

Sim=array(NA,dim=c(T/dt,N)) # creazione matrice vuota delle simulazioni
# Prendiamo il valore numerico dell'ultimo elemento dell'insieme "in sample"
Sim[1,]=as.numeric(tail(SP,1))

dW=array(rnorm(T/dt*N,0,sqrt(dt)),dim=c(T/dt,N))

for (i in 2:(T/dt)){
  Sim[i,]=Sim[i-1,]+
    Sim[i-1,]*mu*dt+
    Sim[i-1,]*sigma*dW[i-1,]
}

dim(Sim)
## [1] 7741 1000
N=1000 # simulazioni
T=len(SPF)*dt # anni da simulare

Sim=np.zeros((int(T/dt),N)) # matrice vuota
# Ultimo elemento dell'insieme "in sample"
Sim[0,:]=SP.tail(1)

dW=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))

for i in range(1,int(T/dt)):
  Sim[i,:]=Sim[i-1,:]+Sim[i-1,:]*mu*dt+Sim[i-1,:]*sigma*dW[i-1,:]

Sim.shape
## (7741, 1000)

5.5 Rappresentazione grafica dei risultati

Ora possiamo rappresentare su uno stesso grafico i valori dell’indice che si sono effettivamente verificati dal 1995 in poi insieme ai valori che sarebbero stati previsti dal modello. Questi valori si possono “riassumere” mediante alcuni indici statistici (per esempio la media, la mediana e alcuni quantili).

# Calcolo della media delle simulazioni (per riga)
M=apply(Sim,1,mean)

# Calcolo dei percentili 20, 50 (mediana) e 80
# L'operatore "t" traspone la matrice Q che avrebbe 3 righe (invece vogliamo 3 colonne)
Q=t(apply(Sim,1,quantile,probs=c(0.2,0.5,0.8)))
# Calcolo della media delle simulazioni (per riga)
M=np.mean(Sim,axis=1)

# Calcolo dei percentili 20, 50 (mediana) e 80
# Si traspone la matrice Q che avrebbe 3 righe (invece vogliamo 3 colonne)
Q=np.transpose(np.quantile(Sim, [0.3,0.5,0.8], axis=1))

I comandi per ottenere il grafico desiderato sono i seguenti. Per non avere una serie “in sample” troppo lunga, iniziamo a rappresentare i dati a partire dal 1980.

# Grafico dell'indice
plot(SP500[time(SP500)>='1980-01-01'],xlab='',ylab='')

# Grafico della media e dei quantili
matlines(time(SPF),cbind(M,Q),type='l',
        col=c('blue','red','orange','darkgreen'),
        lty=1)

# Legenda
legend('topleft',col=c('black','blue','red','orange','darkgreen'),
       legend=c('S&P500','Media simulazioni',
                '20 Percentile','Mediana','80 Percentile'),lty=1,
       bty='n')
grid()

plt.figure()
plt.plot(SP500[SP500.index>='1980-01-01'],color='black',label='SP500')
plt.plot(SPF.index,Q[:,0],color='red',label='10 Percentie')
plt.plot(SPF.index,Q[:,1],color='orange',label='Mediana')
plt.plot(SPF.index,Q[:,2],color='green',label='90 Percentile')
plt.plot(SPF.index,M,color='blue',label='Media simulazioni')
plt.legend()
plt.show()

6 Il processo di Poisson

La distribuzione di Poisson (detta “degli eventi rari”) permette di rappresentare eventi che accadono con bassa frequenza, ma che possono influenzare in modo significativo i mercati finanziari.

Chiamiamo \(d\Pi_t\) la variabile aleatoria di Poisson che può avere valore \(1\) con probabilità \(\phi dt\) oppure valore \(0\).

dt=1/250 # dati giornalieri
phi=10 # frequenza relativa dei salti (ogni anno)
dPi=rpois(250,phi*dt) # estraggo 250 valori dalla Poisson

plot(dPi)

dt=1/250 # dati giornalieri
phi=10 # frequenza relativa dei salti (ogni anno)
dPi=np.random.poisson(phi*dt,250) # estraggo 250 valori dalla Poisson

plt.figure()
plt.plot(dPi,linestyle='',marker='o')
plt.show()

6.1 Simulazione di un processo a salti

Immaginiamo di voler rappresentare una varaibile \(X_t\) che evolve a salti, in modo che la sua variazione \(dX_t\) sia pari al salto

\[dX_t=\gamma d\Pi_t\]

dove \(\gamma\) è l’ampiezza del salto. Proviamo a rappresentare graficamente una variabile che parte da \(X_0=0\) e che salta di un valore \(\gamma=1\) in media due volte all’anno \(\phi=2\).

Se vogliamo simulare \(T=5\) anni di valori giornalieri \(dt=1/250\) i comandi sono i seguenti.

X0=0 # valore iniziale
g=1 # salto
phi=2 # frequenza (annuale) del salto
T=5 # anni da simulare
dt=1/250 # dati giornalieri
N=100 # numero di simulazioni

X=array(NA,dim=c(T/dt,N))
X[1,]=X0
dPi=array(rpois(T/dt*N,phi*dt),dim=c(T/dt,N))

for (i in 2:(T/dt)){
  X[i,]=X[i-1,]+g*dPi[i-1,]
}

matplot(X,type='l',col='lightgray',lty=1)
matlines(apply(X,1,mean),lwd=2)

X0=0 # valore iniziale
g=1 # salto
phi=2 # frequenza (annuale) del salto
T=5 # anni da simulare
dt=1/250 # dati giornalieri
N=100 # numero di simulazioni

X=np.zeros((int(T/dt),N))
X[0,:]=X0
dPi=np.random.poisson(phi*dt,(int(T/dt),N))

for i in range(1,int(T/dt)):
  X[i,:]=X[i-1,:]+g*dPi[i-1,:]

plt.figure()
plt.plot(X,color='lightgray')
plt.plot(np.mean(X,axis=1),color='black')
plt.show()

6.2 Simulazione di processi a diffusione e salto

Immaginiamo, ora, di voler simulare un prcoesso a diffusione e salto nella forma seguente:

\[\frac{dX_t}{X_t}=(\mu -\gamma\phi)dt+\sigma dW_t + \gamma d\Pi_t\]

dove la deriva è stata compensata per avere, in media, un rendimento che è, comunque, uguale a \(\mu\).

Nel nostro caso immaginiamo un processo che cresca, in media, del 7%, con una volatilità di 0.2 e che, mediamente una volta all’anno, possa perdere il 20%. Infine, vogliamo simulare 5 anni di dati giornalieri per un processo che parte dal valore \(X_0=1\).

X0=1 # valore iniziale
g=-0.2 # salto
phi=1 # frequenza (annuale) del salto
mu=0.07 # tasso medio di crescita
sigma=0.2 # diffusione
T=5 # anni da simulare
dt=1/250 # dati giornalieri
N=100 # numero di simulazioni

X=array(NA,dim=c(T/dt,N))
X[1,]=X0
dPi=array(rpois(T/dt*N,phi*dt),dim=c(T/dt,N))
dW=array(rnorm(T/dt*N,0,sqrt(dt)),dim=c(T/dt,N))

for (i in 2:(T/dt)){
  X[i,]=X[i-1,]+
    X[i-1,]*(mu-phi*g)*dt+
    X[i-1,]*sigma*dW[i-1,]+
    X[i-1,]*g*dPi[i-1,]
}

matplot(X,type='l',col='lightgray',lty=1)
matlines(apply(X,1,mean),lwd=2)

X0=1 # valore iniziale
g=-0.2 # salto
phi=1 # frequenza (annuale) del salto
mu=0.07 # tasso medio di crescita
sigma=0.2 # diffusione
T=5 # anni da simulare
dt=1/250 # dati giornalieri
N=100 # numero di simulazioni

X=np.zeros((int(T/dt),N))
X[0,:]=X0
dPi=np.random.poisson(phi*dt,(int(T/dt),N))
dW=np.random.normal(0,np.sqrt(dt),(int(T/dt),N))

for i in range(1,int(T/dt)):
    X[i,:]=X[i-1,:]+X[i-1,:]*(mu-phi*g)*dt+X[i-1,:]*sigma*dW[i-1,:]+X[i-1,:]*g*dPi[i-1,:]

plt.figure()
plt.plot(X,color='lightgray')
plt.plot(np.mean(X,axis=1),color='black')
plt.show()

6.3 Rappresentazione dei rendimenti (logaritmici)

plot(diff(log(X[,1])),type='l')

plt.figure()
plt.plot(np.diff(np.log(X[:,0])))
plt.show()

6.4 Calibrazione del modello a diffusione e salto

Se tutti i parametri del modello a diffusione e salto sono costanti, allora possiamo utilizzare il metodo dei momenti applicandolo ai rendimenti logaritmici.

L’applicazione del lemma di Ito al caso di diffusione e salto porge:

\[d\ln X_{t}=\left(\mu-\phi\gamma-\frac{1}{2}\sigma^{2}\right)dt+\sigma dW_{t}+\ln\left(1+\gamma\right)d\Pi_{t}.\]

I momenti di questa distribuzione sono:

\[M1:=\mathbb{E}_{t}\left[d\ln X_{t}\right]=\left(\mu-\phi\gamma-\frac{1}{2}\sigma^{2}+\phi\ln\left(1+\gamma\right)\right)dt\] \[M2:=\mathbb{V}_{t}\left[d\ln X_{t}\right]=\left(\sigma^{2}+\left(\ln\left(1+\gamma\right)\right)^{2}\phi\right)dt\] \[M3:=\mathbb{E}_{t}\left[\left(d\ln X_{t}\right)^{3}\right]=\left(\ln\left(1+\gamma\right)\right)^{3}\phi dt\] \[M4:=\mathbb{E}_{t}\left[\left(d\ln X_{t}\right)^{4}\right]=\left(\ln\left(1+\gamma\right)\right)^{4}\phi dt\]

Se risolviamo questo sistema di 4 equazioni in 4 incognite, otteniamo

\[\gamma=e^{\frac{M4}{M3}}-1\] \[\phi=\frac{1}{dt}\frac{\left(M3\right)^{4}}{\left(M4\right)^{3}}\] \[\sigma^{2}=\frac{1}{dt}M2-\left(\ln\left(1+\gamma\right)\right)^{2}\phi\] \[\mu=\frac{1}{dt}M1+\phi\gamma+\frac{1}{2}\sigma^{2}-\phi\ln\left(1+\gamma\right)\]

Possiamo iniziare scaricando i dati finanziari dell’indice S&P500.

# Scarichiamo i dati S&P500
library(quantmod)
getSymbols('^GSPC',src='yahoo',return.class='zoo',
           from='1920-01-01')
## [1] "GSPC"
SP=GSPC$GSPC.Adjusted
dt=1/250

# Calcoliamo i rendimenti logaritmici
dlnS=diff(log(SP))

# Calcoliamo i momenti
M1=mean(dlnS)
M2=var(dlnS)
M3=mean(dlnS^3)
M4=mean(dlnS^4)

# Stima dei parametri
g=exp(M4/M3)-1
phi=1/dt*M3^4/M4^3
sigma=sqrt(M2/dt-phi*log(1+g)^2)
mu=M1/dt+phi*g+0.5*sigma^2-phi*log(1+g)
c(mu,sigma,g,phi)
## [1]  0.0783382030  0.1883699169 -0.4725109306  0.0006629961
# Scarichiamo i dati S&P500
S=yf.download('^GSPC',start="1920-01-01")
## <string>:2: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [*********************100%***********************]  1 of 1 completed
SP=S['Close']
dt=1/250

# Calcoliamo i rendimenti logaritmici
dlnS=np.diff(np.log(SP['^GSPC']))

# Calcoliamo i momenti
M1=np.mean(dlnS)
M2=np.var(dlnS)
M3=np.mean(dlnS**3)
M4=np.mean(dlnS**4)

# Stima dei parametri
g=np.exp(M4/M3)-1
phi=1/dt*M3**4/M4**3
sigma=np.sqrt(M2/dt-phi*np.log(1+g)**2)
mu=M1/dt+phi*g+0.5*sigma**2-phi*np.log(1+g)
mu,sigma,g,phi
## (np.float64(0.07837491799957795), np.float64(0.18836260020464987), np.float64(-0.4725121905369317), np.float64(0.0006629591660505729))

6.5 Confronto con una distribuzione normale

Sappiamo che, se \(x\) è distribuita normalmente (con media \(\mu\) e s.q.m. \(\sigma\)), vale che

\[\mathbb{P}\left\{\mu-2\sigma<x<\mu+2\sigma \right\}=0.9545\]

Vediamo nel nostro caso come sono distribuiti gli eventi più alti e più bassi della soglia di 2 volte lo scarto quadratico medio.

Proviamo, dapprima, a rappresentare graficamente i rendimenti logaritmici con la loro media e le due soglie della media aumentata e diminuita dello scarto quadratico medio.

plot(dlnS,type='l')
abline(h=c(M1-2*sqrt(M2),M1,M1+2*sqrt(M2)),
        col=c('red','blue','red'),lty=1,lwd=2)

plt.figure()
plt.plot(dlnS,color='black')
plt.axhline(M1-2*np.sqrt(M2),color='red',linewidth=2)
plt.axhline(M1,color='blue',linewidth=2)
plt.axhline(M1+2*np.sqrt(M2),color='red',linewidth=2)
plt.show()

Vediamo ora la frequenza relativa dei dati dentro le due soglie.

# Isolamento dei salti
salti=dlnS[dlnS<(M1-2*sqrt(M2)) | dlnS>(M1+2*sqrt(M2))]

# Calcolo della media (annuale) dei salti
mean(salti)/dt
## [1] -1.029766
# Calcolo della frequenza relativa dei salti
length(salti)/length(dlnS)
## [1] 0.04654856
# Isolamento dei salti
salti=dlnS[(dlnS<M1-2*np.sqrt(M2)) | (dlnS>M1+2*np.sqrt(M2))]

# Calcolo della media (annuale) dei salti
np.mean(salti)/dt
## np.float64(-1.0297662635329843)
# Calcolo della frequenza relativa dei salti
len(salti)/len(dlnS)
## 0.04654666883857306

7 Processi a ritorno sulla media

7.1 I tassi di interesse

Consideriamo i dati del tasso a 3 mesi dei T-Bill americani (serie FRED ‘DTB3’). Vediaom che questa variabile non è mediamente sempre crecente come lo S&P500, bensì tende a oscillare intorno a un valore medio.

getSymbols('DTB3',src='FRED',return.class='zoo')
## [1] "DTB3"
tasso=na.omit(DTB3)/100 #per non avere dati in percentuale
plot(tasso)
abline(h=mean(tasso),col='blue',lty=2)

tasso=fred.get_series('DTB3')
tasso=tasso[~np.isnan(tasso)]/100

plt.figure()
plt.plot(tasso)
plt.axhline(np.mean(tasso),color='red',linestyle='--')
plt.show()

7.2 Il processo di Vasicek - Stima dei parametri

Cerchiamo di adattare a questi dati il modello di Vasicek nella forma seguente:

\[dr_t=a\left( b-r_t\right) dt + \sigma dW_t\]

Il processo ha volatilità istantanea costante:

\[\mathbb{V}_t\left[dr_t\right]=\sigma^2 dt\]

La stima del parametro \(\sigma\) appare quindi molto semplice.

dr=diff(tasso)
dt=1/250
sigma_r=sd(dr)/sqrt(dt)
sigma_r
## [1] 0.01331758
dr=np.diff(tasso)
dt=1/250
sigma_r=np.std(dr)/np.sqrt(dt)
sigma_r
## np.float64(0.013317209600394241)

La restante parte del processo (quella che contiene i parametri \(a\) e \(b\)) è meno banale da stimare. Dobbiamo ricorrere a una stima dei minimi quadrati (OLS) sul modello dicretizzato:

\[r_t = r_{t-1} + a\left( b-r_{t-1}\right) dt + \epsilon_t\]

ovvero

\[r_t = a\cdot b\cdot dt+ \left( 1 -a\cdot dt \right)\cdot r_{t-1}+ \epsilon_t\] La stima dei parametri avverrà sull’equazione seguente

\[r_t = \beta_0 + \beta_1 r_{t-1} + \epsilon_t\]

Una volta effettuata la stima dei parametri \(\beta_0\) e \(\beta_1\), si potranno ottenere i valori dei parametri originali \(a\) e \(b\):

\[a=\frac{1-\beta_{1}}{dt}\]

\[b=\frac{\beta_{0}}{1-\beta_{1}}\]

Vediamo come procedere sui software.

Nel caso di Python è necessario scaricare un pacchetto statistico per effettuare la regressione dei minimi quadrati.

#Usiamo il comando "lm" (linear model)
regr=lm(tail(tasso,-1) ~ head(tasso,-1))
summary(regr)
## 
## Call:
## lm(formula = tail(tasso, -1) ~ head(tasso, -1))
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0126531 -0.0001998 -0.0000092  0.0001967  0.0134407 
## 
## Coefficients:
##                  Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)     1.829e-05  1.067e-05    1.714   0.0865 .  
## head(tasso, -1) 9.996e-01  2.054e-04 4867.040   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0008422 on 17927 degrees of freedom
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9992 
## F-statistic: 2.369e+07 on 1 and 17927 DF,  p-value: < 2.2e-16
beta0=regr$coefficients[[1]]
beta1=regr$coefficients[[2]]
a=(1-beta1)/dt
a
## [1] 0.1005641
b=beta0/(a*dt)
b
## [1] 0.04546422
# b deve essere prossimo alla media del processo
mean(tasso)
## [1] 0.04195602
import statsmodels.api as stat

#Delle variabili in gioco prendiamo solo i valori (senza le date)
model=stat.OLS(tasso.values[1:],stat.add_constant(tasso.values[:-1]))
results=model.fit()
results.summary()
OLS Regression Results
Dep. Variable: y R-squared: 0.999
Model: OLS Adj. R-squared: 0.999
Method: Least Squares F-statistic: 2.369e+07
Date: lun, 06 ott 2025 Prob (F-statistic): 0.00
Time: 21:19:58 Log-Likelihood: 1.0149e+05
No. Observations: 17929 AIC: -2.030e+05
Df Residuals: 17927 BIC: -2.030e+05
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 1.829e-05 1.07e-05 1.714 0.087 -2.62e-06 3.92e-05
x1 0.9996 0.000 4867.040 0.000 0.999 1.000
Omnibus: 5575.491 Durbin-Watson: 1.717
Prob(Omnibus): 0.000 Jarque-Bera (JB): 963674.652
Skew: 0.249 Prob(JB): 0.00
Kurtosis: 38.913 Cond. No. 32.7


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified. 
beta0=results.params[0]
beta1=results.params[1]
a=(1-beta1)/dt
a
## np.float64(0.1005640884186032)
b=beta0/(a*dt)
b
## np.float64(0.0454642208362537)
# b deve essere prossimo alla media del processo
np.mean(tasso)
## np.float64(0.041956017847183495)

7.3 Il processo di Vasicek - Simulazioni

Adesso possiamo simulare il processo di Vasice con i valori stimati per vedere se il processo si adegua bene ai dati empirici.

La simulazione segue le stesse regole già viste in precedenza. In R i comandi sono i seguenti.

T=length(tasso)*dt
N=100

dW=array(rnorm(T/dt*N,0,sqrt(dt)),dim=c(T/dt,N))
r=array(NA,dim=c(T/dt,N))
r[1,]=tasso[1]

for (i in 2:(T/dt)){
  r[i,]=r[i-1,]+a*(b-r[i-1,])*dt+sigma_r*dW[i-1,]
}

matplot(time(tasso),r,type='l',col='lightgray')
lines(tasso,lwd=2)

T=len(tasso)*dt
N=100

dW=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))
r=np.zeros((int(T/dt),N))
r[0,:]=tasso[0]
## <string>:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
# ciclo
for i in range(1,int(T/dt)):
  r[i,:]=r[i-1,:]+a*(b-r[i-1,:])*dt+sigma_r*dW[i-1,:]

plt.figure()
plt.plot(tasso.index,r,color='lightgray')
plt.plot(tasso,color='black',linewidth=2)
plt.show()

Vediamo così che il processo di Vasicek non garantisce che i tassi di interesse rimangano sempre positivi e, anzi, quando sono molto bassi, ha una probabilità non trascurabile che vadano nel dominio negativo.

7.4 Il processo CIR - Stima dei parametri

Il modello Cox-Ingersoll-Ross restituisce unicamente tassi di interesse positivi e la sua forma algebrica è la seguente:

\[dr_t = a\left( b-r_t\right) dt+\sigma_r \sqrt{r_t}dW_t\]

In questo caso anche la stima di \(\sigma_r\) non può essere ottenuta mediante i metodi visti in precedenza e, inoltre, essendo l’errore eteroschedastico, non si può stimare la dinamica mediante i minimi qudarati ordinari.

Tuttavia, esiste una trasformazione che rende questo processo omoschedastico. Se, infatti, si applica il lemma di Ito alla trasformazione \(y_t=2\sqrt{r_t}\), si ottiene il seguente processo:

\[dy_t=\left(-\frac{a}{2}y_t+\left(2ab-\frac{1}{2}\sigma_{r}^{2}\right)\frac{1}{y_t}\right)dt+\sigma_{r}dW_t\]

Da questa equazione si può stimare immediatamente il termine \(\sigma_r\) osservando che vale

\[\mathbb{V}_t\left[dy_t\right]=\sigma_r^2 dt\]

L’equazione da stimare, quindi, è:

\[y_{t+1}=\underset{\beta_{0}}{\underbrace{\left(1-\frac{a\cdot dt}{2}\right)}}y_t+\underset{\beta_{1}}{\underbrace{\left(2\cdot a\cdot b-\frac{1}{2}\sigma_{r}^{2}\right)dt}}\frac{1}{y_t}+\epsilon_t\]

Dal momento che il termine \(y_t\) copmare al denominatore, bisogna assicurarsi che esso non sia mai nullo. Inoltre, poiché \(y_t\) contiene la radice quadrata del tasso, bisogna anche verificare che il tasso non diventi mai negativo.

I casi di tassi nulli o negativi sul T-Bill trimestrale sono veramente ridotti.

tasso[which(tasso<=0)]
## 2008-12-10 2008-12-18 2008-12-24 2011-09-22 2011-10-05 2011-12-15 2011-12-16 
##      0e+00      0e+00      0e+00      0e+00      0e+00      0e+00      0e+00 
## 2013-09-26 2015-09-18 2015-09-22 2015-09-25 2015-09-30 2015-10-01 2015-10-02 
##      0e+00     -1e-04     -1e-04     -1e-04     -1e-04     -2e-04      0e+00 
## 2015-10-06 2015-10-07 2015-10-08 2015-10-14 2015-10-22 2020-03-25 2020-03-26 
##      0e+00      0e+00     -1e-04      0e+00      0e+00     -4e-04     -5e-04
tasso[tasso<=0]
## 2008-12-10    0.0000
## 2008-12-18    0.0000
## 2008-12-24    0.0000
## 2011-09-22    0.0000
## 2011-10-05    0.0000
## 2011-12-15    0.0000
## 2011-12-16    0.0000
## 2013-09-26    0.0000
## 2015-09-18   -0.0001
## 2015-09-22   -0.0001
## 2015-09-25   -0.0001
## 2015-09-30   -0.0001
## 2015-10-01   -0.0002
## 2015-10-02    0.0000
## 2015-10-06    0.0000
## 2015-10-07    0.0000
## 2015-10-08   -0.0001
## 2015-10-14    0.0000
## 2015-10-22    0.0000
## 2020-03-25   -0.0004
## 2020-03-26   -0.0005
## dtype: float64

Il primo passaggio, quindi, è quello di togliere dal campione i dati negativi o nulli per procedere, poi, con la stima dei parametri.

tasso2=tasso[which(tasso>0)]
y=2*sqrt(tasso2)

sigma_r=sd(diff(y))/sqrt(dt)

Y=tail(y,-1) #definizione var. dipendente
X1=head(y,-1) #definizione primo regressore
X2=head(y,-1)^(-1) #definizione secondo regressore 
regr=lm(Y ~ X1 + X2 -1) #"-1" indica che non c'è la costante di regressione
summary(regr)
## 
## Call:
## lm(formula = Y ~ X1 + X2 - 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.148742 -0.001205 -0.000020  0.001244  0.100869 
## 
## Coefficients:
##     Estimate Std. Error  t value Pr(>|t|)    
## X1 9.998e-01  7.950e-05 12576.91  < 2e-16 ***
## X2 1.981e-05  3.558e-06     5.57 2.59e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004203 on 17906 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 8.515e+07 on 2 and 17906 DF,  p-value: < 2.2e-16
tasso2=tasso[tasso>0]
y=2*np.sqrt(tasso2)

sigma_r=np.std(np.diff(y))/np.sqrt(dt)

Y=y.values[1:]
X1=y.values[:-1]
X2=y.values[:-1]**(-1)

model=stat.OLS(Y,np.column_stack((X1, X2)))
results=model.fit()
results.summary()
OLS Regression Results
Dep. Variable: y R-squared (uncentered): 1.000
Model: OLS Adj. R-squared (uncentered): 1.000
Method: Least Squares F-statistic: 8.515e+07
Date: lun, 06 ott 2025 Prob (F-statistic): 0.00
Time: 21:20:04 Log-Likelihood: 72581.
No. Observations: 17908 AIC: -1.452e+05
Df Residuals: 17906 BIC: -1.451e+05
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
x1 0.9998 7.95e-05 1.26e+04 0.000 1.000 1.000
x2 1.981e-05 3.56e-06 5.570 0.000 1.28e-05 2.68e-05
Omnibus: 11583.395 Durbin-Watson: 1.841
Prob(Omnibus): 0.000 Jarque-Bera (JB): 11978392.236
Skew: -1.667 Prob(JB): 0.00
Kurtosis: 129.658 Cond. No. 23.2


Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.

I valori dei parametri \(a\) e \(b\) sono ottenuti dai valori stimati di \(\beta_0\) e \(\beta_1\).

\[a=2\frac{1-\beta_{0}}{dt}\]

\[b=\frac{\beta_{1}+\frac{1}{2}\sigma_{r}^{2}dt}{4\left(1-\beta_{0}\right)}\]

beta0=regr$coefficients[[1]]
beta1=regr$coefficients[[2]]
a=2*(1-beta0)/dt
a
## [1] 0.07688597
b=(beta1+0.5*sigma_r^2*dt)/(4*(1-beta0))
b
## [1] 0.04659997
beta0=results.params[0]
beta1=results.params[1]
a=2*(1-beta0)/dt
a
## np.float64(0.07688596817900795)
b=(beta1+0.5*sigma_r**2*dt)/(4*(1-beta0))
b
## np.float64(0.04659916741757053)

Il valore a cui il processo tende (\(b\)) è di nuovo simile al valore medio empirico.

In Python vediamo come ottenere gli stessi risultati.

7.5 Il processo CIR - Simulazioni

Per simulare il processo CIR i comandi sono identici a quelli visti per il processo di Vasicek, con l’unica modifica che riguarda la volatilità del processo.

T=length(tasso)*dt
N=100

dW=array(rnorm(T/dt*N,0,sqrt(dt)),dim=c(T/dt,N))
r=array(NA,dim=c(T/dt,N))
r[1,]=tasso[1]

for (i in 2:(T/dt)){
  r[i,]=r[i-1,]+a*(b-r[i-1,])*dt+sigma_r*sqrt(r[i-1,])*dW[i-1,]
}
## Warning in sqrt(r[i - 1, ]): Si è prodotto un NaN
## Warning in sqrt(r[i - 1, ]): Si è prodotto un NaN
## Warning in sqrt(r[i - 1, ]): Si è prodotto un NaN
matplot(time(tasso),r,type='l',col='lightgray')
lines(tasso,lwd=2)

T=len(tasso)*dt
N=100

dW=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))
r=np.zeros((int(T/dt),N))
r[0,:]=tasso[0]
## <string>:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
# ciclo
for i in range(1,int(T/dt)):
  r[i,:]=r[i-1,:]+a*(b-r[i-1,:])*dt+sigma_r*np.sqrt(r[i-1,:])*dW[i-1,:]
## <string>:4: RuntimeWarning: invalid value encountered in sqrt
plt.figure()
plt.plot(tasso.index,r,color='lightgray')
plt.plot(tasso,color='black',linewidth=2)
plt.show()

Possiamo verificare come si muova il tasso “vero” rispetto alla media e ai quantili dei processi simulati.

Q=t(apply(r,1,quantile,probs=c(0.1,0.5,0.9),na.rm=T))
#l'opzione na.rm=T (NA remove = TRUE) è necessaria per ignorare i NAN che si possono produrre nelle simulazioni del processo CIR

matplot(time(tasso),
        cbind(Q,as.numeric(tasso)),
        type='l',col=c('darkgreen','blue','red','black'),lty=1,
        xlab='',ylab='')
legend('topleft',col=c('darkgreen','blue','red','black'),lty=1,
       legend=c('Percentile 0.1','Mediana','Percentile 0.9','Dati storici'),
       bty='n')
grid()

Q=np.nanquantile(r, [0.1,0.5,0.9], axis=1)
#nanquantile è necessario per ignorare i NAN che si possono produrre nelle simulazioni del processo CIR
Q=np.transpose(Q)

plt.figure()
plt.plot(tasso,color='black',label='Dati storici')
plt.plot(tasso.index,Q[:,0],color='red',label='Quantile 0.1')
plt.plot(tasso.index,Q[:,1],color='blue',label='Mediana')
plt.plot(tasso.index,Q[:,2],color='darkgreen',label='Quantile 0.9')
plt.legend(frameon=False)
plt.show()

8 Il portafoglio ottimo

8.1 I dati

Prendiamo tre titoli rischiosi su cui calcoliamo i portafogli ottimi secondo diversi criteri. I titoli sono Amazon, Coca Cola e Google, i cui codici su yahoo finance sono, rispettivamente, AMZN, COKE e GOOG.

Per produrre un grafico con prezzi che siano confrontabili, poniamo pari a 1 il prezzo iniziale di ogni serie, cioè dividiamo i prezzi di ogni colonna per il primo prezzo di ogni colonna.

getSymbols(c('AMZN','COKE','GOOG'),
           src='yahoo',return.class='zoo',
           from='2019-01-01',to='2019-12-31')
## [1] "AMZN" "COKE" "GOOG"
# Si fondono tra loro le serie storiche
S=merge(AMZN=AMZN$AMZN.Adjusted,
        COKE=COKE$COKE.Adjusted,
        GOOG=GOOG$GOOG.Adjusted,all=FALSE)

matplot(time(S),
        S%*%diag(as.numeric(S[1,]^(-1))),type='l',
        col=c('orange','red','blue'),lty=1,
        xlab='',ylab='prezzi')
legend('topleft',legend=c('Amazon','Coke','Google'),
       col=c('orange','red','blue'),lty=1,bty='n')
grid()

S=yf.download(['AMZN','COKE','GOOG'],start="2019-01-01",end='2019-12-31')
## <string>:1: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [                       0%                       ][**********************67%*******                ]  2 of 3 completed[*********************100%***********************]  3 of 3 completed
S=S['Close']

plt.figure()
plt.plot(S['AMZN']/S['AMZN'][1],label='Amazon',color='orange')
plt.plot(S['COKE']/S['COKE'][1],label='Coke',color='red')
plt.plot(S['GOOG']/S['GOOG'][1],label='Google',color='blue')
plt.legend()
plt.show()

8.2 Stima dei parametri

Il modello teorico con i prezzi che seguono moti browniani geometrici con tante fonti di rischio quanti sono i titoli porta a rendimenti logaritmici che si scrivono nel modo seguente:

\[d\ln S_{t}=\left(\mu-\frac{1}{2}\text{diag}\left(\Sigma^{\intercal}\Sigma\right)\right)dt+\Sigma^{\intercal}dW_{t}.\]

Calcolando il primo e il secondo momento si ottiene:

\[\begin{cases} \mathbb{E}_{t}\left[d\ln S_{t}\right]=\left(\mu-\frac{1}{2}\text{diag}\left(\Sigma^{\intercal}\Sigma\right)\right)dt,\\ \mathbb{V}_{t}\left[d\ln S_{t}\right]=\Sigma^{\intercal}\Sigma dt. \end{cases}\]

Dove \(\Sigma^{\intercal}\Sigma\) è la Matrice di Varianze e Covarianze (MVC). Risolvendo questo sistema rispetto ai rendimenti e a MVC si ha:

\[\begin{cases} \mu=\frac{1}{dt}\mathbb{E}_{t}\left[d\ln S_{t}\right]+\frac{1}{2}\text{diag}\left(\Sigma^{\intercal}\Sigma\right),\\ \Sigma^{\intercal}\Sigma=\frac{1}{dt}\mathbb{V}_{t}\left[d\ln S_{t}\right]. \end{cases}\]

La stima dei parametri avviene con i seguenti comandi.

# Calcolo dei rendimenti logaritmici
dlnS=apply(log(S),2,diff)
dt=1/250

# Calcolo MVC
MVC=var(dlnS)/dt
MVC
##            AMZN       COKE       GOOG
## AMZN 0.05191702 0.01216833 0.03363764
## COKE 0.01216833 0.13822214 0.01003392
## GOOG 0.03363764 0.01003392 0.05774797
# Calcolo dei rendimenti medi
mu=apply(dlnS,2,mean)/dt+0.5*diag(MVC)
mu
##      AMZN      COKE      GOOG 
## 0.2082443 0.5457040 0.2738290
# Calcolo dei rendimenti logaritmici
dlnS=np.diff(np.log(S),axis=0)
dt=1/250

# Calcolo MVC
MVC=np.cov(dlnS,rowvar=False)/dt
MVC
## array([[0.05191702, 0.01216827, 0.03363764],
##        [0.01216827, 0.13822202, 0.0100339 ],
##        [0.03363764, 0.0100339 , 0.05774799]])
# Calcolo dei rendimenti medi
mu=np.mean(dlnS,axis=0)/dt+0.5*np.diagonal(MVC)
mu
## array([0.20824434, 0.54570464, 0.27382889])

8.3 Simulazione di più processi stocastici

Abbiamo stimati i parametri \(\mu\) e \(\Sigma^\intercal\Sigma\) per un insieme di processi stocastici. Per poter simulare i processi stocastici, tuttavia, occorre avere a disposizione la matrice \(\Sigma\). Essa non può essere univocamente ottenuta dalla matrice di varianze e covarianze.

Vediamo con un esempio semplice di due titoli influenzati da due fonti di rischio. Il modello si può scrivere nel modo seguente:

\[\frac{dS_{1}\left(t\right)}{S_{1}\left(t\right)}=\mu_{1}dt+adW_{1}\left(t\right)+bdW_{2}\left(t\right)\] \[\frac{dS_{2}\left(t\right)}{S_{2}\left(t\right)}=\mu_{2}dt+mdW_{1}\left(t\right)+ndW_{2}\left(t\right)\]

La matrice di varianze e covarianze, in questo caso, è data da:

\[MVC=\left[\begin{array}{cc} a^{2}+b^{2} & am+bn\\ am+bn & m^{2}+n^{2} \end{array}\right]\]

Se si stima la MVC sul mercato si ottengono i seguenti valori:

\[MVC=\left[\begin{array}{cc} \sigma_{1}^{2} & \rho\sigma_{1}\sigma_{2}\\ \rho\sigma_{1}\sigma_{2} & \sigma_{2}^{2} \end{array}\right]\]

dove \(\rho\) è l’indice di correlazione (di Pearson-Bravais). Per fare in modo che il modello teorico coincida con i dati empirici, dobbiamo risolvere il seguente sistema:

\[\begin{cases} \sigma_{1}^{2}=a^{2}+b^{2}\\ \rho\sigma_{1}\sigma_{2}=am+bn\\ \sigma_{2}^{2}=m^{2}+n^{2} \end{cases}\]

che ha 3 equazioni e 4 incognite (\(a,b,m,n\)). Il problema, dunque, ha infinite soluzioni possibili.

Una soluzione possibile è quella proposta da Cholesky che, in termini finanziari, consiste nell’ipotesi che:

  • il primo titolo dipenda solo dalla prima fonte di rischio;
  • il secondo titolo dipenda dalle prime due fonti di rischio;
  • il terzo titolo dipenda dalle prime tre fonti di rischio;
  • e così via, fino all’ultimo titolo che dipende da tutte le fonti di rischio.

Nel caso del nostro esempio questo coincide con l’ipotesi che si abbia \(b=0\). In questo caso il sistema ha una sola soluzione:

\[\begin{cases} \sigma_{1}=a\\ \rho\sigma_{2}=m\\ \sigma_{2}\sqrt{\left(1-\rho^{2}\right)}=n \end{cases}\]

e, quindi, una volta stimati \(\sigma_1\), \(\sigma_2\) e \(\rho\), si può rappresentare il mercato come:

\[\frac{dS_{1}\left(t\right)}{S_{1}\left(t\right)}=\mu_{1}dt+\sigma_{1}dW_{1}\left(t\right)\]

\[\frac{dS_{2}\left(t\right)}{S_{2}\left(t\right)}=\mu_{2}dt+\rho\sigma_{2}dW_{1}\left(t\right)+\sigma_{2}\sqrt{\left(1-\rho^{2}\right)}dW_{2}\left(t\right)\]

dove la matrice dei termini di diffusione è data da:

\[\Sigma^{\intercal}=\left[\begin{array}{cc} \sigma_{1} & 0\\ \rho\sigma_{2} & \sigma_{2}\sqrt{\left(1-\rho^{2}\right)} \end{array}\right]\]

Questo modo di scomporre la matrice di varianze e covarianze nel prodotto di due matrici quadrate sub-diagonali è stato adottato, per la prima volta, da Cholesky. In molti software, quindi, si usa il comando “cholesky” per effettuare questa scomposizione.

#R restituisce una matrice con zeri nella zona sud-ovest
chol(MVC)
##           AMZN       COKE        GOOG
## AMZN 0.2278531 0.05340426 0.147628655
## COKE 0.0000000 0.36792679 0.005843333
## GOOG 0.0000000 0.00000000 0.189524678
SigmaT=t(MVC)
#Python restituisce una matrice con zeri nella zona nord-est
SigmaT=np.linalg.cholesky(MVC)
SigmaT
## array([[0.22785307, 0.        , 0.        ],
##        [0.05340403, 0.36792667, 0.        ],
##        [0.14762863, 0.00584337, 0.18952475]])

Adesso possiamo simulare contemporaneamente le traiettorie per i tre titoli.

T=5 #anni di simulare
dt=1/250 #frequenza dei dati
N=100 #numero di percorsi da simulare

# Si generano le tre variabili aleatorie
dW1=array(rnorm(T/dt*N,0,sqrt(dt)),
          dim=c(T/dt,N))
dW2=array(rnorm(T/dt*N,0,sqrt(dt)),
          dim=c(T/dt,N))
dW3=array(rnorm(T/dt*N,0,sqrt(dt)),
          dim=c(T/dt,N))

# Si ipotizza di partire dallo stesso valore 1
Sim1=Sim2=Sim3=array(NA,dim=c(T/dt,N))
Sim1[1,]=1
Sim2[1,]=1
Sim3[1,]=1

for (i in 2:(T/dt)){
  Sim1[i,]=Sim1[i-1,]+
    Sim1[i-1,]*mu[1]*dt+
    Sim1[i-1,]*SigmaT[1,]%*%rbind(dW1[i,],dW2[i,],dW3[i,])
  Sim2[i,]=Sim2[i-1,]+
    Sim2[i-1,]*mu[2]*dt+
    Sim2[i-1,]*SigmaT[2,]%*%rbind(dW1[i,],dW2[i,],dW3[i,])
  Sim3[i,]=Sim3[i-1,]+
    Sim3[i-1,]*mu[3]*dt+
    Sim3[i-1,]*SigmaT[3,]%*%rbind(dW1[i,],dW2[i,],dW3[i,])
}

matplot(Sim1,type='l',lty=1,col=rgb(0.8,0.3,0,0.2),
        ylim=c(min(Sim1,Sim2,Sim3),max(Sim1,Sim2,Sim3)))
matlines(Sim2,type='l',lty=1,col=rgb(1,0,0,0.2))
matlines(Sim3,type='l',lty=1,col=rgb(0,0,1,0.2))
grid()

T=5 #anni di simulare
dt=1/250 #frequenza dei dati
N=100 #numero di percorsi da simulare

# Si generano le tre variabili aleatorie
dW1=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))
dW2=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))
dW3=np.random.normal(0, np.sqrt(dt),(int(T/dt),N))

# Si ipotizza di partire dallo stesso valore 1
Sim1=np.zeros((int(T/dt),N))
Sim2=np.zeros((int(T/dt),N))
Sim3=np.zeros((int(T/dt),N))
Sim1[0,:]=1
Sim2[0,:]=1
Sim3[0,:]=1

for i in range(1,int(T/dt)):
  Sim1[i,:]=Sim1[i-1,:]+Sim1[i-1,:]*mu[0]*dt+Sim1[i-1,:]*np.matmul(SigmaT[0,:],np.row_stack((dW1[i,:],dW2[i,:],dW3[i,:])))
  Sim2[i,:]=Sim2[i-1,:]+Sim2[i-1,:]*mu[1]*dt+Sim2[i-1,:]*np.matmul(SigmaT[1,:],np.row_stack((dW1[i,:],dW2[i,:],dW3[i,:])))
  Sim3[i,:]=Sim3[i-1,:]+Sim3[i-1,:]*mu[2]*dt+Sim3[i-1,:]*np.matmul(SigmaT[2,:],np.row_stack((dW1[i,:],dW2[i,:],dW3[i,:])))
## <string>:3: DeprecationWarning: `row_stack` alias is deprecated. Use `np.vstack` directly.
## <string>:4: DeprecationWarning: `row_stack` alias is deprecated. Use `np.vstack` directly.
## <string>:5: DeprecationWarning: `row_stack` alias is deprecated. Use `np.vstack` directly.
plt.figure()
plt.plot(Sim1,color='orange',alpha=0.2)
plt.plot(Sim2,color='red',alpha=0.2)
plt.plot(Sim3,color='blue',alpha=0.2)
plt.grid()
plt.show()

8.4 Il tasso privo di rischio

Nel calcolo del portafoglio compare sempre il tasso privo di rischio. I modelli più semplici lo danno per costante. Prendiamo, quindi, il tasso T-Bill a 3 mesi (nello stesso periodo considerato) e ne calcoliamo la media.

# Scarichiamo la serie T-Bill 3m
getSymbols('DTB3',src='FRED',return.class='zoo')
## [1] "DTB3"
DTB3=na.omit(DTB3)/100

# Prendiamo il periodo desiderato
tasso=window(DTB3,start='2019-01-01',end='2019-12-31')
plot(tasso,xlab='',ylab='')
r=mean(tasso)
r
## [1] 0.0206248
abline(h=r,lty=2,col='blue')
grid()

# Scarichiamo la serie T-Bill 3m
tasso=fred.get_series('DTB3',observation_start='2019-01-01',observation_end='2019-12-31')
tasso=tasso/100
r=np.mean(tasso)

# Grafico
plt.figure()
plt.plot(tasso,color='black')
plt.axhline(y=r,color='blue',linestyle='-')
plt.show()

8.5 Il portafoglio di massima crescita

Il portafoglio di titoli rischiosi che massimizza il rendimento logaritmico atteso (Growth Optimal Portfolio) è

\[\omega_{GOP}=\left(\Sigma^{\intercal}\Sigma\right)^{-1}\left(\mu-r\mathbb{1}\right)\]

I calcoli sono i seguenti.

omega_GOP=solve(MVC)%*%(mu-r)
omega_GOP
##           [,1]
## AMZN 0.5569637
## COKE 3.4991687
## GOOG 3.4522226
omega_GOP=np.matmul(np.linalg.inv(MVC),mu-r)
omega_GOP
## array([0.55696935, 3.49917685, 3.45221699])

8.6 L’ottimizzazione dell’utilità marginale attesa

Questo portafoglio è proporzionale al GOP e il termine di proporzionalità è dato dall’inverso dell’indice di avversione relativa al rischio. Se questo indice viene chiamato \(\delta\), il potafoglio ha la forma:

\[\omega_U=\frac{1}{\delta}\omega_{GOP}.\]

Se consideriamo \(\delta=2.5\) il portafoglio si calcola come segue.

delta=2.5
omega_U=1/delta*omega_GOP
omega_U
##           [,1]
## AMZN 0.2227855
## COKE 1.3996675
## GOOG 1.3808890
delta=2.5
omega_U=1/delta*omega_GOP
omega_U
## array([0.22278774, 1.39967074, 1.3808868 ])

8.7 Il portafoglio sicuro

Per questo portafoglio, che minimizza la probabilità che il rendimento logaritmico del portafoglio scenda sotto una certa soglia \(\mu_E\), occorre definire il quadrato del prezzo di mercato del rischio:

\[\xi^\intercal\xi:=\left(\mu-r\mathbb{1}\right)^\intercal\left(\Sigma^\intercal\Sigma\right)^{-1}\left(\mu-r\mathbb{1}\right)\] In questo caso il portafoglio ottimo è:

\[\omega_P=\sqrt{2\frac{\mu_E-r}{\xi^\intercal\xi}}\omega_{GOP}\]

Ecco i calcoli.

xi2=as.numeric((mu-r)%*%solve(MVC)%*%(mu-r))
mu_E=0.1
omega_P=sqrt(2*(mu_E-r)/xi2)*omega_GOP
omega_P
##           [,1]
## AMZN 0.1322426
## COKE 0.8308248
## GOOG 0.8196781
xi2=np.matmul(mu-r,np.matmul(np.linalg.inv(MVC),mu-r))
mu_E=0.1
omega_P=np.sqrt(2*(mu_E-r)/xi2)*omega_GOP
omega_P
## array([0.13224382, 0.83082582, 0.81967592])

8.8 Il portafoglio media-varianza

Il portafoglio che minimizza la varianza del rendimento logaritmico, dato un rendimento atteso \(\mu_E\) che si vuole raggiungere (in media), è dato da

\[\omega_{MV}=\frac{\mu_E-r}{\xi^\intercal\xi}\omega_{GOP}\]

I comandi sono:

omega_MV=(mu_E-r)/xi2*omega_GOP
omega_MV
##            [,1]
## AMZN 0.01569951
## COKE 0.09863340
## GOOG 0.09731009
omega_MV=(mu_E-r)/xi2*omega_GOP
omega_MV
## array([0.01569963, 0.09863342, 0.09730973])

8.9 Confronto tra portafogli

Si può creare un grafico per confrontare tra loro i portafogli. Da principio creiamo una matrice con tutti i portafogli, aggiungendo, in fondo, anche la quantità di titolo privo di rischio.

# Allineamento di tutti i portafogli
omega=cbind(omega_GOP,omega_U,omega_P,omega_MV)
omega
##           [,1]      [,2]      [,3]       [,4]
## AMZN 0.5569637 0.2227855 0.1322426 0.01569951
## COKE 3.4991687 1.3996675 0.8308248 0.09863340
## GOOG 3.4522226 1.3808890 0.8196781 0.09731009
# Calcolo del titolo privo di rischio
TBILL=1-apply(omega,2,sum)

# Allieamento dei portafogli
omega=rbind(omega,TBILL)
omega
##             [,1]       [,2]       [,3]       [,4]
## AMZN   0.5569637  0.2227855  0.1322426 0.01569951
## COKE   3.4991687  1.3996675  0.8308248 0.09863340
## GOOG   3.4522226  1.3808890  0.8196781 0.09731009
## TBILL -6.5083550 -2.0033420 -0.7827455 0.78835700
# Si sovrappongono i portafoglio
omega=np.vstack((omega_GOP,omega_U,omega_P,omega_MV))
omega=omega.transpose()
omega
## array([[0.55696935, 0.22278774, 0.13224382, 0.01569963],
##        [3.49917685, 1.39967074, 0.83082582, 0.09863342],
##        [3.45221699, 1.3808868 , 0.81967592, 0.09730973]])
# Si aggiunge il titolo privo di rischio
TBILL=1-np.sum(omega,axis=0)
omega=np.vstack([omega, TBILL])
omega
## array([[ 0.55696935,  0.22278774,  0.13224382,  0.01569963],
##        [ 3.49917685,  1.39967074,  0.83082582,  0.09863342],
##        [ 3.45221699,  1.3808868 ,  0.81967592,  0.09730973],
##        [-6.50836319, -2.00334528, -0.78274557,  0.78835722]])

Se confrontiamo i porafogli sotto forma di grafico a barre possiamo dare i seguenti comandi.

barplot(omega,col=c('orange','red','blue','darkgreen'),
        beside=TRUE,
        names.arg=c('GOP','U','P','MV'))
legend('bottomright',legend=c('Amazon','Coke','Google','T-Bill'),
       col=c('orange','red','blue','darkgreen'),lty=1,lwd=10,bty='n')
grid()

plt.figure()
x=np.arange(0,4)
base=0.2
plt.bar(x-0.2,omega[0,:],base,label='Amazon',color='orange')
plt.bar(x,omega[1,:],base,label='Coke',color='red')
plt.bar(x+0.2,omega[2,:],base,label='Google',color='blue')
plt.bar(x+0.4,omega[3,:],base,label='T-Bill',color='green')
plt.xticks(x, ['GOP','U','P','MV'])
## ([<matplotlib.axis.XTick object at 0x717a9b14b450>, <matplotlib.axis.XTick object at 0x717ad82adcd0>, <matplotlib.axis.XTick object at 0x717ad8487c10>, <matplotlib.axis.XTick object at 0x717a9b182e90>], [Text(0, 0, 'GOP'), Text(1, 0, 'U'), Text(2, 0, 'P'), Text(3, 0, 'MV')])
plt.legend()
plt.show()

8.10 Il portafoglio attraverso il tempo

In questo paragrafo vediamo come rappresentare il portafoglio ottimo attraverso il tempo ricalcolando, con una frequenza prefissata, il portafoglio ottimo, aggiornando i dati di volta in volta.

La tecnica che usiamo qui viene chiamata “rolling” e si basa sulla definizine di una “finestra” temporale che, rimanendo costante, si sposta attraverso il tempo.

Se, per esempio, prendiamo un periodo di 90 giorni per la stima dei parametri del modello (medie, varianze e covarianze), allora:

  • la prima stima è calcolata sui dati dal n° 1 al n° 90
  • la seconda stima è calcolata sui dati dal n° 2 al n° 91
  • la terza stima è calcolata sui dati dal n° 3 al n° 92
  • e così via…

Per diversificare l’analisi prendiamo tre nuovi titoli. Caterpillar, Apple e Nike dal mercato USA e, quindi, usiamo come tasso privo di rischio, ancora una volta, il tasso sui T-Bills a 3 mesi (il “Federal Effective Rate” potrebbe essere una alternativa per il tasso privo di rischio).

Il primo passo è quello di scaricare i dati.

# Scarico i prezzi dei titoli rischiosi
getSymbols(c('AAPL','CAT','NKE'),src='yahoo',
           return.class='zoo',from='2019-01-01')
## [1] "AAPL" "CAT"  "NKE"
# Fondo le serie storiche
S=merge(APL=AAPL$AAPL.Adjusted,
        CAT=CAT$CAT.Adjusted,
        NIK=NKE$NKE.Adjusted,
        all=FALSE)

# Calcolo rendimenti logartimici
dlnS=diff(log(S))

# Scarico il tasso privo di rischio
getSymbols('DTB3',src='FRED',return.class='zoo')
## [1] "DTB3"
# Unisco i dati in un'unica matrice (eliminando gli NA)
dati=na.omit(merge(dlnS,DTB3/100,all=FALSE))
# Scarico i prezzi dei titoli rischiosi
S=yf.download(['AAPL','CAT','NKE'],start='2019-01-01')
## <string>:2: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [                       0%                       ][**********************67%*******                ]  2 of 3 completed[*********************100%***********************]  3 of 3 completed
S=S['Close']

# Calcolo rendimenti logartimici
dlnS=np.diff(np.log(S),axis=0)

# Scarico il tasso privo di rischio
DTB3=fred.get_series('DTB3',observation_start='2019-01-01')
DTB3=DTB3/100

Il passo successivo è quello di definire i dati del problema.

dt=1/250 # dati giornalieri
T=dim(dati)[1]*dt # numero periodi del data set (in anni)
n=dim(dati)[2] # numero di titoli (rischiosi + riskfree)

finestra=150 # finestra per il "rolling"
ptf=array(NA,dim=c(T/dt,n)) # portafoglio vuoto
dt=1/250 # dati giornalieri
T=np.shape(dlnS)[0]*dt # numero periodi del data set (in anni)
n=np.shape(dlnS)[1]+1 # numero di titoli (rischiosi + riskfree)

finestra=150 # finestra per il "rolling"
ptf=np.zeros((int(T/dt),n))

Si può procedere, infine, al calcolo dei portafogli “rolling” e alla loro rappresentazione grafica in serie storica. Qui vediamo l’evoluzione del portafoglio media-varianza, ma, ovviamente, lo stesso procedimento può essere applicato a qualsiasi tipo di portafoglio.

mu_E=0.15 # rendimento obiettivo

for (i in finestra:(T/dt)){
  #si prendono i rendimenti rischiosi (tranne la colonna n)
  rend=dati[(i-finestra+1):i,-n]
  #si prendono i tassi privi di rischio (n-esima colonna)
  tasso=dati[(i-finestra+1):i,n]
  MVC=var(rend)/dt
  mu=apply(rend,2,mean)/dt+0.5*diag(MVC)
  r=mean(tasso)
  xi2=as.numeric((mu-r)%*%solve(MVC)%*%(mu-r))
  ptf[i,-n]=(mu_E-r)/xi2*solve(MVC)%*%(mu-r)
  ptf[i,n]=1-sum(ptf[i,-n])
}

matplot(time(dati),ptf,type='l',lty=1,
        col=c('blue','orange','red','darkgreen'))
grid()
legend('topleft',legend=c('APL','CAT','NIK','T-Bill'),
       col=c('blue','orange','red','darkgreen'),lty=1,bty='n')

mu_E=0.15 # rendimento obiettivo

for i in range(finestra,int(T/dt)):
  #si prendono i rendimenti rischiosi (tranne la colonna n)
  rend=dlnS[(i-finestra):i,:]
  #si prendono i tassi privi di rischio (n-esima colonna)
  tasso=DTB3[(i-finestra):i]
  MVC=np.cov(rend,rowvar=False)/dt
  mu=np.mean(rend,axis=0)/dt+0.5*np.diagonal(MVC)
  r=np.mean(tasso)
  xi2=np.matmul(mu-r,np.matmul(np.linalg.inv(MVC),mu-r))
  ptf[i,0:n-1]=(mu_E-r)/xi2*np.matmul(np.linalg.inv(MVC),mu-r)
  ptf[i,n-1]=1-sum(ptf[i,0:n-1])

plt.figure()
plt.plot(ptf[:,0],color='blue',label='APL')
plt.plot(ptf[:,1],color='orange',label='CAT')
plt.plot(ptf[:,2],color='red',label='NIK')
plt.plot(ptf[:,3],color='darkgreen',label='TBill')
plt.legend()
plt.grid()
plt.show()

9 La misurazione del rischio

9.1 Creare una funzione

È possibile creare una funzione personalizzata mediante il comando function in R e il comando def in Python. Tra parentesi quadre bisogna indicare gli input. Questi input possono avere una valore di default (cioè un valore che viene loro attribuito se l’utente non ne indica uno) e questo va assegnato mediante il simbolo =.

Vediamo, per esempio, come scrivere una funzione che, ponendo come input due numeri, ne calcoli sia l asomma sia la differenza.

somdiff=function(a,b=1){
  S=a+b
  D=a-b
  list(somma=S,differenza=D)
}
def somdiff(a,b=1):
  S=a+b
  D=a-b
  return [S,D]

Vediamo come richiamare questa funzione:

  • gli input possono essere richiamati in qualsiasi ordine se si esplicita il loro “nome”
  • se non si esplicita il nome, gli input sono considerati dal software nello stesso ordine in cui sono stati dichiarati
  • se non si introducono gli input che hanno un vlore di default, tale valore viene automaticamente assegnato
somdiff(b=2,a=5)
## $somma
## [1] 7
## 
## $differenza
## [1] 3
somdiff(5)
## $somma
## [1] 6
## 
## $differenza
## [1] 4
somdiff(a=3,b=4)
## $somma
## [1] 7
## 
## $differenza
## [1] -1
somdiff(b=2,a=5)
## [7, 3]
somdiff(5)
## [6, 4]
somdiff(a=3,b=4)
## [7, -1]

Quando si utilizza una funzione si può richiamare solo uno dei sui output. Nel caso di R si può utilizzare la notazione con il simbolo del dollaro oppure la notazione matriciale con le paretensi quadre (che sono il default per Python).

somdiff(b=2,a=5)$somma
## [1] 7
somdiff(5)$differenza
## [1] 4
somdiff(a=3,b=4)[1]
## $somma
## [1] 7
somdiff(b=2,a=5)[0]
## 7
somdiff(5)[1]
## 4

9.2 Calcolo del VaR

Il Value at Risk è l’opposto di un quantile statistico. Il \(VaR_\alpha\), quindi, è l’opposto della perdita che lascia, alla sua sinistra, risultati economici che avvengono con una probabilità (cumulata) pari ad \(\alpha\).

Se \(\gamma\) è l’\(\alpha\) quantile, allora vale la relazione

\[F\left(-\gamma\right)=\alpha\]

da cui, essendo la funzione di ripartizione \(F\) invertibile, vale

\[\gamma=-F^{-1}\left(\alpha\right)\] Possiamo utilizzare i comandi già presenti su R e Python per calcolare i quantili. I metodi per il calcolo dei quantili sono molti. Noi, quindi, supponiamo di guardare semplicemente i dati storici senza alcuna forma di interpolazione tra i quantili. Questo metodo è definito type=1 in R e method='inverted_cdf' in Python.

PL=c(-7,-5,-2,0,1,5,10,12)
quantile(PL,probs=0.125,type=1)
## 12.5% 
##    -7
quantile(PL,probs=0.15,type=1)
## 15% 
##  -5
import numpy as np

PL=[-7,-5,-2,0,1,5,10,12]
np.quantile(PL,0.125,
  method='inverted_cdf')
## np.int64(-7)
np.quantile(PL,0.15,
  method='inverted_cdf')
## np.int64(-5)

9.3 Calcolo del VaR sui dati

Scarichiami i dati dell’indice S&P500 su tutta la serie storica disponibile da Yahoo e ne calcoliamo i rendimenti (logaritmici).

library(quantmod)

getSymbols('^GSPC',src='yahoo',return.class='zoo',from='1920-01-01')
## [1] "GSPC"
S=GSPC$GSPC.Adjusted

dlnS=diff(log(S))

plot(dlnS)

S=yf.download('^GSPC',start="1920-01-01")
## <string>:1: FutureWarning: YF.download() has changed argument auto_adjust default to True
## [*********************100%***********************]  1 of 1 completed
S=S['Close']

dlnS=np.diff(np.log(S['^GSPC']))

plt.figure()
plt.plot(dlnS)
plt.show()

Ora possiamo calcolare un vettore di valori giornalieri del VaR annuale (cioè su 250 osservazioni) a partire dalla 250-esima osservazioni e procedendo “rolling”. Poi rappresentiamo graficamente i valori del vettore VaR usando come valori delle ascisse le date dei prezzi, escludendo i primi 250 valori.

L=250 #lunghezza dell'intervallo temporale
N=length(dlnS) #numero totale di osservazioni
VaR=array(NA,dim=c(N-L+1,1)) #vettore vuoto

for (i in L:N){
  VaR[i-L+1]=-quantile(dlnS[(i-L+1):i],
                       probs=0.01,type=1)
}

plot(tail(time(dlnS),-L+1),VaR,
     type='l',xlab='')

L=250 #lunghezza dell'intervallo temporale
N=len(dlnS) #numero totale di osservazioni
VaR=np.zeros((N-L,1)) #vettore vuoto

for i in range(L,N):
  VaR[i-L]=-np.quantile(dlnS[(i-L):i],0.01,method='inverted_cdf')


plt.figure()
plt.plot(S.index[L+1:],VaR)
plt.show()

Osserviamo così quattro grandi picchi nella misura di rischio VaR:

  1. gli anni della crisi del 1929
  2. il lunedì nero del 1987 https://en.wikipedia.org/wiki/Black_Monday_(1987)
  3. la crisi del 2008
  4. la pandemia del 2020

9.4 L’Expected Shortfall

L’ES a livello di confidenza \(\alpha\), che è una misura di rischio coerente, è l’oppoto della media degli \(\alpha\) risultati peggiori, cioè dei risultati inferiori (o uguali) al VaR \(\alpha\). Formalmente, dunque, si può scrivere:

\[\mathbb{E}\left[\left.x\right|x\leq q_{\alpha}\right]\]

Anche in questo caso, se il quantile non è univoco, ricorriamo al metodo di spezzare la funzione di probabilità. Se esistono \(n\) osservazioni, tutte considerate equiprobabili, e vogliamo calcolare il quantile \(\alpha\) con il prodotto \(n\cdot\alpha\) che non è intero, allora effettuiamo i seguentie passaggi:

  1. si calcola il valore intero di \(n\cdot\alpha\) arrotondato per difetto \(k:=\left\lfloor n\alpha\right\rfloor\)
  2. si calcola la somma delle prime \(k\) osservazioni, pesate per la probabilità \(p:=\frac{1}{n\alpha}\)
  3. si aggiunge ancora l’osservazione \(k+1\) (in ordine crescente \(\overset{\rightarrow}{x}\)) pesata per la probabilità “mancante” data da \(1-p\cdot k\)
  4. si calcola l’opposto del valore così ottenuto

La formuala per il calcolo, dunque, è la seguente:

\[ES_{\alpha}=-\left(p\sum_{i=1}^{k}\overset{\rightarrow}{x_{i}}+\left(1-pk\right)\overset{\rightarrow}{x}_{k+1}\right)\]

Per il calcolo del VaR e dell’ES possiamo definire una funzione che riceva, in ingresso, la serie dei profitti e perdite e che restituisca il VaR e l’ES entrambi allo stesso livello di confidenza.

Nel codice R i utilizza il comando as.numeric per togliere dalla serie storica della variabile \(x\) qualsaisi etichetta (comprese le date)

vares=function(x,alpha){
  x=as.numeric(x)
  VaR=-quantile(x,probs=alpha,type=1)[[1]]
  k=floor(length(x)*alpha)
  p=1/length(x)
  ES=-(p*sum(sort(x)[1:k])+(alpha-p*k)*sort(x)[k+1])/alpha
  array(c(VaR,ES),dim=c(1,2),
        dimnames = list(NULL, c('VaR', 'ES')))
}
def vares(x,alpha):
  VaR=-np.quantile(x,alpha,method='inverted_cdf')
  k=int(len(x)*alpha)
  p=1/len(x)
  ES=-(p*sum(np.sort(x)[0:k])+(alpha-p*k)*np.sort(x)[k])/alpha
  return [VaR,ES]

9.5 Confronto VaR ed ES su dati reali

Sui dati S&P500 mettiamo a confronto, sullo stesso grafico, VaR ed ES mediante la stessa procedura che abbiamo già visto per il solo VaR.

L=250 #lunghezza dell'intervallo temporale
N=length(dlnS) #numero totale di osservazioni
alpha=0.01
VaRES=array(NA,dim=c(N-L+1,2)) #vettore vuoto

for (i in L:N){
  VaRES[(i-L+1),]=vares(dlnS[(i-L+1):i],alpha)
}

matplot(tail(time(dlnS),-L+1),VaRES,
        type='l',xlab='',ylab='',
        col=c('red','blue'),lty=1)
grid()
legend('topright',legend=c('VaR','ES'),
       col=c('red','blue'),lty=1,bty='n')

L=250 #lunghezza dell'intervallo temporale
N=len(dlnS) #numero totale di osservazioni
alpha=0.01
VaRES=np.zeros((N-L,2)) #vettore vuoto

for i in range(L,N):
  VaRES[i-L,:]=vares(dlnS[(i-L):i],alpha)


plt.figure()
plt.plot(S.index[L+1:],VaRES[:,0],
  color='red',label='VaR')
plt.plot(S.index[L+1:],VaRES[:,1],
  color='blue',label='ES')
plt.legend()
plt.show()