Volatilité stochastique¶
Nous considérerons l'EDS suivante modélisant les variations de prix d'un actif:
$$ dS_t = S_t \, ( \mu \, dt + \sigma_t \, dW_t)\qquad t\geq 0,\qquad\text{avec}\qquad S_0 = 100, $$
où $W$ est une mouvement Brownien, $\mu$ représente la tendance, et $(\sigma_t)_t$ la volatilité. Ici, la volatilité sera un processus stochastique défini par
$$ \sigma_t = f(Y_t)\qquad t\geq 0, $$
où $Y$ est la solution d'une EDS et $f$ une fonction. Dans la suite, on considérera un processus d'Ornstein-Uhlenbeck $Y$, c'est à dire que $Y$ vérifie
$$ dY_t = \alpha (m- Y_t)dt + \beta \, dZ_t. $$
On prendra aussi $f(y) = e^y$. On parle alors de modèle expOU pour la volatilité. On supposera aussi que $Z$ est un mouvement Brownien de la forme
$$ Z_t = \rho W_t + \sqrt{1-\rho^2} B_t \qquad t\geq 0, $$
où $\rho$ s'appelle le coefficient de corrélation instantané, $B$ est un mouvement Brownien indépendant de $W$.
Importons les librairies nécessaires et illustrons des trajectoires de la volatilité $(\sigma_t)_t$.
import numpy as np
import numpy.random as npr
import scipy.stats as st
import matplotlib.pyplot as plt
# fonction simulant une trajectoire du processus Y
def OU_simulation(T, n, m, alpha, beta, y0):
"""
T : durée de la simulation.
n : nombre de pas de temps.
alpha : vitesse de retour vers m.
m : valeur moyenne en temps long pour Y
beta : coefficient de diffusion.
y0 : condition initiale pour la volatilité.
"""
h = T / n # pas de temps
Y = np.zeros(n) # vecteur qui enregistrera les differentes valeurs de X
Y[0] = y0 # on enregistre la donnée initiale
for i in range(0, n-1): # schéma d'Euler-Maruyama
Y[i+1] = Y[i] + alpha * (m - Y[i]) * h + beta * np.sqrt(h) * npr.randn()
return Y
T = 20. # temps total de simulation
n = 1000 # nombre de pas de discrétisation
N = 20 # nombre de trajectoires
# initialisation des paramètres pour l'EDS de Y
m = np.log(0.1) # 0.1 est le comportement en temps long pour la volatilité \sigma_t, et donc log(0.1) pour Y
alpha = 1.
beta = np.sqrt(2)
y0 = 0.1 # condition initiale pour Y
t = np.linspace(0, T, n) # grille de temps
for _ in range(N): # boucle pour illustrer les différentes trajectoires
Y = OU_simulation(T, n, m, alpha, beta, y0)
plt.plot(t, np.exp(Y), lw=0.5, label="")
plt.xlabel("$t$", fontsize=15)
plt.ylabel("$\sigma_t$", fontsize=15)
plt.show()
On peut voir les excursions que peut faire la volatilité vers des valeurs assez grandes. Maintenant, illustrons l'effet d'une volatilité stochastique sur les trajectoires du prix de l'actif avec une volatilité stochastique.
# simulation d'une trajectoire du prix de l'actif avec une volatilité stochastique
def simulation_stock_expOU(T, n, mu, m, alpha, beta, rho, s0):
"""
T : durée de la simulation.
n : nombre de pas de temps.
mu: tendance pour le prix de l'actif
alpha : vitesse de retour vers m.
m : moyenne en temps long pour Y.
beta : coefficient de diffusion.
rho: coefficient de correlation instantanée
s0 : prix initial de l'actif.
"""
h = T / n # pas de temps
S = np.zeros(n) # vecteur qui enregistrera les differentes valeurs S_t de l'actif
S[0] = s0 # on enregistre la donnée initiale
Y = m # représentant les valeurs du processus Y_t
for i in range(0, n-1): # schéma d'Euler-Maruyama
V = npr.randn()
S[i+1] = S[i] * ( 1 + mu * h + np.exp(Y) * np.sqrt(h) * V )
Z = rho * V + np.sqrt(1-rho**2) * npr.randn()
Y = Y + alpha * (m - Y) * h + beta * np.sqrt(h) * Z
return S
T = 0.5 # temps total de simulation
n = 1000 # nombre de pas de discrétisation
N = 40 # nombre de trajectoires
# initialisation des paramètres pour l'EDS du prix de l'actif
mu = 0.15
s0 = 100.
# paramètres pour le modèle de volatilité
m = 0.1 # valeur en temps long de la volatilité
alpha = 1.
beta = np.sqrt(2)
rho = -0.2
t = np.linspace(0, T, n) # grille de temps
for _ in range(N):
S = simulation_stock_expOU(T, n, mu, np.log(m), alpha, beta, rho, s0)
plt.plot(t, S, lw=0.5, label="")
plt.xlabel("$t$", fontsize=15)
plt.ylabel("$S_t$", fontsize=15)
plt.show()
Sur cette illustration, on peut voir le large éventail des valeurs possibles de l'actif au temps final $T$. Ci-dessous, on va estimer la densité de $\log(S_T)$, et la comparer à la densité théorique dans le cas d'une volatilité constante $\bar \sigma$. On prendra $\bar \sigma=0.1$, qui est le comportement à long terme de notre volatilité stochastique.
Dans le cas d'une volatilité constante, on rappelle que $\log(S_T)$ suit une loi $\mathcal{N}((\mu - \bar \sigma^2/2)T,\bar \sigma^2 T)$.
# fonction simulant l'évolution de la volatilité et du prix de l'actif
# et renvoyant S_T
def simulation_stock_expOU_time_T(T, n, mu, m, alpha, beta, rho, s0):
"""
T : durée de la simulation.
n : nombre de pas de temps.
mu: tendance pour le prix de l'actif
alpha : vitesse de retour vers m.
m : moyenne en temps long pour Y.
beta : coefficient de diffusion.
rho: coefficient de correlation instantanée
s0 : prix initial de l'actif.
"""
h = T / n # pas de temps
S = s0 # représente les differentes valeurs de S
Y = m # représentant les valeurs du processus Y_t
for i in range(0, n-1): # schéma d'Euler-Maruyama
V = npr.randn()
S = S * ( 1 + mu * h + np.exp(Y) * np.sqrt(h) * V )
Z = rho * V + np.sqrt(1-rho**2) * npr.randn()
Y = Y + alpha * (m - Y) * h + beta * np.sqrt(h) * Z
return S
T = 0.5 # temps total de simulation
n = 1000 # nombre de pas de discrétisation
N = 10000 # nombre de trajectoires pour faire l'histogramme de S_T
# initialisation des paramètres pour l'EDS du prix de l'actif
mu = 0.15
s0 = 100.
# initialisation des paramètres du modèle de volatilité
m = 0.1
alpha = 1.
beta = np.sqrt(2)
rho = -0.2
t = np.linspace(0, T, n) # grille de temps
log_S_T = np.zeros(N) # vecteur qui enregistre les valeurs simulées de log_S_T
for i in range(N):
log_S_T[i] = np.log(simulation_stock_expOU_time_T(T, n, mu, np.log(m), alpha, beta, rho, s0))
log_S_T
array([4.62841559, 4.68141096, 4.69091693, ..., 4.65547508, 4.78497161,
4.63873242])
On peut donc comparer l'histrogramme normalisé de $\log(S_T)$ avec la densité de loi $\mathcal{N}((\mu - \bar \sigma^2/2)T,\bar \sigma^2 T)$.
s_bins = np.linspace(4, 5.20, 80) # définition des classes pour l'histogramme
# initialisation des paramètres pour la loi normale
mean = (mu - m**2/2)*T + np.log(s0)
std = m * np.sqrt(T)
y = st.norm(mean, std).pdf(s_bins)
# comparaison de l'histogramme normalisé de log(S_T) avec la densité de la loi normale
plt.hist(log_S_T, bins=s_bins, density=True, label="volatilité stochastique")
plt.plot(s_bins, y, lw=3, label="volatilité constante")
plt.xlabel("$\log(S_T)$", fontsize=15)
plt.legend()
plt.savefig("/home/cgomez/Nextcloud/cours/cours_calc_sto/Illustrations/stock_price_fat_tail.png")#plt.show()
On peut voir sur ce graphique, que l'histogramme de $\log(S_T)$, obtenu avec une volatilité stochastique, ne correspond pas à la densité de cette quantité dans le cas d'une volatilité constante. En particulier, on voit que l'histogramme est plus étalé sur les côtés que la densité de la loi normale prédite par le modèle à volatilité constante. On parle de queues épaisses (fat tails) pour la loi de $S_T$ avec un modèle de volatilité stochastique. Ce phénomène de fat tail est une carctéristique observée en pratique montrant que le modèle du mouvement Brownien géométrique avec une volatilité constante (que nous avons étudié en cours) n'est pas réaliste. Illustrons ce phénomène avec des données réelles. On va considérer le prix du Bitcoin en euro sur le moins d'octobre observé toutes les 15min. L'intérêt de regarder le prix du Bitcoin est que la cotation a lieu en continu, 7 jours sur 7 et 24h sur 24, ce qui facilite la modélisation mathématique.
Nous allons illustrer le phénomène de fat tail sur le log return, c'est à dire
$$ \log\Big(\frac{S_{t+h}}{S_t}\Big), $$
où $h$ est l'intervalle de temps de $15$ mins. Dans le cas d'une volatilité constante, avec un modèle de mouvement Brownien géométrique pour le prix du Bitcoin, on a
$$ \log\Big(\frac{S_{t+h}}{S_t}\Big) = \sigma(W_{t+h}-W_t) + \Big(\mu-\frac{\sigma^2}{2}\Big)h. $$
On remarque que les observations
$$ U_j = \log\Big(\frac{S_{t_j+h}}{S_{t_j}}\Big),\qquad j=0,\dots, N, $$
sont de loi normale et indépendantes car un mouvement Brownien est un processus à accroisement indépendant. Pour déterminer les paramètres de cette loi normale, on prend l'espérance et la variance de cette relation, et on obtient le système
$$ \mathbb{E}\Big[\log\Big(\frac{S_{t+h}}{S_t}\Big)\Big] = \Big(\mu-\frac{\sigma^2}{2}\Big)h \qquad\text{and}\qquad Var\Big(\log\Big(\frac{S_{t+h}}{S_t}\Big)\Big) = \sigma^2 h, $$
qu'on peut facilement résoudre:
$$ \sigma^2 h =Var\Big(\log\Big(\frac{S_{t+h}}{S_t}\Big)\Big) \qquad\text{and}\qquad \mu h = \mathbb{E}\Big[\log\Big(\frac{S_{t+h}}{S_t}\Big)\Big] + \frac{1}{2}Var\Big(\log\Big(\frac{S_{t+h}}{S_t}\Big)\Big). $$
On peut alors estimer ces paramètres par
$$ \sigma^2 h = \frac{1}{N}\sum_{j=1}^N (U_j - \bar U_N)^2 \qquad\text{and}\qquad \mu h = \bar U_N + \frac{1}{2N}\sum_{j=1}^N (U_j - \bar U_N)^2. $$
où
$$ \bar U_N = \frac{1}{N}\sum_{j=1}^N U_j, $$
Comme nous l'avons vu au Chapitre 5.
# librairie nécessaire à l'importation et la gestion des données
import yfinance as yf
import pandas as pd
# téléchargement des données
ticker = 'BTC-EUR' # code sur yahoo finance de la cotation qui nous intéresse
start_date = '2024-10-01' # date de début pour les données
end_date = '2024-10-31' # date de fin
btceur_data = yf.download(ticker, start=start_date, end=end_date, interval='15m') # téléchargement des données
# Calcul des log returns (creation de la colonne dans la dataframe)
btceur_data['Log_return'] = np.log(btceur_data['Close'] / btceur_data['Close'].shift(1))
n_time = btceur_data.shape[0]-1
[*********************100%***********************] 1 of 1 completed
btceur_data = btceur_data.dropna() # retire les valeurs manquantes
# Illustration de l'évolution des log returns en fonction du temps
plt.figure(figsize=(12, 4))
plt.plot(btceur_data.Log_return)
plt.title('Log Returns of BTC-EUR')
plt.xlabel('Date', fontsize=15)
plt.ylabel('Log Return', fontsize=15)
plt.show()
A partir de ces données, estimons les paramètres de la loi normale prédite par un modèle à volatilité constante.
# estimation des paramètres de la loi normale
sig = np.sqrt(np.var(btceur_data.Log_return))
mu = np.mean(btceur_data.Log_return) + sig**2/2
# création de la grille de valeur pour l'histogramme des log returns
u_min = btceur_data.Log_return.min()
u_max = btceur_data.Log_return.max()
u = np.linspace(u_min, u_max, 100)
# évaluation de la densité de la loi normale sur la grille avec les paramètres estimés
y = st.norm(mu, sig).pdf(u)
# Histogramme des log returns et de la densité de la loi normale
# prédite par le modèle à volatilité constante
plt.hist(btceur_data.Log_return, bins=u, density=True, label="data log returns", alpha=0.3)
plt.plot(u, y, lw=3, label="denisté loi normale")
plt.xlabel("15min log returns")
plt.xlim(-0.01, 0.01)
plt.legend()
plt.show()
On peut observer sur ce graphique que la loi des log returns observés ne correspond pas à une loi normale comme le modèle de mouvement Brownien géométriqueque à volatilité constante prédit. On peut voir sur l'histogramme quelques observations là où les probabilités d'apparition sont très faibles pour une loi normale. Ceux sont ces valeurs qui produisent une variance trop grandes pour la loi normale par rapport aux données, et c'est bien ce phénomène de fat tail, et ces observations qui dévient beaucoup de la moyenne, qui produise cette trop grande variance pour la loi normale vis à vis des données.
L'introduction d'un modèle de volatilité stochastique permet de corriger ce problème, et d'obtenir une loi pour les log returns qui correspond beaucoup mieux aux données observées (voir cours).