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

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

plot(SP500)

# Si caricano i pacchetti necessari
import numpy as np
import matplotlib.pyplot as plt
import yfinance as yf

S=yf.download('^GSPC',start="1920-01-01")
## 
[*********************100%***********************]  1 of 1 completed
SP500=S['Adj Close']

plt.figure()
plt.plot(SP500)
plt.show()

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))

# Pacchetto per gestire le date
import datetime as dt

plt.figure()
plt.plot(SP500)
plt.axvline(dt.datetime(1995,1,1), color = 'red',linestyle='dotted')
plt.text(dt.datetime(1995,1,1),np.max(SP500),'in sample <-',ha='right')
plt.text(dt.datetime(1995,1,1),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

2 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)) # differenze logaritmiche
sigma=np.sqrt(np.var(dlnS)/dt)
mu=np.mean(dlnS)/dt+0.5*sigma**2
mu
## 0.06638414723242678
sigma
## 0.18940631332079627

3 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] 7245 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[-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
## (7245, 1000)

4 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()