Processi a diffusione e salto
2023-10-18
1 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)
# Si caricano i pacchetti necessari
import numpy as np
import matplotlib.pyplot as plt
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()
2 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()
3 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()
4 Rappresentazione dei rendimenti (logaritmici)
plot(diff(log(X[,1])),type='l')
plt.figure()
plt.plot(np.diff(np.log(X[:,0])))
plt.show()
5 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,075229223 0,188893582 -0,462448879 0,000755774# Scarichiamo i dati S&P500
import yfinance as yf
S=yf.download('^GSPC',start="1920-01-01")##
[*********************100%***********************] 1 of 1 completedSP=S['Adj Close']
dt=1/250
# Calcoliamo i rendimenti logaritmici
dlnS=np.diff(np.log(SP))
# 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## (0.07522847518134426, 0.18888962285238553, -0.46244880843876945, 0.0007557746533071398)6 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,032535# Calcolo della frequenza relativa dei salti
length(salti)/length(dlnS)## [1] 0,04654642# 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
# Calcolo della frequenza relativa dei salti## -1.0325355714732578len(salti)/len(dlnS)## 0.04654642174382844