\documentclass[a4paper, 11pt, french]{article}
\usepackage{babel}
\usepackage{fontspec}
\usepackage{amssymb, amsmath}
\usepackage[pagebackref=true,breaklinks=true,colorlinks=true,backref=false,pdfencoding=auto]{hyperref}
\usepackage[capitalize,noabbrev]{cleveref}
\usepackage[backend=biber,style=trad-plain]{biblatex}
\usepackage{url}
\addbibresource{rapport.bib}
\newcommand{\E}{\mathbb{E}}
\newcommand{\V}{\mathbb{V}}
\newcommand\floor[1]{\lfloor#1\rfloor}
\DeclareMathOperator*{\argmax}{arg\,min}

\begin{document}
\section*{Introduction}
\paragraph{Sujet :\\}
Comparer les performances des méthodes de stratification exposées dans \cite{etore:hal-00192540}
et de QMC randomisées exposées dans \cite{tuffin2004randomization} sur les exemples présentés dans \cite{etore:hal-00192540}

Références : \cite{etore:hal-00192540} \cite{tuffin2004randomization}

\section{Quasi Monte Carlo randomisé}
\subsection{Présentation mathématique}
On suppose qu'on cherche à calculer
\begin{equation*}
I=\int_{[0,1]^s}f(x)dx
\end{equation*}
Soit $(\xi^{(n)})_{n\in\mathbb{N}}$ une suite à discrépance faible à valeurs dans $[0,1]^s$.
Au bout de N tirages de cette suite, la méthode de quasi Monte Carlo  nous donnerait une approximation de la valeur de $I$
par la formule :
\begin{equation*}
\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{n=1}^Nf(\xi^{(n)}) = I
\end{equation*}
Le souci de cette méthode est de ne pas pouvoir obtenir d'erreur de l'estimateur facilement. En effet l'inégalité de Koksma–Hlawka nous donne une borne de
cette erreur mais nécessite de pouvoir calculer la variation finie de la fonction $f$ ce qui n'est pas toujours possible en pratique. C'est pourquoi on introduit
la méthode de quasi Monte Carlo randomisé. Elle consiste à rajouter un décalage aléatoire à la suite $(\xi^{(n)})$. Soit $X$ une variable aléatoire
uniformément distribuée sur ${[0,1]^s}$. On s'intéresse dorénavant à la somme aléatoire :
\begin{equation*}
Z=\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\})
\end{equation*}
où \{$X$\} désigne la partie fractionnaire de $X$.
L'idée consiste à faire $K$ tirages de $Z$, ce qui revient à faire $K$ tirages de $X$ vu que la suite $(\xi^{(n)})$ est déterministe, puis de faire
une méthode de Monte Carlo classique sur ces $K$ tirages. On a donc :
\begin{equation}
\label{eqrqmc}\frac{1}{K}\sum_{k=1}^KZ_{k}=\frac{1}{K}\sum_{k=1}^K\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X_{k}\})\overset{p.s}{\to} I
\end{equation}
Il convient maintenant d'estimer l'avantage de cette méthode par rapport à une technique de Monte-Carlo usuelle, en terme de réduction de variance.
On voit dans \eqref{eqrqmc} qu'il nous coûte $NK$ évaluations de la fonction $f$ pour calculer notre estimation de $I$.
Donc la variance pour la méthode de Monte-Carlo de base serait :
\begin{equation*}
\mathrm{Var}(\frac{1}{NK}\sum_{n=1}^{NK}f(X_n))=\frac{1}{NK}\mathrm{Var}(f(X))
\end{equation*}
Il nous faut la comparer avec :
\begin{equation*}
\mathrm{Var}(\frac{1}{K}\sum_{k=1}^K\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X_{k}\}))=\frac{1}{K}\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))
\end{equation*}
On a donc réduction de variance si et seulement si :
\begin{equation*}
\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))<\frac{1}{N}\mathrm{Var}(f(X))
\end{equation*}
Or dans \cite{Tuffin:1997:VRA:268403.268419} il est montré que pour $f$ une fonction à variations bornées on a :
\begin{equation}
\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))=O(N^{-2}(\log N)^{2s})
\end{equation}
Pour cette classe de fonctions on a une réduction de variance par un facteur de $N / (\log N)^{2s}$.

\subsection{Calcul de vitesse et d'intervalle de confiance}
On s'intéresse maintenant à la convergence vers la loi normale de notre estimateur. On rappelle d'abord le cas d'une méthode de Monte-Carlo classique.
Soit $X$ une variable aléatoire, on note $\mu = \E[f(X)]$ son espérance qu'on cherche à calculer et $\sigma$ son écart type. Alors si les $(X_n)_{1\leq n\leq N}$
sont $N$ tirages de la variable $X$, $\bar{X}_N=\frac{1}{N}\sum_{n=1}^Nf(X_n)$ sera notre estimateur de $\mu$ et on aura l'intervalle de confiance
asymptotique usuel pour $\mu$ de la forme : $[\bar{X}_N-c_{\alpha}\frac{\sigma}{\sqrt{N}},\bar{X}_N+c_{\alpha}\frac{\alpha}{\sqrt{N}}]$
avec $\alpha$ le niveau désiré pour l'intervalle et $c_{\alpha}$ le quantile en $1-\frac{\alpha}{2}$ d'une loi gaussienne centrée réduite.

Dans l'article \cite{tuffin2004randomization}, Tuffin suggère d'utiliser la borne de Berry-Esseen pour le cas non asymptotique. Si on introduit les quantités
$\beta = \E[|f(X)-\E[f(X)]|^3]$, $F_n$ la fonction de répartition de $\sqrt{N}\frac{\bar{X}_N-\mu}{\sigma}$ et $\mathcal{N}$ la fonction de répartition d'une
loi normale centrée réduite, alors l'inégalité de Berry-Esseen affirme que :
\begin{equation*}
\forall x, \quad |F_n(x) - \mathcal{N}(x)| \leq C\frac{\beta}{\sigma^3 \sqrt{N}}
\end{equation*}
pour une constante universelle $C$\footnote{La meilleure valeur connue à l'heure actuelle est $C< 0.4748$.}.
Ce qui nous donne $\mathcal{N}(x) - C\frac{\beta}{\sigma ^3 \sqrt{N}} \leq F_n(x) \leq \mathcal{N}(x) + C\frac{\beta}{\sigma ^3 \sqrt{N}} $ et en
soustrayant l'encadrement équivalent de $F_n(-x)$ il vient :
\begin{equation*}
F_n(x)-F_n(-x) \geq \mathcal{N}(x) - \mathcal{N}(-x) -2C\frac{\beta}{\sigma ^3 \sqrt{N}}
\end{equation*}
On veut maintenant égaliser cette borne inférieure avec $1- \alpha$ pour trouver les bornes de l'intervalle de confiance.
Sachant que $\mathcal{N}(x) = 1 - \mathcal{N}(-x)$ on a $1 - \alpha =  1 - 2\mathcal{N}(-x) - 2C\frac{\beta}{\sigma ^3 \sqrt{N}}$ et on trouve en fin de
compte :
 $x = -\mathcal{N}^{-1}(\frac{\alpha}{2} - C\frac{\beta}{\sigma ^3 \sqrt{N}})$.

On a donc l'inégalité suivante :
\begin{equation*}
\mathbb{P}\left(\mathcal{N}^{-1}(\frac{\alpha}{2} -C\frac{\beta}{\sigma ^3 \sqrt{N}}) \leq \frac{\bar{X}_N-\mu}{\sigma / \sqrt{N}} \leq
-\mathcal{N}^{-1}(\frac{\alpha}{2} - C\frac{\beta}{\sigma ^3 \sqrt{N}})\right) \geq 1- \alpha
\end{equation*}
Alors que Tuffin trouve :
\begin{equation*}
\mathbb{P}\left(\mathcal{N}^{-1}(\frac{\alpha}{2}) - C\frac{\beta}{\sigma ^3 \sqrt{N}} \leq \frac{\bar{X}_N-\mu}{\sigma / \sqrt{N}} \leq
-\mathcal{N}^{-1}(\frac{\alpha}{2}) + C\frac{\beta}{\sigma ^3 \sqrt{N}}\right) \geq 1- \alpha
\end{equation*}
ce qui me semble être une erreur.
Si on veut comparer les demi-tailles des intervalles de confiance, on voit que l'intervalle de Tuffin rajoute $\frac{C\beta}{\sigma^2N}$ à
l'intervalle habituel asymptotique tandis que je trouve un intervalle possiblement égal à $\mathbb{R}$ tout entier pour des petites valeurs de $N$.
C'est pourquoi dans les applications ultérieures, j'utilise l'intervalle classique $[\bar{X}_N-c_{\alpha}\frac{\sigma}{\sqrt{N}},
\bar{X}_N+c_{\alpha}\frac{\alpha}{\sqrt{N}}]$.

\section{Échantillonnage stratifié}
\subsection{Présentation mathématique}
Soit $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ et $X$ un vecteur aléatoire à $d$ dimensions.
On cherche à calculer $\mu=\E(f(X))$.
On suppose maintenant qu'on divise $\mathbb{R}^{d}$ en $I$ strates appelées $(A_{i})_{1\le i \le I}$ et de telle sorte que
$\forall i \in {1,...,i}$ on ait $\mathbb{P}(X\in A_{i}) > 0$ . On notera ces probabilités $p_{i}$.
Si on note maintenant $(X_{i})_{1\le i \le I}$ l'ensemble de variables aléatoires suivant la loi conditionnelle de $X$ sachant $X\in A_{i}$,
on a clairement $\E(f(X)) = \sum_{i=1}^{I}p_{i}\E(f(X_{i}))$.
Or si on sait simuler chacune des variables $(X_{i})$, on peut estimer chacune des espérances $\E(f(X_{i}))$ par un estimateur de Monte-Carlo usuel
$\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}f(X_{i}^{j})$ où $N_{i}$ est le nombre de tirages de la variable $X_{i}$. Et donc il vient naturellement :
\begin{equation*}
\mu_{strat}=\sum_{i=1}^{I}p_{i}\sum_{j=1}^{N_{i}}f(X_{i}^{j})=
\frac{1}{N}\sum_{i=1}^{I}\frac{p_{i}}{q_{i}}\sum_{j=1}^{q_{i}N}f(X_{i}^{j}) \overset{p.s}{\to} \mu
\end{equation*}
où $N$ désigne le nombre total de tirages, c'est à dire $N=\sum_{i=1}^{I}N_{i}$, et $q_{i}=N_i / N$ la proportion de tirages dans chacune
des strates.

On veut maintenant regarder dans quels cas cette méthode peut être utile pour réduire la variance.
Si on note $\sigma_{i} = \mathrm{Var}(f(X_{i})) = \mathrm{Var}(f(X)|X\in A_{i})$, la variance de notre estimateur est :
\begin{equation*}
\mathrm{Var}(\mu_{strat})=\sum_{i=1}^{I}\frac{p_{i}^{2}\sigma_{i}^2}{N_{i}}=
\frac{1}{N}\sum_{i=1}^{I}\frac{p_{i}^{2}\sigma_{i}^2}{q_{i}}
\end{equation*}
Tandis qu'avec un estimateur usuel de Monte-Carlo pour le même nombre $N$ total de tirages, c'est à dire $\frac{1}{N}\sum_{j=1}^{N}f(X^{j})$ noté $\mu_{mc}$, on aurait :
\begin{equation*}
\mathrm{Var}(\mu_{mc}) =
\frac{1}{N}\left(\sum_{i=1}^{I}p_{i}(\sigma_{i}^2 + \E^{2}[f(X_{i})])-(\sum_{i=1}^{I}p_{i}\E[f(X_{i})])^{2}\right)
 \geq \frac{1}{N}\sum_{i=1}^{I}p_i\sigma_i^2
\end{equation*}
Pour montrer la dernière inégalité, on applique l'inégalité de Jensen à $\E[f(X_{i})]$ avec la suite des $p_i$ comme probabilité discrète (on rappelle que
$\sum_{i=1}^{I}p_i = 1$) ce qui nous donne :
\begin{equation*}
\frac{1}{N}\sum_{i=1}^{I}p_i\,\E^{2}[f(X_{i})] \geq \frac{1}{N}(\sum_{i=1}^{I}p_i\,\E[f(X_{i})])^2
\end{equation*}
Si on pose $q_i = p_i$, $\mathrm{Var}(\mu_{strat})=\frac{1}{N}\sum_{i=1}^{I}p_i\sigma_i^2$ qui est la borne inférieure trouvée pour $\mathrm{Var}(\mu_{mc})$.
Mais on peut trouver un choix de $q_i$ plus optimal en terme de réduction de variance.
En effet, toujours par l'inégalité de Jensen avec cette fois la suite des $q_i$ vue comme une probabilité discrète, on a :

\begin{multline}
\label{eq:borninf}
\mathrm{Var}(\mu_{strat})=\frac{1}{N}\sum_{i=1}^{I}\frac{p_i^2\sigma_i^2}{q_i}=\frac{1}{N}\sum_{i=1}^I(\frac{p_i\sigma_i}{q_i})^2q_i\\
\geq\frac{1}{N}(\sum_{i=1}^I\frac{p_i\sigma_i}{q_i}q_i)^2=\frac{1}{N}(\sum_{i=1}^Ip_i\sigma_i)^2=\frac{1}{N}\sigma_\star^2
\end{multline}

Et si
\begin{equation}
\label{eqqi}
\forall i \in \{1,\ldots,I\}, \quad q_i=\frac{p_i\sigma_i}{\sum_{j=1}^Ip_j\sigma_j}
\end{equation}
$\textrm{Var}(\mu_{strat})$ atteindra la borne inférieure $\sigma_\star^2/N$. Maintenant qu'on a décrit le principe de la méthode de l'échantillonnage stratifié,
et qu'on a trouvé quelle est la proportion optimale de tirages à faire par strates, intéressons nous à l'algorithme pour implémenter ce modèle.

\subsection{Description de l'algorithme}
On a vu dans la section précédente que les $q_i$ optimaux dépendent de la suite $(\sigma_i)_{1\leq i\leq I}$ qui sont les variances conditionnelles de $X$
sachant $X\in A_{i}$. Mais si ces variances sont inconnues il nous faut les estimer. On ne peut trouver la répartition optimale de tirages à faire par strates
du premier coup. On utilisera donc un algorithme adaptatif. On commence par tirer un nombre $N_1$ assez restreint de variables $(X_{i})_{1\leq i \leq I}$
en choisissant une répartition déterministe, par exemple l'allocation proportionelle $q_i=p_i$ ce qui nous donnera un premier estimateur de  $\mu=\E(f(X))$ mais
surtout les premiers estimateurs de $\sigma_i$. Ensuite on pourra tirer un nombre $N_2 > N_1$ de $(X_i)$ avec de nouveaux $q_i$ donnés par \eqref{eqqi}
ce qui nous donnera un nouvel estimateur de $\mu$ (en utilisant les $N_1 + N_2$ premières valeurs) puis de nouveaux $q_i$ plus précis, et ainsi de suite.

Essayons maintenant de décrire plus précisément comment calculer la proportion de tirages à allouer par strates. En effet, en pratique il y a deux contraintes
en plus de la formule qui nous donne la suite des $q_i$. Tout d'abord le nombre de tirages doit évidemment être entier, et de plus il doit y avoir au moins
un tirage par strate, sinon on ne peut garantir la convergence des variances conditionnelles empiriques $\hat{\sigma}_i$ vers la vraie valeur.
Notons $I$ le nombre de strates, $\hat{\sigma}_i^k$ l'estimateur de la variance de $X$ sachant que $X \in A_i$ à l'étape $k$ , $N_i^k$ le nombre total de tirages effectués de l'étape 0 à l'étape $k$ dans la $i^{\textrm{ème}}$ strate et $M_i^k$ le nombre de tirages à effectuer à
 l'étape $k$ dans la strate $i$ qui est donc $M_i^k = N_i^k - N_i^{k-1}$.

\paragraph{Étape 1 :}
\begin{equation*}
\forall i \in I, \quad M_i^1 = 1+p_i(N_1-I)
\end{equation*}
C'est l'allocation proportionelle modifié de 1, pour garantir au moins un tirage par strate. Si les $M_i^1$ ne sont pas entiers, voir la méthode ci-dessous pour les arrondir.
\paragraph{Étape k :\\}
On commence d'abord par calculer les $\hat{\sigma}_i^{k-1}$ avec les tirages de l'étape précédente.
\begin{equation*}
\forall i \in I \quad \hat{\sigma}_i^{k-1} = \sqrt{\frac{1}{N_i^{k-1}}(\sum_{j=1}^{N_i^{k-1}}(f(X_i^j))^2-(\frac{1}{N_i^{k-1}}\sum_{j=1}^{N_i^{k-1}}(f(X_i^j))^2)}
\end{equation*}
Ensuite on peut calculer les $M_i^k$. On sait qu'il nous faut au moins un tirage par strate on peut donc écrire $M_i^k$ de sous la forme  $m_i^k + 1$ avec
$m_i^k \in \mathbb{N}$ et comme on a $M_i^k = N_i^k - N_i^{k-1}$ il vient $\sum_{i=1}^Im_i^k=N^k-N^{k-1}-I$. Les auteurs de l'article proposent deux méthodes
pour calculer les $m_i^k$, je n'en ai utilisé qu'une seule que je vais présenter ici.

D'après \eqref{eqqi}, on connaît la proportion optimale à faire par strates en fonction des $\hat{\sigma}_i^{k-1}$. Sachant de plus que $\sum_{i=1}^Iq_i=1$,
on voit qu'un choix de $m_i^k$ idéal serait :
\begin{equation*}
\forall i \in I, \quad (m_i^k)^* = \frac{p_i\hat{\sigma}_i^{k-1}}{\sum_{j=1}^Ip_j\hat{\sigma}_j^{k-1}}(N^k-N^{k-1}-I)
\end{equation*}
Pour que les $m_i^k$ soient les entiers au plus près de leurs valeurs optimales, tout en assurant que leur somme ne dépasse pas le nombre total de tirages à faire
à l'étape $k$, Étoré et Jourdain proposent la méthode suivante. On pose d'abord $m_1^k =\floor{(m_1^k)^*}$. Puis pour chaque $i>1$ on pose :
\begin{equation*}
m_i^k=\floor{(m_1^k)^*+\ldots+(m_i^k)^*}-\floor{(m_1^k)^*+\ldots+(m_{i-1}^k)^*}
\end{equation*}
Ce choix des $m_i^k$ assure que notre nombre de tirages par strates sera à une distance d'au plus $1$ de l'allocation optimale.

Étoré et Jourdain prouvent enfin le résultat suivant qui est un théorème centrale limite pour la loi limite de leur estimateur $\mu_{strat}$.
Si $\E[f^2(X)]<\infty$ et si $k/N ^k\to0$ quand $k\to \infty$ alors :
\begin{equation*}
\sqrt{N^k}(\mu_{strat}^k-\mu)\underset{k\to\infty}{\overset{d}{\longrightarrow}}\mathcal{N}(0,\sigma^2_*)
\end{equation*}
avec $\sigma_*=\sum_{i=1}^Ip_i\sigma_i$, ce qui atteint la borne inférieure de $\textrm{Var}(\mu_{strat})$ trouvée en \eqref{eq:borninf}.

\section{Exemples}
\subsection{Calcul de l'espérance d'une loi normale}
Dans ce premier exemple assez simple, on va chercher à calculer $ \mu = \E[X]$ où $X$ suit une loi normale centrée réduite.
\subsubsection{Échantillonnage stratifié}
Pour cet exemple, les strates $A_i$ seront au nombre de $10$ et
\[\forall i \in{\{1,\ldots,10\}}\quad A_i =]y_{i-1}, y_i]
\]
où $y_i$ désigne
le quantile en $i/10$ d'une loi normale, en prenant pour convention $y_0 = - \infty $ et $y_{10} = + \infty $. Ainsi la suite des probabilités $(p_i)_{1\leq i\leq 10}$
est constante à $1/10$. Dans l'article, Étoré et Jourdain ont effectué quatre étapes de leur algorithme avec $N_1 = 300, N_2 = 1300, N_3 = 11300$
et $N_4 = 31300$, j'ai donc pris les mêmes nombres pour ma propre implémentation.
\subsubsection{Quasi Monte Carlo randomisé}
On cherche toujours à calculer $\mu$ mais cette fois avec la méthode de quasi Monte Carlo randomisé dont je rappelle l'estimateur :
\begin{equation*}
\frac{1}{I}\sum_{i=1}^I\frac{1}{N}\sum_{n=1}^N(\{\xi^{(n)}+X_{k}\})
\end{equation*}
Pour ce faire, j'ai utilisé la suite de Sobol à une dimension pour la suite $\xi^{(n)}$ et j'ai fixé $I$ constant à 100 et calculé l'estimateur pour
différents $N \in  \{3; 13; 113; 313\}$ afin que le produit $NI$ nous donne les mêmes nombres de tirages des variables que dans le cas
de l'échantillonnage stratifié. Ainsi on peut comparer les deux méthodes et obtenir le tableau suivant :

\begin{center}
\input{table1.tex}
\end{center}

$\textrm{IC}_{strat}$ et $\textrm{IC}_{rqmc}$ désignent les demi-tailles des intervalles de confiance pour {\mu} obtenues avec les deux méthodes. En comparant ces tailles
on peut constater la plus grande efficacité de la méthode de quasi Monte-Carlo randomisé.

\subsection{Modélisation du prix d'une option asiatique}
Cet exemple est tiré de l'article \cite{glasserman1999asymptotically} de Glasserman, Heidelberger et Shahabuddin.
On considère l'actif $S_t$ solution de l'équation différentielle stochastique suivante :
\begin{equation*}
dS_t=VS_tdW_t+rS_tdt
\end{equation*}
où $r$ représente le rendement de l'actif sans risque, $V$ la volatilité constante de $S_t$, et $W_t$ un processus de Wiener standard. Soit $ T > 0$ la date de maturité
de l'option et $(t_m = \frac{mT}{d})_{1\leq m\leq d}$ une suite de $d$ instants entre $0$ et $T$ où le prix de $S_t$ est connu. L'option asiatique $A$ avec un
prix d'exercice $K$ a pour valeur au temps $0$ :
\begin{equation*}
p=\E[e^{-rT}(\frac{1}{d}\sum_{m=1}^dS_{t_m}-K)^+]
\end{equation*}
qu'on va chercher à déterminer.

On peut exactement simuler la suite des $(S_{t_m})_{1\leq m\leq d}$ en posant $S_{t_0}=S_0$ et :
\begin{equation*}
S_{t_m}=S_{t_{m-1}}\exp([r-\frac{1}{2}V^2](t_m-t_{m-1}) + V\sqrt{t_m-t_{m-1}}X^m) \quad\forall m \in \{1,\ldots,d\}
\end{equation*}
où les $X_m$ sont des normales centrées réduites indépendantes. On cherche donc à calculer une espérance de la forme:
\[p = \E[ g(X) \mathbf{1} _D(X)],\]
où $g$ est une fonction de $\mathbb{R}^d$ dans $\mathbb{R}$ et $D$ le domaine où $g$ est strictement positive. Dans leur article, Glasserman et al. proposent
de d'abord commencer par une méthode d'échantillonage préférentiel, ou importance sampling. Cette méthode consiste à décaler les variables $X$ pour changer
leurs moyennes, afin de donner plus de poids là où le payoff sera important et  ainsi diminuer la variance de notre estimateur par rapport à une méthode
classique de Monte-Carlo. En effet, on remarque que $\forall \mu \in{\mathbb{R}^d}$:
\begin{equation*}
p = \E [g(X+\mu)e^{-\mu'X-(\frac{1}{2})\mu'\mu}\mathbf{1} _D(X+\mu)]
\end{equation*}
Ils affirment alors par un argument heuristique, puis par une preuve plus rigoureuse en utilisant la théorie des grandes déviations que la réduction de variance sera la plus importante pour $\mu^\star$ tel que :
\begin{equation*}
\mu^\star = \argmax_{x\in D} (G(x) - \frac{1}{2}xx')
\end{equation*}
avec $G(x)=\log(g(x))$. Finalement, on estimera $p$ par un estimateur de Monte-Carlo de $\E[ f_{\mu}(X) ]$ où $X\sim \mathcal{N}(0,I_d)$ et
$f_{\mu}(x)=g(x+\mu)e^{-\mu'x-(\frac{1}{2})\mu'\mu}\mathbf{1} _D(x+\mu)$ avec $\mu = \mu^*$.

On peut ensuite utiliser une méthode d'échantillonage stratifié pour calculer cet estimateur, où les différentes strates seront données par
\[
A_i=\{X\in\mathbb{R}^d\textrm{ tels que }u'X\in [y_{i-1},y_i]\}
\]
avec $u$ un vecteur de norme 1 dans $\mathbb{R}^d$ et les $y_i$ seront les quantiles en $i/I$ d'une normale centrée réduite. Il nous faut donc pouvoir tirer la loi conditionnelle de $X$ sachant $u'X\in [a, b]$. Les auteurs de l'article
donnent d'ailleurs une méthode pour simuler ces lois conditionnelle :

Soient $(a, b)$ les bornes de notre intervalle et $u$ un vecteur de $\mathbb{R}^d$ de norme 1.
Notons $F$ la fonction de répartition d'une normale centrée réduite. On simule d'abord $V = F(a) +U(F(b)-F(a))$ avec $U$ une uniforme sur $[0,1]$.
Ensuite en calculant $Z = F^{-1}(V)$ et en simulant $Y \sim \mathcal{N}(0,I_d)$ on obtient une simulation de $X$ sachant $u'X \in [a,b]$ en posant :
\begin{equation*}
X=uZ+Y+u(u'Y)
\end{equation*}
Pour cette méthode d'échantillonage stratifié, Glasserman et ses coauteurs ont pris $u=\mu/(\mu'\mu)$ et ont choisi la répartion proportionnelle du nombre
de tirages par strates (la répartition $q_i = p_i$ donné plus haut lors de la description de la méthode de stratified sampling). Plus précisément, ils ont découpé
$\mathbb{R}$ en $I=100$ strates suivant les quantiles au $1/100^{\textrm{ème}}$ de $\mathcal{N}(0,1)$ notés $y_i$ et ont alloué $N_i = N/I$ tirages par intervalles.

Étoré et Jourdain ont donc refait le même travail avec leur algorithme adaptatif qui converge vers la répartition optimale $q_i=N_i/N$ posés dans \eqref{eqqi},
technique utile dans ce cas car on ne connait pas la suite des $\sigma_i = \mathrm{Var}(f_{\mu}(X)|u'X\in [y_{i-1},y_i])$.

\subsection{Résultats}
Maintenant qu'on a décrit ce qu'est une option asiatique et comment utiliser l'importance sampling pour diminuer au mieux
la variance, on peut comparer les méthodes d'échantillonage stratifié et de quasi Monte-Carlo randomisé.

\begin{center}
\input{table2.tex}
\end{center}
Cette première table représente les différentes valeurs des estimateurs $\mu_{strat}$ et $\mu_{rqmc}$ pour la valeur d'une option d'achat (call) ainsi que les demi-tailles des intervalles de confiance.
$d$ représente la dimension, c'est à dire le nombre d'instants ${t_m}$ entre $0$ et $T$ où le prix du sous-jacent est connu. $d$ sera aussi la dimension
de notre suite à discrépance faible dans le cas du quasi Monte-Carlo randomisé. $K$ est le prix d'exercice de l'option.
On peut constater que l'avantage d'une méthode par rapport à l'autre est assez dépendant du prix d'exercice. L'algorithme de stratified sampling étant supérieur quand le prix de l'option est plux faible.

Regardons maintenant la table équivalente dans le cas de l'option de vente (put). J'obtiens :

\begin{center}
\input{table3.tex}
\end{center}

Pour conclure, on peut constater que la méthode de quasi Monte-Carlo randomisé est — suivant les valeurs des paramètres — entre trois fois plus lente et plus rapide que la méthode d'échantillonnage par strates. En tout cas sur les deux exemples étudiés ici, les deux méthodes sont du même ordre de grandeur. Sachant que quasi Monte-Carlo randomisé est bien plus simple à implémenter, nous recommendons son usage en pratique comme méthode générale pour accélérer les algorithmes de Monte-Carlo.

\section{Choix de l'implémentation informatique}
Dans cette dernière partie je parlerai de la façon dont j'ai implémenté les deux algorithmes présentés précédemment.

\subsection{Quasi Monte-Carlo randomisé}

Pour les deux exemples, j'ai utilisé la suite de Sobol comme suite à
discrépence faible. Cependant, le code est générique et peut prendre
n'importe quelle suite en argument du template. En particulier, j'ai
testé avec la méthode Halton, mais la performance de l'estimateur
était moins bonne. Pour les petites dimensions les suites de Sobol
sont codées avec la GSL \cite{GSL} dans \texttt{\detokenize{low_discrepancy.hpp}}
obtenue sur la page du cours. Malheureusement, l'implémentation de la
GSL est limitée à 40 dimensions. Pour les dimensions supérieures, j'ai
utilisé une implémentation qui fait partie de la librairie NLopt
\cite{nlopt} qui supporte jusqu'à 1111 dimensions. Cependant ces
fonctions ne sont pas publiques, j'ai donc copié directement les
fichiers \texttt{sobolseq.c} and \texttt{soboldata.h} dans le
répertoire source.

Les deux classes principales pour les exemples de quasi Monte-Carlo sont \texttt{\detokenize{monte_carlo}}
et \texttt{\detokenize{quasi_mean}}. J'ai essayé d'abstraire leurs arguments. \texttt{\detokenize{quasi_mean}} a pour template
une variable aléatoire (possiblement quasi aléatoire comme \texttt{\detokenize{quasi_gaussian}}) et une fonction. Cette classe renverra
un tirage de la variable aléatoire $\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X_{k}\})$  qui sera l'argument du template de la
classe \texttt{\detokenize{monte_carlo}}.

\subsubsection{Échantillonnage stratifié}

L'algorithme de l'achantillonnage stratifié adaptatif est principalement implémenté dans la méthode \texttt{update}
de \texttt{\detokenize{stratified_sampling}} qui a un type abstrait. Pour l'exemple de l'option asiatique on doit d'abord résoudre
le problème d'optimisation donné dans la formule 3.2 de l'article d'Etoré-Jourdain, l'implémentation se trouve dans \texttt{opti.cpp}.
J'ai utilisé un algorithme de la librairie NLopt appelé \texttt{COBYLA}\footnote{\url{http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms\#COBYLA_.28Constrained_Optimization_BY_Linear_Approximations.29}}. La classe \texttt{\detokenize{exponential_tilt}} qui est dérivée
de la structure \texttt{\detokenize{unary_function}} est utilisée pour le décalage de nos variables aléatoires décrit plus haut.

\subsection{Option}
Enfin j'ai créé une structure \texttt{\detokenize{asian_option}} qui, en arguments du constructeur, prend les différents
paramètres de l'option, ainsi qu'un booléen qui indique si nous sommes dans le cas d
 un put ou d'un call. L'operator \texttt{\detokenize{()}} nous renverra une évaluation du pay-off de l'option.
Elle est dérivée de \texttt{\detokenize{unary_function}} et peut donc s'utiliser directement dans nos exemples comme argument des templates de
\texttt{\detokenize{quasi_mean}} ou de \texttt{\detokenize{exponential_tilt}}. Ces structures se combinent élégamment entre elles. Par exemple, pour faire un tirage de quasi Mont-Carlo avec $I$ tirages aléatoires, $N$ valeurs de la suite, à discrépence faible, et $d$ la dimension du problème, il suffit d'appeler : \texttt{\detokenize{monte_carlo(I, quasi_mean<asian_option, sobol> (N, d, A)}}

\printbibliography
\end{document}