Copule et chaîne de Markov

2. Copule et chaîne de Markov#

Dans cette feuille, on s’intéresse à la simulation de vecteurs aléatoires.

Dans la première partie on simule un vecteur $(\tau_1, \dots, \tau_n)$ dont la dépendance entre les composantes est donnée par une fonction copule, on parle de dépendance spatiale.
Dans la deuxième partie, on simule un vecteur $(X_0, \dots, X_n)$ dont la dépendance est temporelle, c’est à dire qu’on construit $X_{k+1}$ à partir de $X_k$ et d’un aléa indépendant. C’est un exemple de chaîne de Markov.

2.1. Copule Gaussienne et instants de défaut#

On considère la modélisation suivante: soit $\tau_1, \dots, \tau_n$ des variables aléatoires exponentielles de paramètre $\lambda_i > 0$ qui représentent $n$ instants de défaut (instant de défaillance d’un composant en fiabilité, instant de défaut d’une entreprise en finance, instant de mutation en biologie, etc). On suppose que ces $n$ instants sont corrélés par la copule Gaussienne $C$ de matrice de covariance $\Sigma_{i,i} = 1$ et $\Sigma_{i,j} = \rho \in [0,1]$ pour $i \neq j$.

Dans un premier temps on s’intéresse à la simulation du vecteur $U = (U_1,\dots,U_n)$ correspondant à la copule $C$ i.e. pour tout $i \in \{1,\dots,n\}, U_i \sim \mathcal{U}([0,1])$ et $C(u_1,\dots,u_n) = \mathbf{P}[U_1 \le u_1,\dots,U_n\le u_n]$. On peut montrer que $U$ se représente sous la forme

\[ \forall i \in \{1,\dots,n\}, \quad U_i = \Phi\bigl(\sqrt{\rho} G_0 + \sqrt{1 - \rho} G_i \bigr), \]

où $(G_0,G_1,\dots, G_n)$ est un vecteur normal standard de $\mathbf{R}^{n+1}$ et $\Phi$ est la fonction de répartition gaussienne standard (fonction norm dans le module scipy.stats).

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() 
from numpy.random import default_rng
rng = default_rng()

2.1.1. Question: simulation de la copule gaussienne#

Ecrire une fonction qui prend 3 arguments: size, n et rho et qui renvoie un échantillon de taille size de réalisations $(U_1,\dots,U_n)$ obtenues par la construction précédente.
Dans toute la suite n est la dimension de la copule qui sera par exemple 10, et non le nombre de réalisations considérées.

2.1.2. Question: représentation de la copule gaussienne#

Reproduire le tracé suivant à partir d’échantillons $\big(U_1^{(j)}, U_2^{(j)}\big)_{1 \le j \le M}$ de taille $M = 2000$ de réalisations de la copule gaussienne en dimension 2 et de corrélation $\rho= 0.2$ puis $\rho = 0.8$. Que remarquez-vous?

Compléter ce tracé en vérifiant graphiquement que la loi empirique de chaque composante $\big(U_i^{(j)}\big)_{1 \le j \le M}$ pour $i=1,2$ est proche de la loi uniforme sur $[0,1]$. Vous pouvez augmenter la taille de $M$.

On revient maintenant au problème des instants de défaut dont on rappelle que $\tau_i \sim \mathcal{E}(\lambda_i)$. On s’intéresse à la variable aléatoire

\[ X = \sum_{i = 1}^n \mathbf{1}_{\tau_i \le 1}. \]

C’est donc une variable aléatoire discrète à valeurs dans $\{0,\dots,n \}$.

2.1.3. Question: simulation des instants de défaut#

Ecrire une fonction (avec les arguments adhoc) pour simuler un échantillon de $X$.
On utilise dans la suite $n =20$ instants de défaut et $\lambda_i = 0.4$ pour tout $i \in \{1,\dots,n\}$.

2.1.4. Question: création d’un ensemble de scénarios `samples`#

Créer un vecteur rhos de taille 4 avec les valeurs $(0,0.2,0.6,0.8)$ et créer une matrice samples de 4 lignes et 20000 colonnes, dont chaque ligne sera remplie par des simulations de $X$ pour une valeur de rho donnée. Remplir cette matrice avec réalisations de $X$. Il faut que les éléments de la matrice sample soit de type entier np.int64.

Cette matrice sera utilisée dans la suite pour le calcul de la moyenne et le tracé des histogrammes.

2.1.5. Question: calcul de la moyenne#

Calculer (mathématiquement) $\mathbf{E}[X]$ en fonction de $\rho$ et vérifier ce résultat en calculant la moyenne empirique de chaque ligne de la matrice samples des réalisations de $X$.

2.1.6. Question: représentation graphique de `samples`#

2.2. Simulation d’une chaîne de Markov: Urnes d’Ehrenfest#

On considère $d$ balles ($d > 1$) numérotées de 1 à $d$ et réparties initialement dans deux urnes $A$ et $B$. On note $E = \{1,\dots,d\}$ l’ensemble des balles et on s’intéresse à l’évolution du contenu des urnes après un nombre $n \ge 1$ de changements d’états. Un changement d’état est modélisé de la façon suivante: “on tire un numéro de balle selon la loi uniforme sur $E$ et à un tirage $i$ on déplace la balle numéro $i$ d’une urne à l’autre”. Le contenu des urnes change au cours du temps et on note $A_n$ le contenu de $A$ à l’itération $n$ (de même $B_n$ est le contenu de $B$ au temps $n$). Ce modèle porte le nom d’Urnes d’Ehrenfest.

Dans le programme de test on prendra $d = 20$ c’est à dire $E = \{1,\dots,20\}$ et le contenu initial $A_0 = \{1,\dots,10\}$.

2.2.1. Approche naïve#

On commence par une approche naïve qui consiste à programmer l’évolution du contenu de l’urne. Si on s’intéresse uniquement au nombre de balles dans chaque urne, il est plus commode d’utiliser la propriété de Markov du système, aussi bien pour l’étude mathématique que pour la programmation.

2.2.1.1. Question: pour manipuler le type `list`#

En utilisant le type list initaliser les listes E et A_0, puis programmer l’évolution du contenu de l’urne pas à pas (par exemple en utilisant les méthodes remove et append). Afficher le résultat A_n et B_n après n = 100 itérations.

2.2.1.2. Question: pour manipuler le type `set`#

Reprendre la question précédente avec le type set et les méthodes associées: add et remove.

2.2.1.3. Question facultative#

Cette construction itérative nécessaire dans des dynamiques plus complexes peut être améliorée de la façon suivante. Si on construit un vecteur sample des n réalisations de la loi uniforme sur $E$, on peut en déduire directement la composition des urnes à l’itération $n$. En effet la balle $i$ n’a pas changé d’urne entre 0 et $n$ si et seulement si le numéro $i$ apparait un nombre pair de fois dans sample.

Construire directement A_n et B_n à partir d’un vecteur sample de taille n qui contient les réalisations de la loi uniforme sur $E$.

2.2.2. Modélisation en tant que chaîne de Markov#

On s’intéresse maintenant non plus à la composition de l’urne $A_n$ mais uniquement à sa taille. On note $X_n = \mathrm{Card}(A_n)$. L’évolution de $X_n$ se fait de la façon suivante: si l’urne $A_n$ contient $X_{n}$ balles alors la probabilité de tirer une balle présente dans $A_n$ est $\frac{X_{n}}{d}$. Ainsi avec probabilité $\frac{X_{n}}{d}$, $X_{n+1} = X_n - 1$ (car on déplace la balle dans l’urne $B$), et avec probabilité $\frac{d-X_n}{d}$ on a $X_{n+1} = X_n + 1$. La récurrence aléatoire suivante permet de construire une trajectoire $(X_0, \dots, X_n)$

\[ \forall k \ge 0, \quad X_{k+1} = X_k + 1 - 2 \mathbf{1}_{ \left\{U_{k+1} < \frac{X_k}{d} \right\}}, \qquad X_0 \in \{0,\dots,d\}, \]

avec $(U_k)_{k \ge 1}$ une suite de variables aléatoires indépendantes de loi uniforme sur $[0,1]$.

2.2.2.1. Question: simulation de $N$ trajectoires#

Ecrire une fonction qui permet de simuler $N$ trajectoires indépendantes $(X_0^{(j)},\dots, X_n^{(j)})_{1 \le j \le N}$ de la dynamique précédente, et qui renvoie donc un np.array de dimension (N, n+1) qui contient ces trajectoires. A vous déterminer les arguments de cette fonction.

2.2.2.2. Question: affichage de trajectoires#

Reproduire le graphe suivant qui représente l’évolution de 5 trajectoires de $(X_k)_{0 \le k \le n}$ avec $n=20$.

2.2.2.3. Question: représentation de la loi après $n$ itérations#

Représenter la distribution empirique de $(X^{(j)}_n)_{1 \le j \le N}$ pour $N=100\,000$ et deux valeurs de $n$, $n = 100$ puis $n = 101$. Que remarque-t-on? Etait-ce prévisible?

2.2.2.4. Question: modification de la dynamique#

On modifie un peu la modélisation précédente en considérant la règle suivante: “on tire un numéro de balle selon la loi uniforme sur $E$ et à un tirage $i$ on déplace la balle numéro $i$ d’une urne à l’autre avec probabilité 1/2”.

Refaire les questions précédentes avec ce nouveau modèle. On peut montrer que si $X_{n} \sim \mathcal{B}\bigl(d, \frac{1}{2}\bigr)$ (loi binomiale) alors $X_{n+1} \sim \mathcal{B}\bigl(d, \frac{1}{2}\bigr)$. On dit que la loi binomiale $\mathcal{B}\bigl(d, \frac{1}{2}\bigr)$ est invariante pour la chaine de Markov. De plus pour toute configuration initiale $X_0$ la chaine $(X_n)_{n \ge 0}$ converge en loi vers cette mesure invariante $\mathcal{B}\bigl(d, \frac{1}{2}\bigr)$. On veut illustrer cette convergence.

Comparer cette loi binomiale avec l’histogramme empirique de la loi $X_n$ pour $n$ grand (par exemple $n = 100$, vous pouvez choisir $X_0$ fixé à 10).

2.3. Processus de Poisson#

On considère un processus de Poisson de paramètre (ou intensité) $\lambda > 0$, c’est à dire un processus de comptage associé à un processus ponctuel $(T_n)_{n \ge 1}$ où les variables aléatoires $T_n$ (appelées instants de sauts) sont définies par

\[\begin{equation*} \forall n \ge 1, \quad T_n - T_{n-1} = S_n, \qquad \text{en posant $T_0 = 0$} \end{equation*}\]

avec $(S_n)_{n \ge 1}$ suite i.i.d. de loi exponentielle de paramètre $\lambda > 0$.

Pour tout $t \ge 0$, on définit

\[\begin{equation*} N_t = \sum_{n \ge 0} \mathbf{1}_{T_n \le t}, \end{equation*}\]

et on veut simuler une trajectoire de $(N_t)_{t \in [0,T]}$ pour un horizon $T > 0$ fixé.

2.3.1. Question: simulation de trajectoires#

Ecrire une fonction one_poisson_path(lambd, T) qui renvoie les instants $(T_0, \dots, T_n)$ d’un processus de Poisson d’intensité lambd restreint à $[0, T]$. Par convention, on veut que le dernier point du tableau renvoyé soit $T$.
Reproduire le tracé suivant.

2.3.2. Question: simulation alternative#

On rappelle que si $(N_t)_{t \ge 0}$ est un processus de Poisson d’intensité $\lambda > 0$, alors conditionnellement à l’événement $N_T = n$ les instants de sauts $(T_k)_{k=1,\dots,n}$ (tels que $0 < T_1 < \dots < T_n \le T$) ont même loi que le réordonnement croissant d’un vecteur $(U_1, \dots, U_n) \sim \mathcal{U}([0,T]^n)$.