{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aa734dd-bd14-4f15-aa30-818a7792379c",
   "metadata": {},
   "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee556227-6551-408d-af0a-d9d531990e41",
   "metadata": {},
   "source": [
    "# Modèle de Black-Scholes 1d, Put Bermudéen\n",
    "\n",
    "Soit $(S_t)_{t \\in [0,T]}$ un processus de Black-Scholes (Brownien géométrique) de paramètre $r$, $\\sigma$ et de valeur initiale fixée $x_0 > 0$, c'est à dire $S_t = x_0 e^{(r-\\frac{\\sigma^2}{2}) t + \\sigma W_t}$ où $(W_t)_{t \\in [0,T]}$ est un mouvement Brownien standard.\n",
    "\n",
    "On considère des dates discrètes fixées $t_n = n \\frac{T}{N}$ pour $n = 0, \\dots, N$. La valeur de l'actif aux instants $t_n$ forme une chaine de Markov que l'on note $(X_n)_{n = 0,\\dots,N}$ c'est à dire \n",
    "$$\n",
    "    \\forall n=0, \\dots, N, \\quad X_n = S_{t_n}\n",
    "$$\n",
    "Le put Bermudéen est une option que l'on peut exercer à toute date $t_n$.  Si on exerce en $t_n$ le gain associé (le payoff) est \n",
    "$$\n",
    "    Z_n = \\varphi(n, X_n) = e^{-r n \\frac{T}{N}} (K - X_n)_+.\n",
    "$$\n",
    "\n",
    "On s'intéresse donc au problème d'arrêt optimal à temps discret \n",
    "$$\n",
    "    V_0(x_0) = \\sup_{\\tau \\in \\mathcal{T}_0} \\mathbf{E}\\big[ Z_\\tau \\big]\n",
    "    = \\sup_{\\tau \\in \\mathcal{T}_0} \\mathbf{E}\\big[ \\varphi(\\tau, X_\\tau) \\big],,\n",
    "$$\n",
    "où $\\mathcal{T}_0$ est l'ensemble des temps d'arrêts à valeur dans $\\{0,\\dots,N\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1051121-e67c-408d-bc6d-0ce3b5926536",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 0.1\n",
    "sigma = 0.25\n",
    "x0 = 100\n",
    "K = 110\n",
    "N, T = 10, 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dc270ab-dbfc-4e96-84d6-704101871385",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simu_BS1d(size_path, size_sample): \n",
    "    h = T/size_path\n",
    "    brown_acc = np.sqrt(h)*rng.standard_normal(size=(size_path, size_sample))\n",
    "    sample = np.zeros(shape=(size_path+1, size_sample))\n",
    "    sample[0] = x0\n",
    "    for n in range(1, size_path+1):\n",
    "        sample[n] = sample[n-1] * np.exp((r - 0.5 * sigma**2)*h + sigma*brown_acc[n-1])\n",
    "    return sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2696bde7-dc9c-4db4-a850-acc2dbad24dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def payoff_phi(n, x): \n",
    "    return np.exp(-r*n*T/N) * np.maximum(K-x, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c072d3a-0885-488b-939d-00ddcf3c25e1",
   "metadata": {},
   "source": [
    "# Régression sur une base de fonctions \n",
    "\n",
    "On considère $m$ fonctions $(e_k)_{1 \\le k \\le m}$ et la projection sur le sous-espace engendré par les fonctions $(e_k)_{1 \\le k \\le m}$ c'est à dire la fonction paramètrique \n",
    "$$\n",
    "    \\Phi(x; \\theta) = \\sum_{k=1}^m \\theta_k e_k(x) \\quad \\text{avec $\\theta \\in \\mathbf{R}^m$}.\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9c4356f-ec66-49c2-ad6c-8f3e3801ee73",
   "metadata": {},
   "outputs": [],
   "source": [
    "def base1_ek(x): \n",
    "    return np.array([np.ones_like(x),  x, x**2, x**3])\n",
    "\n",
    "# même base en ajoutant le payoff $(K-x)_+$\n",
    "def base2_ek(x): \n",
    "    return np.array([np.ones_like(x), np.maximum(K-x,0), x, x**2, x**3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7116e935-799d-4df4-b1d3-01251aaca3d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "base1_ek(0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ed61491-a590-4ccc-bb49-cee390e2543a",
   "metadata": {},
   "outputs": [],
   "source": [
    "base1_ek(np.array([0.2, 0.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "952e6ad6-9a02-4803-a3a3-12a0092ac244",
   "metadata": {},
   "outputs": [],
   "source": [
    "def theta_by_regression(payoff, x, ek=base1_ek):\n",
    "    norm = (ek(x) @ ek(x).T)\n",
    "    return np.linalg.inv(norm) @ (ek(x) @ payoff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a594ed39-473d-4202-8838-5882b5406f67",
   "metadata": {},
   "outputs": [],
   "source": [
    "def function_Phi(x, theta, ek=base1_ek): \n",
    "    return np.dot(theta, ek(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e8e0960-9c3b-4ede-952e-2b6f01a88934",
   "metadata": {},
   "source": [
    "# Algorithme de Longstaff-Schwartz\n",
    "\n",
    "On fabrique l'ensemble des scénarios: l'échantillon de $M$ trajectoires $(X^{(j)}_n)_{n=0,\\dots,N}$, $1 \\le j \\le M$ et les payoffs associés $(Z^{(j)}_n)_{n=0,\\dots,N}$, $1 \\le j \\le M$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92248d5c-ef23-4477-a211-5ecb487f6c25",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = int(1e6)\n",
    "sample_X = simu_BS1d(N, M)\n",
    "payoffs_Z = np.empty_like(sample_X)\n",
    "for n in range(0, N+1):\n",
    "    payoffs_Z[n] = payoff_phi(n, sample_X[n])\n",
    "    \n",
    "print(\"Shape of sample_X: \", sample_X.shape)\n",
    "print(\"Shape of payoffs_Z:\", payoffs_Z.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20bce2e2-6a4e-41c6-adfc-f2d07bae7fa6",
   "metadata": {},
   "source": [
    "\\begin{equation} \\tag{$A_{LS}$}\n",
    "\\begin{cases}\n",
    "    \\tau^{(j)}_N = N, & \\text{condition terminale} \\\\\n",
    "    \\tau^{(j)}_{n} = n \\mathbf{1}_{Z^{(j)}_n \\ge \\Phi(X^{(j)}_n; \\theta_n)} \n",
    "    + \\tau^{(j)}_{n+1} \\mathbf{1}_{Z^{(j)}_n < \\Phi(X^{(j)}_n; \\theta_n)}, & 1 \\le n \\le N-1\n",
    "\\end{cases}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f339099a-ef12-4547-8190-89969a3056cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = base1_ek(sample_X[0]).shape[0]   # on récupère le nombre de fonctions utilisées\n",
    "thetas = np.zeros((N, m))            # on va sauver les paramètres \"optimaux\"\n",
    "\n",
    "# tau_opt = N * np.ones(M, dtype=int)\n",
    "payoff_opt = payoffs_Z[N].copy()\n",
    "\n",
    "for n in reversed(range(1, N)):\n",
    "    thetas[n] = theta_by_regression(payoff_opt, sample_X[n])\n",
    "    stop_at_n = payoffs_Z[n] >= function_Phi(sample_X[n], thetas[n]) \n",
    "    # tau_opt[stop_at_n] = n \n",
    "    payoff_opt[stop_at_n] = payoffs_Z[n, stop_at_n].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fdb3647-99c5-4497-9c3a-c041d21d9a50",
   "metadata": {},
   "outputs": [],
   "source": [
    "thetas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e89a388-8a2e-41a1-90e6-7573690d0005",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(90, 130, 1000)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(xx, payoff_phi(N, xx))\n",
    "for n in range(1, N):\n",
    "    ax.plot(xx, np.maximum(payoff_phi(n, xx), function_Phi(xx, thetas[n], base1_ek)), label=fr\"V{n}\")\n",
    "    ax.legend()\n",
    "    ax.set_title(\"Les fonctions valeurs approchées (base1)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "feb89557-9fd9-4839-8f84-5a043ca8c40d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(90, 130, 1000)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(xx, payoff_phi(N, xx))\n",
    "for n in range(1, N):\n",
    "    ax.plot(xx, function_Phi(xx, thetas[n]), label=fr\"C{n}\")\n",
    "    ax.legend()\n",
    "    ax.set_title(\"Les fonctions de continuation approchées (base1)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1295068-fcdd-4d31-92b2-e798c47bc99a",
   "metadata": {},
   "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",
    "    result = { 'mean': mean, 'var': var, \n",
    "               'lower': mean - ci_size, \n",
    "               'upper': mean + ci_size }\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6758ea-4446-4d65-ac7e-27a80fd4a34c",
   "metadata": {},
   "outputs": [],
   "source": [
    "monte_carlo(payoff_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f248c7a5-a16a-4434-8197-25e60c428e59",
   "metadata": {},
   "source": [
    "# Comparaison des 2 bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11f69e46-ae2d-436e-959e-a8236d663d13",
   "metadata": {},
   "outputs": [],
   "source": [
    "bases = [ \n",
    "    { \"name\": \"base1\", \"ek\": base1_ek }, \n",
    "    { \"name\": \"base2\", \"ek\": base2_ek },\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c695b71a-91ef-4dac-9048-2cf2980a256f",
   "metadata": {},
   "outputs": [],
   "source": [
    "for base in bases:\n",
    "    m = base[\"ek\"](sample_X[0]).shape[0]    # on récupère le nombre de fonctions utilisées\n",
    "    thetas = np.zeros((N, m))            # on va sauver les paramètres \"optimaux\"\n",
    "\n",
    "    tau_opt = N * np.ones(M, dtype=int)\n",
    "    payoff_opt = payoffs_Z[N].copy()\n",
    "\n",
    "    for n in reversed(range(1, N)):\n",
    "        thetas[n] = theta_by_regression(payoff_opt, sample_X[n], base[\"ek\"])\n",
    "        stop_at_n = payoffs_Z[n] >= function_Phi(sample_X[n], thetas[n], base[\"ek\"]) \n",
    "        tau_opt[stop_at_n] = n \n",
    "        payoff_opt[stop_at_n] = payoffs_Z[n, stop_at_n].copy()\n",
    "    \n",
    "    base[\"thetas\"] = thetas\n",
    "    base[\"payoff_opt\"] = payoff_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85713019-d6c1-4e36-9969-39ae7d34913e",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.linspace(90, 130, 1000)\n",
    "\n",
    "fig, axs = plt.subplots(ncols=2, figsize=(12,6))\n",
    "for ax, base in zip(axs, bases): \n",
    "    ax.plot(xx, payoff_phi(N, xx))\n",
    "    for n in range(1, N):\n",
    "        ax.plot(xx, np.maximum(payoff_phi(n, xx), \n",
    "                               function_Phi(xx, base[\"thetas\"][n], base[\"ek\"])), \n",
    "                label=fr\"V{n}\")\n",
    "        ax.legend()\n",
    "        name = base[\"name\"]\n",
    "        ax.set_title(f\"Les fonctions valeurs approchées ({name})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526ac04a-3322-4807-b6ee-7cd34693f296",
   "metadata": {},
   "outputs": [],
   "source": [
    "for base in bases: \n",
    "    print(\"Price with\", base[\"name\"]) \n",
    "    print(monte_carlo(base[\"payoff_opt\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6d669fa-ef01-4ba0-8249-2a4c86babcab",
   "metadata": {},
   "source": [
    "# Resimulation\n",
    "\n",
    "..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.14.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}