{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ddaf2c7",
   "metadata": {},
   "source": [
    "# TP Noté : Discrétisation du modèle de Black-Scholes et pricing d'un option à barrière"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7b1b8868",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09c6f71a",
   "metadata": {},
   "source": [
    "## Partie 1 : Discrétisation de l'EDS de Black-Scholes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87e70063",
   "metadata": {},
   "source": [
    "*Question 1*: La solution exacte de cette EDS est : $S_t = S_0 e^{(\\mu - \\frac{\\sigma^2}{2})t + \\sigma W_t}$\n",
    "\n",
    "*Question 2* :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d147570",
   "metadata": {},
   "outputs": [],
   "source": [
    "def S(t, S0, mu, sigma, W):\n",
    "    \"\"\"\n",
    "    Solution exacte de l'EDS de Black-Scholes\n",
    "    \"\"\"\n",
    "    return S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * W)\n",
    "\n",
    "\n",
    "def euler_maruyama(mu, sigma, T, N, X0=0.0):\n",
    "    \"\"\"\n",
    "    Simulation d'une EDS de Black-Scholes par la méthode d'Euler-Maruyama\n",
    "\n",
    "    Paramètres :\n",
    "        mu (float) : drift\n",
    "        sigma (float) : volatilité\n",
    "        T (int) : temps final\n",
    "        N (int) : nombre de pas de temps\n",
    "        X0 (float) : valeur initiale\n",
    "\n",
    "    Retourne :\n",
    "        t (array-like) : tableau des temps\n",
    "        X (array-like) : tableau des valeurs de l'EDS\n",
    "    \"\"\"\n",
    "    dt = T / N\n",
    "\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    X = np.zeros(N + 1)\n",
    "\n",
    "    X[0] = X0\n",
    "\n",
    "    dW = np.random.normal(0, np.sqrt(dt), N)\n",
    "\n",
    "    for i in range(N):\n",
    "        St = S(t[i], X[i], mu, sigma, dW[i])\n",
    "        X[i + 1] = X[i] + mu * St * dt + sigma * St * dW[i]\n",
    "\n",
    "    return t, X\n",
    "\n",
    "\n",
    "def plot_brownien(t, X, B=None):\n",
    "    \"\"\"\n",
    "    Plot la simulation d'Euler-Maruyama\n",
    "\n",
    "    Paramètres :\n",
    "        t (array-like) : tableau des temps\n",
    "        X (array-like) : tableau des valeurs de l'EDS\n",
    "        B (float) : barrière (optionnelle)\n",
    "    \"\"\"\n",
    "    plt.plot(t, X, alpha=0.5, label=\"Euler-Maruyama\")\n",
    "    plt.title(\"Simulation d'Euler-Maruyama pour une EDS\")\n",
    "\n",
    "    if B is not None:\n",
    "        plt.axhline(B, label=\"Barrière\", color=\"red\", linestyle=\"--\")\n",
    "\n",
    "    plt.legend()\n",
    "    plt.xlabel(\"Temps\")\n",
    "    plt.ylabel(\"X(t)\")\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a05f9468",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 59332 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "S0 = 100\n",
    "K = 100\n",
    "sigma = 0.2\n",
    "r = 0.05\n",
    "mu = r\n",
    "T = 1\n",
    "B = 120\n",
    "\n",
    "t, X = euler_maruyama(mu=mu, sigma=sigma, T=T, N=1000, X0=S0)\n",
    "plt.figure(figsize=(10, 6))\n",
    "plot_brownien(t, X)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a605d850",
   "metadata": {},
   "source": [
    "*Question 3*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1cf7a286",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54392 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h_list = [0.01, 0.1, 0.5, 1.0]\n",
    "N_list = [int(T / h) for h in h_list]\n",
    "\n",
    "np.random.seed(333)\n",
    "\n",
    "\n",
    "def plot_convergence(S0, mu, sigma, T):\n",
    "    \"\"\"\n",
    "    Plot la convergence du schéma d'Euler-Maruyama\n",
    "\n",
    "    Paramètres :\n",
    "        S0 (int) : valeur initiale\n",
    "        mu (float) : drift\n",
    "        sigma (float) : volatilité\n",
    "        T (int) : temps final\n",
    "    \"\"\"\n",
    "    errors = []\n",
    "\n",
    "    for N in N_list:\n",
    "        dt = T / N\n",
    "        dW = np.random.normal(0, np.sqrt(dt), N)\n",
    "\n",
    "        exact = S(T, S0, mu, sigma, dW)\n",
    "        _, X = euler_maruyama(mu=mu, sigma=sigma, T=T, N=N, X0=S0)\n",
    "\n",
    "        errors.append(np.max(np.abs(X[1:] - exact)))\n",
    "\n",
    "    plt.plot(np.log(h_list), np.log(errors), \"o-\", label=\"Erreur numérique\")\n",
    "    plt.plot(np.log(h_list), 0.5 * np.log(h_list), \"--\", label=\"Ordre 1/2\")\n",
    "    plt.plot(np.log(h_list), np.log(h_list), \"--\", label=\"Ordre 1\")\n",
    "    plt.plot(np.log(h_list), 2 * np.log(h_list), \"--\", label=\"Ordre 2\")\n",
    "    plt.xlabel(\"log(h)\")\n",
    "    plt.ylabel(\"log(Erreur)\")\n",
    "    plt.title(\"Convergence du schéma d'Euler-Maruyama\")\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plot_convergence(S0, r, sigma, T)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c71cb5db",
   "metadata": {},
   "source": [
    "L'ordre de convergence forte du schéma d'Euler-Maruyama est **${\\frac{1}{2}}$**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b384a5de",
   "metadata": {},
   "source": [
    "## Partie 2 : Pricing d'une option à barrière de type \"Up-and-Out Call\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b3d7e75",
   "metadata": {},
   "source": [
    "*Question 1*: Sous la mesure risque-neutre, l'équation de Black-Scholes s'écrit :\n",
    "\n",
    "$dS_t = r S_t dt + \\sigma S_t d\\tilde{W_t}$\n",
    "- $r$ est le taux d'intérêt\n",
    "- $\\tilde{W_t} = W_t + \\lambda t$ avec $\\lambda = \\frac{\\mu - r}{\\sigma}$\n",
    "\n",
    "La solution s'écrit : $S_t = S_0 e^{(r - \\frac{\\sigma^2}{2})t + \\sigma \\tilde{W_t}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbcbf31b",
   "metadata": {},
   "source": [
    "*Question 2*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ca7a09fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 63688 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "La barrière n'a pas été franchie\n"
     ]
    }
   ],
   "source": [
    "t, X = euler_maruyama(mu=r, sigma=sigma, T=T, N=1000, X0=S0)\n",
    "plt.figure(figsize=(10, 6))\n",
    "plot_brownien(t, X, B=B)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "def is_barrier_breached(X, B):\n",
    "    \"\"\"Renvoie True si la barrière est franchie, False sinon\n",
    "    La barrière est franchie si X >= B\n",
    "\n",
    "    Paramètres:\n",
    "        X (array-like): Trajectoire des valeurs\n",
    "        B (float): Valeur de la barrière\n",
    "\n",
    "    Retourne:\n",
    "        bool: True si la barrière est franchie, False sinon\n",
    "    \"\"\"\n",
    "    return any(X >= B)\n",
    "\n",
    "\n",
    "print(\n",
    "    \"La barrière a été franchie\"\n",
    "    if is_barrier_breached(X, B)\n",
    "    else \"La barrière n'a pas été franchie\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0440a6ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_browniens(trajectories, B):\n",
    "    \"\"\"Trace les trajectoires de Brownien et la barrière\n",
    "\n",
    "    Paramètres:\n",
    "        trajectories (list of tuples): Liste des trajectoires avec le temps et les valeurs\n",
    "        B (float): Valeur de la barrière\n",
    "    \"\"\"\n",
    "    for t, X in trajectories:\n",
    "        col = \"pink\" if is_barrier_breached(X, B) else \"lime\"\n",
    "        plt.plot(t, X, alpha=0.5, color=col)\n",
    "    plt.title(\"Simulation d'Euler-Maruyama pour une EDS\")\n",
    "\n",
    "    plt.axhline(B, label=\"Barrière\", color=\"red\", linestyle=\"--\")\n",
    "\n",
    "    plt.legend()\n",
    "    plt.xlabel(\"Temps\")\n",
    "    plt.ylabel(\"X(t)\")\n",
    "    plt.grid()\n",
    "\n",
    "\n",
    "def payoff(X, B, K):\n",
    "    \"\"\"Calcule le payoff d'une option en fonction des trajectoires.\n",
    "\n",
    "    Paramètres:\n",
    "        X (list): Liste des valeurs de la trajectoire.\n",
    "        B (float): Valeur de la barrière.\n",
    "        K (float): Prix d'exercice de l'option.\n",
    "\n",
    "    Retourne:\n",
    "        float: Valeur du payoff de l'option.\n",
    "    \"\"\"\n",
    "    if not is_barrier_breached(X, B):\n",
    "        return max(X[-1] - K, 0)\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "\n",
    "def call_BS(x):\n",
    "    \"\"\"Calcul du prix d'une option d'achat européenne selon le modèle de Black-Scholes en fonction de x.\n",
    "\n",
    "    Paramètres:\n",
    "        x (float): Le prix de l'actif sous-jacent.\n",
    "\n",
    "    Retourne:\n",
    "        float: Le prix de l'option d'achat européenne.\n",
    "    \"\"\"\n",
    "    d1 = (np.log(x / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))\n",
    "    d2 = d1 - sigma * np.sqrt(T)\n",
    "    return x * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)\n",
    "\n",
    "\n",
    "def compute_payoff_BS():\n",
    "    \"\"\"Calcul du prix d'une option d'achat Up-and-Out selon le modèle de Black-Scholes en fonction de la barrière.\n",
    "\n",
    "    Retourne:\n",
    "        float: Le prix de l'option d'achat Up-and-Out.\n",
    "    \"\"\"\n",
    "    lam = (r + 0.5 * sigma**2) / sigma**2\n",
    "    return (\n",
    "        call_BS(S0)\n",
    "        - call_BS(S0) * (S0 / B) ** (2 * lam)\n",
    "        + (S0 / B) ** (lam - 1)\n",
    "        * (call_BS(B**2 / S0) - (S0 / B) ** 2 * call_BS(B**2 / S0))\n",
    "    )\n",
    "\n",
    "\n",
    "def compute_payoff(trajectories, B, K):\n",
    "    \"\"\"Calcule le payoff d'une option en fonction des trajectoires.\n",
    "\n",
    "    Paramètres:\n",
    "        trajectories (list of tuples): Liste des trajectoires avec le temps et les valeurs.\n",
    "        B (float): Valeur de la barrière.\n",
    "        K (float): Prix d'exercice de l'option.\n",
    "\n",
    "    Retourne:\n",
    "        float: Valeur du payoff de l'option.\n",
    "    \"\"\"\n",
    "    return np.mean([payoff(X, B, K) for (_, X) in trajectories])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2de6a625",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 340452 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Moyenne des payoffs actualisés: 1.4768275985983736\n",
      "Comparaison avec la formule de Black-Scholes: 17.998941120931214\n"
     ]
    }
   ],
   "source": [
    "N_trajectories = 1000\n",
    "trajectories = [\n",
    "    (t, X)\n",
    "    for (t, X) in [\n",
    "        euler_maruyama(mu=r, sigma=sigma, T=T, N=1000, X0=S0)\n",
    "        for _ in range(N_trajectories)\n",
    "    ]\n",
    "]\n",
    "plt.figure(figsize=(10, 6))\n",
    "plot_browniens(trajectories, B=B)\n",
    "plt.show()\n",
    "\n",
    "payoffs = compute_payoff(trajectories, B, K)\n",
    "payoffs_BS = compute_payoff_BS()\n",
    "print(\"Moyenne des payoffs actualisés:\", np.mean(payoffs))\n",
    "print(\"Comparaison avec la formule de Black-Scholes:\", payoffs_BS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5671e9e8",
   "metadata": {},
   "source": [
    "## Partie 3 : Analyse d'erreur et convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2dcb9c95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54452 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_list = [10, 25, 50, 100, 200, 500, 1000]\n",
    "\n",
    "np.random.seed(333)\n",
    "\n",
    "\n",
    "def plot_payoff_errors():\n",
    "    \"\"\"Trace l'erreur de convergence du payoff actualisé en fonction de N.\"\"\"\n",
    "    errors = []\n",
    "\n",
    "    for N in N_list:\n",
    "        trajectories = [\n",
    "            (t, X)\n",
    "            for (t, X) in [\n",
    "                euler_maruyama(mu=r, sigma=sigma, T=T, N=N, X0=S0)\n",
    "                for _ in range(N_trajectories)\n",
    "            ]\n",
    "        ]\n",
    "        payoff_BS = compute_payoff_BS()\n",
    "        payoffs = compute_payoff(trajectories, B, K)\n",
    "\n",
    "        errors.append(np.max(np.abs(payoffs - payoff_BS)))\n",
    "\n",
    "    plt.plot(np.log(N_list), np.log(errors), \"o-\", label=\"Erreur numérique\")\n",
    "    plt.plot(np.log(N_list), 0.5 * np.log(N_list), \"--\", label=\"Ordre 1/2\")\n",
    "    plt.plot(np.log(N_list), np.log(N_list), \"--\", label=\"Ordre 1\")\n",
    "    plt.plot(np.log(N_list), 2 * np.log(N_list), \"--\", label=\"Ordre 2\")\n",
    "    plt.xlabel(\"log(h)\")\n",
    "    plt.ylabel(\"log(Erreur)\")\n",
    "    plt.title(\"Convergence de l'erreur du payoff actualisé\")\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plot_payoff_errors()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8005cd64",
   "metadata": {},
   "source": [
    "*Question 2*: On observe que plus le pas de discrétisation est petit, plus l’approximation du payoff se rapproche de sa courbe réelle. Cela permet de mieux capturer les franchissements de barrière pouvant se produire entre deux instants de discrétisation, réduisant ainsi les erreurs liées à la non-détection de ces événements."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "studies",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
