{
"cells": [
{
"cell_type": "markdown",
"id": "5f78b14c-8c80-4f65-8e80-5f71f0e802cf",
"metadata": {},
"source": [
"# Options basket, pytorch et réduction de variance"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "2997c1ad",
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import time\n",
"import os\n",
"import pickle\n",
"import numpy as np \n",
"import scipy.stats as sps \n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns \n",
"sns.set_theme()\n",
"import pandas as pd\n",
"from tqdm import tqdm\n",
"\n",
"from numpy.random import default_rng, SeedSequence\n",
"sq = SeedSequence()\n",
"seed = sq.entropy # on sauve la graine pour reproduire les résultats\n",
"rng = default_rng(sq)"
]
},
{
"cell_type": "markdown",
"id": "a21391bb-67a2-441d-a36a-66add60ca264",
"metadata": {},
"source": [
"## Retour sur Monte Carlo avec erreur préscrite\n",
"\n",
"Pour un échantillon $(X_i)_{1 \\leq i \\leq n}$ de taille $n$ l’estimateur\n",
"Monte Carlo se contruit à partir de l’estimation de la moyenne empirique\n",
"$\\hat m_n$ et de l’estimation de la variance empirique $\\hat \\sigma^2_n$\n",
"définies par $$\n",
" \\hat m_n = \\frac{1}{n} \\sum_{i=1}^n X_i \\quad \\text{et} \\quad \n",
" \\hat \\sigma_n^2 = \\frac{n}{n-1} \\big( \\frac{1}{n} \\sum_{i=1}^n X_i^2 - \\hat m_n^2 \\big) \n",
"$$ La taille de l’intervalle de confiance de niveau $\\alpha$"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "a5cfa0ca",
"metadata": {},
"outputs": [],
"source": [
"def monte_carlo_adaptive(sampling_function, epsilon: float, \n",
" batch_size: int = 100000, \n",
" proba: float = 0.95) -> dict:\n",
" \"\"\"\n",
" Effectue une estimation Monte Carlo adaptative jusqu'à ce que l'intervalle\n",
" de confiance de la moyenne soit inférieur à epsilon.\n",
"\n",
" Args:\n",
" - sampling_function (callable): Une fonction qui produit des échantillons \n",
" de taille donnée en argument.\n",
" - epsilon (float): La taille maximum tolérée pour l'IC.\n",
" - proba (float): Niveau de confiance pour l'intervalle (défaut: 0.95).\n",
" - batch_size (int): Taille du batch pour les simulations.\n",
"\n",
" Returns:\n",
" - dict: Un dictionnaire contenant les valeurs suivantes :\n",
" - \"mean\" (float): Moyenne de l'échantillon final.\n",
" - \"var\" (float): Variance de l'échantillon final.\n",
" - \"ci_size\" (float): Taille de l'intervalle de confiance.\n",
" - \"samples_size\" (int): Nombre de tirages effectués.\n",
" - \"time (s)\" (float): Temps d'execution en secondes.\n",
" \"\"\"\n",
" alpha = 1 - proba\n",
" quantile = sps.norm.ppf(1 - alpha / 2)\n",
"\n",
" sum_, sum2_, size_ = 0, 0, 0\n",
" start = time.time()\n",
" with tqdm(total=None, desc=\"Adaptive Monte Carlo\") as pbar:\n",
" while True:\n",
" # Génération d'un nouveau batch d'échantillons\n",
" samples = sampling_function(batch_size)\n",
" sum_ += samples.sum().item()\n",
" sum2_ += (samples**2).sum().item()\n",
" size_ += batch_size\n",
"\n",
" mean = sum_ / size_\n",
" var = size_ / (size_-1) * (sum2_ / size_ - mean**2)\n",
" ci_size = 2 * quantile * math.sqrt(var / size_)\n",
" \n",
" # Mise à jour de tqdm avec la taille de l'IC actuelle\n",
" pbar.set_postfix({\"IC_size\": ci_size, \"mean\": mean, \"var\": var})\n",
" pbar.update(batch_size)\n",
"\n",
" if ci_size <= epsilon:\n",
" break\n",
" stop = time.time()\n",
" return {\n",
" \"mean\": mean,\n",
" \"var\": var,\n",
" \"ci_size\": ci_size,\n",
" \"samples_size\": size_,\n",
" \"time (s)\": stop-start,\n",
" }"
]
},
{
"cell_type": "markdown",
"id": "8b1fa86d-45d7-4612-96fa-ec0da3fa4db8",
"metadata": {},
"source": [
"## Options basket: pricing Monte Carlo\n",
"\n",
"Dans cette première partie on na fait pas de la programmation orientée\n",
"objet mais on s’en approche en rassemblant tous les paramètres dans un\n",
"dictionnaire.\n",
"\n",
"> **Note**\n",
">\n",
"> Il est conseillé d’éviter les variables globales (pour prendre de\n",
"> bonnes habitudes de programmation et réutiliser facilement son code),\n",
"> donc on va regrouper les différents paramètre dans un dictionnaire."
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "2868dd4e",
"metadata": {},
"outputs": [],
"source": [
"def model_bs(d, r, S0, sigma, correlation, T):\n",
" \"\"\"\n",
" Initialise les paramètres du modèle multidimensionnel de Black-Scholes.\n",
"\n",
" Args:\n",
" - d (int): Dimension, représentant le nombre d'actifs dans le panier.\n",
" - r (float): Taux d'intérêt sans risque.\n",
" - S0 (np.ndarray): Vecteur des prix initiaux des actifs, de taille (d,).\n",
" - sigma (np.ndarray): Vecteur des volatilités des actifs, de taille (d,).\n",
" - correlation (np.ndarray): Matrice de corrélation, de taille (d, d).\n",
" - T (float): Maturité de l'option (temps jusqu'à l'échéance).\n",
"\n",
" Returns:\n",
" - dict: Un dictionnaire contenant les paramètres du modèle.\n",
" \"\"\"\n",
" return {\n",
" \"d\": d,\n",
" \"r\": r,\n",
" \"S0\": S0,\n",
" \"sigma\": sigma,\n",
" \"mu\": r - 0.5 * sigma**2,\n",
" \"correlation\": correlation,\n",
" \"correlation_cholesky\": np.linalg.cholesky(correlation),\n",
" \"T\": T,\n",
" \"actualization\": math.exp(-r * T)\n",
" }\n",
"\n",
"d = 40\n",
"rho = 0.3\n",
"bs = model_bs(d=d, r=0.1, S0=np.full((d), 100), \n",
" sigma=np.full((d), 0.3), \n",
" correlation=np.full((d,d), rho) + (1-rho)*np.eye(d),\n",
" T=1)"
]
},
{
"cell_type": "markdown",
"id": "42375da0-5d08-4b2d-a2a6-9334188a15ef",
"metadata": {},
"source": [
"### Code `numpy`\n",
"\n",
"On considère $d \\ge 2$ actifs financiers dont la loi à l’instant $T > 0$\n",
"est modélisée par une loi log-normale c’est à dire $$\n",
" \\forall i \\in \\{1,\\dots,d\\}, \\quad\n",
" S^i_T = S^i_0 \\exp\\Bigl( \\bigl(r-\\frac{\\sigma_i^2}{2}\\bigr) T + \\sigma_i \\sqrt{T} \\tilde G_i \\Bigr)\n",
"$$ où le vecteur $(\\tilde G_1,\\dots, \\tilde G_d)$ est gaussien centré de\n",
"matrice de covariance $\\Sigma$ et les constantes $r > 0$, $\\sigma_i > 0$\n",
"sont fixées. Il s’agit d’actifs financiers $(S^i_t)_{t \\in [0,T]}$,\n",
"$1 \\le i \\le d$, modélisés par un processus de Black-Scholes\n",
"multidimensionnel. On introduit la matrice $L$ triangulaire inférieure\n",
"obtenue par la décomposition de Cholesky de la matrice\n",
"$\\Sigma = L L^\\top$.\n",
"\n",
"A l’aide de cette matrice $L$, on définit la fonction\n",
"$\\Phi:\\mathbf{R}^d \\to \\mathbf{R}^d$ telle que $$\n",
" (S^1_T, \\dots, S^d_T) = \\Phi(G_1, \\dots, G_d) \\quad \\text{ou encore} \\quad S^i_T = \\Phi_i(G_1, \\dots, G_d)\n",
"$$ où $(G_1, \\dots, G_d) \\sim \\mathcal{N}(0, I_d)$ (l’égalité précédente\n",
"est à considérer en loi)."
]
},
{
"cell_type": "code",
"execution_count": 80,
"id": "39e0f974",
"metadata": {},
"outputs": [],
"source": [
"def phi(bs, Gn):\n",
" mu_T = bs[\"mu\"] * bs[\"T\"]\n",
" sig_T = bs[\"sigma\"] * math.sqrt(bs[\"T\"]) * bs[\"correlation_cholesky\"] \n",
" ST = bs[\"S0\"] * np.exp(mu_T + np.einsum('ij,pj->pi', sig_T, Gn))\n",
" return ST"
]
},
{
"cell_type": "markdown",
"id": "a1b2f9df-be5a-439d-8fd1-c41b8b6695d4",
"metadata": {},
"source": [
"La fonction `phi` correspond au modèle financier choisi (ici\n",
"Black-Scholes) et doit être modifié si on change de modèle. Il est utile\n",
"d’écrire cette fonction comme transformation d’un vecteur\n",
"$(G_n)_{n=1,\\dots,d}$ pour faciliter l’implémentation des méthodes de\n",
"réduction de variance.\n",
"\n",
"La fonction suivante `sampling_payoffs` regroupe\n",
"\n",
"- la simulation du vecteur $(G_n)_{n=1,\\dots,d}$ (en fait un nombre\n",
" `size` de vecteurs de $\\mathbf{R}^d$)\n",
"- l’application du modèle via l’appel de `phi`\n",
"- la transformation via la fonction payoff\n",
" $g:\\mathbf{R}^d \\to \\mathbf{R}$ $$\n",
" g(S_T) = \\Bigl( \\frac{1}{d} \\sum_{i=1}^d S^i_T - K \\Big)_+\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "a42002fa",
"metadata": {},
"outputs": [],
"source": [
"def sampling_payoffs(K, size, bs, rng): \n",
" Gn = rng.standard_normal((size, bs[\"d\"]))\n",
" samples = phi(bs, Gn)\n",
" payoffs = np.maximum(np.mean(samples, axis=1)-K, 0.) \n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "markdown",
"id": "4f1c0de4-1384-409c-9316-ffe63976ac76",
"metadata": {},
"source": [
"On regroupe dans une fonction `run` l’appel de la fonction\n",
"`monte_carlo_adaptive` pour différentes valeurs de\n",
"$K \\in \\{80,90,100,110,120\\}$. La fonction prend comme argument la\n",
"fonction `sampling_payoffs` (et ses arguments) car c’est la fonction qui\n",
"sera modifiée dans la suite."
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "17a820c7",
"metadata": {},
"outputs": [],
"source": [
"def run(name, function, epsilon, batch_size, **kwargs): \n",
"\n",
" result = {} \n",
" for K in [80, 90, 100, 110, 120]:\n",
" result[K] = monte_carlo_adaptive(\n",
" lambda size: function(K, size, **kwargs), \n",
" epsilon=epsilon, \n",
" batch_size=batch_size)\n",
" result_df = pd.DataFrame(result).T\n",
" result_df[\"K\"] = result_df.index\n",
" result_df[\"method\"] = name \n",
"\n",
" return result_df"
]
},
{
"cell_type": "markdown",
"id": "a9ba1a9f-824b-4be2-bfc0-30f44ed29ce5",
"metadata": {},
"source": [
"On vérifie que le temps d’execution est proportionnel au nombre de\n",
"d’échantillons utilisé dans l’estimateur Monte Carlo.\n",
"\n",
"> **Note**\n",
">\n",
"> Avant de lancer des gros calculs, pensez à vérifier que votre\n",
"> algorithme se comporte comme attendu. Ici l’estimateur Monte Carlo est\n",
"> linéaire donc le temps d’execution doit être proportionnel au au\n",
"> nombre de simulations. Si ce n’est pas le cas il y a un problème\n",
"> d’implémentation ou de gestion de la mémoire)"
]
},
{
"cell_type": "code",
"execution_count": 83,
"id": "e54bb71d",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 44000000it [00:36, 1189714.25it/s, IC_size=0.0099, mean=27.8, var=281]\n",
"Adaptive Monte Carlo: 39000000it [00:32, 1184472.25it/s, IC_size=0.00992, mean=19.4, var=250]\n",
"Adaptive Monte Carlo: 30000000it [00:25, 1195204.81it/s, IC_size=0.00991, mean=12.3, var=192]\n",
"Adaptive Monte Carlo: 19000000it [00:15, 1198813.95it/s, IC_size=0.00996, mean=6.96, var=123]\n",
"Adaptive Monte Carlo: 11000000it [00:09, 1181993.16it/s, IC_size=0.0096, mean=3.56, var=66] \n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.766534 \n",
" 280.616704 \n",
" 0.009899 \n",
" 44000000.0 \n",
" 36.985466 \n",
" 80 \n",
" numpy \n",
" \n",
" \n",
" 90 \n",
" 19.391415 \n",
" 249.953203 \n",
" 0.009924 \n",
" 39000000.0 \n",
" 32.926979 \n",
" 90 \n",
" numpy \n",
" \n",
" \n",
" 100 \n",
" 12.264333 \n",
" 191.765775 \n",
" 0.009911 \n",
" 30000000.0 \n",
" 25.101238 \n",
" 100 \n",
" numpy \n",
" \n",
" \n",
" 110 \n",
" 6.964153 \n",
" 122.694511 \n",
" 0.009961 \n",
" 19000000.0 \n",
" 15.849783 \n",
" 110 \n",
" numpy \n",
" \n",
" \n",
" 120 \n",
" 3.555150 \n",
" 66.036318 \n",
" 0.009604 \n",
" 11000000.0 \n",
" 9.307451 \n",
" 120 \n",
" numpy \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax1 = plt.subplots(layout=\"tight\")\n",
"\n",
"color = 'C0'\n",
"ax1.set_xlabel(r'strike $K$')\n",
"ax1.set_ylabel('samples size', color=color)\n",
"ax1.plot(result_numpy[\"K\"], result_numpy[\"samples_size\"], \n",
" marker=\"o\", ls=':', lw=2, color=color)\n",
"ax1.tick_params(axis='y', labelcolor=color)\n",
"\n",
"ax2 = ax1.twinx() \n",
"\n",
"color = 'C1'\n",
"ax2.set_ylabel('time (s)', color=color) \n",
"ax2.plot(result_numpy[\"K\"], result_numpy[\"time (s)\"], \n",
" marker=\"o\", ls=':', lw=2, color=color)\n",
"ax2.tick_params(axis='y', labelcolor=color)\n",
"ax2.grid(False);"
]
},
{
"cell_type": "markdown",
"id": "a4387c9c-430c-4321-b638-8e6d24f6cb3a",
"metadata": {},
"source": [
"### Passage du code en `pytorch`\n",
"\n",
"PyTorch est une bibliothèque open-source de calcul scientifique et\n",
"d’apprentissage automatique largement utilisée dans la recherche et\n",
"l’industrie. Elle est particulièrement conçue pour le développement et\n",
"l’entraînement de modèles d’apprentissage profond, mais elle offre\n",
"également des outils puissants pour la simulation et les probabilités\n",
"numériques.\n",
"\n",
"- PyTorch repose sur la manipulation de tenseurs, des structures\n",
" multidimensionnelles similaires à des matrices, qui permettent des\n",
" opérations rapides et efficaces sur de grandes quantités de données.\n",
"- PyTorch peut exécuter des calculs sur GPU, offrant ainsi des\n",
" performances considérablement améliorées pour les tâches exigeantes\n",
" en ressources.\n",
"- Différentiation automatique (Autograd): L’un des points forts de\n",
" PyTorch est son système de différenciation automatique, qui\n",
" simplifie le calcul des gradients pour optimiser des fonctions\n",
" complexes. Cela est essentiel pour l’entraînement des modèles\n",
" d’apprentissage profond.\n",
"\n",
"Pour passer d’un code NumPy à un code PyTorch il faut remplacer les\n",
"objets `np.ndarray` par des objets `torch.tensor`. Une nouveauté est la\n",
"notion de `device`.\n",
"\n",
"> **Important**\n",
">\n",
"> En PyTorch, un device désigne l’endroit où les calculs et les données\n",
"> sont stockés et exécutés. En pratique il y a essentiellement 2 types\n",
"> de device:\n",
">\n",
"> - CPU (Central Processing Unit) : l’unité de traitement standard de\n",
"> l’ordinateur.\n",
"> - GPU (Graphics Processing Unit) : une unité de traitement\n",
"> spécialisée, capable d’exécuter des opérations massivement\n",
"> parallèles.\n",
"\n",
"Le type de données stockées dans le `torch.tensor` doit être\n",
"**compatible avec le device utilisé**. Dans la suite on utilise le type\n",
"`torch.float32` pour que le code suivant puisse être executé sur un GPU\n",
"Nvidia d’il y a quelques années ou sur le GPU intégré de la puce ARM\n",
"d’apple."
]
},
{
"cell_type": "code",
"execution_count": 85,
"id": "e3a9ded1",
"metadata": {},
"outputs": [],
"source": [
"import torch\n",
"dtype = torch.float32\n",
"device = torch.device(\"cpu\")\n",
"rng_torch = torch.Generator(device=device)\n",
"\n",
"bs_torch = { \n",
" key: torch.tensor(item, dtype=dtype, device=device) \n",
" if type(item) is np.ndarray else item \n",
" for key, item in bs.items() \n",
"}"
]
},
{
"cell_type": "markdown",
"id": "78906369-6991-413d-90e1-99dbe8f5213a",
"metadata": {},
"source": [
"Il faut remplacer la fonction `phi` par la version `pytorch` tout\n",
"simplement en remplaçant les appels `np.exp` et `np.einsum` par les\n",
"versions `torch.exp` et `torch.einsum`."
]
},
{
"cell_type": "code",
"execution_count": 86,
"id": "1156dc67",
"metadata": {},
"outputs": [],
"source": [
"def phi_torch(bs, Gn):\n",
" mu_T = bs[\"mu\"] * bs[\"T\"]\n",
" sig_T = bs[\"sigma\"] * math.sqrt(bs[\"T\"]) * bs[\"correlation_cholesky\"] \n",
" ST = bs[\"S0\"] * torch.exp(mu_T + torch.einsum('ij,pj->pi', sig_T, Gn)) \n",
" return ST"
]
},
{
"cell_type": "markdown",
"id": "d3c85c9c-f7b0-433f-9605-44c0341e9e81",
"metadata": {},
"source": [
"De même on modifie la fonction `sampling_payoffs`."
]
},
{
"cell_type": "code",
"execution_count": 87,
"id": "0bee5beb",
"metadata": {},
"outputs": [],
"source": [
"def sampling_payoffs_torch(K, size, bs, device, rng): \n",
" Gn = torch.randn((size, bs[\"d\"]), dtype=dtype, \n",
" device=device, generator=rng)\n",
" samples = phi_torch(bs, Gn)\n",
" payoffs = torch.nn.functional.relu(torch.mean(samples, axis=1) - K)\n",
" #payoffs = torch.maximum(torch.mean(samples, axis=1) - K, \n",
" # torch.full((size,), 0.0, device=device)) \n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "markdown",
"id": "aa80f0c0-06cc-4833-a50c-1e5619ee3f47",
"metadata": {},
"source": [
"> **Warning**\n",
">\n",
"> Dans certains cas, la fonction `pytorch` a le même nom que la fonction\n",
"> `numpy` mais les appels autorisés peuvent être différents par exemple\n",
"> la fonction `np.maximum` peut comparer un objet `np.ndarray` et un\n",
"> scalaire alors que la version pytorch `torch.maximum` fair la\n",
"> comparaison uniquement entre 2 objets `torch.tensor`.\n",
">\n",
"> Une alternative ici est d’utiliser la fonction `ReLU` (qui est\n",
"> simplement la fonction partie positive) dans `torch.nn.functional`."
]
},
{
"cell_type": "code",
"execution_count": 88,
"id": "5819d660",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 44000000it [00:06, 6375267.95it/s, IC_size=0.0099, mean=27.8, var=281]\n",
"Adaptive Monte Carlo: 39000000it [00:06, 6327924.80it/s, IC_size=0.00993, mean=19.4, var=250]\n",
"Adaptive Monte Carlo: 30000000it [00:04, 6461682.90it/s, IC_size=0.00991, mean=12.3, var=192]\n",
"Adaptive Monte Carlo: 19000000it [00:02, 6461918.87it/s, IC_size=0.00996, mean=6.96, var=123]\n",
"Adaptive Monte Carlo: 11000000it [00:01, 6525178.40it/s, IC_size=0.0096, mean=3.55, var=65.9]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.769639 \n",
" 280.706259 \n",
" 0.009901 \n",
" 44000000.0 \n",
" 6.903001 \n",
" 80 \n",
" torch_cpu \n",
" \n",
" \n",
" 90 \n",
" 19.399929 \n",
" 250.090608 \n",
" 0.009926 \n",
" 39000000.0 \n",
" 6.164299 \n",
" 90 \n",
" torch_cpu \n",
" \n",
" \n",
" 100 \n",
" 12.259754 \n",
" 191.608462 \n",
" 0.009907 \n",
" 30000000.0 \n",
" 4.644058 \n",
" 100 \n",
" torch_cpu \n",
" \n",
" \n",
" 110 \n",
" 6.963109 \n",
" 122.697598 \n",
" 0.009961 \n",
" 19000000.0 \n",
" 2.941267 \n",
" 110 \n",
" torch_cpu \n",
" \n",
" \n",
" 120 \n",
" 3.552353 \n",
" 65.912922 \n",
" 0.009595 \n",
" 11000000.0 \n",
" 1.687111 \n",
" 120 \n",
" torch_cpu \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.767202 \n",
" 280.529912 \n",
" 0.009898 \n",
" 44000000.0 \n",
" 0.726342 \n",
" 80 \n",
" torch_cuda \n",
" \n",
" \n",
" 90 \n",
" 19.395393 \n",
" 250.044280 \n",
" 0.009926 \n",
" 39000000.0 \n",
" 0.340875 \n",
" 90 \n",
" torch_cuda \n",
" \n",
" \n",
" 100 \n",
" 12.265392 \n",
" 191.688457 \n",
" 0.009909 \n",
" 30000000.0 \n",
" 0.232198 \n",
" 100 \n",
" torch_cuda \n",
" \n",
" \n",
" 110 \n",
" 6.965617 \n",
" 122.798718 \n",
" 0.009965 \n",
" 19000000.0 \n",
" 0.137718 \n",
" 110 \n",
" torch_cuda \n",
" \n",
" \n",
" 120 \n",
" 3.552220 \n",
" 66.027713 \n",
" 0.009604 \n",
" 11000000.0 \n",
" 0.080104 \n",
" 120 \n",
" torch_cuda \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"result = pd.concat([ result_numpy, result_torch_cpu, result_torch_mps ])\n",
"ax = sns.barplot(data=result, x=\"K\", y=\"time (s)\", hue=\"method\")\n",
"ax.set_title(fr\"Temps d'execution en fonction du strike $K$ pour \" + \n",
" fr\"atteindre une précision de $\\epsilon$ = {epsilon}\");"
]
},
{
"cell_type": "markdown",
"id": "b2665367-327c-4e34-8fe4-acd388534e21",
"metadata": {},
"source": [
"## Méthodes de réduction de variance\n",
"\n",
"On rappelle que le modèle est codé à travers la fonction `phi`. Cette\n",
"fonction ne sera pas modifiée dans la suite. Le payoff est donné par une\n",
"fonction $g:\\mathbf{R}^d \\to \\mathbf{R}_+$ et le code est fait dans la\n",
"fonction `sampling_payoffs`. On va implémenter les méthodes des\n",
"réduction de variance vues en cours en modifiant la fonction\n",
"`sampling_payoffs`.\n",
"\n",
"### Variables antithétiques\n",
"\n",
"Pour utiliser les variables antithétiques on utilise le fait que\n",
"$G \\sim \\mathcal{N}(0, \\operatorname{Id}_d)$ et $-G$ ont même loi et\n",
"donc que $$\n",
" \\mathbf{E} \\big[ (g \\circ \\Phi)(G) \\big] = \\frac{1}{2} \n",
" \\mathbf{E} \\big[ (g \\circ \\Phi)(G) + (g \\circ \\Phi)(-G) \\big]. \n",
"$$ Cela suggère d’utiliser la représentation de droite pour écrire un\n",
"estimateur Monte Carlo de la forme $$\n",
" \\tilde I_n = \\frac{1}{2 n} \\sum_{i=1}^n \\big( g \\circ \\Phi(G_i) + g \\circ \\Phi(-G_i) \\big),\n",
"$$ où $(G_i)_{i \\ge 1}$ est une suite *i.i.d.* de loi\n",
"$\\mathcal{N}(0, \\operatorname{Id}_d)$.\n",
"\n",
"Voici le code modifié pour implémenter les variables antithétiques."
]
},
{
"cell_type": "code",
"execution_count": 93,
"id": "28b183d2",
"metadata": {},
"outputs": [],
"source": [
"def sampling_payoffs_antithetic(K, size, bs, rng): \n",
" Gn = rng.standard_normal((size, bs[\"d\"]))\n",
" payoffs_1 = np.maximum(np.mean(phi(bs, Gn), axis=1) - K, 0) \n",
" payoffs_2 = np.maximum(np.mean(phi(bs, -Gn), axis=1) - K, 0) \n",
" payoffs = 0.5 * (payoffs_1 + payoffs_2)\n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "markdown",
"id": "d47e15a7-49ed-419a-85c6-18288b0fe2d2",
"metadata": {},
"source": [
"> **Warning**\n",
">\n",
"> Dans la suite comme le nombre d’échantillons va diminuer (car on\n",
"> utilise des méthodes de réduction de variance) on change le\n",
"> `batch_size` dans le Monte Carlo adaptatif: on passe de $1\\,000\\,000$\n",
"> à $100\\,000$."
]
},
{
"cell_type": "code",
"execution_count": 94,
"id": "9c76b825",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 1200000it [00:01, 696174.60it/s, IC_size=0.00984, mean=27.8, var=7.57]\n",
"Adaptive Monte Carlo: 2500000it [00:03, 687148.39it/s, IC_size=0.00983, mean=19.4, var=15.7]\n",
"Adaptive Monte Carlo: 4600000it [00:06, 684916.69it/s, IC_size=0.00992, mean=12.3, var=29.5]\n",
"Adaptive Monte Carlo: 5700000it [00:08, 681781.83it/s, IC_size=0.01, mean=6.96, var=37.1] \n",
"Adaptive Monte Carlo: 4100000it [00:06, 680311.51it/s, IC_size=0.00999, mean=3.55, var=26.6]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.763688 \n",
" 7.565365 \n",
" 0.009842 \n",
" 1200000.0 \n",
" 1.724977 \n",
" 80 \n",
" antithetic \n",
" \n",
" \n",
" 90 \n",
" 19.399543 \n",
" 15.724141 \n",
" 0.009831 \n",
" 2500000.0 \n",
" 3.639759 \n",
" 90 \n",
" antithetic \n",
" \n",
" \n",
" 100 \n",
" 12.264132 \n",
" 29.457681 \n",
" 0.009920 \n",
" 4600000.0 \n",
" 6.717031 \n",
" 100 \n",
" antithetic \n",
" \n",
" \n",
" 110 \n",
" 6.962649 \n",
" 37.074154 \n",
" 0.009997 \n",
" 5700000.0 \n",
" 8.361540 \n",
" 110 \n",
" antithetic \n",
" \n",
" \n",
" 120 \n",
" 3.548517 \n",
" 26.633407 \n",
" 0.009991 \n",
" 4100000.0 \n",
" 6.027799 \n",
" 120 \n",
" antithetic \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"result = pd.concat([ result_numpy, result_antithetic ])\n",
"fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6.4,6.4), \n",
" sharex=True, layout=\"tight\")\n",
"sns.barplot(data=result, x=\"K\", y=\"var\", hue=\"method\", ax=ax1)\n",
"sns.barplot(data=result, x=\"K\", y=\"time (s)\", hue=\"method\", ax=ax2)\n",
"fig.suptitle(fr\"Variance et temps d'execution en fonction du strike $K$, \"\n",
" +fr\"précision demandée $\\epsilon$ = {epsilon}\");"
]
},
{
"cell_type": "markdown",
"id": "dfc495a5-bc21-475b-ad64-e44d9232054e",
"metadata": {},
"source": [
"#### Passage en `pytorch`\n",
"\n",
"Pas de difficulté particulière ici, si suffit de modifier la fonction\n",
"`sampling_payoffs_torch` pour évaluer 2 fois $g \\circ \\Phi$ (via `relu`,\n",
"`torch.mean` et `phi_torch`)."
]
},
{
"cell_type": "code",
"execution_count": 112,
"id": "25f38a46",
"metadata": {},
"outputs": [],
"source": [
"from torch.nn.functional import relu \n",
"def sampling_payoffs_torch_antithetic(K, size, bs, device, rng): \n",
" Gn = torch.randn((size, bs[\"d\"]), dtype=dtype, \n",
" device=device, generator=rng)\n",
" #samples = phi_torch(bs, Gn)\n",
" payoffs_1 = relu(torch.mean(phi_torch(bs, Gn), axis=1) - K)\n",
" payoffs_2 = relu(torch.mean(phi_torch(bs, -Gn), axis=1) - K)\n",
" payoffs = 0.5 * (payoffs_1 + payoffs_2)\n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "code",
"execution_count": 113,
"id": "3531b1e6",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 1180000it [00:00, 9147256.66it/s, IC_size=0.00993, mean=27.8, var=7.57]\n",
"Adaptive Monte Carlo: 2420000it [00:00, 11499575.92it/s, IC_size=0.00997, mean=19.4, var=15.7]\n",
"Adaptive Monte Carlo: 4540000it [00:00, 15031852.39it/s, IC_size=0.00998, mean=12.3, var=29.4]\n",
"Adaptive Monte Carlo: 5720000it [00:00, 18833446.41it/s, IC_size=0.00998, mean=6.96, var=37.1]\n",
"Adaptive Monte Carlo: 4120000it [00:00, 19881923.06it/s, IC_size=0.00998, mean=3.55, var=26.7]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.765727 \n",
" 7.574819 \n",
" 0.009932 \n",
" 1180000.0 \n",
" 0.132232 \n",
" 80 \n",
" antithetic_cuda \n",
" \n",
" \n",
" 90 \n",
" 19.393324 \n",
" 15.652379 \n",
" 0.009969 \n",
" 2420000.0 \n",
" 0.211281 \n",
" 90 \n",
" antithetic_cuda \n",
" \n",
" \n",
" 100 \n",
" 12.263614 \n",
" 29.434453 \n",
" 0.009981 \n",
" 4540000.0 \n",
" 0.303142 \n",
" 100 \n",
" antithetic_cuda \n",
" \n",
" \n",
" 110 \n",
" 6.962103 \n",
" 37.108212 \n",
" 0.009984 \n",
" 5720000.0 \n",
" 0.304725 \n",
" 110 \n",
" antithetic_cuda \n",
" \n",
" \n",
" 120 \n",
" 3.553585 \n",
" 26.681542 \n",
" 0.009976 \n",
" 4120000.0 \n",
" 0.208037 \n",
" 120 \n",
" antithetic_cuda \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"result = pd.concat([ result_numpy, result_torch_cpu, result_torch_mps, \n",
" result_antithetic, result_antithetic_mps ])\n",
"ax = sns.barplot(data=result, x=\"K\", y=\"time (s)\", hue=\"method\")\n",
"ax.set_title(fr\"Temps d'execution en fonction du strike $K$ pour \" + \n",
" fr\"atteindre une précision de $\\epsilon$ = {epsilon}\");"
]
},
{
"cell_type": "markdown",
"id": "f597ae0f-c0af-4d07-911f-54783c1f0532",
"metadata": {},
"source": [
"### Variable de contrôle\n",
"\n",
"On rappelle que le prix d’un Call dans un modèle de Black-Scholes en\n",
"dimension 1, le prix est donnée par une formule fermée. Pour une option\n",
"Basket (en dimension $d \\ge 2$) on approche le prix par Monte Carlo mais\n",
"on peut utiliser des approximations pour trouver un produit\n",
"unidimensionnel proche du produit Basket. Ces approximations servent de\n",
"variables de contrôles: **on ne rajoute pas une erreur, on retire de la\n",
"variance**…\n",
"\n",
"On rappelle que, en posant $\\mu_i = r - \\frac{1}{2}\\sigma_i^2$, $$\n",
" X = \\biggl(\\frac{1}{d} \\sum_{i=1}^d S^i_0 e^{\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i} - K\\biggr)_+\n",
"$$ et en introduisant $a^i_0 = \\frac{S^i_0}{\\sum_{j=1}^d S^j_0}$ (t.q.\n",
"$\\sum a^i_0 = 1$) et $\\bar S_0 = \\frac{1}{d} \\sum_{i=1}^d S^i_0$ on a $$\n",
" X = \\biggl(\\bar S_0 \\sum_{i=1}^d a^i_0 e^{\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i} - K\\biggr)_+.\n",
"$$ La variable de contrôle proposée est obtenue en échangeant\n",
"l’exponentielle et la moyenne pondérée par les poids\n",
"$\\big(a^i_0\\big)_{i=1,\\dots,d}$: $$\n",
" Y = \\bigl(\\bar S_0 e^Z - K\\bigr)_+\n",
" \\quad \\text{avec} \\quad \n",
" Z = \\sum_{i=1}^d a^i_0 \\big(\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i\\big) \n",
"$$ La variable aléatoire $Z$ suit une loi gaussienne\n",
"$Z \\sim \\mathcal{N}(m T, s^2 T)$ avec $$\n",
" m = \\sum_{i=1}^d a^i_0 \\mu_i\n",
" \\quad \\text{et} \\quad\n",
" s^2 = \\sum_{i=1}^d \\Big( \\sum_{j=1}^d a^i_0 \\sigma_i L_{ij} \\Big)^2. \n",
"$$ Ainsi l’espérance de la variable de contrôle $Y$ est connue par la\n",
"formule de Black-Scholes, car elle correspond au prix d’un call de\n",
"strike $K$ d’un actif Black-Scholes de dimension 1, de valeur initiale\n",
"$\\bar S_0$, de taux $\\rho = m+\\frac{1}{2} s^2$ et de volatilité $s$ (à\n",
"un facteur d’actualisation près… attention à ça). On a donc $$\n",
" e^{-\\rho T} \\mathbf{E} \\big[ Y \\big] = P_{\\text{BS}}\\big(\\bar S_0, \\rho, s, T, K\\big),\n",
"$$ où $$\n",
" P_{\\text{BS}}\\big(x, r, \\sigma, T, K\\big) = x F_{\\mathcal{N}(0,1)}(d_1) - K e^{-r T} F_{\\mathcal{N}(0,1)}(d_2),\n",
"$$ avec $F_{\\mathcal{N}(0,1)}$ est la fonction de répartition de la loi\n",
"normale centrée réduite et la notation $$\n",
" d_1 = \\frac{1}{\\sigma \\sqrt{T}} \\Big( \\log\\big( \\frac{x}{K} \\big) \n",
" + \\big(r + \\frac{\\sigma^2}{2}\\big) T \\Big)\n",
" \\quad \\text{et} \\quad\n",
" d_2 = d_1 - \\sigma \\sqrt{T}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "3b52ec13",
"metadata": {},
"outputs": [],
"source": [
"def black_scholes_call(S0: float, K: float, T: float, \n",
" r: float, sigma: float) -> float:\n",
" \"\"\"\n",
" Calcule le prix d'un Call européen dans le modèle de Black-Scholes.\n",
"\n",
" Args:\n",
" - S0 (float): Prix actuel de l'actif sous-jacent.\n",
" - K (float): Prix d'exercice du Call.\n",
" - T (float): Temps jusqu'à l'échéance en années.\n",
" - r (float): Taux d'intérêt sans risque.\n",
" - sigma (float): Volatilité du sous-jacent.\n",
"\n",
" Returns:\n",
" - float: Prix du Call européen.\n",
" \"\"\"\n",
" # Calcul des paramètres d1 et d2\n",
" d1 = (math.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*math.sqrt(T))\n",
" d2 = d1 - sigma*math.sqrt(T)\n",
" \n",
" # Calcul du prix du Call\n",
" call_price = S0*sps.norm.cdf(d1) - K*math.exp(-r*T)*sps.norm.cdf(d2)\n",
" return call_price"
]
},
{
"cell_type": "markdown",
"id": "6279cf4b-a9eb-48d2-873f-d97fd0ce8de4",
"metadata": {},
"source": [
"> **Note**\n",
">\n",
"> On code une fonction `payoffs_control` qui renvoie des réalisation de\n",
"> la variable de contrôle centrée $\\tilde Y = Y - \\mathbf{E}[Y]$. Il est\n",
"> important de noter que\n",
">\n",
"> - l’espérance est calculée en utilisant la formule fermée obtenue\n",
"> car $Z \\sim \\mathcal{N}(m T, s^2 T)$\n",
"> - les réalisations de $Y$ (et donc de $Z$) sont obtenues par\n",
"> transformation de $G$ et donc la fonction `payoffs_control` doit\n",
"> prendre comme argument l’objet `Gn`"
]
},
{
"cell_type": "code",
"execution_count": 100,
"id": "5d9d6519",
"metadata": {},
"outputs": [],
"source": [
"def payoffs_control(bs, K, Gn): \n",
" weight = bs[\"S0\"] / bs[\"S0\"].sum()\n",
" m = (weight * bs[\"mu\"]).sum()\n",
" s2 = (((weight * bs[\"sigma\"]) @ bs[\"correlation_cholesky\"])**2).sum()\n",
" rho = m + 0.5*s2\n",
"\n",
" mu_T = bs[\"mu\"] * bs[\"T\"]\n",
" sig_T = bs[\"sigma\"] * math.sqrt(bs[\"T\"]) * bs[\"correlation_cholesky\"] \n",
" Z = np.sum(weight * (mu_T + np.einsum('ij,pj->pi', sig_T, Gn)), axis=1)\n",
" Y = np.maximum(bs[\"S0\"].mean() * np.exp(Z) - K, 0)\n",
" Y_mean = np.exp(rho*bs[\"T\"]) * \\\n",
" black_scholes_call(bs[\"S0\"].mean(), K, bs[\"T\"], \n",
" r=rho, sigma=math.sqrt(s2))\n",
" return Y - Y_mean"
]
},
{
"cell_type": "markdown",
"id": "288b7880-7a5a-4973-b684-b5a902c15bac",
"metadata": {},
"source": [
"Etant donné cette variable aléatoire\n",
"$\\tilde Y = Y - \\mathbf{E}[Y] = \\psi(G)$ centrée, l’estimateur Monte\n",
"Carlo s’écrit $$\n",
" \\bar I_n = \\frac{1}{n} \\sum_{i=1}^n \\big( (g \\circ \\Phi)(G_i) - \\psi(G_i) \\big).\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 101,
"id": "d76fab43",
"metadata": {},
"outputs": [],
"source": [
"def sampling_payoffs_control_variate(K, size, bs, rng): \n",
" Gn = rng.standard_normal((size, bs[\"d\"]))\n",
" payoffs_X = np.maximum(np.mean(phi(bs, Gn), axis=1) - K, 0) \n",
" payoffs_Y = payoffs_control(bs, K, Gn)\n",
" payoffs = payoffs_X - payoffs_Y\n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "markdown",
"id": "8a17eb2f-1e40-456f-9573-87d5d1831c4b",
"metadata": {},
"source": [
"> **Warning**\n",
">\n",
"> On réduit encore le `batch_size` et on passe à $20\\,000$."
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "3e1d3bf1",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 180000it [00:00, 793527.46it/s, IC_size=0.00961, mean=27.8, var=1.08]\n",
"Adaptive Monte Carlo: 300000it [00:00, 766177.62it/s, IC_size=0.00972, mean=19.4, var=1.85]\n",
"Adaptive Monte Carlo: 440000it [00:00, 784494.61it/s, IC_size=0.00996, mean=12.3, var=2.84]\n",
"Adaptive Monte Carlo: 500000it [00:00, 787215.67it/s, IC_size=0.00998, mean=6.96, var=3.24]\n",
"Adaptive Monte Carlo: 420000it [00:00, 796634.97it/s, IC_size=0.00998, mean=3.56, var=2.72]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.767146 \n",
" 1.080913 \n",
" 0.009606 \n",
" 180000.0 \n",
" 0.227993 \n",
" 80 \n",
" var_cont \n",
" \n",
" \n",
" 90 \n",
" 19.392709 \n",
" 1.846307 \n",
" 0.009725 \n",
" 300000.0 \n",
" 0.392570 \n",
" 90 \n",
" var_cont \n",
" \n",
" \n",
" 100 \n",
" 12.268541 \n",
" 2.838352 \n",
" 0.009956 \n",
" 440000.0 \n",
" 0.561966 \n",
" 100 \n",
" var_cont \n",
" \n",
" \n",
" 110 \n",
" 6.964850 \n",
" 3.243183 \n",
" 0.009983 \n",
" 500000.0 \n",
" 0.636577 \n",
" 110 \n",
" var_cont \n",
" \n",
" \n",
" 120 \n",
" 3.557704 \n",
" 2.724161 \n",
" 0.009983 \n",
" 420000.0 \n",
" 0.528792 \n",
" 120 \n",
" var_cont \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.768241 \n",
" 1.076860 \n",
" 0.009588 \n",
" 180000.0 \n",
" 0.019557 \n",
" 80 \n",
" var_cont_cuda \n",
" \n",
" \n",
" 90 \n",
" 19.393873 \n",
" 1.854623 \n",
" 0.009746 \n",
" 300000.0 \n",
" 0.030391 \n",
" 90 \n",
" var_cont_cuda \n",
" \n",
" \n",
" 100 \n",
" 12.262587 \n",
" 2.838412 \n",
" 0.009956 \n",
" 440000.0 \n",
" 0.043680 \n",
" 100 \n",
" var_cont_cuda \n",
" \n",
" \n",
" 110 \n",
" 6.963644 \n",
" 3.239082 \n",
" 0.009977 \n",
" 500000.0 \n",
" 0.053451 \n",
" 110 \n",
" var_cont_cuda \n",
" \n",
" \n",
" 120 \n",
" 3.558166 \n",
" 2.727846 \n",
" 0.009990 \n",
" 420000.0 \n",
" 0.041741 \n",
" 120 \n",
" var_cont_cuda \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" lambda \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.768880 \n",
" 0.667401 \n",
" 0.009244 \n",
" 120000.0 \n",
" 0.157157 \n",
" 1.038273 \n",
" 80 \n",
" var_cont_opt \n",
" \n",
" \n",
" 90 \n",
" 19.392490 \n",
" 1.085898 \n",
" 0.009628 \n",
" 180000.0 \n",
" 0.230559 \n",
" 1.062720 \n",
" 90 \n",
" var_cont_opt \n",
" \n",
" \n",
" 100 \n",
" 12.264128 \n",
" 1.513434 \n",
" 0.009844 \n",
" 240000.0 \n",
" 0.304000 \n",
" 1.090117 \n",
" 100 \n",
" var_cont_opt \n",
" \n",
" \n",
" 110 \n",
" 6.966865 \n",
" 1.636613 \n",
" 0.009835 \n",
" 260000.0 \n",
" 0.330056 \n",
" 1.128270 \n",
" 110 \n",
" var_cont_opt \n",
" \n",
" \n",
" 120 \n",
" 3.555140 \n",
" 1.345922 \n",
" 0.009696 \n",
" 220000.0 \n",
" 0.278473 \n",
" 1.170608 \n",
" 120 \n",
" var_cont_opt \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"result = pd.concat([ result_antithetic, result_antithetic_mps, \n",
" result_varcont, result_varcont_opt ])\n",
"ax = sns.barplot(data=result, x=\"K\", y=\"time (s)\", hue=\"method\")\n",
"ax.set_title(fr\"Temps d'execution en fonction du strike $K$ pour \" + \n",
" fr\"atteindre une précision de $\\epsilon$ = {epsilon}\");"
]
},
{
"cell_type": "markdown",
"id": "aff696b2-cfe2-4d99-b9af-ba29fef340b2",
"metadata": {},
"source": [
"### Combinaison des 2 méthodes\n",
"\n",
"Il est possible de combiner plusieurs méthodes de réduction de variance.\n",
"Ici on fait une variable de contrôle (optimale) sur des variables\n",
"antithétiques."
]
},
{
"cell_type": "code",
"execution_count": 109,
"id": "62176a17",
"metadata": {},
"outputs": [],
"source": [
"def sampling_payoffs_antithetic_control_variate_opt(K, size, bs, rng): \n",
" Gn = rng.standard_normal((size, bs[\"d\"]))\n",
" payoffs_X1 = np.maximum(np.mean(phi(bs, Gn), axis=1) - K, 0) \n",
" payoffs_X2 = np.maximum(np.mean(phi(bs, -Gn), axis=1) - K, 0)\n",
" payoffs_X = 0.5 * (payoffs_X1 + payoffs_X2)\n",
" payoffs_Y1 = payoffs_control(bs, K, Gn)\n",
" payoffs_Y2 = payoffs_control(bs, K, -Gn)\n",
" payoffs_Y = 0.5 * (payoffs_Y1 + payoffs_Y2)\n",
"\n",
" cntxt[\"lambda_numerator\"] += (payoffs_X * payoffs_Y).sum() \n",
" cntxt[\"lambda_denominator\"] += (payoffs_Y**2).sum() \n",
" cntxt[\"lambda\"] = cntxt[\"lambda_numerator\"] / cntxt[\"lambda_denominator\"]\n",
"\n",
" payoffs = payoffs_X - cntxt[\"lambda\"] * payoffs_Y\n",
" return bs[\"actualization\"] * payoffs"
]
},
{
"cell_type": "code",
"execution_count": 110,
"id": "afdce2a0",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Adaptive Monte Carlo: 80000it [00:00, 428014.23it/s, IC_size=0.00992, mean=27.8, var=0.512]\n",
"Adaptive Monte Carlo: 100000it [00:00, 417425.51it/s, IC_size=0.00901, mean=19.4, var=0.528]\n",
"Adaptive Monte Carlo: 100000it [00:00, 424307.82it/s, IC_size=0.00902, mean=12.3, var=0.53]\n",
"Adaptive Monte Carlo: 60000it [00:00, 413794.52it/s, IC_size=0.00884, mean=6.97, var=0.305]\n",
"Adaptive Monte Carlo: 100000it [00:00, 424733.62it/s, IC_size=0.00901, mean=3.55, var=0.528]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" var \n",
" ci_size \n",
" samples_size \n",
" time (s) \n",
" lambda \n",
" K \n",
" method \n",
" \n",
" \n",
" \n",
" \n",
" 80 \n",
" 27.766779 \n",
" 0.511902 \n",
" 0.009916 \n",
" 80000.0 \n",
" 0.188013 \n",
" 0.972516 \n",
" 80 \n",
" antithetic_var_cont_opt \n",
" \n",
" \n",
" 90 \n",
" 19.398933 \n",
" 0.528058 \n",
" 0.009008 \n",
" 100000.0 \n",
" 0.240708 \n",
" 0.931822 \n",
" 90 \n",
" antithetic_var_cont_opt \n",
" \n",
" \n",
" 100 \n",
" 12.266264 \n",
" 0.529686 \n",
" 0.009022 \n",
" 100000.0 \n",
" 0.236727 \n",
" 0.952754 \n",
" 100 \n",
" antithetic_var_cont_opt \n",
" \n",
" \n",
" 110 \n",
" 6.965330 \n",
" 0.304993 \n",
" 0.008838 \n",
" 60000.0 \n",
" 0.146040 \n",
" 1.060539 \n",
" 110 \n",
" antithetic_var_cont_opt \n",
" \n",
" \n",
" 120 \n",
" 3.552981 \n",
" 0.528052 \n",
" 0.009008 \n",
" 100000.0 \n",
" 0.236170 \n",
" 1.140471 \n",
" 120 \n",
" antithetic_var_cont_opt \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"result = pd.concat([ result_varcont, result_varcont_opt, \n",
" result_antithetic_varcont_opt ])\n",
"fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6.4,6.4), \n",
" sharex=True, layout=\"tight\")\n",
"sns.barplot(data=result, x=\"K\", y=\"var\", hue=\"method\", ax=ax1)\n",
"sns.barplot(data=result, x=\"K\", y=\"time (s)\", hue=\"method\", ax=ax2)\n",
"fig.suptitle(fr\"Variance et temps d'execution en fonction du strike $K$, \"\n",
" +fr\"précision demandée $\\epsilon$ = {epsilon}\");"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}