ipynb/john_stachurski-euro_option.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "id": "8ed62ea1",
            "metadata": {
                "id": "8ed62ea1"
            },
            "source": [
                "## Pricing a European Call Option: CPU vs GPU and CUDA vs JAX."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "db3b91a2",
            "metadata": {
                "id": "db3b91a2"
            },
            "source": [
                "#### Author: [John Stachurski](http://johnstachurski.net/)\n",
                "#### Date: Nov 2022"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "820988bd",
            "metadata": {},
            "source": [
                "In this notebook we implement a simple task on the CPU and the GPU.  The task is estimating the risk-neutral price of a call option.  All our code is accelerated by a JIT-compiler, meaning that it runs fast on a single thread.  The main aim is to parallelize the work across multiple threads in the most efficient way.  We start on the CPU, which only has a few threads, and then shift to the GPU.\n",
                "\n",
                "On the GPU we try a traditional implmentation using Numba-CUDA and a JAX-based implementation.  We will see that the JAX-based implementation is must faster than all others.\n",
                "\n",
                "While the Numba-CUDA implementation could be improved by thinking more carefully about the underlying hardware, JAX takes care of efficient parallization automatically."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "f02daa8f",
            "metadata": {},
            "source": [
                "If you want to run this lecture and don't have a configured GPU, you can try it on [Google Colab](https://colab.research.google.com/).\n",
                "\n",
                "If you have access to your head of department's credit card, I recommend Colab Pro (faster GPUs)."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "af728dbd",
            "metadata": {},
            "source": [
                "We will use the following imports:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "id": "2447bf5f",
            "metadata": {
                "id": "2447bf5f"
            },
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "import matplotlib.pyplot as plt\n",
                "from collections import namedtuple"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "d16fe980",
            "metadata": {
                "id": "d16fe980"
            },
            "source": [
                "We are going to price a **European call option** under the assumption of risk neutrality.  \n",
                "\n",
                "Recall that the risk-neutral price of an asset that pays $G_n$ in $n$ periods (one payoff $n$ periods in the future) is \n",
                "\n",
                "$$ P = \\beta^n \\mathbb E G_n $$\n",
                "\n",
                "where $\\beta$ is a discount factor.\n",
                "\n",
                "For a European call option on with expiry at time $n$ and strike price $K$, the payoff is \n",
                "\n",
                "$$ G_n = \\max\\{S_n - K, 0\\} $$\n",
                "\n",
                "where $S_n$ is the price of the stock at time $n$.\n",
                "\n",
                "Thus, the risk-neutral price is\n",
                "\n",
                "$$ P = \\beta^n \\mathbb E \\max\\{ S_n - K, 0 \\} $$\n",
                "\n",
                "(For example, if the call option is to buy stock in Amazon at strike price $K$, the owner has the right to buy 1 share in Amazon at price $K$ after $n$ days.  The price is the expectation of the return $\\max\\{S_n - K, 0\\}$, discounted to current value.)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "5236b451",
            "metadata": {
                "id": "5236b451"
            },
            "source": [
                "## Lognormal Example\n",
                "\n",
                "Let's start with an easy case where $S_n$ has a known distribution.\n",
                "\n",
                "Suppose in particular that $S_n$ has the [log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution) distribution with parameters $\\mu$ and $\\sigma$.  Let $f$ denote the density of this distribution.  \n",
                "\n",
                "Then\n",
                "\n",
                "$$ P = \\beta^n \\int_0^\\infty g(x) dx $$\n",
                "\n",
                "where\n",
                "\n",
                "$$g(x) := \\max\\{x - K, 0\\} f(x)$$ \n",
                "\n",
                "Let's plot $g$ over the interval $[0, 400]$ when"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "id": "c7003806",
            "metadata": {
                "id": "c7003806"
            },
            "outputs": [],
            "source": [
                "μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "id": "8e909043",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/",
                    "height": 265
                },
                "id": "8e909043",
                "outputId": "1851ede9-9162-4ad4-f57f-8b1ba14baddd"
            },
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAf7klEQVR4nO3de3Bc533e8e8Pi8viDhAESRAACZKiJVGySEnULbblOGPXlOyYUZKJZbu+tPao8lhJPZ1MK08q1246aRy7vrRWwmFtNU6VVGkbxWZsurLj+B5HJiVRNCmZF/EiglcQ1wWBBXaBt3/sWXAFLoC9nN2zu3g+Mxjsnj0456cz0MMX73nP+5pzDhERKX9VQRcgIiL+UKCLiFQIBbqISIVQoIuIVAgFuohIhagO6sQrV650fX19QZ1eRKQsPffcc5edc53pPgss0Pv6+ti/f39QpxcRKUtmdnqhz9TlIiJSIRToIiIVQoEuIlIhAutDFxEptlgsRn9/P9FoNOhSlhQOh+np6aGmpibjn8ko0M1sB/AlIAR8xTn3x/M+/1XgG8BJb9PTzrn/mHEVIiJF0N/fT3NzM319fZhZ0OUsyDnH4OAg/f39bNiwIeOfWzLQzSwEPA68DegH9pnZHufcS/N2/bFz7p3ZFC0iUkzRaLTkwxzAzOjo6GBgYCCrn8ukD/1O4Lhz7oRzbhp4CtiZQ40iIoEr9TBPyqXOTAK9GziT8r7f2zbfPWb2opl928xuWqDAh8xsv5ntz/ZfnkrwvZcv8pfPniY2Mxt0KSJSgTLpQ0/3z8T8SdSfB9Y758bN7H7g68Dma37Iud3AboDt27cvq4nYf3DkEh/+WuJBqleHJvjEfTcGXJGIVJpMWuj9QG/K+x7gXOoOzrkx59y493ovUGNmK32rsgJ84e+PsXFlI+/aupav/Pgkl8ZK/y67iJSXTAJ9H7DZzDaYWS3wILAndQczW2Neh4+Z3ekdd9DvYsvV8UvjvHhmhPfetY6Pv3UzM7OOp184G3RZIhKQl19+mXvvvZdbbrmFz372s1x33XW+HHfJLhfnXNzMHgGeITFs8Qnn3GEze9j7fBfw28BHzSwOTAIPOq1tN+fHxxL3C3bcvIae9ga29rTyncMXePjNmwKuTGT5+vTfHealc2O+HnPL2hb+w6+nvYU4Jx6P8773vY+vfvWr3HrrrXz0ox/l5ptv9uX8GY1D97pR9s7btivl9ZeBL/tSUQX62SuD9K6op6e9AYA3X7+KL//DMUYmpmlrqA24OhEppqeffpqtW7dy6623ArBlyxZWrVrly7H1pGgRPP/qCPe+7uothTdtXsl//d4x9p0a5m1bVgdYmcjytVRLulAOHjzItm3b5t4fOnSIHTt2+HJszeVSYIPjU1wen2JLV8vcttd3t1JdZRw4MxxgZSIShI6ODo4ePQrAgQMHePLJJ9m6dasvx1YLvcCOXIgAcP2a5rlt4ZoQN3Q1c+DMSEBViUhQ3v/+9/OOd7yDO+64g3vuuYe+vj42btzoy7HVQi+wX6YJdIBtvW0cPDPK7KzuHYssJ+FwmGeffZZ9+/bR29vLAw884NuxFegFduRChI7GWjqb6l6zfWtPG5GpOCcuXwmoMhEJwhe+8AVuuukmtm3bxqlTp3jsscd8O7a6XArs2KUIr1vdfM28DDd6fepHL0a4blVTEKWJSAAee+wxX0M8lVroBfbq0CTrOxqu2X7dqibMrvaxi4jkS4FeQJPTM1wen6J3xbWBHq4J0dfRyNGLCnQR8YcCvYD6hycA6GmvT/v561Y3cUSBLlJU5fIQey51KtAL6IwX6Ola6ADXr27m1OUrRGMzxSxLZNkKh8MMDg6WfKgnVywKh8NZ/ZxuihbQq4NeoLenD/RNq5qYdXBmaILNq5vT7iMi/unp6aG/vz/rlYCCkFxTNBsK9AI6MzxJfU2IlU3p52vp62gE4OTlKwp0kSKoqanJao3OcqMulwLqH56gu71+waWk+lYmAv3UoMaii0j+FOgFdHFsiq7WhfvAWutrWNFYy8nLE0WsSkQqlQK9gC6ORVndsvhNjb6OBk7paVER8YECvUBmZx2XIlOsbqlbdL++lY3qchERXyjQC2TwyjQzs27JFvqGjkbOj0aZnNbQRRHJjwK9QC56i0Cval6iy8W7MXp6SK10EcmPAr1AkoG+ZJeLN3TxlG6MikieFOgFcnFsCoA1i4xygavTApwdmSx4TSJS2RToBXJhLIoZrGxavIXe1lBDU1313LwvIiK5UqAXyEAkSkdjLTWhxS+xmdHTXs+ZIbXQRSQ/CvQCuTw+TUfj4q3zpJ72erXQRSRvCvQCGboyzYrG9HO4zNfT3sDZ4cmSnwFOREqbAr1Ahq5M07HApFzz9bTXE5mKMzYZL3BVIlLJFOgFcnl8io4sWuhwdf50EZFcKNALYDo+SyQaZ0UWfeiA+tFFJC8K9AIYnpgGyLjLJbkARv+wRrqISO4U6AUwOO4FeoZdLi311TTXVSvQRSQvCvQCGLySeEo001EuZka3hi6KSJ4U6AUwdCXZ5ZJZHzokboyqhS4i+cgo0M1sh5kdMbPjZvboIvvdYWYzZvbb/pVYfrLtcgFY2xbmnOZzEZE8LBnoZhYCHgfuA7YA7zGzLQvs9xngGb+LLDeDV6YIVRmt9TUZ/0xXaz1j0ThXpjQWXURyk0kL/U7guHPuhHNuGngK2Jlmv98F/ga45GN9ZWnoyjTtDbVUVaVfHDqdtW2JWRnPj6qVLiK5ySTQu4EzKe/7vW1zzKwbeADYtdiBzOwhM9tvZvsHBgayrbVsDI5PZ9XdAokWOsC5kWghShKRZSCTQE/XzJw/6cgXgX/nnFt0HTXn3G7n3Hbn3PbOzs4MSyw/g1nM45LU1aoWuojkpzqDffqB3pT3PcC5eftsB54yM4CVwP1mFnfOfd2PIsvN8MQ0N65pyepn1rSGMVMLXURyl0mg7wM2m9kG4CzwIPDe1B2ccxuSr83sz4FvLtcwBxidiNHakPkNUYCaUBWdTXVqoYtIzpYMdOdc3MweITF6JQQ84Zw7bGYPe58v2m++3DjnGJ2M0ZbFCJekrrZ6zo+qhS4iucmkhY5zbi+wd962tEHunPtQ/mWVryvTM8RnXVZDFpPWtoY5ejFSgKpEZDnQk6I+G52MAYm1QrPV1ZpooWuhCxHJhQLdZyPeTIs5tdDbwkxMz2ihCxHJiQLdZ8kWemt9dsMWIWUsum6MikgOFOg+G53Io8tFT4uKSB4U6D4bmWuh53JTVE+LikjuFOg+y+emaGdzHaEqUwtdRHKiQPfZyESMmpBRXxPK+mdDVcbq5jrOq4UuIjlQoPtsdDJGa30t3jQIWetqq9dNURHJiQLdZ6OT0zl1tyR1tYb1tKiI5ESB7rNECz33QF/bpoeLRCQ3CnSfjUzkNo9LUldrmOn4LIPeuqQiIplSoPtsZCK/Fnry4SLdGBWRbCnQfTY2mf3UuamSS9HpxqiIZEuB7qP4zCyRqThtOTz2n3S1ha5AF5HsKNB9NBZNTKrVWp/RrMRpdTTWUhuq0kgXEcmaAt1HyZkW2xpyb6FXVRlrWsOcU6CLSJYU6D4azWMel1RdrWEuqA9dRLKkQPfR3MRcedwUhUSga4IuEcmWAt1HY16gt4TzDPS2ei6ORZmZ1cNFIpI5BbqPkjdFW/K4KQqJtUXjs47L41N+lCUiy4QC3UeRqE8t9Ll50dWPLiKZU6D7aGwyTk3IqKvO77JeXblI/egikjkFuo8i0RjN4Zqcp85NSq5cpEAXkWwo0H0UicZpDufXfw6J1Y7qqqv0tKiIZEWB7qNINJZ3/zmAmc1NoysikikFuo/8aqGDNxZdDxeJSBYU6D7yN9DrNYWuiGRFge6j5E1RP6xtC3MpEiU+M+vL8USk8inQfeRnC31Na5hZB5cierhIRDKjQPfJ7KxjfDruXwt9buii+tFFJDMKdJ+MT8dxDlr86kNPrlykfnQRyVBGgW5mO8zsiJkdN7NH03y+08wOmtkBM9tvZm/0v9TSlpyYy8+boqAWuohkbsn0MbMQ8DjwNqAf2Gdme5xzL6Xs9j1gj3POmdktwP8GbihEwaUq4k3M5VeXS0u4msbakFroIpKxTFrodwLHnXMnnHPTwFPAztQdnHPjzrnkXK+NwLKb9zUZ6H48WASJh4u62uq5oIeLRCRDmQR6N3Am5X2/t+01zOwBM/sl8C3gX6Y7kJk95HXJ7B8YGMil3pKVnGnRry4XSDxcpC4XEclUJoGebqapa1rgzrm/dc7dAPwG8IfpDuSc2+2c2+6c297Z2ZlVoaXuapeLv4GutUVFJFOZBHo/0Jvyvgc4t9DOzrkfAZvMbGWetZWVqy10f7pcIHFj9PL4FNNxPVwkIkvLJND3AZvNbIOZ1QIPAntSdzCz68ybM9bMbgNqgUG/iy1lYwVooa9tC+McXBxTK11ElrZk+jjn4mb2CPAMEAKecM4dNrOHvc93Ab8FfMDMYsAk8O6Um6TLQiQapzZURbgm5NsxU1cu6l3R4NtxRaQyZdScdM7tBfbO27Yr5fVngM/4W1p5Sczj4l/rHBItdIALaqGLSAb0pKhPxnycxyVpzVwLXYEuIktToPvEz5kWk5rqqmkOV2vooohkRIHuk0g0Tku9vy10SEzSpRa6iGRCge6TSDRGc52/LXRITNKlFrqIZEKB7hM/50JP1dWqtUVFJDMKdJ8kAr0ALfTWMENXponGZnw/tohUFgW6D2ZmHeNThWqhe0MX1UoXkSUo0H0wPuX/U6JJa9u8oYvqRxeRJSjQfZCcx8WvqXNTJVvo5zXSRUSWoED3wdhk4VroWrlIRDKlQPdBIWZaTKqvDdHeUKNpdEVkSQp0H8ytVlSAB4sgMQWAboqKyFIU6D6ITBWuhQ7Q3Rbm3Ii6XERkcQp0HxRitaJUPe0N9A9PssxmJBaRLCnQfVD4QK9nfCrO6GSsIMcXkcqgQPdBJBqntrqKumr/FrdI1dOeGOnSP6xuFxFZmALdB4mJuQrTOodElwtA//BEwc4hIuVPge6DQj32n9Q7F+hqoYvIwhToPohE4zQVMNBb6qtprqtWoIvIohToPhiPxgsyF3qSmdHdXq8uFxFZlALdB2PRWEFb6JDoRz8zpBa6iCxMge6DQi1ukarHa6FrLLqILESB7oPxqXhBR7kA9K5o4Mr0DCMTGosuIukp0PPkXHJxi8L1oYPGoovI0hToeZqMzTAz64rQh54MdN0YFZH0FOh5KvRj/0k9GosuIktQoOcpGehNBe5Db62voTlcrRa6iCxIgZ6nQi4/N1+vN+uiiEg6CvQ8JReILnQfOiSHLirQRSQ9BXqeitWHDol+9FeHNBZdRNLLKNDNbIeZHTGz42b2aJrP32dmB72vfzSzrf6XWprGi9SHDtC3soHJ2AwD41MFP5eIlJ8lA93MQsDjwH3AFuA9ZrZl3m4ngTc7524B/hDY7XehpWqsgAtEz7e+oxGA04O6MSoi18qkhX4ncNw5d8I5Nw08BexM3cE594/OuWHv7T8BPf6WWbrm+tCL0ULvSAxdPHX5SsHPJSLlJ5NA7wbOpLzv97Yt5MPAt/MpqpxEonEaa0OEqqzg5+puq6e6ytRCF5G0MmlWpkuqtHflzOwtJAL9jQt8/hDwEMC6desyLLG0jUcL/9h/UnWoiu72ek4NqoUuItfKpIXeD/SmvO8Bzs3fycxuAb4C7HTODaY7kHNut3Nuu3Nue2dnZy71lpzIVOGnzk21vqNRLXQRSSuTQN8HbDazDWZWCzwI7EndwczWAU8D73fOHfW/zNJVjKlzU/V1NHBq8IqGLorINZZMIudc3MweAZ4BQsATzrnDZvaw9/ku4JNAB/CnZgYQd85tL1zZpaPYgb6+o5FINM7wRIwVjbVFO6+IlL6Mksg5txfYO2/brpTXHwE+4m9p5WF8Kk53W33Rzpcc6XJ68IoCXUReQ0+K5ikSjRVlyGKSxqKLyEIU6HkqdpdL74p6zNBIFxG5hgI9DzOzjonpmaKOcqmrDrG2tV4tdBG5hgI9D+NzE3MVZxx60vqOBk7oaVERmUeBnofIlDePSxH70AE2dTZx4tK4hi6KyGso0PNQzKlzU123qonIVJyBiGZdFJGrFOh5KObiFqk2dTYBcHxgvKjnFZHSpkDPQ6SIU+emum5VItBfuaRAF5GrFOh5KNYC0fOtbqmjqa6aVwZ0Y1RErlKg5yEZ6C1F7nIxMzZ1NnJcLXQRSaFAz0NQfeiQ6Ed/RX3oIpJCgZ6HSDRGqMqorwkV/dybVjVxfjQ694+KiIgCPQ/j3mP/3gyTRZUc6XJCrXQR8SjQ8xCJxot+QzQpOdJF/egikqRAz0NkqnjLz823vqOB6ipToIvIHAV6HiLRWNEf+0+qCVWxYWUjRy9GAjm/iJQeBXoexqeKO3XufDd2tfDyeQW6iCQo0PMQicYDGbKYdENXM2dHJhnznlgVkeVNgZ6HYi9uMd+NXS0A/FKtdBFBgZ6X8WicprpgbooC3LjGC/QLY4HVICKlQ4Geo2hshumZ2UBb6Ktb6mhrqOHl8wp0EVGg5yz5hGaQgW5m3LhGN0ZFJEGBnqOgFreY74auZo5ciDAzq9WLRJY7BXqOxuemzg2uDx0SN0YnYzO8OqRFo0WWOwV6jpJDBYs9de58yRuj6kcXEQV6jsYmE4He2hBsC33z6iaqq4xDZ0cDrUNEgqdAz9HoZLKFHmygh2tC3NDVzIv9I4HWISLBU6DnaK7LpT7YQAfY2tPGwf5RZnVjVGRZU6DnaGwyTqjKaKwt/uIW823taSMSjXNyUGuMiixnCvQcjU7GaAlocYv5tva2AXBQ3S4iy5oCPUdj0VhJdLdAYrGLhtoQL57RjVGR5SyjQDezHWZ2xMyOm9mjaT6/wcx+ZmZTZvb7/pdZesYmY4HfEE0KVRk3d7fqxqjIMrdkoJtZCHgcuA/YArzHzLbM220I+D3gc75XWKLGonFaS6SFDrCtt43D58aYis8EXYqIBCSTFvqdwHHn3Ann3DTwFLAzdQfn3CXn3D5g2UzMPToZo6U+2IeKUm1f3850fJaD/ep2EVmuMgn0buBMyvt+b1vWzOwhM9tvZvsHBgZyOUTJKKUuF4A7+lYA8POTQwFXIiJBySTQ0w3jyGnAs3Nut3Nuu3Nue2dnZy6HKBmldFMUoL2xls2rmhToIstYJoHeD/SmvO8BzhWmnPIwFZ8hGpstqT50gDs2rOD508OaeVFkmcok0PcBm81sg5nVAg8CewpbVmkbm0zMtBj0xFzz3bVhBZGpuCbqElmmlkwk51zczB4BngFCwBPOucNm9rD3+S4zWwPsB1qAWTP7OLDFOVeRyVJKj/2nSvajP3tyiJu7WwOuRkSKLaMmpnNuL7B33rZdKa8vkOiKWRaSMy2WWqCvbatnw8pGfnr8Mh9+44agyxGRItOTojkolZkW03nT5pX87JVBjUcXWYYU6DkY81Yrai2hcehJ927uZDI2w3Onh4MuRUSKTIGeg7ESbqHfvamD6irjR0cvB12KiBSZAj0HoyXahw7QVFfN7evb+dHR8n5wS0Syp0DPwVg0Rm11FeGa4OdCT+fN13fy0vkxzo9OBl2KiBSRAj0HY5PxkuxuSXr7TWsA+M7hiwFXIiLFpEDPwdhkrCRviCZt6mxi86om/t+hC0GXIiJFpEDPwVg0RnMJt9ABdty8hmdPDjJ0ZTroUkSkSBToORiemKa9obQD/e03rWHWwXdfUitdZLlQoOdg+EqM9sbaoMtY1E1rW+jraODrLyzredRElhUFeg4SLfTSDnQz47du6+FnJwY5MzQRdDkiUgQK9CxFYzNMTM+wosRb6AC/eXsPZvD082eDLkVEikCBnqWRicRDRW0l3ocO0N1Wz69s6uD/Pn+GWc2RLlLxFOhZGp5IjBop9S6XpN/Z3suZoUl+qCdHRSqeAj1L5Rbo97++i9UtdXz1JyeDLkVECkyBnqXhK4kul/bG0u9yAagJVfGBe/r4yfHLHLkQCbocESkgBXqWyq2FDvC+u9YRrqli1w9fCboUESkgBXqWhr0nL8vhpmhSW0MtH7ynj68fOMuxi2qli1QqBXqWhidiNNaGqKsuzZkWF/LwmzfRWFvN5797NOhSRKRAFOhZGpmYpq2MuluS2htr+cibNvDtQxf4+cmhoMsRkQJQoGdpaGK6LB4qSuehezfS3VbPv//6L5iOzwZdjoj4TIGepeGJWFn1n6dqqK3m0++6iaMXx9n9I90gFak0CvQsDV8p3xY6wFu3rOYdt3Txxb8/xvOvaiFpkUqiQM9SOUzMtZQ/euD1rGkN87t/9YLmSxepIAr0LERjM0SicTrKuIUO0Fpfw5ffexuXx6f4yNf2EY3NBF2SiPhAgZ6FQa8129lcF3Al+dvW28aXHtzGC2dGePjJ5xTqIhVAgZ6FgcgUUBmBDrDj5i7+6IHX88OjA3zof/yckQl1v4iUMwV6Fiot0AHec+c6vvjubTx3eph3/refcOjsaNAliUiOFOhZSAb6yqbKCXSAndu6+et/dQ8zs44H/vSnfO6ZI+qCESlDCvQsXB5PBHpHU3nfFE3ntnXtfOv33sSvb13Ll79/nLd87gf8+U9PKthFyogCPQsDkSnaGmrKbh6XTK1orOXzv7ONpx66m572ej71dy9x93/+Hp/8xiGef3WYGa16JFLSqjPZycx2AF8CQsBXnHN/PO9z8z6/H5gAPuSce97nWgM3EJmquO6WdO7e2MH/efhXePbEIE8++yp/ve8Mf/Gz07Q11PCG61Zy+7p2blrbwpa1LTSHy/OpWZFKtGSgm1kIeBx4G9AP7DOzPc65l1J2uw/Y7H3dBfyZ972iDIxP0bkMAj3pro0d3LWxg7FojO//8hI/PnaZnxy7zLcOnp/bZ1VzHd3t9XS3Jb5WNNbS1lBDa33ie1NdNeGaKuqqQ9TVVBGuCVFXXUVtqIpEO0BE/JJJC/1O4Lhz7gSAmT0F7ARSA30n8BfOOQf8k5m1mVmXc+78tYfLzw+PDvCfvvnS0jsWwLFL4/zmrd2BnDtILeEadm7rZue2xH/7pUiUw+fGeOncGKcHr3B2ZJJDZ0f5zuGLTM9kNumXGVRXGVWW+ApVGVUGVVVGyIwq733ytRkYds0xXvP+mnPYop+n25jTMUSy9O47evnImzb6ftxMAr0bOJPyvp9rW9/p9ukGXhPoZvYQ8BDAunXrsq0VgKa6ajavbsrpZ/P1+u5WPvZr1wVy7lKyqjnMquvDvOX6Va/Z7pxjYnqGkckYIxPTjE7EiEzFmYrPEo3NMBWfZcr7Ho3NEJ91zM46Zp1jZhZmXfK1Y9bB7Kxjxtvm5nXfu3kb5vfuX7N/mv+OpY4xf4NLexSR7BWq6zaTQE/XKJn/m53JPjjndgO7AbZv357T/x23r2/n9vW35/KjUmBmRmNdNY111XS31Qddjsiyk8kol36gN+V9D3Auh31ERKSAMgn0fcBmM9tgZrXAg8CeefvsAT5gCXcDo4XoPxcRkYUt2eXinIub2SPAMySGLT7hnDtsZg97n+8C9pIYsnicxLDFf1G4kkVEJJ2MxqE75/aSCO3UbbtSXjvgY/6WJiIi2dCToiIiFUKBLiJSIRToIiIVQoEuIlIhbP7TckU7sdkAcDrHH18JXPaxHL+Ual1QurWpruyoruxUYl3rnXOd6T4ILNDzYWb7nXPbg65jvlKtC0q3NtWVHdWVneVWl7pcREQqhAJdRKRClGug7w66gAWUal1QurWpruyoruwsq7rKsg9dRESuVa4tdBERmUeBLiJSIcou0M1sh5kdMbPjZvZowLWcMrNfmNkBM9vvbVthZt81s2Pe9/Yi1PGEmV0ys0Mp2xasw8w+4V2/I2b29iLX9SkzO+tdswNmdn8AdfWa2ffN7GUzO2xm/9rbHug1W6SuQK+ZmYXN7Odm9qJX16e97UFfr4XqCvx3zDtXyMxeMLNveu8Lf72cc2XzRWL63leAjUAt8CKwJcB6TgEr5237E+BR7/WjwGeKUMe9wG3AoaXqALZ4160O2OBdz1AR6/oU8Ptp9i1mXV3Abd7rZuCod/5Ar9kidQV6zUisSNbkva4BngXuLoHrtVBdgf+Oeef7N8BfAd/03hf8epVbC31uwWrn3DSQXLC6lOwEvua9/hrwG4U+oXPuR8BQhnXsBJ5yzk05506SmMP+ziLWtZBi1nXeOfe89zoCvExiDdxAr9kidS2kWHU559y497bG+3IEf70WqmshRfsdM7Me4B3AV+adv6DXq9wCfaHFqIPigO+Y2XOWWAAbYLXzVmvyvq9a8KcLa6E6SuEaPmJmB70umeSfnYHUZWZ9wK0kWnclc83m1QUBXzOv++AAcAn4rnOuJK7XAnVB8L9jXwT+LTCbsq3g16vcAj2jxaiL6A3OuduA+4CPmdm9AdaSqaCv4Z8Bm4BtwHngv3jbi16XmTUBfwN83Dk3ttiuabYVrLY0dQV+zZxzM865bSTWC77TzG5eZPeg6wr0epnZO4FLzrnnMv2RNNtyqqvcAr2kFqN2zp3zvl8C/pbEn0kXzawLwPt+KaDyFqoj0GvonLvo/U84C/x3rv5pWdS6zKyGRGj+pXPuaW9z4NcsXV2lcs28WkaAHwA7KIHrla6uErhebwDeZWanSHQL/5qZPUkRrle5BXomC1YXhZk1mllz8jXwz4BDXj0f9Hb7IPCNIOpbpI49wINmVmdmG4DNwM+LVVTyF9rzAIlrVtS6zMyArwIvO+c+n/JRoNdsobqCvmZm1mlmbd7reuCtwC8J/nqlrSvo6+Wc+4Rzrsc510cio/7BOffPKcb1KtQd3kJ9kViM+iiJO8F/EGAdG0ncmX4ROJysBegAvgcc876vKEIt/4vEn5YxEv/af3ixOoA/8K7fEeC+Itf1P4FfAAe9X+SuAOp6I4k/aQ8CB7yv+4O+ZovUFeg1A24BXvDOfwj45FK/6wHXFfjvWMr5fpWro1wKfr306L+ISIUoty4XERFZgAJdRKRCKNBFRCqEAl1EpEIo0EVEKoQCXUSkQijQRUQqxP8Hp3TC7g2q2ssAAAAASUVORK5CYII=\n",
                        "text/plain": [
                            "<Figure size 432x288 with 1 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "from scipy.stats import lognorm\n",
                "\n",
                "def f(x):\n",
                "    \"Lognormal density.\"\n",
                "    return lognorm.pdf(x, σ, scale=np.exp(μ))\n",
                "\n",
                "def g(x):\n",
                "    return np.maximum(x - K, 0) * f(x)\n",
                "\n",
                "x_grid = np.linspace(0, 400, 1_000)\n",
                "y_grid = g(x_grid) \n",
                "\n",
                "fig, ax = plt.subplots()\n",
                "ax.plot(x_grid, y_grid, label=\"$g$\")\n",
                "ax.legend()\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "57358ce7",
            "metadata": {
                "id": "57358ce7"
            },
            "source": [
                "### Integration with Quadrature"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "2dda4e74",
            "metadata": {
                "id": "2dda4e74"
            },
            "source": [
                "We can integrate numerically using quadrature to obtain the price:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "id": "04ca1571",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "04ca1571",
                "lines_to_next_cell": 2,
                "outputId": "d69eccb2-10a0-4596-a2ed-e15988155f50"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "The numerical integration based option price is 15.188\n"
                    ]
                }
            ],
            "source": [
                "from scipy.integrate import quad\n",
                "\n",
                "integral, error = quad(g, 0, 1_000)\n",
                "P = β**n * integral \n",
                "print(f\"The numerical integration based option price is {P:.3f}\")"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "01bc94d2",
            "metadata": {
                "id": "01bc94d2"
            },
            "source": [
                "### Integration via Monte Carlo"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "ef33d190",
            "metadata": {
                "id": "ef33d190"
            },
            "source": [
                "We can get a similar result using Monte Carlo to compute the expectation term in the option price, rather than `quad`."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "3f1bf35b",
            "metadata": {
                "id": "3f1bf35b"
            },
            "source": [
                "We use the fact that if $S_n^1, \\ldots, S_n^M$ are independent draws from the lognormal distribution specified above, then, by the law of large numbers,\n",
                "\n",
                "$$ \\mathbb E \\max\\{ S_n - K, 0 \\} \n",
                "    \\approx\n",
                "    \\frac{1}{M} \\sum_{m=1}^M \\max \\{S_n^m - K, 0 \\}\n",
                "    $$\n",
                "  "
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "id": "cabb4a3c",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "cabb4a3c",
                "outputId": "f9286a55-de80-432b-e64f-633d550dea4a"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "The Monte Carlo option price is 15.195\n"
                    ]
                }
            ],
            "source": [
                "M = 10_000_000\n",
                "S = np.exp(μ + σ * np.random.randn(M))\n",
                "payoffs = np.maximum(S - K, 0)\n",
                "P = β**n * payoffs.mean()\n",
                "print(f\"The Monte Carlo option price is {P:.3f}\")"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "574b8ecd",
            "metadata": {
                "id": "574b8ecd"
            },
            "source": [
                "## A Model with Stochastic Volatility\n",
                "\n",
                "In the last example we had a known distribution for $S_n$.\n",
                "\n",
                "Now let's investigate a more realistic model, where we begin by specifying the underlying dynamics."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "4d25b29d",
            "metadata": {
                "id": "d2904228"
            },
            "source": [
                "One common model for $\\{S_t\\}$ is\n",
                "\n",
                "$$ \\ln \\frac{S_{t+1}}{S_t} = \\mu + \\sigma \\xi_{t+1} $$\n",
                "\n",
                "where $\\{ \\xi_t \\}$ is IID and standard normal.  However, its predictions are counterfactual because volatility is not stationary but rather changes over time.  \n",
                "\n",
                "Here is an improved version:\n",
                "\n",
                "$$ \\ln \\frac{S_{t+1}}{S_t} = \\mu + \\sigma_t \\xi_{t+1} $$\n",
                "\n",
                "where \n",
                "\n",
                "$$ \n",
                "    \\sigma_t = \\exp(h_t), \n",
                "    \\quad\n",
                "        h_{t+1} = \\rho h_t + \\nu \\eta_{t+1}\n",
                "$$\n",
                "\n",
                "Here $\\{\\eta_t\\}$ is also IID and standard normal."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "d2904228",
            "metadata": {
                "id": "d2904228"
            },
            "source": [
                "Putting this together, with $s_t := \\ln S_t$, the dynamics become\n",
                "\n",
                "$$\n",
                "\\begin{aligned}\n",
                "s_{t+1} & = s_t + \\mu + \\exp(h_t) \\xi_{t+1}\\\\\n",
                "h_{t+1} & = \\rho h_t + \\nu \\eta_{t+1}\n",
                "\\end{aligned}\n",
                "$$"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "e1291080",
            "metadata": {},
            "source": [
                "### Default Specification"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "5e4b5ce5",
            "metadata": {},
            "source": [
                "We use a namedtuple for storing parameters and initial conditions"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 6,
            "id": "e447d71b",
            "metadata": {
                "id": "e447d71b",
                "lines_to_next_cell": 1
            },
            "outputs": [],
            "source": [
                "Model = namedtuple(\"Model\", \n",
                "                   (\"β\", \"μ\", \"K\", \"n\", \"ρ\", \"ν\", \"S0\", \"h0\"))"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "60dba19a",
            "metadata": {},
            "source": [
                "Here's a function to create an instance of `Model` with a set of default parameters."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "id": "ded4af65",
            "metadata": {
                "id": "ded4af65",
                "lines_to_next_cell": 1
            },
            "outputs": [],
            "source": [
                "def create_model(β=0.95,     # discount factor\n",
                "                 μ=0.0001,   # mean price growth rate\n",
                "                 K=100,      # strike price\n",
                "                 n=10,       # expiry\n",
                "                 ρ=0.1,      # volatility correlation\n",
                "                 ν=0.001,    # volatility of volatility\n",
                "                 S0=10,      # initial S\n",
                "                 h0=0):      # inital h\n",
                "\n",
                "    return Model(β=β, μ=μ, K=K, n=n, ρ=ρ, ν=ν, S0=S0, h0=h0)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "9ff267d1",
            "metadata": {},
            "source": [
                "Let's create a default instance to use below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "id": "75afe6de",
            "metadata": {
                "id": "75afe6de"
            },
            "outputs": [],
            "source": [
                "default_model = create_model()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "322139c0",
            "metadata": {
                "id": "322139c0"
            },
            "source": [
                "### Visualizing Paths\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "99f53167",
            "metadata": {
                "id": "99f53167"
            },
            "source": [
                "\n",
                "Here is a function to simulate a path for prices."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "id": "1695150b",
            "metadata": {
                "id": "1695150b"
            },
            "outputs": [],
            "source": [
                "from numpy.random import randn\n",
                "\n",
                "def simulate_asset_price_path(model):\n",
                "    \n",
                "    β, μ, K, n, ρ, ν, S0, h0 = model\n",
                "    s = np.empty(n+1)\n",
                "    s[0] = np.log(S0)\n",
                "\n",
                "    h = h0\n",
                "    for t in range(n):\n",
                "        s[t+1] = s[t] + μ + np.exp(h) * randn()\n",
                "        h = ρ * h + ν * randn()\n",
                "        \n",
                "    return np.exp(s)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "81c7da30",
            "metadata": {
                "id": "81c7da30"
            },
            "source": [
                "Here we plot some paths and the log of the paths."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 10,
            "id": "0ccd29c0",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/",
                    "height": 297
                },
                "id": "0ccd29c0",
                "outputId": "821b719b-337e-4866-c9f8-ece1081737d5"
            },
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACBJ0lEQVR4nO29d5hlV30lunblnFNXV1fnlrpbqVErk4WQAAkBNskE2diWB2NgZpjPDzx+bzxje8yb8TcPbDMeAwaTQSRLBAmEQAmUuluhpW51TlVdXTnnW7XfH+su731P3Upd4d6q2uv7znfPPffcc/bZ55y99i8bay0CAgICAgLSDRmpbkBAQEBAQEAyBIIKCAgICEhLBIIKCAgICEhLBIIKCAgICEhLBIIKCAgICEhLBIIKCAgICEhLBIIKCJgDjDGnjTFvSHU7ZsNKaWdAwFwQCCogYIXCGPMvxpi/SnU7AgKWCoGgAgICAgLSEoGgAgLmCWNMrjHmM8aY8/HlM8aYXO/3PzXGtMR/+wNjjDXGbJvmWA8bY/7GGPO0MabXGHOvMabC+/27xpgL8d8eNcbsjm+/G8D7APypMWbAGPMj77BXGWNeiP/nO8aYvPh/qowxPzbG9BhjuowxjxljwhgQkLYID2dAwPzxnwFcD+AqAFcCuBbAnwOAMeY2AP8RwBsAbAPwmjkc74MAPgSgHkAMwN95v90PYDuAGgAHAHwDAKy1n4+v/w9rbZG19g7vP+8CcBuAzQCuAPC78e2fANAEoBpALYA/AxBynQWkLQJBBQTMH+8D8N+stW3W2nYA/xXAB+K/vQvAl621L1lrh+K/zYavWWtftNYOAvi/AbzLGJMJANbaL1lr+621owD+AsCVxpjSWY73d9ba89baLgA/AokUAMYBrAOw0Vo7bq19zIZknAFpjEBQAQHzRz2AM973M/Ft+u2c95u/Ph38fc4AyAZQZYzJNMZ82hhzwhjTB+B0fJ+qWY53wVsfAlAUX/+fAI4D+Lkx5qQx5pNzaFtAQMoQCCogYP44D2Cj970xvg0AWgA0eL9tmMPx/H0aQUmnA8DvALgTVBeWAtgU38fEP+cl/cQlsU9Ya7cAuAPAfzTG3DyfYwQELCcCQQUEzB/fAvDnxphqY0wVgP8HwNfjv90D4PeMMTuNMQXx32bD+40xu+L7/zcA37PWTgAoBjAKoBNAAYD/HvlfK4Atc220MeZ2Y8w2Y4wB0AdgIr4EBKQlAkEFBMwffwVgH4AXABwEnRf+CgCstfeDTg6/AtVpT8T/MzrD8b4G4F9A1VwegI/Ft38VVPk1AzgE4MnI//4ZwK64V96/zqHd2wH8AsBAvF3/21r78Bz+FxCQEphgIw0IWDoYY3YCeBFArrU2luT3hwF83Vr7xeVuW0BAuiNIUAEBiwxjzNuNMTnGmHIA/y+AHyUjp4CAgJkRCCogYPHxRwDaAZwAbTwfTm1zAgJWJoKKLyAgICAgLREkqICAgICAtERWqhswG6qqquymTZtS3YyAgICAgCXC/v37O6y11dHtaU9QmzZtwr59+1LdjICAgICAJYIx5kyy7UHFFxAQEBCQlkh7CSogIGD5YS0QiwEZGYAxbgkIWE4EggoIWOMYHQVaW7lcuODWx8cT9zPGEdZyfvpLZubM67P9vtD/BZJeXgSCCghYI7AW6OlJJKILF4DubrdPfj5QWwu84hVAcTH/Mzk59TPZtov5nJiYfT9/mZiYur6cKC4GKircUl7u1vPylrctawGBoAICViHGx4G2tkQiam2ltARQEqioAOrrgT17gLo6ElNJycqSEqx1RDcdgc1lfS77xmJAXx/Q1QUcPw709ye2paAgkbz8JT9/ZfVruiAQVEDACoa1HDSj6rnOTv4GALm5JJ8rrnBEVFMD5OSktu2LAdnGMlLg7jU2RumzqytxOXsWOHjQ9T9A6Wo6yauoKJDXdAgEFRCwQhCLAe3tU6Wi4WG3T3k5Ceiyy/hZVweUlYUBcCmQk8M+rq2d+lssRnVqlLzOnwcOHUpUTWZnTy95rTSJdrERCCogIA0xMDCViDo63MCWnU0paNcuR0S1tZSWFgprgcOHgebmRC++qEffUv+m7zk5QGXlyhqos7KAqiouUUxMAL29U8mrowM4epS/+8eRtBWVuoCpHpbz+X4x//G/LwcCQQUEpAGsBU6dAp58ksQwOOh+Ky0l+Vx6qSOj8vKlUWv19wM//jFw5Ai91gDnsJBKlJYCO3Zw2byZA/dyYniYNr3sbLaloODiB+nMTEc0UUxO8h5EyaurCzh5cqpnZSqh6//whzlZWgoEggoISCGspcH9kUeApiZ6ie3YkSgV5ecvTzuefx544AGqp269FbjuukQSlEOCCCu6Ptv3i/1tYICSxfPPA888Q4lqyxbgkkuA7dudNLFY/dDdTalVS2srJR4fWVlUv5WUkLCSfeblzZ/EMjL4/9JSEnG0bQMDJKuhIbdNkwe//6Lf57LPxX4vKJjfNc4HgaACAlIAaymlPPoo7RKlpcDttwNXXbX80kFfH/CjHwHHjgGNjcCdd1KlFkUqHRJe8QoS5+nT7LejR4GXX+Zv69eTrETscyUF39PxwgWgpYUODv39wMgInSBycrhE46HGxijlNjdTJTc+TukqN5f3zxj+bzYSm49K1hhOYIqL5919KxZpX25j7969NuTiC1gtsJZG8kcf5cy8vBx41auAK690KrXlbMtzz1FqmpwE3vAG4NprV4atx1r239GjXJqbuS2ZKlCSR3Mz1ainTpGI5GAyMkKSycqiNFBUxKWwkJ+lpXQ0KSkhMfX3u8WXZMbG6MavDByZmU4alLu6SEykl5s7PXnpMzs7pV09BcPD7MumJn6+4x0Ll/KNMfuttXuj24MEFRCwDJicBF58EXjsMXriVVUBb387cPnlqZFIenspNR0/DmzaBLz1rcltIukKY6gCrasDXv1qEtChQ8CBA7yu9nZnx1MwcCxe0zgvz5HQxo2MBauv52RBZKTP2VzxYzFHVn19ieSl73193G9y0pFYfz/blJHBgd4PWs7JcSSWk0PSTEZeamNx8fyeIZ0rFpt+GR/n58gIpczz5526s7vbqWCLi4EbbqDKdSkQCCogYAkxMcGYmEcfpe2gpgb47d+m910qiMlaDuI//znX3/IWYO/elSE1jY7Sdbu3l59tbZSEzp2jeq6zk6Q0McFZ/uioU7VJstqzh+rCbdu4baHqVHnZlZdPv4+1HOiTkZe/PjDAto+OcunuJlEYQ0IYHeV1jY3xvIWFibFueXn89MktJ4cSmLWOdETU07Wzry+xfTqHr7KsqnLZM5ZS5RgIKiBgCTAxQfXZY49xMK2rA979bnripYoMenqA++6jN9jmzZSaZhpYU42uLvbh0aOcvXd0cBDXMjbmBubSUkqjGzbw2rZto3SUleVUgUeOUCV1//2OsC65hBLkUtr9jKEKLD9/em+3yUkSb1MTpZWmJuDMGba3tZX3bmiIBDIywmMWFdHmVlbm1Jidnc7Tz89hmJ9PYikudhJYfj73FXl2dro+Xb+eUmVDA5cNG0hI2dnLm5Mw2KACAhYRsRgllMcf5yx0/XrgNa+ht1mqiMlaYN8+4MEH2YY3vpFSRDpKTb297LsnnwROnKBEJCkgL4+D6rp1JJXGRq7X1c09G4M8Ao8e5fHHx5fWK9DHyIgL3j1/nktLi4txGxqilDQywnbJ4UIkIzLOzKQKs72d/ykqAqqr2fatW9kvGRkkGxHa4KDLMtLezn4eHXVxZ6Wl7McNG0jsjY0kJKkSS0qWVuKfzgYVCCpgSWCtU0kMD/NF0rq/+CoSLcvhVr3YGB8nCfzmN5yNNjaSmLZsSS0RdHcD995L77etWyk1lZamrj1CLMZBWV50hw7RjfzMGUoTBQXsw8sv54Ape9NipmiKxegwIemqr4/3av1652gxH69ASUFtbc6BoKWF3zs6HCmMjnLfrCwSb0YGPwsKSEByoPDtUbm53ObbxtrbaUM8coSkNzrK7YWFzs09M5Pvma+qy87mPnl5XM/Kcjap8XGnCpSjR2Ym9ykqogTm2+kqKmgDXKgkHggq4KIg76TpSGa6bSMjM2eazs11KgY/KBXgixMlLS2lpcvv7TYTRkcZm/PEE7yOzZv5wm7alFpishZ4+mngF7/gQHPbbXRhX+42KYN6Wxtn7/qUvUjbjOGgd8UVwI03kphKS5evvb5XoFSBQKIqcONGvgvNzbR7SQJqayNZdHfz+ZeKzRinEisocGo+2Yas5cDvu6YXFfE5LyvjovXycpJOsme/t5c2zsceo72zrY0EI9JrbKRqWWSfk8Pfx8YcIUXXh4dJ2L29zk42MMB7NjDA517U8Q//wNRaC0EgqIB/Q0cHH7jZSEbrcyEaf/FfxmTbNLMTZPzu7k6++KlfpI6YjsCWK2v0yAgJ4Ikn2EfbtpGYGhuX/tyzobOTtqYzZ6j2ueMODm5LjaGhRBJqa+MyNub2KSri4NjdzYGvsBDYvRu4+moOouniUj0wQKJ65hnawTo62N6xscSsGllZHPCLihKfbxFOZiZ/l3osPz+ReKLrM12/tezjjg5KnXL17uri77J1jYw4e5K1PKbewYYGOujs3HnxUo8cLrq7Sczbty88xVYgqDUOpdJ55BEOXFHk5MyPZLQstTRjLQeG6chrYCBx/9zcmaWvhRrDh4dpH3nqKQ4EO3ZQlbd+/cUdT2UctExOUo1yMf06Ocl2/fKX/P9ttzG+arEJe3ycA1NUKvLvRUEB1XG1tbSPxGIcTI8e5cBZUUGJ7sor00PlCPBZa2+nOvT0ab4vmjh1dPAzN9fFRBUV8btPKtnZiVJPlITmUjMqFiPBdHbyvP76yIjbr7iYz11Dg3NqEFFMTDDw+tln2eeDg05VB3B93ToS1c6dvEepRCCoNQpr+aA++igHiJISxi3U1ycSzXJnL1gsjI3NLH35LrXG8Po1YPhLcfHU2BD/e18fX/aDBynxbdhAtUZpafL9pztOdHsy6TQzkwOGUh3pc6aUMh0dtDWdO0d11O23L9z9d3KSfRgloq4uJ0VkZbGtKuGhz6Ii9tnzz1MC6eriJGj3bhJTY2PqnTSsZb8pO8ULL7D/pKqTGkxZJWQ7ysx0AbeFhXyeKivZD9XVLrFrfj73y852jh5az8riMQcGnAqtr4/PslTe6p/iYt57SWUTEzzPTTclz6QeRV8fHXcOHOD1xmLORb2/n9dUVkZ76aZNXPcnTv7zq3V/2xvfuHAJPRDUGoMyUj/6KNUBZWUuY8FKJaP5Qq6305FXtOBcMoyOOnvD5CQH38ZG5+mVlZW4aJY6l+3JthnD2bKymPtSSUlJImnV1ZEgn3oK+NWvOPi96U20NUw3+E9OUk00ODj7Z09Poj1FpTx8IqqoSPTuGh9nCqLnnqM7u7Uc9Pbs4Uw9lTWorKXN6Omnqbo7dIj93NeXmKpITgpyGsjI4IC+bh2JYmQkUQUux4dYLDHjt0hJ/aNPpUbSfhkZLouFzu2nlJK9R+mWCgocOZSVcbJUV+eCkZMRiq6/t5d90N7OZ0GegUpSC/C7yLa4OPGa/GdW6+99b/Ks7fNBIKg1gmjGgspKEtPll6eXc0E6YHzcSV++CiQriwPPs89yEAMoLb3ylRyY9XIuRzzIwMDUEu0quzE46FRmW7ZwRl1VxQFmfDw58fi1o3zIflFQwMG4oICDX00Nl+rq6cnFWtpDnn2Wz97oKP8rFd58bB0aRM+dc0tnpyOPuS6Zmbyvhw9T6j12jBqEnh4O2pmZvN6qKjq2lJc7+05fH/u0pMSVNNm5kxOT7Gz+5i8jIy7BbHMzP9vbKTXKoUDSsi99yYMuI8ORi6r3ysEiqlKXfUvPRlcX73V+Pp+BnTtJVjp2MkLJyuK5Tp6k5Kj+3biRE57BQZfdorSU1797Nyca0TFkeHhxvG4DQa1yTExQRfHYYy5jwatfnbqMBSsVPT2Mw3n2Wb6gV11FYlqONECTkxxshodnlmz6+igBvPQSB7/ycqf6EQoKODmpqiKprlvH9cJCR0D+p5xLens5uMreIVfjZJ8jIxz4jx5lv2Vn0yZ32WUc7ETiMx1jctK5ZZ87x0/N5HNyaFupreWAKknFX4aH+by3tSXabEQyfv2s4mISbWMj34+JCZL++fNuglJaSoeXK65wnntK6moM74FqN+l8XV2JfZ+fz74vL3cDen8/+1XSuPL4jYwkkm9pqasjVVLCY4lw5WEnL7vh4eRSUWkpQwq2bOG67vl0S2cn1X8vvMD2lJdz3CgoYHuPH3dqwe3b+S6MjlI1ev488NGPLvz9CAS1ShGLcTB9/HE+qOvW0Wh/ySWp1/GnO2QUV3bss2f5ouflMf3Pddc5r7O52pQudlssNnPNJbkqT05Sqhsc5ED6mtdwMNRvAwMcnLu7nb1IyM93qjnNxCcnuY9IabZ6Q5OTzotMx1aQZ3X17Orj0VHnvqx0Orpu2XNUmE8hBZpgDQ66mKLOThf0OjTkXKQB5/BTUcFrratju3p6nCTa38/7INdueYb6JTKsdR6tg4Nse3a2s2PKOaGszNmIFCArm5JgDNsj1Vl1tSMhZXHwSc/3wpN6z09aW1DAc+nZURLbkyc5aeju5v/kpAEkHs8vZyJPw8xMXmdXF4+VmelSGkmqkgt7RgavYcMG4K/+is/iQpA2BGWMOQ2gH8AEgFiyRvkIBJUcY2PA/v0uMHTDBkpM27atDGIaHKQBfWzMvbCVlUurhtRg/MILJPTnnuMMsK/PGcA1WGVlcbDyJQ0tyVyBjZmfrWmmbfn5U6WcrCzg17+mTTEvjzn0du2a+XrHx3l9R49yFnz6NAeZ1lZnlzCG/V5fT2lh2zZKQRs3sh1SPfkqPGUv2L2bKqWyMqfC8j9VzqKlxcUM9fW5JKkasKuqnLOJ4mxEJiLOri6XdXx83KnpFHRaXk4V1MaNPN74OP9z4QInHrI5ZmSw7Rs2uMwLup9yClHpDRG2tc6NvK+P5KGsD4DLSp6b6zI+KPyiuNglptXvvhQZLV+id3dign3R2+tSGHV1sV8GBpz9KBrIq//HYtxvaIjbCgoS4wj9MUJqRS3KAyiikh1OSWvlYKHr//u/59izEKQbQe211nbMZf9AUIlI18DQucBaurjv20fbgGawegQzMhK9oWT7mA9xKcamq4ufnZ0857PPcmbZ1uZUIwUFlDgvvZSD8q5dbJOkEM30e3vdoJqRwUFHqjPN0GtqOLNdjHvgz3InJji4/+hHHDQvvRS4+WZn2FZi1I4ODuZRNZcGHcAl+iwpIeFJZdTby/0HBty5c3OdQ4Y8yxSzdNVVfO6iquOhIaeqO3eOpDY6ynNkZ5OQiot5nJwcHlMhBLL9+GUsFMiqDAsigspKtkteZ3V1JK9z50hGZ844SWt0lORQUwNccw0l4y1bHCl1d1PqOHqU5Nvezn7QuTIz2ceKG4xmXJANSjYjkZEcHuQN5ydq1T2ZTvUJJKow/VhESYcZGY6sJbllZrJvRY7Z2Tyf4tHUp/6kR5KUUizJ6cNaXpNUm0piqzboPCMjwJe+xH5dCAJBrXAMD9Nb68kn+VBs307nh3QIDJ0NIyOUlvbt4wCQl0fj+d69nH1rcNXS1uZUFMBU4vIDGnt7p5bGlsrk/Hl+Hx52RvEtW5iH7tWv5mA7lxicyUkO0mpnR4db9+NScnNJXJIKtD5befbhYZfB4MQJNxBMTroBNyuLbc/Pd/YoqZ78QFjNlpNJftHquOPjjkDGxng8JWLt7+d9mJzkdVx6Kcmgqsq55Y+NuUmAMilo4Paza8uGo4Gwt9flh9O1aKIi1+3iYldKQtLspk1cKiv5/7NnubS18bhdXY4MiospFe7cSVVcbq6zdR09yhin5mbeV2U9yctziwJwlWRVqsfKSi5KN5STw7ZLpXrhAvvDGB5n40ZKadu3sz0K2h0fd6meTp3i5OnkSbavs9ORoV96Q+pT3a+JCXdPJcFrgiTbVXGxa3NmprNXyb6Yn5+Y4VyOMbpGTYIk3arPlRaqqgr4y7/kBG8hSCeCOgWgG4AF8E/W2s8n2eduAHcDQGNj49VnkkWWrhEMDlJaeuYZvoSXXsrBtb4+1S2bHc3NJKUXX+QLuX49Semyy2bPGNDfz8H6xAm+EE1NTuXjD8gaSPLynPQjnbwCGa+8Erj+ehJT1BNNRvqODmfUn8uiWevISCJh6dN3Yc/MdA4LIq+sLO6r61OA7rZt3P/sWTq8dHQ4w7m8uzSzrajgomNLVeYb0UVkvqOF1iUV+m7EBQVsR3a2cyI4dMjZm6zl4CTvMc20Fd8jO0lRkZPyRISy4+ga5LhQVpYYO6T2yB4iB4yJCfa37HVyTJBTifpaxDgxwWvt73f9oWdHklhdHSd5DQ1OYpfkWFrK/pivVDw8zD576SVOOi5ccJOOoqLE0hnKzwe4RLgNDSRiBTj397MPenqcxKN7oevr72cfiGgkpflqRVUCHhnhe6JcmDU1DAO48Ua+L8nUtTrexASXCxdI9GfO0Aa1YcP8+iiKdCKoemvteWNMDYAHAXzUWvvodPuvVQmqv582h/37+dDt3k2JaS6BeanE2BgJad8+SjDZ2fSIuvrqqaQ6MuLccf2lu9tVKhWky5druGaf+o9mdDk5HOx27uR5d+xIVBVOTDg1lIjPJ7yLgT8AaJmcTEwdNTjoJESVilDgcHW1q5w6NOTUiwUFVKft2sVrkh1EedykIlMtIX0me6UllSjhZ3RdwaAnTwI//jHJUapN2SREaiUlHMgqK13CU6kaNZDqGkVgiiXypZHCQufBpqBZSdAauGVvkmrML+wntZdsUZIaRLraR5JQdbWzs23Y4PaP2nKmW9ekanCQUk53t3PUUI0qOX+ovcoKoWSxfX3OflZURDKShGVMoqpWqj09JyqVUVTEa/TzXo6MODWdAn+lmpUTie/RqElddnYimVVWUlKXS/1Mz7v6+fd/f+GZKNKGoBJObsxfABiw1v7tdPusNYLq6SExHTjAF+yKK+jmvNBAuKVGWxtJ6fnn+aLU1FBauuIKvgjW8qXz41va293/lWNPUoHUC8rE7Mf/9Pa6OjbZ2dx3xw7OPktLeX7p3c+fd3YkeW5J/dXYSG/HXbs44GrwS7ZEK41GP8fGpnrljYw4g/uFC87IrwF9cNCla/JzvEndVFHhXI1loC4v57rsIypL7hNNlISKipLb8IaGGFT70ku00e3bR9IeGnJOGzKOV1RQssjNZZtHR11eucnJREcTIJFMVTBPg6kkAs3iVV4CSCwrIRWaJDFldNB1l5a6QNJYzKkoAfdcRDNxKHO3VJuy9Uw3DOra/ImT71KuZ9d3VNBzofb69qHSUrZN976z00mFeXncp6aGz2NNjfOii8Y1JbNdDQ3xWdOi4GFlNgcSJX5JXZKE1V6A927DBhcDtWUL2yRHi8LCxQ1fSQuCMsYUAsiw1vbH1x8E8N+stQ9M95+1QlCdnfQse/55PvCKv5lLkKOyA0h1M1uWgOFhN+BFS0jr+1zq68RidHZ45hkObJmZfJj37uWAcP68I6OmJhckmp/viqDV1Tn9eEsL91PRNt9tODPTlcrOzWX7t22jynP7dh5zYoLH0DmlK9cMUh54qlDqD9oyBkddcP3vwOx9MjbmZsxKdCuDuEjLHxSkjvKrpsoGIiN4tDyD3LHXraNEsGkTpVO5ZqsUgq/OHB52QauHD1O1qEwVctXWTH39evZpXZ3zBDx3jqSitvhZFioreR8uuYT/yctLVEv5bun9/c7mowwKJSWOeKVClEu3VJe633v2cLDs6qK0d/Kki9lS6ZbiYrZBEl1UQlbGcD9HoyTInBzet2PHuJw65SQZqVNVcqKkxKUL6uhw+QglydfXcxIklW53t5usKOuEnmVNxjQ5yMnhvd2yhUtNzdRnb2yMnplShXfErfrFxZTKFAtVWDj1ObXW1aU6eZLXKeeSCxe49Pe7rOiS9GUXVBBzSQnwnvcsXLOTLgS1BcAP41+zAHzTWvvXM/1ntRNUWxvVKS++yAHz6qtpL8nKSiSWmUhnpuwAmonKgOobVeVCLFdSH5mZjriiBDYxwRfipZfYhvJyDk4VFXxJzp3jy6rZdHU1X1YFL2oGqXILTU2UcHTedeu4n/LstbU5V99LLuFguHkzf/clsubmRGOvCtqVl3Mg6O/n+bS/KpbKrqWBR27mspX4i9Qavpvw6KhL6Nnb62Kd/OwA6nO5REu9JU8wDUqKxfHVWbK9yIFBfSepTd6QIjDVEpL9R/sBjhRyc91vkjYaG3mP5I2ooGEV9ZPNCnCDl6TTWMyp0tQnyqCg6y4q4nOwfTszm+zcyXMfO0bS7Otj27dt4z61tY7Qs7I4aJ48yUFeDhWbN7tBvLJy6iAuyc1PceXnbuztnSoljY87CVbenZdfzue4s9M5Cpw75949SeT19fzv2JibcPX0uL6vq3PxUw0NPIff5pERko5PvgD7bvNmpyE4c4bHFxlu2uRIqbp6fnYzldQZHOQxz5xx78eLL5K8lA7KD4nQBO7rX2d+z4UgLQjqYrAaCKqnhyq7o0cpNdTU8GV7+mnGp0xO8oFdt84NRsnge2j5ObsEDUq+Oso3nEq1MTzsBuGMDDfjlp5ddhQ/qFIZBuQZNzbm4jw0kKrwmgzeZWXcR7pxqWB6elzku85fUcGX3BhnTwG4TbPorKzE9sghQnYGvwCbAhclrcjzyXfvzc1NjFERUWdl8Xx1dRxwCgvdSymVVGurK8etMgxyQwecyqeiwqkes7NdHI1m0OXl7voUOCxnAt+uoPsmW4JUXn5OOKmKfAnQJ1Q9Q3IyEKnJFhNNbOrXJhJpaiDXc+p7EsoBoqSEJHP55ZxMbNnC5379eh7n7Fk6ERw+zL7LyuL+u3eTEESecv0+doxtzszkoC5CWr/+4tVMg4N8915+mW2RdKg4OEnk003+qqp4TWVlbEN/v4v1kgqwpMSVTG9o4HMw33IiLS3Ms6i4vQsXXNaLggL3/Eazgqj6sK8q1aRUkx7fVV559vzsH3pejOFzeeIESUsSrqjjy19mcuKFYDqCylrYYQOmQ2cnve/uvZeBoUpDogdDAYPbtwPXXssHWPENfiS7nwBSsxxl6Z6cTCw0Nj7uHjDpwDV4KWmqX2zMt1/47qZ+lm29fN3dTjWjh9lf1+Dsu7r6A6R/HbID1dSwDX19nDVq0FTcisivt5dkoGu01r10yqIs76TBQTdrF1kq7Y8fwa9ieCIA2UHa2zljlHv36KgjJz+7tUhfxebkHl1ezoFTEwDN3qWCqqmhF+PmzZz19vc7Nc3QkMsPJwO97FVSxfX0OHdmRfTLviF7jF8mXM+M4lt8pwFNCHRvRDLROCzAEVlpqfMa1ECu1DdjY05C7e8nATQ0sE9iMQ60hw/zWrKz+ezv2uXqCcVi7IcXX6T3m6Taq68mganY3sVgcpIEcvw4Ce/8eSeV79njgnaj5TBUpl21j9Q/nZ2ujAXA66mvp/ZDEtJ8MnxPTDhtgUjz+HE+97KfyoYp9bv6u6nJvVe+OhqYalvz7VZ+Tkn/Hfb3UTiCpHDfq08qyiUtBR8kqMWBZr+/+AXwk5/wRezqcjcyI8NFn2sW79s1NMj5CSElIfmDvgYULVEXUJGMH1GuGbM8pvLznfPByIhrh++5pLQ5MgrLeK0ZWkkJByoZ4BU46NtYNOsXCUgtIzdkOQpodqzBUIThByfKS0s6cBGDSEoGZnm7qZ25uS6xqhYRrtSk8naSoVh9KzfgwUG33U/gKaN7UZEbOOX26xNpdH/1vZw2MjJcAG1Rkcu3p6wKIk8NQoAbRHQvZBcoLnYeXiK1/n7nhNLYyH3kJZmRwTb68U1SP+o4Fy6QOE6dcil4Jibc/aiu5gC/cyclm5ISTsruv59EMzLisjfs2cP0TNdf7+7N6dPONjYywn7atYtE3th48QPg0BDbfewYB3vZ2Roa2F7Z2aZTh1lLwpA96tw590xWVTnJSPkCZ2qnninZq7R0dpKETp50noFSjyrhr2KWlMBX7vh6lyX96vw+ifjBvvIm9Z1TfOKJjin6ruckWi5Ek+1vfhN4/esv7h4JQcW3BLCWD9dPfgL87GcuaaYIQTdZEpMCDv3obKlvNABJilB+LP+hEHnoAdEArQe0pITfx8f5EMp2oEU2CT+7wNiY85Dr6XGEpCh+pZDZsYMDnNyg5TLrE65UDgoKHRtzhLRtG9vX1eWCJDVgSBL0Z+y6Zv2u9vqDtO/VJEJSAkzZQvQSinDkiq0XU15hUo8orx3gylxrBikVqdSDUo1KnTo56dyZJWXofrS2OrtFX5+TqHwPNNUg0rXpvqpia3k5B2wZz5W4UyW8T5zgIH/qFO9pdzf/t3MncOutHPD9kg/t7c7LUMZ7ORT4diyA55FTizIInDjBgV/OKMrAIWlfz6cIwC/joHdD4QE6h9SJmkxIjetLbsqwrn21yMPw1Cm2TRm5CwrYV9u2UUqK1tXSxETu48rNeOKEKzki6VE2WT8LQ3SRaj2ZrVjvlt57EYhPEnrWABc2oWS/VVXOQ1JjiNRwvtTre5hqDNGia9GkKDvbqd3Pn3caGtlEi4pcCMTAgFOxa/nSl0LJ91Q3499gLV+A73+flUtPnOBN8jMn+0ZxzTz84E4/bkaG7Shh6aHVDNv3PNMictJ3axN/y8py2zTw+TpnXypT+zUTs9YNOCIhEaeIVtel4yoQ0lcB6By+JClJTseQ3l8pcDRgKzWLP1DK9iYiVHS9JD7ZQ+R5JSlUhnq9mKrr1NjIY5w/T/f+48edw0ZpKWfHl1/OF1CeYS+/zP3kXFFezsHvFa8gGWRmcub9wguUDE6edNkTRG7KKKD2i/TlECD7nspBXHIJPyXpSFI8edLZaaTOGxlxmadf9SragWYzmouwjhzh9R09yn7ftYsepWVlTgXofw4NURV66BDfBUmmfkofqaE1EfPvi4hYnnQafDVpUMoe9U80nkdSmO65npFonj+9u2qTn2RVMUNKtaTnvKKCEpKcR2SL1LvoB8nKbqqgYPW34r+ysth+BZv7UrmeacCpa1VDq7FxakydJload3wVv9Tqfn5ALSJxaT38WmkAf2ts5POycyeffQUqNzVx3Dt50qnbreV79O53k0AXgkBQC8DkJD3WvvtdJus8ftzNrP18ciISzVTy8vjSaQakgMy8PKd+kVeRXt7xcXcM/yWSm7EISIO4DOp+8GL0lvqqRBGMP9Py7U3axycTEZ0v/kd11SIrSYe+vUnnThaZrv4FnKODsiRv3eokBOWOU0Ck0utInSbSldTpB1n6RCgSNMaVapArumbp1dW8FkkXkiylTpOjS34+2yHPst5e7qc+1eSjtJT3f8MGF48zOsr/nD7N42dkuIh+uWm3tzsVTyzmjucPjllZPL6yrhcVUdq95RY6HSQjptFRZ8dSDaNTp1xCVan/ZG+TO7zOk5XFfvO9Fo1x6kHtMznJ42lfOXn40oHUmgDPUVLi8uxt2UKJTXYz2SNPnSKJnj7t+khpgfx3xPc48x1YZKPVRFCEWVTkskvIjiiVuwZ3pdHq7nbPLcB7qslDTg6Pffw4Jdpjx1xaLNltSktps5K3rtqvHI/RGkv+O6OYKUnSGza4iZdsaMqi4YefnDvHicfJky67RUaGmxz4kqHU4D097jnLyWHbGhp47sZG/mf37pmrPc8FgaDmiYkJetl997vAI4/wZZBUo8Ha93yRPlgu2bKTKNGkZvayY6g2zbp1TsyWHUPVTH31gEpv+0krgUSi8B0WRFayMWhGLuO94kO0Xck4RVp+vi/NWuV6rJm+Zm9AIkHn5fFllwrBJw1felDQquwyfhYFka5sX37GaV+dIiOt+k+DkLwdNUEwhu2W3l+2GblZqw8UM+IHY+q6fS85STklJc67T5KDtS7N0iteQQmssJBE8PTTHCC6u9n29evpJPO617lsAoJUg83NlMaee44DsyRLGcsHB9mm8nKn9gNcNnb1m0qMd3Wxf+Uwo/7yHU90f+R1qEmB3gFfdbt1q6ubVFiYGI/W18djyiGitNSlwHrySZLN0BD3kQZBJKZnt6QkMWegnjXFEOm8ule+1O6Tm1Reeu6Ki10mh23bXDySVOGyhfpLUVFi8Kwyjp8/z3t05IibOMmJJS+P90ZeoWVlLtO6tSTFnTu5qFru6Kh792Ur1nVpnJA6trAw0Vtw/XruMzbG/j16lCSp8h/r13MSs327k3xaWpxq8+RJN9EqKnJeuYWFrl2aOAHARz6ySjNJzAXLSVCxGPDQQ8A99wAPP+x08v5gJbWMdLkykCtaXp42fmaBqOdLlOAkOUjU10Cv4/jqQ19V6KsstC7ClHpE0oLvej4x4UpF19Q4m5GvKvBdmZUIs63NSY7qB82s5NXnS3K+/UR1aTT4yb6gFwFw5/MzMSgaXv2koE4Vg/OdQdR2vUSyi/g6fkGqJc365cmmwbegwKmIKirYl5IkFePS3MzjZ2XxRb/mGuC1r2X2jEOHGN/25JMcHCSV1NZyoL75Zue9mUzKsZaD3sGDlN77+91AX1bGgeSZZ3hv8vNdGinZTJT7TrYtkZiu3Z9U6bok+QPORV5B1Mq9p4Bnzbzz8qh6rK7mQH3kCG1So6MuUakIRmToY2yM/djWlhgDp4z0knA12VKqpLIyJ60DbqLjx+jIJb611U2ofPW37zGra5bTgSZZmmhqX78ch6SLoaFEV+3SUr5bO3dyciIJZ3CQz8KZM85ZY9Mmtkepu6Re1jVFM/vX1PBYcuBQYHtTE48hrYCe5cJC/u+yy6gm3rbNFSk8dcotUo9KpdzYSCLzNTnRRc/YZZclDwaeDwJBTYOREeCHPwS++lXObpVHy4cGY80wMzKmes35Kivf+0VEJPhBn9Gu13afpHR+nxB8+49/bhGnXlbf0cAX+TXLl+TlX5sknFjMDQ7yqMvJcRKJYoekevOdEaIxPJLIdL3yPpMu3Pda8uMv1H6pZERCGkyj1Uf9LA6awfr3WV6UfryY71KvnG6SOKWirazki63CdyralpvLwVTb5CwgG4RUY35GccD1Y309Bw8Rtrz32ts5qMqRYt06DmQVFc4OMDLi1JGxGAdMqYx1z3T9svWIiJQLDnADv6+S1jUrf1x2Ngfba68lQcp1/uWXGd/n2+0U+6NSGL4xProoUzlAknr0URJcbq4ri6F4t40bOWDK61J57+Ss09fnaj+dP+9qOSl2UHFCvt1VE0D/+RO5+ymufPW279UmD876etrprr2WNr+6usTn7je/oRbm1CkeW84WqrelbBz+ojRX/iTUX9S20VGSkwoVNjWxv2R60HOsa9Z/NC7JEUJt8u3Vc4G1wB//8SrJJHExWAqC6u0F/s//YYDZqVPTJwtNRiQzdZf29yUspa9XAbiTJ/kyZWc7Y/fwcGL8iGZjegGU8VkSmvbR4AK4AVcvWZQ01bZo/i4/iBNwEpakMcW5bNjgFrmXy4AtCRJw+cBOneLL0tzsXGd91Y1PQJq5iowkpfieT1ok6fjHkwSpa9Fg7JMQkGiL8CcBGnx8FZAGBv/+a/KQl+c8+LKznZTte2T6EqRfA0nBn7IdAs4moXpJGrDq6zlItbTQEeLMGfavJhWyvfmSqyAJOzrY+M41ChWQ95xi4ZTHUAO8PNw0CZLKVddYXs5nvKjIqRGV2VyZHhobEx16fNuQ1GeqJ/XUUyQ9a0nOcgWX+s9Xw0n1qFgzqX1lL8vLm1pWxHe91j3zNR7J3Kyjz4Cvsvbz5PnaD6noANcfSg8ktaTvcTob9JzKsUMSpmyAZWWO4HJyuE9rKxff2cnP6uHn/FOmdf+90afeNU1O5PXb3Q380z9x4rIQrHmCOn6cTP/YY9NnapgN/gAobxg9dIoUV/lnzVb0Uoh0iov5Qh08yNmeUq34pQeUiFEOAZJ6NMvV7xpAfbfwvLyppOVLM/5/gMTEmb5kFQ3W8yW/6Evq24eiThd+jIafW85PhRM9vv7rG7z9AdcnJ18C8gcTeS9mZyfGj/i2QZGwAnqVZ+3CBUoHTU2JEqQfU+araf1sDT5hqw/VLrVFhvjaWh7Xz86ge6z0SZIGrXUDi9832l+eoH4NqLo6PpO1te7aBgbcoKXM4RpsVBtJkysl7i0t5f/OnKF0ojINtbUuM7iPsTGnsurocNKz1M6+550vWfuQF6lqeWnSoePrmVW/a7CXXc5/xqKDrT9h8dXUfvaFqPu6Eh5L6lGOQTkSSPWt4HhNDFTcT0TkZxJXdndNYBT/6E+KfElO5gO9JyrPsmMHiUapoE6ccFWTjaGUrXsJOOcYPQMi9qysxDyDsjn19zu1uQLXFdZSWgp89KOUlheCNUlQQ0ML041qANAgoIc5mXrO/x514dai2bNmvb7XnYhDx4ouPvRCSbpSYKWvAgPci+BnJYjFnMOFPODk2lpb64qqaZYZrSXkz0R9z6JkqkzADd7+QK51kYjfV/41++oM/zxC1Ibg3wffgcU/h+/U4huAAUoqQ0POvqKMDkNDzkFAgbPqB//aNUHw1XoakI1x7shyptBg44cR+BKCJGTFl0n60b3UIFJU5OoINTZysGhsdGSWDP39LmPBiRNOKh8ZcV53kpykalR5iLo6tl22HWOcW3durns+/HPJTVwSTlUVB9XGRucRKZWxkrwq+/3hw65CsAhH/e2HWOj+6/7q2dcEpbbWPeMbNrDPKiudum06KWZwkKrM/ftdXj09f5rYjIy4CU1xsQtaVqXf5mYXbyjvR71PmqhFQzl80lSWkqoqp8EYHOSEQYl4NblU1nsFtesZ8wnYJ+KxMbZHOSpbW50ThBIT19ayv0Ry4+MuGXJwkrhIgurrm1vF1LnCH2CjM2df4oiq0URsUgP5Khktvi0LcKokuX/KhbSvz3nzSZfs2xB8xw0d1ydBrUutJq8rDbyaxftSgLzYNPD4SSP92ZaCgSXlaEDXf/2sCFFJy+87X6WnNvsvrQhXKk/NwmX49mfnvoSlQU2DvPTtkoo04AHOHpaZ6VSxUoPJO1DS7sCAy9iushq+Q4CIx5+4qD1Ru4ZUhrIfyLvQz47hZ8nQNWvQjt5rf5IyMZFYXl0elwomVmYLSca6hzq+n+FD/SV1rCY6creWyi8aYD40lBheIXuYnDXUL8pgIecMlfxQ3kZV6c3IcPck6vzgv1s+YUqS1z2Wik7Pn/qxu9s5cAAuDEGTgcpKtvHYMVdUUUmR/dIV0TFE9zIvL3Gy5qu25R2oJK16h+Q1Kw9BuX5HJybJ8vH5E2BlMm9uduTpJ4yOxVycVGenK2OSmZkYm3fZZXQEiaaImi/WJEEBZPef/YyzxMXEdFLSdNuBxJmdJDHplP1BKhrUGD2n7yjhewz64r+vBorm3BL51dQ4CUL2EXnmyFlEbdbMbfNm6psVf+Env5Rqp7PTldY+cYJqB5XlVkJRucpKJaL0SL7RXgO0yEAu64pFUfyL+lMzUXlraZCKxZz0ohm4BgTfKcO30Ykw/cmE7Gbqe5Gpr0Ly9ff+PZX0I+cWX3XqBxNLjaTcdL4dUffSL7Qn70k/X5pP8CIZn8h9yTM6SfC1A3l5lJb8yr2qTyU1nd+fKj45MsL/NjTweSkpcSolqf3kTCKJUpMjvy91X3wnIT+9kiQgX60NTJXA1T9+7Sk/Ji9qu5XWQDY3qb1U5l4xQsqEokzmW7cmxrtFa3IpuHwu9iZhfNw5gUgSlQ1bhJSbO/3/reX/mpsdIZ0/n2jm0PMM8Pr9BLlFRXwGioocgUoDo2fnD/+Q0vBCsGYJykdTE/CxjzFfXrTyqK8qiHrkReETkD5nU2n5MyW9cH76Ec0G/ZfTt2FJdecPInppVcBO6hUNGNnZidkT5DZdWcmBpb2d51OpaSW3bG93MydJIWVlnB3u3s19lS3h/HkOPLJP+K63vlOBbD+ZmU59o7Q4vnrQl0D9gFvARetHiVv3wO9vwScinwjU17o+tUGEFpW8fEQ9uTTp8F2UNWAq7Y3sEprx+yl7srOd5xWQ6EYNsB3RXIJKuaR0VHLaKChwqp+WFicVqQ9lS/D7UFKkJj+6DoD942fD9snPJ3RJJXJ2kOu1+l72lspKlylCpOF740nSzMpyGT/Kyx35aZIRJVZf9RmNCfTbmexTExERGOBUasXFTqqemCC5qmaT+l0SptI7yQFFAa2bN9NONN8yGNPBt6FFMTTEe+9LR5J+5HmqybGefbVJAeX+UlycvA1jYzxHUxNw3XXzz9IeRSAoD9ayguif/znT2wwMOP2v3C5VW0YShe+GqgfEt4kkIy3/5fBVDr6eWQlYZbOQ55EGCalzZDvq6Ej0/vJn7BkZJBE9WEpFo0BIDRQ5OXzx6+o4iO3bx0Dk/n7nQabZ6eSki7dobXVxNuqvKDHINia72MSEy8bs2yZ8e42vApJEIKmqtzd5jjrf60/E7xO5Jhnq96hqzXfM8J1JNGuWGs9/of0JggZkSX6SNOWivW6dyxaumbpichTv4qsWdf/98AD1jS8t6VoV5+Wn1/FVQHK60f0eG3ODsWximhxVV7tyL6oFZYxLlHr6tHOe0DXp3kTtKFGpUL9p4pRsImEt21hW5moaVVXxeJrs+DFKk5Oc4Jw/79TdIkN5y1VXOycRSX7KWqF3rKeHMWbKoZmZyfdH4QSDg04joImXMjgoOaykLz0/fX2uaoBUvbrWvDy2RXawhgauy4NTnyJwPWft7S5PniaEgMt2kZ3NcymPnrWu4KOeXT0XIjU9o/KSXLdu4dkgFoJAUNNgaIg59f7n/2Skvtxp9SLX1nL2o3QsetmUp86POI+qV5LNwH3HCJ+wNGuMzsY1K9TAJVuCZqq+6kMz7f7+RA+1qMv15GRiwK10+HqwNdgBiYHEkgx8EohKaZKs/FLmIp/qag4+O3ZwUQ0kGeVV7PDECX4q44FPMP6kQAOurwLypS45KfhFG30p2Z+BGsOXNiuL7fCL4ilztcqbawYvNdDoqFNd9fY6KVhee8ozqO+yzSjo07cp+qoqf2LjX6PieiTlSOrQBERqITljGOPi1qRyksOAJDhfHab+0eA7NOQK2HV3sw3V1VQxXXqpK5vS3Oyqu54+zeuSalEkHbU/6vnys3PLO0wTNzm0+BWR/cnR6Kir2itJy09h5L9T8kzTs5mVRVJYv54DtbwHFWcl0lW4gEri6LnSsfWeqn2+qtyP21KNMo0Rap8vpfoqyokJ91zLA3TdOvbj6dMumW52tku1pAmiyEnqUJ+QZlILpgKBoGaBtXzB7r8f+MIXGDCoAUieVOvXMwuAgtJyc90sWkZFv4SFZo9+RgTZVnzSkNrFLzEhaUTEoIdfg6oeYA0AIhSfkKSyiM5s/UwTEvOzstxArHgIlQaQK7wxTo2Umck2d3S4AcF3U5ckJQcKGVGj9hrZSTQQ+mmcorEm8kSSFKF2K5VM1J6mAcpfNMOWTUFLXp5T06gthw5xOXXKDQQic9+V33fk8LN/6N5IjSLHEt8m6KsE5eWpewY4VbCu2U81oz6UJClpzbdBSurTuURufoVkTSymgwbx4mKXxfzgQRKQXxBSk4isrES3/vx8l2FfnmGyZclxQo46isGSylODuDQOvjOQ/q8B3u8DpQtTNhE/NMEYR9YqxOgnW5XdSXnxVE5DEws5cyhrg19Dy5fi5OAkLYhIWNkzjh7lcvo03yN/QqHJqk9eUkVqEiG1LeAmrcqZuGUL02xdcw1d0ReqgltqBIKaB0ZHGRj5k5+w4ODx44kzsoICivlXXMGHuLCQv+kBjsVcPi5JKv6MEXBxIb47clYWZ6gqudzZ6VQnvpeRBjn/gdWg6pOcBm8NnH4gou89pYHPz4bsz+Sj7t5RV2/fW9H/nwhRg4gGWv3mG/dlUyks5Cxv+3bO0JUOprGRxmc/iebEBF/uo0c5oZALcH09PYy2bOHxJFUmyzwt9ZgPX82mSUFnpyvhrfo8kko08x0fdwOYBqusrMQ4Et/pRIOjAiAlQfmOMnJxV4yS70os0peNS8Tme9n5TjlRGw3gvD4Vy7dxI8lFGRzGxhivd+qUi5eSSlzPmdTiqopbW+via5R+JyODM/err+Z5xsfZl88/z/ero4NtXreOCXOvuYbtGh3l+U+f5mdPjwt/kKej3qGaGpKJMiMUFLAdL7/MRW7xkmJkb1JeQdloCwqcus73BPVd2XXdfsZ1kYAmJ/498b1gNbHwJw6SiH3nId+rLiqpaRzy3yuRvTwke3uddmX9eleiRY4rfrmNi4EmwZIqF4JAUBeJtjbmPHvgAZZdPn06UbIqLuYN37uX4rcGmPx8V1gsI4MPi2Jp/EVSg6/q82d3ihrXf/2sBb6Xl/+yiAT9LNgS8zMyXCbwc+d4TA2qftoX37UbcMSoF0Eze83kojFZY2OcXcpO4UsRetEA/r+ykmR//fWc9dXWOlvHXGEt75XIqrmZ20pLnUpx82Z3PT5kP9AgM9PLNjLCuJyDBzloT046Y35REe9VS4sbYFSeG3BBjyMj7JfMTEcOjY2J2a2VX03uv5pMqF8Vx+ZLm0pkWlfHRR6PkiKU+1BZxaPea3qeJL0rNML3HvQzzW/dyn3PnGF/NDU5Q3x1tZP+5RSiLP6ypVVWcuCsr2fbT52imr2piW1SVgQ/4a1f7FGlIkTevo3Yz7KhwVySiNSkGRmJbt2SXvLzeV5lWFBcoJxL/MB3vxin7/zk2ws18dO98tX2ep98hw3fFq3nU8+o78giBw09/9LG+HWm1Ed6D9Ufkuj8eCg/I4bvaKK+9zUMfuWFz36Wqt6FIBDUAhGLcRb21FMkqgMHnIeMyKqkhC/tjTdylujXqFE109palwCyvNxJEvKak/3j5ElXTkCzfg0sQGKOOBUTlEFXA6DsZIp4V/ySn7VZM9LsbBfXI+OssmZXVXEQUeodSXtyE1a7e3r43+JiDpDV1WyXPH4uXHB2u6ws9seuXcCVV/Icmg37nxc7MxscdCqUEyecLWzLFkpXO3YsPMFlfz+N7AcP8lkwhhLI7t3sq95eVwywpcU528iranjYxZgoQr+01EkkGzcyxqSoiKT07LM8T1+f8yKTmk6qI3leytlDak8tBQW8T8rpd/q0c5n20/3I41DZwquqnA1OpT50HUriOz7uglGVAUIhApL4fFugT5B6hyQ9Sh0octF/facRfzIxMuIyqCvtkX9M3ytUx/MdlDSpEqH7Kcfk3LBunUtZpmFzeNi9Y11dLvOLiFLnlOpNkp0vLfkTVj9zuq/V0OJnjSkuZnt27KCWQa7tsilKihOJdnW5iU9LS2I2ez+0wg+RAJxDjV8QVE5ERUUM5QmZJNII3d0cLJ580s2mfZdeY/jCbtvGUshbt/LmK/5DUDYAPzhUajpJD8oU4Ged7ukhWZ49y5m0PLRUSj03N9HICjgVlHT7IioRTWsr/yMjq8pQxGJOCigt5Uxpxw43OCrIc3SU7RIptLe7OCpldr7qKj7IExNOB9/ZOdVl3Idvb4iSlz6j2/ygRIDtUNmBI0fcPVAGAdmNpPLwbTP+oDKTRNfVxefg4EHe54wM3v/LLychZmfz+fAJ68IFN/vXLFfPkG9vknSt7BZS1WzbRiL0M3kb47w9Zc85cybRE02DW02NC7TcutVl1BBxqjyJSqbLSULqWDnsaIKi50r2p4kJF0rgx9LJM9InimTenOPjfE6U4UPPsq9qGxpyjhFjYy6Br5LLKgOHvA6VTNcP5Jak4wdxa6AWGfvlcpT9XpKT3i/9V0Sqa/CDbKNqdBGmbIN+zTJJWL5ziY7T2spJn2yjUudKi1FcnLx/laewuNjdH+XU0z1dv54T7Pp6TnAk+SsGSwSu9/bDHw7JYlPdjKSYnOSLe+AAVROnT/P7+fPOA8hXA15/PR9OvSjKFq6Zq9RfyoMm+5a8uURgUkmIGBT0qBdGXlt6EP1BWyoHzc40S1SAYlYWSejNb3ZurIWFvLbf/IYDBkAJKS8vMeBP3mFSjchF9vhxV7VT9iF5fwFO0pT6KfqZbNtM+RQ184+SV1YW+/zkSRJVayuPVVjoAlKl+ol6WsrmpAFEn/4ilZgq0x496lRaivSvr3dSs6QRSb7yRFNgswYBnVv1i2prXSHHZBgbc4OOJGQ5sIgcAKeK8kMrqqqcdN/V5Qoq9vU5qUdtbm11FWj9LPcauOTQIMlkcNAF8F52GaVn5Qj0bSLJ1LDW8pwvv0z78NNPc12lOCoqXL0l1TxTG/xja1IldbacKdTfvn1Lx7CW7R4a4j0dGaGd7NprXRmMsjI3Jmhi6Ace+yqy3l4+I5K6/MKPImE9t34ohjxGZdeSU40y4EvqlDNRRYWLI1N5GpGxPCZ96VQ5ITWZkRSuhMIVFc4UIE/m/n7gE59YRRV1jTG3AfgsgEwAX7TWfnqm/dOZoHz099Pgu38/JRup6ZqaXEoXP/5JSSSVgbi01LmnKgeaHCwA55mlPFtSE/r2hljMBc2qENzEhNNnK25Ddq6GBs7CVeK5pITqywce4Dn37OF/WlvdDLmz05V8rqigMft1r+PgK1f8KJRl4MgRLk1N3F5WRrLauZMD+HzUebJHzERsCliVWqOz09lU5PVXWOgGXFVJlbQo12Zfty/VS7LaRgpK9tVQSiHU1+eCp1XXJzc3Md+hpOeSEvbH9ddzEK+r4z09coTSjILMq6u5X1UVr/XsWe6ne5WR4dIDKU2TZuCSeJV2SJmx5XkozzSRjJKwatCVnenyy6miVaJkpSGSykyZIgYHSQKqb6U2VlQ4G5QkVKkVRVyZmVTdHTvG+9ne7hyVXv1q4A1vcOpHP0WQ8v9NN/MXRGLyftXETfcmM9MNyhUVPH9vL0n2bW9zsUwLcRQYHXVZV86c4b1UqifdE9/71p9wivxEVn6ZHN8rULYtxcH5mTSA5HGNIlxJe6ojt2ULx46SEpo0FOR/sUgLgjLGZAI4CuAWAE0AngHwXmvtoen+s1IISrCWD9iBA3wR29v5sBw/7h42P7hPMyEg0cjqF0yTW65yd/nxEZK6lJlAgYa1tTzXyZOc9cnO1NjodOn+rEk2MGW4PnyY/9m5E7jpJpcVW5mxm5sZ5HzkCF+cq64CbriBv82GgQG26eWX2b5YjNcooty8eWa355nQ18cX/fhxHlu1merr+VJt3MhBRkHMGoQkkcmOcfass4HU11MFtmULiU0DkUhFg7AcEbTI/ihnidZWN1seH+d9q652UqckXU1Y/Ng13zbU1sZrPHvWeXoqBkqebFVVLlGooBm3T4gictnHfO82DVa6XnlZ1tfzHlVVzc+Rxc/eoPpNLS3sCzlglJW50AN50vm2rpoaDoivfS0HyPlmZlAMYJS8/HVlj9f+foJnqV3PneN9yMtz3qV+JnCl5pKrvTFT+17PjBIx+8QgVamvWvZzXKp9vgu+HxPoh4mcOeMkLKUGU+xXfb1Tx/q2MN/ZYmQk0cEmqqn53Oc4UVgI0oWgbgDwF9baW+PfPwUA1tq/me4/K42gfAwPc8a4fz9nfgqe1CAgLxjAzeyBRE86eQ1JR+6ndvFT92gwk1iuwUWBohLN5fKrB1HQ4OOrBmXDqKhgOhOVudYit9YjR0gGAAeNvXtJBP6+/uK7OY+OkkwOHyZpjY3xBdqxg2Q1WwxHLMYXUKSkxJ7Fxa6U95YtHEBU7tyXGPRd3ox+kLQSm8o9OSvLBW7v3EmpUSoS5QJUn7W2cvCVa7RedKlTNRGQKjgry1UerqlJVLfl5jpvLBGcDNfl5S5jh+wRsZgbLDXA+emnNFhq39FRZwvZtInX19jI8/hl2VUiQs4FKi++bp2z5fk57qZb97cND/P5OXSIfTI25vrUdwzIzma/1NY61Zs82fyltNStX2wCUw3uycjLL/Ny6hSdpsbG2BcTE86W6Gee8aUYP4RB75tf2iNZctcoogH+QjR7TZSwurvdBKyjIzFcYt06vi8qMulnY9E90GRG8WqacH3pS3QMWgjShaB+G8Bt1to/iH//AIDrrLV/Mt1/VjJBCda6FCVKuqnkiyqn7i9SNSn2ys+LJvWfZjKa7ejBV6YEP2uCJDY90LEYZ3fSK2vm6meuUOzU6dPU94+NUcJZty7R+0pGXwUfNjfzd5W9VsqoaH/oBVD2a3l59fY6A7B+v/RSqpEuvZRt7OhwhHTmjHPDVpLOigq+eD4RyYtRkLpUfaAMA378iv+pMAF5KyowOpo4Vg4L8vwqKkoMtFUBu6wsnrOy0qm/VBNKNgo/rk02A5GWVFLKKOKnhtKMXM4zksZFpNHME8opV1Y21f6zcSMnHDt3so+7uzkZOXkysVR4bS0nApJSo/axWIz9pz707W5Sl/b3O1tQTg6ftb17gVe+ksdUGY7pFv/+6t5ECcxf5ktgUgX75zx7Fnj4Yd67LVtc0lT1vSZjemf8dxVw6mRJXH6Arx8nl4zwfRKJZgGZbdvQEN/V48d5HyWJZ2ezDUr0W1fnwkmix9CxP/QhPpcLQboQ1DsB3BohqGuttR+N7Hc3gLsBoLGx8eozZ84sWxvTAZqJR4lLL7If5Cu1n2wKra0uf53cVRWr5AfJaqArLHQuwb5ko5mVBsIXX+RLeMklwB138MFVye6JiURj8/79XLq6+Lv2lTpMwcdSXfkJQhVDJpuapB7NXn3PJ1+lIvuKXn4tsjXJ3dovwjebGtFaV97ET76pEgwarEQCyvmWl+c84WQH9B0dNAApx6LuldReknR8DzM/A7u1yVVK8qwsK+O1d3VRKj12zHlVSiW1axdjzjZuTJyJ+58qKT8dJiedU4dsJyJVuSDLuC8bnKDSFVVVzruvqsq5qx86BDzxBO26g4P8vaGBi0+oWnJyZiaw6AQFmBuBjYywLQcPcrLmX4Mk3qIivh9tbbTbvuMdvJZkz5f6orPTOSToUxNTITc3MZO8/znfrOgzQQ4oBw9SItS1jI+7oohXX02b6ELU7zMhXQhqTan4FhPy/vMJyy9UBnBQqq3lC+yrCfr6XKLJjg4nomvAkl3JN7pKKpIEMTzsBmXl9xI5KsZEXkGSxBRMOT7OAWn3bgbkVlY69cb4eKKU09npbBOq+KkM5rIDyE4no69fBsJfNHgobsOP4YiuZ2Qkque07nsLlpQ4e1FJibNnKBaqqcklB5Wn2tVXOweU4WFnAJdDgzwJZUcSSYuU5GIcVQcrdkhxNnKCULLd8XGXMFYZtUtKONg3NztX+KIi57K+dWtipo7poHvmS0Fya1cmfMVW5ee71Ds7d/L+b9+ePBdcsmBr1SQS8SldmFypVYIlmbrPX+SQMpMEJql4ctIRhsrXV1e7+DY5cESlTU3MqquBN74xMRN4VBXnw88zKdWiMkH46kV/qC4p4T2trnbkLvKazrtzrpiY4LP89NNMJH3kiFMrV1byub72Wi5zsTnPBelCUFmgk8TNAJpBJ4nfsda+NN1/AkHNjPHx5NKWb18qLU10cMjJ4QP3wguu9PzIiJNKsrPdgOAb6OVYoCJyxcUudsYPIvQTtmZnOyeB0VGeY/t2DlTr1yfm9JNzgqRFDUZ+8lCpOZTNe3TUDSBy5a6ocMZ4P1OAVKOzQUmCZfOQw4IyaCtrNOBeWtlijCG5Hj7sCtn5rvyKjQHcACupStKfZsq+ZKF+ljqvt5fEeOAAB/TmZl6fXIj96qt+6Q5B7s5RD8P6ehLa1q38n2936O93KkRf4lZ8VlQiysoiCUvC6upy1y114IYN7FfFz/X2cp/1610GkNpa5z3a3Oy8MmXLNcbZxnzX8Og9nYnAcnI4ED/9NKU2xcqp/xSeMRs6OnjvMzJIaHI/Xyj8wFs5Wyi9mrLWSEOisibRpazs4qSf4WE+Z488ws+zZ52qurAQ+Pu/p416IUgLgoo35M0APgO6mX/JWvvXM+0fCGr+sJaDSZS0OjqSB8UqUFSuv0VFVOVdeSUHED/+RzPRn/+cBvSyMnpUKZHo0FDirFQzwP5+zq6bmlxJgMJCvjCytUkKq6hwnoZKIeTnIpM3lWpaSSpR0G9xscvYoTLVUvFJdaOKvCKuiQk3sPo1dVTSQfYDZQTQwJCb63LOtbS4OlfKv+dLR6oKq0qkO3a42CEt02WZ1qxWtp+mJucS39jI2XRDA7c1NbE/lDdwfJznVq45X30lb7KWFlfTK1nwtNLj+Pkj5Q7upxHyg22j6yMjvF+nT3MQl70pP5/366qraG/avXv6OkRCLEZJ21e9dnS4598vryHbm7KUy2lF91/GfnmlaQK1a5dTB+fnJybx1Xmi7x3AY33vezzPzTdTlRrFdP+dy3dJmi+/zEmQ1MiKXysoSEyPJijsIEpcFRVOQ+PHz2ldsWaCYi/b29l3X/4yJxsLQdoQ1HwRCGrxEIs5tYxft0lEMTnJwe/ZZzmATExwRv2KV3BAjRqVjx8H/vVf+RLccgtnUXPRi3d0AA89RBvD+LhLBKuyDdGYo2Sfit3woYwYLS3O4UA56pTQ088oLsKQVKnS1yqpojQ3fq5B2cYU1Kmie8qHmJdHUmxs5Eu7dSuvbWKCjh0K4FV8zaZNJKxLLuFAKMi+o5REZ886NZ/cvJNVNfahY0ilePasc9nOy2MbtSgPIMBrO3KEnxrk/YwhignSerJt0d81YCq1k4Kf/QSrgItRU3yf8jKqZtG6dS5mMArVRBJhyXMW4LmUmby4mM+J4qoU36REscqX6EOpwBoa2O+bN08f96e2/OAHvNd79gBveUvyAOSFQllSDh/mPZN0K5Xqhg18Ppubef/PneP7oYzxfs0qfwJSVcW+Xr+e76UkzfLyxOw38iZcKAJBBcwLQ0NU/x04wJc5O5szyj17nGEd4Atx3318ObZuZeDibLNfwS+uN18oYe50BNbd7WoTKbBYrvTyjvK9IWMx52KvgnLV1XxZpdrs73fZvOXmrzyF69e7mkKzDUSSchS4rJl/XR37UBkcNPutrnaD4qZNF+8+La9Rn7B07sxMtl+EFc0cP1+MjroYvKNHXcqjujp3fBXVVP9K8lOMkYLbFSwtSaysjM/gli0kjHXrXF7LKGR/PXGCSZ9ffNFdc0UFJ15797LflTnFGD4T0gBImujqYvvm6rloLT38HnmEx373uxce0DodlGvv8GGq7g8fdnk55QyifJeAk/7llq4wFb1DftCx7Lz+IrXzYjlMBIIKuChYyxnXgQMkrNFRvth79lAFWFLCffbvB372M5LNW9+68OzGi4nhYQ6Shw9T6ovFOOgpw7lsRwUFTr2nWfiFC06tU1TkiEjLxZKFj44O5xhw9iz7VIS0efPcCf9iMDjobH/KgKLrVUZ9qVdnWpfnYlMTVYVSQebluYTG27bNvWrr8DBJWtlYuroSEycrdZJc+ouLeR9VnmXrVn5vaSEpHT/O9lRV8feyMh5D91gSuQKR/XvstznquXj2bGLOSRFWfb0b5A8fBn74QxLYu95FgvYzQPiLn1l+pkVes1LDRSU+TaAGB12CAEnNV13Fd7emJjmp+9W7JfVqke0J4H/Ly4H3vMelLbtYBIIKWDDGx/myHTjAwcMY6uv37OEA1N0NfP/7HBSuvhq49daFexQtNsbGEgODFaRaW5votZeTM3WgUjn0pYQGu6U+z3RQglYFGPtSaVTFKjd8efINDfEYqvxaX89+9Us6zIfwopKoKgXLrtnZ6Wxtzc2U9JVNXOorhRps28bndOdOl47JzwQfDSVob3d2n/Jyl3VBCZT9Ehj6n2+/89NpKXPJM89wgN++feEDulSOvurNX5eHrdDV5epinTvnEkPv3MnJZEPD9KmaFAohL8gLF1z14vZ2ZjOfKRxhLggEFbCoUPmH557jTLSwkLOyK67gjPXXv+aL+Y53uODFdIOvv29ro/pJZFRVtfAibKsRw8NUlx0+zMFOue3Ud7W1HBx9UpuO5GYbevzchgo89+sV+eu9vRx45S04MuKkCOWyk2u+ijRq8QNlFUtWWOgGZoVlyMNQ8LNDaFEJGz+3odz5a2pIoiMjVC3eeKMj4qh9dKYlmo1lLpBKfHjYeRoePsyJpuIHlfhW4R8ipWT2Xh//7t/x/i8EgaAClgSTkxywDhygimpy0mV0UD6/172O+fzCgL/yYG2iG/jZsy7Dh9zAt2yZv6rTj7NLRl7+d4Uc+NnBYzFHSk1NzvFCwcXRYn6yIUoC6+1NrG7rlytRGimfABQ6IXJTILJCBFS5WtkjtD487CQNpXJSFvgNG4A3vYnSnZI567/+MZJtU6YZ3/U82bo+5YQSRSzmYhD7+pxjkIpSKvWRn4bJT6g7OkrJdDrv07kiEFTAkmNwkDEkzz7Ll9EYl3Ln8sspTS1WXEhAcigebaYceLPlx9P6hQskJcUvKQ/hjh2UllIx4RgYoITuF4nctInP165djijVD1FS87OXtLc7ddWFC07FK9uQVI7RVFUKIfCD2v3CilGVpdTc4+Ouxpa8M41x8XYqixFN8jtXKEzDr6GWbD26TeQ6Nsa+ff55xtgpfZRUh/Kw9V3tAeCP/mgVlduYLwJBrTxYy1nts89yIDl3jvr5+nrg935v4UF9KxkqROhnstaijNFSRc2VVKL52RYLWVl00hAp+W7wy4nRUaoTX3iBg7q1VCldcQVVZYvlGWctybilhc+rCkvKLpmZSfvRNdfwWfZjhqKfUbWYKgvLVlRezuN+73uUqOTmr3CFhgZKMRs2uKwl/n1WSjCfdGbyHlUyYb/0SHTxS5H4deu0PTfXFR+9/HK2saTEVeFeCAJBBaQEY2OcjT3+OPDgg3zY9+wBfvd3ObisBrWfMlv4RJNsXct0On2lXfJtK36Jivlsm+/+ybYp0e1i9pOCl30VVLIEvf56Xx/7rKyMA+MVVyzcKD9XWEvCURbw55/nva6ooHv6VVdN9Uy01gUFK+2TT169vW4iMT5O8pXzxPbtzkNPJVmqqly5l2TxV4rLm4l8pAL14Wdp8Ys6+ovc0i9ccE4Wra3cVldHB4u9e2eOCZsLAkEFpBxtbcDXvkZ39IwMevrdeCNjSBobly5G5GKgRLXTEY3/XWXZo1DclRY/MDpZRd6VQNa+sX0+JKOMGjNBThG+6klJcS+9lLP3VHk3CrEYbavPPEN7XFYWM19ccw3VnnNpn0qti7y6uoDHHqPDkdzys7NdPJ+CwtU/9fUkLaWhkiOGX3AwNzfRU1GOIPosL2ffXswzJ4/Aw4epKfnoR0nYC0EgqIC0wZkzwOc/T6cKwNXLkXuyAgGli58u0eZ0STfnuu5/V0Z2v3hfMqgWUjKSiX6PuvqmM8bHp9bJUkYFn2Rmy0enYptRkpnLtqhjQrqjtZXJVJ9/npqCujoS1eWXX1x4xYsvAvfey3649Vb2iSSvri5XKVupt/zKzdFlrlkrkjljzGd9fJxZMhaq/g0EFZBWGB0FfvlLGrpV1bejw7nmAi5xqhZ5ZkXzks22Ppf9VBtpNklnMVVey42xsUQC8helBBKUZ28+ZDObHWS1YnSUttZnniFp5eYy5OKaa+avirxwAfj2tzk5uP12qsN9WJuYLHYuBSKj6xfzn5nW3/e+hWc1DwQVsCIQTcdz5gzdnAEOfg0Niel4FureutowOjo9CfX3J+5bWJhYY8hPqruQNEdrFXIOeuYZ2l0nJqi+vuYaV/BxLhgaovPEyZMsaXHrrUtTgymdEAgqYMVicDAxf1xLC2duxtA9t7HR2bGWMi1QumBkxBW6iy5+KhqA/eGTj78Ecl86DA7SprRvHydchYVMunz11XMLtZicBH7xC+A3v+Gz/a538RirFYGgAlYNxsZcYlGlulEgYnl5YpbuqqqVZdcAnBpnOklIKYWEkpKp5FNZyb5It1RTaw3WukS1R49ym1zVt26d3Unh4EEmY87PZ867+vqlb3MqEAgqYNViYoK6f6kE/YJqBQWJhLVu3fKqS5SvLerVNtP3aPJPY1wRuigRlZevbLvYWkJvL5MqHzhAL9DyckpUe/bMLB21tADf+Q7/c8cdtG+tNgSCClgzUMClrxaUHSs7m+7AUgk2NMys6lJg7XwIxv8+m9ebsnH7jgZlZVNJaC06H6xWTEzQRXvfPubCy8xkFoxrrpnelX5oCPjud+nFd/31rL+2muxSgaAC1jQGBqbasVRsra6O6WaUG84nG6W0mQ7KXeYv0ezdybbNJfo/YPWjvZ1E9dxzfNZqaxn4esUVUydOk5OuknVFBSdYyrJeW7uyn6VAUAEBHkZH6eLul4tXMtD5kE1u7sqzcQWkH8bGaG/at4+Tp5wcktQ115B8fLz0Egnt/Hmnys7M5CRLhKWCmytFygoEFRAQEJDmsJYTp337GLgbi1FS2ruXakBfSlI9LtWi0iL7ZVaWK4Mi0qqsTM+MJYGgAgICAlYQhoacq3pXFx1+9uwhWZWXJ/+P4gh90mppcV6uKsTpL+XlqdcCBIIKCAgIWIGwlkG7+/YxPZjqXs2VVJS8VnWwVJFYZTOysxNz9qn+01yP/973Ll0miRVsVgsICAhY/TCGMVNbt5JYDhxgPr65/E+I2rEmJ105DRUsbGtjqiWABCUPUnmTTleUcimdMwJBBQQEBKwQlJQAr33t0hx7fJzxhL49q73dBYiXlk5VDy51SqxlIyhjzF8A+EMA7fFNf2at/elynT8gICAgYHpkZzMusKHBbRsbcwUctRw+7H6vqADe/e6pEtpiYbklqP/PWvu3y3zOgICAgICLQE4Og9o3bnTbhocTSWsp67gFFV9AQEBAwJyRn8/qvlu2LP25ltsj/k+MMS8YY75kjJnGURIwxtxtjNlnjNnX3t4+3W4BAQEBAasYi+pmboz5BYC6JD/9ZwBPAugAYAH8JYB11toPzXbM4GYeEBAQsLqRVnFQxphNAH5srb1sDvu2AzizwFNWgeQY4BD6JBGhP6Yi9MlUhD5JxGL1x0Zr7ZT6w8vpxbfOWtsS//p2AC/O5X/JGn0R596XjJ3XMkKfJCL0x1SEPpmK0CeJWOr+WE4nif9hjLkKVPGdBvBHy3jugICAgIAVhmUjKGvtB5brXAEBAQEBKx9pmNd2SfD5VDcgDRH6JBGhP6Yi9MlUhD5JxJL2R9oniw0ICAgIWJtYKxJUQEBAQMAKQyCogICAgIC0xKomKGPMbcaYI8aY48aYT6a6PamGMWaDMeZXxpjDxpiXjDEfT3Wb0gXGmExjzLPGmB+nui3pAGNMmTHme8aYl+PPyw2pblMqYYz5D/F35kVjzLeMMdMUn1i9iGcAajPGvOhtqzDGPGiMORb/nDZD0MVg1RKUMSYTwOcAvAnALgDvNcbsSm2rUo4YgE9Ya3cCuB7AR0Kf/Bs+DuDwrHutHXwWwAPW2ksBXIk13DfGmPUAPgZgbzy5QCaA96S2VSnBvwC4LbLtkwAestZuB/BQ/PuiYdUSFIBrARy31p601o4B+DaAO1PcppTCWttirT0QX+8HB531qW1V6mGMaQDwFgBfTHVb0gHGmBIArwbwzwBgrR2z1vaktFGpRxaAfGNMFoACAOdT3J5lh7X2UQBdkc13AvhKfP0rAN62mOdczQS1HsA573sTwmD8b4inm9oD4KkUNyUd8BkAfwpgMsXtSBdsAeu2fTmu9vyiMaYw1Y1KFay1zQD+FsBZAC0Aeq21P09tq9IGtcoQFP+sWcyDr2aCMkm2BZ96AMaYIgDfB/DvrbV9qW5PKmGMuR1Am7V2f6rbkkbIAvAKAP9ord0DYBCLrLpZSYjbVe4EsBlAPYBCY8z7U9uqtYHVTFBNADZ43xuwBsXyKIwx2SA5fcNa+4NUtycNcBOAtxpjToNq4NcbY76e2ialHE0Amqy1kq6/BxLWWsUbAJyy1rZba8cB/ADAjSluU7qg1RizDmC+VQBti3nw1UxQzwDYbozZbIzJAY2a96W4TSmFMcaAdoXD1tr/ler2pAOstZ+y1jZYazeBz8gvrbVrenZsrb0A4Jwx5pL4ppsBHEphk1KNswCuN8YUxN+hm7GGnUYiuA/AXfH1uwDcu5gHX7UVda21MWPMnwD4Geh18yVr7UspblaqcROADwA4aIx5Lr7tz6y1P01dkwLSFB8F8I345O4kgN9LcXtSBmvtU8aY7wE4AHrCPos1mPLIGPMtAK8FUGWMaQLwXwB8GsA9xpjfB4n8nYt6zpDqKCAgICAgHbGaVXwBAQEBASsYgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAgICAtISgaACAtIYxhhrjNmW6nYEBKQCgaACAtIExpiHjTF/kOp2BASkCwJBBQQEBASkJQJBBQQsMowxp40xnzLGHDLGdBtjvmyMyTPGlBtjfmyMaY9v/7ExpiH+n78G8CoA/2CMGTDG/IN3yDcYY47F//O5eFVXGGO2GWMeMcb0GmM6jDHfScHlBgQsGQJBBQQsDd4H4FYAWwHsAPDn4Pv2ZQAbATQCGAbwDwBgrf3PAB4D8CfW2iJr7Z94x7odwDUArgTwrvhxAeAvAfwcQDmABgB/v7SXFBCwvAgEFRCwNPgHa+05a20XgL8G8F5rbae19vvW2iFrbX98+2vmcKxPW2t7rLVnAfwKwFXx7eMg2dVba0estY8vwXUEBKQMgaACApYG57z1MwDqjTEFxph/MsacMcb0AXgUQJkxJnOWY13w1ocAFMXX/xSAAfC0MeYlY8yHFqvxAQHpgEBQAQFLgw3eeiOA8wA+AeASANdZa0sAvDr+u4l/2vmcwFp7wVr7h9baegB/BOB/B5f0gNWEQFABAUuDjxhjGowxFQD+DMB3ABSDdqee+Pb/EvlPK4Atcz2BMeadcrIA0A0S3MSCWx4QkCYIBBUQsDT4JujAcDK+/BWAzwDIB9AB4EkAD0T+81kAvx331vu7OZzjGgBPGWMGANwH4OPW2lOL0/yAgNTDWDsvrUJAQMAsMMacBvAH1tpfpLotAQErGUGCCggICAhIS8yJoIwxZcaY7xljXjbGHDbG3GCMqTDGPBgPIHzQGFPu7f8pY8xxY8wRY8yt3varjTEH47/9nQIOAwICAgICopirBPVZAA9Yay8FgwUPA/gkgIestdsBPBT/DmPMLgDvAbAbwG2gZ5HcaP8RwN0AtseX2xbpOgIC0gbW2k1BvRcQsHDMSlDGGLnD/jMAWGvHrLU9AO4E8JX4bl8B8Lb4+p0Avm2tHY0bbI8DuNYYsw5AibX2CUvD11e9/wQEBAQEBCQgaw77bAHQDuDLxpgrAewH8HEAtdbaFgCw1rYYY2ri+68HPZSEpvi28fh6dPuMqKqqsps2bZpDMwMCAgICViL279/fYa2tjm6fC0FlAXgFgI9aa58yxnwWcXXeNEhmV7IzbJ96AGPuBlWBaGxsxL59++bQzICAgICAlQhjzJlk2+dig2oC0GStfSr+/XsgYbXG1XaIf7Z5+/tR9A1gFH1TfD26fQqstZ+31u611u6trp5CqgEBAQEBawCzEpS19gKAc8aYS+KbbgZwCAwMvCu+7S4A98bX7wPwHmNMrjFmM+gM8XRcHdhvjLk+7r33Qe8/AQEBAQFzQEcHcGqNhGPPRcUHAB8F8A1jTA4YFf97ILndY4z5fQBnAbwTAKy1Lxlj7gFJLAbgI9ZapV/5MIB/AaPp748vAQEBAQFzwOQk8K1vAZ2dwKWXAm96E1BamupWLR3SPpPE3r17bbBBBQQEBACHDgH33ANcdhlw5AhgDPD61wPXXgtkrOC0C8aY/dbavdHtc5WgAgICAgJSCGuBxx8HKiuBd7wD6O0FfvIT4IEHgBdeAG6/HaivT3UrFxcrmHMDAgIC1g5OnQLOnwduuonSUnk58L73Ae98J9DXB3zhCySr0dFUt3TxECSogICAgBWAxx4DiouBK65w24wBdu8Gtm4FHnoIeOopqgHf/GbaqFY6ggQVEBAQkOZobqYEdcMNQFYSsSIvD3jLW4APfQjIzwe+/W0uvb3L39bFRCCogICAgDTH44+TeK6+eub9NmwA7r4buOUW4MQJ4HOfA558kt5/KxGBoAICAgLSGO3twOHD9NTLzZ19/8xM2qn++I+BjRtpl/rCF2i/WixYCxw9Cnzta0tr8wo2qICAgIA0xq9/DWRnA9ddN7//lZcDv/M7tEndfz9J6rrrgNe9bm5ElwzWAseOAQ8/TMIrKwO6u4G6uos73mwIBBUQEBCQpujtpQv5NdcABQXz//9iOVEkI6a3vhW48kpKbEuFQFABAQEBaYrf/IafN964sOPIieLKK4Ef/YgOFHPJRJEqYhICQQUEBASkIYaGgAMH6Fa+WOmMGhroRPHkkySdz30ueSYKa4Hjx7lPc/PyE5MQCCogICAgDfHUU0AsRoeHxYScKHbtAn76UzpRPP88cMcdwLp16UFMQiCogICAgDTD6Cjw9NNUwy1VxSHfieKnPwU+/WnarIqLgaqq1BKTENzMAwICAtIM+/cDw8PAK1+59OfKyQEKC1nG48gRqhZvuQV4xStSS05AkKACAgIC0gqxGPDEE8DmzcD69UtzjmQ2pv/0n5iI9qc/Bb77XXoPvvnNqS3nEQgqICAgII3wwgtAfz/w9rcv/rGtZYaJhx8GmpqS25iiThSvex3jp1JRziMQVEBAQECaYHKSgbn19ZSgFgvJiOmOO4CrrpqqxpMTxe7dLOfxs5+RNO+4Y/nLeQSCCggICEgTHD7MarnvehcdFhaK+RBTFGVldKI4fNhlorj2WrqlX2wmivkiEFRAQEBAGsAvSLjQUhkLISYfxtAdfcsWZqJ4+mkS1pvexDYuBonOhEBQAQEBAWmAEyeAlhbgzjsv3t4TJabS0osjpiiimSi+8x3gkkuW3okiEFRAQEBAGuDxx4GSksSChHPFUhFTFMpE8dRTwK9+RSeKu+5aOm/DQFABAQEBKUZTE3D6NHDrrfMjFGuBkydJTOfOLR0x+cjMZG7AXbuYK3CpMpkD8yAoY0wmgH0Amq21txtjKgB8B8AmAKcBvMta2x3f91MAfh/ABICPWWt/Ft9+NYB/AZAP4KcAPm6ttYt1MQEBAQErEXMtSCikgpiiKCujim8pMR8J6uMADgMoiX//JICHrLWfNsZ8Mv79/zLG7ALwHgC7AdQD+IUxZoe1dgLAPwK4G8CTIEHdBuD+RbmSgICAgBWItjbg5ZeB176WWR1mQ3Mz8+eJmG6/HdizJ/VZH5YCczLFGWMaALwFwBe9zXcC+Ep8/SsA3uZt/7a1dtRaewrAcQDXGmPWASix1j4Rl5q+6v0nICAgYAqefx74/OeBnp5Ut2Tp8Otfk5iuvXb2fXt7ga9/nZ+33w587GPA3r2rk5yAuefi+wyAPwXgV7avtda2AED8sya+fT2Ac95+TfFt6+Pr0e0BAQEBU/Dii8C//ivrEN1zD1MArTb09AAHD1K1N1tBwokJpiCanAR+93dXNzEJsxKUMeZ2AG3W2v1zPGYyz3g7w/Zk57zbGLPPGLOvvb19jqcNCAhYLXj5ZeAHPwAaG4Hf+i2S1AMPpLpVi4/f/IaxRDfcMPu+v/gFnSnuvBOoqFj6tqUD5mKDugnAW40xbwaQB6DEGPN1AK3GmHXW2pa4+q4tvn8TgA3e/xsAnI9vb0iyfQqstZ8H8HkA2Lt3b3CiCAhYQzh+nJJCfT0zGeTmAhcuUBXW2HhxbtjpiMFBV5CwpGTmfQ8fZgLZ666j99xawawSlLX2U9baBmvtJtD54ZfW2vcDuA/AXfHd7gJwb3z9PgDvMcbkGmM2A9gO4Om4GrDfGHO9McYA+KD3n4CAgACcPs1y5NXVwPve51Lq3HwzsHEjg0Tb2mY8xIrBU09RbTdbQcKuLqo6168H3vjGZWla2mAh+Wk/DeAWY8wxALfEv8Na+xKAewAcAvAAgI/EPfgA4MOgo8VxACcQPPgCAgLiOHcO+OY3WUjvAx+g27WQkQH89m+TsO65hwX9VjJUkHDnThYHnA6xGKXJjAzgne9c/TanKOZFUNbah621t8fXO621N1trt8c/u7z9/tpau9Vae4m19n5v+z5r7WXx3/4kxEAFBAQAtDF9/etAURHwwQ+ygF4UxcUkqc5O4L77GAu0UrFvHzAyMntBwgceYPqjt7+dcUdrDaGibkBAQErR1gZ87WvM93bXXSSi6bBpE9V9L71ECWQlQgUJt2yZuXzFwYMksle+EtixY/nal04IBBUQEJAydHYCX/0qkJVFcppL4tGbbmKi0p/9jGrBlYbnnwcGBoBXvWr6fdrbaW/buJHlLdYqAkEFBASkBN3dwFe+QlXdBz84d9dpY6jyKi2lfWZwcGnbuZhQQcL16ykNJsPYGK8rO5sqzVRUsk0XrOFLDwgISBX6+ig5jY/TIaK6en7/z8tjUb+hIcZLTU7O/p90wKFD9Mp75SuT11KyllVs29sZ/zWTunMtIBBUQEDAsmJggJLT0BDw/vdffDbsdeuYrPTECeCRRxa3jUsBFSSsqpq+IOGzz1IF+NrX0ka11hEIKiAgYNkwNETJqa+PcU4LrSO0Zw8zeD/6KAN80xnHjzPgeDrp6cIF4Kc/BbZundk+tZYQCCogIGBZMDJCb72uLuC972VWiIXCGFZ6rakBvv99JlFNVzz+OO1ml18+9bfRUdqd8vOBd7xjbdudfIRuCAgIWHKMjQHf+AZdyt/97sVVX2Vn0x41OZm+SWXPngXOnGGhv2iwrbWM6+ruZjBushiwtYpAUAEBAUuK8XHgW99iotPf+i1g+/bFP0dlJfC2t7FW0s9/vvjHXygef5zZyvfsmfrbM88wruvmmxdHqlxNCAQVEBCwZIjFgO98hzn23v72pU10unMns4I//TSDXNMFra3A0aNM9BotSNjczHiuSy6hdBWQiEBQAQEBS4KJCdqFjh9ncb3lyEL+hjdQCvnRj+iqnQ6YriDh8DDtTsXFlP6SOU6sdQSCCggIWHRMTjID9+HDwJvexIJ8y4HMTNpxsrNpjxobW57zTofubhZe3Ls3MfmttcAPfwj097O9/m8BDoGgAgICFhXWUoI5eJASzXXXLe/5lVS2oyP1SWWnK0j4m99Q7XfrrQt3tV/NCAQVEBCwaLAWuP9+Bpy+5jWzZ+teKmzezBx2L75IJ4RUYGCA/XDllYkZIc6cAR56CNi9G7jmmtS0baUgEFRAQMCiwFrgwQfppHDjjcyGkEooC/jPfkYPwuVGsoKEg4PA977HmldvfWuwO82GQFABAQGLgkceoerqmmuAW25J/eCrpLLFxXRGGBpavnOPjJCod+2iCzxAu9z3v0/niHe+01ULDpgegaACAgIWjF//Gnj4YaYdevObU09OQn4+g3gHBphUdrnsUfv2MTuEr+J89FHg5En2z8XmH1xrCAQVEBCwIDz9NFV7l12Wnmqr+np6Eh4/TpJYaoyPA08+CWzbxoS2AInpkUdI4MmCdQOSIxBUQEDARePAASY4vfRSqtPSNYfc1VczDuvhh5n9fCnx3HOU2CQ99fVRtVddnV7S5UpAmj5OAQEB6Y6DB+lOvm0b3bqjOebSCcYwWLi6emmTyk5O0g7X0MBquJOTdIoYH6fdKZpJImBmBIIKCAiYNw4fZqDpxo1M/pqVleoWzY6cHNqjYjE6TUxMLP45XnqJwbmvehVJ8Ze/ZKLYO+6Yf1HGgDkQlDFmgzHmV8aYw8aYl4wxH49vrzDGPGiMORb/LPf+8yljzHFjzBFjzK3e9quNMQfjv/2dMUHYDQhYaTh2jFLB+vUsm5GdneoWzR1VVcCdd9Lt/MEHF/fYKkhYXU339qNH+X3v3uQlNgJmx1wkqBiAT1hrdwK4HsBHjDG7AHwSwEPW2u0AHop/R/y39wDYDeA2AP/bGCPh/x8B3A1ge3y5bRGvJSAgYIlx6hSTv9bUsODgSnSV3r0buP56OjK89NLiHffYMSaGfeUrqUL84Q/pJHFbGOUuGrMSlLW2xVp7IL7eD+AwgPUA7gTwlfhuXwHwtvj6nQC+ba0dtdaeAnAcwLXGmHUASqy1T1hrLYCvev8JCAhIc5w9y7IZFRXABz4A5OWlukUXj1tuATZsAO69lymRFgOPPw6UlTGr+ne/S4nqne9cGerPdMW8bFDGmE0A9gB4CkCttbYFIIkBqInvth7AOe9vTfFt6+Pr0e0BAQFpjvPnWXCwuBj44AdZ22glQ0lls7IWJ6ns2bNcbryRaYyam6lKrKhYnPauVcyZoIwxRQC+D+DfW2v7Zto1yTY7w/Zk57rbGLPPGLOvPV1y5gcErFG0trJUe34+yamoKNUtWhyUlLCAYns78OMfLyyI97HHWAk3J4cpjm64gZJUwMIwJ4IyxmSD5PQNa+0P4ptb42o7xD/b4tubAGzw/t4A4Hx8e0OS7VNgrf28tXavtXZvdXB9CQhIGTo6gK9+lY4Qd90FlJamukWLi61bmTPwhReA/fsv7hgXLtD+dOmlTJTb0MAs7gELx1y8+AyAfwZw2Fr7v7yf7gNwV3z9LgD3etvfY4zJNcZsBp0hno6rAfuNMdfHj/lB7z8BAQFphq4u4CtxK/MHP8gEp6sRr341Y7nuv5+qzPni17+mqvDUKac6TOeYsJWEuUhQNwH4AIDXG2Oeiy9vBvBpALcYY44BuCX+HdbalwDcA+AQgAcAfMRaq4iDDwP4Iug4cQLA/Yt5MQEBAYuD3l5KTrEYyamqKtUtWjoYA7zjHVRd3nMPk7nOFV1dLOkxOcn1t7999UmZqYSxqazmNQfs3bvX7tu3L9XNCAhYE7CWrtcPPMDsB3fdxVx2awHNzcCXvgRs2QL8zu/MLSXRj3/Mviouplrv5puXvp2rEcaY/dbavdHtIZNEQEAAAKCnB/jmNxmEW1ICfOhDa4ecAAYe33or7UmPPTb7/gMDdC3v7QUuuQR43euWvo1rDcFDPyBgjWNykkGrv/oVpYbbbgOuvTZ9E78uJa65Bjh3jn3R0EBpajo8+iidK66/nt6Aa7G/lhqBoAIC1jDOn2fC15YWpud5y1vWtg3FGObNu3CBSWX/6I8oTUYxPMy4sKIiBi37Jd0DFg+B8wMC1iDGxlgK/QtfoKrqXe9iXr21TE6CksqOj1PdmSyp7Ne/TpvVe98LbN68/G1cKwgEFRCwxnD0KPC5zwFPPME6SR/5CEuTh9TNDtXVLL549izwi18k/nb2LIlr507gbW9LSfPSAiMjfJYmJ5fuHEHFFxCwRjAwwFifl17iAPyhDwGNjaluVfrisstIRk88wbx9u3ZxUP7MZ/j7xz++tkh9bIz9cfo0Y77On6fX5913L50zTSCogIBVDmtZ+fbBBxnX9PrXAzfdFIJJ54I3vpGqvHvvBWprgZ//nLWw3vAGeu6tZsRiLEty6hSX5maqOzMz6fH46ldTvVlTM/uxLhaBoAICVjHa2+kEcfYssGkTHQAqK1PdqpWDrCzao/7pn4AvfpHSQ309E8GuNulpYoJSkQjp3DmSlDG85htuICFt2LB8lYEDQQUErELEYozReewxDiZ33glcddXqG1SXA6WlzDTx9a8Do6O0223fnupWLRyTk/RWPHWKxHvmjMvqXldHl/tNm1g1OVWlVQJBBQSsMpw+zQwHHR3AFVcw+LSwMNWtWtnYto3xYT/9qSvnvtJgLdDW5mxIp0/TpgbQJnnllZSQNm1Kn3IqgaACAlYJhodpZzpwgIld3/9+DqwBC4O1lDSef57q0csuS3WL5gZrmR9QKrvTp4HBQf5WXk6nDxFSusZxBYIKCFjhUP68++8nSd10E0tIZGenumUrF+PjwMmTdKM+ehTo73dJZdM5Y0RPjyOkU6fYboDBxtu2OUIqK0thI+eBQFABASsYPT3AT37C/HHr1zOrQV1dqlu1MtHby348epTkFIsBubmsGbVjB+1OqVKVTk7So+7YMeD4carm6upoHxsb48SkpYXPA8B2iow2b2Zl35WolgwEFRCwAhHNn/emN9Gonc6z+3SDtXSdlpR04QK3l5cDe/eSlDZuTJ07/sAAyejYMeDECZLS+DgrGw8O0glGBcezskhY27ZRBXnppfS8Ky9fmcQkBIIKCFhh8PPnXXIJ8OY3hxRFc8XoaKLqbnCQA3hjI3DLLSSlqqrUDOqTkyRMkVJzMyWjsTGSEkCyjMVoM3rb2yg1FxZyW2srSfbQIdaoAigB1tUB69a5z6qqlRMDFwgqIACcTaf7THNsjBLTk08ySem73sV0O+ne7lSjpwc4coSEdPo0433y8qiy27GDUocIYLkxOEhCOn6c7WtvB/r6KBFlZLBdpaW8342Nbqmrm55kYjEep6WFy4ULLGc/Ps7fs7IYXOuTVm1tetosA0EFrHrEYnzp+/poZ0i2PjTEASEzM3FJtm267Yuxb2Ym7QXFxYnEc/QobU29vVTl3Xxz6mJT0h2y10hKamvj9qoq4LrrSEobNqRGirCWEvCxY8xI8fLLvKcjIySOkhKSRU0N1YsipMrKuU9EsrJIOuvWuW2Tk0BnJ8lKpHXoEIkL4LGrqhJJq65uKnGPjDB8oaODJNjRwZyFS2WbCwQVsKIxMUFPJZFNMgKSa60PzUxLSlj3p7CQL/HEBBd/3V/87WNjc9v/YpJp5udzoCop4WDW1sYBa6Xmz7OWE4WRES7Dw2492Xd/WyxGMs7PZ3yO/6n1jAyquM6d4zI6ym0bNzIObMeO1GXQGBqiDengQRJCayufy4wMetNVVdERY9MmR0iL7fadkcFYp+pq4PLLuc1aviM+aZ0+TXf6sTE3acvN5efkJIksJ4efmkwNDgaCCliDmJwk+fhkEyWgwUG+aD7y8jiwl5bSUKz1khK3LFeqFl1HlMCmI7RYjLPSCxeAZ54Bnn6aqpmNG6nmue8+Epe/lJYuj5pvcpID/3zIxf+erGyFj5wc3jstpaW8vqwsdwwNqENDXO/s5NLby+cgO5uDZm0t7TPDw5Sizp2bntz0mZe3OE4mkpKefZb38OhRtm98nG2rqQFe+UqqGDdu5AQpFdKwtbwnxrDf8vNJjP39fK8GBriuYN7CQvZTZSULOW7b5t6vpUIgqICUYHKSL8BMqrf+/qnkk5PjyKamxq37BJSbm5prmg4ZGVyy5vi2FRfTyJ2ZCbznPcxcMD7OmXdrKwe/l15y++fmTiWtmprZ+yEW4z0YGHADkv99aCiRbEZHZ79OkUt+viMZ/7u/RLfNpnKbmCDRHD1Km9LoKK9z+3YO8nK7Hh1l24eH3WdbGz+Hh2eWaPPypiew/PzpCWxoiCq7F1+kPam7m/esoIAS0mWXsZ1yUtCz0NPDJfqcz4asLB5bpFFQwPudbKIyPj5VLdfRQWL3Jw3FxWzb3r2UtKqquBQXU6JqbU20az3xBPvyj/946QjW2Pn2zDJj7969dt++faluRsAcYS0HCEk+/f2J6/ocGJj6UmZnTyWbKAFN9xKuBExMsG9GRxnh39rKAaOri0tPD783NfE6r76aXnoaHDUY5eWx73p6OBB2dXGwaW/nQBKLObVYcbEjgpwc/m9wkPdhYICDV3QRoebm8nxFRVyKi7n490ZLQQHv32LcG2tJiCKXri6SkuJ/MjMZ26PYpPLy+R1bUuDQkFMP9/SwT3p63PPpSxLDw+xT9Y0x3F/3Tc9zZib7qrqahFldzT6T04P+uxTwVcqxmHvWRkd5b3JyuFRVUfKpryfBi4zmSzJyxqitXbjkaYzZb63dO2X7chOUMeY2AJ8FkAngi9baT8+0fyCo9IHsPRpQ/YG1p8dJPiMjTl0l1VVGBl+SrKzElzU/n4NdebkbTP2Z9UzrWVnTv+waiDTzHx52g40+h4a4DA5yYM/MTGxjdnbiiye1nAhgYoLHFfH6x9TgqvOPjnLw13+tdQQd/ZTkoevTNVrrBqBk2/w+j8US1YiAc86QKk2kV1jI9ZwcXrMIMTPTkRaQOMD665mZiRKHlsJCHlPSguxQ4+Psb0k6vtQsCW5szO03MUGyrKx0s3pdi9RU4+NTVaX+ot/977of6kPdE61L0tIz4KucBwddv+blOVtSZSWv3Xd60fOu58qfMKjPRSDa1+9vIPE5kPahrc1NTrq72Ze6pslJHkvn1z33J3giVJ1H99P/XedTf/j3T334t3/LtEkLwXQEtawqPmNMJoDPAbgFQBOAZ4wx91lrDy1nO5YCelFmsjPogdcD4H8uZBvAY4+NcUDUw6PAPg0GMlL7A2Us5iLR/UE2quYZGnKzc12vzq8XQC+DXozcXC5S38hzTYOaMW6AP3mS7Rodde0dHU1uo/H7Vm3xBxf1R7K5l//C67vW1S/J7qN/Dn8Q0/mjRKlt/uKfKzrA+4PE+DhVKBoQdE3JSC16bdFnRX0OuP4bHSUh+G3yB0kNmpmZ7jrUD+of9ZU/8OueiDC1JGtztP3Jnmu/X/U9ivnOr/1jzvRf//n2P7Xue2GqP8+eTfxf9H5F16PXNN2EZab2+W2a7jrneuyLlVWuvBL4m7+5uP/OhuW2QV0L4Li19iQAGGO+DeBOAEtGUCtVHRQQsJYwl0F5sc81n32j/xFJBwBf+crSEdRyJ0ZZD+Cc970pvi0Bxpi7jTH7jDH72pXLIyAgICAg7dDdvXTHXm4JKpk8M2U+Y639PIDPA7RBLeSElZWuKqRvT5BqRUimf00mfUX3mU7NM90sLdkx/RgD6dWjx/HVAsm2+9t8L6XosfzjRNUW/nqy9kevb7rrjf7XP2dUBRdV5Ux3vuixk11XsnMmw3Sz4pmk7flI4ostBUgdOlO2C53zYo3VMz2vUdVyVP2V7FhqR1TlFT32dOeO3r+59uls93C+z0SyY8/lWVhuzU1UtevbKYGp6lJ9zvQeJVNDaruP//7fF+cakmG5CaoJwAbvewOA80t5wo6OpTx6QEBAQMBSYblVfM8A2G6M2WyMyQHwHgD3LXMbAgICAgJWAJZVgrLWxowxfwLgZ6Cb+ZestS/N8reAgICAgDWItA/UNca0AzizwMNUAQjKvkSEPklE6I+pCH0yFaFPErFY/bHRWlsd3Zj2BLUYMMbsSxYEtpYR+iQRoT+mIvTJVIQ+ScRS90eovxkQEBAQkJYIBBUQEBAQkJZYKwT1+VQ3IA0R+iQRoT+mIvTJVIQ+ScSS9seasEEFBAQEBKw8rBUJKiAgICBghSEQVEBAQEBAWmJVE5Qx5jZjzBFjzHFjzCdT3Z5UwxizwRjzK2PMYWPMS8aYj6e6TekCY0ymMeZZY8yPU92WdIAxpswY8z1jzMvx5+WGVLcplTDG/If4O/OiMeZbxpgUFGlPLYwxXzLGtBljXvS2VRhjHjTGHIt/zqN85OxYtQTl1Z56E4BdAN5rjFlgWa0VjxiAT1hrdwK4HsBHQp/8Gz4O4HCqG5FG+CyAB6y1lwK4Emu4b4wx6wF8DMBea+1lYBac96S2VSnBvwC4LbLtkwAestZuB/BQ/PuiYdUSFLzaU9baMQCqPbVmYa1tsdYeiK/3g4POlHInaw3GmAYAbwHwxVS3JR1gjCkB8GoA/wwA1toxa21PShuVemQByDfGZAEowBInuU5HWGsfBdAV2XwngK/E178C4G2Lec7VTFBzqj21VmGM2QRgD4CnUtyUdMBnAPwpgMlZ9lsr2AKgHcCX42rPLxpjClPdqFTBWtsM4G8BnAXQAqDXWvvz1LYqbVBrrW0BOAEGULOYB1/NBDWn2lNrEcaYIgDfB/DvrbV9qW5PKmGMuR1Am7V2f6rbkkbIAvAKAP9ord0DYBCLrLpZSYjbVe4EsBlAPYBCY8z7U9uqtYHVTFDLXntqJcAYkw2S0zestT9IdXvSADcBeKsx5jSoBn69MebrqW1SytEEoMlaK+n6eyBhrVW8AcApa227tXYcwA8A3JjiNqULWo0x6wAg/tm2mAdfzQQVak9FYIwxoF3hsLX2f6W6PekAa+2nrLUN1tpN4DPyS2vtmp4dW2svADhnjLkkvulmAIdS2KRU4yyA640xBfF36GasYaeRCO4DcFd8/S4A9y7mwZe7ou6yIdSeSoqbAHwAwEFjzHPxbX9mrf1p6poUkKb4KIBvxCd3JwH8XorbkzJYa58yxnwPwAHQE/ZZrMGUR8aYbwF4LYAqY0wTgP8C4NMA7jHG/D5I5O9c1HOGVEcBAQEBAemI1aziCwgICAhYwQgEFRAQEBCQlggEFRAQEBCQlggEFRAQEBCQlggEFRAQEBCQlggEFRAQEBCQlggEFRAQEBCQlvj/AUl1km+2GmlEAAAAAElFTkSuQmCC\n",
                        "text/plain": [
                            "<Figure size 432x288 with 2 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "fig, axes = plt.subplots(2, 1)\n",
                "\n",
                "titles = 'log paths', 'paths'\n",
                "transforms = np.log, lambda x: x\n",
                "for ax, transform, title in zip(axes, transforms, titles):\n",
                "    for i in range(50):\n",
                "        path = simulate_asset_price_path(default_model)\n",
                "        ax.plot(transform(path), 'b-', alpha=0.5)\n",
                "    ax.set_title(title)\n",
                "    \n",
                "fig.tight_layout()\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "723b16c7",
            "metadata": {
                "id": "723b16c7"
            },
            "source": [
                "### Monte Carlo Computation via Numba"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "b3e57768",
            "metadata": {
                "id": "b3e57768"
            },
            "source": [
                "Now we compute the price of the option $P$ by Monte Carlo, averaging over realizations $S_n^1, \\ldots, S_n^M$ of $S_n$ and appealing to the law of large numbers:\n",
                "\n",
                "$$\n",
                "\\begin{aligned}\n",
                "P & = \\beta^n \\mathbb E \\max\\{ S_n - K, 0 \\} \n",
                "  \\\\\n",
                "  &\\approx\n",
                "    \\beta^n \\frac{1}{M} \\sum_{m=1}^M \\max \\{S_n^m - K, 0 \\}\n",
                "\\end{aligned}\n",
                "$$\n",
                "    \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 11,
            "id": "153ec2db",
            "metadata": {
                "id": "153ec2db"
            },
            "outputs": [],
            "source": [
                "from numba import njit\n",
                "default_M = 20_000_000"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 12,
            "id": "3da12c7a",
            "metadata": {
                "id": "3da12c7a"
            },
            "outputs": [],
            "source": [
                "@njit\n",
                "def compute_call_price(model, M=default_M):\n",
                "    \"\"\"\n",
                "    Generate M observations of S_n and average to estimate\n",
                "    the option price.\n",
                "    \"\"\"\n",
                "\n",
                "    # Set up\n",
                "    β, μ, K, n, ρ, ν, S0, h0 = model\n",
                "    current_sum = 0.0\n",
                "    \n",
                "    # For each sample path\n",
                "    for m in range(M):\n",
                "        s = np.log(S0)\n",
                "        h = h0\n",
                "        \n",
                "        # Simulate forward in time\n",
                "        for t in range(n):\n",
                "            s = s + μ + np.exp(h) * randn()\n",
                "            h = ρ * h + ν * randn()\n",
                "            \n",
                "        # And add the value max{S_n - K, 0} to current_sum\n",
                "        current_sum += np.maximum(np.exp(s) - K, 0)\n",
                "        \n",
                "    return β**n * current_sum / M"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "6382e6ed",
            "metadata": {},
            "source": [
                "Let's time it."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 13,
            "id": "ab8ec64a",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "ab8ec64a",
                "lines_to_next_cell": 2,
                "outputId": "082d3534-d083-4dd5-93fd-92665f79823f"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 10.1 s, sys: 7.08 ms, total: 10.1 s\n",
                        "Wall time: 10.1 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "838.4871579285684"
                        ]
                    },
                    "execution_count": 13,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price(default_model)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "36634042",
            "metadata": {},
            "source": [
                "We'll call it again to exclude compile time."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 14,
            "id": "20afabf5",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "ab8ec64a",
                "lines_to_next_cell": 2,
                "outputId": "082d3534-d083-4dd5-93fd-92665f79823f"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 10.3 s, sys: 52.6 ms, total: 10.4 s\n",
                        "Wall time: 10.2 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "869.2102128128225"
                        ]
                    },
                    "execution_count": 14,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price(default_model)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "e99ee41f",
            "metadata": {
                "id": "e99ee41f"
            },
            "source": [
                "### Acceleration via `prange`\n",
                "\n",
                "Let's use `prange` from Numba to parallelize this codeand make it faster.\n",
                "\n",
                "Note that, at this stage, parallelization targets host chips (onboard CPUs), not the GPU."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 15,
            "id": "ea4d17c9",
            "metadata": {
                "id": "ea4d17c9",
                "lines_to_next_cell": 1
            },
            "outputs": [],
            "source": [
                "from numba import prange"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 16,
            "id": "6980ff3a",
            "metadata": {
                "id": "6980ff3a"
            },
            "outputs": [],
            "source": [
                "@njit(parallel=True)\n",
                "def compute_call_price_parallel(model, M=default_M):\n",
                "\n",
                "    β, μ, K, n, ρ, ν, S0, h0 = model\n",
                "\n",
                "    current_sum = 0.0\n",
                "    # For each sample path\n",
                "    for m in prange(M):      # Note use of prange\n",
                "        s = np.log(S0)\n",
                "        h = h0\n",
                "        # Simulate forward in time\n",
                "        for t in range(n):\n",
                "            s = s + μ + np.exp(h) * randn()\n",
                "            h = ρ * h + ν * randn()\n",
                "        # And add the value max{S_n - K, 0} to current_sum\n",
                "        current_sum += np.maximum(np.exp(s) - K, 0)\n",
                "        \n",
                "    return β**n * current_sum / M"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 17,
            "id": "0a98e90b",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "0a98e90b",
                "outputId": "da705bad-0917-4f05-d6e5-ec0f6104f26c"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 15.2 s, sys: 36 ms, total: 15.3 s\n",
                        "Wall time: 2.5 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "856.8948485567332"
                        ]
                    },
                    "execution_count": 17,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time\n",
                "compute_call_price_parallel(default_model)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 18,
            "id": "d8140f11",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "d8140f11",
                "lines_to_next_cell": 2,
                "outputId": "b957919d-2724-4e76-b2ef-af9e5cc4c30a"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 14.8 s, sys: 19.4 ms, total: 14.8 s\n",
                        "Wall time: 2.03 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "873.4848984973764"
                        ]
                    },
                    "execution_count": 18,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time\n",
                "compute_call_price_parallel(default_model)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "d26d83a8",
            "metadata": {
                "id": "d26d83a8"
            },
            "source": [
                "### A GPU Implementation Using CUDA"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "353726be",
            "metadata": {},
            "source": [
                "Next we will use a GPU implmentation using CUDA, which is an interface to GPU programming designed by Nvidia.  Rather than using CUDA directly, which requires C, we use a Numba-based front end --- which is just as fast, if not faster."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "9eaf0df5",
            "metadata": {},
            "source": [
                "If you are using Colab, you need to first need to make Numba aware of the location of certain libraries.  Uncomment the next two cells."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 19,
            "id": "c45c423c",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "1_lbXMhOTRra",
                "outputId": "4e455ec2-deea-4ca5-e24a-cb3dfcd6ded8"
            },
            "outputs": [],
            "source": [
                "#!find / -iname 'libdevice'\n",
                "#!find / -iname 'libnvvm.so'"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "cecf954d",
            "metadata": {},
            "source": [
                "Change `cuda-xxx` to the version displayed."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 20,
            "id": "14b93d5c",
            "metadata": {
                "id": "26bAoSYNTXFf"
            },
            "outputs": [],
            "source": [
                "#import os\n",
                "#os.environ['NUMBAPRO_LIBDEVICE'] = \"/usr/local/cuda-11.2/nvvm/libdevice\"\n",
                "#os.environ['NUMBAPRO_NVVM'] = \"/usr/local/cuda-11.2/nvvm/lib64/libnvvm.so\""
            ]
        },
        {
            "cell_type": "markdown",
            "id": "fcee7eb6",
            "metadata": {},
            "source": [
                "We use the following libraries."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 21,
            "id": "dc42806e",
            "metadata": {
                "id": "dc42806e"
            },
            "outputs": [],
            "source": [
                "from math import log, exp\n",
                "from numba import cuda, objmode, float32, int32\n",
                "from numba.cuda.random import xoroshiro128p_normal_float32 as xnf32\n",
                "from numba.cuda.random import create_xoroshiro128p_states as x128_states"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "a2bceb31",
            "metadata": {},
            "source": [
                "GPUs typically operate on 32 bit floats so we use the following function to convert as required."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 22,
            "id": "2a4dd2f3",
            "metadata": {},
            "outputs": [],
            "source": [
                "def convert_to_32_bit(scalar):\n",
                "    \"Convert scalar to 32 bits.  Must be a float or integer.\"\n",
                "    if isinstance(scalar, float):\n",
                "        return float32(scalar)\n",
                "    return(int32(scalar))    "
            ]
        },
        {
            "cell_type": "markdown",
            "id": "421340fa",
            "metadata": {},
            "source": [
                "Here's a function to update the log price of the asset and the volatility term by one period."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 23,
            "id": "2ef635b3",
            "metadata": {
                "id": "2ef635b3"
            },
            "outputs": [],
            "source": [
                "@cuda.jit(device=True)\n",
                "def update_price(params, s, h, ξ_1, ξ_2):\n",
                "    \"Update the state given standard normal draws ξ_1 and ξ_2.\"\n",
                "    β, μ, K, n, ρ, ν = params\n",
                "    s = s + μ + exp(h) * ξ_1\n",
                "    h = ρ * h + ν * ξ_2\n",
                "    return s, h"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "c0bc33b5",
            "metadata": {},
            "source": [
                "Now we write a **kernel** function, which describes the task of a single thread on the GPU."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 24,
            "id": "797f1f05",
            "metadata": {
                "id": "797f1f05"
            },
            "outputs": [],
            "source": [
                "@cuda.jit\n",
                "def price_generation_kernel(params, init_vals, price_sample, rng_states):\n",
                "    \"\"\"\n",
                "    Kernel function for the GPU.  Produces one Monte Carlo observation of S_n.\n",
                "    \n",
                "    This draw of S_n will be written into a slot in the array `price_sample` \n",
                "    specific to each thread.\n",
                "    \"\"\"\n",
                "\n",
                "    # Get position and unpack parameter values\n",
                "    pos = cuda.grid(1)\n",
                "    β, μ, K, n, ρ, ν = params\n",
                "    S0, h0 = init_vals\n",
                "    \n",
                "    s, h = log(S0), h0\n",
                "    for t in range(n):\n",
                "\n",
                "        # Draw shocks and then update price\n",
                "        ξ_1 = xnf32(rng_states, pos)\n",
                "        ξ_2 = xnf32(rng_states, pos)\n",
                "        s, h = update_price(params, s, h, ξ_1, ξ_2)\n",
                "\n",
                "    price_sample[pos] = exp(s)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "a1e98800",
            "metadata": {},
            "source": [
                "Next we create a host function that allocates memory on the device, launches the kernel, collects the output and returns an estimate of the option price."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 25,
            "id": "c0804638",
            "metadata": {
                "id": "c0804638"
            },
            "outputs": [],
            "source": [
                "def compute_call_price_cuda(model,\n",
                "                            M=default_M,\n",
                "                            seed=1234, \n",
                "                            threads_per_block=132):\n",
                "    \"\"\"\n",
                "    This function launches the kernel and returns a \n",
                "    Monte Carlo estimate of the price\n",
                "    \"\"\"\n",
                "    # Set up\n",
                "    nthreads = threads_per_block\n",
                "    nblocks = M // nthreads\n",
                "    β, μ, K, n, ρ, ν, S0, h0 = model\n",
                "    params = tuple(map(convert_to_32_bit, (β, μ, K, n, ρ, ν)))\n",
                "    init_vals = tuple(map(convert_to_32_bit, (S0, h0)))\n",
                "\n",
                "    # Generate the random states\n",
                "    rng_states = x128_states(M, seed=seed)\n",
                "\n",
                "    # Allocate memory\n",
                "    price_sample = cuda.to_device(np.zeros(M))\n",
                "\n",
                "    # Launch kernel and compute prices\n",
                "    price_generation_kernel[nblocks, nthreads](\n",
                "        params, init_vals, price_sample, rng_states\n",
                "    )\n",
                "\n",
                "    # Compute average payoff and discount\n",
                "    price_sample = price_sample.copy_to_host()\n",
                "    payoffs = np.maximum(price_sample - K, 0)\n",
                "    return β**n * np.mean(payoffs)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 26,
            "id": "474d80fd",
            "metadata": {
                "colab": {
                    "base_uri": "https://localhost:8080/"
                },
                "id": "474d80fd",
                "outputId": "66221def-87da-4783-b989-54b2aca8308a"
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 5.02 s, sys: 180 ms, total: 5.2 s\n",
                        "Wall time: 5.24 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "849.0907210556882"
                        ]
                    },
                    "execution_count": 26,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time\n",
                "compute_call_price_cuda(default_model)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "11b27cba",
            "metadata": {},
            "source": [
                "This result is disappointing.  The main issue is that the kernel performs a very small task, so the overhead of parallelization exceeds the speed gains.\n",
                "\n",
                "We could try to reoptimize this --- or we could just use JAX and let Google's engineers do the job for us."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "33bd155d",
            "metadata": {
                "id": "a74eshKnWVyP"
            },
            "source": [
                "## JAX Version"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 27,
            "id": "3e74526f",
            "metadata": {},
            "outputs": [],
            "source": [
                "import jax\n",
                "import jax.numpy as jnp"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0fdba2f5",
            "metadata": {},
            "source": [
                "Our JAX-based function keeps track of all $M$ prices at once and shifts them forward in time over $n$ steps.  Then it computes the sample mean and the option price."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 28,
            "id": "d1c39b5a",
            "metadata": {},
            "outputs": [],
            "source": [
                "\n",
                "def compute_call_price_jax(model=default_model,\n",
                "                           M=default_M,\n",
                "                           key=jax.random.PRNGKey(1)):\n",
                "    \"\"\"\n",
                "    Estimate the price of the call option using Monte Carlo.\n",
                "    \"\"\"\n",
                "    \n",
                "    β, μ, K, n, ρ, ν, S0, h0 = model\n",
                "\n",
                "    s = jnp.full(M, jnp.log(S0))\n",
                "    h = jnp.full(M, h0)\n",
                "    for t in range(n):\n",
                "        key, subkey = jax.random.split(key)\n",
                "        Z = jax.random.normal(subkey, (2, M))\n",
                "        s = s + μ + jnp.exp(h) * Z[0, :]\n",
                "        h = ρ * h + ν * Z[1, :]\n",
                "    expectation = jnp.mean(jnp.maximum(jnp.exp(s) - K, 0))\n",
                "        \n",
                "    return β**n * expectation"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 29,
            "id": "e642fe26",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 696 ms, sys: 36.4 ms, total: 733 ms\n",
                        "Wall time: 957 ms\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "DeviceArray(883.57184, dtype=float32)"
                        ]
                    },
                    "execution_count": 29,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price_jax().block_until_ready()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 30,
            "id": "406cecf2",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 23.1 ms, sys: 8.01 ms, total: 31.1 ms\n",
                        "Wall time: 44.7 ms\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "DeviceArray(883.57184, dtype=float32)"
                        ]
                    },
                    "execution_count": 30,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price_jax().block_until_ready()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "86c6e930",
            "metadata": {},
            "source": [
                "We can do even better by using JAX to JIT-compile the function."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 31,
            "id": "32820138",
            "metadata": {},
            "outputs": [],
            "source": [
                "compute_call_price_jax_jit = jax.jit(compute_call_price_jax)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "3514085b",
            "metadata": {},
            "source": [
                "Let's time it."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 32,
            "id": "8bd16be7",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 1.23 s, sys: 55.9 ms, total: 1.29 s\n",
                        "Wall time: 1.64 s\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "DeviceArray(883.57184, dtype=float32)"
                        ]
                    },
                    "execution_count": 32,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price_jax_jit().block_until_ready()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0bf64a88",
            "metadata": {},
            "source": [
                "Now here's the timing after compilation is complete:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 33,
            "id": "b926f388",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU times: user 3.32 ms, sys: 21 µs, total: 3.34 ms\n",
                        "Wall time: 21.4 ms\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "DeviceArray(883.57184, dtype=float32)"
                        ]
                    },
                    "execution_count": 33,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "%%time \n",
                "compute_call_price_jax_jit().block_until_ready()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "c85d8e7f",
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "accelerator": "GPU",
        "colab": {
            "provenance": []
        },
        "gpuClass": "premium",
        "kernelspec": {
            "display_name": "Python 3 (ipykernel)",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.9.12"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 5
}