Computational_Finance_A2_V1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0paxbpsJQBG8"
   },
   "source": [
    "# Assignment 2\n",
    "\n",
    "**Deadline**:  December 8, 23.59.\n",
    "\n",
    "Your notebook should not give errors when executed with `Run All`. Please submit your answers via [Canvas](https://canvas.uva.nl/courses/12430/assignments/78394).\n",
    "\n",
    "|**Name**|**Student ID**|**Email**|\n",
    "|:-------|:-------------|:--------|\n",
    "|Jason Spano|*******|******|\n",
    "\n",
    "\n",
    "****Hand in the following:****\n",
    "* Your notebook. N.B. **click on `Kernel`-> `Restart & Run All`** before submitting, see notes.\n",
    "* A (printed) pdf version of your notebook.\n",
    "\n",
    "****NOTES****:\n",
    "* This is an individual assignment. Teamwork is **not** permitted, including during lectures.\n",
    "* The assignment is a partial stand-in for a final examination, so the usual rules regarding plagiarism and fraud apply, with all attendant consequences. Code found on the internet or elsewhere is not acceptable as a solution.\n",
    "* Efficient code is required for full marks.\n",
    "* Make sure that any function you write has a docstring, and comments where appropriate.\n",
    "* Before submitting your work, **click on `Kernel`-> `Restart & Run All`** and verify that your notebook produces the desired results and does not error.\n",
    "* Some questions require you to write code to obtain a numerical result (e.g., an option price). In that case, don't just give the function, but also the result of calling it with the given parameter values (i.e., the numerial value that it returns). If your function uses random numbers, then set the seed to 0 before calling it. This makes it much easier to grade the assignments (at least as long as the answer is correct).\n",
    "\n",
    "\n",
    "**Declaration of Originality**:\n",
    "\n",
    "By submitting these answers, I declare that\n",
    "1. I have read and understood the notes above.\n",
    "2. These solutions are solely my own work.\n",
    "3. I have not made these solutions available to any other student."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EGURrplgQBG_"
   },
   "source": [
    "**Question A**:\n",
    "Consider an $N$-period model with two assets: a riskless bond with value $B_t=e^{rt}$, and a stock that evolves according to a binomial tree, such that $S_t=S_{t-\\delta t}u$ or $S_t=S_{t-\\delta t}d$ in the good and bad states of the world, respectively. Here $u=1/d=\\exp(\\sigma\\sqrt{\\delta t})$ and $\\delta t= T/N$ as usual. We would like to price a *chooser* option, which has the property that, after $M<N$ periods, the holder can choose whether the option is a European call or put option, both maturing after $N$ periods, and with the same strike price $K$ (See Hull, section 25.7 (ed. 8) or 26.8 (ed. 9/10)). Denoting their prices after $i$ periods as $C_{i}$ and $P_{i}$, respectively, the payoff of the chooser option after $M$ periods is $V_{M}=\\max \\{C_{M},P_{M}\\}\\,$. Our starting point is the function `calltree` (or `calltree_numpy`) from Week 4. **For each of the following subquestions, your answer should include the function itself, as well as the result of calling it with the following parameter values**: $S_0 = 9$, $K = 11$, $T = 3/12$, $r = 0.025$, $\\sigma = 0.3$, $q = 0.01$, $N = 500$, $M=200$.\n",
    "\n",
    "   1. We want to allow for a dividend yield, $q$. Construct a calltree pricer (vectorised) which accepts an additional input argument `q`. That is, the function becomes `calltree_numpy(S0, K, T, r, sigma, q, N)`. In the function change the risk-neutral probability to $(e^{(r-q)\\delta t}-d)/(u-d)$.\n",
    "\n",
    "   * Building on your function from the previous question write a function `chooser(S0, K, T, r, sigma, q, N, M)` that prices a European chooser option. Along with the call price tree, this function should also compute a put price tree stored in a matrix `P` (of the same size as the call price matrix). Following computation of `C` and `P`, create a $(M+1) \\times (M+1)$ matrix `CH`, which stores the tree for the chooser option. At point in time `M`, the last column of this tree will be populated with $V_M=\\max \\{C_{M},P_{M}\\}\\,$. Based on the `CH` tree perform the backwards induction to determine the price of the chooser.\n",
    "\n",
    "   * Now consider a chooser option for which after $M$ periods, the holder will choose between an *American* call or an *American* put (above we were choosing between a *European* call and a *European* put). Modify the function from above so that before populating the chooser tree `CH` the first loop also determines whether to exercise or keep the respective *American* put and *American* call. Call the function that evaluates the price of this particular derivative `chooser_american(S0, K, T, r, sigma, q, N, M)`.\n",
    "\n",
    "   * For the given parameters values ($S_0$, $K$, $T$, $r$, $\\sigma$, $q$, $N$) plot the prices of `chooser` and `chooser_american` as a function of $M$. That is, produce a graph that has the prices of the two options on the y-axis and $M$ on the x-axis. Describe and explain (in words) the price convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Mm2bVDN8QBHC"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats.mstats as stats\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "T1_ELdy9QBHN"
   },
   "source": [
    "Question A1. We want to allow for a dividend yield, $q$. Construct a calltree pricer (vectorised) which accepts an additional input argument `q`. That is, the function becomes `calltree_numpy(S0, K, T, r, sigma, q, N)`. In the function change the risk-neutral probability to $(e^{(r-q)\\delta t}-d)/(u-d)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "s9pFSdwyQBHP"
   },
   "outputs": [],
   "source": [
    "#def calltree with addition of dividend term (q)\n",
    "def calltree_numpy(S0, K, T, r, sigma, q, N): \n",
    "    \n",
    "    ## define delta(T)\n",
    "    deltaT = T / float(N)\n",
    "    \n",
    "    ## define the up and down movementsi n the tree \n",
    "    u = np.exp(sigma*np.sqrt(deltaT))\n",
    "    d = 1.0 / u\n",
    "    \n",
    "    #risk neutral probability adjusted according to specification (of the question)\n",
    "    p = (np.exp((r-q)*deltaT)-d) / (u-d)\n",
    "    \n",
    "    ## code the up and down movements in the tree \n",
    "    piu = np.exp(-r*deltaT) * p\n",
    "    pid = np.exp(-r*deltaT) * (1-p)\n",
    "    \n",
    "    ## create a stock price tree based off the movement of the stock price under (u) and (d)\n",
    "    C = np.zeros((N+1, N+1))\n",
    "    S = S0 * u**np.arange(N+1) * d**(2*np.arange(N+1)[:, np.newaxis])\n",
    "    \n",
    "    #keep only the upper triangular part of the matrix (where the payoff is != 0)\n",
    "    S = np.triu(S) \n",
    "    \n",
    "    ## apply the options payoff function\n",
    "    C[:, N] = np.maximum(0, S[:, N]-K) \n",
    "    \n",
    "    ## backwards induction after the final node: \n",
    "    for j in range(N-1, -1, -1):\n",
    "        C[:j+1, j] = piu * C[:j+1, j+1] + pid * C[1:j+2, j+1]\n",
    "        \n",
    "    return C[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "jM9xkR-SQBHU",
    "outputId": "f25ad19f-60cd-44a7-fa7c-28c930006e5d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06554359536314218"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S0=9.; K=11.; T=3/12.; r=.025; sigma=.3; q=.01; N=500; #run function with given values\n",
    "calltree_numpy(S0, K, T, r, sigma, q, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "11q4Vet7QBHd",
    "outputId": "c816e5ec-9a5e-491b-d524-9b3445c9dce8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.48 ms ± 564 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit calltree_numpy(S0, K, T, r, sigma, q, N) #ipython magic for timing things"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EMkaN-pEQBHj"
   },
   "source": [
    "Question A2. Building on your function from the previous question write a function `chooser(S0, K, T, r, sigma, q, N, M)` that prices a European chooser option. Along with the call price tree, this function should also compute a put price tree stored in a matrix `P` (of the same size as the call price matrix). Following computation of `C` and `P`, create a $(M+1) \\times (M+1)$ matrix `CH`, which stores the tree for the chooser option. At point in time `M`, the last column of this tree will be populated with $V_M=\\max \\{C_{M},P_{M}\\}\\,$. Based on the `CH` tree perform the backwards induction to determine the price of the chooser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-4GDd_vZQBHl"
   },
   "outputs": [],
   "source": [
    "## chooser function with additional dividend term (q)\n",
    "def chooser(S0, K, T, r, sigma, q, N, M): \n",
    "    \n",
    "    ## define time intervals delta_t\n",
    "    deltaT = T / float(N)\n",
    "    \n",
    "    ## define up and down movements in the tree\n",
    "    u = np.exp(sigma*np.sqrt(deltaT))\n",
    "    d = 1.0 / u\n",
    "    \n",
    "    #risk neutral probability adjusted according to specification\n",
    "    p = (np.exp((r-q)*deltaT)-d) / (u-d) \n",
    "    piu = np.exp(-r*deltaT) * p\n",
    "    pid = np.exp(-r*deltaT) * (1-p)\n",
    "    \n",
    "    # create matrix P/C\n",
    "    C = np.zeros((N+1, N+1))\n",
    "    P = np.zeros((N+1, N+1)) \n",
    "    \n",
    "    S = S0 * u**np.arange(N+1) * d**(2*np.arange(N+1)[:, np.newaxis])\n",
    "    \n",
    "    ## keep only the upper triangular part (not necessary but speed enhancement)\n",
    "    S = np.triu(S) \n",
    "    \n",
    "    ## calculate payoffs \n",
    "    C[:, N] = np.maximum(0, S[:, N]-K) \n",
    "    P[:, N] = np.maximum(0, K-S[:, N]) \n",
    "    \n",
    "    #for loop to populate put and call matrices\n",
    "    for j in range(N-1, -1, -1):\n",
    "        C[:j+1, j] = piu * C[:j+1, j+1] + pid * C[1:j+2, j+1]\n",
    "        P[:j+1, j] = piu * P[:j+1, j+1] + pid * P[1:j+2, j+1]\n",
    "        \n",
    "    #return C[0, 0], P[0, 0] #check to see prices for put and call\n",
    "    \n",
    "    #create chooser matrix\n",
    "    CH = np.zeros((M+1, M+1)) \n",
    "    \n",
    "    #populate last column with maximum of put/call\n",
    "    CH[:, M] = np.maximum(C[:M+1, M], P[:M+1, M])\n",
    "    \n",
    "    #backwards induction\n",
    "    for i in range(M-1, -1, -1):\n",
    "        CH[:i+1, i] = piu * CH[:i+1, i+1] + pid * CH[1:i+2, i+1]\n",
    "        \n",
    "    ## return price     \n",
    "    return CH[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "bisilVExQBHu",
    "outputId": "3c0e00b0-4597-4deb-d35a-604c04b2c08f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.025914426146793"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#run function with given values\n",
    "S0=9.; K=11.; T=3/12.; r=.025; sigma=.3; q=.01; N=500; M=200; \n",
    "chooser(S0, K, T, r, sigma, q, N, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "jgpFwep-QBH0",
    "outputId": "9a7e9917-1927-4bc3-f76c-83c1c78fedd4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13.8 ms ± 1.14 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit chooser(S0, K, T, r, sigma, q, N, M) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2WOLP67yQBH3"
   },
   "source": [
    "Question A3.  Now consider a chooser option for which after 𝑀 periods, the holder will choose between an American call or an American put (above we were choosing between a European call and a European put). Modify the function from above so that before populating the chooser tree CH the first loop also determines whether to exercise or keep the respective American put and American call. Call the function that evaluates the price of this particular derivative chooser_american(S0, K, T, r, sigma, q, N, M)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "y4yJ3TPtQBH4"
   },
   "outputs": [],
   "source": [
    "def chooser_american(S0, K, T, r, sigma, q, N, M): #\n",
    "    ## variables as defined above: \n",
    "    deltaT = T / float(N)\n",
    "    u = np.exp(sigma*np.sqrt(deltaT))\n",
    "    d = 1.0 / u\n",
    "    p = (np.exp((r-q)*deltaT)-d) / (u-d) \n",
    "    piu = np.exp(-r*deltaT) * p\n",
    "    pid = np.exp(-r*deltaT) * (1-p)\n",
    "    \n",
    "    C = np.zeros((N+1, N+1))\n",
    "    P = np.zeros((N+1, N+1)) \n",
    "    S = S0 * u**np.arange(N+1) * d**(2*np.arange(N+1)[:, np.newaxis])\n",
    "    S = np.triu(S) \n",
    "    \n",
    "    C[:, N] = np.maximum(0, S[:, N]-K)\n",
    "    P[:, N] = np.maximum(0, K-S[:, N]) \n",
    "    \n",
    "    ## loop for the payoffs \n",
    "    for j in range(N-1, -1, -1):\n",
    "        C[:j+1, j] = piu * C[:j+1, j+1] + pid * C[1:j+2, j+1] \n",
    "        C[:j+1, j] = np.maximum(C[:j+1, j], S[:j+1,j]-K)\n",
    "        \n",
    "        #take the maximum of the remaining value or value at n-1\n",
    "        P[:j+1, j] = piu * P[:j+1, j+1] + pid * P[1:j+2, j+1]\n",
    "        #take the maximum of the remaining value or value at n-1\n",
    "        P[:j+1, j] = np.maximum(P[:j+1, j], K-S[:j+1, j] )\n",
    "        \n",
    "\n",
    "    #create chooser matrix\n",
    "    CH = np.zeros((M+1, M+1)) \n",
    "    \n",
    "    #populate last column with maximum of put/call\n",
    "    CH[:, M] = np.maximum(C[:M+1, M], P[:M+1, M])\n",
    "    \n",
    "    #backwards induction\n",
    "    for i in range(M-1, -1, -1):\n",
    "        CH[:i+1, i] = piu * CH[:i+1, i+1] + pid * CH[1:i+2, i+1]\n",
    "    return CH[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "El-jug_IQBH7",
    "outputId": "aa939667-7400-4fe4-c6e9-d87ea9051acb"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0418690863961153"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S0=9.; K=11.; T=3/12.; r=.025; sigma=.3; q=.01; N=500; M=200; #run function with given values\n",
    "chooser_american(S0, K, T, r, sigma, q, N, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "6vm1FrKmQBIA",
    "outputId": "b32677ae-c362-4997-8351-2c2366244444"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.4 ms ± 1.35 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit chooser_american(S0, K, T, r, sigma, q, N, M) #ipython magic for timing things"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DQ21U4lfQBIF"
   },
   "source": [
    "Question A4. For the given parameters values ($S_0$, $K$, $T$, $r$, $\\sigma$, $q$, $N$) plot the prices of `chooser` and `chooser_american` as a function of $M$. That is, produce a graph that has the prices of the two options on the y-axis and $M$ on the x-axis. Describe and explain (in words) the price convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Y1Qx8_5AQBIG"
   },
   "outputs": [],
   "source": [
    "S0=9.; K=11.; T=3/12.; r=.025; sigma=.3; q=.01; N=500\n",
    "M = range(0, 499)\n",
    "\n",
    "## create empty lists to hold the prices of the choosers \n",
    "Chooser_EU = []\n",
    "Chooser_USA = []\n",
    "\n",
    "for j in M:\n",
    "    \n",
    "    ## apply valuation function for european chooser up until M=499 (must be < 500)\n",
    "    x1 = chooser(S0, K, T, r, sigma, q, N, j)\n",
    "    \n",
    "    #calculate prices for american chooser up until M=499 (should be less than 500)\n",
    "    x2 = chooser_american(S0, K, T, r, sigma, q, N, j)\n",
    "    \n",
    "    ## append values to lists \n",
    "    Chooser_EU.append(x1)\n",
    "    Chooser_USA.append(x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 299
    },
    "colab_type": "code",
    "id": "1ubDP5wzQBII",
    "outputId": "8369a8bb-e680-4711-f1d1-7a182c120a48"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEaCAYAAAAPGBBTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZzN9f7A8dc7EUWRNlnSguxEEZE2RW50W7i3Qsul0i2/9rTo3va62m6rlitUKhRFSYUoSwhjKEnKmLHvy1hm3r8/3t+ZjuPMmOEsc2bez8fjPM6c7/mc7/l8v+fM930+u6gqzjnnXF4OSnQGnHPOFW0eKJxzzuXLA4Vzzrl8eaBwzjmXLw8Uzjnn8uWBwjnnXL48UIQRkTdFREXk2UTnJRIRaRfkr12i8xJPIjIxOO5It+cTnb9kIyIdROQzEVklIrtEZKWIjBaRS0PSPByc34MTmdd4EpFSInKTiMwQkc0iskVEfhCRm0Wk1H7us2JwLk+L8NxEEZl4wBmPsRLzBSgIESkHXBE8vEpE7lbV3YnMUwSzgTOBBYnOSALMA3pH2J4R74wkMxEZANwODAduAVYAxwJ/AT4UkeaqOjeBWUwIESkNfAJcALwEPAAocBHwHNBRRLrsxzWhItAfSMP+f0PdfECZjhMPFHu6FDgcGAt0xL4gnyU0R4Hg14yo6iZgWqLzkyCbVTVmxy4ih6jqjljtvygQkauxIHGnqg4Ie/ojEXkBWB//nMXHPj7j+7H/+y6qOipk+3gR+RYLIvcD/4pWflQ1OX7wqarfghswDlgHHA1sAz6MkOZh7FfGqUH6rcAfwLXB89cAPwFbgAnAyRH28Q9gLpAJrAHeAo4MS6PAY8C9wG9AFtAUaBc81y4s/aXAd8H7bgJmAJeEPH8LMDU4vg1YsLk4bB81g333Bv6N/VLfAHwKVNvHuXsFWAkcHLb9EOzC83zwuDzw3+Cc7Qhe8xVw6j72PxGYUoDPcCkwKMJ2BR6O8Dk2CD7HLcCo4DkB/g/4GdgZnIeXgMPz+Izux34tbge+BZpEeP+/Bud8W3BOPwJqhKXpBnwDrA7y8yPQI49jeRS4NfhubAYmAfULcH4WACkF/H/IOUe1gDFBnn4HHgIOCktbB/g4OLbtwbFeFGGfFwXfw+3ARuziWycszYXA98HzW4LP4aGwNI2B0cF3azv23W8TlmZQ8LmcGexvO/BCHsd6SJD3Mfmcj7HB+x0SPG4XnJ/Lgvdaj/3vvQtUDvufCr/1DPleTyzsuSzMZxONW1R3lsw34HjsYvxq8Pg97EJeKY8PKCX4R70g+FAVeDz4QnbBqrDSgelhr38S2AUMANoD1wLLgelAqZB0GmyfHHwRL8KqB3K+nO1C0v4z2PZxkPZC4D7g1pA0/wGuB84Lnn8peE2HkDQ5X+qlwfF3AHpgwWzSPs5fy+C1HcO2XxZsbxY8fgMLDtcDbbEA9x+g5T72PxGYgpWCw28Skm4phQsUvwL9gHNzzmnwOWpwji7EgsaW4LM4KGyfy7CLVBegK3ZRW0tI4AduDNK+jf1i7QosxC7yFULS9cOqItoD52PBehdwY4RjWYoFuEuAy4N9LSYsUEf4jivweAH/J3LO0XzgjiBPLwTbrg3b72pgCXA1VoX1Bfb/FPr9uijYNj7I99+DPK8GqgZpTsJ+QLwbpD8X++HyVMh+TsN+oE0Jjr0jFjR25HzPgnSDsCD6O/Y/0g5okcextgqOq3c+5yPnczwzeNwu5DvwvyC//wzec0KQ5hDsO55zfWgZ3I4O+V5P3I9zWaDPJmrXx2jvMFlvwD1hX4ILg8fh/6Q5H1D3kG2VgN3YBeLwkO23BmlPCB7XDD7w8F9HrYN0XUK2KRZoyoWlzflytgseHx58MUcW4lgPwi6wXxL8ig7JnxIWFIA7g+3H72O/i4D3w7Z9AiwIeTwfeHY/Pp+JRP5lpsDlIemWUrhAcVtYuiOxHwiDwrZfHaS/JGyfa4DDws7hLuCR4HF57Jfx22H7q4mVVvru4zN6A5gb4Vh+AUqHbLs82N4qn3PYgn1cDPP4rl8btj0F+DLk8X+C7/8pIdtKYUFzdsi2mUG+Dw7ZdmJwvp4NO47D88nX11igLRP2fguBT0K2DQr21bkAx9o1SHthPmkuCtJcGfa/+EVYuquC7eeF/V/dkMf3euJ+nMsCfTbRunmvpz91B35R1anB46+wC3X3PNJ/nvOHqq4HVgHT1NoQcvwU3FcP7i/ALgDvisjBOTesNLEJ+4Ud6gtV3b6PfLfCLkYD80skIs2CXi4rsS/iriA/dSIkHxP2OCW4r7GPvAwFOotIheA9j8RKJYND0vwA9BSRfiLSvJA9SeYCp0e4fV2IfYT7OOxxS+xX4NCw7cOw83Z22Paxqro154GqLsWqCs4MNp2JBfPwzzwN+37kfuYiUktE3heR5djnswu4gcif0XhV3RXyuKCf0f4I/z7MD3uftth3f3HOBlXNAt4HmojI4SJyGFYS+EBDGoNV9TesRJZzXudgxz1MRC4XkWNC3zjocHI2VnWXHXI+BfufDf8f2k3B2hnlANJ8GPb4IyCbP78DhbHPcxmWfl+fTVR4oABE5HSgHjAy6MpWEagAjATOFJHaEV4W3uC3M49tAGWD+5wv/WL+vBDk3A4HKoe9viC9eXJek5ZXAhGpjl1Mj8SKxq2wC+wXIXkLtS7scU7jX6S0oYYEaS4PHncDSmPVCDn+CbwOXIcFjVUi8pyIHLqPfQNsUdWZEW4H0vgafo6PjLQ9uLitDXk+x8oI+1wJVA3+zvnMv2Lvz7whwecnIuWxKpnGWLtUG+wzehsLXOH25zNaFtyfkE+aSCK9V+j7HEnk7+oK7OJaKbhJPumOBAgukBdi16YhwAoRmS4iOYHkSOwX9oPsfT5vASqJSOh1bVVwod2XnHNTM580OedtWdj2Pb4DqppzLahK4RXkXIba12cTFd7ryfQI7u8JbuG6Y13lDtTa4L49kXuWrA17rAXY55rgvir2ayKSi4AjsCJzbkAp4MW5wFT1NxH5Dqum+V9wP1FVl4Wk2YK1n9wnIidgQeVJLKhGOveFlQmUCd0QlGzyzHbY45x/vOOA1JB9HIxd1MM/o2Mj7PNYrH2JkPQ9Q/cXYnNwfyZ2IWqjqlPC3jcqVDVdRBZi9d79orVf7JwdF2H7cdj5XYf9wtZ80uWeV1WdAEwQkUOwatl/A2NEpCbWwJsNvMyeJVVCXp8d+rCAxzATK9Vfgv2QieQSrBoxvIvrHt8BESmDXdCXU3gFOZdxV+JLFMGH2g2r/jknwm0OcI2IFKRoui/jsS95jTx+Gf+2H/v8Hmto7ZVPmpyAkFtVEZSSWu/H++3LEKBdMCDwTPL4ZwZQ1d/VumimYL2PouH3CPvqVIjXT8N+lXUL294V+2E1KWx7x6BaBYDgYtYS69kD9vlsxuqcI33mPwfpIn1GlYDOhch7QTwONBCR2yM9KSJNRaSwVReTgJbBsefspxR2zn5U1c1B9dws4IrQ6sbgx0Ir9j6vqOoOVf0GeBo4DDgx2M9krOQ1O9I5LWTec98LeBH7PPc658G2DlivqfDutVeGPb4Cu7bmfAdy0pcrQFb2eS4LsI+o8xKFXUQqA3eo6sTwJ0XkdeBVrOFqwoG8kar+KiJPAS+JSB3sS5GJtWFcALwZ/JoqzD43i8h9wH9FZARWzbMZaAJkqup/sWqP3cDgYLBVFawv+B9E/8fCh9g/3FCsa9+I0CdFZCrWQyUFC3BnY//07xRg3xVEpGWE7etDLrjDgLdF5Dmsbrox9mu+QFR1XTAq/z4R2Yp1iayLdUedwt51wtuBL0XkGayK6F/YL9Pngv1tEpG7gJdF5GisbWsjVgI8GytxvYcFlE1Buv7YhfEBrMR4REHzX4DjGxqMEB4gImdin9cKrIrsYqx7d3Psu1FQz2HneHyQ901Y763awT5zPIidv89E5BWsbe1f2PkYACAiN2L19GOxKp6jsBJoOn+WmG/HuiGPE5G3sKqao7A2kFKqem8h8h7q39ixfygiL2OflfJnb6YvsO9BuPoi8j/su1cb6zI9SVVz2s5WYiWmbiIyD+ux9ZuqhpdOoeDnMr6i3TqebDdgVPBhHJrH80dgfd8H6Z69DcLHCywFhoZtaxekPT9s+zXYL9et2MVyIdYVs1pIGgUejZCfnH22C9t+OVYq2h4cz3SgU8jzV2KNp5lYFUg3rFfI0pA0NYnQOyOv98znnH4UpH8vwnNPYeMDNgbHn0JIN9589jmRvHs9fRaS7iCsL/nvwec2DjiZvHs97dWdlMjjKF4m73EU/bA2okzs126kcRQdsR8am4LPaDHW/lAvJM25wbnZjnXbvTUnnxHe99GwbTmfXc8CfkYdsYv2aqwUsxL7X/jLvs5R+Pcm2FYH6+G2MTgPBR1HMYqQcRRYKXQUFiR2BOf+I/Yea1EXuzCvCtKlYT9AOoblM62Q14ODgT5Y+9nW4DYTa/8IPw/tgvPz1+C9NmA/0t4DjgpL2wUbw7Ir9HMi73EU+Z7Lwnw20bhJsHPnXCGJiAKPqWo02q9ckgmqVycAF6jqVwnOTkyV+DYK55xz+fNA4ZxzLl9e9eSccy5fXqJwzjmXr2LXPfaoo47SmjVrJjobzjmXVGbNmrVGVY+O9FyxCxQ1a9Zk5sz9GnPjnHMlloj8ntdzXvXknHMuXx4onHPO5csDhXPOuXwVuzaKSHbt2kVaWhqZmZmJzopLoLJly1KtWjVKly6d6Kw4l1RKRKBIS0ujQoUK1KxZk+hMAuuSjaqydu1a0tLSOPHEExOdHeeSSomoesrMzKRy5coeJEowEaFy5cpeqnRuP5SIQAF4kHD+HXBuP5WYQOGcc8XWhg3wzDPopwVZHrzwSkQbhXPOFVu//MLODpdQ5tef+KnSmdRdV5gFHQvGSxTOOZesJk9GmzalzK8/8QY3MP76D2LyNiUyUBx3HIhE73ZcpKXQI3j99depUqUKTZo0yb2lpKTE9mBjZMWKFXTr1o2TTz6ZevXq0bFjRxYtWsTSpUtp0CBay1875/I0aBC0bYts3cpNvMKSu1/n1meqx+StSmTV08qVidnfvHnzePTRR7n++usL/R5ZWVmUKlVq3wnjQFW59NJL6dGjB8OGDQNgzpw5rFy5kurVY/NFLUzeVJWDDiqRv4FcSbB7N9x7LwwYAEBbJlH972155cnYvaX/N8VRSkoKTZo02Wt7y5YtWbp0KQDLly+nefPmAFxxxRXcfvvtnHPOOTzxxBP89NNPtG3blvr163P++eezZs0aALp160bXrl1p0aIFJ5xwAmPGjMnd92+//Ubnzp1p3rw5Z5xxBj///DMAw4cPp2XLljRu3JizzjqL1atXA3DppZfywAMP0KZNG4477ji++mrvFR4nTJhA6dKlufHGG3O3NWnShDZt2gAW1P7xj39Qv3592rdvz/bt2wF49tlnadCgAQ0aNOD555/PfW2k7Vu3buXiiy+mcePGNGjQgA8+sCL10KFDOeOMM2jSpAm9e/cmKyuLpUuXUrduXW6++WZOO+00li1bth+fjnNJYP16uPRSGDCA3RUqUvOgP9jYqC3PPWe1GzET7UW4E31r1qyZhluwYMEejyH6t4I48sgjtWHDhtq4cWNt3Lixvv7665qdna1VqlTR7OxsVVUdO3as9uzZU1VV69Spow8++KCqqmZmZmq9evV09uzZqqr65JNPar9+/VRV9dRTT9V7771XVVUnT56sp59+uqqq7ty5U88991xdvHixqqqOGTMmd99r1qzJzdfDDz+sL730kqqqnnLKKfrMM8+oquqIESNy04d64YUXtG/fvhGP8bffftNSpUrpjz/+qKqqV1xxhQ4ZMkRnzpypDRo00C1btujmzZtzjyWv7cOHD9cbbrghd78bNmzQBQsWaKdOnXTnzp2qqnrTTTfpO++8o7/99puKiE6dOnWfn0H4d8G5pPHrr6oNGqiC7mjWUqscukHLllX95Zfo7B6YqXlcV0tk1VMiLFu2jGOOOYZ58+btsf2XX37hxBNPzO3jP2/ePBo2bEhmZibr1q3joYceAuCTTz7hrLPOomnTpgDUq1eP0aNHs337dtasWUP//v1zt69fvz73NampqVx22WUA7N69O/dX/6BBg/jggw/YsWMHK1as4PHHH2fbtm1s3LiR//u//8tNX7FixUIf64knnphbcmrWrBlLly5l7dq1XHrppRx22GEA/PWvf2Xy5Mm51Vjh2y+66CLuvPNO7rnnHjp16kSbNm0YMmQIs2bN4vTTTwdg+/btHHPMMbRt25YTTjiBli1bFjqvziWFH36A9u1hwwYyr7qOtnNfYsX2cnz/PZxySuzf3gNFnMybN49TTz11r+0pKSk0bNgw9/HMmTPp3bs3qamptGjRgoMPto9owYIFe6RLSUmhXr16zJ8/n1q1alG2bFkAZs+eTePGjQGYO3cujz322F5tIoMHD2bGjBl88803lC9fPrc6KzU1lWbNmuW2hcybNy9iw3T9+vUZPnx4nsd6yCGH5P5dqlQptm/fjuax5G5e22vXrs2sWbMYO3Ys9913H+3bt6dSpUr06NGDJ554Yo+0S5cuzQ00zhU7H30EV14JQPYz/+HmlL78ML8UAwZAvH4beRtFnKSkpEQMFOvWraNcuXIALFy4kDFjxtCwYUNSUlJo1KhRbrqqVauyYMECAJYsWcKQIUPo3r07c+fO5Y8//iAzM5OtW7fSv3//3BJBlSpVGDduHNnZ2bl5UFVSUlJo1aoV5cuXZ8SIEXz//fc0bNiQ+fPn79GGMm/evD3ykOPcc89lx44dvPHGG7nbfvjhByZNmpTn8bdt25ZPPvmEbdu2sXXrVj7++GPatGmT5/b09HQOPfRQrr76au68805mz57Neeedx/Dhw1m1alXuufv99zzXWnEuuWVlwSOP5AYJxozhse138L/Bpfjb3+D22+OXlRJZojj22Oj2fDr22H2nSUlJYdKkSXz++eeATScxefJkLrzwQl588UWuvPJKGjRoQOXKlTn22GNJSUnhjDPOyH39Nddcw9ixY2nYsCHlypXj7bffpnLlysydO5errrqKdu3asWnTJvr160fr1q0BuO6665gwYQJ169alXLlyNGjQgKFDh9KjRw86d+7M8OHD6dixIyeddBKHHXYYKSkptGjRIvc958+fH7FEISJ8/PHH9O3blyeffJKyZctSs2bNPRqow5122mn07Nkz95huuOGG3Gq0SNvHjRvHXXfdxUEHHUTp0qV59dVXqVevHo8++ijt27cnOzub0qVL8/LLL3NcQfsnO5cstmyBm2+GIUOgUiWYPJkPU+vz0ENw9tnwzjvxzY7kVfRPVs2bN9fwpVAXLlxI3bp1E5Sj2Grbti1vvPEGderUSXRWkkJx/i64YmL1aujcGaZOhSZNYMwYlmQeT5MmcOSR8PXXcPLJ0X9bEZmlqs0jPRfzqicRqS4iE0RkoYikishtEdKIiLwoIotFZJ6InBby3NPB6xYGaXxmtxC//vortWrVSnQ2nHPR8Ntv0KiRBYnLL4eJE9lW8Xhat4atW+H992MTJPYlHm0Uu4E7VLUu0BLoIyL1wtJ0AGoFt17AqwAi0gpoDTQCGgCnA2fHIc9JY/ny5T64zLni4Lvv4KSTYMUK6NcP3n+f7WWOoEMH2zRiBJx5ZmKyFvMrjKpmqOrs4O/NwEKgaliyzsDgoDvvNKCiiFQBFCgLlAEOAUoDUR5X7ZxzCaRqjQ5nnWWj5gYPhsceQ0sdzL/+Bd9+C3feabVRiRLXn6IiUhNoCkwPe6oqEDqcNg2oqqpTgQlARnAbp6oLY59T55yLg6wsKz307AkVK1pUuOYaAN56C556Ctq1g6efjvHI632IW68nESkPjAD6quqm8KcjvERF5BSgLlAt2DZeRNqq6rdh++6FVVlRo0aN6GbcOediYfdum47js8+gXj34+GOoXRuAceOgTx9o0AC++CKxQQLiVKIQkdJYkHhXVUdGSJIGhM4mVw1IBy4FpqnqFlXdAnyOtXPsQVUHqmpzVW1+9NFHR/8AnHMumrZtgzp1LEh06AATJ+YGifR0uPZa6+H07rsQMn41YeLR60mAt4CFqvpsHslGA92D3k8tgY2qmgH8AZwtIgcHweZsrI3DOeeS05IlcNhhdn/LLfDJJxD8wM3KgnPPhYwMq3qKMN41IeJR9dQauAZIEZE5wbZ+QA0AVX0NGAt0BBYD24Brg3TDgXOBFKxh+wtV/TQOeXbOuegbP97mbAJ46SWrXwpkZcFf/gI//wz/+x907JigPEYQ80ChqlOI3AYRmkaBPhG2ZwG9Y5Q155yLn//+F269FSpUgE8/tSHWIZ5+Gj7/HG68Mbc9u8gokVN4OOdcXN10E7z2GtSqBaNHQ9i8byNHWuenxo2toFFE1ijLVTJHaiVqLdTALbfcwgknnBCjgzPff/997tTjseBLoTpXANu32/iI116D886zQXVhQeLbb+Gqq6BmTZg+vegFCSipJYpErYWKrTg3ceJEdu7cyebNm6lQoUJ084KtMNeqVStatWoV9X1D0V4KFXw5VFdELFkCbdvC8uXQuze8+CKUKbNHkvR0uO46KF/eShVFoYdTJP6fFGf9+/fngQceoF69eqSmpuZuv+KKK7jllls466yzOOGEE5gyZQrdu3endu3ae6wnkdfSpuHLpl5xxRVMmTIFgPT0dC677DKaNm3KqaeeyowZM5J2KVTw5VBdEpg4ERo2tCDx/PNWoggLEqrQpQv8+iu88goEkykXTXktfZest4IshZqotVDnz5+vTZs21ezsbO3Tp4+++eabuc/VqVNHBwwYoKqqDz74oNauXVvT09N1x44dWrFiRc3MzMx3adPQZVNVbXnUDRs26K5du7RRo0b66aefqqrq1q1bddOmTUm5FKqqHvByqL4Uqou5t9+2a0K5cqrB/1247GzVLl0s2csvxzl/ecCXQi0a7r//fh555BFEhLp16zJ//nwAMjMz2bBhA3379gWgXLlyXH/99VSpUgWAQw89lDJlyjB8+PCIS5uGL5uamZnJrl27OOKIIxg+fDh169alU6dOufsCGDBgQNIthQrw9ddf+3KormhShXvugWeegerVrWdTsNpkeLL//MeGT/TubbeizgNFnEyfPp1x48YxZ84c+vTpQ2ZmZu7qcampqZx22mm5depz587lpptuAiAtLY3jjz8eEclzadNZs2btsWxqamoq9erZBL1z5szZ6+KZrEuhPvTQQ6iqL4fqip6tW6FrVxgzBs44w3o25bGi2fjxcPfdVjP1wgtFs/E6nLdRxEm/fv347LPPWLp0KUuXLmXu3Lm5JYqUlJTcda5hzyVI586dm/t3fkubhi5ZGvr4uOOO26MtZPXq1Um7FCrgy6G6omfpUpv/e8wY+PvfrWdTHkFi6lTo1AmqVIEZM4pu43W4khkoCrJ2aRT3N378eHbs2MF5550X8pJj2bp1K+vWrSMlJSX3Ap2Zmcn27dupVKkSsOfF+rrrriM7O5u6devSpEkTnnrqKUQk30DRs2dPVq5cSf369WnSpAlTp06lR48evPjii7Rp04ZFixbtsRRqaKDY11Ko48eP5+STT6Z+/fo8/PDDHH/88Xmeg9ClUFu0aJG75Gle23OWgm3SpAmPPfYYDzzwAMAey6E2atSICy64gIyMjHzPv3MxM3WqBYmUFHj8cZuc6eDIFTW//26TxJYvD19+CWXLxjerB8KXQnUlin8XXNR88IENoc7OhmHDbEW6PGRnQ5s28P33tkpdt25xzGcBJXQpVOecK1ZU4d//tqt9hQpWh5RPkAD4618tSDz/fNEMEvvijdnOOVdQW7dCr17w3nu2hsTEibkzv+bl0Udh1CibKPaf/4xPNqOtxJQoilsVmys8/w64A7Jsmc38+t57tuDQvHn7DBKffw4PPgitWtlqdck6WUCSZrtwypYty9q1a/1CUYKpKmvXrqVsMrUguqJj1qw/Gxn697f5NvbRr3X6dJs2vEoV6xIbDGFKSiWi6qlatWqkpaXlTlPhSqayZctSrVq1fSd0LtSoUdbtdds2+PBDuOKKfb5k0SJrizj8cOstm8xBAkpIoChdujQnnnhiorPhnEsm2dnw7LNw111QrhzMmRNxpHW4bdusG2xGBnz0ERSHS0+JCBTOOVcoW7fCbbfZeqQnnWRVTwWczubKK214xQsvWNVTcVAi2iicc67AMjLg4ostSFxyCSxeXOAgcfvtNkD7vvuSt4dTJB4onHMux7x50Lo1TJpkjdajRtniZAXw3nvw3HPQubP1dCrgy5KCVz055xxYUaBrV6t2KmCjdY5p0/5cpW7oUGvSKE68ROGcK9lyGq2DqfiZO7dQQWLOHLjoIjjqKJg92+ZyKm68ROGcK7m2bYNbb7X2iNq1rRX6yCML/PKMDOvhpAqffQbBXJ7FTsxLFCJSXUQmiMhCEUkVkdsipBEReVFEFovIPBE5Ldh+jojMCbllikiXWOfZOVcCZGTAhRdakOjcGebPL1SQULXqprlzbTnsFi1imNcEi0eJYjdwh6rOFpEKwCwRGa+qC0LSdABqBbcWwKtAC1WdADQBEJEjgcXAl3HIs3OuOJs3z3o2paVZo3X//oVufe7aFSZMsNXqevSIUT6LiJgHClXNADKCvzeLyEKgKhAaKDoDg4N1W6eJSEURqRK8NsflwOequi3WeXbOFWOjR9tVfudOGDHCpnYtpMces8F0ffpAsHJwsRbXxmwRqQk0BaaHPVUVWBbyOC3YFqob8H4e++0lIjNFZKZP0+GciygrC554wqqZDj3UShX7ESQ+/RQeeADatoWnn07eif4KI26HKCLlgRFAX1XdFP50hJfkzuAnIlWAhsC4SPtW1YGq2lxVmx+9j9kcnXMl0JYt0L079OsHTZvaILr69Qu9mylTbOLY6tVtZthkn8OpoOISKESkNBYk3lXVkRGSpAHVQx5XA9JDHl8JfKyqu2KXS+dcsbRsGbRrZyPirrnGBj3sR/eklBSrsTrqKOscVVKCBMSn15MAbwELVfXZPJKNBroHvZ9aAhvD2if+Rh7VTopK31UAACAASURBVM45l6fp06FZM5uracAAeOcdKFOm0LtZt84arLdssXhTNbxivJiLR6+n1sA1QIqIzAm29QNqAKjqa8BYoCPWq2kbcG3Oi4N2jerApDjk1TlXHKjaEOlrr4WDD4avvoLzztvvXXXpAj/+CIMGwbnnRjerySAevZ6mELkNIjSNAn3yeG4pezdsO+dcZDt3WlvEgAFQo4bN21Sz5n7tStWqmyZPtvWui3s32Lz4yGznXPGxZo21Q3zxBZx/PnzyCRx22H7tShWefNK6wd58s615XVKVgI5dzrkSYeFCaNnSgsTdd9v9fgYJgI8/toJJ69ZWmtjHyqfFmpconHPJb8wYWzFoxw4YNszqiw7Al1/avIC1atno69Klo5TPJOUlCudc8srKsvqhTp2sv+qcOQccJGbNson+qlf3IJHDSxTOueS0YQPcdJOVIM44w6qaDnD61rVr4e9/t4LJBx+UvG6wefFA4ZxLPosXw+WX29StvXvDSy9ZN9gDsHs3nHMOLFpkbeBt2kQpr8WABwrnXHL56itrj1i/Ht58E6677oDXHd2922qvUlJsTF7nzlHKazHhbRTOueSQlWXdj9q3t1XpZs2C668/4CCRnQ333QfjxsH991vvWrcnL1E454q+jRttJbrBg6FRI2tlLsQiQ/l59VVbU6JzZ3jkkQOOO8WSlyicc0XbokVwwQUWJK691ubSiFKQGDLEBtK1bm3jJjxIROaBwjlXdI0bB2edBT/8YO0Rb78dtQUgxo+3INGsGYwa5UEiPx4onHNFT3Y2PPMMdOxoLc0zZ1p7RJQsWWLdYCtVsoJK5cpR23Wx5G0UzrmiZf16+6n/3nv2c//rr+GII6K2+40boVUr2LTJFh+qVy9quy62vEThnCs65s+3yfzeew9uvNGqnKIYJDZssJqslSvhm2+gefOo7bpY8xKFc65oGDkSevWCrVttMEP37lHd/datNuQiNdV237p1VHdfrHmJwjmXWNnZ8MADNgtfhQq2zmiUg4Qq3HWX9Wy6//6o777Y8xKFcy5x1qyx1YDGjrUqp+HDo1rVlOO++2y8xA032FgJVzheonDOJcYPP8CZZ9pkfvfdZ3N7xyBIvPQSPPWUTSr7yitR332J4IHCORd/AwfaDHybNlkp4vHHYzKQ4YsvoG9fOPtsCxI+Zfj+8aon51z8ZGVZ/c+gQXD66TaI4dRTY/JW8+bZRH9Vq8KHH0ZtMHeJ5IHCORcf6ek2oV9qqg2ee/55KF8+Jm81f751g61QwWq4jjkmJm9TYsS86klEqovIBBFZKCKpInJbhDQiIi+KyGIRmScip4U8V0NEvgxev0BEasY6z865KPv8c6hTx4ZEv/aaTccRoyCRs1RFqVLw7bceJKIhHm0Uu4E7VLUu0BLoIyLhYyE7ALWCWy/g1ZDnBgPPBK8/A1gV+yw756IiKwvuvdem4jj+eJv1tXfvmL3dhg1w9dWQlmYT/jVsGLO3KlFiXvWkqhlARvD3ZhFZCFQFFoQk6wwMVlUFpolIRRGpAlQCDlbV8cHrt8Q6v865KFm92ubunjoVLrsMXn4Zjj02Zm+nChdfDNOnWxNIp04xe6sSJ669noJqo6bA9LCnqgLLQh6nBdtqAxtEZKSI/Cgiz4hIqQj77SUiM0Vk5urVq2OTeedcwU2ZYlVNU6fCgAG2rnUMg0R2tgWG77+33k09esTsrUqkuAUKESkPjAD6quqm8KcjvESxEk8b4E7gdOAkoOdeCVUHqmpzVW1+9NFHRzXfzrlC2LULnnjCFpwuXRq++w5uv/2A17POT1aWrWk0diw89JBNEeWiKy6BQkRKY0HiXVUdGSFJGlA95HE1ID3Y/qOqLlHV3cAnwGkRXu+cS7Tly6FLF+jXD847z7oetWoV07dUtQLLyy9b20T//r6uRCzEo9eTAG8BC1X12TySjQa6B72fWgIbg7aNH4BKIpJTTDiXPds2nHNFwcSJ0LKl/ax/9FFbcCgOpfuBA+GeeywuDR4ctTWNXJh4jKNoDVwDpIjInGBbP6AGgKq+BowFOgKLgW3AtcFzWSJyJ/B1EHBmAW/EIc/OuYLIyrL5MR56CCpWtKqmGJcicgwZAn36WC3X5597SSKW4tHraQqR2yBC0yjQJ4/nxgONYpA159yBSE+Hf/zDShHnnGNTs8ZgrqZIvvwSbrrJ1pP48EOfmiPWvKDmnCu8L7+0qqavvoJ//9tWAYpTkEhNtd62xx0H775r9y62PFA45wpO1RqrL77YfsaPGwcPPhi3t//jD4tPIvbWJ58ct7cu0XyuJ+dcwaxaBX/5C8yYAX/9qy3wEMf5MX77Ddq2tWaR2bM9SMSTlyicc/v26adw0kkwZw688IJNDR7HIJGWZgvgrV9vTSIxmnDW5cEDhXMub1u2WNeiSy6B6tWtV9Ott8a1i9HmzXDVVTYs4803oV27uL21C3jVk3MusjlzbBRbaqqtIfHMM9YFNo6ys62269tvbcxEt25xfXsX8EDhnNtTdratH3r33TZX9/Dh1iYR54EKOfM3TZpk2fnHP+L69i6EBwrn3J/S022AwujR0KIFjBoV08n88rJ7t83Z9Pnn8MgjcPPNcc+CC1HgNgoRqSUib4vIy7HMkHMuQT79FM44w/qdPvywzfyagCCRnW09bt96ywLE/ff7qOtEK0xj9hDgI2w2V0SkgYgMjkmunHPxs3On/Xy/5BIbNPf11wmdXe+JJ+DJJ62X00sveZAoCgoTKA5S1c+BLABVnQ80iEmunHPxMW+e9TV9/XXo1ct6NbVunbDsPP00PPCATUL7wQceJIqKwgSKdBE5EVsnImdW2HIxyZVzLrZ27LCf7U2awMaNMHKkBYs492oK9frrNuj7L3+BoUM9SBQlhWnM7gu8CRwnItcCFwHzY5Ir51zspKZag/XkyXDRRfD221ClSkKzNGLEnzPB/u9/cNhhCc2OC1PgEoWqLsWCw63YSnOTsOnDnXPJQNUq/Vu1glmzbGDC2LEJDxLffGPtEaeeagGjcuWEZsdFUOAShYg0xxYfGg4Mj12WnHNRt3KlDZ776itrg3j3XTjhhETniokToX17qFHD1rs+/PBE58hFUpg2isEEDdkAInKUiHSKfpacc1GTnW0r/JxyijVUP/mkDXMuAkFiwgTo3NmyMmOGB4mirDCBIlNVM3MeqOoa4N/Rz5JzLiqWLbPuQ927Q716FijuuadIrBc6dSr87W82r+DYsXGdX9Dth8J8Y5aISIewbWWimRnnXBSo2gLSTZrAmDHw2GM2D0bTponOGQA//WQlidKl4aOPoE6dROfI7Uthej39E/hcRK4BpmFjKH6NSa6cc/tnxQro3dum4GjYEKZMgbp1E52rXMuWWVt6VpY1lzTyRY6Twj5LFCLSP/jzRGxU9gjgKGAO8LfYZc05V2CqMGiQBYXx420KjjlzilSQWLoUmjWD7dth2jQPEsmkICWKL4P7vkB9oDSwINjWEe8B5Vxi/fqrlSK+/tp+rr/0UpGpZsrxyy/QoYMFiSJWyHEFsM8ShapOFZGDgDmqWh9oiDVi/wq03NfrRaS6iEwQkYUikioit0VIIyLyoogsFpF5InJayHNZIjInuI0u1NE5V5xlZcHzz9tP82+/heees0EJRSxILFkCl10Ga9bAJ59YqcIllwK1UahqtoicBzyuqjuA2cGtIHYDd6jqbBGpAMwSkfGquiAkTQegVnBrAbwa3ANsV9UmBXwv50qG+fNtgYZp0+Dss63xukaNROdqL6tW2VyDv/8O778P552X6By5/VGYXk9zRKR/ULooMFXNUNXZwd+bgYVA1bBknYHBaqYBFUUkscNFnSuKtm61CZGaNrWf6m+8YaPWimCQ2LwZzjoLfv7ZJvjr5KOuklZhLvrVgW7Y5ICjROQREbmiMG8mIjWBpsD0sKeqAstCHqfxZzApKyIzRWSaiHTJY7+9gjQzV69eXZgsOZc8vvjCejI98QRceaVNw3HDDYnOVUQbNkDz5rB4sY2T6Ngx0TlyB6LA3WNV9UoAETkEa9RuCJyBrVGxTyJSHusx1VdVN4U/Hektg/saqpouIicB34hIiqru0S1XVQcCAwGaN2+uOFecrFoFffta3U3Vqraw0AUXFNnpVVetsobrX3+Fzz6zrLrkVpDusT1EZI2IrBORd4AyqjpbVd9R1bsK8iYiUhoLEu+q6sgISdKwEkuOakA6gKrm3C8BJmIlEueKv1274JVXoHZt+PhjW8N6yRKbHKmIBokVK2x57dRUi2tekigeClL19CBwAXAq8AfweGHeIFi34i1sQsFn80g2Guge9H5qCWxU1QwRqRSUYBCRo4DW/Nk117ni6/vvbc3qPn3gtNNszounnoIyRXcyhHXr4PLLYeZMazq5olAV064oK0jV0yZV/TH4+0ERCW9f2JfW2HTkKSIyJ9jWD6gBoKqvAWOxMRmLgW3AtUG6usDrIpKNBbUnw3pLOVe8rF8Pd91la0RUqmQr+HTrBqVKJTpn+dq2zQo6c+bYuL+rr050jlw0FSRQVBGRXlhvpZ+wAXcFpqpTiNwGEZpGgT4Rtn+PtYU4V7zt3AlvvQX33mtX3ZtvthJEEqzgs3UrnHkmpKTA8OE2ZsIVLwUJFP2BRsBV2EW7vIiMBeYC81T1/Rjmz7ni76uvrBQxZ46NiRgwIGlGpa1fb2MjUlJsNdVLL010jlws7DNQBD2KcolINSxwNMSqizxQOLc//vjDejN9/DEcfzwMG2YV+0VgGvCCWLnSGq5TUizrHiSKr8LMHguAqqZhvZTGRj87zpUA69fDf/4DTz9tc23ff7/dypVLdM4KbMUKa7ieNcvWuO7aNdE5crFU6EDhnNtPWVm2TvWjj0J6ul1dH38cTjop0TkrlHXrbJR1Soq1uV91VaJz5GLNA4VzsaZqCwjdcQcsWmRrVo8cad1fk8ymTZb9RYus4dqrm0qG5KgMdS5ZzZxprb1/+Qvs2AEffgiTJydlkFizxoZ0/PyzTcvhQaLk8EDhXCykplol/umnw8KFNgX44sXWWF1ER1Xn548/rCSxbBl8+SVceGGic+TiyauenIumtDR48EEbKHfwwfDQQ9azqVKlROdsvy1aBF26QEYGfPopnH9+onPk4s0DhXPRsHy5DZB75RUbRX3jjbYcaeXKic7ZAVm4EDp3hrVrrU2ifftE58glggcK5w7E8uXwzDM2udHOndCzp3V1PfHEROfsgP38s838un27tb2ffXaic+QSxQOFc/tj5Urr2jpwoDVS9+gBDzwAJ5+c6JxFxU8/2bQcqjBhQpFbXdXFmQcK5wpjyRKbYmPgQGuU7toV+veHU05JdM6iZtYsaNsWDjkEfvihWB2a208eKJwriJQUa4P48EP7md2zp83PVLt2onMWVV9/bW0SlSrBd9/BCSckOkeuKPBA4VxeVGHKFAsQY8ZA2bI2q+tdd9lKc8XMyJHQvbstvz12rAcJ9ycfR+FcuF27rItPs2ZWB/PDD9bNNS0Nnn++WAaJIUNsKo569Wxp7po1E50jV5R4icK5HBs2wOuvWzBYsQJq1bLeTN26Qfnyic5dzLz2mi2k16oVjBgBxxyT6By5osYDhXMzZtj4hw8/tL6g7dtbwLj44iK/styBULVJbO++G849F0aNKtbx0B0ADxSuZNq61X4+v/aarUddvjxccw3cdBM0aZLo3MVcVpYN93jqKZuGauRIG0juXCT+1XAly9y5FhyGDoUtW6x66cUXrRdThQqJzl1c7N4NvXvbFOE9e9p9Ek4/5eLIA4Ur/hYtgg8+sGXYFiyw3ktdu8INN9hMdyXoKrlrl83bNHYs3HMPPPFEiTp8t588ULji6fffrc1h2DCYPdu2tWkDL71kjdNJPgfT/tiyxWY8nzHD2utvvdWDhCsYDxSu+MjIgI8+suAwdaptO+MMePZZm967WrXE5i+BMjIsSPz8M7z3Hvztb4nOkUsmMQ8UIlIdGAwcB2QDA1X1hbA0ArwAdAS2AT1VdXbI84cDC4GPVfWWWOfZJZE1a6xRetgwmDTJuvI0bmx1KldemXTLjMbC/PnWYL16NXz2GXTokOgcuWQTjxLFbuAOVZ0tIhWAWSIyXlUXhKTpANQKbi2AV4P7HI8Ak+KQV5cMfv8dPv/c+nOOH29deE491eZc6trV/nYATJxoNW0iMG6cNck4V1gxDxSqmgFkBH9vFpGFQFUgNFB0BgarqgLTRKSiiFRR1QwRaQYcC3wBNI91fl0RtGWLTaXx5ZcWIH76ybafdJINAujaFRo18gr3MMOHw7XX2kDyjz6Chg0TnSOXrOLaRiEiNYGmwPSwp6oCy0IepwFVRWQlMAC4Bjgvn/32AnoB1KhRI3oZdomxbRt8/73Nbz1hgk2hsXu3TWd69tnWt7NDB5uQz4NDRC+/DLfdZmtcjxxZoptnXBTELVCISHlgBNBXVTeFPx3hJQrcDIxV1WWSzwVBVQcCAwGaN2+u0cmxi5vly60rzowZNmXptGnWj7NUKVtz+q674JxzrN7k0EMTndsibfdum5bqiSdsgPmIET7a2h24uAQKESmNBYl3VXVkhCRpQPWQx9WAdOBMoI2I3AyUB8qIyBZVvTfWeXYxoAp//GFTds+dCzNnWnBIT7fnS5e2UdF9+1pgOOusEjMILhq2bLGB5UOH2kC6N98s1jOQuDiKR68nAd4CFqrqs3kkGw3cIiLDsEbsjUHbxlUh++kJNPcgkQR274alS+GXX+z2008wb54FiE0hhclatSwgtGhh3VgbN7bBcK7QVq2ypppJk6xN/+GHE50jV5zEo0TRGmtjSBGROcG2fkANAFV9DRiLdY1djHWPvTYO+XL7Y/t2KwFkZNh9zi3n8bJl8NtvFixyHH64NTZfdZXdN2oEDRrYdnfAFi+Gjh0tNr/zjk1Z5Vw0iXU0Kj6aN2+uM2fOLPwLt22DF4LhHSJ73hKx7aCD9ryVKrX/20SsC+nu3XvfduywOovNm+0+5xb6ePNmG6+Qnm5TcYc75BA4/nioUsW62JxyipUWcm7HHOONzjEybZq162dl2dpKbdokOkcuWYnILFWN2LPUR2bn2LIF+vVLdC4S76CDrPUz51ahgt3XqWNzUR9//J9BIefvSpU8ECTAsGFw9dX2EUyc6GMLXex4oMhx9NGQmWkNrqE3SMw2VcjOtltW1p9/7882VZtDOtKtTJk9A0K5cn7RL+Kys+Ff/4J//xtatrShJRUrJjpXrjjzQJFDxKpQnCvCtm+3da2HD7fSxFtvWax3LpY8UDiXJFauhAsvtJ7FTz5pg9K98OfiwQOFc0ngxx9tAN3mzfDpp9CpU6Jz5EqSgxKdAedc/t5/34aZlCkDs2Z5kHDx54HCuSJq926bveTvf7cxiXPnQv36ic6VK4m86sm5ImjDBrj8cvj6a5sD8YUXvK+FSxwPFM4VMQsW2CC6tDR44w24/npvtHaJ5YHCuSLko4+s22uFCjbTeosW+36Nc7HmbRTOFQE7dsCdd9rqrY0bQ2qqBwlXdHigcC7B0tOt6+uAATZN+JQpcOyxic6Vc3/yqifnEujbb60UsWEDDBli1U7OFTVeonAuAbKy4OmnbTmOQw+F2bM9SLiiy0sUzsXZqlXQqxeMGgVdutiKdIcdluhcOZc3L1E4F0eTJ9uMr198Ac8+CyNHepBwRZ8HCufiQBWeeALOO8/GREyaBP/3fz4+wiUHr3pyLsZWrbL2h/Hj4bLL4M03ff0Il1y8ROFcDI0bBw0bWpfXl1+2AXUeJFyy8UDhXAxkZkLfvnDRRTYm4vvv4eabvarJJSevenIuylJSoFs3m7PplltskSFvsHbJLOYlChGpLiITRGShiKSKyG0R0oiIvCgii0VknoicFmw/QURmicic4LU3xjq/zu2vrCzrydSsma1GN3Ys/Pe/HiRc8otHiWI3cIeqzhaRCsAsERmvqgtC0nQAagW3FsCrwX0G0EpVd4hIeWC+iIxW1fQ45Nu5AluyBG64ASZMsIWFBg2CypUTnSvnoiPmJQpVzVDV2cHfm4GFQNWwZJ2BwWqmARVFpIqq7lTVHUGaQ+KRX+cKQxVeeQVOO81GV7/xBowe7UHCFS9xvfCKSE2gKTA97KmqwLKQx2nBtpyqq3nB809FKk2ISC8RmSkiM1evXh2LrDu3l+XL4eKLoU8fOP10CxQ33OAN1q74iVugCKqORgB9VXVT+NMRXqIAqrpMVRsBpwA9RGSveTVVdaCqNlfV5kcffXS0s+7cHlThtdegTh2ranrpJfjySzjppETnzLnYiEugEJHSWJB4V1VHRkiSBlQPeVwN2KPkEJQkUoE2scqnc/uyZImtPnfTTXDmmTB/vpUovBThirN49HoS4C1goao+m0ey0UD3oPdTS2CjqmaISDURKRfspxLQGvg51nl2LtyuXfD441C/Pnz3nZUixo2Dk09OdM6ci7149HpqDVwDpIjInGBbP6AGgKq+BowFOgKLgW3AtUG6usAAEVGseuo/qpoShzw7l+vHH23d6h9/tCk4nn8eqlVLdK6ci5+YBwpVnULkNojQNAr0ibB9PNAoRllzLl/r1sEDD8Drr8NRR8GIEfDXvyY6V87Fn3c3dS6MKrz3HtStCwMH2ujqn3/2IOFKLp/Cw7kQs2fb9N/ffmtdXseNgyZNEp0r5xLLSxTOAWvXwo03QvPmNkfTq6/C1KkeJJwDL1G4Em7HDhtZ/eijsHEj3HYb9O/vU4E7F8pLFK5EUoUPPrB2iNtvt5LE7Nnw3HMeJJwL54HClTjffgstWthU4BUqWDvEuHHQyPvXOReRBwpXYsybB507w9lnQ3o6/O9/Vopo3z7ROXOuaPNA4Yq9hQuha1do3BgmTrT2iEWLoGdPKFUq0blzrujzxmxXbM2ZY9NuDB9uiwfdf7+1Rxx5ZKJz5lxy8UDhip3vvrMAMXYsHH443HuvBYijjkp0zpxLTh4oXLGgCuPHw2OPWWP1UUfZ3zff7L2YnDtQHihcUtu9G0aNgieegFmzoGpVm7Tvhht8rWrnosUDhUtKa9bAm2/aCOo//rDpvt94A665Bg45JNG5c6548UDhksrMmbYWxLBhNqr6nHPg2Wet2+vB/m12Lib8X8sVeRs22Cjqt9+GGTOsSun66639oX79ROfOueLPA4UrkrKyrHF60CD45BMrPdSrBy++CN27wxFHJDqHzpUcHihckaFqa1APGQJDh0JGho15+Mc/oEcPaNbM16Z2LhE8ULiEysqy6bw/+cRuv/5qbQ0dO1pwuPhib5x2LtE8ULi4274dvvrKAsOnn8Lq1VCmDJx3Htx9N3TpAscck+hcOudyeKBwMZedbfMtffutBYgvvoBt26yd4eKLrcfSRRfZKGrnXNHjgcJF3a5d8OOPMHmyBYcpU2DdOnuuWjWbjK9LF5vFtUyZhGbVOVcAMQ8UIlIdGAwcB2QDA1X1hbA0ArwAdAS2AT1VdbaINAFeBQ4HsoDHVPWDWOfZFc7atTB3rgWEyZOtzWHrVnvulFMsKLRpA23bwokneoO0c8kmHiWK3cAdwYW/AjBLRMar6oKQNB2AWsGtBRYcWmBBo7uq/iIixwevHaeqG+KQbxdm3TpITf3ztmCB3a9cac+LQMOGcO21FhjatIEqVRKbZ+fcgYt5oFDVDCAj+HuziCwEqgKhgaIzMFhVFZgmIhVFpIqqLgrZT7qIrAKOBjxQRJkqbNpkC/osX/7n/fLl8NNPFhBWrPgzffnyNq6hY0e7b9DAVo2rVClxx+Cci424tlGISE2gKTA97KmqwLKQx2nBtoyQ154BlAF+jbDfXkAvgBo1akQzy0lB1doFtm/f87Zt256PN260Uc45t5Ur9wwMOdVFoY44AmrXtsbmevVsJHT9+lC9Ohzky145VyLELVCISHlgBNBXVTeFPx3hJRry2irAEKCHqmbvlVB1IDAQoHnz5hr+fEGsWwdnnWUXXdvnnn9H2nagz0fjNVlZFgSy9zor+StfHo4+2mZbbdoUOnWC44+3xzn3Var4DKzOuTgFChEpjQWJd1V1ZIQkaUD1kMfVgPTgtYcDY4AHVHVarPJ48MFWfWLv+ed96N+Rth3o8wf6moMOgnLl9rwdeuje28qVs9JBxYp27xPoOecKKh69ngR4C1ioqs/mkWw0cIuIDMMasTeqaoaIlAE+xtovPoplPg8/HD78MJbv4JxzySkevytbA9cAKSIyJ9jWD6gBoKqvAWOxrrGLsZ5O1wbprgTaApVFpGewraeq5uzHOedcjMWj19MUIrdBhKZRoE+E7UOBoTHKmnPOuQLwfivOOefy5YHCOedcvjxQOOecy5cHCuecc/nyQOGccy5fHiicc87lS1T3a8aLIktEVgO/H8AujgLWRCk7ycKPuWTwYy4Z9veYT1DVoyM9UewCxYESkZmq2jzR+YgnP+aSwY+5ZIjFMXvVk3POuXx5oHDOOZcvDxR7G5joDCSAH3PJ4MdcMkT9mL2NwjnnXL68ROGccy5fHiicc87lywNFQEQuEpGfRWSxiNyb6PxEi4i8LSKrRGR+yLYjRWS8iPwS3FcKtouIvBicg3kiclricr7/RKS6iEwQkYUikioitwXbi+1xi0hZEZkhInODY/5XsP1EEZkeHPMHwWJgiMghwePFwfM1E5n/AyEipUTkRxH5LHhcEo55qYikiMgcEZkZbIvZ99sDBfZFA14GOgD1gL+JSL3E5ipqBgEXhW27F/haVWsBXwePwY6/VnDrBbwapzxG227gDlWtC7QE+gSfZ3E+7h3AuaraGGgCXCQiLYGngOeCY14PXB+kvx5Yr6qnAM8F6ZLVbcDCkMcl4ZgBzlHVJiFjJmL3/VbVEn8DzgTGhTy+D7gv0fmK4vHVBOaHPP4ZqBL8XQX4Ofj7deBvkdIl8w0YBVxQUo4bOBSYjS0r3tvZlAAAA2FJREFUvAY4ONie+z0HxgFnBn8fHKSTROd9P461WnBRPBf4DFskrVgfc5D/pcBRYdti9v32EoWpCiwLeZwWbCuujlXVDIDg/phge7E7D0H1QlNgOsX8uIMqmDnAKmA88CuwQVV3B0lCjyv3mIPnNwKV45vjqHgeuBvIDh5XpvgfM4ACX4rILBHpFWyL2fc7HmtmJ4NIS7WWxH7Dxeo8iEh5YATQV1U3ieS5Im+xOG5VzQKaiEhF4GOgbqRkwX3SH7OIdAJWqeosEWmXszlC0mJzzCFaq2q6iBwDjBeRn/JJe8DH7SUKkwZUD3lcDUhPUF7iYaWIVAEI7lcF24vNeRCR0liQeFdVRwabi/1xA6jqBmAi1j5TUURyfhCGHlfuMQfPHwGsi29OD1hr4BIRWQoMw6qfnqd4HzMAqpoe3K/CfhScQQy/3x4ozA9AraC3RBmgGzA6wXmKpdFAj+DvHlgdfs727kEviZbAxpyibDIRKzq8BSxU1WdDniq2xy0iRwclCUSkHHA+1sA7Abg8SBZ+zDnn4nLgGw0qsJOFqt6nqtVUtSb2P/uNql5FMT5mABE5TEQq5PwNtAfmE8vvd6IbZYrKDegILMLqde9PdH6ieFzvAxnALuyXxfVYvezXwC/B/ZFBWsF6f/0KpADNE53//Tzms7Ci9TxgTnDrWJyPG2gE/Bgc83zgoWD7ScAMYDHwEXBIsL1s8Hhx8PxJiT6GAzz+dsBnJeGYg+ObG9xSc65Xsfx++xQezjnn8uVVT8455/LlgcI551y+PFA455zLlwcK55xz+fJA4ZxzLl8eKJxzzuXLA4Vzzrl8eaBwLsZEpLeIqIicHbLtlmDb+YnMm3MF4YHCudhrhI2YrgsgIodiI+RXYyNlnSvSPFA4F3sNsalUTg0e34pNJZGtqisTlivnCsgDhXOxVxf4EDhVRI4AugLfY3MyOVfkeaBwLoZEpDqwVlWXYAvJ3A38F6iNVUc5V+R5oHAuthrxZzvEZmz98new6ihvn3BJwVe4cy62QgPCM1jpIktEGgKDEpYr5wrBA4VzsdUQW2kPVf0sZHs9YEFCcuRcIfl6FM455/LlbRTOOefy5YHCOedcvjxQOOecy5cHCuecc/nyQOGccy5fHiicc87lywOFc865fP0/zx7rgCKfxZsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plot the prices of the two options over time \n",
    "plt.plot(Chooser_EU, color='blue')\n",
    "plt.plot(Chooser_USA, color='red')\n",
    "plt.title('American vs European Chooser Option', fontsize=16)\n",
    "plt.xlabel('$M$')\n",
    "plt.ylabel('$Price$')\n",
    "\n",
    "#create legend with corresponding colours\n",
    "eu_patch = mpatches.Patch(color='blue', label='$European$ Chooser')\n",
    "us_patch = mpatches.Patch(color='red', label='$American$ Chooser')\n",
    "plt.legend(handles=[eu_patch, us_patch])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EDFtBcOcQBIK"
   },
   "source": [
    "To begin with, the American chooser option has a clear premium over the European one, due to the fact that the American one has the additional flexibility of being exercisable before the actual exercise date, as opposed to a European one. As time goes by, the holder of an American option has fewer and fewer opportunities to exercise the option early, hence both prices converge towards slightly over 2.08 as M approaches 500(=N--> exercise date)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1yZzGAqaQBIL"
   },
   "source": [
    "#  **Question B**:\n",
    "\n",
    " We wish to obtain the no-arbitrage price $X_{0}$ of an index-linked certificate. The derivative has payoff $X_{T}$, defined from the path of an equity index $\\{S_{t}\\}_{0\\leq t\\leq T}$ as follows:\n",
    " $$\n",
    " X_{T}=S_{0}\\prod\\nolimits_{t=1}^{T}(1+\\max (R_{t},R)),\\qquad R_{t}=\\frac{S_{t}-S_{t-1}}{S_{t-1}},\\qquad R=e^{r}-1,\n",
    " $$\n",
    " where $r\\equiv \\log (1+R)$ is the continuously compounded interest rate, so that $R$ is the corresponding annually compounded interest rate. Time is measured in years, and $T$ is a positive integer, representing the number of years after which the payoff is received. The idea behind this derivative is that the investor participates in the upside potential of the stock market, but receives a guaranteed minimum return. The value $X_{0}$ of this product cannot be evaluated analytically. \n",
    "\n",
    " We assume a Black-Scholes economy and use the following parameter values: $S_0=12$, $T=10$, $R=0.025$, $\\sigma =0.3$, and $q =0.01$. \n",
    "     \n",
    "   1. Using `asianmc_vec` (and `bmsim_vec`) as a starting point, modify the `payoff` of the derivative accordingly. Call the function `certificatemc(S0, R, T, sigma, q, numsim=10000)`, which will return the Monte Carlo estimate of $X_{0}$, along with a 95% confidence interval. **Notice that both the interval and the CI should be stated explicitly in your answers, based on a random seed of 0**.\n",
    "   \n",
    "   2. Consider a second derivative with payoff $F_{T}=\\max (S_{0}(1+R)^{T},S_{T})$. Show (with pen and paper, or preferably inside the notebook using Markdown and $\\LaTeX$) that the no-arbitrage price $F_{0}$ is\n",
    " $$  \n",
    " F_{0}=S_{0}(1+e^{-q T}\\Phi (d_{1})-\\Phi (d_{2})),\n",
    " $$\n",
    " where\n",
    " $$\n",
    " d_{1}=\\frac{(-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}},\\qquad\n",
    " d_{2}=d_{1}-\\sigma \\sqrt{T}=\\frac{(-q -\\frac{1}{2}\\sigma ^{2})T}{\\sigma\n",
    " \\sqrt{T}}.\n",
    " $$\n",
    " Hint: The payoff $F_T$ can be replicated with a position in cash bonds and a European call option with a non-standard strike price; re-arrange $F_{T}=\\max (S_{0}(1+R)^{T},S_{T})$ to understand this.\n",
    " \n",
    " Building on `blackscholes(S0, K, T, r, sigma)` (see week 4), write a function `F0(S0, R, T, sigma, q)` that evaluates the price of the derivative. \n",
    "**The resulting value for the given parameter values should be given explicitly in your answer**.\n",
    "\n",
    "3.   Write a function `certificatecv(S0, R, T, sigma, q, numsim=10000)` that improves upon your code from the previous question by using $F_{T}$ as a control variate.\n",
    "**Report the improved estimate $\\hat{X}_{0}$ along with an apropriately constructed confidence interval, based on a random seed of 0.** Then, compare the CI with the one based on the crude MC estimate, and discuss your findings; are they in line with the theory? Has the computational cost increased?\n",
    "\n",
    "4. Test if the confidence interval returned by `certificatecv` has correct coverage, by simulating 1000 independent ones (with `numsim=1000`), and recording in a hit series $I_j$, $j\\in\\{1, \\ldots, 1000\\}$, whether an interval contains the true price ($I_j=1$) or not ($I_j=0$). The test can be conducted by regressing $\\{I_j-0.95\\}_{j=1}^{1000}$ on an intercept and testing whether that is zero. The 'true' price should be obtained from `certificatecv` with `numsim=100000`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DhiVFzNuQBIL"
   },
   "source": [
    "Question B1. Using asianmc_vec (and bmsim_vec) as a starting point, modify the payoff of the derivative accordingly. Call the function certificatemc(S0, R, T, sigma, q, numsim=10000), which will return the Monte Carlo estimate of 𝑋0, along with a 95% confidence interval. Notice that both the interval and the CI should be stated explicitly in your answers, based on a random seed of 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "dvDSIzQSQBIM"
   },
   "outputs": [],
   "source": [
    "S0=12; T=10; R=0.025; sigma=0.3; q=0.01; N=T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "EpLJ3slUQBIP"
   },
   "outputs": [],
   "source": [
    "## Note new input: numsim --  the number of paths.\n",
    "def bmsim_vec(T, N, X0=0, mu=0, sigma=1, numsim=1): \n",
    "    \n",
    "#Simulate `numsim` quantity of Brownian motion paths.\n",
    "    deltaT = float(T)/N\n",
    "    tvec = np.linspace(0, T, N+1)\n",
    "    z = np.random.randn(numsim, N+1) #(N+1)->(numsim, N+1)\n",
    "    dX = mu*deltaT + sigma*np.sqrt(deltaT)*z\n",
    "    \n",
    "    dX[:, 0] = 0. #dX[0]->dX[:, 0]\n",
    "    X = np.cumsum(dX, axis=1) #cumsum(dX)->cumsum(dX, axis=1)\n",
    "    X += X0\n",
    "    return tvec, X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WuJKaZ02QBIS"
   },
   "outputs": [],
   "source": [
    "def certificatemc(S0, R, T, sigma, q, numsim=10000):\n",
    "    \n",
    "    #Monte Carlo price of an arithmetic average Asian call, adapted to index-linked certificate.\n",
    "    X0 = np.log(S0)\n",
    "    r = np.log(1+R) # define r as above\n",
    "    nu = r-q-.5*sigma**2\n",
    "    \n",
    "    #simulate all paths at once:\n",
    "    # simulate the stock price\n",
    "    _, X = bmsim_vec(T, N, X0, nu, sigma, numsim)\n",
    "    S = np.exp(X)\n",
    "    S_tlag = S[:,:N]\n",
    "    S_t = S[:,1:]\n",
    "    \n",
    "    #calculate return of index\n",
    "    Rt = (S_t - S_tlag) / S_tlag \n",
    "    \n",
    "    # take maximum of index return or guaranteed rate\n",
    "    R_max = np.maximum(Rt,R) \n",
    "    \n",
    "    #calculate return for each year\n",
    "    R_max = np.cumprod(1+R_max,axis=1) \n",
    "    \n",
    "    ## payyoffs\n",
    "    payoffs = S0 * R_max[:,N-1] \n",
    "    \n",
    "    ## discount the payoffs\n",
    "    g = np.exp(-r*T)*payoffs \n",
    "    \n",
    "    ## get mean and standard_deviation for the confidence intervals\n",
    "    C = g.mean(); s = g.std()\n",
    "    zq = norm.ppf(0.975)\n",
    "    Cl = C - zq/np.sqrt(numsim)*s\n",
    "    Cu = C + zq/np.sqrt(numsim)*s\n",
    "    return C, Cl, Cu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "_SuaQclBQBIW",
    "outputId": "54f04c86-16f7-47cc-e457-660badad4558"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35.232976937666535, 34.80365527102579, 35.66229860430728)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "certificatemc(S0, R, T, sigma, q, numsim=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "wgAEeqy1QBIZ",
    "outputId": "bc510b80-c7ff-4894-912d-88d6e290c334"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.39 ms ± 236 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit certificatemc(S0, R, T, sigma, q, numsim=10000) #ipython magic for timing things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "jQIv5PcjQBId",
    "outputId": "4637a438-5e53-4cba-9ea5-7bca6ee1f694"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The Monte Carlo estimation of X0 is 35.233, and the estimation of the 95% confidence interval is [34.804, 35.662].'"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "C, Cl, Cu = certificatemc(S0, R, T, sigma, q, numsim=10000)\n",
    "'The Monte Carlo estimation of X0 is %4.3f' %C+', and the estimation of the 95%% confidence interval is [%4.3f'%Cl+', %4.3f'%Cu+'].'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zpritCgSQBIh"
   },
   "source": [
    "Question B2. Consider a second derivative with payoff $F_{T}=\\max (S_{0}(1+R)^{T},S_{T})$. Show (with pen and paper, or preferably inside the notebook using Markdown and $\\LaTeX$) that the no-arbitrage price $F_{0}$ is\n",
    " $$  \n",
    " F_{0}=S_{0}(1+e^{-q T}\\Phi (d_{1})-\\Phi (d_{2})),\n",
    " $$\n",
    " where\n",
    " $$\n",
    " d_{1}=\\frac{(-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}},\\qquad\n",
    " d_{2}=d_{1}-\\sigma \\sqrt{T}=\\frac{(-q -\\frac{1}{2}\\sigma ^{2})T}{\\sigma\n",
    " \\sqrt{T}}.\n",
    " $$\n",
    " Hint: The payoff $F_T$ can be replicated with a position in cash bonds and a European call option with a non-standard strike price; re-arrange $F_{T}=\\max (S_{0}(1+R)^{T},S_{T})$ to understand this.\n",
    " \n",
    " Building on `blackscholes(S0, K, T, r, sigma)` (see week 4), write a function `F0(S0, R, T, sigma, q)` that evaluates the price of the derivative. \n",
    "**The resulting value for the given parameter values should be given explicitly in your answer**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4b6MrzJgQBIi"
   },
   "source": [
    "Solution:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pa6eMkKIQBIj"
   },
   "source": [
    "$F_{T}$can be seen as a position in a cash bond plus a European call option, and $S_{0}(1+R)^{T}$ can be seen as $K$ in the Black-Scholes model, hence we can replace $K$ by $S_{0}(1+R)^{T}$. This leads to:\n",
    "$$\n",
    "F_{T}=\\max (S_{0}(1+R)^{T},S_{T}) \n",
    "$$\n",
    "being rearranged as:\n",
    "$$\n",
    "F_{T}=S_{0}(1+R)^T+\\max(0,S_{T}-S_{0}(1+R^T)\n",
    "$$\n",
    "$F_{0}$ then becomes:\n",
    "\n",
    "$$\n",
    "F_{0}=S_{0}+S_{0}e^{-q T}\\Phi (d_{1})-S_{0}(1+R)^{T}e^{-rT}\\Phi (d_{2})\n",
    "$$\n",
    "since $(1+R)^{T}e^{-rT}=1$:\n",
    "$$\n",
    "F_{0}=S_{0}+S_{0}e^{-q T}\\Phi (d_{1})-S_{0}\\Phi (d_{2})\n",
    "$$ \n",
    "rearranging we get:\n",
    "<br>\n",
    "$$\n",
    "F_{0}=S_{0}(1+e^{-q T}\\Phi (d_{1})-\\Phi (d_{2}))\n",
    "$$\n",
    "<br>\n",
    "We also replace $K$ by $S_{0}(1+R)^{T}$ in $d_{1}$ and $d_{2}$. $d_{1}$ goes from:\n",
    "$$\n",
    "d_{1}=\\frac{log(\\frac{S_{0}}{K})+(r-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}}\\qquad\n",
    "$$\n",
    "To:\n",
    "$$\n",
    "d_{1}=\\frac{log(\\frac{S_{0}}{S_{0}(1+R)^{T}})+(r-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}}\\qquad\n",
    "$$\n",
    "<br>\n",
    "since $log(\\frac{S_{0}}{S_{0}(1+R)^{T}})=log(\\frac{1}{(1+R)^{T}})=-log(1+R)^T=-rT$:\n",
    "<br>\n",
    "$$\n",
    "d_{1}=\\frac{-rT+(r-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}}\\qquad\n",
    "$$\n",
    "Rearranging:\n",
    "<br>\n",
    "$$\n",
    "d_{1}=\\frac{(-q +\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}},\\qquad\n",
    "$$\n",
    "Hence, $d_{2}$ becomes:\n",
    "$$\n",
    "d_{2}=d_{1}-\\sigma \\sqrt{T}=\\frac{(-q -\\frac{1}{2}\\sigma ^{2})T}{\\sigma \\sqrt{T}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "c0TFriUWQBIk"
   },
   "outputs": [],
   "source": [
    "#F0 function for derivative based on blackscholes function\n",
    "def F0(S0, R, T, sigma, q): \n",
    "    d1 = (-q+0.5*(sigma**2))*T/(sigma*np.sqrt(T))    \n",
    "    d2 = d1-sigma*np.sqrt(T) \n",
    "    F0 = S0*(1+np.exp(-q*T)*norm.cdf(d1)-norm.cdf(d2))\n",
    "    return F0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "RxYabzUqQBIn",
    "outputId": "ee0347c0-ffc4-4b8a-b43d-9e2dcd19fa44"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.619120903425273"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#F0(12, 0.025, 10, 0.3, 0.01)\n",
    "F0(S0, R, T, sigma, q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 51
    },
    "colab_type": "code",
    "id": "vBhQZD1PQBIp",
    "outputId": "012b151d-55e0-4a2c-80e3-718738b28a2d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "201 µs ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit F0(S0, R, T, sigma, q) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "GWSliGNZQBIs",
    "outputId": "d633dd54-a308-46b1-afd7-975310fd743d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The price of the derivative based on the Black-Scholes model is 15.619'"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F01=F0(S0, R, T, sigma, q)\n",
    "'The price of the derivative based on the Black-Scholes model is %4.3f'%F01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3RmCxWgqQBIu"
   },
   "source": [
    "Question B3. Write a function `certificatecv(S0, R, T, sigma, q, numsim=10000)` that improves upon your code from the previous question by using $F_{T}$ as a control variate.\n",
    "**Report the improved estimate $\\hat{X}_{0}$ along with an appropriately constructed confidence interval, based on a random seed of 0.** Then, compare the CI with the one based on the crude MC estimate, and discuss your findings; are they in line with the theory? Has the computational cost increased?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yS4pA5BCQBIv"
   },
   "outputs": [],
   "source": [
    "## Monte Carlo price of an arithmetic average Asian call using control variates.\n",
    "def certificatecv(S0, R, T, sigma, q, numsim=10000):\n",
    "    ## variables as defined above \n",
    "    X0 = np.log(S0)\n",
    "    r = np.log(1+R) \n",
    "    nu = r-q-.5*sigma**2\n",
    "    #simulate all paths at once:\n",
    "    \n",
    "    _, X = bmsim_vec(T, N, X0, nu, sigma, numsim)\n",
    "    S = np.exp(X)\n",
    "    S_1 = S[:,:N]\n",
    "    S_2 = S[:,1:]\n",
    "    Rt = (S_2 - S_1) / S_1 \n",
    "    R_max = np.maximum(Rt,R) \n",
    "    R_max = np.cumprod(1+R_max,axis=1)\n",
    "    mu = F0(S0, R, T, sigma, q)\n",
    "    payoffs = S0 * R_max[:,N-1]\n",
    "    \n",
    "    ## create control variate as per financial theory \n",
    "    control_Ft = np.maximum(S0*(1+R)**T, S[:,N])\n",
    "    c = (np.cov(payoffs, control)/np.var(control))[0, 1]\n",
    "    \n",
    "    \n",
    "    g = np.exp(-r*T)*payoffs - c*(np.exp(-r*T)*control-mu)\n",
    "    C = g.mean(); s = g.std()\n",
    "    zq = norm.ppf(0.975)\n",
    "    Cl = C - zq/np.sqrt(numsim)*s\n",
    "    Cu = C + zq/np.sqrt(numsim)*s\n",
    "    return C, Cl, Cu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "colab_type": "code",
    "id": "3pKtRqP1QBIx",
    "outputId": "8bdfc13c-fa75-498e-c206-a5bec7ce1ba6"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35.136318890302036,\n",
       " 0.5325508901919989,\n",
       " 34.870043445206036,\n",
       " 35.402594335398035)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "C, Cl, Cu = certificatecv(S0, R, T, sigma, q, numsim=10000)\n",
    "C, Cu - Cl, Cl, Cu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WNZaiPnWQBI2"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "f8lcQx6xQBI5"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Vkpnc1bQQBI8"
   },
   "source": [
    "Question B4. Test if the confidence interval returned by `certificatecv` has correct coverage, by simulating 1000 independent ones (with `numsim=1000`), and recording in a hit series $I_j$, $j\\in\\{1, \\ldots, 1000\\}$, whether an interval contains the true price ($I_j=1$) or not ($I_j=0$). The test can be conducted by regressing $\\{I_j-0.95\\}_{j=1}^{1000}$ on an intercept and testing whether that is zero. The 'true' price should be obtained from `certificatecv` with `numsim=100000`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "k0CjF5LyQBI9"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aLFpaYz_QBJA"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Assignment 2.ipynb",
   "provenance": []
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}