9879HW10_Chen_Weiyi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>\n",
    "### Student: Chen Weiyi, weiyi.alan.chen@gmail.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MTH9879 Homework 10\n",
    "\n",
    "Assigned: April 28, 2015.\n",
    "Due: May 12, 2015 by 6pm. \n",
    "\n",
    "Late homework **will not be accepted**.\n",
    "\n",
    "$$\n",
    "\\newcommand{\\supp}{\\mathrm{supp}}\n",
    "\\newcommand{\\E}{\\mathbb{E}}\n",
    "\\newcommand{\\Eof}[1]{\\mathbb{E}\\left[ #1 \\right]}\n",
    "\\def\\Cov{{ \\mbox{Cov} }}\n",
    "\\def\\Var{{ \\mbox{Var} }}\n",
    "\\newcommand{\\1}{\\mathbf{1} }\n",
    "\\newcommand{\\PP}{\\mathbb{P} }\n",
    "%\\newcommand{\\Pr}{\\mathrm{Pr} }\n",
    "\\newcommand{\\QQ}{\\mathbb{Q} }\n",
    "\\newcommand{\\RR}{\\mathbb{R} }\n",
    "\\newcommand{\\DD}{\\mathbb{D} }\n",
    "\\newcommand{\\HH}{\\mathbb{H} }\n",
    "\\newcommand{\\spn}{\\mathrm{span} }\n",
    "\\newcommand{\\cov}{\\mathrm{cov} }\n",
    "\\newcommand{\\sgn}{\\mathrm{sgn} }\n",
    "\\newcommand{\\HS}{\\mathcal{L}_{\\mathrm{HS}} }\n",
    "%\\newcommand{\\HS}{\\mathrm{HS} }\n",
    "\\newcommand{\\trace}{\\mathrm{trace} }\n",
    "\\newcommand{\\LL}{\\mathcal{L} }\n",
    "%\\newcommand{\\LL}{\\mathrm{L} }\n",
    "\\newcommand{\\s}{\\mathcal{S} }\n",
    "\\newcommand{\\ee}{\\mathcal{E} }\n",
    "\\newcommand{\\ff}{\\mathcal{F} }\n",
    "\\newcommand{\\hh}{\\mathcal{H} }\n",
    "\\newcommand{\\bb}{\\mathcal{B} }\n",
    "\\newcommand{\\dd}{\\mathcal{D} }\n",
    "\\newcommand{\\g}{\\mathcal{G} }\n",
    "\\newcommand{\\p}{\\partial}\n",
    "\\newcommand{\\half}{\\frac{1}{2} }\n",
    "\\newcommand{\\T}{\\mathcal{T} }\n",
    "\\newcommand{\\bi}{\\begin{itemize}}\n",
    "\\newcommand{\\ei}{\\end{itemize}}\n",
    "\\newcommand{\\beq}{\\begin{equation}}\n",
    "\\newcommand{\\eeq}{\\end{equation}}\n",
    "\\newcommand{\\beas}{\\begin{eqnarray*}}\n",
    "\\newcommand{\\eeas}{\\end{eqnarray*}}\n",
    "\\newcommand{\\cO}{\\mathcal{O}}\n",
    "\\newcommand{\\cC}{\\mathcal{C}}\n",
    "\\newcommand{\\cF}{\\mathcal{F}}\n",
    "\\newcommand{\\cL}{\\mathcal{L}}\n",
    "\\newcommand{\\BS}{\\text{BS}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color = \"red\">Homework is to be done by each student individually.  To receive full credit, you must email a completed copy of this iPython notebook to Yu Gan (yugan323@gmail.com), Fubo Shi (fubo.shi.baruch@gmail.com), and Tai-Ho Wang (tai-ho.wang@baruch.cuny.edu) by the due date and time.  All R-code must run correctly and solutions must be written up neatly in Markdown/LaTeX format.\n",
    "\n",
    "<font color=\"blue\">If you encounter problems with Markdown/LaTeX or iPython notebook, please contact your TAs Yu Gan and/or Fubo Shi.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up R environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical solution of Fredholm equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. (12 points)\n",
    "Consider the Fredholm equation (7) in Lecture 10:\n",
    "\n",
    "$$\n",
    "\\int_0^T{G}(|t-s|)\\,dx_s =\\text{const.}\\qquad\\text{for all $t\\in [0,T]$.}\n",
    "$$\n",
    "\n",
    "Define $t_i:={i}\\frac T{N}$. One way to discretize this integral equation is as follows:\n",
    "\n",
    "$$\n",
    "\\sum_{j=1}^N G_{ij}\\,v_j=\\text{const.}\n",
    "$$\n",
    "\n",
    "where $v_j=x_{t_{j}}-x_{t_{j-1}}$ and\n",
    "\n",
    "<a name=\"eq:1\"></a>(1)\n",
    "$$\n",
    "G_{ij}=\n",
    "\\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt.\n",
    "$$\n",
    "\n",
    "(a) Derive an explicit expression for the $G_{ij}$ in the case\n",
    "$$\n",
    "G(\\tau)=\\frac {1}{(1+\\tau)^2}.\n",
    "$$\n",
    "You will need to compute separate expressions for the cases $i \\neq j$ and $i=j$.\n",
    "\n",
    "(b) With $N=512$, use matrix algebra to solve equation [(1)](#eq:1) numerically with this decay kernel.\n",
    "\n",
    "(c) Plot your results - you should end up with something like the top figure in Figure 2 of Lecture 10.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>\n",
    "\n",
    "### Solution\n",
    "**(a)** When $i < j$ , substitute the expression of $G(\\tau)$ to \n",
    "\n",
    "$$\n",
    "G_{ij} = \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt.\n",
    "$$\n",
    "\n",
    "we derive\n",
    "\n",
    "\\begin{equation}\n",
    "G_{ij} = \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,\\frac{1}{(1+s-t)^2}\\,ds\\,dt \n",
    "= \\ln\\left(\\frac{1+t_{j-1}-t_{i-1}}{1+t_{j-1}-t_i} \\cdot \\frac{1+t_j-t_i}{1+t_j-t_{i-1}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "guaranteed positive inside $\\ln$, similarly when $i > j$,\n",
    "\n",
    "\\begin{equation}\n",
    "G_{ij} = \\ln\\left(\\frac{1+t_{i-1}-t_{j-1}}{1+t_{i-1}-t_j} \\cdot \\frac{1+t_i-t_j}{1+t_i-t_{j-1}}\\right)\n",
    "\\end{equation}\n",
    "\n",
    "When $i = j$, we need to separate the inner integral: \n",
    "\n",
    "\\begin{equation}\n",
    "G_{ij} = \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{i-1}}^{t}\\,\\frac{1}{(1+t-s)^2}\\,ds\\,dt + \\int_{t_{i-1}}^{t_i}\\,\\int_{t}^{t_i}\\,\\frac{1}{(1+s-t)^2}\\,ds\\,dt = 2(t_i-t_{i-1}) - 2\\ln(1+t_i-t_{i-1})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAD8GlDQ1BJQ0MgUHJvZmlsZQAAOI2N\nVd1v21QUP4lvXKQWP6Cxjg4Vi69VU1u5GxqtxgZJk6XpQhq5zdgqpMl1bhpT1za2021Vn/YCbwz4\nA4CyBx6QeEIaDMT2su0BtElTQRXVJKQ9dNpAaJP2gqpwrq9Tu13GuJGvfznndz7v0TVAx1ea45hJ\nGWDe8l01n5GPn5iWO1YhCc9BJ/RAp6Z7TrpcLgIuxoVH1sNfIcHeNwfa6/9zdVappwMknkJsVz19\nHvFpgJSpO64PIN5G+fAp30Hc8TziHS4miFhheJbjLMMzHB8POFPqKGKWi6TXtSriJcT9MzH5bAzz\nHIK1I08t6hq6zHpRdu2aYdJYuk9Q/881bzZa8Xrx6fLmJo/iu4/VXnfH1BB/rmu5ScQvI77m+Bkm\nfxXxvcZcJY14L0DymZp7pML5yTcW61PvIN6JuGr4halQvmjNlCa4bXJ5zj6qhpxrujeKPYMXEd+q\n00KR5yNAlWZzrF+Ie+uNsdC/MO4tTOZafhbroyXuR3Df08bLiHsQf+ja6gTPWVimZl7l/oUrjl8O\ncxDWLbNU5D6JRL2gxkDu16fGuC054OMhclsyXTOOFEL+kmMGs4i5kfNuQ62EnBuam8tzP+Q+tSqh\nz9SuqpZlvR1EfBiOJTSgYMMM7jpYsAEyqJCHDL4dcFFTAwNMlFDUUpQYiadhDmXteeWAw3HEmA2s\n15k1RmnP4RHuhBybdBOF7MfnICmSQ2SYjIBM3iRvkcMki9IRcnDTthyLz2Ld2fTzPjTQK+Mdg8y5\nnkZfFO+se9LQr3/09xZr+5GcaSufeAfAww60mAPx+q8u/bAr8rFCLrx7s+vqEkw8qb+p26n11Aru\nq6m1iJH6PbWGv1VIY25mkNE8PkaQhxfLIF7DZXx80HD/A3l2jLclYs061xNpWCfoB6WHJTjbH0mV\n35Q/lRXlC+W8cndbl9t2SfhU+Fb4UfhO+F74GWThknBZ+Em4InwjXIyd1ePnY/Psg3pb1TJNu15T\nMKWMtFt6ScpKL0ivSMXIn9QtDUlj0h7U7N48t3i8eC0GnMC91dX2sTivgloDTgUVeEGHLTizbf5D\na9JLhkhh29QOs1luMcScmBXTIIt7xRFxSBxnuJWfuAd1I7jntkyd/pgKaIwVr3MgmDo2q8x6IdB5\nQH162mcX7ajtnHGN2bov71OU1+U0fqqoXLD0wX5ZM005UHmySz3qLtDqILDvIL+iH6jB9y2x83ok\n898GOPQX3lk3Itl0A+BrD6D7tUjWh3fis58BXDigN9yF8M5PJH4B8Gr79/F/XRm8m241mw/wvur4\nBGDj42bzn+Vmc+NL9L8GcMn8F1kAcXgSteGGAAAysElEQVR4Ae3dC5AcVdnw8WeTEELuhCRgLlgJ\nqIGES0SRoIWA4MvFUiGAEoQQ9ANBYoVriVpyS4kXIMVVBZRXgcAHUVAhSiIERQGFKpBICHwhRLkF\nQhJyQXLZbH/9dPYMs5vZ2Wd2e/r06fl31WZmuk+fy++czDOn+8xuUxRvwoYAAggggAACmQr0yLQ0\nCkMAAQQQQACBRIAAzEBAAAEEEEDAgwAB2AM6RSKAAAIIIEAAZgwggAACCCDgQYAA7AGdIhFAAAEE\nECAAMwYQQAABBBDwIEAA9oBOkQgggAACCBCAGQMIIIAAAgh4ECAAe0CnSAQQQAABBAjAjAEEEEAA\nAQQ8CBCAPaBTJAIIIIAAAgRgxgACCCCAAAIeBAjAHtApEgEEEEAAAQIwYwABBBBAAAEPAgRgD+gU\niQACCCCAAAGYMYAAAggggIAHAQKwB3SKRAABBBBAgADMGEAAAQQQQMCDAAHYAzpFIoAAAgggQABm\nDCCAAAIIIOBBgADsAZ0iEUAAAQQQIAAzBhBAAAEEEPAgQAD2gE6RCCCAAAIIEIAZAwgggAACCHgQ\nIAB7QKdIBBBAAAEECMCMAQQQQAABBDwIEIA9oFMkAggggAACBGDGAAIIIIAAAh4ECMAe0CkSAQQQ\nQAABAjBjAAEEEEAAAQ8CBGAP6BSJAAIIIIAAAZgxgAACCCCAgAcBArAHdIpEAAEEEECAAMwYQAAB\nBBBAwIMAAdgDOkUigAACCCBAAGYMIIAAAggg4EGAAOwBnSIRQAABBBAgADMGEEAAAQQQ8CBAAPaA\nTpEIIIAAAggQgBkDCCCAAAIIeBAgAHtAp0gEEEAAAQQIwIwBBBBAAAEEPAgQgD2gUyQCCCCAAAIE\nYMYAAggggAACHgQIwB7QKRIBBBBAAAECMGMAAQQQQAABDwIEYA/oFIkAAggggAABmDGAAAIIIICA\nBwECsAd0ikQAAQQQQIAAzBhAAAEEEEDAgwAB2AM6RSKAAAIIIEAAZgwggAACCCDgQYAA7AGdIhFA\nAAEEECAAMwYQQAABBBDwIEAA9oBOkQgggAACCBCAGQMIIIAAAgh4ECAAe0CnSAQQQAABBHo1EsGc\nOXOkubm5kZpMWxFAAAEEqggMHz5cDj300Cop6neoKYq3+mWfn5x//etfy1VXXSVTp07NT6WoCQII\nIICAV4Frr71W7rjjDtl3330zr0fDzIB15nvKKafIGWeckTkyBSKAAAII5FPgxRdflJaWFi+V4x6w\nF3YKRQABBBBodAECcKOPANqPAAIIIOBFgADshZ1CEUAAAQQaXYAA3OgjgPYjgAACCHgRIAB7YadQ\nBBBAAIFGFyAAN/oIoP0IIIAAAl4ECMBe2CkUAQQQQKDRBQjAjT4CaD8CCCCAgBcBArAXdgpFAAEE\nEGh0AQJwo48A2o8AAggg4EWAAGxgj9asMaQiCQIIIICAD4FoyxaJ1q/3UXS3yiQAG/hapv0fiVat\nMqQkCQIIIIBA5gL/b4m0XHBR5sV2t0ACcHcFOR8BBBBAwK+A/lG/pia/dehC6QRgC1qAHWtpFmkQ\nQAABBPwJEID92VMyAggggEAaAoH+WXsCsKXzdQYcaAdbmkcaBBBAIHiBAK9UEoCDH3U0AAEEEEAg\nRAECsKXXmAFblEiDAAII+BFgEZYfd0pFAAEEEEAgRAFmwNZe4x6wVYp0CCCAQLYCzICz9aY0BBBA\nAAEEQhZgBmzpvQBX11maRRoEEECgEAKBXqEkABdi9NEIBBBAoMEFApwoEYAtY5ZV0BYl0iCAAAII\n1CBAAK4Bi6QIIIAAAjkUSBZh5bBenVSJANwJUHKYGbBFiTQIIIAAAjUIEIBrwCIpAggggEAOBfga\nUg47Jc0qBbrKLk0C8kIAAQQQSE+AGXB6luSEAAIIIOBDINAJEgHYMlgCXN5uaRZpEEAAAQT8CRCA\n/dlTMgIIIIBAWgIBTpQIwJbOZxW0RYk0CCCAgB8BFmH5cadUBBBAAAEEQhRgBmzpNWbAFiXSIIAA\nAn4EmAH7cadUBBBAAAEEQhRgBmzttUCXuVubRzoEEEAgWIFA358JwMGOOCqOAAIIIBCyAAHY0nsB\nLm+3NIs0CCCAQGEEAnyfJgAXZvTREAQQQKBBBViEVeCOZxV0gTuXpiGAAAJ+BJgB+3GnVAQQQACB\ntASYAaclmcN8mAHnsFOoEgIIIBC2ADPgsPuP2iOAAAII8DWkgo+BQDu44L1C8xBAAIFgBZgBB9t1\nVBwBBBBAoCTA15BKFMV6EmDHFqsDaA0CCCBQPAFmwMXrU1qEAAIINJYAq6AL3N+sgi5w59I0BBAI\nXoAAHHwX0gAEEEAAAQQyE+AStIWaGbBFiTQIIICAH4FAv6VCAPYzXCgVAQQQQKDBBQjA1gEQ6Ccs\na/NIhwACCAQtEOC3VQjAQY84Ko8AAgggEKoAAdjScwF+srI0izQIIIBAIQRYBV2IbqQRCCCAAAKh\nCRCAQ+uxGurLKugasEiKAAIIIGAR4BK0RYk0CCCAAAL5FQh0kSwB2DKkmAFblEiDAAIIIFCDAAG4\nBiySIoAAAgjkVCDAxbIEYOtYCvQSh7V5pEMAAQQQyFaAAJytN6UhgAACCKQtwCrotEVzlF+AlzZy\npEdVEEAAAQQqCDADroDCLgQQQACBgASYAQfUWbVWNVkFXetJpEcAAQQQyEQg0DU6zIAzGR0UggAC\nCCCAQFsBAnBbj8qv+B5wZRf2IoAAAnkRCHCtDgE4L4OHeiCAAAIINJQAAdja3YHeY7A2j3QIIIBA\nqAJR/P7cxAw41O6j3ggggAACCGQrwAzY4h3gJytLs0iDAAIIFEKAryEVohtpBAIIIIBAaAKB3iJk\nBmwZaKyCtiiRBgEEEECgBgECcA1YJEUAAQQQyKlAgLcKCcCWscQM2KJEGgQQQACBGgQIwDVgkRQB\nBBBAIIcCLMJKp1Oam5tl9erV6WSWZi6B3uRPk4C8EEAAAQTSE8jFDHjTpk3y7W9/W0aPHi29e/eW\nIUOGSL9+/WTChAly6623ptdackIAAQQQKJ5AoDPgXnnoienTp8vy5cvlgQcekLFjxybBd+3atbJo\n0SKZMWOGbNiwQc4880x/VQ3w5r4/LEpGAAEEELAI5GIGPG/ePPnZz34me++9t/Tv3z/5lWKDBg2S\nSZMmyTXXXCP33XefpS2kQQABBBBoRIFAbxHmIgDrpeYFCxZUHDb333+/DBs2rOKxzHayCjozagpC\nAAEEuiQQ4JXKXFyCvuyyy2TKlCkya9Ys2W233WTgwIGyZs0aef7550UXZc2dO7dL/cFJCCCAAAII\n5FUgFwF44sSJ8vTTT8vjjz8uy5YtS+4H66xX7/sedNBB5r9ycdNNN8ns2bMrWi9ZsiQJ7l//+tcr\nHq+6kxlwVR4OIoAAAl4FWITVPf4+ffrIIYcc0q1MTj/9dNGfSts555yTBPZKx9iHAAIIIIBA1gK5\nuAecdaO7VF6gN/m71FZOQgABBEISYAbc9d666qqrZPPmzR1mMG7cOPniF7/Y4XEOIIAAAgggEJpA\nLu4B633f66+/XqZOnZp8B7g9Yi5WQbevFK8RQAABBPIhEOgVylwE4Ouuu05aWlqSnxtuuCEfHUot\nEEAAAQTCEQjwa0i5uQf8wx/+UPS3X61fvz5/Hd4UVynQT1j5w6RGCCCAAAIqkIsZsFZEfwPWHXfc\noU/ZEEAAAQQQsAtEcVJmwHavoFLyPeCguovKIoAAAiEI5OYSdAhY1BEBBBBAIIcCgX4NiQBsHUvc\nA7ZKkQ4BBBBAwCBAADYgkQQBBBBAIMcCgU6QCMCWMRXgzX1Ls0iDAAIIIOBPgADsz56SEUAAAQTS\nEghwokQAtnQ+q6AtSqRBAAEE/AiwCMuPO6UigAACCCAQogAzYEuvMQO2KJEGAQQQ8CPADNiPO6Ui\ngAACCCAQogAzYGuvBbrM3do80iGAAALBCgT6/kwADnbEUXEEEEAAgZAFCMCW3gtwebulWaRBAAEE\nCiMQ4Ps0Abgwo4+GIIAAAg0qwCKsAnc8q6AL3Lk0DQEEEPAjwAzYjzulIoAAAgikJcAMOC3JHObD\nDDiHnUKVEEAAgbAFmAGH3X/UHgEEEECAryEVfAwE2sEF7xWahwACCAQrwAw42K6j4ggggAACJQG+\nhlSiKNaTADu2WB1AaxBAAIHiCTADLl6f0iIEEECgsQSSVdDhNZkAbOkzVkFblEiDAAII+BHga0h+\n3CkVAQQQQACBEAWYAVt6jRmwRYk0CCCAgB+BQL+lQgD2M1woFQEEEECgwQUIwNYBEOgnLGvzSIcA\nAggELRDgt1UIwEGPOCqPAAIIIBCqAAHY0nMBfrKyNIs0CCCAQCEEWAVdiG6kEQgggAACoQkQgEPr\nsRrqyyroGrBIigACCCBgEeAStEWJNAgggAAC+RUIdJEsAdgypJgBW5RIgwACCCBQgwABuAYskiKA\nAAII5FQgwMWyBGDrWAr0Eoe1eaRDAAEEEMhWgACcrTelIYAAAgikLcAq6LRFc5RfgJc2cqRHVRBA\nAAEEKggwA66Awi4EEEAAgYAEmAEH1Fm1VpVV0LWKkR4BBBDITiDQNTrMgLMbIpSEAAIIIIBASYAA\nXKKo8oQZcBUcDiGAAAI5EAhwrQ4BOAfjhioggAACCDSeAAHY2ueB3mOwNo90CCCAQLACLMIKtuuo\nOAIIIIAAApkLMAO2kAd4b8HSLNIggAAChRBgBlyIbqQRCCCAAAKhCQR6i5AZsGGgNcUz4CjQDjY0\njyQIIIAAAh4ECMAe0CkSAQQQQCBlgQBvFRKALWNAO5YZsEWKNAgggAACRgECsBGKZAgggAACORVg\nEVZOOyatajEDTkuSfBBAAAEEYgFmwAwDBBBAAIGwBZgBh91/VWsf4M39qu3hIAIIIICAdwFmwN67\ngAoggAACCHRLINBbhARgS6+zCtqiRBoEEEDAn0CAVyoJwP6GCyUjgAACCDSwAAHY0vnMgC1KpEEA\nAQT8CLAIy487pSKAAAIIIBCiADNga68FepPf2jzSIYAAAsEKMAMOtuuoOAIIIIAAApkLMAO2kAe4\nus7SLNIggAAChRAI9AolAbgQo49GIIAAAg0uEOBEiQBsGbOsgrYokQYBBBBAoAYBAnANWCRFAAEE\nEMihAIuwctgpaVWJGXBakuSDAAIIINAqwAyYoYAAAgggELYAM+Cw+6/T2ge6yq7TdpEAAQQQQMCL\nADNgL+wUigACCCCQmkCgEyQCsGUEBLi83dIs0iCAAAII+BMgAPuzp2QEEEAAgbQEApwoEYAtnZ+s\ngrYkJA0CCCCAQOYCLMLKnJwCEUAAAQQQCFaAGbCl6/gesEWJNAgggIAfgUBnwL060tq0aZM8/fTT\n8txzz8nSpUtl9913l/3220/22GMP6dWrw9M6yo79CCCAAAIIIFAmsM0MOIo/Sdx5550yfvx4Oeec\nc+Thhx+W9957T+bOnSsnnXSSjBo1Sq666irRAN1QW6DL3Buqj2gsAgg0pkCg789tprIbNmyQyZMn\ny+GHHy6PP/64DB06dJvOXL16tdx4441y6KGHyt133y0jRozYJg07EEAAAQQQQKC6QJsArJeW77rr\nLhkwYECHZ+24447yne98R84991xpaWnpMF2hDgS4vL1Q/jQGAQQQ6EwgwPfpNpegNQCXB9/XXntN\nNm/eLOvWrZNrrrlG7r333hLBDjvsIP369Su95gkCCCCAAAJeBIq2COuxxx5LLkW/8MILctlll8lT\nTz2V3PddtWqVfPWrX/Vi7K1QVkF7o6dgBBBAoKgCbWbA5Y28/fbb5ZZbbpGdd945udf7q1/9SnTf\nnDlzypPxHAEEEEAAAb8Cgc6AOwzAa9askWHDhsmjjz4qw4cPlwkTJsjGjRtl4MCBfqF9lM4M2Ic6\nZSKAAAKFFmizCKu8pUcffbTMmDFDtmzZIqeeeqosWrRIpk6dmizAKk9Xr+crVqwQXfDFd47rJUy+\nCCCAQEEEAv0aUocz4ClTpsjVV18tl156qXzrW9+S3r17y09+8hM5+eSTU++xU045RRYvXpzkq/ec\nNfiPHj1adtllFzn77LOThWCpF1prhoF2cK3NJD0CCCCAQDYCbQKw/sKNBQsWlEr+7Gc/KyeccIL0\n6NEj+U1YhxxySOmYXpp+4403Sq+78+Rf//qXvPvuu0kWV1xxhYwbN05ef/110YVgy5YtE93HhgAC\nCCCAQIcCAX4Nqc0laJ3lPvLII3Lttdcmgfeoo46SQYMGldr76quvyt///vfkK0l6T/hHP/pR6Vha\nTx588EF58cUXk69DDRkyRGbOnJl85/h73/tep0XcdNNNMnv27IrplixZImPGjKl4rNOdAXZsp20i\nAQIIIICAV4E2Abhnz57JJWedfV588cUyffr05Jdt6EpoDb6DBw+WSZMmyQ033CB77bVXqhXX2a7+\nVq0DDjhAVq5cWfo+8sKFC2XixImmsk4//XTRn0qb/lrN5cuXVzrEPgQQQACBkAUCXQXdJgA7fw2E\nN998c/KjQcv9MQZdDV2PTX/H9O9//3u5/PLLRVdf9+nTJ/l91JdcckkS7B966KF6FGvPk1XQditS\nIoAAAlkLRHGBAV6prBiAy+10IZT+1HM777zzRH9009++tXbt2uT5EUccIeeff770798/ec0/CCCA\nAAIIFEWgwwCsf/Hotttu26adujCrHvd+XUEjR44U/dFNL0fnYmMGnItuoBIIIIBARYFAv6XSYQA+\n5phjZP/990/aqn+iUO8L6++D1oVZbAgggAACCCDQPYEOA/DYsWNFf8o3fX3llVfKwQcfXL67MZ4H\n+gmrMTqHViKAQMMLBHgPuM33gDvrwJdffjlZJNVZOo4jgAACCCCAQHWBDmfAOtPVP8DgNv0lHa+8\n8kqyOtnta5jHAD9ZNUzf0FAEEECgSF9D0t6cPHly8p1f17P6O5n1ErT+gQY2BBBAAAEEciNQtACs\nvzWqy785Kje9klJFmuJ8uAecEibZIIAAAgioQE33gCFDAAEEEEAgdwKBTpAIwJaRxPeALUqkQQAB\nBBCoQYAAXAMWSRFAAAEEcioQ4GJZArB1LAV6icPaPNIhgAACCGQrQADO1pvSEEAAAQTSFgh0FTQB\n2DIQAry0YWkWaRBAAAEE/AkQgP3ZUzICCCCAQBoCzIDTUMxpHqyCzmnHUC0EEEAgFgh0jQ4zYEYv\nAggggAACHgQIwBZ0ZsAWJdIggAAC/gQCXKtDAPY3XCgZAQQQQKCBBQjA1s4P9B6DtXmkQwABBIIV\nYBFWsF1HxRFAAAEEEMhcgBmwhTzAewuWZpEGAQQQKIRAMgMOryUE4PD6jBojgAACCJQLBHqLkABc\n3okdPWcVdEcy7EcAAQQQ6KIAAbiLcJyGAAIIIJAjgQBvFRKALeOHGbBFiTQIIIAAAjUIEIBrwCIp\nAggggEAOBfgaUg47Jc0qBXqTP00C8kIAAQQQSE+AGXB6luSEAAIIIOBDgBmwD/WMygzw5n5GMhSD\nAAIIINBFAWbAXYTjNAQQQACBnAgEeouQAGwZP6yCtiiRBgEEEPAnEOCVSgKwv+FCyQgggAACDSxA\nALZ0PjNgixJpEEAAAT8CLMLy406pCCCAAAIIhCjADNjaa4He5Lc2j3QIIIBAsALMgIPtOiqOAAII\nIIBA5gLMgC3kAa6uszSLNAgggEAhBAK9QkkALsTooxEIIIBAgwsEOFEiAFvGLKugLUqkQQABBBCo\nQYAAXAMWSRFAAAEEcijAIqwcdkpaVWIGnJYk+SCAAAIItAowA2YoIIAAAggELRDFM+Am7gEH3YfV\nKx/oKrvqjeIoAggggIAvAWbAvuQpFwEEEEAgHYFAJ0gEYEv3B3hpw9Is0iCAAAII+BMgAPuzp2QE\nEEAAgbQEApwoEYAtnc8qaIsSaRBAAAE/AnwNyY87pSKAAAIIIBCiADNgS68xA7YokQYBBBDwI8AM\n2I87pSKAAAIIIBCiADNga68Fuszd2jzSIYAAAsEKBPr+TAAOdsRRcQQQQACBkAUIwJbeC3B5u6VZ\npEEAAQQKIxDg+zQBuDCjj4YggAACDSrAIqwCdzyroAvcuTQNAQQQ8CPADNiPO6UigAACCKQlwAw4\nLckc5sMMOIedQpUQQACBsAWYAYfdf9QeAQQQQICvIRV8DATawQXvFZqHAAIIBCvADDjYrqPiCCCA\nAAIlAb6GVKIo1pMAO7ZYHUBrEEAAgeIJMAMuXp/SIgQQQKCxBFgFXeD+ZhV0gTuXpiGAQPACBODg\nu5AGIIAAAgggkJkAl6At1MyALUqkQQABBPwIBPotFQKwn+FCqQgggAACDS5AALYOgEA/YVmbRzoE\nEEAgaIEAv61CAA56xFF5BBBAAIFQBQjAlp4L8JOVpVmkQQABBAohwCroQnRjh41omXWttPzfezo8\nzgEEEEAAgewFWmZeIfL2SpEAJ0rMgC3jJe7YpgM+IfLWW5bUpEEAAQQQyEgg+s8rIs3NGZWWbjEE\n4Fo8N4fZybU0kbQIIIBAUAKbNokEukiWAGwZaXppQztYO5oNAQQQQCA/Aps356cuNdaEAGwF0yAc\ncEdbm0k6BBBAICgBNzHiHnBQ3VZbZeMZcOQ6urYzSY0AAgggUC+BgCdGzICtg0I/XRGArVqkQwAB\nBLIR0PdlvoaUjbWXUtw94EBX2nkxo1AEEEAgC4FN3APOgtl/GQF3tH88aoAAAgikLxDpJWhmwOnD\n5iZHd3M/4HsNubGkIggggECaAi4Ap5lnRnlxD7gWaO4B16JFWgQQQKD+AgG/LxOALcNDZ8AtLSzC\nsliRBgEEEMhSwF2ZdFcqsyy7m2URgGsBdB1dyzmkRQABBBCon8DGjfXLu845E4BrAQ74UkctzSQt\nAgggEIJApN9M2W67rVcomQGH0GXdqCMz4G7gcSoCCCCQsoBOinr3TjnT7LJjBmyxdp+s+BqSRYs0\nCCCAQDYCLgDzNaRsvH2Wwq+i9KlP2QgggEA7AQ3A22+/9RJ0u0MhvMzdDLg5vqa/evXqfNm5GfDi\nF6T5E5+SSFdEsyGAAAIIeBOIVq2SLftNkqZevbzVobsF56Lmm+JPMZdcconcdttt8tprr8W/1CSS\nvn37ypgxY+S8886TadOmdbed6Zw/eJA0jRwhWz40XpqOOkJEV9/17CmiAbpH/FlGL4O4SyG6zwVu\nt09rkeyPH+Ok0kO/3hQ/0UfdXHoN8Jpv+Xl63JXh8tZHV6Yec8c1re4vf61pXb76qMfc+Zre5anP\nddPzdXN56PHyTY8PGSKyZs3Wve5Diat3edr25bQ/Fr9uistJPtjE+erzSpuOizbH4tfJPldmax2j\nLVukSffppvV0bdi65/19rl6aRp+3np84aVr32j26dM7P5aeP7piWpccHDhBZu+79FC4Pl1Yf3T59\n1PP1PH3u8nBpXf1det2vm56jx9ymr91WnrZ8vzteaV95uXq+y0PT6k95Pcpf63NNq49uc0Zunzse\nP0bxjKVpw4atKd1+l4f2m57b/jyXrz7qMf0pr4/uj/u99P9GX2veuq/cSPe7vPW5bppOf7Tc9lv7\ntHrc7YvP0THbpOfqpo/xMTdO3TF9rfuTMa71KU+v5+lCIm237m/fdj1eXq/yPtJj7rWm0fPdaz2m\n5bqy9LXLW9O0d2mtoyYruZaXmxxo/UfT6ubK3Pqqbf4uTXm6Nvta89B97euo+Wod3TF97N9PZN36\n98tQryE7bv3RAOzSuroE8piLADx9+nRZvny5PPDAAzJ27Fjp16+frF27VhYtWiQzZsyQDfF/1jPP\nPLNT0tmzZ8tvf/vbiumefvppGTFiRMVjne7UAaIdHN8Dbpp0gPS44Nz4DXZgfOkjvvmv40iP6aaP\n7rk7x712x9x+fXQDXI+5/xRbc9r6r0tbLY/yY64O7tGdr6/Ly3PH9bHSf9jy/wDt/5O5c+Nyo149\npam59Q1O6+F+3H8el7ZSObrPbXpeXL+m1ke3u/Qf3B3XA64MbY87R/eVbfGRrZumiY8lb4j6XLey\nfcnr1ry3HorfUNvllZSnB8vSlZ5X2pdkGv/TM/6v1Rz/ijxXbnm+rXVwSZPHSnm5fe5cd57LU090\nx5y5O8dl7o63P8e9bp/eneceKx13+9yjptXnbix1kndTj/jNM2oNspXSunz1UTeXZuurUh+WXsbl\nug9vSR30QHke7ny3Lz7sgqbr7/LXyXM9p7z8snNLeZflU9oXp3N5JY/6WuvTWgd3THe1yd+Vp+Wo\no3ut6VrP1aelc1x/677y4+37wB0vr3/5PldW+f/z9nm4/F0e7rXmo5u+dvv00bVBH93xrc+2ptP9\nGkBdfu6xPK3mU1aPKH5dctQnW1okWvqyyLJ/i7y90uUe3GMuAvC8efPk8ccfl1122aUEOGjQIJk0\naZJcc801cvHFF5sC8Be+8AU57LDDSnmUP7nnnnuSQF6+r6bnboAOH5YE4ZrOLWji5D9EAG2rVM9K\n+7QpHe0PoJkNXcWu9lv789q/bmjUHDW+Ur80feTDsmXefJGVq94P5jmqs6UquQjAEyZMkAULFsiJ\nJ564TZ3vv/9+GTZs2Db7K+3QmbP+VNqGDh3avXvLLgDvtFOl7NmHAAIIIJCxQFP8vh41vfj+lYGM\ny+9ucbkIwJdddplMmTJFZs2aJbvttlt8dXdgfGtxjTz//PPx7ZFmmTt3bnfb2f3zWwNw07Ch3c+L\nHBBAAAEEui+g78fusrc+BrblIgBPnDhR9B6tXoZetmxZcj9YZ7163/eggw6KfT3Dug7W1dkjPhBY\nF1NdBBBAoKACo0aJDB4ssj5eoBXglosArG59+vSRQw45JL+E8Qy454J50jR6dH7rSM0QQACBBhLo\nEX8bpSX+Fk30i/8NstVl32EIsv7ZVNrNgH3PxLNpLaUggAAC4QjoamndAnx/JgBv7brO/9V7wK6j\nO09NCgQQQACBLATcBCmLslIugwBsAdUOJgBbpEiDAAIIZCugE6Py7xJnW3q3SiMAW/laV0Fbk5MO\nAQQQQCADgQAvPTsVArCT6OxRf6kLl6A7U+I4AgggkK1A8it9W38VZ7Yld7s0ArCVMJkBWxOTDgEE\nEEAgEwFmwJkw+yvE3eRnBuyvDygZAQQQqCSg78v6R20C3JgBWzuNe8BWKdIhgAAC2Qm4GbB7zK7k\nbpdEALYQaseyCtoiRRoEEEAgW4GAr0wSgI1DJflzZwF+wjI2j2QIIIBAmALJJWgWYYXZeZZa66+i\n1hmw/g1LNgQQQACB3Ah4/1sB3ZBgBmzFC/SL3tbmkQ4BBBAIUqA0Aw6v9gRga59xD9gqRToEEEAg\nO4GAbw0SgK3DhBmwVYp0CCCAQHYCOgPW9+cANwKwpdNYBW1RIg0CCCCQvYC+PxOAs3fPtERmwJly\nUxgCCCBgEnBfQwrwUjQzYEsPMwO2KJEGAQQQyF4gfn8O9WuiBGDrcGEGbJUiHQIIIJCdgJsBZ1di\naiURgC2UOgNubuavIVmsSIMAAghkKaDvz4FOkAjA1oGiHRzwJy1rM0mHAAIIBCXAKuiguqtrleWP\nMXTNjbMQQACBegq4GXA9y6hT3syArbDMgK1SpEMAAQSyEwj4yiQB2DJM9BMWM2CLFGkQQACBbAXc\nDFgfA9sIwNYOYwZslSIdAgggkJ2AuwdMAM7OPNOSmAFnyk1hCCCAgFnAzYDNJ+QnITNga18wA7ZK\nkQ4BBBDITkD/TKy+PzMDzs4805KSGTBfQ8rUnMIQQAABi0CAgdc1ixmwk+jsMV6E1RTwarvOmsdx\nBBBAIEiBHq2LZAOsPAHY2ml6iYMNAQQQQCBfAsyA89UfqdeG+7+pk5IhAgggkIoAq6BTYcxvJvoJ\nK+BPWfmFpWYIIIBANwUCfm/mErS17/U+AxsCCCCAQL4EmAHnqz/qUpuAP2XVxYNMEUAAgTwI6Htz\noL+pkBmwdQCxAtoqRToEEEAgO4GA35sJwNZhEnAnW5tIOgQQQCA4AZ0BB/otFQKwdbRxCdoqRToE\nEEAgOwGdHOkl6AA3ArC105gBW6VIhwACCGQn4GbAAU6SCMCWYRLo7xm1NI00CCCAQNACzICD7r7O\nK6+frJgBd+5ECgQQQCBrAX1/ZhV01uoZlscMOENsikIAAQRqEHAzYC5B14AWUlLt2AA7NyRi6ooA\nAgh0SUDfm1kF3SW6cE4iAIfTV9QUAQQaR0BnwFtYBV3sDtc/+syGAAIIIJAvAQ3AzIDz1Sep14YZ\ncOqkZIgAAgh0W0Dfm1mE1W3G/Gagn65YBZ3f/qFmCCDQuALMgBug75kBN0An00QEEAhOgBlwcF1W\nW4W1g5kB12ZGagQQQCALAWbAWSh7LIPvAXvEp2gEEECgioBOkLZsCfKrovwqyir9WjrEDLhEwRME\nEEAgVwLMgHPVHXWpTJMGYTYEEEAAgXwJ6Hszfw0pX32Sem0IwKmTkiECCCDQbQE3Aw7wPZpL0Nbe\nZxGWVYp0CCCAQGYCydVJfhFHZt7ZF8QirOzNKREBBBCwCujslxmwVSvAdMyAA+w0qowAAg0hoME3\nwFkwl6CtozPAT1fWppEOAQQQCFpAJ0gBvkcTgK2jrgeroK1UpEMAAQQyF2AGnDl5NgXqJ6sAP11l\ng0MpCCCAgGeBQG8RMgO2jBv9ZBVoB1uaRxoEEEAgaAEuQQfdfdUrr7NfAnB1I44igAACvgT0/ZlL\n0L7061wuM+A6A5M9Aggg0A2BQG8Tcgm6G33OqQgggAACORBgBpyDTqhXFTZvFlm7rl65ky8CCCCA\nQHcEVq7sztnezu3lreSACu5xxP+I6A8bAggggEDuBHq981bu6mSpEJegLUqkQQABBBBAIGUBAnDK\noGSHAAIIIICARYAAbFEiDQIIIIAAAikLEIBTBiU7BBBAAAEELAIEYIsSaRBAAAEEEEhZgACcMijZ\nIYAAAgggYBEgAFuUSIMAAggggEDKAgTglEHJDgEEEEAAAYsAAdiiRBoEEEAAAQRSFiAApwxKdggg\ngAACCFgECMAWJdIggAACCCCQskBTFG8p55nL7J555hk5+uijZeLEiTXV78knn5SNGzdKz549azqv\nURLr8Hnvvfekb9++jdLkmtupPr1792YMdSDHGOoApmw3Y6gMo8JTF8YOPPDACker71q6dKnMnz9f\nRo4cWT1hPY7GFWerIjB9+vTo2WefrZKisQ+9/fbb0eTJkxsboZPWn3322dHChQs7SdW4h1esWBEd\nd9xxjQtgaPlZZ50VPffcc4aUjZnkzTffjE444YTgGs8l6Hp8qiFPBBBAAAEEOhEgAHcCxGEEEEAA\nAQTqIUAArocqeSKAAAIIINCJAAG4EyAOI4AAAgggUA8BAnA9VMkTAQQQQACBTgQIwJ0AcRgBBBBA\nAIF6CDTM94C7irdy5UoZMGBA8j3OruZR5PNaWlok/iqSDB8+vMjN7FbbdAwNHDhQtttuu27lU9ST\nGUOd9yxjqLpRqGOIAFy9XzmKAAIIIIBAXQS4BF0XVjJFAAEEEECgugABuLoPRxFAAAEEEKiLAAG4\nLqxkigACCCCAQHUBAnB1H44igAACCCBQFwECcF1YyRQBBBBAAIHqAgTg6j4cRQABBBBAoC4CBOC6\nsJIpAggggAAC1QUIwNV9OIoAAp4E4j/u6qlkis27QHNzsxRhfBCA45G2evVqif+Ys3zoQx+Svfba\nSx577LGK48+aruLJge+84oorZO+995YxY8aIPq+0bd68WS644AL52Mc+lvxcdNFFsmnTpkpJC7ev\n1rGxdu1a+eAHPyh/+tOfCmfRUYMsY0jPvfrqq2WPPfZI/i+effbZsmXLlo6yLNR+6xh69dVX5eST\nT5Z9991XjjzySPnzn/9cKIfOGvPKK68k/3eWLl3aYVLrWOswg6wOxJ8iGn47/vjjo8svvzyKf51Z\ntGDBgmjnnXeO/vvf/27jYk23zYmB77j77rujT37yk9E777wTvfHGG9E+++wTzZ07d5tW3XzzzdEx\nxxwTxUE3+fn85z8f6b5G2GodG9OmTYsGDx4czZ8/vxF4IusYeuihh6KJEydG69ati+IPdNGUKVOi\n22+/vSGMrGPoa1/7WvT9738/MXnyySejsWPHJlaNgHTLLbdEu+22WxT/WtdoyZIlFZtsHWsVT854\nJzPg+JPOH//4RznrrLOkqalJDj74YBk1apT89a9/3eYzkDXdNicGvkPb/ZWvfEUGDRoku+yyi5x4\n4oly7733btOqODDLj3/84+R3HuvvPd5zzz3lb3/72zbpirijlrFx3333SZ8+feTDH/5wESkqtsk6\nhuIPbMn/xX79+kn8IVjuuOMOOemkkyrmWbSd1jG0YsWK0u+m199Tv3z58oa4SqBX0+LgKvGHf4k/\nvHbY/dax1mEGGR5o+ACsl302btwoQ4YMKbFrkHnrrbdKr/WJNV2bkwry4j//+Y984AMfKLVGfd58\n883Sa/fk4x//uMSfTpOX7777rsyePVs+97nPucOFfaxlbOi4iq+2yA9+8IPCelRqmHUMaTq9xDpi\nxAjZaaed5Nhjj22I2xi1jKGZM2eKflA57rjj5PDDD5cbb7xRtt9++0rshdrXu3dvefDBBzv94God\na3nAafgArH9lRD9tl2877LCDrF+/vnyXWNO1OakgL9q3vW/fvqIBtqNNP6l++ctfFg3IkydP7ihZ\nYfa399GGVRpDuv+MM86QSy+9NPnrSPq6Ubb2Rh2NIf2Acueddyb3xhcvXizxZcbkTbfoTu19tL0d\njSG9qhRfKU3uk48cOVIeeeQR0UVJbFsF2lt2NNby4NXwAXjo0KGiC2LKN32tn8DLN2u68nOK8rx9\n2yv5uLZq8NVZiy6c0RlwI2ztfbTNlYzuuusueemllxKS+++/X+J76vLEE0/IsmXLkn1F/qe9USUf\nbb9eWtQFRuPHj0+upkydOlXuueeeItMkbWvvozsrGWmgvfDCC0XHkl5J0WA8b968irfMCo/WQQPb\nW1Zy7ODUzHc3fADW//D6SVMve7lN3xB33XVX9zJ5tKZrc1JBXug98X//+9+l1qjP6NGjS6/dE31z\n0JmvBt/f/OY3pftU7nhRH61jI15YJP3795d4AU3y8/rrryezvaeeeqqoNKV2WceQjiu9r+k2vexY\n7WqLSxf6o3UM6aVq/f+l6y1069Gjh+y3337y8ssvh06QWv2tYy21AruTUcaLvnJZ3GmnnRZNnz49\nWUk4Z86caNy4cckqXq2sroqOFz0k9a6WLpcNS6lSf/jDH6L4K0jRa6+9FsX/0aPdd9890tWXui1a\ntChauHBh8nzWrFlJujiwRPFloORHV7M2wlZtbJSPoXKL/fffv2FWQVvHUDzbTcaQfgtBx86BBx4Y\nXXfddeVshX1uHUOHHXZYFM+AE4f4qzjRjjvuGK1ataqwLpUaNmzYsDaroMvfh6qNtUp5+dyn9xIa\nftOgMmHChCi+7Jwscdc3TLfpV5IeeOCB5GW1dC59ER/161nuazPxAqzo4osvLjUz/p5mdOqppyav\n4++16m9OaPNz1FFHldIW+Um1sVE+hsoNGikAW8dQfBUl+TCs4yxe+Jd8DSme8ZWzFfa5dQz94x//\niPT/lX4o1h/9ak6jbe0DcPn7ULWxljenJq1Qd2bQRTpXl/fHHdtpk6zpOs0osAR6L0VXWzbCisuu\ndk2jjg2rl3UMvffee8nCovLL0dYyQk9nHUNqOXDgwNCbW7f6W8da3SpgyJgAbEAiCQIIIIAAAmkL\nNPwirLRByQ8BBBBAAAGLAAHYokQaBBBAAAEEUhYgAKcMSnYIIIAAAghYBAjAFiXSIIAAAgggkLIA\nAThlULJDAAEEEEDAIkAAtiiRBgEEEEAAgZQFCMApg5IdAggggAACFgECsEWJNAgggAACCKQsQABO\nGZTsEEAAAQQQsAgQgC1KpEEAAQQQQCBlAQJwyqBkhwACCCCAgEWAAGxRIg0CCCCAAAIpCxCAUwYl\nOwQQQAABBCwCBGCLEmkQQAABBBBIWYAAnDIo2SGAAAIIIGARIABblEiDAAIIIIBAygIE4JRByQ4B\nBBBAAAGLAAHYokQaBBBAAAEEUhbolXJ+ZIcAAnUUuO222+Qb3/hGUsLGjRsliiLp06dP8vrKK6+U\nl156SbbbbjuZOXNmHWtB1gggkIZAU/wfOEojI/JAAIFsBc4991zZtGmTXH/99aWC33nnHWlqapJB\ngwaV9vEEAQTyKcAl6Hz2C7VCoEsCP//5z+XWW29Nzv3MZz4jN998s4wZM0bGjRsnjz32mGjQHj58\nuJx44omybt26JN2KFSvk2GOPlcGDB8s+++wjf/nLXyqWrbPvXXfdVXbaaSc5/vjjZfXq1RXTsRMB\nBGwCBGCbE6kQCELgrbfekrfffjupq16OnjNnjjz66KNy9NFHy2GHHSb9+/eXZ555Rl599VWZP39+\nku60005LZsyLFy+WGTNmyLRp07Zp64YNG+Sss86S3/3ud8ll7nfffVd++tOfbpOOHQggYBcgANut\nSIlAcALf/OY3ZdSoUfKlL31Jmpub5YILLpARI0bIpz/9aXniiSdk1apVMnfuXLnwwgulb9++Mnny\n5OT4s88+26atelm7paVFHn74YdF7z/fdd59cdNFFbdLwAgEEahMgANfmRWoEghLQYKvbDjvskATW\nAQMGJK979+6dBGSdCWtwPfTQQ+UjH/lI8rNkyZLkcnWSsPWf7bffXu6++2755S9/KSNHjkxm1C+8\n8EJ5Ep4jgECNAgTgGsFIjkBIAr16Vf+ig94b1gVbCxculDfeeCP50dmvXpYu33T2+9GPflT++c9/\nJj8DBw4srcYuT8dzBBCwCxCA7VakRKBwAjoT1sVaN9xwQ3KJefny5bLnnnuK3g8u3/S+8oQJE5J7\nx+PHj5cjjzyy/DDPEUCgCwLVPx53IUNOQQCBsAT0Xq7eI/7FL34hPXv2lPPPP1/23nvvNo3QldPf\n/e535VOf+pT069dP1q9fnyzwapOIFwggUJMA3wOuiYvECBRXQL+ONHTo0OSecLVWarphw4ZVS8Ix\nBBAwCBCADUgkQQABBBBAIG0B7gGnLUp+CCCAAAIIGAQIwAYkkiCAAAIIIJC2AAE4bVHyQwABBBBA\nwCBAADYgkQQBBBBAAIG0BQjAaYuSHwIIIIAAAgYBArABiSQIIIAAAgikLUAATluU/BBAAAEEEDAI\nEIANSCRBAAEEEEAgbQECcNqi5IcAAggggIBBgABsQCIJAggggAACaQsQgNMWJT8EEEAAAQQMAgRg\nAxJJEEAAAQQQSFuAAJy2KPkhgAACCCBgECAAG5BIggACCCCAQNoCBOC0RckPAQQQQAABgwAB2IBE\nEgQQQAABBNIW+P/1x9X20c28pwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# (b) With N=512, use matrix algebra to solve equation (1) numerically with this decay kernel.\n",
    "\n",
    "N <- 512\n",
    "\n",
    "G <- function(i,j,dt){\n",
    "    # function to calculate G_ij\n",
    "    if (i == j) {\n",
    "        return(2*dt - 2*log(1 + dt))\n",
    "    } \n",
    "    else if (i > j) {\n",
    "        return(log((1+(i-j)*dt)*(1+(i-j)*dt) / ((1+(i-j-1)*dt)*(1+(i-j+1)*dt))))\n",
    "    }\n",
    "    else {\n",
    "        return(log((1+(j-i)*dt)*(1+(j-i)*dt) / ((1+(j-i-1)*dt)*(1+(j-i+1)*dt))))\n",
    "    }\n",
    "}\n",
    "\n",
    "# To solve equation (1), we need to construct mat_G and ls_const\n",
    "vec_const = rep(1/10000,N)\n",
    "mat_G = matrix(0,N,N)\n",
    "for (i in 1:N) {\n",
    "    for (j in 1:N) {\n",
    "        mat_G[i,j] = G(i,j,1/N)    \n",
    "    }\n",
    "}\n",
    "\n",
    "# Solve\n",
    "vec_v = solve(mat_G) %*% vec_const\n",
    "\n",
    "# (c) Plot your results \n",
    "plot(seq(1/N, 1, 1/N), vec_v, type = \"l\", col = \"red\", xlab = \"Time s\", ylab = \"u(s)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimal interval VWAP strategy in the square-root model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. (12 points)\n",
    "\n",
    "Consider equation (11) of Lecture 10 for the cost of liquidation using interval VWAPs:\n",
    "\n",
    "$$\n",
    "\\cC=\\sum_{j=1}^m\\,C_{jj}+\\sum_{i< j}^m\\,C_{ji}.\n",
    "$$\n",
    "\n",
    "(a) With $m=3$ (that is 3 interval VWAP slices), and assuming the square-root process where\n",
    "\n",
    "$$\n",
    "h(v)=\\frac  3 4 \\sigma\\,\\sqrt{\\frac{v}{V}} =\\frac  3 4 \\sigma\\,\\sqrt{\\frac{\\dot x}{V}}, \n",
    "\\quad \\mbox{ and } \\quad\n",
    "G(\\tau)=1/\\sqrt{\\tau},\n",
    "$$\n",
    "\n",
    "compute explicit expressions for the matrix elements $C_{ji}$ in terms of $x_1$, $x_2$, $t_1$, and $t_2$ (noting that $t_0=0$ and $t_3=T$).\n",
    "\n",
    "(b) Minimize the resulting expression for the expected cost of liquidation numerically to find optimal choices of $x_1$ and $x_2$ when $t_1=T/3$, and $t_2=2 \\,T/3$.\n",
    "\n",
    "(c) What is the percentage saving from trading with the 3-slice strategy relative to a VWAP?\n",
    "\n",
    "(d) Explain to your non-mathematical manager why this strategy has lower expected cost."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>\n",
    "\n",
    "#### Solution\n",
    "**(a)** According to the lecture notes,\n",
    "\n",
    "\\begin{equation}\n",
    "C_{ii}=\\frac{x_i}{t_{i}-t_{i-1}}\\,h(v_i)\\,\\int_{t_{i-1}}^{t_i}\\,dt\\,\\int_{t_{i-1}}^{t}\\,\\frac{1}{\\sqrt{t-s}}\\,ds=\n",
    "\\frac{4}{3}\\,x_i\\,h(v_i)\\,\\sqrt{t_i-t_{i-1}} = \\frac{\\sigma}{\\sqrt{V}} x_i^{3/2}\n",
    "\\end{equation}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "C_{ji}&=&\\frac{x_j}{t_{j}-t_{j-1}}\\,h(v_i)\\,\\int_{t_{j-1}}^{t_j}\\,dt\\,\\int_{t_{i-1}}^{t_i}\\,\\frac{1}{\\sqrt{t-s}}\\,ds\\\\&=&\n",
    "\\frac{4}{3}\\,x_j\\,h(v_i)\\,\\frac{\\left(t_{j}-t_{i-1}\\right)^{3/2}\n",
    "-\\left(t_{j}-t_{i}\\right)^{3/2}\n",
    "-\\left(t_{j-1}-t_{i-1}\\right)^{3/2}\n",
    "+\\left(t_{j-1}-t_{i}\\right)^{3/2}\n",
    "}{t_{j}-t_{j-1}}\\\\&=&\n",
    "\\frac{\\sigma}{\\sqrt{V}}x_j\\sqrt{x_i}\\,\\frac{\\left(t_{j}-t_{i-1}\\right)^{3/2} - \\left(t_{j}-t_{i}\\right)^{3/2} - \\left(t_{j-1}-t_{i-1}\\right)^{3/2} + \\left(t_{j-1}-t_{i}\\right)^{3/2} }{(t_{j}-t_{j-1})\\sqrt{t_{i}-t_{i-1}}}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "**(b)** When $t_1=T/3$, and $t_2=2 \\,T/3$:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "C_{11} &=& \\frac{\\sigma}{\\sqrt{V}} x_1^{3/2} \\\\\n",
    "C_{21} &=& (2\\sqrt{2}-2)\\frac{\\sigma}{\\sqrt{V}}x_2\\sqrt{x_1}\\\\\n",
    "C_{22} &=& \\frac{\\sigma}{\\sqrt{V}} x_2^{3/2} \\\\\n",
    "C_{31} &=& (3\\sqrt{3}-4\\sqrt{2}+1)\\frac{\\sigma}{\\sqrt{V}}x_3\\sqrt{x_1}\\\\\n",
    "C_{32} &=& (2\\sqrt{2}-2)\\frac{\\sigma}{\\sqrt{V}}x_3\\sqrt{x_2}\\\\\n",
    "C_{33} &=& \\frac{\\sigma}{\\sqrt{V}} x_3^{3/2} \\\\\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "$par\n",
       "[1] 6.471808e+03 5.276104e-09\n",
       "\n",
       "$value\n",
       "[1] 883282.3\n",
       "\n",
       "$counts\n",
       "function gradient \n",
       "     249       NA \n",
       "\n",
       "$convergence\n",
       "[1] 0\n",
       "\n",
       "$message\n",
       "NULL\n",
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "cost <- function (x){\n",
    "    # @summary: cost function\n",
    "    # @param c11, c22, c33: cost parameters for i=j\n",
    "    # @param c21, c32, c31: cost parameters for i<j\n",
    "    c11 = x[1]^1.5\n",
    "    c22 = x[2]^1.5\n",
    "    c33 = (10000 - x[1] - x[2])^1.5\n",
    "    c21 = (2*sqrt(2)-2) * x[2] * sqrt(x[1])\n",
    "    c32 = (2*sqrt(2)-2) * (10000 - x[1] - x[2]) * sqrt(x[2])\n",
    "    c31 = (3*sqrt(3)-4*sqrt(2)+1) * (10000 - x[1] - x[2]) * sqrt(x[1])\n",
    "    return (c11 + c22 + c33 + c21 + c32 + c31)\n",
    "}\n",
    "\n",
    "# Optimization\n",
    "optim(c(2000,2000), cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>Therefore, the distrubtion of total shares at three layers should be 64.72%, 0% and 35.28% respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>\n",
    "**(c)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 11.67179\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# What is the percentage saving from trading with the 3-slice strategy relative to a VWAP?\n",
    "\n",
    "cost_vwap = 10000*sqrt(10000)\n",
    "cost_3layer = cost(c(6472,0))\n",
    "print((cost_vwap-cost_3layer)/cost_vwap*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkblue'>\n",
    "\n",
    "The percentage saving from trading with the 3-slice strategy relative to a VWAP is 11.67%.\n",
    "\n",
    "**(d)** This strategy has lower expected cost is mainly because there exists a temporary component in the market impact. \n",
    "\n",
    "Tranditionally, we sell all shares continuously then the market impact will be accumelated when time goes by. But in this strategy, there is a cool-down / wait for the market impact between selling different layers, then the expected cost in slippage will lower."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}