{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DiscreteDP Example: Option Pricing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Daisuke Oyama**\n",
    "\n",
    "*Faculty of Economics, University of Tokyo*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Miranda and Fackler, <i>Applied Computational Economics and Finance</i>, 2002,\n",
    "Section 7.6.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "import quantecon as qe\n",
    "from quantecon.markov import DiscreteDP, backward_induction, sa_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = 0.5       # Time expiration (years)\n",
    "vol = 0.2     # Annual volatility\n",
    "r = 0.05      # Annual interest rate\n",
    "strike = 2.1  # Strike price\n",
    "p0 = 2        # Current price\n",
    "N = 100       # Number of periods to expiration\n",
    "\n",
    "# Time length of a period\n",
    "tau = T/N\n",
    "\n",
    "# Discount factor\n",
    "beta = np.exp(-r*tau)\n",
    "\n",
    "# Up-jump factor\n",
    "u = np.exp(vol*np.sqrt(tau))\n",
    "\n",
    "# Up-jump probability\n",
    "q = 1/2 + np.sqrt(tau)*(r - (vol**2)/2)/(2*vol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We follow the state-action pairs formulation approach.\n",
    "We let the state space consist of the possible values of the asset price\n",
    "and the state indicating that \"the option has been exercised\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Possible price values\n",
    "ps = u**np.arange(-N, N+1) * p0\n",
    "\n",
    "# Number of states\n",
    "n = len(ps) + 1  # State n-1: \"the option has been exercised\"\n",
    "\n",
    "# Number of actions\n",
    "m = 2  # 0: hold, 1: exercise\n",
    "\n",
    "# Number of feasible state-action pairs\n",
    "L = n*m - 1  # At state n-1, there is only one action \"do nothing\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Arrays of state and action indices\n",
    "s_indices, a_indices = sa_indices(n, m)\n",
    "s_indices, a_indices = s_indices[:-1], a_indices[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reward vector\n",
    "R = np.empty((n, m))\n",
    "R[:, 0] = 0\n",
    "R[:-1, 1] = strike - ps\n",
    "R = R.ravel()[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Transition probability array\n",
    "Q = sparse.lil_matrix((L, n))\n",
    "for i in range(L-1):\n",
    "    if a_indices[i] == 0:\n",
    "        Q[i, min(s_indices[i]+1, len(ps)-1)] = q\n",
    "        Q[i, max(s_indices[i]-1, 0)] = 1 - q\n",
    "    else:\n",
    "        Q[i, n-1] = 1\n",
    "Q[L-1, n-1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a DiscreteDP\n",
    "ddp = DiscreteDP(R, Q, beta, s_indices, a_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The backward induction algorithm for finite horizon dynamic programs is offered\n",
    "as the function `backward_induction`.\n",
    "(By default, the terminal value function is set to the vector of zeros.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vs, sigmas = backward_induction(ddp, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the returns, `vs` is an $(N+1) \\times n$ array that contains the optimal value functions,\n",
    "where `vs[0]` is the value vector at the current period (i.e., period $0$)\n",
    "for different prices\n",
    "(with the value `vs[0, -1] = 0` for the state \"the option has been exercised\" included),\n",
    "and `sigmas` is an $N \\times n$ array that contins the optimal policy functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v = vs[0]\n",
    "max_exercise_price = ps[sigmas[::-1].sum(-1)-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtgAAAEWCAYAAABc2eBcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeYVOX5xvHvsw0W2KXsLr0jFhRBg4rYe4lRgybBiNgi\nEk0iJtbEjiXRqLGiqIhYfzH2iiUqiqACAQUrVap0WPouPL8/zlkdcWFn2Zk9M7P357rOxZQzZ+4z\nO5x55p33vK+5OyIiIiIikhhZUQcQEREREckkKrBFRERERBJIBbaIiIiISAKpwBYRERERSSAV2CIi\nIiIiCaQCW0REREQkgVRgS51mZq+Z2elR59iSmc0ys8OjziEiqSNVj1dbY2ZTzezgJGz3YDObm+jt\nJpKZ/dXMHow6R7KY2TVm9ljUOVKZCmxJmLAoXGdmq83sOzMbYWaN4nhcRzNzM8upYr1uZvaima00\ns1Ize8fM+lQj308OCO5+jLs/Eu824nyey8xsdCW3F5vZRjPbLZHPJyLVp+PVT56vj5n9N8y60sxe\nMrNu1Xj8CDO7PvY2d9/V3d9NeNiqs7iZrQn/thXLJbWZwd1vdPffJXKb4XuiLGafvjCzkxL5HJI4\nKrAl0X7h7o2APYFewBWJ2KiZdQHGAJ8BnYDWwHPAG2a2byKeI4EeA/qYWactbu8HfObuUyLIJCI/\npeMVEGZ6A3iBIGsnYDIwxsw6R5mtBnq4e6OY5eZEbryqL1hJ9H8V+wQMBh4zsxYRZdkuZpYddYba\noAJbksLd5wGvAbvBT7s8bNE6U9HauyL8Vl7ZB9A1wFh3/5u7L3P3Une/E3gU+Ee4zYqWpYFmNt/M\nFpjZReF9RwN/BX4TPsfk8PZ3zex34eUsM7vCzGab2SIzG2lmjbfY9ulm9q2ZLTGzv21l3+cC/wVO\n2+KuAcDIcHtdwtaipeG2HjezJpVtb8uWoS1/HjWz1mb2jJktNrOZZvanyrYjIpWry8er0M3ASHe/\nI8y6zN2vAMaF+/L9cceCrg9Lwtfo1PC+gcCpwCVh3pe2fB3D1/BpM3ssbCX/zMx2NLPLw/xzzOzI\nmNf8zLCFttTMZpjZuVX9HeNhZq+a2a0x158ys+Ex188Kn3e5mY0ysw4x97mZnW9m3wDfhLftamZv\nmtkyC34J+WvM/j4WXq4f7vdSM1thZp9UFMVm1tjMHgr//vPM7HqLswB191FAKdAlJuM5ZjYtzPOi\nmbUOb//JLy9bvJ/OMLMPzOyf4b7PNLNjYtbtZGbvhX+PN4HiLV7Xp81soQW/fow2s11j7hthZkPD\n134N8OfwtcqOWadvxfs8U6jAlqQws3bAscD/4lj9wPDfJuE387GVrHME8HQlt/8b2M/M8mNuOwTo\nChwJXGpmh7v768CN/PDtv0cl2zojXA4BOgONgLu3WGd/YCfgMOAqM9tlK/v0CDEFtpntBPQEnqi4\nCbiJoLVoF6Ad4QdZdZhZFvASQWtTmzDXYDM7qrrbEqmr6vLxyswaAH22kfeImOstCQqrNsDpwDAz\n28ndhwGPAzeHeX9RybYAfkHwJaMpwWs9iqAOaQNcB9wfs+4i4DigEDgTuN3M9tzKdqvjLOA0Mzs0\n/IKwN3ABgJmdQPDFpi9QArwPPLnF408E9gG6mVkB8BbwOsGxfAfg7Uqe83SgMcFxvggYBKwL7xsB\nlIeP3YPgfVBl1xIL/BzIAz4PbzuU4HPl10ArYDbwVFXbirEP8BXB3/hm4CEzs/C+J4AJ4X1Dwn2K\n9RrB+7g5MJHg/RDrt8ANQAFwF7CUYF8rnEbYAJUpVGBLoj1vZiuAD4D3CD4kEqEYWFDJ7QsI3sfN\nYm671t3XuPtnwMPAKXE+x6nAbe4+w91XA5cD/ezHPwVe6+7r3H0yQVFb2QcfBD8Ht7Af+lwOAF5z\n98UA7j7N3d909w3hbbcBB8WZM9ZeQIm7X+fuG919BvAAQXcUEdk2Ha+CLFnbyFu8xW1Xhset94BX\nCIq5eL3v7qPcvZygoC8B/u7uZQSFYEcLf8lz91fcfboH3iPownJANZ5rYthaXLEcFW53IfB7gkaQ\nO4AB7l4aPmYQcJO7fxFmvBHoGduKHd6/zN3XEXwBWOjut7r7+rD1/6NKspQRFNY7uPsmd5/g7qvC\nVuxjgcHhe2ARcDvbPn7/OnzPrgZeBG509xXhfacCw919ortvIHhP7GtmHeN8zWa7+wPuvil8fVoR\nfI61J/isqfjbjyZo2Pmeuw8P938DQWNRDwt/UQm94O5j3H2zu68Pt98fwMyaAUfxQwNURlCBLYl2\nors3cfcO7n5eeBBKhCUE/9m31ArYDCyPuW1OzOXZBC0L8Wgdrh/72Bwgtn/bwpjLawlajX7C3dcS\nfIAMCFsATiXm27mZtQh/mpxnZqsI+m1v+UEWjw5A69gPEoIWmLTqkycSER2vgiybt5F3Sey67r5m\nO/MCfBdzeR2wJCzmKq5TkdHMjjGzcWFXhxUEhWh1jpF7hn/bimVUzH0vAdnAV+7+QcztHYA7Yo6l\nywh+bWwTs07s36sdMD2OLI8StNY/ZUF3oJvNLDd8vlxgQcxz3k/QCrw1/w73pyFB15AB9kP3mR+9\nJ8IvXku3yL8t379fws8wCP4eran8bw8EfarN7O9mNj38PJsV3hX794p93SD4zPuFmTUk+JL2vrtX\n9iUvbanAltqyBmgQc71lzGWP4/FvAb+q5PZfE/R1XBtzW7uYy+2B+XE+z3yCA17sY8v58YdCdTwS\n5juC4Gex2G/8N4Z5urt7IcE3efvJFgLbeu3mADO3+CApcPdjtzOziNSh41VYNI3dRt7YLg9Nw4Jo\ne/LGzczqAc8A/wRauHsT4FW2foysrhuAL4BWZhb7i8Ec4Nwtjqf57v5hzDq+xfpVngTq7mXufq27\ndyPojnMcwa+ac4ANQHHM8xW6+67b2l7MdmcRdM2o6JLzo/dE+LcqAuYRvKdh6+/rbVlA5X/7Cr8F\nTgAOJ+gK07EiQmzcLbLPI3jf9SXoHvJonFnShgpsqS2TCH6+zDWzXsDJMfctJmhB2daB6lqCkTlu\nMLNmZlZgZn8kOEhdusW6V5pZg/AkizOB/wtv/47gJ8itve+fBC4MT+ZoxA99IMurs6Mx3gdWAMOA\np9x9Y8x9BQQ/8a00szbAxdvYziTg2HC/WxKcOV7hY6DUzC41s/ywJWE3M9trOzOLSN07Xl0GnG5m\nfwqzNrXgxOp9w3350b6ZWZ6ZHUBQKFb03f6OOIrNOOUB9Qhe63ILTrY7ctsPiY+ZHUjwOg8g6Ed8\nV3gMBrgPuLziBD0LTkCs7ItHhZcJivTBZlYvfO32qeQ5DzGz7hac1LeKoMvI5rDF9g3gVjMrtODE\n1S5mFld3QTNrCxwNTA1vehI408x6hl9SbgQ+cvdZYVfEeUD/8HPiLGJOjtwWd58NjOeHv/3+/FDU\nQ/B5toGgtbwB8Xe1GglcAnQHno3zMWlDBbbUlisJ/jMvJzhgf9/XKmzNuYFgSKgVZtZ7ywe7+zcE\nJ+z0IPj5aQFwEnCUu4/ZYvX3gGkELS//dPc3wtsrPgiWmtnESjIOJ/gWPRqYCawH/ljtPf0hsxMc\nQDrw05M3riUYGmwlQT/GbR1cHiXoPzmL4GBc8QFM+PPqcQQnUM4k+Dn3QYJWBBHZPnXqeBV2kziK\noDVxAcHP/3sA+4f7UmEhwWsyn+AktkHu/mV430MEJ/6tMLPntydHTJ5S4E8EJ1kuJ2ghfbGam5ls\nPx4H+19mVkhwLP6Du89z9/fD3A+bmbn7cwSjvDwVdnWYAhyztScIcx5BUGwuJBhZ5JBKVm0J/Ieg\nuP6C4G9e0WI7gB9OVFwerldZd50KFSPLrAY+IRgO8towz1sE791nCP6OXfhxf+5zCBpzlgK7ArEt\n81X5LcFJkMuAq/nxZ9pIgvfMvHA/xsW5zecIPh+f2+JXnYxgQQ0gkv7CEzlmArk1aHUWEUm6dDte\nWTAj42Pu3jbqLJI5zGw6Qbect6LOkmhqwRYRERGRWmXBLJROMG9ExlGBLSIimFk7C6bz/tzMpprZ\nBZWss7OZjTWzDRZOiiIiUl1m9i4wFDjf3TdHHCcp1EVEREQws1ZAK3efaMEEGhMIhrH7PGad5gR9\nJk8kGLbrn9GkFRFJbWrBFhER3H2Bu08ML5cSnIzVZot1Frn7JwSjIIiIyFbkVL1KasnJyfEddtiB\nRo0qnd9DRCSlTZgwYYm7l0SdY1vCE/D2ACqblS6exw8EBgI0bNjwZzvvvHPCsomI1KYJEyasJphC\nvsIwdx9W1ePSrsDOzs6mrKyMd955h4KCgqjjiIhUi5nNrnqt6IRjKj9DMH3zqu3ZRvjhMwygV69e\nPn78+AQmFBGpPWb2lbv3qu7j0q6LSKdOnZg1axaDBw+uemUREYlbOH3zM8Dj7p5xEz+IiNSWtCuw\nGzVqxGWXXcbw4cN5/vkajWUvIiIhMzOCSTe+cPfbos4jIpLO0q7ABrj66qvZc889GTp0aNRRREQy\nxX7AacChZjYpXI41s0FmNgjAzFqa2Vzgz8AVZjY3nB1PRERipF0fbIC8vDxeeukliouLo44iIpIR\nwumyrYp1FgLVmslvzvK1rFxXRuP83JrEExFJK2nZgg3QunVr8vLyWL58OW+++WbUcUREpBIr1pYx\ncfbyqGOIiNSqtC2wK1x44YWccMIJfPXVV1WvLCIitcqAj2YuizqGiEitSlqBHee0u2Zmd5rZNDP7\n1Mz2rO7z3HjjjeTn59O/f3/KyjT3gYhIKsnPy+aTWSqwRaRuSWYLdjnwF3fvBvQGzjezbluscwzQ\nNVwGEsxLXy2tW7fmgQceYPz48Vx33XU1zSwiIgnUMC+HT+euYN3GTVFHERGpNUkrsOOZdhc4ARjp\ngXFAEzNrVd3n6tu3L2eccQY33ngjH374YY2zi4hIYjSsl0PZJud/c9QPW0TqjloZRWQb0+62AebE\nXJ8b3rZgi8d/P+1u+/btK32OO+64AzOjY8eOiYgsIiIJ0KBeNhsNPp65jD5dNPKTiNQNST/JMVHT\n7rp7L3fvVVJSUuk6hYWFDB8+nNatW+PuNUgsIiKJkm3GLi0L1Q9bROqUpBbYcUy7Ow9oF3O9bXjb\ndluyZAmHH364ZnkUEUkRe3dqxoTZy9lYvjnqKCIitSKZo4jEM+3ui8CAcDSR3sBKd1+wlXXjUlhY\nyIoVKzjnnHNYuHBhTTYlIiIJsE+nZqwv28yU+SujjiIiUiuS2YJd5bS7wKvADGAa8ABwXk2fNC8v\nj8cee4zVq1dz1llnqbuIiEjE9urUDAj6YYuI1AXJHEXkA3c3d9/d3XuGy6vufp+73xeu4+5+vrt3\ncffu7j4+Ec+9yy67cMstt/Daa69x3333JWKTIiKynYob1aNzSUMV2CJSZ6T9TI5bc/7553PUUUdx\n1113UV5eHnUcEZE6bZ9Ozfhk1jI2bdaviiKS+TK2wDYzRo4cydixY8nJqZXRCEVEZCv27tSM0vXl\nfLlwuwaTEhFJKxlbYAM0b96cxo0bs2HDBkaNGhV1HBGROmvvTkWA+mGLSN2Q0QV2hRtuuIFjjz1W\nszyKiESkTZN82jTJV4EtInVCnSiwL7roItq3b89pp51GaWlp1HFEROqkfTo34+OZyzS6k4hkvDpR\nYBcWFvLoo48ya9YsLrzwwqjjiIjUSb07FbF0zUamLVoddRQRkaSqEwU2wP77789ll13GQw89pFke\nRUQisHc4HvZH6iYiIhmuTg2vcfXVVzNv3jy6du0adRQRkTqnQ1EDWhTW46OZy+jfu0PUcUREkqZO\nFdh5eXmMGDHi++vuTjCju4iIJJuZsU+nIsbNWKrjr4hktDrTRSTWunXrOO2007j//vujjiIiUqfs\n3akZi0o3MHvp2qijiIgkTZ0ssOvVq8eiRYv485//zFdffRV1HBGROqN354p+2EsjTiIikjx1ssDO\nysri4YcfJj8/n/79+1NWVhZ1JBGROqFLSSOKGubx0Qyd6CgimatOFtgArVu3ZtiwYYwfP54hQ4ZE\nHUdEpE4wM/bu1EwjiYhIRquzBTbASSedxOmnn87QoUNZuXJl1HFEROqEfTo1Y96Kdcxdrn7YIpKZ\n6nSBDXDnnXfyv//9j8aNG0cdRUSkTti7UxGAuomISMaq8wV2YWEhbdu2ZfPmzbz55ptRxxERyXg7\ntyygcX6uTnQUkYxV5wvsCsOHD+fII4/ULI8iIkmWlWXs06kZ49SCLSIZSgV2aMCAAey5556cc845\nLFy4MOo4IiIZbZ/ORXy7bC3zV6yLOoqISMKpwA7l5eXx2GOPsXr1as4++2zcPepIIiIZS+Nhi0gm\nU4EdY5ddduHmm2/m1Vdf5b777os6johIxtq5ZSGF9XN0oqOIZKScqAOkmvPPP59Jkyax8847Rx1F\nRCRjZWcZe3cqYtwMtWCLSOZRgb2FrKwsHnrooahjiIhkvN6dm/HWF9+xcOV6WjauH3UcEZGEUReR\nrdi8eTOXX34511xzTdRRREQyUu/O4XjY6octIhlGBfZWZGVlsWDBAoYMGcLYsWOjjiMiknF2aVVI\nQf0cDdcnIhlHBfY23HnnnbRv357+/ftTWloadRwRkaQxs3Zm9o6ZfW5mU83sgkrWMTO708ymmdmn\nZrZnTZ4zO8vYu2MzPlI/bBHJMCqwt6GwsJCRI0cyc+ZMLrzwwqjjiIgkUznwF3fvBvQGzjezblus\ncwzQNVwGAkNr+qS9OxcxY8kaFq1aX9NNiYikDBXYVTjggAO49NJLefTRR5k5c2bUcUREksLdF7j7\nxPByKfAF0GaL1U4ARnpgHNDEzFrV5Hkr+mGPVSu2iGQQFdhxuPbaa/nf//5Hp06doo4iIpJ0ZtYR\n2AP4aIu72gBzYq7P5adFOGY20MzGm9n4xYsXb/O5urUupKBeDh/NVD9sEUlJxRXHs3AZGM+DVGDH\nIS8vj27dgl9KR48erVkeRSRjmVkj4BlgsLuv2p5tuPswd+/l7r1KSkq2uW4wHnYzxk1XC7aIpKQl\nFcezcBkWz4NUYFfDm2++yUEHHcT9998fdRQRkYQzs1yC4vpxd3+2klXmAe1irrcNb6uRin7Y36kf\ntohkCBXY1XDYYYdx5JFH8pe//IWvv/466jgiIgljZgY8BHzh7rdtZbUXgQHhaCK9gZXuvqCmz71v\nl6AftmZ1FJFMoQK7GrKysnj44YepX78+/fv3p6ysLOpIIiKJsh9wGnComU0Kl2PNbJCZDQrXeRWY\nAUwDHgDOS8QT79KqkML6OSqwRSRjaKr0amrdujX3338/v/rVr7j0un9w25Aroo4kIlJj7v4BYFWs\n48D5iX7uoB92EWPVD1tEMoRasLfDySefTN/BN/B8WU8mzF4edRwRkbTXu3MzZi1dy4KV66KOIiJS\nYyqwt9PDN11MqyYNuPD/JrF6Q3nUcURE0lrFeNjqJiIimUAF9nYqrJ/L7b/pybfL1nDs3x6KOo6I\nSFrr1qqQxvm5jJuu8bBFJP2pwK6BvTs1Yxfm8m1OW4Y8/GLUcURE0lZWOB62ZnQUkUygAruGnr7q\nDLJXzefBSWuYOn1O1Q8QEZFK7du5iG+XrWXeCvXDFpH0plFEaqhRg/rc038vzn12Br+57SU+vWsQ\nWVn63iIiUl0V/bD/8MREihvV+8n9p+zdjkN3blHbsUREqk2VYAIc3acnhzQrZXVBB4Y8+V7UcURE\n0tLOLQs4eteWrC/bzNzl6360jJu+lHvemR51RBGRuKgFO0EevKQ/J935X578cj39F6+mS0mjqCOJ\niKSVrCzjvtN+Vul9/3j9Sx4YPYM1G8ppWE8fXSKS2tSCnSA5OdkM+90B1M/NZuDwMaxdvyHqSCIi\nGWPfzkWUb3bGa+4BEUkDSSuwzWy4mS0ysylbuf9gM1sZMyXvVcnKUltaFNbnvL2aMn15Of2GPBJ1\nHBGRjNGrY1Nys40Ppy+JOoqISJWS2YI9Aji6inXed/ee4XJdErPUmoHH7EWLtbOYXNaSR14ZHXUc\nEZGM0CAvh57tmjBO06mLSBpIWoHt7qOBOjljwLNX9MPWLufq12by3dIVUccREckI+3Yp5rN5K1m1\nvizqKCIi2xR1H+w+Zvapmb1mZrtubSUzG2hm481s/OLFi2sz33Zp07yIK45ojzdsxi+HPB51HBGR\njLBv5yI2O3w8o0623YhIGomywJ4ItHf33YG7gOe3tqK7D3P3Xu7eq6SkpNYC1sTvTjiE7tkLmV+/\nI6OmLIg6johI2tujfRPycrI026OIpLzICmx3X+Xuq8PLrwK5ZlYcVZ5kePa6s9m1dSGXPzeFRaXr\no44jIpLW6udm06tDU8aqH7aIpLjICmwza2lmFl7eO8ySUUfNvJws/vWbnpSu28DRVz3G5s2bo44k\nIpLW9u1cxOcLVrF8zcaoo4iIbFUyh+l7EhgL7GRmc83sbDMbZGaDwlVOBqaY2WTgTqCfu3uy8kSl\na4sC9i9YyrJ6rRj0T/XHFhGpiX27BNOpfzQzo9pjRCTDJG06LHc/pYr77wbuTtbzp5IHLu7P7n+6\nn1FlLXn7kykcttduUUcSEUlLu7dtQoO8bP711je8NmVh3I/rUNSQPx+xYxKTiYj8QPPN1oKcnGye\nHPxzjr93HOc+/CGfdu9Kg/r1oo4lIpJ28nKyOK13B0ZNXcjkOfENg7pm4yZemDSfX/2sLe2aNUhy\nQhERFdi1pseOHRjQbQKPzmjDH4a+wvAL+0YdSUQkLV1+7C5cfuwuca8/bVEph982mjHTltBv7/ZJ\nTCYiEoh6HOw6ZcjAvhzQNpd3F9VjwmyN4yoiUhu6lDSieUE9PtToIyJSS1Rg17J7f3cIrZvk84fH\nxmuWRxGRWmBm9OlSxIfTl5KB59KLSApSgV3LCurncvVRHZm/aoNmeRQRqSV9uhSzZPUGvlm0Ouoo\nIlIHqMCOwBE9O9PN5jG/fkeuf/jFqOOIiGS8iuH9Ppy2JOIkIlIXqMCOyNNXnUH2qvk8MGkNU6fP\niTqOiEhGa9esAe2a5asftojUChXYEWnUoD539+8FufX4zW0vaZZHEZEk69O5mHEzlrJps/phi0hy\nqcCO0DF99uDQolJWF3Rg5NhZUccREclofXYoYtX6cqbOXxl1FBHJcCqwI/bQpQM4oGsxf3/9K6Yv\n1sk3IiLJ8n0/bHUTEZEkU4Edsaws45+/6kEOzvE3Psva9RuijiQikpGaF9Sna/NGKrBFJOk0k2MK\naFFYn77t1zNyehH9hjzCizcMjDqSiEhG6tOliCc/nkPfe8fE/Zjc7CyuO2E3dmpZkMRkIpJJ1IKd\nIq47py8t1s5kcllLHnnl/ajjiIhkpH57t6fPDkU0rJcT9zJh9nKenTg36ugikkbUgp1Cnr3iFPa7\n7mWufm05R/fuTouiJlFHEhHJKLu0KmTEmXtX6zG/vn8sY6Zr/GwRiZ9asFNIm+ZFXHFEe7xhM06/\n4+Wo44iICLBfl2Kmzl/FirUbo44iImlCBXaK+d0Jh3Bc5zy+3NiUUVMXRh1HRNKMBfqb2VXh9fZm\nVr0mW/mR/bsW4Q5jdXKkiMRJBXYKuv3sI9itTSGX/mcyn8+cF3UcEUkv9wL7AqeE10uBe+J5oJkN\nN7NFZjZlK/c3NbPnzOxTM/vYzHZLTOTUtnvbJjTMy+YDTbMuInFSgZ2C8nKyuPmXu7K8dC2//ucL\nmuVRRKpjH3c/H1gP4O7Lgbw4HzsCOHob9/8VmOTuuwMDgDtqkDNt5GZnsU/nIg3vJyJxU4Gdorq1\nbcahzVayuqADg259Iuo4IpI+yswsG3AAMysB4vqW7u6jgWXbWKUb8N9w3S+BjmbWomZx08N+OxQz\nc8ka5q9YF3UUEUkDKrBT2AOXnEaDlbMYtbABb39S6S+2IiJbuhN4DmhuZjcAHwA3Jmjbk4G+AGG/\n7g5A2y1XMrOBZjbezMYvXrw4QU8drf12CGaBHKNuIiJ1TXHF8Sxc4pqsRAV2CsvJyeaJwT/HN23k\n3Ic/1CyPIlIld38cuAS4CVgAnOjuTydo838HmpjZJOCPwP+ATZVkGObuvdy9V0lJSYKeOlo7tSig\nuFGeuomI1D1LKo5n4TIsngepwE5xPXfqyOnd8igvbMO/3vwq6jgikuLMrDcwz93vcfe7gXlmtk8i\ntu3uq9z9THfvSdAHuwSYkYhtpzozY98uxXwwbQnuHnUcEUlxcU00Y2bHAUMIfg7MAQxwdy9MYjYJ\nDRnYl7X/nsSDH87hqN3b8rMOzaKOJCKpayiwZ8z11ZXctl3MrAmw1t03Ar8DRrv7qppuN13sv0MR\nL02ez+S5K+lU1DCh287PyyYvR21eIpki3pkc/0XQ7+4z11f3SFxz/K6M+XoRv73zDUb/9WjN8igi\nW2Oxx2l332xm8TamPAkcTNDncC5wNZAbbuc+YBfgETNzYCpwdoKzp7Q+XYoBOPGeMQnfdpsm+bx3\n8cHkZKvIFskE8RbYc4ApKq6jU1A/l991z2XImAb0vf4Jxtx+XtSRRCQ1zTCzPxG0WgOcR5zdONz9\nlCruHwvsWLN46atdswbc139P5q9Yn9DtTlu8mic++pYp81fRs50aT0QyQbwF9iXAq2b2HvD9mXbu\nfltSUkmlfnfCITwzdihf1OvAkIdf5Mozj486koiknkEEI4lcQTBU39tAXGe9S9WO3q1Vwre5dPUG\nnvjoWz74ZrEKbJEMEe9vUTcAa4H6QEHMIrXs6avOJHvlPB6cvJap0+dEHUdEUoy7L3L3fu7e3N1b\nuPtv3X1R1Llk64oa1aNbq0LNFCmSQeJtwW7t7nViStxU16hBfe4+bS8GPTuTM4a+xce3nIGZRR1L\nRCJmZpe4+81mdhfhJDOx3P1PEcSSOO3ftZgRY2axdmM5DfLi/WgWkVQVbwv2q2Z2ZFKTSNyO6bMH\n/XdrxOKc5jw2bnbUcUQkNXwR/jsemFDJIils/x2K2bhpM5/MWh51FBFJgHi/Jv8euMjMNgBlaJi+\nyF1/+hHMefgTbnjlC7o0KqdP9x2ijiQiEXL3l8Ip0ru7+0VR55Hq2atjM/KysxgzbQkH7ZgZk/OI\n1GVxtWDO0NSTAAAgAElEQVS7e4G7Z7l7vrsXhtdVXEfIzLj5pO5sWLeaAUPf0SyPIoK7bwL2izqH\nVF9+XjY/69CUD75RP2yRTBBXgW1mB1a2JDucbFvLxvn03ymb8sLW9BsyMuo4IpIaJpnZi2Z2mpn1\nrViiDiVV279rMZ8vWMWS1WowEUl38fbBvjhmuRJ4CbgmSZmkGoYM7EuLtTOZXNaCR155P+o4IhK9\n+sBS4FDgF+FyXKSJJC777RBMZPPh9KURJxGRmoqrD7a7/yL2upm1I5jdUVLAs1ecwn7XvczVry3n\n6N7dNcujSN12sburn0Ea6t6mMYX1cxjzzRKO79E66jgiUgPbOxbQXIIpcyUFtGlexN8Ob8cNY9cw\n5JUvuHvAvlFHEpFaZma/AIYDZWa2Gfi1u38YcSyphuwso0+XYt7+8jtufv3LWnveQ3Zuzl4dm9Xa\n84nUBXEV2FuMq5oF9AQmJiuUVN85Jx7K8vpfcu+70zluyoKkzDYmIintBuAAd//SzPYBbgYOijiT\nVNPxPVvzzleLeOD9uGa3r7Hyzc5/v1zE64N1WpVIIsXbgj0+5nI58KS7j0lCHqmBwYfvyH8/X8Af\nR47l+XP3Ytcu7aKOJCK1p9zdvwRw94/MTLPtpqFju7fi2O6110ByzzvTuGXUVywu3UBJQb1ae16R\nTBfvMH2PxCyPq7hOTXk5WQzu3ZSNm+E3t73E5s2bo44kIrWnuZn9uWKp5LrITxzYNRhze4ymaRdJ\nqG0W2Gb27/Dfz8zs05jlMzP7tHYiSnUc3acnhzZbxeqCDgy69Ymo44hI7XkAKIhZtrwu8hO7ti6k\naYNc3tf42yIJVVUXkQvCfzXEUxp58NLT6P6H+xhV3pq3Pp7C4XvvFnUkEUkyd7826gySfrKyjD47\nFPPBtMW4O2YWdSSRjLDNFmx3XxD+O9vdZwPLgdKYRVJQdnY2Twz+OV6+gT8+/gkby9VVREREKnfA\nDsV8t2oD0xatjjqKSMaIdybHc81sIfApMCFcxm/7URKlnjt15II+JazLb86db38TdRwREUlR+3cN\nJrgZrW4iIgkT70yOFwG7uXtHd+8ULp2TGUxq7i+/PoyTf9aWe9+dxisffRF1HBGpBWbWKZ7bRCq0\nbdqATsUN+eCbxVFHEckY8RbY04G11dmwmQ03s0VmNmUr95uZ3Wlm08ITJ/eszvYlPlf/ohu5Zas5\n/7GP+W7piqjjiEjyPVPJbf+p9RSSVg7oWsxHM5epS6FIgsRbYF8OfGhm94dF8Z1mdmcVjxkBHL2N\n+48BuobLQGBonFmkGgrq53LRfiV4flP6Xv9k1HFEJEnMbGczOwlobGZ9Y5YzgPoRx5MUt/8Oxazd\nuImJ3y6POopIRoh3opn7gf8CnwFxfb1199Fm1nEbq5wAjHR3B8aZWRMza1VxYqUkzjm/PJRnxw3l\ni3rtGfLwi1x55vFRRxKRxNuJYMSnJsAvYm4vBc6JJJGkjd5disjOMl6fspAWhan7faxV4/rUz82O\nOoZIleItsHPdPdETFbQB5sRcnxve9pMC28wGErRy0759+wTHqBv+fdUZ9LzoUR6cXEjf6XM0y6NI\nhnH3F4AXzGxfdx8bdR5JL4X1c9mzfRNGfDiLER/OijrOVh3QtZhHz94n6hgiVYq3wH4tLHJfAjZU\n3Ojuy5KSagvuPgwYBtCrVy+vjefMNAUN87n7tL34/XOzueqVr/nPH9tqvFORzLTUzN4GWrj7bma2\nO3C8u18fdTBJbbf9uicTZqduF5E3P/+OUVMXsnpDOY3qxVu+iEQj3nfoKeG/l8fc5kBNRhKZB8Q2\no7YNb5MkOabPHlzjTbjmpc95dNxsBuzbMepIIpJ4DwAXE3Ttw90/NbMnABXYsk3tmjWgXbMGUcfY\nquaF9XjlswWMm76Uw7u1iDqOyDbFdZJjzNB8nRI4TN+LwIBwNJHewEr1v06+0/t0ZP8uzbjm+U95\n6+NKB3gRkfTWwN0/3uK28kiSiCTQzzo0JT83m/c1nKCkgXgnmmlgZleY2bDwelcz2+b06Wb2JDAW\n2MnM5prZ2WY2yMwGhau8CswAphG0uJy33XshcTMzLj6oFeUb1jJoxFjWrNtQ9YNEJJ0sMbMuBL8y\nYmYnU8m5LSLppl5ONr07N9OEOJIW4h2m72FgI9AnvD6PKn5udPdT3L2Vu+e6e1t3f8jd73P3+8L7\n3d3Pd/cu7t7d3TUzZC3psWNHTt8ll/LC1pxy/cio44hIYp1P0D1kZzObBwwGfh9tJJHEOKBrCTOX\nrGHOsmpNzSFS6+ItsLu4+81AGYC7rwV0hlwaG3LuSbRYM4PJZS145NUPoo4jIgni7jPc/XCgBNjZ\n3fd391kRxxJJiAN3LAHgfbViS4qLt8DeaGb5/PCTYxdiRhOR9PTslb+Ftcu59vWZrFq3Meo4IpIA\nZnaBmRUSzL57u5lNNLMjo84lkghdShrSunF99cOWlBdvgX018DrQzsweB94GLklaKqkVbZoXMeTn\nO0DDZgx5+Yuo44hIYpzl7quAI4Ei4DTg79FGEkkMM+OAriWMmbaE8k2a1l1SV5UFtgWDJX8J9AXO\nAJ4Eern7u0lNJrViwDH78fuDu/D0hLk8/p5GFRHJABXd944lmC13KurSJxnkgB2LWbW+nE/nrYw6\nishWVVlgh1OZv+ruS939FXd/2d3V+SmDXHDYjjRlNX99bgqfz5gbdRwRqZkJZvYGQYE9yswKADX1\nScbYr0sxZvD+1ypFJHXFO9HMRDPby90/SWoaiUReThY3/mJHBj07k9/c9hKT7zyXrKx4ew+JSIo5\nG+gJzHD3tWbWDDgz4kwiCdO0YR67t23C3e98w8Mfzow6To20a9qAZ8/rQ262PnMzTbwF9j5AfzOb\nBawh+LnR3X33ZAWT2nXMfntw6LuTeKe0PYNufZJhF58adSQR2T77ApPcfY2Z9Qf2BO6IOJNIQv31\nmJ159bP0Ht594ar1jJr6HZPmrGCvjs2ijiMJFm+BfVRSU0hKePCyAXT/w1BGlbfh7U+mctheu0Yd\nSUSqbyjQw8x6AH8BHgRGAgdFmkokgfbpXMQ+nYuijlEjK9eV8dYXbzL668UqsDPQNn+TMLP6ZjYY\nuBg4Gpjn7rMrllpJKLUmOzubJwb/HNtUxnVvzmZjubptiqSh8vDcmROAu939HqAgngea2XAzW2Rm\nlZ7xbGaNzewlM5tsZlPNTF1PRLZT4/xcerZrwntfa8jBTFRVp59HgF7AZ8AxwK1JTySR6rlTJ+49\ncz9mr3LuePvrqOOISPWVmtnlQH/gFTPLAnLjfOwIgsaUrTkf+NzdewAHA7eaWV4NsorUaQd2LeGz\neStZtkZzUWSaqgrsbu7e393vB04GDqiFTBKxY3dvw8l7tuHed6ZplkeR9PMbgonAznb3hUBb4JZ4\nHujuo4Fl21oFKAiHb20Urltes7giddeBOxbjjibOyUBVFdhlFRfcXQfROuQvh7Rn8+qlXP3adBYu\nXRF1HBGJk7svdPfb3P398KYOBCeqJ8LdwC7AfIJfNi9w95/0JTOzgWY23szGL16swkFka3Zv24Qm\nDXIZrSEHU1lxxfEsXAbG86CqCuweZrYqXEqB3Ssum9mqmmeWVNWqpBlXHN4Oz2/KSdc/GXUcEakG\nM9vDzG4JR366DkjUVK1HAZOA1gRDAd4dTsv+I+4+zN17uXuvkpKSBD21SObJzjL236GY979ZTHDq\nhKSgJRXHs3AZFs+Dtllgu3u2uxeGS4G758Rc/slBVTLLOSceSjeby7x67Rky4qWo44jINpjZjmZ2\ntZl9CdwFfAuYux/i7ncn6GnOBJ71wDRgJrBzgrYtUicduGMJi0o38OXC0qijSAJpZHPZpqevPpPs\nlfN4cNIavl2saWlFUtiXwKHAce6+v7vfBWxK8HN8CxwGYGYtgJ2AGQl+DpE65cCuwa88ozWaSEZR\ngS3b1KhBPvcO2Jt6DQu48qWv9BOWSOrqCywA3jGzB8zsMIJJweJmZk8CY4GdzGyumZ1tZoPMbFC4\nyhCgj5l9BrwNXOru6jwqUgMtG9dnpxYFGq4vw8Q70YzUYUft25MraMJVL0zlnjc+4w9HaQJPkVTj\n7s8Dz5tZQ4IxsAcDzc1sKPCcu78RxzZOqeL++cCRicgrIj84cMdiHvpgJgfe/E7UUSJ15n4dOXO/\nTlHHSAgV2BKX03p3YMQbE7jlzens3CyLw/faLepIIlIJd18DPAE8YWZNgV8BlwJVFtgiEo3+vTuw\nYm0Z5Zvr7q/EE79dzogPZ6nAlrrFzLit354cf+9HDBoxlsm7daVhfr2oY4nINrj7cmBYuIhIiupQ\n1JBbftUj6hiRGjl2Fle9MJVZS9bQsbhh1HFqTH2wJW49d+rE6d1yKC9ozSnXj4w6joiIiGSIg3YM\nTvbMlL7oKrClWoacezLN18xgclkLzfIoIiIiCdGhqCEdixqowJa667krTiFr3Qr+NW45pevLqn6A\niIiISBUO2rGEsdOXsr4s0SOM1j4V2FJtbVoU89TgY1hZlsW1L30edRwRERHJAAftVMK6sk2Mn7U8\n6ig1pgJbtss+XUr4/cFd+M+Eudww8pWo44iIiEia6925iLycLN79alHUUWpMBbZst/MP6kT2ynkM\nm1jK5zPmRh1HRERE0liDvBz26dQsI/phq8CW7dagfj3u7r8XZOfx69teZvPmzVFHEhERkTR20I4l\nfLNoNfNWrIs6So2owJYaOWa/PTi0aCWrG7Vj0G1PRh1HRERE0ljFcH2j07wVWxPNSI09eOkAuv9h\nKKPK2zDu81n07tYx6kgiIiKShnZo3ojWjetz5fNTGPJy+g6koAJbaiw7O5snBv+cUx+dyvVvz+HZ\nHduTl6MfR0RERKR6zIybTtqdD75JjRbsK7bzcSqwJSF67tSJW/vlM+ixCdz00mSu/uUeUUcSERGR\nNHTQjiXfdxWJ2vYW2GpmlIQ5ereW7NsCho+dyyOvaJZHERERqZtUYEtC/fO3+8CapVz9+gwWLl0R\ndRwRERGRWqcCWxKqTYtirjyiPZ7fhL7Xa1QRERERqXtUYEvC/e7EQ+lmc5lfrz1DRrwUdRwRERGR\nWqUCW5Li6avPJHvVPB77ajPfrVofdRwRERGRWqMCW5KiUYN8Xv7br7CcPC7+z6e4e9SRRERERGqF\nCmxJml3aNOVvx+7C6K8Xc8Hdz0YdR0RERKRWqMCWpDp1n/Y0WDmLF2Zn8fYnU6OOIyIiIpJ0KrAl\nqbKysnhi8LF42QbOHTGWtes3Rh1JREREJKlUYEvS9dypE6d3y6G8oBX9rh8ZdRwRERGRpFKBLbVi\nyLkn02LNdCZvbM6/35kQdRwRERGRpFGBLbXm2St+S9M8586PV1K6vizqOCIiIiJJoQJbak2bFsU8\ndM4BzF+xjr89MynqOCIiIiJJkdQC28yONrOvzGyamV1Wyf0Hm9lKM5sULlclM49E72cdmtF350a8\n+Nkihox4Oeo4IiIiIgmXtALbzLKBe4BjgG7AKWbWrZJV33f3nuFyXbLySOq45ld7k71iLg9OWsPU\nGXOjjiMiIiKSUMlswd4bmObuM9x9I/AUcEISn0/SRKMG+dxz2t6QnctvbnuZzZs3Rx1JREREJGGS\nWWC3AebEXJ8b3ralPmb2qZm9Zma7VrYhMxtoZuPNbPzixYuTkVVq2dH77cEhzVawulE7Bt32ZNRx\nRERERBIm6pMcJwLt3X134C7g+cpWcvdh7t7L3XuVlJTUakBJnocuO50GK2fy1uICpi0qjTqOiIiI\nSEIks8CeB7SLud42vO177r7K3VeHl18Fcs2sOImZJIVkZ2cz6rrTaNIonwuemsTGcnUVERERkfSX\nzAL7E6CrmXUyszygH/Bi7Apm1tLMLLy8d5hnaRIzSYppV9KYm/p2Z+r8VZxx69NRxxERERGpsaQV\n2O5eDvwBGAV8Afzb3aea2SAzGxSudjIwxcwmA3cC/dzdk5VJUtNRu7ak5dqZjFmazyOvfBB1HBER\nEZEaSWofbHd/1d13dPcu7n5DeNt97n5fePlud9/V3Xu4e293/zCZeSR1PfO3frBmGVe/PoPvlq2M\nOo5InWRmw81skZlN2cr9F8fMWzDFzDaZWbPazikikuqiPslRBAhmebzi8LZ4fhN+eb1GFRGJyAjg\n6K3d6e63VMxbAFwOvOfuy2ornIhIulCBLSnjnF8eTjfmMD+vHf94fFTUcUTqHHcfDcRbMJ8C6Nuw\niEglVGBLSnn6mrNonrOOJ6cZ361aH3UcEamEmTUgaOl+Ziv3a+4CEckUxRXHs3AZGM+DVGBLSmnU\nIJ8nLzia9eWb+PNTEzXLo0hq+gUwZmvdQzR3gYhkkCUVx7NwGRbPg1RgS8rpUtKI8/ZtyZgZy/n9\nbU9FHUdEfqof6h4iIrJVKrAlJZ1/1O40WDGD1xfm8/YnU6OOIyIhM2sMHAS8EHUWEZFUpQJbUlJ2\ndjZPXngcXraec0eMZe36jVFHEsl4ZvYkMBbYyczmmtnZW8xdAPBL4A13XxNNShGR1KcCW1JWj506\ncfouuZQXtKLf9SOjjiOS8dz9FHdv5e657t7W3R+KnbsgXGeEu/eLMqeISKpTgS0pbcigk2m+ehqf\nlrVk3PQlUccRERERqZIKbEl5r15/Fu2LGvCXpz+ldH1Z1HFEREREtkkFtqS84iYF3P6bPViwch39\nbn426jgiIiIi26QCW9LCzzo0pRtzmbq2EUNGvBx1HBEREZGtUoEtaeOpKweQvWIuD05aw9QZc6OO\nIyIiIlIpFdiSNho1yOfu0/bCs3P49e2vUVa+KepIIiIiIj+hAlvSyjH77ckxRctZ07A1Z96u/tgi\nIiKSelRgS9oZetlZdMtbxgdLG/D2F99FHUdERETkR1RgS9rJysri2St/y66tC7nw/yYxc9GqqCOJ\niIiIfE8FtqSl+rnZ3HLCjpSWlvKrf73O+jL1xxYREZHUoAJb0la3Di3YacXHLNnckLPvfwd3jzqS\niIiIiApsSW9P3fpXNk96jjFzN3DLa1OjjiMiIiKiAlvSW5MmTRhxyW9Z/emb3Dt6Ni9Mmhd1JBER\nEanjVGBL2jv00EM5fdc8spfO4KKnJzNh9rKoI4mIiEgdpgJbMsKN1w/hg3+cQZsm+fzukfF8/V1p\n1JFERESkjlKBLRmhXr16tGpWwF0n78LG9evo/+BHzFqyJupYIiIiUgepwJaMMmbU83wz7I+UrlnH\nqQ9+xLwV66KOJCIiInWMCmzJKAMHDuSAHl2Z/8TlrFi7gVMfGMeiVeujjiUiIiJ1iApsyShZWVmM\nGDGCrJXzqP/RQywq3cApD4xjvlqyRUREpJaowJaM065dO4YOHcrEN57h8NyvWLRqAycN/ZBpi3Ti\no4iIiCSfCmzJSP369eMf//gHF53+S546tzdlm5yT7xvLxG+XRx1NREREMpwKbMlYl1xyCe3bt6db\nq0KeOvtnNM7P5dQHPuKdrxZFHU1EREQymApsyWibNm3i5z//Obde91f+M6gPnUsacvaIT7j33Wls\n3uxRxxMREZEMpAJbMlp2djbdunXjvvvu4+PRb/Hvc/fl2O6tuPn1rxj46ARWriuLOqKIiIhkGBXY\nkvFuuOEGunfvztlnn83aVcu565Q9uOq4brz71SJOuPsDvliwKuqIIiIikkFUYEvGq1evHo8//jjL\nly/nnHPOAeCs/Tvx1MDerCvbxIn3jOGhD2aqy4iIiIgkhApsqRO6d+/OTTfdxPjx41mwYAEAvTo2\n4+U/HsD+OxQz5OXP6ffAOL75TkP5iYiISM2owJY6Y/DgwUyZMoXWrVt/f1tJQT0ePL0XN5+8O18u\nWMXRd7zPdS99rr7ZIiIist1UYEudkZWVRZMmTdi4cSP33nsv5eXlAJgZv+7VjncuOpjf7NWOhz+c\nyaH/fJcnP/6WTeo2IiIiItWkAlvqnDfeeIPzzz+fm2666Ue3FzWqx42/7M5Lf9ifziUNufzZzzji\ntvf4v0++ZUP5pojSioiISLpRgS11znHHHccpp5zCtddeyyeffPKT+3dr05h/n7sv9566Jw3qZXPp\nM59x4M3vMGz0dFZvKI8gsYiIiKQTc0+vn8B79erl48ePjzqGpLnly5fTo0cP8vPzmThxIg0bNqx0\nPXfng2lLGPrudD6cvpSGedkc070VffdoQ+/ORWRlWS0nl3RnZhPcvVfUOWqLjtkiks6295idk4ww\nIqmuadOmPPLIIxx22GFceuml3H333ZWuZ2Yc0LWEA7qWMHnOCh7/aDavfraQ/0yYS6vG9TmhZxuO\n2a0lu7VpTLaKbREREUEFttRhhxxyCLfffjsHHXRQXOv3aNeEHu2acO3xu/HmF9/x3MS5PPD+DO57\nbzpNG+Sy3w7FHLhjCQd0LaZV4/wkpxcREZFUpQJb6rQLLrjg+8tlZWXk5uZW+Zj8vGyO79Ga43u0\nZunqDbz/zRJGf7OY979ZwsufBmNst2pcn93aNGb3No3p3rYxu7VpTFHDPMzUyi0iIpLpklpgm9nR\nwB1ANvCgu/99i/stvP9YYC1whrtPTGYmkcqcd955zJ8/n+eee65aRXBRo3qcuEcbTtyjDe7OlwtL\nGTNtCZ/NW8ln81by5ufffb9uQf0cOhU3pENRQzoVNaB9UUNaFNajeUF9WhTWo3F+rgpwiZSZDQeO\nAxa5+25bWedg4F9ALrDE3eP7CUhEpA5JWoFtZtnAPcARwFzgEzN70d0/j1ntGKBruOwDDA3/FalV\nXbp0YejQoQwfPpyzzz57u7ZhZuzSqpBdWhV+f9uq9WVMnbeKzxesYvbSNcxcsobJc1bwyqfz2XKI\n7bycLEoa1aMwP5eC+jkU1s+hoH5wuVG9Hy7Xz80mLyeLejlZ3/8bLMHtOVlGdpaRZUZWlpFtRpbx\nw+Ws4HrFOt+va6jAlxHA3cDIyu40sybAvcDR7v6tmTWvxWwiImkjmS3YewPT3H0GgJk9BZwAxBbY\nJwAjPRjKZJyZNTGzVu6+IIm5RH7iwgsv5JVXXuH111/f7gK7MoX1c9m3SxH7din60e0byzczf8U6\nFpVuYFHpehat2sB3petZXLqBVevKKV1fxvwV6yndUErp+nJK15fX+qQ3sbW2/eh228rtsetX/uB4\n1t/a80ryuftoM+u4jVV+Czzr7t+G6y+qjVwiIukmmQV2G2BOzPW5/LR1urJ12gA/KrDNbCAwMLy6\nwcymJDZqxikGlkQdIoVt8/VRKy6g91BVavL6dEhkkFq2I5BrZu8CBcAd7v6T1m4ds+vk/x/tc91Q\nF/d5NzOLHWt0mLsPq+pBaXGSY7gjwwDMbHxdGkN2e+g12ja9PlXTa7Rtdfj1yQF+BhwG5ANjzWyc\nu38du1JdP2Zrn+sG7XPdsL37nMwCex7QLuZ62/C26q4jIiKpYS6w1N3XAGvMbDTQA/h62w8TEalb\nkjlV+idAVzPrZGZ5QD/gxS3WeREYYIHewEr1vxYRSVkvAPubWY6ZNSDo9vdFxJlERFJO0lqw3b3c\nzP4AjCIYpm+4u081s0Hh/fcBrxIM0TeNYJi+M+PYdJX9XkSvURX0+lRNr9G2ZeTrY2ZPAgcDxWY2\nF7iaYDg+3P0+d//CzF4HPgU2Ewy/WlX/6ox8raqgfa4btM91w3btswUDeIiIiIiISCIks4uIiIiI\niEidowJbRERERCSB0qrANrOjzewrM5tmZpdFnSfVmNlwM1tUB8ecjYuZtTOzd8zsczObamYXRJ0p\nlZhZfTP72Mwmh6/PtVFnSkVmlm1m/zOzl6POkkqqOj6HJ7PfGd7/qZntGUXORIpjn3c2s7FmtsHM\nLooiY6LFsc+nhn/fz8zsQzPrEUXORIpjn08I93mSmY03s/2jyJlI8dZbZraXmZWb2cm1mS/R4vgb\nH2xmK8O/8SQzu6rKjbp7WiwEJ0pOBzoDecBkoFvUuVJpAQ4E9gSmRJ0lFRegFbBneLmAYGgxvYd+\neH0MaBRezgU+AnpHnSvVFuDPwBPAy1FnSZUlnuMzwQntr4Xvs97AR1HnroV9bg7sBdwAXBR15lra\n5z5A0/DyMXXk79yIH85p2x34Murcyd7nmPX+SzBgxclR507y3/jg6h7z06kF+/up1919I1Ax9bqE\n3H00sCzqHKnK3Re4+8TwcinB8GJtok2VOjywOryaGy46CzqGmbUFfg48GHWWFBPP8fkEYGT4PhsH\nNDGzVrUdNIGq3Gd3X+TunwBlUQRMgnj2+UN3Xx5eHUcwv0U6i2efV3tYhQENSf/jZrz11h+BZ4BF\ntRkuCZJSX6ZTgb21adVFqs3MOgJ7ELTSSijs/jCJ4ID5prvr9fmxfwGXEAxRJz+I5/icacfwTNuf\neFR3n88m+NUincW1z2b2SzP7EngFOKuWsiVLlftsZm2AXwJDazFXssT7vu4TdgV6zcx2rWqj6VRg\niySEmTUi+NY92N1XRZ0nlbj7JnfvSdDqtLeZ7RZ1plRhZscBi9x9QtRZRFKdmR1CUGBfGnWW2uDu\nz7n7zsCJwJCo89SCfwGXuntdaWyYCLR3992Bu4Dnq3pAOhXYmlZdaszMcgmK68fd/dmo86Qqd18B\nvAMcHXWWFLIfcLyZzSL4CfFQM3ss2kgpI57jc6YdwzNtf+IR1z6b2e4E3ahOcPeltZQtWar1dw67\nanY2s+JkB0uiePa5F/BUeDw8GbjXzE6snXgJV+X+uvuqii6U7v4qkFvV3zidCux4pl4X2SozM+Ah\n4At3vy3qPKnGzErMrEl4OR/4//buP9bquo7j+PMFleaQGlJpWlL2BzILBowcSLsyXc02ooURAWZR\nG7QK3HJlK+ZabaxcZDZHRYrWVikoslqKk/yVy4i4V8ycGtZaoZksp/0EefXH53Pn8XbPvefioXMP\n9/XYzi7n++vz/nJ239/3+Xw/9/s5H3i4s1GNHrYvs32a7SmU/LPT9vIOhzVatJKftwMX1aeJnA08\nY9HE/7wAAAZHSURBVHv//zvQNhqL16Rhz1nSG4GbgBW2H+lAjO3Wyjm/pV5fqE/HOQ7o5i8Ww56z\n7TfZnlLz4Rbg47aH7dUdpVr5jE9u+IznUOrnIT/jozZVeru5ydTrHQ5rVNEg0xzb/m5noxpV5gEr\ngL11nDHA5+q30ShPWblO0nhK8rjBdh5FF8Nqlp8lrarrN1KeNHAB8BjwD+DDnYq3HVo5Z0knA78C\nJgKHJa2lPJ2gK4emtfg5rwNOovRoAhyyPbtTMb9ULZ7z+yhfHg8C/wSWNPzRY9dp8ZyPGS2e72Jg\ntaRDlM/4A8N9xpkqPSIiIiKijbppiEhERERExKiXAjsiIiIioo1SYEdEREREtFEK7IiIiIiINkqB\nHRERERHRRimwo2tIWiTJkqYeheNOa7Luckl/ktQr6UFJC5tst1DSZ9sZV0REJ0k6qea+XklPNOTC\nXkn3HYX2Zki6YIT79NTrwkcHHMeSPj3Mvk1z/zD7vV7SliONOcaGFNjRTZYC99af7bQIGCrJbqjT\nh18IXCPpRb83kl5me7vt9W2OKyKiY2w/bXtGzX8bqbmwvuYehSZnUJ6VPlIPAu9veL8U6Gthv+Fy\n//+o+f7PthfXRUcacxzjUmBHV5A0ATgHWEmZZal/+SmS7m7oYZ4vabykzfX9XkmX1G3PkHSrpN2S\n7pE0VdJcYCHw1XqMM5rFYPu3wCHKRD6bJW2UdD/wFUkXS/pmbed1km6W1Fdfc+vy5ZJ+Wdv5Vp3Q\nJSKi60h6rv7skXSXpFsk7ZO0XtKymuv29udUlZlit0raVV/zBhzvFcAXgSU1Ry6RNEnSNkkPSPqF\nyhTsg/kDcHzNvQLeBfy04dgfq2321RhOGCz3S7pT0uy6z2SVacCp+X27pJ3AHZKm1OvLYDE/Kuk1\ndb9xkh7rfx9jS9fM5Bhj3nuAW20/IulpSbNs7wY+CNxm+8u1YD2B0qNwqu2zAFSn/wa+Dayy/aik\ntwNX214gaTvwY9tbhgqg7nMYeKouOg2Ya/t5SRc3bPoN4C7b760xTZB0JrAEmGf7oKSrgWXA9S/1\nPyYiosOmA2cCB4B9wCbbcyStAT4JrAWupPSA36synfptdR8AbP9H0jpgtu1PAEi6Cthje5GkBZR8\nOaNJDFsodxn3AL8G/t2w7ibb36nH/BKw0vZVA3N/qc2bmgm8zfYBSVOGiHkqJbd/HTgP6LP91OCH\njGNZCuzoFkspCRrgh/X9bmAXZdjGy4Fttnsl7QPeXJPzT4AdtQd8LnBjQxI9rsW2L5G0HHiWOgVu\nPcaNtp8fZPsFwEUAdf0zklYAs4Bddd9XAn9p+ewjIkavXbb3A0j6HbCjLt8LnFv/fR4wrSH/TpQ0\nwfZzQxz3HMo05NjeqTImfGKTqeZvAH4ETAV+QMn3/c6qhfWrgQmU4n6kbrd9oIXtrgFuoRTYHwGu\nPYK24hiQAjtGPUmTKEXrWyUZGA9Y0qW275b0DuDdwGZJX7N9vaTpwDuBVZSxeWuBv9WxhCO1wfYV\ngyz/+0hOA7jO9mVH0H5ExGjW2Ft8uOH9YV6oM8YBZ9v+19EIwPYTkg4C5wNreHGBvRlYZLuv3m3s\naXKYQ7wwdPb4Aetayve2/yjpydrjPofSmx1jUMZgRzdYDHzP9um2p9h+A/A4MF/S6cCT9fbfJmCm\npMnAONtbgc8DM2uPx+OSLgRQMb0e/1ngxDbGewewurYzXtKr6rLFkl5bl0+qsUdEjAU7KMNFgPL0\njUG2GZiL76EWqJJ6gL826b3utw74zCB3Fk8E9tc7nY0F78D2fk+50wjlutOKwa4fm4Dv0/wuZ4wB\nKbCjGywFbh6wbGtd3gP0SdpDGeN8JXAqcKekXkqS6+81XgaslNQH/IYyrhvKkJNLJe3REH/kOAJr\ngHMl7aUMY5lm+yFKsb9D0gPA7cApbWgrIqIbfAqYXf9g8SHK3cWBfkYZRtIraQlwOTCr5sz1wIeG\nasD2fba3DbLqC8D9wM+BhxuWD8z9VwCr6/VkcovnNTBmgO2UoSgZHjKGyXanY4iIiIg4JtQnkWyw\nPb/TsUTnZAx2RERERBuoTDi2moy9HvPSgx0RERER0UYZgx0RERER0UYpsCMiIiIi2igFdkRERERE\nG6XAjoiIiIhooxTYERERERFt9F9A4zHsrd3W/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11cdb9208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "axes[0].plot([0, strike], [strike, 0], 'k--')\n",
    "axes[0].plot(ps, v[:-1])\n",
    "axes[0].set_xlim(0, strike*2)\n",
    "axes[0].set_xticks(np.linspace(0, 4, 5, endpoint=True))\n",
    "axes[0].set_ylim(0, strike)\n",
    "axes[0].set_yticks(np.linspace(0, 2, 5, endpoint=True))\n",
    "axes[0].set_xlabel('Asset Price')\n",
    "axes[0].set_ylabel('Premium')\n",
    "axes[0].set_title('Put Option Value')\n",
    "\n",
    "axes[1].plot(np.linspace(0, T, N), max_exercise_price)\n",
    "axes[1].set_xlim(0, T)\n",
    "axes[1].set_ylim(1.6, strike)\n",
    "axes[1].set_xlabel('Time to Maturity')\n",
    "axes[1].set_ylabel('Asset Price')\n",
    "axes[1].set_title('Put Option Optimal Exercise Boundary')\n",
    "axes[1].tick_params(right='on')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The figure in the right panel looks different from\n",
    "the corresponding figure in Miranda and Fackler (Figure 7.4(b), p.180)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}