{ "cells": [ { "cell_type": "markdown", "id": "7a149ebe-0121-407d-ab82-22dbdae6833e", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Calcul de sensibilités" ] }, { "cell_type": "code", "execution_count": 1, "id": "1bb383e8-5124-4169-b672-77bd10f1b6d1", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:26.569687Z", "iopub.status.busy": "2025-12-20T09:23:26.569479Z", "iopub.status.idle": "2025-12-20T09:23:27.113919Z", "shell.execute_reply": "2025-12-20T09:23:27.113442Z" }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "sns.set_theme() \n", "from numpy.random import default_rng\n", "rng = default_rng()\n", "import pandas as pd" ] }, { "cell_type": "markdown", "id": "e93be7c7-b92f-44c1-9dde-c4e06aaa199c", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Problème de charge sinistre\n", "\n", "On définit la _charge sinistre totale_ (sur une période $T$) par la variable aléatoire positive\n", "\n", "$$\n", " S = \\sum_{i = 1}^N X_i\n", "$$\n", "où $N$ est une variable aléatoire à valeurs dans $\\mathbf{N}$ représentant le nombre de sinistres sur la période $T$, et pour $i \\ge 1$, $X_i$ est une variable aléatoire à valeurs dans $\\mathbf{R}_+$ représentant le coût du i-ème sinistre, avec la convention selon laquelle la somme est nulle si $N = 0$. Les $(X_i)_{i \\ge 1}$ sont supposées indépendantes et identiquement distribuées, et indépendantes de $N$ (indépendance fréquences - coûts).\n", "\n", "Une modélisation classique est de considérer \n", "\n", "- $N$ de loi de Poisson de paramètre $\\lambda > 0$, \n", "- $X_1$ de loi log-normale de paramètres $\\mu > 0$, $\\sigma^2 > 0$, c'est à dire $X_1 = \\exp(G_1)$ avec $G_1 \\sim \\mathcal{N}(\\mu, \\sigma^2)$. \n", "\n", "Le but est d'estimer la **probabilité de dépassement** c'est à dire calculer la probabilité que la charge sinistre totale dépasse un seuil $K$:\n", "\n", "$$\n", " p = \\mathbf{P}\\bigl[S > K\\bigr] \\quad \\text{pour $K$ grand}\n", "$$\n", "\n", "Dans la suite on prend $\\lambda = 10$, $\\mu = 0.1$ et $\\sigma = 0.3$ et on considère plusieurs valeurs du seuil $K$." ] }, { "cell_type": "markdown", "id": "0c00bdf6-02a1-4f80-90bd-41b25134f0ff", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "On utilisera la notation $S^{(\\lambda)}$ pour indiquer la dépendance de variable aléatoire $S = \\sum_{i=1}^{N} X_i$ en le paramètre $\\lambda > 0$ (paramètre de la loi de Poisson sous-jacente). On s'intéresse à la sensibilité de la probabilité $p$ en fonction de lambda c'est à dire \n", "$$\n", " \\frac{\\partial}{\\partial \\lambda} p(\\lambda) = \\frac{\\partial}{\\partial \\lambda} \\mathbf{P} \\bigl[ S^{\\lambda} > K \\bigl].\n", "$$" ] }, { "cell_type": "markdown", "id": "965ce11d-62ad-48be-9ebb-881ef2c6de88", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Méthode des différences finies\n", "Implémenter l'estimateur Monte Carlo basé sur les différences finies d'ordre 2\n", " \n", "$$\n", " \\frac{\\partial}{\\partial \\lambda} p(\\lambda) = \\frac{p(\\lambda + h) - p(\\lambda- h)}{2h} + \\mathcal{O}(h^2)\n", "$$\n", "Comme vu en cours, il y a plusieurs façon d'implémenter l'estimateur Monte Carlo dans ce cadre biaisé. \n", "\n", "Le **premier estimateur** naïf $J^{(1)}_{n,h}(\\lambda)$ est basé sur des réalisations indépendantes de $S^{(\\lambda+h)}$ et $S^{(\\lambda-h)}$ et n'est pas efficace: la variance explose lorsque $h$ tend vers 0. Ainsi on pose\n", "\n", "$$\n", " J^{(1)}_{n, h}(\\lambda) = \\frac{1}{2 h n} \\bigl( \\sum_{k = 1}^n \\mathbf{1}_{\\{S^{(\\lambda+h)}_k > K\\}} - \\sum_{k = 1}^n \\mathbf{1}_{\\{\\tilde S^{(\\lambda-h)}_k > K\\}} \\bigr),\n", "$$\n", "où $(S^{(\\lambda+h)}_k)_{k \\ge 1}$ et $(\\tilde S^{(\\lambda-h)}_k)_{k \\ge 1}$ sont des suites indépendantes de variables aléatoires _i.i.d._.\n", "\n", "Le **deuxième estimateur** $J^{(2)}_{n,h}(\\lambda)$ utilise des réalisations fortements corrélées de la loi de Poisson au sens suivant: on utilise la même réalisation uniforme $U$ pour constuire deux réalisations $N^{(\\lambda+h)}$ et $N^{(\\lambda-h)}$ en utilisant la méthode de l'inverse de la fonction de répartition. Dans ce deuxième estimateur, les lois log-normales sont indépendantes. On a donc\n", "\n", "$$\n", " J^{(2)}_{n, h}(\\lambda) = \\frac{1}{2 h n} \\sum_{k = 1}^n \\bigl(\\mathbf{1}_{\\{S^{(\\lambda+h)}_k > K\\}} - \\mathbf{1}_{\\{\\bar S^{(\\lambda-h)}_k > K\\}} \\bigr),\n", "$$\n", "où pour $k \\ge 1$, $S^{(\\lambda+h)}_k = \\sum_{i = 1}^{G(\\lambda+h, U_k)} X_{i,k}$ et $\\bar S^{(\\lambda-h)}_k = \\sum_{i = 1}^{G(\\lambda-h, U_k)} \\bar X_{i,k}$ avec $G(\\lambda, u)$ l'inverse généralisée de la loi de Poisson de paramètre $\\lambda$, $(U_k)_{k \\ge 1}$ suite _i.i.d._ uniforme sur $[0,1]$ indépendante de $(X_{i,k})_{i\\ge1, k\\ge 1}$ et $(\\bar X_{i,k})_{i\\ge1, k\\ge 1}$ deux suites (doublement indicées) _i.i.d._ de loi log-normale (de paramètres $\\mu$ et $\\sigma$ inchangés).\n", "\n", "Un **troisième estimateur** $J^{(3)}_{n,h}(\\lambda)$ utilise des réalisations fortements corrélées de la loi de Poisson et des variables aléatoires log-normales communes. \n", "\n", "$$\n", " J^{(3)}_{n, h}(\\lambda) = \\frac{1}{2 h n} \\sum_{k = 1}^n \\bigl(\\mathbf{1}_{\\{S^{(\\lambda+h)}_k > K\\}} - \\mathbf{1}_{\\{S^{(\\lambda-h)}_k > K\\}} \\bigr),\n", "$$\n", "où pour $k \\ge 1$, $S^{(\\lambda+h)}_k = \\sum_{i = 1}^{G(\\lambda+h, U_k)} X_{i,k}$ et $S^{(\\lambda-h)}_k = \\sum_{i = 1}^{G(\\lambda-h, U_k)} X_{i,k}$ avec $G(\\lambda, u)$ l'inverse généralisée de la loi de Poisson de paramètre $\\lambda$, $(U_k)_{k \\ge 1}$ suite _i.i.d._ uniforme sur $[0,1]$ indépendante de $(X_{i,k})_{i\\ge1, k\\ge 1}$ une suite (doublement indicée) _i.i.d._ de loi log-normale." ] }, { "cell_type": "markdown", "id": "6299f26f-2fdf-4f27-84bb-6b7186050017", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: plusieurs estimateurs des différences finies \n", "\n", "On fixe les paramètres $\\lambda = 10$, $\\mu = 0.1$, $\\sigma = 0.3$ et $K = 20$.\n", "Programmer ces 3 estimateurs pour différentes valeurs de $h$ (par exemple, $h=1$, 0.5, 0.1 et 0.01), et donner le résultat des estimateurs Monte Carlo avec $n = 50\\,000$.\n", "\n", "Que se passe-t-il lorsque $h$ tend vers 0? Comparez le comportement pour ces 3 estimateurs. Il est très important de bien interpréter ces tableaux de résultats et de conclure qu'il faut utiliser l'estimateur $J^{(3)}_{n, h}(\\lambda)$ et en aucun cas l'estimateur $J^{(1)}_{n,h}(\\lambda)$.\n", "\n", "_Remarque_: on considère ici uniquement l'étude de l'erreur statistique dûe à la méthode de Monte Carlo. On ne considère pas l'erreur de biais qui décroît lorsque $h$ tend vers 0 et qui est peut-être non négligeable pour $h = 1$. Les IC construits ici sont biaisés et on ne peut pas affirmer que la vraie valeur est dans l'IC à 95% (au moins pour les grandes valeurs de $h$). " ] }, { "cell_type": "code", "execution_count": 2, "id": "677c88cf-29b1-4339-a6cc-87b6fa46b3dd", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.115201Z", "iopub.status.busy": "2025-12-20T09:23:27.115127Z", "iopub.status.idle": "2025-12-20T09:23:27.116822Z", "shell.execute_reply": "2025-12-20T09:23:27.116530Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "def monte_carlo(sample, proba = 0.95):\n", " mean = np.mean(sample)\n", " var = np.var(sample, ddof=1)\n", " alpha = 1 - proba \n", " quantile = stats.norm.ppf(1 - alpha/2) # fonction quantile \n", " ci_size = quantile * np.sqrt(var / sample.size)\n", " return (mean, var, mean - ci_size, mean + ci_size)" ] }, { "cell_type": "code", "execution_count": 3, "id": "fee8116c-f36d-4cbc-a725-857eeef2d44d", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.117731Z", "iopub.status.busy": "2025-12-20T09:23:27.117672Z", "iopub.status.idle": "2025-12-20T09:23:27.119178Z", "shell.execute_reply": "2025-12-20T09:23:27.118929Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "def simu_S(size, mu, sigma, lambd): \n", " sample_N = rng.poisson(size=size, lam=lambd)\n", " sample_S = np.empty(size)\n", " for k, Nk in enumerate(sample_N):\n", " sample_S[k] = np.sum(rng.lognormal(size=Nk, mean=mu, sigma = sigma))\n", " return sample_S" ] }, { "cell_type": "code", "execution_count": 4, "id": "5edc11d4-8d65-483c-9db9-e1a128934e06", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.119943Z", "iopub.status.busy": "2025-12-20T09:23:27.119896Z", "iopub.status.idle": "2025-12-20T09:23:27.121171Z", "shell.execute_reply": "2025-12-20T09:23:27.120904Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "lambd, mu, sigma = 10, 0.1, 0.3\n", "K = 20" ] }, { "cell_type": "code", "execution_count": 5, "id": "ff9cb886-0eba-4ee9-b180-409760b528d8", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.121985Z", "iopub.status.busy": "2025-12-20T09:23:27.121943Z", "iopub.status.idle": "2025-12-20T09:23:27.123424Z", "shell.execute_reply": "2025-12-20T09:23:27.123196Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "def J1(n, h): \n", " sample_Sph = simu_S(n, mu, sigma, lambd+h)\n", " sample_Smh = simu_S(n, mu, sigma, lambd-h)\n", " xph = (sample_Sph > K).astype(int)\n", " xmh = (sample_Smh > K).astype(int)\n", " return monte_carlo((xph - xmh)/(2*h))" ] }, { "cell_type": "code", "execution_count": 6, "id": "d71e1fcd-2412-42b5-93e5-8299a0c24749", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.124336Z", "iopub.status.busy": "2025-12-20T09:23:27.124294Z", "iopub.status.idle": "2025-12-20T09:23:27.125847Z", "shell.execute_reply": "2025-12-20T09:23:27.125474Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "# on crée une fonction pour ne pas répéter ces lignes de code\n", "def result_estimator(J, n=50000, hs=[1, 0.5, 0.1, 0.01]):\n", " result = [ J(n, h) for h in hs ]\n", " df = pd.DataFrame(result, \n", " columns=['mean', 'var', 'lower', 'upper'], \n", " index=hs)\n", " return df" ] }, { "cell_type": "code", "execution_count": 7, "id": "acd60d27-5f1b-4905-abd1-a6c9ede8bbf2", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.126649Z", "iopub.status.busy": "2025-12-20T09:23:27.126604Z", "iopub.status.idle": "2025-12-20T09:23:27.636394Z", "shell.execute_reply": "2025-12-20T09:23:27.636078Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meanvarlowerupper
1.000.017650.0127340.0166610.018639
0.500.017640.0444900.0157910.019489
0.100.010601.0279080.0017130.019487
0.010.04300108.950330-0.0484910.134491
\n", "
" ], "text/plain": [ " mean var lower upper\n", "1.00 0.01765 0.012734 0.016661 0.018639\n", "0.50 0.01764 0.044490 0.015791 0.019489\n", "0.10 0.01060 1.027908 0.001713 0.019487\n", "0.01 0.04300 108.950330 -0.048491 0.134491" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result_estimator(J1)" ] }, { "cell_type": "code", "execution_count": 8, "id": "3d60ff8d-b9c4-435e-a9d4-650b1fb57880", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.637327Z", "iopub.status.busy": "2025-12-20T09:23:27.637277Z", "iopub.status.idle": "2025-12-20T09:23:27.639101Z", "shell.execute_reply": "2025-12-20T09:23:27.638813Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "def J2(n, h): \n", " sample_U = rng.random(size=n)\n", " sample_Nph = stats.poisson(mu = lambd+h).ppf(sample_U)\n", " sample_Nmh = stats.poisson(mu = lambd-h).ppf(sample_U)\n", " sample_Sph = np.empty(n)\n", " sample_Smh = np.empty(n)\n", " for k, (Nph, Nmh) in enumerate(zip(sample_Nph, sample_Nmh)):\n", " sample_Sph[k] = np.sum(rng.lognormal(size=int(Nph), mean=mu, sigma = sigma))\n", " sample_Smh[k] = np.sum(rng.lognormal(size=int(Nmh), mean=mu, sigma = sigma))\n", " xph = (sample_Sph > K).astype(int)\n", " xmh = (sample_Smh > K).astype(int)\n", " return monte_carlo((xph - xmh)/(2*h))" ] }, { "cell_type": "code", "execution_count": 9, "id": "88489b47-47b7-4fcd-8847-7db7d5a666f1", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:27.639902Z", "iopub.status.busy": "2025-12-20T09:23:27.639858Z", "iopub.status.idle": "2025-12-20T09:23:28.333419Z", "shell.execute_reply": "2025-12-20T09:23:28.333045Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meanvarlowerupper
1.000.017840.0089620.0170100.018670
0.500.017940.0234590.0165970.019283
0.100.012600.4598500.0066560.018544
0.010.0000042.800856-0.0573440.057344
\n", "
" ], "text/plain": [ " mean var lower upper\n", "1.00 0.01784 0.008962 0.017010 0.018670\n", "0.50 0.01794 0.023459 0.016597 0.019283\n", "0.10 0.01260 0.459850 0.006656 0.018544\n", "0.01 0.00000 42.800856 -0.057344 0.057344" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result_estimator(J2)" ] }, { "cell_type": "code", "execution_count": 10, "id": "024afcc3-7876-4030-845a-4246ad447661", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:28.334369Z", "iopub.status.busy": "2025-12-20T09:23:28.334317Z", "iopub.status.idle": "2025-12-20T09:23:28.336324Z", "shell.execute_reply": "2025-12-20T09:23:28.335987Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [], "source": [ "def J3(n, h): \n", " sample_U = rng.random(size=n)\n", " sample_Nph = stats.poisson(mu = lambd+h).ppf(sample_U).astype(int)\n", " sample_Nmh = stats.poisson(mu = lambd-h).ppf(sample_U).astype(int)\n", " sample_Sph = np.empty(n)\n", " sample_Smh = np.empty(n)\n", " for k, (Nph, Nmh) in enumerate(zip(sample_Nph, sample_Nmh)):\n", " max_N = max(Nph, Nmh)\n", " sample_X = rng.lognormal(size=max_N, mean=mu, sigma = sigma)\n", " sample_Sph[k] = np.sum(sample_X[:Nph])\n", " sample_Smh[k] = np.sum(sample_X[:Nmh])\n", " xph = (sample_Sph > K).astype(int)\n", " xmh = (sample_Smh > K).astype(int)\n", " return monte_carlo((xph - xmh)/(2*h))" ] }, { "cell_type": "code", "execution_count": 11, "id": "b7eb6916-542b-4c57-bba5-e2d5ac8d9981", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-12-20T09:23:28.337129Z", "iopub.status.busy": "2025-12-20T09:23:28.337080Z", "iopub.status.idle": "2025-12-20T09:23:28.978392Z", "shell.execute_reply": "2025-12-20T09:23:28.978029Z" }, "slideshow": { "slide_type": "" }, "tags": [ "correction" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meanvarlowerupper
1.000.016970.0081970.0161760.017764
0.500.017140.0168470.0160020.018278
0.100.018300.0911670.0156530.020947
0.010.013000.6498440.0059340.020066
\n", "
" ], "text/plain": [ " mean var lower upper\n", "1.00 0.01697 0.008197 0.016176 0.017764\n", "0.50 0.01714 0.016847 0.016002 0.018278\n", "0.10 0.01830 0.091167 0.015653 0.020947\n", "0.01 0.01300 0.649844 0.005934 0.020066" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result_estimator(J3)" ] } ], "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" }, "toc-autonumbering": true }, "nbformat": 4, "nbformat_minor": 5 }