{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44f28e3d-80ab-4a08-a928-96bfc77b58b3",
   "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": "code",
   "execution_count": null,
   "id": "b036e38d-b7e4-4481-8c9c-ca50495e6e7e",
   "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": "893bfd44-77e9-43dd-9265-465a9dee44e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b1fb65d-744d-48a0-89ca-a8d201d8ed06",
   "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": "2cf7b92a-d84e-4a9b-8962-d8548e144c54",
   "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": "de5bb5f8-4049-4472-a429-e817a406caa7",
   "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": "99fe6991-c8e1-4cbb-beed-7944afdec885",
   "metadata": {},
   "outputs": [],
   "source": [
    "def payoff_phi(n, x): \n",
    "    return np.exp(-r*n*T/N) * np.maximum(K-x, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7466320f-07de-4253-8b9f-08e0f65c21ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 2**20\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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c69c71c7-3c56-4f6e-a84b-09d91f2788ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_X.shape "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e10dc1c3-dd43-40b2-a284-f03426531d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "payoffs_Z.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa8f6f27-9e6e-4789-8e90-67e455c45aeb",
   "metadata": {},
   "source": [
    "# Feedforward Neural Network \n",
    "\n",
    "On considère la fonction paramétrique, pour un $L \\ge 2$ fixé, obtenue par composition de fonctions \n",
    "$$\n",
    "    \\Phi(x; \\theta) = A_L \\circ \\big( \\sigma \\circ A_{L-1} \\big) \\circ \\cdots \\circ \\big( \\sigma \\circ A_1 \\big)(x)\n",
    "$$\n",
    "avec $\\theta = \\big( W_\\ell, \\beta_\\ell \\big)_{1 \\le \\ell \\le L}$ et pour tout $\\ell \\in \\{1,\\dots,L\\}$ la fonction affine $A_\\ell$ s'écrit $A_\\ell(x) = W_\\ell x + \\beta_l$. La fonction $\\sigma$ est la fonction ReLU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf01bfc4-6af5-40ab-9923-6f9bad63d51c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NeuralNetwork(nn.Module):\n",
    "    def __init__(self, input_size, layer_sizes, output_size):\n",
    "        super().__init__()\n",
    "        layers = [ nn.Linear(input_size, layer_sizes[0]),\n",
    "                   nn.ReLU() ]\n",
    "        for (ls_in, ls_out) in zip(layer_sizes, layer_sizes[1:]):\n",
    "            layers.append(nn.Linear(ls_in, ls_out))\n",
    "            layers.append(nn.ReLU())\n",
    "        layers.append(nn.Linear(layer_sizes[-1], output_size))\n",
    "        self.linear_relu_stack = nn.Sequential(*layers)\n",
    "\n",
    "    def forward(self, x):\n",
    "        out = self.linear_relu_stack(x)\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53395bec-2819-4aaa-91e3-d8c1337ca191",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = NeuralNetwork(input_size=1, layer_sizes=[16, 16, 32, 8], output_size=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34179829-983d-427f-a669-4d8d796b613c",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd36492e-b763-4563-b2df-9a65e444dee0",
   "metadata": {},
   "source": [
    "De façon plus simple on peut utiliser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "411b36a7-2757-4fa1-91f9-ddc409ea2a47",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NeuralNetwork2(nn.Module):\n",
    "    def __init__(self, input_size, hidden_size, output_size):\n",
    "        super().__init__()\n",
    "        self.linear_relu_stack = nn.Sequential(\n",
    "            nn.Linear(input_size, hidden_size),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(hidden_size, hidden_size),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(hidden_size, output_size)\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        out = self.linear_relu_stack(x)\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c0abb1-4603-491f-85e0-f99c5a266e13",
   "metadata": {},
   "outputs": [],
   "source": [
    "#phi = NeuralNetwork2(input_size=1, hidden_size=16, output_size=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a1985c9-3448-44b0-8698-47407e79117f",
   "metadata": {},
   "source": [
    "# Etape préliminaire: fonction de continuation en $N-1$\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$. \n",
    "\n",
    "Et on stocke nos données: $(xj, yj)_{1 \\le j \\le M}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b25bd96-655f-456d-9d9b-5617328bdf77",
   "metadata": {},
   "outputs": [],
   "source": [
    "#device = \"mps\" if torch.mps.is_available() else \"cpu\"\n",
    "device = \"cpu\"\n",
    "print(f\"Using {device} device\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f814077f-2f39-4f85-98d3-b29e9429506f",
   "metadata": {},
   "outputs": [],
   "source": [
    "xj = torch.from_numpy(sample_X[N-1].astype(np.float32))[:,None].to(device)\n",
    "yj = torch.from_numpy(payoffs_Z[N].astype(np.float32))[:,None].to(device)\n",
    "print(\"Shape of xj:\", xj.shape)\n",
    "print(\"Shape of yj:\", yj.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71330012-d9b1-4025-8d8b-038ac649c0c8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xj.min(), xj.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "652c9648-fccf-48a3-be39-db05f2801ff5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "yj.min(), yj.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7bc14e6-5610-4057-995a-b4533c4153ee",
   "metadata": {},
   "source": [
    "## Sans normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807c48f0-d84a-4197-864d-aace6b829e4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 20\n",
    "batch_size = 2*1024\n",
    "learning_rate = 1e-3\n",
    "\n",
    "data_size = xj.shape[0]\n",
    "n_upd = data_size // batch_size \n",
    "\n",
    "phi = NeuralNetwork(input_size=1, layer_sizes=[32, 32], output_size=1).to(device)\n",
    "optimizer = torch.optim.Adam(phi.parameters(), lr=learning_rate)\n",
    "with tqdm(range(n_epochs)) as nrange:\n",
    "    for n in nrange:\n",
    "        indexes = torch.randperm(data_size)   # on pourrait resimuler...\n",
    "        for k in range(n_upd):\n",
    "            idx = indexes[ k*batch_size : (k+1)*batch_size ] # à la main mais on peut utiliser dataloader\n",
    "            loss = ((yj[idx] - phi(xj[idx]))**2).mean()  # cf torch.MSEloss avec reduction = 'mean' \n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()  # theta_{k+1} = theta_k - ...\n",
    "        nrange.set_postfix(loss=loss.item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c216dedd-9c9f-41d8-b645-c1c9bb178961",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi.cpu()\n",
    "xx = torch.linspace(90, 130, 1000)[:,None].float()\n",
    "yy = phi(xx)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "with torch.no_grad():\n",
    "    ax.plot(xx, payoff_phi(N, xx))\n",
    "    ax.plot(xx, yy)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adb8c6c3-9a78-4e81-972f-121c456d5a05",
   "metadata": {},
   "source": [
    "## Avec normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b19278f3-f262-46d5-8c7b-ad1dc4ec91bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "cst_norm = xj.mean(), xj.std(), yj.mean(), yj.std()\n",
    "x_mean, x_std, y_mean, y_std = cst_norm\n",
    "\n",
    "def normalize_x(xi):\n",
    "    return (xi-x_mean)/x_std\n",
    "def normalize(xi, yi):\n",
    "    return (xi-x_mean)/x_std, (yi-y_mean)/y_std \n",
    "def unnormalize(xi, yi):\n",
    "    return xi*x_std + x_mean, yi*y_std + y_mean "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc9da3b-5f37-4ec0-a879-315d409bc343",
   "metadata": {},
   "outputs": [],
   "source": [
    "xj_, yj_ = normalize(xj, yj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3259d210-a397-4b63-b529-98b797c79d83",
   "metadata": {},
   "outputs": [],
   "source": [
    "xj_.min(), xj_.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "687cfd4f-ce65-4de3-a145-3937d29dd2ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 20\n",
    "batch_size = 2*1024\n",
    "learning_rate = 1e-3\n",
    "\n",
    "data_size = xj.shape[0]\n",
    "n_upd = data_size // batch_size \n",
    "\n",
    "phi = NeuralNetwork(input_size=1, layer_sizes=[32, 32], output_size=1).to(device)\n",
    "optimizer = torch.optim.Adam(phi.parameters(), lr=learning_rate)\n",
    "with tqdm(range(n_epochs)) as nrange:\n",
    "    for n in nrange:\n",
    "        indexes = torch.randperm(data_size)\n",
    "        for k in range(n_upd):\n",
    "            idx = indexes[k*batch_size:(k+1)*batch_size] \n",
    "            loss = ((yj_[idx] - phi(xj_[idx]))**2).mean()\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "        nrange.set_postfix(loss=loss.item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0e61ff-f5d9-4f7a-8fbc-5af2fb5b7c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx_, _ = normalize(torch.linspace(90, 130, 1000)[:,None], 0)\n",
    "yy_ = phi(xx_)\n",
    "xx, yy = unnormalize(xx_, yy_)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "with torch.no_grad():\n",
    "    ax.plot(xx, payoff_phi(N, xx))\n",
    "    ax.plot(xx, yy)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc24876b-f850-4887-8302-acbe73e92524",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf3aa7b8-7309-461a-bcdd-cfb9acd0f2ff",
   "metadata": {},
   "source": [
    "## Création d'une fonction `learning`, approche **données**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0edcf1b-f12e-4d20-b17b-d9eb5bace03e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def theta_by_learning(xj, yj, layer_sizes, n_epochs, batch_size, \n",
    "                      learning_rate=1e-3, \n",
    "                      theta_init=None, device=device):\n",
    "    data_size, input_size = xj.shape\n",
    "    n_upd = data_size // batch_size \n",
    "    \n",
    "    phi = NeuralNetwork(input_size, layer_sizes, output_size=1).to(device)\n",
    "    if theta_init is not None:\n",
    "        phi.load_state_dict(theta_init.copy())\n",
    "    optimizer = torch.optim.Adam(phi.parameters(), lr=learning_rate)\n",
    "    \n",
    "    with tqdm(range(n_epochs)) as nrange:\n",
    "        for n in nrange:\n",
    "            indexes = torch.randperm(data_size)\n",
    "            for k in range(n_upd):\n",
    "                idx = indexes[ k*batch_size : (k+1)*batch_size ] \n",
    "                loss = ((yj[idx] - phi(xj[idx]))**2).mean()\n",
    "                optimizer.zero_grad()\n",
    "                loss.backward()\n",
    "                optimizer.step()\n",
    "            nrange.set_postfix(loss=loss.item())\n",
    "    return phi.state_dict(), phi(xj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ba2ed60-135a-425d-8a12-e134bf637f4c",
   "metadata": {},
   "source": [
    "# Algorithme de Longstaff-Schwartz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6867aee9-c24d-4c10-8f39-dc961ee9d439",
   "metadata": {},
   "outputs": [],
   "source": [
    "xj_path = torch.from_numpy(sample_X.astype(np.float32)).to(device).float()\n",
    "yj_path = torch.from_numpy(payoffs_Z.astype(np.float32)).to(device).float()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9833bc0-00c6-4ac4-bde1-b62f642a043d",
   "metadata": {},
   "source": [
    "## Sans normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1af2fba-b863-4b20-91d3-019c672a605f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "layer_sizes = [16, 16]\n",
    "phi = NeuralNetwork(1, layer_sizes, 1).to(device)\n",
    "theta_n = phi.state_dict()\n",
    "thetas = []\n",
    "\n",
    "payoff_opt = yj_path[N].clone()\n",
    "for n in reversed(range(1, N)):\n",
    "    print(\"learning at timestep:\", n)\n",
    "    theta_n, cont_value_n =  theta_by_learning(xj_path[n][:,None], payoff_opt[:,None], \n",
    "                                               layer_sizes, \n",
    "                                               20 if n==N-1 else 2, \n",
    "                                               2*1024, theta_init = theta_n, \n",
    "                                               device=device)\n",
    "    thetas.insert(0, theta_n)\n",
    "    stop_at_n = yj_path[n,:] >= cont_value_n.flatten()\n",
    "    payoff_opt[stop_at_n] = yj_path[n, stop_at_n].clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1120ae00-13d8-48a9-8e23-00696fabcb88",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi.cpu()\n",
    "def continuation(n, xx):\n",
    "    phi.load_state_dict(thetas[n-1])\n",
    "    return phi(xx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c898bdaa-676c-463a-a386-09239beb0089",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = torch.linspace(90, 130, 1000)[:,None].float()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "with torch.no_grad():\n",
    "    ax.plot(xx, payoff_phi(N, xx))\n",
    "    for n in range(1, N):\n",
    "        ax.plot(xx, continuation(n, xx), label=fr\"C{n}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1af464d6-80a7-49e7-83e6-db2898a9a566",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = payoff_opt\n",
    "monte_carlo(sample.detach().cpu().numpy().flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "102ccf1e-11d5-4103-8fe2-f4dec57ae065",
   "metadata": {},
   "source": [
    "## Avec normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3378495-d2b1-43c3-8f51-d176a2d63d58",
   "metadata": {},
   "outputs": [],
   "source": [
    "cst_norm = xj_path.mean(), xj_path.std(), yj_path.mean(), yj_path.std()\n",
    "x_mean, x_std, y_mean, y_std = cst_norm\n",
    "\n",
    "def normalize_x(xi):\n",
    "    return (xi-x_mean)/x_std\n",
    "def normalize(xi, yi):\n",
    "    return (xi-x_mean)/x_std, (yi-y_mean)/y_std \n",
    "def unnormalize(xi, yi):\n",
    "    return xi*x_std + x_mean, yi*y_std + y_mean "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eecdf068-88b6-4b8b-a789-9bfa6ad9b9ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "xj_path, yj_path = normalize(xj_path, yj_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9bd09d5-764c-4117-86d0-bf4302b888c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "layer_sizes = [32, 16]\n",
    "phi = NeuralNetwork(1, layer_sizes, 1).to(device)\n",
    "theta_n = phi.state_dict()\n",
    "thetas = []\n",
    "\n",
    "payoff_opt = yj_path[N].clone()\n",
    "for n in reversed(range(1, N)):\n",
    "    print(\"learning at timestep:\", n)\n",
    "    theta_n, cont_value_n =  theta_by_learning(xj_path[n][:,None], payoff_opt[:,None], \n",
    "                                               layer_sizes, \n",
    "                                               40 if n == N-1 else 5, \n",
    "                                               2*1024, theta_init = theta_n)\n",
    "    thetas.insert(0, theta_n)\n",
    "    stop_at_n = yj_path[n,:] >= cont_value_n.flatten()\n",
    "    payoff_opt[stop_at_n] = yj_path[n, stop_at_n].clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b1a0a25-c916-4ea3-8ad2-1015010ce5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def continuation(n, xx):\n",
    "    xx_ = normalize_x(xx)\n",
    "    phi.load_state_dict(thetas[n-1])\n",
    "    xx, yy = unnormalize(xx_ , phi(xx_))\n",
    "    return yy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98248008-6afb-463c-91cf-e9a818ffdbd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = torch.linspace(90, 130, 1000)[:,None].float()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "with torch.no_grad():\n",
    "    ax.plot(xx, payoff_phi(N, xx))\n",
    "    for n in range(1, N):\n",
    "        ax.plot(xx, continuation(n, xx), label=fr\"C{n}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64320f40-e641-4d16-ba0b-dc48ba93805d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = y_mean + payoff_opt*y_std\n",
    "monte_carlo(sample.detach().numpy().flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b6f8e64-2da2-45b8-909f-7b86a7786bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#thetas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a584158e-4762-461d-a073-5e8ebb27a18a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Recalcul d'un prix sur de nouvelles trajectoires "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d601d2c-9ffb-465e-9f77-c83e76e32c6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simu_data_torch(size_sample, device=device): \n",
    "    h = T/N\n",
    "    brown_acc = np.sqrt(h)*torch.randn(size=(N, size_sample), device=device)\n",
    "    sample_X = torch.zeros(size=(N+1, size_sample))\n",
    "    sample_X[0] = x0\n",
    "    for n in range(1, N+1):\n",
    "        sample_X[n] = sample_X[n-1] * torch.exp((r - 0.5 * sigma**2)*h + sigma*brown_acc[n-1])\n",
    "\n",
    "    payoffs_Z = torch.empty_like(sample_X)\n",
    "    for n in range(0, N+1):\n",
    "        payoffs_Z[n] = np.exp(-r*n*h) * torch.relu(K-sample_X[n])\n",
    "    return normalize(sample_X, payoffs_Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d3e95b5-b9e3-4b00-9039-60ec96012bbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 2**22\n",
    "xj_path, yj_path = simu_data_torch(sample_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "996301b0-8a58-472f-a692-e5b71de50957",
   "metadata": {},
   "outputs": [],
   "source": [
    "payoff_opt = yj_path[N].clone()\n",
    "for n in reversed(range(1, N)):\n",
    "    phi.load_state_dict(thetas[n-1])\n",
    "    with torch.no_grad():\n",
    "        cont_at_n = phi(xj_path[n,:][:,None]).flatten()\n",
    "    stop_at_n = yj_path[n,:] >= cont_at_n\n",
    "    payoff_opt[stop_at_n] = yj_path[n, stop_at_n].clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d35d693c-2053-4853-b974-689f8dd8ae47",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample = y_mean + payoff_opt*y_std\n",
    "monte_carlo(sample.detach().numpy().flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef4e7d9-13e6-4d53-8143-a7490ab973c7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}