{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training \"inspecteur modèles\" BNP Paribas - Monte-Carlo methods - January 7, 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook 4 - Examples of Nested Monte-Carlo and empirical regression methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1. Nested Monte-Carlo for expected loss\n",
    "\n",
    "\n",
    "#### $\\blacktriangleright$ Setting\n",
    "We consider a Black-Scholes model at discrete dates $t_i$:\n",
    "\n",
    "$$\n",
    "S_{t_i} = S_0 \\, e^{\\sigma \\, W_{t_i} - \\frac 12 \\sigma^2 t_i},\n",
    "\\qquad t_0 < t_1 < \\dots < t_n,\n",
    "$$\n",
    "\n",
    "which we use to generate option prices (note we consider zero interest rate $r=0$). \n",
    "\n",
    "In particular, the price at time $t_1$ of a Put option with maturity $t_2$ and strike price $K$ is given by \n",
    "\n",
    "$$\n",
    "P(t_1, S_{t_1}) = \\mathbb E \\bigl[ (K - S_{t_2})^+ \\big| S_{t_1} \\bigr]\n",
    "$$\n",
    "\n",
    "where $x^+$ denotes the positive part of $x$; the explicit Black-Scholes formula is of course available for the function $P(t,S)$\n",
    "\n",
    "$$\n",
    "P(t,S) = K \\mathcal N \\biggl(\\frac{\\ln \\frac K {S}}{\\sigma \\sqrt{t_2 - t}} + \\frac 12 \\sigma \\sqrt{t_2 - t} \\biggr)\n",
    "- S \\, \\mathcal N \\biggl(\\frac{\\ln \\frac K {S}}{\\sigma \\sqrt{t_2 - t}} - \\frac 12 \\sigma \\sqrt{t_2 - t} \\biggl),\n",
    "\\qquad t \\le t_2, \\ S > 0\n",
    "$$\n",
    "\n",
    "where $\\mathcal N$ is the standard gaussian cdf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We wish to apply Monte-Carlo simulation to the evaluation of the _nested_ expectation\n",
    "\n",
    "$$\n",
    "I = \\mathbb E \\Bigl[\\Bigl(P(t_1, S_{t_1}) - P(t_0, S_0) \\Bigr)^+ \\Bigr]\n",
    "$$\n",
    "\n",
    "which represents the _expected loss_ (or exposure) at time $t_1$ of a short position on the put option of maturity $t_2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as sps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  1. Benchmark value (semi-explicit in this case)\n",
    "\n",
    "In the framework of the model above, it is possible to evaluate the nested expectation $I$ using numerical evaluation of a 1D integral with a deterministic quadrature formula, since\n",
    "\n",
    "$$\n",
    "I = \\int_{\\mathbb R} \\Bigl(P \\bigl(t_1, S_0 \\, e^{\\sigma \\, \\sqrt{t_1} \\, y \\, - \\, \\frac 12 \\sigma^2 t_1} \\bigr) - P(t_0, S_0) \\Bigr)^+\n",
    "\\frac{e^{-\\, y^2/2}}{\\sqrt{2 \\pi}} \\, dy\n",
    "$$\n",
    "\n",
    "and the integrand function is explicit.\n",
    "\n",
    "We can numericall evaluate the integral on the right hand side using the function `scipy.integrate.quad` (we suggest to check out the input and output formats of this function on the online documentation).\n",
    "\n",
    "We will use this value as a benchmark value for $I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Benchmark value for the expected loss at time t_1: 5.296\n"
     ]
    }
   ],
   "source": [
    "# Option parameters\n",
    "t_2 = 2.\n",
    "K = 100.\n",
    "S_0 = 100.\n",
    "\n",
    "# Model parameter\n",
    "sigma = 0.3\n",
    "\n",
    "# Intermediate date\n",
    "t_1 = t_2/2\n",
    "\n",
    "def putPrice(t, S, T, K, sigma):\n",
    "    \"\"\"\n",
    "    Black-Scholes price at time t for put option with maturity T.\n",
    "    \"\"\"\n",
    "    deltaT = T - t\n",
    "    sigmaSqrtDeltaT = sigma * np.sqrt(deltaT)\n",
    "    \n",
    "    d = np.log(K/S) / sigmaSqrtDeltaT + sigmaSqrtDeltaT/2\n",
    "    \n",
    "    return K * sps.norm.cdf(d) - S * sps.norm.cdf(d - sigmaSqrtDeltaT)\n",
    "\n",
    "##############################################################\n",
    "# Evaluation of the benchmark value with scipy.integrate.quad \n",
    "##############################################################\n",
    "    \n",
    "##############################################################\n",
    "# TO DO: evaluate the initial price P(0, S_0)\n",
    "\n",
    "price_at_zero = putPrice(t=0, S=S_0, T=t_2, K=K, sigma=sigma)\n",
    "\n",
    "def to_integrate(y):\n",
    "    \n",
    "    S_t_1 = S_0 * np.exp(sigma*np.sqrt(t_1) * y - 0.5*sigma*sigma*t_1)\n",
    "    \n",
    "    loss_t_1 = putPrice(t=t_1, S=S_t_1, T=t_2, K=K, sigma=sigma) - price_at_zero\n",
    "\n",
    "    positivePart = np.maximum(loss_t_1, 0.)\n",
    "\n",
    "    density = np.exp(-0.5*y*y) / np.sqrt(2*np.pi)\n",
    "\n",
    "    return positivePart * density\n",
    "    \n",
    "###################################################\n",
    "# To Do:\n",
    "# + Import the function quad from scipy.integrate\n",
    "# + Evalue the value of the integral I\n",
    "\n",
    "from scipy.integrate import quad\n",
    "\n",
    "benchmark_value = quad(to_integrate, a= -10., b =10.)[0]\n",
    "\n",
    "print(\"\\n Benchmark value for the expected loss at time t_1: %1.3f\" %benchmark_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Nested Monte-Carlo estimator\n",
    "\n",
    "In the definition of $I$, we can recognize a nested expectation of the form\n",
    "\n",
    "$$\n",
    "\\mathbb E \\bigl[ g \\bigl( E[f(X,Y)|X] \\bigr)  \\bigr]\n",
    "$$\n",
    "\n",
    "where $f,g$ are two given functions and $X$ and $Y$ are two independent random variables, since\n",
    "\n",
    "$$\n",
    "I = \\mathbb E \\biggl[ \n",
    "\\biggl( \\mathbb E \\Bigl[ \\Bigl(K - S_{t_1} e^{\\sigma \\, (W_{t_2} - W_{t_1}) - \\frac 12 \\sigma^2 (t_2 - t_1)} \\Bigr)^+ \\Big| S_{t_1} \\Bigr] - P(t_0, S_0) \\biggr)^+\n",
    "\\biggr]\n",
    "$$\n",
    "\n",
    "and the Brownian increment $W_{t_2} - W_{t_1}$ is independent of $S_{t_1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider\n",
    "    \n",
    "+ $M$ i.i.d. samples $S_m$ of $S_{t_1}$\n",
    "\n",
    "    \n",
    "+ $M \\times N$ i.i.d. samples of the Brownian increment , denoted $\\Delta W_{m,n}$\n",
    "\n",
    "The idea behind the __Nested Monte-Carlo estimator__ is the following: for every value $S_m$ sampled at $t_1$, we use the vector of $N$ samples $(\\Delta W_{m,n})_{n = 1, \\dots, N}$ of the Brownian increment to approximate the _inner_ conditional expectation in the expression above.\n",
    "Then, we have to take the outer mean with respect to the simulated values of $S_{t_1}$.\n",
    "\n",
    "The resulting Nested Monte-Carlo estimator is\n",
    "\n",
    "$$\n",
    "I_{M,N} = \\frac 1 M \\sum_{m=1}^M\n",
    "\\biggl(\n",
    "\\frac 1 N \\sum_{n=1}^N\n",
    "\\Bigl(K - S_m \\, e^{\\sigma \\, \\Delta W_{m,n} \\, - \\, \\frac 12 \\, \\sigma^2 (t_2 - t_1)} \\Bigr)^+\n",
    "- P(t_0, S_0)\n",
    "\\biggr)^+.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\blacktriangleright$ Implementation: \n",
    "\n",
    "In the cell below, complete the function `NestedEstimator` so that it outputs one single sample of the estimator $I_{M,N}$.\n",
    "\n",
    "Then, generate a sample of $k$ i.i.d. values of the estimator $I_{M,N}$, choosing\n",
    "\n",
    "+ first $M=N$ (and for example $M=1000$)\n",
    "\n",
    "+ then a smaller value of $N$, of the order of $\\sqrt M$.\n",
    "\n",
    "In both case, give an estimate of the RMSE \n",
    "\n",
    "$$\n",
    "\\sqrt{ \\mathbb E \\Bigl[ \\bigl(I - I_{M,N} \\bigr)^2 \\Bigr]}\n",
    "$$\n",
    "\n",
    "and plot the histogram of the $k$ simulated values of the estimator $I_{M,N}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NestedEstimator(t_1, t_2, K, S_0, sigma, M, N):\n",
    "    \"\"\"\n",
    "    This function outputs one sample of the \n",
    "    Nested MC estimator I_{M,N}\n",
    "    \"\"\"\n",
    "    #######################################\n",
    "    # M i.i.d. samples of S_t_1\n",
    "    G = np.random.normal(loc=0, scale=1, size=M)\n",
    "    \n",
    "    S_1 = S_0 * np.exp(sigma*np.sqrt(t_1)*G - 0.5*sigma*sigma*t_1)\n",
    "    \n",
    "    #######################################\n",
    "    # To Do: simulate M*N i.i.d. samples \n",
    "    # of the Brownian increment Delta W\n",
    "    Delta_W = np.sqrt(t_2 - t_1) * np.random.normal(loc=0, scale=1, size=(M,N))\n",
    "    \n",
    "    ##############################################\n",
    "    # We evaluate the NestedMonte-Carlo estimator\n",
    "    ##############################################\n",
    "    \n",
    "    # This variable will contain the outer sum over {m = 1, ..., M}\n",
    "    outerSum = 0.\n",
    "    \n",
    "    for m in range(M):\n",
    "        # We work with one-dimensional arrays of size N\n",
    "        exponential_factor = np.exp(sigma * Delta_W[m] - 0.5*sigma*sigma*(t_2 - t_1))\n",
    "        \n",
    "        S_2 = S_1[m] * exponential_factor\n",
    "        \n",
    "        ##################################################\n",
    "        ## To Do: implement the inner Monte-Carlo mean 1/N sum_{n = 1, ..., N}\n",
    "        innerMean = 1/N * np.sum(  np.maximum(K - S_2, 0.) )\n",
    "        \n",
    "        ##################################################\n",
    "        ## To do: update the outer sum_{m = 1, ..., M}\n",
    "        outerSum += np.maximum(innerMean - price_at_zero, 0.)\n",
    "    \n",
    "    Nested_estimator = outerSum / M\n",
    "    \n",
    "    return Nested_estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\blacktriangleright$ Simulation of the estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**Nested estimator \n",
      "\n",
      "M = 1000 outer samples and N = 1000 inner samples \n",
      "\n",
      "k = 100 samples of the Nested estimator \n",
      "\n",
      "Estimateur I_{M,N} (1 sample) = 5.295 \n",
      "\n",
      "RMSE = 0.259 \n",
      "\n",
      "Time = 10.69 \n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAD8CAYAAADZu7i7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW9//HPNyNJgAxkgCRAEFBABoGIEwUEZRJBLfSiKFoHsApYqL2CqAxeEe1VKWJxZNCqFFCvgForDvirFRUQBLFARJSAEIaAgYSEJOv3RxIaQkIOcMgO4f16njyes/c6e33PIUc+rLX23uacEwAAALwT4HUBAAAAZzsCGQAAgMcIZAAAAB4jkAEAAHiMQAYAAOAxAhkAAIDHCGQAAAAeI5ABAAB4jEAGAADgsSCvOo6NjXUpKSledY8z0IYNGyRJ5513nseVAN5ZuXLlbudcnNd1APAvzwJZSkqKVqxY4VX3OAN169ZNkvTJJ594WgfgJTP70esaAPgfU5YAAAAeI5ABAAB4jEAGAADgMQIZAACAxwhkAAAAHqs0kJnZLDPLMLN1lbS70MwKzGyg/8oDAACo+Xy57MUcSTMkvVxRAzMLlPSYpPf9UxYAoLpYtWpVr6CgoAnOufpiZgU4UYVmtiM/P39Shw4dKsxJlQYy59ynZpZSSbORkt6QdOEJlQgAqNZWrVrVKzQ0dEZKSkpeWFhYZkBAgPO6JuBMUlhYaDk5OZFbtmyZsWrVqhEVhbJT/peOmSVJulbSsz60HWZmK8xsxa5du061awDAaRYUFDQhJSUlLyIiIocwBpy4gIAAFxERkZOSkpIXFBQ0oaJ2/rhS/zRJ9znnCszsuA2dc89Lel6SUlNT+WKfIS6b+pG27cvxugzt2LxHkpQy9p3T3ldSVJg+G9v9tPcDVHfOufphYWGZXtcBnOnCwsIOFU/7l8sfgSxV0rziMBYrqa+Z5Tvn/s8Px0Y1sG1fjrZMvcrrMtRt+Z8kSZ9UQS1VEfqAM0QAI2PAqSv+HlU4M3nKgcw516TksZnNkbSEMAYAAOA7Xy578bqkzyWdZ2bpZnabmd1pZnee/vIAADh5Ztaxsp8lS5bU8brOU3XXXXclJSQktPW6Dpw8X86yvN7XgznnbjmlagAA8KOlS5f+u+RxdnZ2QP/+/c8dNWrUz/37999fsr19+/beL5LFWc8fa8gAAKiWevTocbDk8f79+wMkqWnTprmlt1ckOzvbwsPDWT+HKsEF/gAAZ73HH388zsw6/vOf/wzv2LHjebVq1erw6KOPxi9cuLCumXVct25daOn27dq1azFgwIAmpbctXry4Tslro6Oj2w0ZMqRRVlZWhX/PPvbYY3FhYWHtS4JiiU8//TTczDr+4x//iJCkV155Jeriiy8+Nzo6ul2dOnUuaN++fYvFixcfd5q15P0cPnz4qO2xsbHtRo0alVh62+zZs6NbtWrVMjQ0tEN8fHzbkSNHJuXn5x/384L/EcgAACg2ZMiQcwYMGJC5YMGCTX379v3F19ctWrSozrXXXts8OTk57+WXX06bNGlS+vvvvx994403Nq7oNTfddFPm4cOHbd68eVGlt7/22msx9evXz7viiisOStLmzZtD+vfvn/nSSy/9MHfu3M1t27Y9eO2115776aefhp/8Oy0yY8aMerfffvs5F1100YHXX3897Z577tkxa9as+DFjxiRW/mr4E1OWAAAUGzFixI4//vGPu0ueL1y4MNiX191///3Jl156adbbb7/9Q8m2hISE/BtuuKHZ2rVrt7dp0ya37GsSExPzL7744qwFCxZEDx8+fG/J9sWLF0dfffXVmQEBRWMmEyZMyCjZV1BQoKuvvvqX9evXhz///POxXbp0+ekk36ry8/M1adKkpMGDB++aPXv21uLNvwQGBrrJkycnT5w4cUdMTEzhyR4fJ4ZABgA4IbfeemvDdevWnfLozMlo3bp19qxZs7ZW3vLkXHfddfsrb3W03bt3B65fvz78qaee2lJ6irBPnz5ZkrR8+fLw8gKZJA0cOHDvfffd13jPnj2B9erVK/jwww8jtm/fHjJkyJAjAW3jxo0h9957b9Ly5cvr7N69O9i5omVtwcHBp7S+7auvvgrbvXt38H/9139llq67V69eWePGjQtYvXp1WPfu3Stdawf/YMoSAIBiycnJJ7x4aufOnUHOOf3+979PCQkJ6VjyExkZ2b6wsFBbt24Nqei1N9544z7nnF5//fUoqWi6Mjk5Obdr167ZknT48GFdddVVzdeuXRv+wAMPbHvnnXc2LFu27LuLL744Kzc39/i3x6lERkZGkCQNGDDg3NJ1d+jQ4XxJ+vHHHyusG/7HCBlQjqSoME+u1s8tm3AmOJ0jVF4re1eCWrVqOUkqG372798fJClXkmJjY/Mlady4cdt69+59zLqzlJSUvIr6i42NLejcufMvCxcujL7rrrv2LFmyJHrQoEF7SvavWrUqLC0trdbixYs39uvXL6tke05OTkBgYGCFI2S1atUqlKRDhw4FBAcHF0pSYWGhDhw4EFiq73xJevrpp7e0bdv2mEt/tGzZstxRPZweBDKgHF6FIm7ZBFQvJWFq7dq1tTp27HhIkr799tvQ9PT00Hbt2h2UpISEhIKWLVtmp6WlhXbp0iX7RPsYNGjQ3lGjRqW89tprURkZGcE33XTTkenKgwcPBkj/CVjFtYSuW7cuol27dgcqOmajRo3yJGn16tW1fvWrX2VL0nvvvVe7dLDs1KlTTnR0dP6PP/4YMmLEiD0VHQtVg0AGAEAFWrdundu8efOciRMnJgcFBSkvL8/+9Kc/NYiMjDxqavPRRx9NHzhwYPOBAwfaddddl1m7du3CzZs3h7z77rtRM2bM2HruuedWOEp2ww037BszZowbM2ZMoyZNmhy66KKLjoxWXXjhhTn16tXLHz16dKMHH3xw2969e4OmTJmSmJCQUOHxJKlnz54HYmJi8keMGNFo/Pjx23fu3Bn89NNPJ4SFhR0JdsHBwZo0aVL6mDFjGu/duzeoV69e+4ODg7Vp06bQRYsWRX3yySebgoKICVWFNWQAABzHvHnzNterV+/wHXfc0eSRRx5JnDhx4rbk5OSjpvMGDBiQ9c4772zYvn178LBhw84ZPHhwsxkzZtRPSUnJTUhIOO66tMjIyMKuXbvu37VrV/C11167t/S+OnXqFM6bNy+tsLBQQ4cObTZ16tTE+++/f3v79u2Pu9g+PDzczZs3L+3w4cN28803N33uuefin3vuuS3h4eFHnTU5cuTIPXPnzv1+1apVETfffHPTm266qemcOXNiL7zwwoMlZ3mialjJ2RpVLTU11a1YscKTvnFiUsa+oy1Tr/K6DHXr1k2S9Mknn3hax+lUXT5rVF9mttI5l1pV/a1Zs2ZLu3btdlfeEkBl1qxZE9uuXbuU8vYRfwEAADxGIAMAAPAYgQwAAMBjBDIAAACPEcgAAAA8RiADAADwGIEMAADAY1yC9wxx2dSPtG3fMbcaqxJJUWGe9AsAwNmCQHaG2LYvhwuGAgBQQ1U6ZWlms8wsw8zWVbB/iJl9U/zzLzNr5/8yAQAAai5f1pDNkdT7OPt/kNTVOddW0sOSnvdDXQAA+MWYMWMSzaxj586dm5fd17t373M6dep0nj/7++abb0LHjBmTuHv37kB/HdOXOn/961+nmFnHSy+99Jj3eeDAAYuIiGhvZh2nT59er/S+goICPfnkk7Ht27dvUbt27fahoaEdmjdvfv6DDz6YsH///tO61vzjjz8OHzNmTGLZ7WPGjEmMjo6usgGeiuqoSpVOWTrnPjWzlOPs/1epp8slJZ96WQCA6qzTI0vbZGTlhlR1v/F1QvO+HH/F2pN57WeffVZ32bJl4V27ds32d12lffvtt7WeeuqpBsOHD98dGxtbcDr7Kis8PLzwyy+/rLt169aghg0bHrmp+fz586PKa19QUKB+/fqd89FHH0UNHTo0Y/z48dtDQkLcypUrw1988cX47du3h7z00ktbT1e9n3/+ecRTTz3V4Mknn9xeevvdd9+969prr913uvr1tY6q5O81ZLdJes/PxwQAVDMZWbkhW6ZetbKq+00Z+07Hk3ldZGRkQUJCQt7DDz/coGvXrt/7u67qokmTJocOHDgQ+Morr0Tff//9u0q2/+1vf4vp0aPHvsWLF8eUbj916tT4999/P/rNN9/ceM0112SVbO/fv3/Wfffdl/HBBx/Ursr6SzRt2vRw06ZND3vRtz8cOHDAateu7U7kNX4bijSzy1UUyO47TpthZrbCzFbs2rWromYAAPiVmbl77733548++ijqyy+/PO6p45s2bQrp16/fOZGRkReEhYW179y5c/M1a9aElm4zbty4+o0aNWodGhraoV69eu1+9atfNf/pp5+ClixZUueGG25oJkktWrRoY2Ydk5KS2pzIsdPS0oK7du3arFatWh2SkpLaPPnkk7En8l6vueaavW+88caR4JWZmRmwbNmyyMGDB+8t23bmzJkJV1555b7SYaxEeHi4GzBgwDHbS/vqq69qdevWrVlERET7iIiI9n369Dnnp59+OjLYk5uba8OGDUtu0KBBm5CQkA7x8fFtr7zyyqaHDh2y6dOn1xs/fnwjSTKzjmbWsWRatuyU5ZIlS+qYWce33367To8ePZqGhYW1b9y4ces333yzbn5+voYPH54cHR3dLj4+vu3EiRMTSte4dOnSiO7duzeLj49vGxYW1r5FixatZs6ceeTzOV4dkrRo0aI6bdu2bVHyZ33jjTc2Kj2VW1LbG2+8Ubd79+7NwsPD2996662Njve5lccvgczM2kp6UdIA59yeito55553zqU651Lj4uL80TUAAD659dZbMxs3bnxo8uTJDSpqs3PnzsAuXbqc9/3339d64oknfpw1a9bmnJycgF69ep134MABk6QZM2bUmz59eoO77rpr55tvvrnxySef/LFJkya5WVlZgZdeeunBhx56KF2S5s6d+/3SpUv/PX/+/DRfj11YWKj+/fs327hxY9i0adO2TJkyZeuzzz4bv2rVKp9HqoYOHbr366+/rr1p06YQSXr11Vej69atm9+7d++jwlVaWlrwtm3bQnr27Ln/xD9Nad26daHdu3dvkZubG/Dss8/+8Mwzz/ywcePGsL59+zYvLCyUJI0fP77+W2+9FXP//fdvf/vttzdOmTJla926dQvy8/M1cODA/XfcccdOSVq6dOm/ly5d+u+ZM2f+eLw+R44c2fjSSy898Oqrr36flJSUN3To0KY333xzo6ysrICXXnrph6uuuipz0qRJyR9++GFEyWs2b94ccskllxx4+umnf/zb3/6WdvXVV2eOGjUq5bnnnouRdNw6Vq5cWevXv/5185iYmPy5c+d+f999921/++23Y/r169e0bG133XVXSps2bbLnzZuXdscdd+w+0c/zlKcszayRpDcl3eSc23iqxwMA4HQIDAzU6NGjd4wePTrlm2++CW3btm1u2TZTpkxJyMnJCVy9evX6hISEAkm64oorDpxzzjltnn766dhx48bt+vLLLyM6d+78y9ixY49M9dx8881H1ju1aNHikCRddNFF2eedd17eiRx7wYIFkd999134hx9++O/u3bsflKRLLrkku1WrVm1SUlKOqbc8HTp0ONS8efOcl19+Ofrhhx/euWDBguh+/fplBgYefY7Bjz/+GCJJKSkpeeUeqBIPPPBAYr169fI//vjjTbVq1XKS1LFjx5wLLrig9fz58yMHDx68f+XKlRHXXHPN3pEjRx4ZrLn99tszJal27dr5Je+pR48eB33pc9CgQXsefvjhnZLUuHHjvNTU1PO///77WsuXL98oSQMGDPhl8eLFMQsXLowqOeawYcMyS15fWFioPn36ZG3bti1k9uzZscOHD9+bmJhYYR0TJkxITExMzFu6dGlaUFBRZIqJicm/4447zlm6dGnEFVdccaR9v379Mv/85z+f9Bo0Xy578bqkzyWdZ2bpZnabmd1pZncWN3lIUj1JfzGz1Wa24mSLAQDgdPrd7363p379+nkVjZItW7asbufOnX+JiYkpOHz4sA4fPqyoqKiC888/P3vlypURknTBBRdkL1u2LHL06NGJH3/8cXh+fn55hzqpY3/xxRcR9erVyy8JY5J07rnn5rVq1cqnwFLiuuuu2/vmm2/G7Ny5M/Bf//pX3SFDhhwzXVnCzE7k0Ed89tlndfr27ZsZGBjoSt5PixYtchMTE3O/+uqrcElq06ZN9oIFC2IfeOCBhC+++CKsZOTsZPXs2fPIKN/555+fK0ldunQ5si0wMFANGzbM3b59+5ETTnbt2hV4yy23NExMTGwTEhLSMSQkpOPrr78e+8MPP9SqrL/Vq1dH9OnTZ19JGJOkW265JTMwMNAtW7bsqFHLq6+++pROQvDlLMvrK9l/u6TbT6UIAACqQnBwsEaNGrVj/PjxDTdu3HjMaEZmZmbQmjVrIkJCQo45eSAgICBLku65557dWVlZAXPnzo2bNm1ag6ioqPyhQ4fueuKJJ7aX/ov7ZI69Y8eOoJiYmGMWs8fGxuYfOHDA58to3HzzzXunTp2a9OCDDzaIj48/3KNHj4NlL2HRuHHjPEnasmXLSZ0tm5mZGTRz5sz6M2fOrF92X3p6eogkTZ069eeAgADNnj07/pFHHkmOj48/PGLEiB0PPvhgxsn0GRMTcyT9lozKRUVFHZWIg4ODXW5u7pGUOXjw4JTVq1fX/sMf/rC9devWh6KiogpmzJgR98EHH5R75mlpu3fvDk5ISDjqzyMoKEhRUVEFe/fuPeoPOzEx0bdkXgGu1A8AOKuMGjVq9xNPPNFg8uTJxwSJyMjI/O7du+dMmDDh53L2FUhFozATJkzImDBhQkZaWlrwrFmz6j322GNJSUlJh//7v/+7wjPWfDl2/fr18/fu3Rtcdv/u3buDSgKIL1q0aJHXpk2bgy+99FLC8OHDd5TXplmzZoeTk5NzP/jgg7pjxow54TVPkZGRBb169cq88847j3ltQkJCvlR0YsC0adO2T5s2bfvatWtDp0+fHvfQQw81bNmy5aGBAwf+cqJ9nqjs7GxbtmxZ1JQpU34q/Wczffp0n4YFY2NjD2dkZByVlfLz87Vv377A0uFQKjpx5FRq5ebiAICzSlhYmLv77rt3LFiwIDYjI+Oo0aEuXbpkbdq0KaxDhw45Xbp0yS79065du2PWcDVr1uzwlClTdjRs2DB3/fr1tSQpNDTUSVJ2dnbAiR67U6dOB/fs2RP00UcfHVmUvmnTppD169eHn+j7HDVq1M7LL79832233VbhyXZ33nlnxj/+8Y/oxYsX1ym7Lzs72xYtWnTM9hKXXnrpLxs2bAjr3Llzdtn3U3rtXIk2bdrkPvfcc+khISFu3bp1YZIUEhJS8lmd3LxpJXJycgIKCgoUGhp6ZK40MzMzYOnSpUeNjlVUR/v27Q+899570aWnpV9++eXogoIC69q16wF/1soIGQDgrDNmzJjd06ZNa/D1119HXHjhhUf+Yh0/fvzON954I6Zz587n3XnnnRkNGzbM+/nnn4OXLVtWp3PnzgeGDx++94YbbmgcHR2df8kllxyMiooq+PDDD+v89NNPoT169MiSpNatWx+SpBkzZsQNGTJkb+3atQs7deqU48uxf/Ob3+yfOHFizpAhQ86ZOHHitrCwsML/+Z//SSw7GuOL22+/PbNkAX1Fxo4dm/HPf/6z9qBBg5oPHTo0o1evXr+EhIS4r7/+OuzFF1+Mv/LKK/f379+/3EtfPPLII9svu+yylpdffnmzW265ZU9cXFz+1q1bg5cuXVr3t7/97Z5+/fplXXnllU3bt2+f3aFDh+zw8PDC+fPnRxcUFFj37t2zJOn8888/VHyshJ49e/4SFRVVUF7wPVn16tUraN26dfbjjz+eGBkZWRAQEKD//d//rV+7du2CgwcPHgnMFdUxYcKEny+55JJWPXv2bDZ8+PBdW7duDZ48eXJy586dfym9oN8fCGQAgBMWXyc072Qv0nqq/frjOHXq1CkcPnz4zsceeyyp9PYGDRrkf/755/++9957k8aPH98wKysrMC4u7nBqauqBjh07ZkvSxRdffGDOnDlxf/3rX+Py8vKsYcOGuU888cSPN9100z6paBH+Qw89lP7CCy/Ez507Nz4hISFv27Zta305dkBAgBYvXpx26623Nr7nnntSYmJiDo8ePfrnDz/8sG55U5mnKjAwUEuWLNn85z//OXbu3Lmxf/3rX+MKCgqscePGhwYNGrTn/vvvr3CtV9u2bXM//fTTf48bNy7x97//fePc3NyA+Pj4vM6dO2e1bNnyUMln9dZbb8XMnDkzobCw0Jo2bZozZ86c77t06ZItSb169TowfPjwnS+88EL8o48+mpSamnrgyy+/3ODP9/j6669vvv322xv/7ne/axIVFZV/2223ZWRnZwfOnj37yPW3KqojNTX10MKFCzc9+OCDSUOHDm0aERFR0L9//73PPPNMuj9rlCRz7pSmPE9aamqqW7GCEzJ9lTL2HW2ZepXXZXiqW7dukqRPPvnE0zpOJ/6cURkzW+mcS62q/tasWbOlXbt2J7y+CMCx1qxZE9uuXbuU8vaxhgwAAMBjBDIAAACPEcgAAAA8RiADAADwGIEMAHA8hYWFhaflGlHA2aT4e1ThvaMIZACACpnZjpycnErv+Qfg+HJycmqZWbl3TZAIZACA48jPz5+0ZcuWkIMHD4YxUgacuMLCQjt48GDYli1bQvLz8ydV1I4LwwIAKtShQ4f3V61aNeL777+f4JyrL/4hD5yoQjPbkZ+fP6lDhw7vV9SIQAYAOK7iv0Qq/IsEwKnjXzoAAAAeI5ABAAB4jEAGAADgMQIZAACAxwhkAAAAHiOQAQAAeIxABgAA4LFKA5mZzTKzDDNbV8F+M7PpZpZmZt+YWQf/lwkAAFBz+TJCNkdS7+Ps7yOpefHPMEkzT70sAACAs0elgcw596mkvcdpMkDSy67IcklRZtbAXwUCAADUdP5YQ5YkaWup5+nF2wAAAOADfwQyK2ebK7eh2TAzW2FmK3bt2uWHrgEAAM58/ghk6ZIalnqeLGl7eQ2dc88751Kdc6lxcXF+6BoAAODM549AtkjS0OKzLS+WtN8597MfjgsAAHBWCKqsgZm9LqmbpFgzS5c0QVKwJDnnnpX0rqS+ktIkZUv67ekqFgAAoCaqNJA5566vZL+TdLffKgIAADjLcKV+AAAAjxHIAAAAPEYgAwAA8BiBDAAAwGMEMgAAAI8RyAAAADxGIAMAAPAYgQwAAMBjBDIAAACPEcgAAAA8RiADAADwGIEMAADAYwQyAAAAjxHIAAAAPEYgAwAA8BiBDAAAwGMEMgAAAI8RyAAAADxGIAMAAPAYgQwAAMBjBDIAAACP+RTIzKy3mW0wszQzG1vO/kZm9rGZfW1m35hZX/+XCgAAUDNVGsjMLFDSM5L6SGol6Xoza1Wm2QOS5jvn2ksaLOkv/i4UAACgpvJlhKyTpDTn3GbnXJ6keZIGlGnjJNUtfhwpabv/SgQAAKjZfAlkSZK2lnqeXryttImSbjSzdEnvShpZ3oHMbJiZrTCzFbt27TqJcgEAAGoeXwKZlbPNlXl+vaQ5zrlkSX0lvWJmxxzbOfe8cy7VOZcaFxd34tUCAADUQL4EsnRJDUs9T9axU5K3SZovSc65zyXVkhTrjwIBAABqOl8C2VeSmptZEzMLUdGi/UVl2vwkqYckmVlLFQUy5iQBAAB8EFRZA+dcvpmNkPS+pEBJs5xz35rZZEkrnHOLJP1B0gtmNlpF05m3OOfKTmsCqERSVJhSxr7jSb+fje1e5f0CAIpUGsgkyTn3rooW65fe9lCpx+slXebf0oCzj1ehyIsQCAD4D67UDwAA4DECGQAAgMcIZAAAAB4jkAEAAHiMQAYAAOAxAhkAAIDHCGQAAAAeI5ABAAB4jEAGAADgMQIZAACAxwhkAAAAHiOQAQAAeIxABgAA4DECGQAAgMcIZAAAAB4jkAEAAHiMQAYAAOAxAhkAAIDHCGQAAAAeI5ABAAB4zKdAZma9zWyDmaWZ2dgK2vzGzNab2bdm9pp/ywQAAKi5giprYGaBkp6RdKWkdElfmdki59z6Um2aSxon6TLnXKaZxZ+ugiXpsqkfadu+nNPZRYWSosL02djunvQNAABqpkoDmaROktKcc5slyczmSRogaX2pNndIesY5lylJzrkMfxda2rZ9Odoy9arT2UWFUsa+40m/AACg5vJlyjJJ0tZSz9OLt5V2rqRzzewzM1tuZr39VSAAAEBN58sImZWzzZVznOaSuklKlvT/zKy1c27fUQcyGyZpmCQ1atTohIsFAACoiXwZIUuX1LDU82RJ28tp87Zz7rBz7gdJG1QU0I7inHveOZfqnEuNi4s72ZoBAABqFF8C2VeSmptZEzMLkTRY0qIybf5P0uWSZGaxKprC3OzPQgEAAGqqSgOZcy5f0ghJ70v6TtJ859y3ZjbZzPoXN3tf0h4zWy/pY0l/dM7tOV1FAwAA1CS+rCGTc+5dSe+W2fZQqcdO0pjiHwAAAJwArtQPAADgMQIZAACAxwhkAAAAHvNpDRn+IykqzJOr9SdFhVV5nwAAoGoQyE4Q97EEAAD+xpQlAACAxwhkAAAAHiOQAQAAeIxABgAA4DECGQAAgMcIZAAAAB4jkAEAAHiMQAYAAOAxAhkAAIDHCGQAAAAeI5ABAAB4jEAGAADgMQIZAACAxwhkAAAAHiOQAQAAeIxABgAA4DGfApmZ9TazDWaWZmZjj9NuoJk5M0v1X4kAAAA1W6WBzMwCJT0jqY+kVpKuN7NW5bSrI2mUpC/8XSQAAEBN5ssIWSdJac65zc65PEnzJA0op93Dkh6XdMiP9QEAANR4vgSyJElbSz1PL952hJm1l9TQObfEj7UBAACcFXwJZFbONndkp1mApKck/aHSA5kNM7MVZrZi165dvlcJAABQg/kSyNIlNSz1PFnS9lLP60hqLekTM9si6WJJi8pb2O+ce945l+qcS42Lizv5qgEAAGoQXwKW0LPWAAAJEklEQVTZV5Kam1kTMwuRNFjSopKdzrn9zrlY51yKcy5F0nJJ/Z1zK05LxQAAADVMpYHMOZcvaYSk9yV9J2m+c+5bM5tsZv1Pd4EAAAA1XZAvjZxz70p6t8y2hypo2+3UywIAADh7cKV+AAAAjxHIAAAAPEYgAwAA8BiBDAAAwGMEMgAAAI8RyAAAADxGIAMAAPAYgQwAAMBjBDIAAACPEcgAAAA8RiADAADwGIEMAADAYwQyAAAAjxHIAAAAPEYgAwAA8BiBDAAAwGMEMgAAAI8RyAAAADxGIAMAAPBYkNcFAPBeUlSYUsa+41nfn43t7knfAFBdEMgAeBqIvAqCAFCd+DRlaWa9zWyDmaWZ2dhy9o8xs/Vm9o2ZfWhmjf1fKgAAQM1UaSAzs0BJz0jqI6mVpOvNrFWZZl9LSnXOtZW0UNLj/i4UAACgpvJlhKyTpDTn3GbnXJ6keZIGlG7gnPvYOZdd/HS5pGT/lgkAAFBz+RLIkiRtLfU8vXhbRW6T9F55O8xsmJmtMLMVu3bt8r1KAACAGsyXQGblbHPlNjS7UVKqpD+Vt98597xzLtU5lxoXF+d7lQAAADWYL2dZpktqWOp5sqTtZRuZ2RWSxkvq6pzL9U95AAAANZ8vI2RfSWpuZk3MLETSYEmLSjcws/aSnpPU3zmX4f8yAQAAaq5KA5lzLl/SCEnvS/pO0nzn3LdmNtnM+hc3+5Ok2pIWmNlqM1tUweEAAABQhk8XhnXOvSvp3TLbHir1+Ao/1wUAAHDW4F6WAAAAHuPWSQDOSpdN/Ujb9uV40jf37wRQFoEMwFlp274cbZl6lSd9c/9OAGUxZQkAAOAxAhkAAIDHCGQAAAAeI5ABAAB4jEAGAADgMQIZAACAxwhkAAAAHiOQAQAAeIxABgAA4DGu1A/AU0lRYZ5cuT4pKqzK+wSAihDIAHiKezoCAFOWAAAAniOQAQAAeIxABgAA4DECGQAAgMcIZAAAAB4jkAEAAHjMp0BmZr3NbIOZpZnZ2HL2h5rZ34r3f2FmKf4uFAAAoKaqNJCZWaCkZyT1kdRK0vVm1qpMs9skZTrnmkl6StJj/i4UAACgpvJlhKyTpDTn3GbnXJ6keZIGlGkzQNLc4scLJfUwM/NfmQAAADWXL4EsSdLWUs/Ti7eV28Y5ly9pv6R6/igQAACgpvPl1knljXS5k2gjMxsmaVjx0wNmtsGH/ssvqvJJ0VhJu0/2+NXAmVz/aa39NA++nsmfu3Rm139W1e7D/8Mq0vikXwmg2vIlkKVLaljqebKk7RW0STezIEmRkvaWPZBz7nlJz59cqSfGzFY451Kroq/T4Uyun9q9cybXT+0Azma+TFl+Jam5mTUxsxBJgyUtKtNmkaSbix8PlPSRc+6YETIAAAAcq9IRMudcvpmNkPS+pEBJs5xz35rZZEkrnHOLJL0k6RUzS1PRyNjg01k0AABATeLLlKWcc+9KerfMtodKPT4kaZB/SztlVTI1ehqdyfVTu3fO5PqpHcBZy5hZBAAA8Ba3TgIAAPBYjQhkZhZoZl+b2ZJy9jUys4+L939jZn29qLEiZrbFzNaa2WozW1HOfjOz6cW3pfrGzDp4UWd5fKh9SHHN35jZv8ysnRd1lqey2ku1u9DMCsxsYFXWdzy+1G5m3Yr3f2tmy6q6xuPx4fcm0swWm9ma4vp/60Wd5TGzKDNbaGb/NrPvzOySMvur7fcVQPXm0xqyM8A9kr6TVLecfQ9Imu+cm1l8y6d3JaVUYW2+uNw5V9E1jPpIal78c5GkmcX/rS6OV/sPkro65zLNrI+K1tmcKbWX3DbsMRWd0FLdVFi7mUVJ+ouk3s65n8wsvmpL88nxPvu7Ja13zl1tZnGSNpjZq8V3CvHanyX93Tk3sPis8/Ay+6v79xVANXXGj5CZWbKkqyS9WEETp/8EtUgdew216m6ApJddkeWSosysgddF+cI59y/nXGbx0+UquobdmWSkpDckZXhdyAm6QdKbzrmfJMk5d6bV7yTVKb79Wm0Vnbmd721JkpnVldRFRWeVyzmX55zbV6bZGft9BeCtMz6QSZom6b8lFVawf6KkG80sXUWjYyOrqC5fOUn/MLOVxXcyKMuXW1d5pbLaS7tN0ntVUJOvjlu7mSVJulbSs1VeWeUq+9zPlRRtZp8UtxlaxfVVprL6Z0hqqaJ/PK2VdI9zrqLvd1U6R9IuSbOLl0C8aGYRZdpU5+8rgGrsjA5kZtZPUoZzbuVxml0vaY5zLllSXxVdL606ve/LnHMdVDTVcbeZdSmz36fbUnmkstolSWZ2uYoC2X1VWVwlKqt9mqT7nHMFVV9apSqrPUhSRxWNHPeS9KCZnVvFNR5PZfX3krRaUqKkCyTNKB6d8lqQpA6SZjrn2ks6KGlsmTbV+fsKoBqrTsHkZFwmqb+ZbZE0T1J3M/trmTa3SZovSc65zyXVUtF956oF59z24v9mSHpLUqcyTXy5dZUnfKhdZtZWRdPJA5xze6q2wor5UHuqpHnFv1sDJf3FzK6p0iIr4OPvzN+dcweL12l9KqnanFDhQ/2/VdGUq3POpaloLWKLqq2yXOmS0p1zXxQ/X6iigFa2TbX8vgKo3s7oQOacG+ecS3bOpajo7gAfOeduLNPsJ0k9JMnMWqookO2q0kIrYGYRZlan5LGknpLWlWm2SNLQ4rO3Lpa03zn3cxWXegxfajezRpLelHSTc25j1VdZPl9qd841cc6lFP9uLZR0l3Pu/6q82DJ8/J15W9KvzCzIzMJVtKj8u6qttHw+1l/6O5sg6TxJm6uyzvI453ZI2mpm5xVv6iFpfZlm1fL7CqD6qylnWR7Fjr6t0x8kvWBmo1U0dXBLNbrPZoKkt4rWLitI0mvOub+b2Z2S5Jx7VkXr3vpKSpOUraLRg+rAl9ofklRPRaNLkpRfTW7A7Evt1VWltTvnvjOzv0v6RkVrK190zpUNPV7x5bN/WNIcM1uroinA+453NmwVGynp1eIzLDdL+u0Z8n0FUM1xpX4AAACPndFTlgAAADUBgQwAAMBjBDIAAACPEcgAAAA8RiADAADwGIEMAADAYwQyAAAAjxHIAAAAPPb/AbSXtFm2tiQEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "################################################\n",
    "# Parameters M et N of the Nested MC estimator\n",
    "M = int(1.e3)\n",
    "N = M\n",
    "\n",
    "################################################\n",
    "# We generate a sample of k i.i.d. values tirages\n",
    "# of the estimator \n",
    "k = 100\n",
    "\n",
    "sample_nested_estim = np.zeros(k)\n",
    "\n",
    "from time import time\n",
    "time0 = time()\n",
    "\n",
    "##########################################\n",
    "## To Do: generate the samples of I_{M,N}\n",
    "for i in range(k):\n",
    "    sample_nested_estim[i] = NestedEstimator(t_1, t_2, K, S_0, sigma, M, N)\n",
    "\n",
    "time1 = time()\n",
    "\n",
    "###########################################\n",
    "# To Do: estimate the RMSE from the sample\n",
    "# of k values of I_{M,N} generated above\n",
    "\n",
    "MSE = np.mean( (sample_nested_estim - benchmark_value)**2 )\n",
    "\n",
    "RMSE = np.sqrt(MSE)\n",
    "\n",
    "#######################################\n",
    "# Some printing of the results\n",
    "print(\"**Nested estimator \\n\")\n",
    "print(\"M = %1.0f outer samples and N = %1.0f inner samples \\n\" %(M,N) )\n",
    "print(\"k = %1.0f samples of the Nested estimator \\n\" %k )\n",
    "print(\"Estimateur I_{M,N} (1 sample) = %1.3f \\n\" %sample_nested_estim[0] )\n",
    "print(\"RMSE = %1.3f \\n\" %RMSE )\n",
    "print(\"Time = %1.2f \\n\" %(time1-time0) )\n",
    "\n",
    "#####################################\n",
    "# To Do: histogram of the estimator\n",
    "\n",
    "plt.hist( sample_nested_estim , density=True, histtype='step', bins=int(np.sqrt(k)), label=\"Nested MC estimator\")\n",
    "\n",
    "plt.axvline(benchmark_value, linewidth=1.5, color='k', label=\"True value\")\n",
    "\n",
    "plt.legend(loc=9, fontsize=15, bbox_to_anchor=(1.4, 1.0), ncol=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\blacktriangleright$ Remark (pricing measure and historical measure)\n",
    "\n",
    "The estimator above can also be applied to the evaluation of the (proper) expected loss\n",
    "\n",
    "$$\n",
    "\\mathbb E^{\\mathbb P^{\\mathrm{hist}}} \\Bigl[\\Bigl(P(t_1, S_{t_1}) - P(t_0, S_0) \\Bigr)^+ \\Bigr]\n",
    "$$\n",
    "\n",
    "where $\\mathbb P^{\\mathrm{hist}}$ is the physical measure of the market, under which we assume that the underlying spot price follows the dynamics\n",
    "\n",
    "$$\n",
    "S_{t_i} = S_0 \\, e^{\\sigma \\, W_{t_i} \\, + \\, (\\mu - \\frac 12 \\sigma^2 ) t_i},\n",
    "\\hspace{30mm} (1)\n",
    "$$\n",
    "\n",
    "where $\\mu$ is a drift parameter.\n",
    "\n",
    "Note that we still have a nested expectation problem, the difference being that now the outer expectation $\\mathbb E^{\\mathbb P^{\\mathrm{hist}}}$ and the inner expectation $P(t_1, S_{t_1}) = \\mathbb E \\bigl[ (K - S_{t_2})^+ \\big| S_{t_1} \\bigr]$ are computed according to different probability measures.\n",
    "\n",
    "$\\blacktriangleright$ __Question__: How would you modify the estimator and the code above so to take the historical model $(1)$ into account? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 2. Dynamic programming with empirical regression: the example of the Bermudan Put"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\blacktriangleright$ Setting\n",
    "\n",
    "In the framework of the discrete-time Black-Scholes model of the previous exercise\n",
    "\n",
    "$$\n",
    "S_i = S_0 \\, e^{\\sigma \\, W_i - \\frac 12 \\sigma^2 t_i},\n",
    "\\qquad t_0 < t_1 < \\dots < t_n\n",
    "$$\n",
    "\n",
    "we wish to evaluate the price of a Bermudan put option of maturity $T=t_n$, which \n",
    "\n",
    "+ can be exerced by its owner at any date $t_i$ previous to $T$\n",
    "\n",
    "\n",
    "+ if exercised at $t_i$, generates the gain $(K - S_i)^+$.\n",
    "\n",
    "This pricing problem for such an option is an example within the more general family of _optimal stopping problems_ (the owner looks for an optimal exercise time)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\blacktriangleright$ The recursive dynamic programming equation\n",
    "\n",
    "The theory of optimal stopping allow to represent the price $Y_i$ of the Bermudan option at any date $t_i$ in terms of the corresponding _dynamical programming equation_: starting from the final value $Y_n = (K - S_n)^+$, we have to evaluate recursively\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "&Y_n = (K - S_n)^+\n",
    "\\\\\n",
    "&Y_i = \\max \\left( (K - S_i)^+, \\mathbb E \\bigl[Y_{i+1} \\big| S_i \\bigr] \\right)\n",
    "\\qquad i = n-1, \\dots, 0.\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "The equation above\n",
    "\n",
    "+ Requires to evaluate the conditional expectation $\\mathbb E[Y_{i+1} | S_i]$ ($Y_{i+1}$ being known from the evaluation already performed at date $t_{i+1}$).\n",
    "\n",
    "\n",
    "+ is arguably still simpler to implement than the representation of the option price as\n",
    "$$\n",
    "Y_0 = \\sup_{\\tau \\in \\mathcal{T}_n } {\\mathbb E}\\bigl[ (K - S_\\tau)^+ \\bigr],\n",
    "$$\n",
    "where $\\mathcal{T}_n$ is the set of all possible (random) stopping times with values in $\\{0,t_1,\\dots, t_n\\}$.\n",
    "\n",
    "\n",
    "Note that the dynamic programming for optimal stopping is a specific example of more general recursive equations of the form\n",
    "\n",
    "$$\n",
    "Y_i =  f\\left( t_i, S_i, \\mathbb E \\bigl[Y_{i+1} | S_i \\bigr] \\right)\n",
    "$$\n",
    "\n",
    "where $f$ is a given function, that appear when simulating _backward_ stochastic differential equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\blacktriangleright$ The pricing problem in a simple two-period framework\n",
    "\n",
    "When $n=2$, the Bermudan put price is given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Y_0 &= \\max \\left\\{ (K-S_0)^+, \\mathbb E \\left[ Y_1 \\right] \\right\\}\n",
    "\\\\\n",
    "&= \\max \\left\\{ (K-S_0)^+, \\mathbb E \\left[ \\max \\Bigl( (K - S_1)^+, \\mathbb E \\left[ (K - S_2)^+ |S_1 \\right] \\Bigr) \\right] \\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The problem is analogous to one considered in the previous exercise: we are in front of a nested expectation of the form $\\mathbb E \\bigl[ g \\bigl(X, \\mathbb E[f(X,Y)|X] \\bigr)  \\bigr]$ (now for a different expression of the outer function $g$). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Benchmark price\n",
    "\n",
    "Exactly as in the previous exercise, in the present two-period setting we have a semi-explicit formula for the Bermudan option price, which we can use a reference value for other estimators.\n",
    "\n",
    "Here again, we exploit the explicit knowledge of the function (the European put price)\n",
    "\n",
    "$$\n",
    "\\mathbb E \\left[ (K - S_2)^+ |S_1 \\right] = P(t_1, S_1)\n",
    "$$\n",
    "\n",
    "which allows to define the estimator\n",
    "\n",
    "$$\n",
    "Y_0^{M} = \\max \\Bigl\\{\n",
    "(K-x_0)^+, \\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, P(t_1, S_1^m) \\right) \\Bigr\\}\n",
    "$$\n",
    "\n",
    "based on a sample $(S_{1}^m)_{m = 1, \\dots, M}$ of i.i.d. values of $S_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cells below:\n",
    "\n",
    "+ For future use of empirical regression and for visualisation purposes, we rewrite the put price as a function of the Brownian motion $W_1$ at time $t_1$ rather than the spot price $S_1$ (the one being a function of the other): we define\n",
    "$$\n",
    "v_1(W_1) := P(t_1, S_1) = P \\Bigl( t_1, S_0 \\, e^{\\sigma \\, W_1 \\, - \\, \\frac 12 \\sigma^2 t_1} \\Bigr)\n",
    "$$\n",
    "so that the estimator above reads\n",
    "$$\n",
    "Y_0^{M} = \\max \\Bigl\\{\n",
    "(K-x_0)^+, \\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, v_1(W_1^m) \\right) \\Bigr\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Problem parameters\n",
    "K = 1.2\n",
    "S_0 = 1.\n",
    "sigma = 0.2\n",
    " \n",
    "t_2 = 2.\n",
    "t_1 = t_2 / 2.\n",
    " \n",
    "M = int(5e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Benchmark price= 0.2455 +/- 0.0039 \n",
      "\n",
      "Erreur relative (TCL) = 0.016 \n",
      "\n",
      "Time: 0.0030 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "def putBlackScholes(S_0, t_1, t_2, K, W_1, sigma):\n",
    "    \"\"\"\n",
    "    Black Scholes price at t_1 of a put option with maturity t_2,\n",
    "    as a function of the underlying Brownian motion\n",
    "    at time t_1.\n",
    "     \n",
    "    W_1: numpy array\n",
    "    \"\"\"\n",
    "    S_1 = S_0 * np.exp(sigma*W_1 - 0.5*sigma*sigma*t_1)\n",
    "    \n",
    "    price = putPrice(t=t_1, S=S_1, T=t_2, K=K, sigma=sigma)\n",
    "        \n",
    "    return price\n",
    "\n",
    "###################################################\n",
    "# The function below is useful to obtain the values\n",
    "# of S at time t once the values of W_t are given.\n",
    "#\n",
    "# GBM stands for \"Geometric Brownian Motion\"\n",
    "###################################################\n",
    "def GBM(S_0, t, sigma, W):\n",
    "    return S_0 * np.exp(-0.5*sigma*sigma*t + sigma * W)\n",
    " \n",
    "time0 = time()\n",
    " \n",
    "#################################################\n",
    "# To Do: generate a sample of M values of W_1\n",
    "# and the corresponding sample of values of S_1\n",
    "\n",
    "W_1 = np.sqrt(t_1) * np.random.randn(M) \n",
    "S_1 = GBM(S_0, t_1, sigma, W_1)  \n",
    "\n",
    "################################################\n",
    "\n",
    "time0_1 = time()\n",
    " \n",
    "timeSimulations = time0_1 - time0\n",
    " \n",
    "###################################################\n",
    "### 1. Prix Benchmark\n",
    "###################################################\n",
    " \n",
    "time1 = time()\n",
    " \n",
    "v_1 = putBlackScholes(S_0, t_1, t_2, K, W_1, sigma)\n",
    " \n",
    "################################################\n",
    "## TO DO: evalute the benchmark price Y_0^M\n",
    "\n",
    "def pos_part(x):\n",
    "    return np.maximum(x, 0)\n",
    "\n",
    "benchmark_price = np.mean(  np.maximum( pos_part(K -  S_1), v_1) ) \n",
    "\n",
    "benchmark_price = np.maximum(benchmark_price, 0)\n",
    "\n",
    "var_TCL = np.var(np.maximum( pos_part(K -  S_1), v_1))  # update here\n",
    "################################################\n",
    " \n",
    "radiusIC = 1.96 * np.sqrt(var_TCL/M)\n",
    " \n",
    "time2 = time()\n",
    " \n",
    "print(\"Benchmark price= %1.4f +/- %1.4f \\n\" %(benchmark_price, radiusIC) )\n",
    "print(\"Erreur relative (TCL) = %1.3f \\n\" %(radiusIC / benchmark_price) )\n",
    "print(\"Time: %1.4f \\n\" %(time2 - time1 + timeSimulations) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Price by empirical regression\n",
    "\n",
    "$\\blacktriangleright$ In this question, we proceed as if we did not know explicitly the conditional expectation $\\mathbb E \\left[ (K - S_2)^+ \\big|S_1 \\right]$ (which is the realistic case in practice, for more general stochastic models and more complex portfolios).\n",
    "\n",
    "$\\blacktriangleright$ We therefore have to _approximate_ the conditional expectation, which we can do applying __least square linear regression__ to a simulated i.i.d. sample of the couple $(W_1, (K - S_2)^+)$.\n",
    "\n",
    "\n",
    "$\\blacktriangleright$ We will therefore need to simulate a sample $(W_1^m, S_2^m)_{1 \\le m \\le M}$ of i.i.d. values of the couple $(W_1, S_2)$.\n",
    "\n",
    "One can note that, in this part, we work with $M$ samples in total - as opposed to the $M \\times N$ samples of the Nested Monte-Carlo procedure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Choice of the regression space\n",
    "\n",
    "+ We are going to build an empirical approximation $\\tilde v_1$ of the function\n",
    "$$\n",
    "x \\mapsto v_1(x) = \n",
    "P \\left(T_1, e^{\\sigma \\, x \\, - \\, \\frac12 \\sigma^2 T_1 } \\right)\n",
    "$$\n",
    "on the state space of the Brownian motion $W_1$. \n",
    "\n",
    "\n",
    "+ We choose as an approximation space the space generated by the indicator functions of disjoint intervals $I_k$\n",
    "$$\n",
    "\\varphi_k(x) = 1_{I_k}(x),\n",
    "\\quad\n",
    "\\quad k = 1, \\dots, n,\n",
    "$$\n",
    "where\n",
    "$$\n",
    "I_k = \\left[-a + (k-1) \\delta, -a + k \\delta \\right[, \\qquad \\delta = \\frac{2a}n\n",
    "$$\n",
    "so that the union of the $I_k$'s is the (large) interval $[-a,a]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The regression problem for $v(x) = \\mathbb E \\left[ (K - S_2)^+ \\big|W_1 = x \\right]$\n",
    "\n",
    "+ As stated above, we approximate the function $v_1$ using the elements of the vector space\n",
    "$$\n",
    "\\Phi = \\text{Vect}(\\varphi_1,\\dots, \\varphi_n)\n",
    "= \\left\\{\\sum_{k=1}^n \\alpha_k 1_{I_k}(\\cdot): \\alpha_1, \\dots, \\alpha_n \\in \\mathbb{R} \\right\\}.\n",
    "$$\n",
    "that is, using piece-wise constant functions of the form $\\sum_{k=1}^n \\alpha_k 1_{I_k}(\\cdot)$.\n",
    "\n",
    "\n",
    "+ The __least square regression coefficients__ for the function $v(x) = \\mathbb E \\left[ (K - S_2)^+ |W_1 = x \\right]$ are defined by the minimization problem\n",
    "$$ \\label{e:regrEmp}\n",
    "\\begin{aligned}\n",
    "\\alpha^* = \\underset{\\alpha\\in \\mathbb{R}^n}{\\rm arg\\min}\n",
    "\\ \\frac 1 M \\sum_{m = 1}^M\n",
    "\\left(\n",
    "\\bigl(K - S_2^m \\bigr)^+ - \\sum_{k=1}^n \\alpha_k 1_{I_k}(W_1^m)\n",
    "\\right)^2.\n",
    "\\end{aligned}\n",
    "$$\n",
    "Once the $\\alpha^*_k$ are known, the regression function (that is, the empirical approximation of $v_1(\\cdot)$) is given by\n",
    "$$\n",
    "\\tilde v_1(\\cdot) = \\sum_{k=1}^{n} \\alpha_k^* \\ 1_{I_k}(\\cdot)\n",
    "$$\n",
    "\n",
    "\n",
    "+ Finally, once the regression function $\\tilde v_1$ has been evaluated, the new estimator of the Bermudan put is\n",
    "$$\n",
    "\\tilde{Y}_0^{M} = \n",
    "\\max \\Bigl\\{\n",
    "(K - x_0)^+, \\frac 1M \\sum_{m=1}^M \\max \\left( (K - S_{1}^m)^+, \\tilde{v}_1(W_1^m) \\right)\n",
    "\\Bigr\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\blacktriangleright$ __Question__: do we have an explicit expression for the regression coefficients $\\alpha^*_k$? \n",
    "\n",
    "__Answer__: in the case of disjoint intervals $I_k$, the expression of the optimal regression coefficients $\\alpha^*$ is relatively simple: denoting\n",
    "\n",
    "$$\n",
    "N_k = \\#\\{m: \\ W_1^m \\in I_k \\}\n",
    "$$\n",
    "\n",
    "the number of samples  $W_1^m$ contained in the cell $I_k$, we have\n",
    "\n",
    "$$\n",
    "\\alpha_k^* = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\\frac 1 {N_k}\n",
    "\\sum_{m = 1}^M\n",
    "\\left(K - S_2^m \\right)^+ 1_{\\{W_1^m \\in I_k\\}}\n",
    "& \\quad \\mbox{if } N_k \\neq 0\n",
    "\\\\\n",
    "0 & \\quad \\mbox{otherwise}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "The default choice $\\alpha_k^*=0$ corresponds to the vector of coefficients $\\alpha^*$ of minimal norm in $\\mathbb R^n$.\n",
    "\n",
    "In the cell below, evaluate\n",
    "\n",
    "+ the coefficients $\\alpha^*_k$ in the Python function `EmpiricalRegression`\n",
    "+ the Bermudan price estimator $\\tilde Y_0^{M}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Price by empirical regression = 0.2483\n",
      "Time: 0.0020 \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x9074c50>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEVCAYAAAAb/KWvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW5x/HPwwAJQtgEKhJCuAIiRAIaQaUqrgW8grYuQBWpetEqrdbaq6IFl1ZqXXCjrbi2vQUX3LBiXeqC0qIgRiQgsggY0CJrANnzu3/8JjKECRmSmZxZvu/Xa16ZOefMzDMYn/zmd37necw5h4iIpJd6QQcgIiLxp+QuIpKGlNxFRNKQkruISBpSchcRSUNK7iIiaUjJXUQkDSm5i4ikISV3EZE0VD+oN27VqpXLz88P6u1FRFLSRx99tMY517q64wJL7vn5+cyePTuotxcRSUlmtjyW4zQtIyKShpTcRUTSkJK7iEgaUnIXEUlDSu4iImlIyV1EJA0puYuIpKGUS+6ffQa//jVs2xZ0JCIiySuwi5hq6uWX4Te/gWeegYkT4aSTgo5IJP1s376ddevWsWnTJnbv3h10OGkrFAqRk5NDy5YtycrKiutrp1xy/9WvoGdPuPxy6NcPLrsM7roLmjcPOjKR9LB9+3ZWrFhBixYtyM/Pp0GDBphZ0GGlHeccO3fupKysjBUrVpCXlxfXBJ9y0zIAp58O8+b5RP/EE3DEEfD880FHJZIe1q1bR4sWLWjVqhUNGzZUYk8QM6Nhw4a0atWKFi1asG7duri+fkomd4CDDoLf/x4+/BAOOQR+9CN/+/rroCMTSW2bNm2iadOmQYeRUZo2bcqmTZvi+popm9wrHHWUT/DjxsErr0C3bvDnP4NzQUcmkpp2795NgwYNgg4jozRo0CDu5zZSPrkDNGgAN9wAn3zik/uIETBgAKxYEXRkIqlJUzF1KxH/3mmR3CscfjhMnw4PPgjvvw8FBX5FjUbxIpJpUm61DNdcA8XFVe6uB4wCLi2AhQthw+Uw93qf+LOzKx08bBiMHJnIaEVEApFWI/dIjbKhsBC6dIayMpg1C1atgu8G8cXFMGlSkCGKiCRM6o3c77sv5kMNOBTYudyvh3/zTTitAzz6KHS4uF+iIhQRCVzajtwjdegAr78Of/oTzJwJRx4JX30VMYoXEUkzGZHcAcz8Va2ffgpHHw0LP/f3V60KOjIRkfjLmOReIT8f/vlP6NQJNmzwK2qeeiroqEQk2YwbNw4z48EHH4y6f+nSpWRlZdG7d29cEi7Jy7jkDlCvHuS2g6Iiv4pm6FC44AJYuzboyEQkWRQWFgIwb968qPuvv/56duzYwfjx45PyuoCMTO4VDmoE770Hv/0tvPCCn4v/xz+CjkpEkkGPHj0AKCkp2WffjBkzmDJlCueffz59+/at69BiEtNqGTPrD9wPhIBHnXO/i3LM+cAt+POUnzjnhsUxzoSpXx9Gj4aBA+Gii/yVrVdcAXffDY0bBx2dSHKp5jKTpNCz5wEtqqtSbm4uLVu23Ce5O+e49tprycrK4s4776z9GyVItSN3MwsBE4ABQDdgqJl1q3RMZ+BGoK9zrjtwTQJiTaiePf1a+F/+Eh5+GHr18o9FJHP16NGDDRs2UFpa+t22yZMn8+GHH3LNNdeQn58PwPTp0xk0aBDt2rXDzHjyySeDCThCLCP33sBi59xSADN7ChgMzI845n+ACc659QDOudXxDjQhiot9UfiwbOBu4OYesOAz+LY3LMuHDnl+tU2VdKWrZIh4jIhTSWFhIe+88w4lJSXk5uaybds2Ro8eTZs2bRg9evR3x23evJmCggKGDx/O8OHDA4x4j1jm3NsBX0Y8Lg1vi9QF6GJmM8xsZngaJ7kNG+aH61E0bw7HFEHrNrBsGXxcDFurauunK11F0lbFvHvFSdX77ruP5cuXc/vtt+9VFnngwIHccccdnHvuudSrlxynMmMZuUcbs1Ze91Mf6Az0A3KB98yswDm3Ya8XMhsJjATIy8s74GDjauTI/Y626+PnoD6ZDGf/FMrnw4QJcOGFlUbxESN/EUkvkStmVq9ezbhx4ygoKODSSy8NOLLqxfInphRoH/E4F6h86U8p8JJzbqdz7gtgIT7Z78U5N9E5V+ScK2rdunVNY65TQ4f6UsI9e8Lw4X7Av2FD9c8TkdTXvXt3QqEQJSUljB07lrKyMu69915CoVDQoVUrluQ+C+hsZh3NrCEwBJha6ZgXgZMBzKwVfppmaTwDDVKHDvD2237J5LPP+kQ/Y0bQUYlIomVnZ9OlSxfmzp3LI488wplnnsnpp58edFgxqTa5O+d24avovgYsAJ5xzpWY2W1mNih82GvAWjObD7wN/Mo5l1aXBIVCfsnkjBn+/oknwi23qFa8SLorLCxk+/btmBl333130OHELKZ17s65acC0StvGRNx3wLXhW1rr0wc+/hhGjYJbb4XBzXyD7sql4kUkPUyePJnJkycHHcYBS47TuimmaVP4y1/gr3+FzZtg9mx4/vmgoxKRoGzevJni4mKKi4spLy9nxYoVFBcXsyLAXp9K7rVw4YW+Pk2jbPjRj+Cqq2BbVUsmRSRtzZ49m169etGrVy+2bt3K2LFj6dWrF2PGjKn+yQmSes06kkyjRtDrKPjlMXDPPb5369NPQ9euQUcmInWlX79+SVcZUiP3OKhnvhbNK6/4+vBHH+2nbUREgqLkHkcDB/oLVo85Bi6+GH7yE9iyJeioRCQTKbnHWbt2vlfrzTfDn//sE32UiqEiIgml5J4A9evD7bf7vq1r10Lv3j7Ri4jUFSX3BDrtND9N07s3jBgBl14K334bdFQikgmU3BOsbVt44w0/TfPEE3DssfD550FHJSLpTsm9DlRM00yb5lfTFBX5GjUiIomi5F6H+vf3pQsKCuD8833Lsh07go5KRNKRknsda98e3nkHrr4a7r8fTj4ZVq4MOioRSTdK7gFo2NC3K3v6aZg71/drfeutoKMSkXSi5B6g88/3TbhbtYLTT4e77lIJYRGJDyX3gHXtCh984AuP/e//wrnnQllZ0FGJSKpTck8COTl+iuaee+Cll3zN+M8+CzoqEUllSu5JwgyuvdaXLqi4qvWFF4KOSkQASktLGT9+PK+++ioA5eXl3HLLLZSWlgYcWdVU8jceiouhX7/av86wYfQbOZI5c/w0zQ9/CDfe6NfIp0A/XpG0VF5ezimnnELDhg1ZsmQJAwYM4Pjjj+fWW2/lkksuCTq8Kim519awYfF5neJi/3PkSHJzYfp038pv3Di/629/gxYt4vNWIhK7+fPns2jRItauXcuiRYs4++yzeeGFFxg6dCh5eXlBh1clJffaGjnS32qr0sg/KwsmTvRXs/7sZ7665Isv+gugRKTu5OfnM2vWLFq2bEmfPn1YvHgxq1atolOnTkGHtl+ac09iZnD55fD2274u/LHHah5epK41adKEoqKi7x43btyYzp07Y2YBRlW9mJK7mfU3s4VmttjMboiyf4SZfWNmxeHbZfEPNXP17eubcHfv7ufhb7kFysuDjkpEklm1yd3MQsAEYADQDRhqZt2iHPq0c65n+PZonOPMeO3awbvv+tLBt97qT7hu2hR0VCLpa9y4cZgZDz74YNT9S5cuJSsri969eydd/1SIbc69N7DYObcUwMyeAgYD8xMZWEaqZtVNNvA4cHsnWPIifHYIFBwJjbIrHThsWHzOA4hksMLCQgDmzZsXdf/111/Pjh07GD9+fFJO0cQyLdMO+DLicWl4W2U/MrO5ZjbFzNrHJbpMMmwY9OxZ7WEG5LaDHj1g+w746CNYvyHigOJimDQpYWGKZIoePXoAUBKlT+aMGTOYMmUK559/Pn379q3r0GISy8g92p+kyt9BXgYmO+e2m9kVwJ+BU/Z5IbORwEggqZcQBeIAV920AA5ZDIMGwaISeOAB+OlPic96e5GqXHPNnmW7yapnT1+Zr5Zyc3Np2bLlPsndOce1115LVlYWd955Z63fJ1FiGbmXApEj8VxgVeQBzrm1zrnt4YePAEdHeyHn3ETnXJFzrqh169Y1iVcidOoEM2fCGWfAlVfCVVdBefJN/YmkrB49erBhw4a9rkSdPHkyH374Iddccw35+fnfbZ8+fTqDBg2iXbt2mBlPPvlk3QccIZaR+yygs5l1BFYCQ4C9rtwxs7bOua/CDwcBC+IapVSpaVOYOhVuuAHuvhv+p7lfVdMg6MAkPcVhRJxKCgsLeeeddygpKSE3N5dt27YxevRo2rRpw+jRo/c6dvPmzRQUFDB8+HCGDx8eUMR7VJvcnXO7zGwU8BoQAh53zpWY2W3AbOfcVODnZjYI2AWsA0YkMGapJBTy5YILCmDjT2DOHGjxOXTpEnRkIqmtYt593rx5/OAHP+C+++5j+fLlPPzwwzRt2nSvYwcOHMjAgQMBGDFiRF2Huo+YrlB1zk0DplXaNibi/o3AjfENTQ7UxRfDxvtg3jxfWfLZZ+G004KOSiR1Ra6YWb16NePGjaOgoIBLL7004MiqpytU00yzZnDU0X5dfP/+voSBiNRM9+7dCYVClJSUMHbsWMrKyrj33nsJpUAlP9WWSUONsuFfr8MFF/jyBZ995qdtUuD3USSpZGdn06VLF+bOncucOXM488wzOf3004MOKyYauaeppk3h5Zd90bHx4+Gcc2Dz5qCjEkk9hYWFbN++HTPj7rvvDjqcmCm5p7H69f369wkT4JVX4MQTYeXKoKMSSS2TJ0/GOcfOnTvp2rVr0OHETMk9A1x5Jfz977Boke/w9PHHQUckkn42b95McXExxcXFlJeXs2LFCoqLi1mxYkUg8Si5Z4gBA2DGDD/vfsIJfiQvIvEze/ZsevXqRa9evdi6dStjx46lV69ejBkzpvonJ4BOqGaQHj3ggw/grLN82YIHHvBXtYpI7fXr1y+pqkMquaej/VSXbAt80Ajmt4C1o2Dx7+Cw//KNQfah6pIiKUvJPd3E0NM1FIKC7rB4CZSWwrZtcMQREIqcpIvo6SoiqUfJPd3EWF3SgM7AtPvhF7+APjt9jZrv6rmpuqRIStMJ1Qx39dXw3HN+oH7ccX5FjYikPiV34ZxzfBPujRvh+ON9GWERSW1K7gLAscfCv/7la9OccgqsWRN0RBKkZFr1kQkS8e+t5C7f6dzZJ/gjj4R5JbByVfXPkfQTCoXYuXNn0GFklJ07d8a9GJmSu+ylTRt46y04uKWff7/pJtAgLrPk5ORQVlYWdBgZpaysjJycnLi+ppK77KNxY9/4o21buOMO+MlPQAO5zNGyZUvWr1/PmjVr2LFjh6ZoEsQ5x44dO1izZg3r16+nZcuWcX19LYWUqMx8J6dbr4CxY+E///HNP5o0CToySbSsrCzy8vJYt24dy5YtY/fu3UGHlLZCoRA5OTnk5eWRlZUV19dWcpcqGTBmDBx6qK8Lf/LJviZNmzZBRyaJlpWVRdu2bWnbtm3QoUgNaVpGqnXZZfDii1BSAn37wtKlQUckItVRcpeYnHWWP9G6bp1fC6+ywSLJTcldYnbssfD++9CwIZx0kr/wSUSSU0zJ3cz6m9lCM1tsZjfs57hzzcyZWVH8QpRkcsQRfi18+/a+AfeUKUFHJCLRVJvczSwETAAGAN2AoWbWLcpxOcDPgQ/iHaQkl9xceO89KCqC88+Hhx8OOiIRqSyWkXtvYLFzbqlzbgfwFDA4ynG3A78HtsUxPklSLVvCG2/AwIFwxRXwm9/oYieRZBLLUsh2wJcRj0uBPpEHmFkvoL1z7u9mdl1VL2RmI4GRAHl5eQcerdSt/TT9ADgImOpg4ffgP7+GxX+CTp38Esq9qOmHSJ2LJblH69Hz3RjNzOoB44ER1b2Qc24iMBGgqKhI47xkFkPTD4B6Bl27QoMGvvHHrp1weFe/HVDTD5GAxJLcS4H2EY9zgciSUjlAAfCO+V5thwBTzWyQc252vAKVOhZj0w/wf/0Pc/Ds72D0aBhY5K9mPegg1PRDJCCxzLnPAjqbWUczawgMAaZW7HTObXTOtXLO5Tvn8oGZgBJ7hjGDG2/0J1dffRXOOAM2bAg6KpHMVW1yd87tAkYBrwELgGeccyVmdpuZDUp0gJJaRo6Ep5+GDz/0g/YdO4KOSCQzxVRbxjk3DZhWaduYKo7tV/uwJJWdd55v+nHOOfBxORQWQnbQQYlkGF2hKglxxhnw5puwc5cvVbBgQdARiWQWJXdJmOOOg549/fr3E0+Ejz4KOiKRzKHkLgnVpDH06uUbgJx8Mrz7btARiWQGJXdJuEaNYMYMX7agf3+YNq3654hI7Si5S51o1w6mT4du3WDwYHjmmaAjEklvSu5SZ1q18jXhjz0Whg6Fxx4LOiKR9KXkLnWqWTN47TW/muayy+C++4KOSCQ9KblLnTvoIN+270c/gl/8Am6/XRUlReJNyV0CkZUFTz0Fw4f7JtzXX68ELxJPMV2hKpII9evDE09AkyZw112waRNMmAD1NOQQqTUldwlUvXrw0EOQkwN33glbtsDjj/vELyI1p/+FJHBmMG6cT/A33wzffguTJvlG3CJSM0ruknjVdHQCXxP+JuDiw2Dxc/BZW+jeHUKRUzTq6CQSMyV3SawYOzpVyM2FeiH4/HP49FM4sgBCIdTRSeQAKblLYh1AR6cKhwJv/R+cNgJ674ZpL0Pzs/slIjqRtKXkLknpwgv9evghQ+DUU2Fmtu/TKiKx0aIzSVo//CG89BLMn+9nZdTVSSR2Su6S1AYMgFdegW3b4ONi+PLLoCMSSQ1K7pL0TjkFehT6kfuJJ8IXXwQdkUjyU3KXlNCsKfQshLIyOOEEv5pGRKqm5C4pIycH3n4bdu70I/h584KOSCR5xZTczay/mS00s8VmdkOU/VeY2admVmxm75tZt/iHKgI9evhWfaGQvy5qzpygIxJJTtUmdzMLAROAAUA3YGiU5D3JOXekc64n8Hvg3rhHKhLWtavv6tS4sZ+Pnzkz6IhEkk8sI/fewGLn3FLn3A7gKWBw5AHOubKIh40BFW+VhDrsMHjvPd/d6fTTfbIXkT1iSe7tgMgFaKXhbXsxs6vMbAl+5P7zaC9kZiPNbLaZzf7mm29qEq/Id/LyfFJv39433n7jjaAjEkkesVyhalG27TMyd85NACaY2TDgZuDiKMdMBCYCFBUVaXQvByZKAbJDgeIW8MkK2PoDWNsdDj64mtdRATLJALGM3EuB9hGPc4FV+zn+KeDs2gQlso9hw6Bnz6i7Gjbwuxo39ito9vulsLjY1xMWSXOxjNxnAZ3NrCOwEhgC7FXqz8w6O+cWhR+eCSxCJJ6qKUDWAOi0EQYOhA8+gD+Phh//OMqB1ZQeFkkX1SZ359wuMxsFvAaEgMedcyVmdhsw2zk3FRhlZqcBO4H1RJmSEUm0Zs3gtddg0CC46CJfsuDSS4OOSiQYMVWFdM5NA6ZV2jYm4v7VcY5LpEaaNPG1aH74Q7jsMp/gr7oq6KhE6p6uUJW006gRvPgiDB4Mo0bB3XcHHZFI3VNyl7SUlQXPPgsXXAC/+hXcdhs4rc+SDKJmHZK2GjSAv/3Nj+THjvWNt8cRfW2vSLpRcpe0FgrBY49BdjbceSdc2g46dVKCl/Sn5C5pr149+MMf/Ah+5Xgo3w2ddocbb4ukKc25S0Ywg3vugQ558NXXMHw47NoVdFQiiaORu2QMM+jYEeqF/EWqW7fC5Mn+5KtIutHIXTJOhzy47z544QU4+2yf5EXSjZK7ZKSrr4ZHHvFXtA4YAJs2BR2RSHwpuUvGuuwyv1Ty/fd9Tfj164OOSCR+lNwlow0dCs89Bx9/DCefDKtXBx2RSHwouUvGGzwY/v53+Pxz33j7yy+rf45IstNqGck8UZp+nA58dTjM/RRWdoZWhX5dfJXU8EOSnEbukln20/SjWTO/a/duP02zZUsVr6GGH5ICNHKXzFJN048coMV8f4J121J49VXo3bvSQWr4ISlAI3eRSrp1g/feg+bN4dRT4e23g45I5MApuYtE8V//5RN8Xp5fB//yy0FHJHJglNxFqnDooTB9Ohx5JJxzjl8TL5IqNOcush8HHwxvveWXS154ob/QaVTQQYnEQCN3kWrk5MC0aT7B/+xnsGw5qKmTJDsld5EYZGfDlCkwYgQsWwaLF0F5edBRiVQtpuRuZv3NbKGZLTazG6Lsv9bM5pvZXDP7p5l1iH+oIsGqX993dWqfCytX+WmaHTuCjkokumqTu5mFgAnAAKAbMNTMulU67GOgyDnXA5gC/D7egYokg3r14LDD/GqayZNh0CDYvDnoqET2FcvIvTew2Dm31Dm3A3gKGBx5gHPubefct+GHM4Hc+IYpklzy2vtR/Btv+LXwa9YEHZHI3mJJ7u2AyFJKpeFtVbkUeDXaDjMbaWazzWz2N998E3uUIknokkt8w4+5c6FvX1i+POiIRPaIJblHaxQfdbGAmV0IFAF3RdvvnJvonCtyzhW1bt069ihFktSgQX70vno1HHccfPJJ0BGJeLEk91KgfcTjXGBV5YPM7DTgJmCQc257fMITSX7f/76/mrVePV8y+K23go5IJLbkPgvobGYdzawhMASYGnmAmfUCHsYndrU7kIxTUAD//je0bw/9+8NTTwUdkWS6apO7c24X/qK814AFwDPOuRIzu83MBoUPuwtoAjxrZsVmNrWKlxNJW+3b+xH8ccf5Dk/33ANOVztJQGIqP+CcmwZMq7RtTMT90+Icl0hKatHCN90ePhyuu86fZB0/HkKhoCOTTKPaMiI1EaWbU4Vs4GngjlwofRAWTIIjukEo2vdkdXSSBFFyFzlQw4ZVe4gBnQ7zZQsWL/Z/C44sgIYNIw4qLvY/ldwlAZTcRQ5UNd2cIuUCs1+EgcPge1/BK6/4ZiCAOjpJQqlwmEiCnX22rwu/dSscfzz8859BRySZQMldpA4UFcEHH/gVNT/4AUycGHREku6U3EXqSIcOMGMGnHEGXH45LF6ipZKSOEruInWoaVOYOtU3/SgthXnzoKws6KgkHSm5i9Sx+vXhgQegc2dYt85f9LRkSdBRSbpRchcJSLtDoUchfPUV9O6tmjQSX1oKKRKgFsuKWXl4Pz6dB1tPhdLDoF1u9FKs+6WLoaQSJXeRoIQvhmoEHHUULFjgT7Ju3gxduvgqkzHRxVAShZK7SFAiLoaqD3Qvh9tug1tvhaL28NxzkJcXw+voYiiJQnPuIkmiXj245RZ48UX4/HM4+mhd8CQ1p+QukmQGD4ZZs6BNG78m/ne/g/LyoKOSVKPkLpKEunTxV7Sedx7ceKNP+OvXBx2VpBIld5Ek1aQJTJ4MDz7oa8QfdRTMnh10VJIqlNxFkpgZjBrlOzzt3u0Lj91/v8oWSPWU3EVSQJ8+8PHHvj/rNdf4SpPr1gUdlSQzJXeRFHHwwfDSS75t36uvQmEhvPtu0FFJslJyF0khZn7k/q9/+S5PJ58MX3wB5ZqmkUpiSu5m1t/MFprZYjO7Icr+E81sjpntMrNz4x+miEQqKvLTNCNGwPIVUPyxXxsvUqHa5G5mIWACMADoBgw1s26VDlsBjAAmxTtAEYmuSRN4/HHftu/brdCrF/zxjzrZKl4sI/fewGLn3FLn3A7gKWBw5AHOuWXOubmALrUQqWNtWsMxRfD978OVV8KAAb5WvGS2WGrLtAO+jHhcCvRJTDgiUhNZC4r5R89+rOoMS96AZfkQ6gSHHHKAFSZVXTJtxJLco/1u1OiLn5mNBEYC5MVUEUlEqhWuLmn4GvEtW8BnC2HhQvhmtb/aNTs7htdRdcm0EktyLwXaRzzOBVbV5M2ccxOBiQBFRUWaGRSJh4jqkuBLCBeWw4QJvnQB82DcOD9lEwrt53VUXTKtxDLnPgvobGYdzawhMASYmtiwRKQ26tXzfVrnzfNz8T//uf85d27QkUldqTa5O+d2AaOA14AFwDPOuRIzu83MBgGY2TFmVgqcBzxsZiWJDFpEYpOf7y94+utfYfFiX5/muut8QxBJbzGtc3fOTXPOdXHOHeac+2142xjn3NTw/VnOuVznXGPn3MHOue6JDFpEYmcGF17o5+AvuQTuuQeOOAKeflrLJtOZrlAVyRAtW8LEif7q1latYMgQOOUU+PTToCOTRFByF8kwxx3nSwf/6U9+Dr5nT/jpT2HHzqAjk3hSchfJQKEQXH45LFoEV10Fjz7qm4OsWAFbtwYdncSDkrtIBmvZEh54wK+qad4cln7h18U/9hjs2hV0dFIbSu4iwuGHw5EF0LMQ2rWDyy6DI4+EZ55R/9ZUpeQuIt9p3hz+/W94/nm/yuaCC3zd+OeeU5JPNUruIrIXMzjnHL+KZtIk2LEDzj3XJ/lJkzRdkyqU3EUkqlAIhg6FkhL4v//zI/cf/xi6doWHH9aJ12QXS20ZEckUxcX71JipD/wYGNYa1nSHFcth0xUwZ5Sfnz+0HTRsUOl1VF0ycEruIuKFq0tWxYDWrfwFUBs3wJdfwrLlvhNUmzaQ2w5yclB1ySRhLqDrj4uKitzs2bMDeW8RiY/PPoOHHoInn4QtW6B3b3hxYz8O+boY69mz9m+gbwD7MLOPnHNF1R2nOXcRqbGuXX1yX7kS7r8fNm2CsQuH8f7mnnz+OZRtqmHzB/DfACapc2dNaeQuInHjHMyY4U+4TpkC27b5ImXDh/taNvn5B/BiFXP/77wT/0BTmEbuIlLnzHzd+L/+Fb7+Gh55BFq08E1DOnaEvn39SH9Vjdr9yIFQcheRhGjWzF/pOmMGLF0Kd9wBZWW+iUi7dj7R33OPr28j8afkLiIJ17GjH71/+inMnw+33+5PwF53na9l07Ur/OpX8NZbsH170NGmByV3EalTRxwBN9/sz5d+8YUvXJaX50/InnoqHHwwnHUWlJbC5i0qe1BTWucuIoHJz/fTND/7mW/99/bb8I9/wOuvw+Il/phz28AJJ+y5FRZCw4aBhp0SlNxFJCk0aeJH7Ged5R9vOw42rIf/Phbeew9efNFvz86Go4+GPn2gqMjf79TJNwWXPZTcRSQpZWfBIYf4C6TAr7CZMQNmzvS3CRP2zM83bQqWx61sAAAMEUlEQVQ9evhRfY8eUFAA3br5KpeZSsldRJJXRK2bQ4HzwjcaQHkf+HaLv3Bq02bYMhc2/xt274YdQDF++uagg6CkcBjfnD2STp3gsMP8Cd7s7KA+VN2IKbmbWX/gfiAEPOqc+12l/VnAX4CjgbXABc65ZfENVUQySjW1buqZn8pp0gTahrc5/IVT326BLd/6FTm5a4rZ/D5c+O7eZQwOPRQ6dPC3vDzIzYX27f32tm39t4YGlQuipZBqr1A1sxDwOXA6UArMAoY65+ZHHHMl0MM5d4WZDQHOcc5dsL/X1RWqIlIn+vXDAWuefYclS/jutmwZLF/uf5aW+rr1lbVq5Yuife97/mercOG0gw/2LQpbtPA/mzf3t2bN/DcCs8R9nFivUI1l5N4bWOycWxp+4aeAwcD8iGMGA7eE708BHjIzc0HVNhARiWBA69b+duyx++4vL4c1a3yly1Wr4Kuv/M+vv4bVq/1tzhxYuxbWrdv/e9Wv788B5OT4W8W3i8aN/RTRQQfBRRfBSScl5KPuiSOGY9oBX0Y8LgX6VHWMc26XmW0EDgbWxCNIEZFaiVKnPlI9oE34dnS0A+rjJ/0P9fVzdu7yHal27Yy4H77t3g27d8GubbB7S/hx+Fa+G3aXQ/Mve8JJ9yXgg+4dcnWifcGoPCKP5RjMbCQwEiAvLy+GtxYRqaVq5u4PlJlvTtKwAdCohi9yRDwjii6W5F4KtI94nAtULvtTcUypmdUHmgH7fHlxzk0EJoKfc69JwCIiB2TkyIysCR/Lsv9ZQGcz62hmDYEhwNRKx0wFLg7fPxd4S/PtIiLBqXbkHp5DHwW8hl8K+bhzrsTMbgNmO+emAo8BfzWzxfgR+5BEBi0iIvsX0zp359w0YFqlbWMi7m8jfG2BiIgET9UYRETSkJK7iEgaUnIXEUlDSu4iImlIyV1EJA1VWzgsYW9s9g2wvIZPb0X6lDbQZ0k+6fI5QJ8lWdXms3RwzrWu7qDAknttmNnsWKqipQJ9luSTLp8D9FmSVV18Fk3LiIikISV3EZE0lKrJfWLQAcSRPkvySZfPAfosySrhnyUl59xFRGT/UnXkLiIi+5Gyyd3MbjezuWZWbGavm9mhQcdUU2Z2l5l9Fv48L5hZ86BjqgkzO8/MSsys3MxSclWDmfU3s4VmttjMbgg6npoys8fNbLWZzQs6ltows/Zm9raZLQj/bl0ddEw1ZWbZZvahmX0S/iy3JvT9UnVaxsyaOufKwvd/DnRzzl0RcFg1YmZn4Gvg7zKzOwGcc9cHHNYBM7MjgHLgYeA651xKdUCPpRl8qjCzE4HNwF+ccwVBx1NTZtYWaOucm2NmOcBHwNkp+t/EgMbOuc1m1gB4H7jaOTczEe+XsiP3isQe1pgobf1ShXPudefcrvDDmfhuVynHObfAObcw6Dhq4btm8M65HUBFM/iU45ybTpRuaKnGOfeVc25O+P4mYAG+Z3PKcd7m8MMG4VvC8lbKJncAM/utmX0J/BgYU93xKeIS4NWgg8hQ0ZrBp2QiSUdmlg/0Aj4INpKaM7OQmRUDq4E3nHMJ+yxJndzN7E0zmxflNhjAOXeTc6498DdgVLDR7l91nyV8zE3ALvznSUqxfI4UFlOjd6l7ZtYEeA64ptK39pTinNvtnOuJ/3be28wSNmUWUyemoDjnTovx0EnAK8DYBIZTK9V9FjO7GPhv4NRk7j97AP9NUlEszeCljoXnp58D/uacez7oeOLBObfBzN4B+gMJOemd1CP3/TGzzhEPBwGfBRVLbZlZf+B6YJBz7tug48lgsTSDlzoUPgn5GLDAOXdv0PHUhpm1rlgJZ2aNgNNIYN5K5dUyzwGH41dnLAeucM6tDDaqmgk3Fs8C1oY3zUzFlT9mdg7wINAa2AAUO+d+EGxUB8bMBgL3sacZ/G8DDqlGzGwy0A9fffA/wFjn3GOBBlUDZvZ94D3gU/z/6wCjw32dU4qZ9QD+jP/dqgc845y7LWHvl6rJXUREqpay0zIiIlI1JXcRkTSk5C4ikoaU3EVE0pCSu4hIGlJyFxFJQ0ruIiJpSMldkpKZjTAzFy4WlfLMrKeZvWtmZeHPdXbQMUl6U3IXSbBwnfgp+Lo11wMXAVFr3YcbnrjwlbKV970Q3vfzKPseCDdJOSzO4UuKSurCYSJpoiNwGPAL59wfqzl2Q/hn08iNZtYeOCv8sFmlfY2Bi4FpzrkltQ9X0oGSu0jitQn/3LDfo7yN4Z9NK22/EtiKrz9Ued9F4W0P1TRAST+alpFaM7Nzw9MFZ0bZd0J43yXhxx3M7KFwT8wt4TnoN83s+Bje5xYz26cYkpn1C79Hv4ht3zOzP5nZSjPbEe6JeqOZ1av03MZmdqeZLTGzbWa21sxmmtm5McTT3cxeMrMNZvZt+Hn/XemYJ4EZ4YdPhONctp+XrUjuORGvkQVcii86tYroiX8R8Fp1MUvm0Mhd4uHvQBkwFF9XP9IQYDtQUYf7GODk8OPl+KqFlwJvmVmRc67Wta3NrBW+XWE2MBGfEPsCdwAdgMiKm38Ix/0HfF3tpkBPoA9+nryq9+gC/AvYia8iWQaMAKaa2QXOuWfDhz4MrAB+HY7lPXxv06pEG7kPwf87PQg8ELkv3Cv1SHwvTlUBlD2cc7rpVusb8CSwCWgUsS2Ebyf2QsS2g6I8t2X4uEcito3Ad0HKj9h2i/+V3ef5/cLH9gs/fhhYg2+sHHncHfiysV0itq0HJtTg807Bd83qHrEtB1gKrATqR2z/fji+ETG8bqPwseMjtn0IvBbxvq9E7Hs6/O/eNOjfAd2S66ZpGYmXyUAT9pz0A9+MoHV4HwAuohmJmTUys4Px04MfAkfXNohwc4fzgGnATjNrVXHDT1sY/ptDhQ1An/AJy1jfI4TvoDPNOVdSsd35Bs5/BA4FjqpJ/M65rfhvA03D79UH/23ngfAhZRH72gLnAH9xKdx6ThJDyV3i5U386HtIxLYh+FHlyxUbzKyh+cbmK4Bv8SPsb4AzgeZxiKM10AJ/kvGbSrd3wse0iTj+l8ARwHIzKzazu8ysuj8yrYHGRO+iMz/8s2ONovc2smfOfRSwhD1N0zexZ1pmJNCASidSzexn4c+yy8xuqUUcksI05y5x4ZzbbWbPApeZWTNgG35U+WJ4NFrhfnxSmoA/0bgeP1VyI3654H7fportoYj7FQOWp4FHqzh+aUTcz5vZ+/hvHKcBlwC/NLObnHPjqoknmoom27WZ/94INDWzNvhvITc45yq6EG0K76uP/3f8p3NuQaXnlwI34z+LZCgld4mnScBV+KS+Eb8ee3KlY4bgpxH2uhDHzGJpN7Y+fGxz51zkssL8iPvf4KcuGjrn3owlaOfcanyfzsfM7CD8SeFbzexu59zOKE/5BtgCdI2yr2Lbsljeuwob8aPzkfgpmici9lWM3M/BT/9cWfnJzrkXwK9iqkUMkuI0LSNx45z7Fz6pDQnf1gBvVDqsnEq/d2Z2AnBsDG+xOPzzuznz8Aj2u9UvzrndwLPAIDM7pvILmFlOeGkhZhYKf8uI/AzfAgvx0x2NowURfo9XgQFmdkTEazcBfopfnTMnhs9TlY34k8yX4/8QbozYVzHnfhV+tdHL+z5dRCN3ib+ngOvwI84nnXO7Ku1/CbjYzDYDxfj57suAEiLWdlfhdfwfj0fNrCv+op5hUY67Eb+C5j0zewyYG37t7sC5+KWDy8LbVprZC8AnwDqgVzieVyt9O6jsZuAMYLqZPcSepZAdgQuifO4DsZE9f8AqX5i0Cf//7UnA9RHTNSJ7UXKXeJsE3ID/3ZoUZf/V+Pn4HwI/wXe1PxefpPvt74Wdc7vCBbcm4JdFrgUeAabjT+hWHPdNeJXJzcBg4H/wq2IWAbcDX4cP/RafPE/Dn9DNwq9JvwP4fTWxLAxfeDUOuBZoiP8DMcg59/f9PTcGFSP1N6LMp28K/9yGn0oSicqc03UPIukofHXsMufcLQGHIgHQyF0kzYTPQ9THryKqb2bZwK5aThVJitEJVZH0czP+fMSFwE3h+zcHGpHUOU3LiIikIY3cRUTSkJK7iEgaUnIXEUlDSu4iImlIyV1EJA0puYuIpCEldxGRNKTkLiKShv4fnx/LI70hfCoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "##################################################\n",
    "### 2. Bermudan put price by empirical regression \n",
    "##################################################\n",
    "\n",
    "#################################################\n",
    "# First of all: \n",
    "# + Generate a sample of M values of W_2\n",
    "#   and the corresponding values of S_1\n",
    "# + Taking into account that the values of W_1\n",
    "#   have already been simulated above\n",
    "#################################################\n",
    "W_2 = W_1 + np.sqrt(t_2 - t_1) * np.random.randn(M) \n",
    "S_2 = GBM(S_0, t_2, sigma, W_2) \n",
    "################################################\n",
    "\n",
    "###############################\n",
    "## We generate the cells I_k\n",
    "###############################\n",
    "a = 3. * np.sqrt(t_1)\n",
    " \n",
    "#################################################\n",
    "# n = the dimension of the approximation space\n",
    "# \n",
    "# It can be shown that n = M*{1/3} is a pertinent\n",
    "# choice for a problem in dimension d = 1, as here\n",
    "#################################################\n",
    "n = int(M**(1./3))\n",
    " \n",
    "partition_points = np.linspace(-a, a, n+1)\n",
    " \n",
    "####################################################\n",
    "## TO DO: evaluate\n",
    "## + The regression coefficients alpha\n",
    "## + The vector v_1_tilde(W1) of size M \n",
    "##   inside the array v_1_tilde\n",
    "####################################################\n",
    "def EmpiricalRegression(W_1, W_2, partition_points, t_2, K, S_0, sigma):\n",
    "    step = partition_points[1] - partition_points[0]\n",
    "    \n",
    "    alpha = np.zeros(partition_points.size - 1)\n",
    "    \n",
    "    v_1_tilde = np.zeros(M)\n",
    "     \n",
    "    for k in range(n):\n",
    "        \n",
    "        leftPoint = partition_points[k]\n",
    "        \n",
    "        insideCell = np.logical_and( leftPoint <= W_1, W_1 < leftPoint+step )         \n",
    "        \n",
    "        GBM_t_2 = GBM(S_0, t_2, sigma, W_2[insideCell])\n",
    "        \n",
    "        #############################################\n",
    "        ## TO DO: evaluate the coefficients alpha[k]\n",
    "        #############################################\n",
    "        numberInsideCell = np.sum(insideCell)\n",
    "        \n",
    "        if numberInsideCell != 0:\n",
    "            alpha[k] = np.sum( pos_part(K - GBM_t_2) ) / numberInsideCell\n",
    "        \n",
    "        #############################################\n",
    "        ## TO DO: evaluate the empirical approximation\n",
    "        ## v_1_tilde(W1) inside the array v_1_tilde\n",
    "        #############################################\n",
    "        v_1_tilde[insideCell] = alpha[k] \n",
    "     \n",
    "    return alpha, v_1_tilde\n",
    " \n",
    "time3 = time()\n",
    " \n",
    "alpha, v_1_tilde = EmpiricalRegression(W_1, W_2, partition_points, t_2, K, S_0, sigma) \n",
    "\n",
    "################################################\n",
    "# TO DO: evaluate the new estimator Y_tilde_0^M \n",
    "# of the Bermudan put price\n",
    "\n",
    "price_EmpRegr = np.mean(  np.maximum( pos_part(K -  S_1), v_1_tilde) ) \n",
    "\n",
    "price_EmpRegr = np.maximum(price_EmpRegr, 0)\n",
    "\n",
    "################################################\n",
    " \n",
    "time4 = time()\n",
    " \n",
    "print(\"Price by empirical regression = %1.4f\" %price_EmpRegr)\n",
    "print(\"Time: %1.4f \\n\" %(time4 - time3 + timeSimulations))\n",
    " \n",
    "#######################################################\n",
    "## We can also plot the true function v_1 and its\n",
    "## empirical approximation v_1_tilde for comparison\n",
    "#######################################################\n",
    "x = np.linspace(-a, a, 100)\n",
    " \n",
    "v_1_exact = putBlackScholes(S_0, t_1, t_2, K, x, sigma)\n",
    " \n",
    "plt.plot(x, v_1_exact, color=\"b\", label=\"$v_1$\")\n",
    " \n",
    "plt.step(partition_points, np.append(alpha, alpha[-1]), where=\"post\", color=\"r\", label=r\"$\\tilde{v}_1$\")\n",
    " \n",
    "plt.xlabel(\"values of $W_1$\", fontsize=17)\n",
    "plt.legend(loc=\"best\", fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Nested Monte-Carlo price\n",
    "\n",
    "In the spirit of Exercise 1, we can also propose a Nested Monte-Carlo estimator that does not require the knowledge of the inner conditional expectation $\\mathbb E \\left[ (K - S_2)^+ \\big|W_1 = x \\right]$.\n",
    "\n",
    "This third estimator is\n",
    "\n",
    "$$\n",
    "\\hat{Y}_0^{M} = \\max \\Bigl\\{ (K-x_0)^+,\n",
    "\\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, \\hat{v}_{1}^m \\right)\n",
    "\\Bigr\\}\n",
    "$$\n",
    "\n",
    "where now, for every $m$, the estimate $\\hat{v}_{1}^m$ of $v_1(W_m^1)$ is obtained from a sample of $N$ i.i.d. values $(\\Delta W^{m, n})_{1 \\le n \\le N}$ of the Brownian increment between $t_1$ and $t_2$:\n",
    "\n",
    "$$\n",
    "\\hat{v}_{1}^m = \\frac1{N} \\sum_{n = 1}^N \\Bigl(K - S_0 \\exp\n",
    "\\Bigl(-\\frac12 \\sigma^2 T_2 \\ + \\sigma (W_1^m + \\Delta W^{m, n} \\Bigr)\n",
    "\\Bigr)^+.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\blacktriangleright$ Simulate the Nested estimator $\\hat{Y}_0^{M}$ in the cell below, and compare with the previous methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nested MC price = 0.2479\n",
      "Time: 0.1940 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "##################################\n",
    "## 3. Nested Monte-Carlo price\n",
    "##################################\n",
    " \n",
    "###################################################\n",
    "# We now have to generate N samples of S_2 for every\n",
    "# point in the sample of M values of W_1.\n",
    "#\n",
    "# We will therefore need to repeat N gaussian iid samples\n",
    "# of the Browian increment W_2 - W_1\n",
    "# for every value of W_1 previously simulated\n",
    "###################################################\n",
    "N = int(np.sqrt(M))\n",
    "\n",
    "# This variable will contain the sum over m in {1,...,M}\n",
    "sum_1 = 0.\n",
    " \n",
    "time7 = time()\n",
    " \n",
    "for m in range(M):\n",
    "    w1 = W_1[m]\n",
    "    s1 = GBM(S_0=S_0, t=t_1, sigma=sigma, W=w1)\n",
    "     \n",
    "    ###################################################\n",
    "    ## To Do: \n",
    "    ## - Simulation of N Brownian increments\n",
    "    ## - Generation of W_2, therefore of S_2, conditionally to W1 = w1\n",
    "    ## - update the sum over m with the current contribution\n",
    "    ###################################################\n",
    "    G = np.random.randn(N)\n",
    "    \n",
    "    W2 = w1 + np.sqrt(t_2 - t_1) * G\n",
    "    \n",
    "    S2 = GBM(S_0=S_0, t=t_2, sigma=sigma, W=W2)\n",
    "    \n",
    "    v_1_hat = np.mean( pos_part(K - S2) )\n",
    "    \n",
    "    sum_1 += np.maximum( v_1_hat, pos_part(K - s1) )\n",
    "    ###################################################\n",
    "\n",
    "Y_0_hat = np.maximum( np.maximum(K-S_0, 0.), sum_1 / M )\n",
    "    \n",
    "time8 = time()\n",
    " \n",
    "print(\"Nested MC price = %1.4f\" %Y_0_hat)\n",
    "print(\"Time: %1.4f\" %(time8 - time7), \"\\n\")"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}