Introduction aux Processus Stochastiques (PRSTO)

P. Carmona

Pour avancer dans le notebook et exécuter les cellules il faut taper Shift+Enter ou utiliser la barre d'outils ci-dessus et choisir Cell, Run Cell and select Below

Consignes pour ce TP temps réel noté

Vous répondrez aux questions en modifiant ce notebook. En insérant des cellules de type Markdown pour le texte et des cellules de type code pour le code.

Ensuite vous sauvez ce notebook sous le nom Prenom_Nom_tpidla.ipynb et vous le déposez sur Moodle

La ruine du joueur

$$ \def\etp#1{{\left ( {#1} \right )}} \def\prob#1{{{\mathbb{P}}\left ( {#1} \right )}} \def\esp#1{{{\mathbb{E}}\left [ {#1} \right ]}} \def\ens#1{{\left\{#1\right\}}} \def\etc#1{{\left [ {#1} \right ]}} $$

Soit $S_n$ une marche aléatoire simple issue de $x$. On peut réaliser $S_n=x+\xi_1 + \cdots + \xi_n$ avec $(\xi_i)_{i\ge 1} $ iid de loi $P(\xi_1=\pm 1 ) = 1/2$. Etant donnés des entiers $a\le x\le b$ on pose $T_x = \inf\{n\ge 0 : S_n = x\}$ et en considérant la martingale $S_n$ arrêtée en $T=T_a \wedge T_b$ on montre que $E_x[S_T]=E_x[S_0]$ et on en déduit (ce résultat a été montré en cours et on ne vous demande pas de le redémontrer) $$ \mathbb{P}_x(T_a < T_b) = \frac{b-x}{b-a}\,.$$ Illustrer ce résultat par des simulations en traçant, pour $a=-10,b=10$ une approximation de la fonction $ x\to \mathbb{P}_x(T_a < T_b)$.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rnd
%matplotlib inline

IDLA

On étudie le modèle de croissance/exploration suivant. Une particule part de l'origine et se comporte comme une marche aléatoire simple jusqu'à ce qu'elle rencontre une zone inexplorée, et la elle se fixe sur le premier site inexploré. Si $A_n$ désigne l'ensemble aléatoire des sites de $\mathbb{Z}$ explorés au temps $n$, on a $A_0=\ens{0}$ et $A_1=\ens{0,1}$ si la particule fait un premier pas à droite, $A_1=\ens{-1,0}$ sinon. Ainsi étant donné $A_n$, on fait partir une marche aléatoire de $0$ et on pose $T= \inf\ens{k\ge 1 : S_n \notin A_n}$ et on pose $A_{n+1}=A_n \cup \ens{S_T}$. On pose $A_n=\etc{G_n,D_n}$ et on remarque que $D_n -G_n=n$ (pourquoi ?). L'objet de notre étude est $X_n = G_n + D_n$.

Ecrire une fonction spdla(n) qui prend pour paramètre l'instant $n$ et renvoie la valeur de $X_n$ obtenue par la simulation décrite précédemment. Attention: vous devez attendre le premier instant où la marche sort de l'intervalle $A_n$ déjà exploré. Vous devez donc utiliser une boucle while (tant que).

In [ ]:

En utilisant les résultats précédents sur la ruine du joueur on peut montrer que en fait $(X_n,n\in\mathbb{N})$ est une chaîne de Markov inhomogène i.e. une suite $(X_n)_{n\in\N}$ de variables aléatoires entières dont les lois sont données par la condition initiale $X_0=0$ p.s. et la récurrence $$ p_{n,i}=\prob{X_{n+1}=i+1 \mid X_n =i} = 1-\prob{X_{n+1}=i-1 \mid X_n =i}=\frac{n+2-i}{2(n+2)}$$

En utilisant la relation précédente écrire une fonction mhdla(h) qui prend pour paramètre l'instant $n$ et renvoie la valeur de $X_n$ obtenue. C'est donc un programme de simulation analogue à celui permettant de simuler une chaine de Markov homogène : la différence est que la chaîne est inhomogène, et que donc nous vous conseillons de coder d'abord une fonction pidla(n,i) qui calcule la probabilité $p_{n,i}$

On peut tirer une variable de Bernoulli de parametre $p$ en écrivant

In [2]:
p=0.3
rnd.binomial(1,p)
Out[2]:
0

Pour $n=10$ et $n=20$, (pourquoi ne prend on pas $n=1000$ ?) générer un grand échantillon de taille $N=500$ de variables aléatoires de loi $X_n$ en utilisant les deux fonctions précédentes. Tracer des histogrammes de ces deux échantillons sur un même graphique. Comparer numériquement les moyennes empirique et variances empirique de ces échantillons. Commentez.

Comme superposer des histogrammes n'est pas évident nous vous conseillons d'utiliser des estimations de la densité (kernel density estimator). A titre d'exemple nous utilisons cette méthode ci-dessous pour deux echantillons de taille 500 de variables aléatoire suivants une loi normale centrée réduite. Il vous appartient d'adater cet exemple pour comparer vos echantillon obtenus avec spdla() et mhdla()

In [3]:
from scipy.stats import kde
 
# create data
x = np.random.normal(size=500)
x.sort()
y = np.random.normal(size=500)
y.sort()
# Evaluate a gaussian kde 

k = kde.gaussian_kde(x)
plt.plot(x,k(x),label="x")
kp = kde.gaussian_kde(y)
plt.plot(y,kp(y),label="y")
plt.legend()
Out[3]:
<matplotlib.legend.Legend at 0x7fa629834ef0>

IDLA Limite

On peut montrer (facilement par récurrence) que $\esp{X_n}=0$ et que si $a_n=\esp{X_n^2}$, alors $a_{n+1}= 1 + \frac{n}{n+2} a_n$. On en déduit que $\lim \frac{a_n}{n} = 1/3$ (un exercice de prépa pas si évident qu'il n'y paraît).

On peut en fait montrer la convergence en loi (admise) $$ \frac{X_n}{\sqrt{n}} \xrightarrow[n\to \infty]{\mathcal{L}} \mathcal{N}(0,1/3)\,.$$ Pour illustrer la convergence en loi, choisir $n$ assez grand, fabriquer un $N$ échantillon, avec $N$ assez grand, de la loi de $X_n$, $$ X_n^{(1)}, X_n^{(2)}, \ldots, X_n^{(N)}$$ Comparer la moyenne, la variance de cet échantillon avec la moyenne et la variance de la loi cible.

Enfin, superposer l'histogramme de l'échantillon avec la densité de la loi cible.

In [ ]:
In [ ]: