Simulation numérique de quelques EDS¶

Nous présentons ici le schéma d'Euler-Maruyama appliqué à deux EDS standards.

In [1]:
# import des librairies nécessaires aux simulations numériques

import numpy as np
import numpy.random as npr
import scipy.stats as st
import matplotlib.pyplot as plt

Processus d'Ornstein-Uhlenbeck¶

Etant donné un espace de probabilité filtré et $ W $ un mouvement Brownien, un processus d'Ornstein-Uhlenbeck est la solution de l'EDS

$$ dX_t = \theta (m - X_t) \, dt + \sigma \, dW_t \qquad X_0 = x_0 $$

où :

  • $ m $ est la valeur vers laquelle le processus tend à revenir (valeur moyenne en temps long),
  • $ \theta >0 $ est le paramètre de vitesse de retour vers $m$ (vitesse de retour à la moyenne),
  • $ \sigma $ représente l'intensité des fluctuations aléatoires (appélé coefficient de diffusion).

Ce processus permet de modéliser des phénomènes physiques, économiques, financiers, etc...

In [2]:
#Discrétisation de l'EDS par un schéma d'Euler-Maruyama

def OU_simulation(T, n, m, theta, sigma, x0):
    """
    T : durée de la simulation.
    n : nombre de pas de temps.
    theta : vitesse de retour vers m.
    m : paramètre vers lequel le processus revient.
    sigma : coefficient de diffusion.
    x0 : condition initiale du processus.
    """
    
    h = T / n                 # pas de temps
    X = np.zeros(n)           # vecteur qui enregistrera les differentes valeurs de X 
    X[0] = x0                 # on enregistre la donnée initiale

    for i in range(0, n-1): # Schéma d'Euler-Maruyama
        X[i+1] = X[i] + theta * (m - X[i]) * h + sigma * np.sqrt(h) * npr.randn() 

    return X
In [3]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 0.       # initialisation des paramètres de l'EDS
theta = 1.
sigma = 1.5

# Simulation d'une trajectoire de processus
x0 = 10.       # condition initiale
X = OU_simulation(T, n, m, theta, sigma, x0)

t = np.linspace(0, T, n) # grille de temps

#illustration
plt.plot(t, X, label="$x_0=0")

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)
plt.show()
No description has been provided for this image

On peut dessiner des trajectoires partant de conditions initiales différentes.

In [4]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 0.
theta = 1.
sigma = 1.5

# Simulation de trois trajectoires de processus pour différentes conditions initiales
x0 = 0.
X0 = OU_simulation(T, n, m, theta, sigma, x0)

x01 = 10.
X1 = OU_simulation(T, n, m, theta, sigma, x01)

x02 = -10.
X2 = OU_simulation(T, n, m, theta, sigma, x02)

t = np.linspace(0, T, n)

plt.plot(t, X0, label="$x_0=%s$" %x0)
plt.plot(t, X1, label="$x_0=%s$" %x01)
plt.plot(t, X2, label="$x_0=%s$" %x02)

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)

plt.legend()
plt.show()
No description has been provided for this image

On voit bien que, malgré des points de départ différents, les trajectoires finissent par tourner autour de $m=0$. Pour une condition initiale donnée on peut allustrer plusieurs trajectoires.

In [5]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 0.
theta = 1.
sigma = 1.5

t = np.linspace(0, T, n)
x0 = 10.

N = 20 # nombre de trajectoires

for _ in range(N):
    X = OU_simulation(T, n, m, theta, sigma, x0)
    plt.plot(t, X, lw=0.5, label="")

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)

plt.show()
No description has been provided for this image

On peut aussi évaluer numériquement le comportement en temps long de ce processus, c'est à dire la loi de $X_t$.

In [6]:
# version de la fonction précédente n'engeristrant pas toute la trajectoire,
# mais ne renvoyant uniquement que la position final au temps T

def OU_simulation_time_T(T, n, m, theta, sigma, x0):
    """
    T : durée de la simulation.
    n : nombre de pas de temps.
    theta : vitesse de retour vers m.
    m : paramètre vers lequel le processus revient.
    sigma : coefficient de diffusion.
    x0 : condition initiale du processus.
    """
    
    h = T / n                 # pas de temps
    X = x0                    # initialisation en la condition initiales  

    for i in range(0, n-1):
        X = X + theta * (m - X) * h + sigma * np.sqrt(h) * npr.randn() 

    return X
In [7]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 0.
theta = 1.
sigma = 1.5

t = np.linspace(0, T, n)
x0 = 10.

N = 20000 # nombre de trajectoires

XT = np.zeros(N) # variable qui correspondra au processus X au temps final X_T. 

for j in range(N):
    XT[j] = OU_simulation_time_T(T, n, m, theta, sigma, x0)
In [8]:
x = np.linspace(-4, 4, 100) # définition des classes pour l'histogramme
plt.hist(XT, bins=x, density=True, label="histogramme $X_t$") 
plt.plot(x, st.norm(m, 1/np.sqrt(theta)).pdf(x), lw=3, label="densité loi normale")
plt.legend()
plt.show()
No description has been provided for this image

L'approximation avec la loi $\mathcal{N}(m,1/\sqrt{\theta})$ est plutôt bonne, tout en sachant que $X_t$ ne suit pas exactement cette loi. C'est $\lim_{t\to +\infty} X_t$ qui suit cette loi (voir Exercice 4.4.2).

Modèle CIR (Cox-Ingersoll-Ross)¶

Ce modèle correspond à l'EDS

$$ d X_t = \theta (m - X_t) dt + \sigma \sqrt{X_t} \, dW_t \qquad X_0 = x_0, $$

où les paramètres ont la même interprétation que précédemment. Ce modèle est utilisé en finance pour modéliser l'évoltion des taux d'intérêt, ou une volatilité non constante (volatilité stochastique/stochastic volatility). Nous reviendrons sur ce sujet dans le dernier chapitre de ce cours.

On peut faire le même type de graphique que précédemment.

In [9]:
def CIR_simulation(T, n, m, theta, sigma, x0):
    """
    T : durée de la simulation.
    n : nombre de pas de temps.
    theta : vitesse de retour vers m.
    m : paramètre vers lequel le processus revient.
    sigma : coefficient de diffusion.
    x0 : condition initiale du processus.
    """
    
    h = T / n                 # pas de temps
    X = np.zeros(n)           # vecteur qui enregistrera les differentes valeurs de X 
    X[0] = x0                 # on enregistre la donnée initiale

    for i in range(0, n-1):
        X[i+1] = X[i] + theta * (m - X[i]) * h + sigma * np.sqrt(max(X[i],0.) * h) * npr.randn() 

    return X
In [10]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 1.      # initialisation des paramètres de l'EDS
theta = 1.
sigma = 1.5

# Simulation d'une trajectoire de processus
x0 = 10.       # condition initiale
X = CIR_simulation(T, n, m, theta, sigma, x0)

t = np.linspace(0, T, n) # grille de temps

# illustration
plt.plot(t, X, label="$x_0=0")

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)
plt.show()
No description has been provided for this image

Pour 3 trajectoires partant de 3 conditions initiales différentes.

In [11]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 1.
theta = 1.
sigma = 1.5

# Simulation de trois trajectoires de processus pour différentes conditions initiales
x0 = 0.
X0 = CIR_simulation(T, n, m, theta, sigma, x0)

x01 = 5.
X1 = CIR_simulation(T, n, m, theta, sigma, x01)

x02 = 10.
X2 = CIR_simulation(T, n, m, theta, sigma, x02)

t = np.linspace(0, T, n)

plt.plot(t, X0, label="$x_0=%s$" %x0)
plt.plot(t, X1, label="$x_0=%s$" %x01)
plt.plot(t, X2, label="$x_0=%s$" %x02)

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)

plt.legend()
plt.show()
No description has been provided for this image

Pour plusiseurs trajectoires partant de la même condition initiale.

In [12]:
T = 10.    # temps total de simulation
n = 1000   # nombre de pas de discrétisation

m = 1.
theta = 1.
sigma = 1.5

t = np.linspace(0, T, n)
x0 = 5.

N = 20 # nombre de trajectoires

for _ in range(N):
    X = CIR_simulation(T, n, m, theta, sigma, x0)
    plt.plot(t, X, lw=0.5, label="")

plt.xlabel("$t$", fontsize=15)
plt.ylabel("$X_t$", fontsize=15)

plt.show()
No description has been provided for this image