{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Bayesian Test for Cointegration\n",
    "\n",
    "## Python notebook implementation\n",
    "\n",
    "The code in this notebook is adapted from the original MATLAB implementation by Chris Bracegirdle for the paper [*Bayesian Conditional Cointegration*](http://icml.cc/2012/papers/570.pdf) presented at [*ICML 2012*](http://icml.cc/2012/).\n",
    "\n",
    "<a class=\"typeform-share button\" href=\"https://bayesiancointegration.typeform.com/to/gRS0pq\" data-mode=\"popup\" style=\"display:inline-block;text-decoration:none;background-color:#267DDD;color:white;cursor:pointer;font-family:Helvetica,Arial,sans-serif;font-size:20px;line-height:50px;text-align:center;margin:0;height:50px;padding:0px 33px;border-radius:25px;max-width:100%;white-space:nowrap;overflow:hidden;text-overflow:ellipsis;font-weight:bold;-webkit-font-smoothing:antialiased;-moz-osx-font-smoothing:grayscale;\" target=\"_blank\">Contact me</a><script>(function(){var qs,js,q,s,d=document,gi=d.getElementById,ce=d.createElement,gt=d.getElementsByTagName,id=\"typef_orm_share\",b=\"https://s3-eu-west-1.amazonaws.com/share.typeform.com/\";if(!gi.call(d,id)){js=ce.call(d,\"script\");js.id=id;js.src=b+\"share.js\";q=gt.call(d,\"script\")[0];q.parentNode.insertBefore(js,q)}})()</script>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "' Bayesian Cointegration \\nImplementation of a Bayesian test for cointegration\\n\\nWritten by Chris Bracegirdle\\n(c) Chris Bracegirdle 2015. All rights reserved.'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from random import random, normalvariate\n",
    "from numpy import sum, empty, polyfit, inf, array, isinf, where, nan, isnan, isreal\n",
    "from scipy import log, sqrt, exp, std, pi\n",
    "from scipy.misc import logsumexp\n",
    "from scipy.stats import norm \n",
    "\n",
    "\"\"\" Bayesian Cointegration \n",
    "Implementation of a Bayesian test for cointegration\n",
    "\n",
    "Written by Chris Bracegirdle\n",
    "(c) Chris Bracegirdle 2015. All rights reserved.\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some notebook-specific requirements here: we're going to show some plots, so let's show them inline. And for debugging, never swallow an overflow with just a warning, please numpy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from numpy import seterr\n",
    "seterr(over='raise')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a rather contrived function to randomly generate two time series, `x` and `y`, after making a random decision as to whether they are cointegrated, and generating according to the corresponding generating function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def GenerateData(T):\n",
    "    cointegrated = random() > 0.5\n",
    "    phi = random() * 2 - 1 if cointegrated else 1\n",
    "\n",
    "    std_eta = exp(normalvariate(0,1))\n",
    "    std_x = exp(normalvariate(0,1))\n",
    "    intercept = normalvariate(0,5)\n",
    "    slope = normalvariate(1,5)\n",
    "    #print \"Intercept\", intercept\n",
    "    #print \"Slope\", slope\n",
    "    \n",
    "    epsilon = empty([T])\n",
    "    x = empty([T])\n",
    "    y = empty([T])\n",
    "    epsilon[0] = normalvariate(0,std_eta)\n",
    "    x[0] = normalvariate(0,std_x)\n",
    "    y[0] = intercept + slope * x[0] + epsilon[0]\n",
    "    for t in range(1,T):\n",
    "        epsilon[t] = phi * epsilon[t-1] + normalvariate(0, std_eta)\n",
    "        x[t] = x[t-1] + normalvariate(0 ,std_x)\n",
    "        y[t] = intercept + slope * x[t] + epsilon[t]\n",
    "    return cointegrated,x,y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what some randomly-generated data look like. Give it a whirl!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cointegrated? False\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEACAYAAABGYoqtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VNX9+PH3yT7ZJvtKAggJBEE2RUWoERV3EFfQilYt\ntlpt/Vm1Vr8VqK3aaqu1dS0KWJWqRUFFFIsRULCy7yRASEL2PTPZJsv5/XEyNxOSIdskk8B5PU8e\nZ+69c+9nYrife3YhpUTTNE3TOuLh7gA0TdO0gUsnCU3TNM0pnSQ0TdM0p3SS0DRN05zSSULTNE1z\nSicJTdM0zaleJwkhxBAhxHohxD4hxB4hxAMt20OFEF8KIQ4JIb4QQpgdPvOYECJDCHFACDGztzFo\nmqZpfUP0dpyEECIGiJFS7hRCBALbgNnAT4BSKeWfhBCPAqFSyt8IIcYA7wDnAEOAr4AkqQdsaJqm\nDTi9LklIKQuklDtbXluBA6ib/2xgWcthy4BrW17PAlZIKRullMeADGBKb+PQNE3TXM+lbRJCiGHA\nBGALEC2lLASVSIColsPigRyHj+W2bNM0TdMGGJcliZaqpg+BX7aUKE6sPtLVSZqmaYOMlytOIoTw\nQiWIt6WUq1o2FwohoqWUhS3tFkUt23OBBIePD2nZ1tF5dWLRNE3rASmlcMV5XFWSeBPYL6V80WHb\nauCOlte3A6scts8VQvgIIYYDI4H/OTuxlFL/SMmTTz7p9hgGyo/+Xejfhf5dnPzHlXpdkhBCXADc\nCuwRQuxAVSv9FngWeF8IcSeQBdwEIKXcL4R4H9gPNAD3Sld/K03TNM0lep0kpJTfAp5Odl/i5DNP\nA0/39tqapmla39IjrgeJ1NRUd4cwYOjfRSv9u2ilfxd9o9eD6fqSEELXRGmapnWTEAI5wBquNU3T\ntFOQThKapmmaUzpJaJqmaU7pJKFpmqY5pZOEpmma5pROEpqmaZpTOklomqZpTukkoWmapjmlk4Sm\naZrmlE4SmqZpmlM6SWiapmlO6SShaZqmOaWThKZpmuaUThKapmmaUzpJaJqmaU7pJNHHTsX1MJpl\ns7tD0DStn+gk0YcamxvxWOxBSU2Ju0NxiU3Zm5j+1nQ8F3vS1Nzk7nA0TesHLkkSQoglQohCIcRu\nh22hQogvhRCHhBBfCCHMDvseE0JkCCEOCCFmuiKGgej7498DUGgtdHMkrjH9relsyt4EQL41383R\naJrWH1xVkngLuOyEbb8BvpJSjgLWA48BCCHGADcBKcAVwMtCCJcsszfQfJbxGcApU5JwlFOZ027b\n15lfc7zquBui0TStr7gkSUgpNwHlJ2yeDSxreb0MuLbl9SxghZSyUUp5DMgAprgijoFmyY4ljAof\nRXFNsbHt++Pfc6jkkBuj6h5bkw2Ac944B4DVc1czLXEah8sOtzluY9ZGZiyfwewVs/s9Rk3T+k5f\ntklESSkLAaSUBUBUy/Z4wPExNLdl2ymlvrGe8tpyfjT0RxRVFxnbF29YzL/3/duNkXXd4/99HN+n\nfJFSsjVvKwDXjLqG2aNm87/c/7U5dvmu5QBsz99OfWN9v8eqaVrf6M+G61Ovm89JFFUXERUQxciw\nkcZTd1NzE5uyN5Fblevm6DpXUVfBc5ufA6C0trTNvlHhozhacbTNtjBTGAAB3gE8v/n5/glS07Q+\n59WH5y4UQkRLKQuFEDGA/XE6F0hwOG5Iy7YOLVy40HidmppKamqq6yPtAwXWAqIDoxkdMZplu5bh\nsciDyXGTqaqvIs+a5+7wOlVoLSQhOIFAn0CyK7MBuHHMjQBEB0ZTaC1ESskNH9zAhzd+CEByeDIX\nDbuIFXtX8Nvpv3Vb7Jp2uklLSyMtLa1Pzu3KJCFafuxWA3cAzwK3A6sctr8jhPgrqpppJNC27sKB\nY5IYTAqsBcQExpASkcL+4v0AbM3byrTEaaSXprs5us4V1xQTGRBJhH8EOZU5jAwbyVMzngIgOiCa\n7fnb+ST9E1YeWInHYg/umngXD53/EBNiJhhVU5qm9Y8TH6AXLVrksnO7qgvsu8B3QLIQIlsI8RPg\nGeBSIcQh4OKW90gp9wPvA/uBNcC98hQccVZYXUhMQAzDQ4e32T73zLnkWfIory0fsIPSmmUzz29+\nnkj/SBKDE8mpyqGmoYYA7wAAogKikMg2jdQ7CnaQEJxAuCmcstoyd4WuaZqLuap30y1Syjgppa+U\nMlFK+ZaUslxKeYmUcpSUcqaUssLh+KellCOllClSyi9dEcNAY69u8vJQhbXx0eMBVVWTaE4k9vlY\n7l59tztDNHx19CvEIsFHBz4CoLSmlI8Pfky4fzgJ5gSyK7OptlXj7+0PgMnb1O4c2/O3M33odMJM\nYW2ShL13lKZpg5Mecd0HpJT839f/R2xgrLGtobkBgHBTOFEBUdQ31bMha4O7QjQ0NDVw6duXAvDC\n9y8AUNdYB4C3hzcJwQmtJQmfAONzUxOmAq3J757J9xDoE4jZz4zVZqWxuREpJb5P+WK1WfvzK2ma\n5kI6SfSB7Mps/L39uefsewD4+vav+eDGDwDw9/Y3prQI9Al0W4x2OVWtvZHttX41DTUAeApPEs2J\n7C/eT7NsxtvD2zj22zu/BSDIN6jNZz2EByF+IZTXllPfpLrCVtQZhUhN0wYZnSRcZPWh1UYD9Y6C\nHaQOS8XH0weA1GGpjIkcA0BsUCw+nj54CA+qG6rdFq9dobWQIB91o/fz8gMw4rpvyn0kmBPYXbib\nJtlERwPj7YmjrqnO2BbuH05RdZGRHBZ8sqBPv4OmaX1HJwkXmb1iNou/WQyo+vmJMRPbHSOflCSa\nE/lk3ifkPJhDVX1Vf4fZzu7C3UYCsz/51zTUMDVhKmOjxjIkeAiTYiex6SebOvy8vc3Fse9BmCmM\ni5ZdRPJLyQB8fvjzvvwKmqb1IZ0kXMB+g7S3QWRWZJIcnuz0eJO3iVC/ULckCSkl2/O3G+9/9tnP\n8PXyZetPt2KptwBQbas2ejJ5eXixbcE2Lki8oMPzeXl4seWuLTw38zljW7gpnOKaYiw2Sx9+E03T\n+oNOEi5gb5gtr1PTVxVVFxEdEH3Sz/h5+dHY3NjvvX++y/mOya9PpqFJNaQnhyfz2tWvEeQbZNzU\naxpqjJ5MnfHy8OLcIecSFRBlbLOPvtY0bfDTScIF7Mlh2a5lrD602piS42SEEAT7BhtP7/1l8/HN\nAPg8pdpL7LEG+QRxuOwwn2d83q4n08nYq5schZvC27wP9QvVXWE1bZDSScIFymvLGRE6AlBtEzsL\ndnaaJACCfIL6vcopqyKLq5KuAtQkhFablRC/ECIDIvH39ueT9E/YW7SXxODELp0vPqj93IxhpjCC\nfYMBNfCuvK6ch7982HVfQtO0fqOThAuU15UTHxzPb6e1zlcUGRDZ6eeCfYP7PUlkV2Vz18S7CPQJ\nJKsyi1C/UDyEB14eXqyau4oDJQdYe2Qts0d3PuV3zoM5/Hnmn9ttD/cPJz4onlmjZhm/k7I6PQpb\n0wajvpzg75RXWlPK1ryt1DTUEOoX2mbcg73768m4I0kUVRcRHRiN1WZl1N9HcUboGca+5PBk0o6l\nATAhZkKn5xoSPKTD7WGmMOKC4lg1dxVSSn71xa+IMEW4JH5N0/rXoE8Sf/v+b+wv3s+rV7/a79d+\nY/sbPPbfxwC4Y8IdXa7Ht3NHknDsuQRq8JtdQnACXh5eNDY3GmMmeuKKkVcwKnwUoNpeRoaNJNw/\nvJNPaZo2EA366qbnNz/Pa9tec8u19xTtMV57Ck8amxu79Xm3JImGagJ8ApgcOxlQpSE7IQRvXPMG\nd028q1fXMPuZmRjbOk5k/lnzjak+NE0bXAZ9knDnzee7nO+M12W1ZVTb1Ejlw/cfdvaRNvo7Sfzp\n2z9xtPwoAd4BbF2gpvO298yyu2PCHfxz1j9del2Tt4nahlqXnlPTtP6hk0QPr5lTmdOm+6rVZjWq\nm0aEjejSeYJ9gympKcHvKT8+Tf+0T2J19OhXjwIYcYb4hfT5NUGNCdElCU0bnAZ8knhxy4s8su4R\np/vdcfMJ/GMg096axvkJ5xPpH8n5Q87njWve4P4p93Psl8e6fJ5g32D+ve/f1DfVs+7Iur4L+AT2\nNomRYSP75XomLxO1jbokoWmD0YBPEn/Z8hf+/F37bpagqnhsTTbMvmaKq4v7LaYm2UR2ZTZTh0wl\n/6F8vr3zW4aGDMXb05uhIUO7fJ5g32D2FO1h5oiZZJRl9GHEGCOsATw9PAH47JbPSP9F36+SZ/I2\nUdNQw7z/zNOD6jRtkBnwScI+Q+mKvSva7fs682tiAmOorK8k6rnOB6+52tSEqXh6eHY4O2pX2L/b\n/VPu7/MkUVzTPolGBUSRFJ7Up9cFVZLYkLWBFXtXcLDkYIfH/Gf/fwbsSn2adjob8EliX/E+AOb9\nZ167fUfLj3J23NnG+/5YBdWxB5PjtXvCPj/SRcMuIqcyp83TvqsVWAsYHz2efffu67NrOBPgE0C+\nNR+A749/325/bUMtN3xwA4dKDvV3aJqmdWLAJ4mTsdqshPqFAqoLqn2q6750oPiA8bq74yJOZK+n\nD/AJID44nqPlR9vsF4sExyqO9eoadgXWAmICY4xpwfuTvcQEsODTBe3WwLaXLlz1XTVNcx23JQkh\nxOVCiINCiHQhxKM9OYfFZmFs1FiWzFpCk2zC9If2ay+7mr1k4wqOI5bHR49nZ8FO431WRRYAewr3\n8HXm172+VqG1kJjAmF6fpyfs8ziZfc0ApJe2bQexjzdxR5LYcnwLuwp29ft1NW2wcEuSEEJ4AH8H\nLgPOBOYJIUZ35bM5lTlGtcTWvK2E+IVw58Q7jf0HSw6y5fgW1wfd4tP0T7l9/O38d/5/e32umSNm\n0vw7VQ+fEpHSpl0iq1Ilifkfz2fG8hm9vtauwl0nXeOiL9mTxIPnPcj1Kde3SxJ7i/YCuGUt7POX\nnG90DdY0rT13lSSmABlSyiwpZQOwAuhwRrmFFy5s8370P0Yz/a3pbMrexMbsje3WiZ69YjbnLzmf\nrXlbuX/N/U4DOFp+lPrG7lVP2ZpsvLPnHQK8A5gxvPc3bsBo9A73D29TDWMfCW1fAjSnMqf9h7th\ny/EtTE+c3qtz9JR9HewQvxDGRY1rlySKqosYah5KXqmVjRthzRpocNI8U18PxcVQVgZ1deq4Q4fU\n++6ytwHpBnNNc85dczfFA453veOoxNHOk6lPsvCbhcb7moYaahpqmP6WuuHZmmxkdNAx6KMDH/H3\nH/7OS1e+1Gb70fKjZJRmcPk7l/O7H/2ORRctchpkejrk5IDJBJ+XvMKYeFU99Oupv6a6Wt2gQlw0\nHi3UL9SodimvLSezItPYNyR4CIkvJPL2nLe5fOTlhJvC2/WosrfZNzSAtzc0N8PRo1BZqbYVVhcR\nGxTrklilhIoKCA3t2vH2NokI/wh8GqP4cM/H7ItWsW3dCjssVkqKY3g1zcp/j0JNDQgBqalgs8He\nvZCZCeecA3v2qGvX1ak4hIDAQPU6MhIWL4bwcCgthYAAqKqC8ROaSUmR2Oo9CWptHuF41XEAiqtL\nWbwYDh6Eyy+HvDzw94foaPj6a3WNOXNg1Ch47z1IS4NbblHxfPml+g5eXnDWWTBjBowdC56evfsd\nV1bC55+D1QoXXACRQ0vYcnwLE/yvZtUq9bc5bBiUl8PFF0N8PMTEwLFj6v99Xh7k5kJ+vkqsYWHq\n9xkbq94XFcHZZ6vvpmknM+An+PP0XAgtD8CmZ9LgN233/7+HmvDZB9POzmJL5AKOVh2GYHh+yTEY\nBTt2NhMX60F5OYRG1fLQ2sf4OON9AJavyqFsJRw4oJ5Es7PVjSUuTt10MjMhORnqbI1svfJePNaN\nwaN2OkmRw2lqAh8fdVMYNQrOOANGjoTdu9XN+Te/UdscWa3q/B39wwwzhVFUXQRAyj9SKKwuNPad\nE3cOx6uOk3Ysjds+uo2Nt29mdOB5+PhAUBCsXAkPPaS+Q309NDWpn5gY9V3q6uDYtcU8+etIzjlL\n3XSio9WNTUo4flwlkm3bYOdOdWO88EJ1UxFCnXfrVvW72LdP/Rw4oG5OY8aoG2J1tbpOaSlYLOqm\ndcMNcO654O/vickjiEXzLyejOBOv69K54mkVu58fFF1sJSk2hnMutvL6NSqWL79UCbqpCa6/Xv0u\njx6FIUNg4kQVV2OjulZIiHq9cSP8/vfq/4vZDBkZ6qa44OvZNHpV4Ll8Iw2/9eDsPRu4cuw0YiZU\n4i9C2XO0iJGZMH06fPqp+kxZmUpGM2ao7/bgg5CVpb7z/Pnwr3/BU09BUpL6nLe3+n//2GMwYoSK\n8aKLVBxRUepvxGxWN/8dO1r/X1laBu3v2AGXXAIffKCuk5UF55+v/l4efhj8r3qRvKSnCH1RcvHF\nKiHu3q0eYObMUYmhvl4lSG9vFUNCgvr/39QES5fCW29BYaF67+enkuukSXDPPSrhHTqk/hYCA1UC\n8vJSMQcH9/zfr5QqthOTZnOz2ldcrGK22dR33bJFXT8/X/3eRo8GX9+eX/90kZaWRlpaWp+cW/RH\nt9F2FxXiPGChlPLylve/AaSU8tkTjpO5uZL4N9RdtfiBeiL/pv5i/jDtRa4ZfQW+tcM5Y5gXXl5w\n3Xtz+fzIJ9Q11RBrGkp+bRYxb1qpLA4gPh4OJ98HU142zj/WexY3Nq5i1ChITFT/sKqq1I0yLAym\nTlX/mAqsBcQ+r57Cbxp9K388+19ERak//FWroKAAdu1S/8jGjVPbly9Xf+BXX61uCFu2qBttfDzc\neKO62U2frv6hnHUWrDuyntn/uZg7qg6zNFiNhD7X9hjn8HPSQ/7BlzXPEl5+GaWhX8CSTQRXXkBT\nE9R5FTBqSDT/+Ltg7Fj1BOzpqW4UAB4ealS6+WkzfzbXkZYm2LdPPUnW1KgbQFmZuimcf76Ka/58\ndbM6dEg9jY4dC8OHq6fT6GgV97hx8Oab6qbi7a22BwSo0kV9vXoq37hR/Xh7w7RpcPvtMO3iSsa+\nFY/lMYtRGpr25jTGRo3FYrPwznXvuPzvLfLPkZTUlPDOnBXc+tFcAB6yNrFm77fkT/o5Nr9cqv+v\nvJOzdI2UqqTxww8qkRYXq9/b4cMqIQQFweTJKoE3Nanfm82m/vZWrlQllGnTWhM8qIeLX7zzF5YV\nPETFryRmc/trdrdE0NCgSmhpabBsmfr7TUlRiaWoCPbvVzfvqiqVKEJC1P/3adNUovPwaE0imZmw\nebMq5Qmh4jGZ1L5Dh9SDg7+/SkRms/p7+OYb9XcXGKgeLkym1usFBqokmJOj3p9/vkqEYWHqQUyX\ngDonhEBK6ZLfkrtKEj8AI4UQQ4F8YC7QfiAE6h9K8++aiXouimafCmP7YzPub1flEh4QRF1TDQD5\ntarhd8eeWsJNAXh7w7XvFrHKoWoqMt7CfTeWEuwbjLenuqtGRak/REeF1kKSw5PJrsxmREQiIxym\nZprXYdTwt7+p6oJvv1U3zv/7P7j0UvjuO/WUvGsXvPOO+geZmwvC+wJ4HErMX0JL3k4wjcajOIGi\n/GQYCo3hu6EZlr1XzW1T1T8UsSiWP9y0itSUWU5/2cXVxUQERPDAA4IHHnB6WBtPPtn5Mb/8ZefH\nVFWppBVg9BY2Y/Yzk12ZzdCQodzzyT18m/Mtl5xxCfkF+V0LrptC/EIoqSkxEgTAnY8c5KLyKl74\nPpa0Y64bnyGEKkFcdFH7fTU16mbo7Ab3iJPZZwID4dxJJpatAS9TNdC263VPbpje3qq0M3GiKiVV\nVakEZj+XPfE0N7ferF95BV56ST1IREaq5JaerkpJ110HTz+tEoeHh0psFotKgpEt62+tWwdHjqjz\n/uEPKiGZTOr3UlOjksmIEa0POKCuu2YNrF+vXm/frq5x112qpLp9O6xerR5kgoPhvPPUA4wzdXWq\npOnhgtbYigr1QNnQoO5T6emqJCqEKmGfdZb6tx8VpR42w8Ndc93+5pYkIaVsEkL8AvgS1Xi+REp5\nwNnxQghC/EJYd2Qd3h7efPHjLzoc5ezYQFpRV0FUQBQ2WWP80ZXVF7Y53mKzEPHnCP50yZ94+ALn\ny2sWVheSEJxAQnACQ81dm3bD2xtmzVI/ji6+WP04qqkBb29fZvxrGp9m32tsv+eWOC45AzZlJzP9\nLahsVjdR/9AqJM2Iln4H+0v2ci0nSRI1xUT6d75SXl/oqKpiQswEdhbsJNGcyOvbXwcgNjCW9bXr\n+yQG+1gaR9vzt+Pl4aXadxDYmmxdWiiqN/z9e/5Ze8+vAmtBlyeQ7I4T/z/Z/3l5eMDQlj/5l1+m\nV268sePtAQHqJ7KDP9HgYJg7V/2ASl6bNqlS7B//qEq4Z52lSrSFhfDEE6qKMztbJaCkJJWoPvtM\nJZn8fJX44uJU7UFCgko6I0aoUpynZ2s13iWXqJLw99+r38HmzfDVV6pqMz5e1Q4kJ6skkJenkufM\nmeocX3yhahMaG9WDYEODqraMilLf4eyzVaIcO7bt97XZ4P331bVj3NNjvR23tUlIKdcCo7p6fFRA\nFD/+6McAXDS8g8c0WhtIvT1UVgjxC6G2oZa6xjomvDqBQ6WH2HHPDoaHDGdT9iYe/OJBQCWhN3e8\n2aYrraOy2jLC/cN55uJniPB3/Qpr9pvHyptWkvDXBGNQYFyQqm84b8h5vHPdO9y68lYAbvzgRqYl\nTmPDHRsAeHz949ww5ganXVxLakq6tJxqf5kQrZLEJWdcYmybmjCVh9c9TENTg1GqcxUvj/Z/5rd9\ndBuvXf0aQT5BBPoEYrVZCTOFufS6rmSxqcaLwurCPkkSg4UQqtprupOOeoWFqgo4Lk6VGHbuVNVj\nL7yg2pqSktSNODtbJY3MTFU9un+/Kv15eKjqwOxs+MlPVAKZNUuVZM45B95+WyW0AwdU4nHsCOHo\noYfab6uvV7UGHh4qgVx0kUoqkyap0th336mOEVLCbbep6j9PT3V8fb2qBj7rLFXi+uYbVYV55Ii6\nf8ybp9q1PD3Vca404Buu7W4ac1Ob9Rs6Yi9JjI0ay5bjW7A12Xhv73uU1JRwqFRVKZh9VXXHhJgJ\nHCk/AkB9Yz13rb6L+ePn4+XhxVdHv2J0xGhjsFtlXSVmXzPDQ4f34TdU62L7efkZSWJ4iLqel4cX\nN4650UgSAJuyN7UZV5BZnklcUByBPoH8/pvfc6T8CEuvXQqo6iZ3lSQ6Mj5mPO/tfY/H1z8OqGlJ\nRkeMJtg3mMLqQqfLovaUvRsxgPUxK0NfGEppbSlV9VUE+bYmidSlqSyZtYRz4s9x6fVdwT4t/d++\n/xtTE6a6OZqBKzoaFixofT9zZvtj/Pxan+CnOvwq/9wyj6i9FHX0qKoiOrENCODMM7sfm69va2eW\n3/0OHnhAJYuMDPjoI1XKePtt1bHEYlEdRKqrVVtMRYUqYbz7rioBXX21SgZBQSoJrlunSkLFxapK\nzZUGTZJIMCcAsOOeHU6PsVcrLLt2GeH+4QT8MYBF37Tt4mpfQyHIN8joH59ryQXUP8RQUyiXvn0p\nN4y5gQ9u/ACAqvoqY0BYX4sLiqOyvhJQs6fadfR07ThocPaK2dQ31SOflPxrz79IL003kkRJTUmf\nlIB6akjwEHKrcvn44Mf4ePqw/nZVzWT2M1NZV+nyJFHdUM1LV7zE/Z/fT4BPgFGyyKrIYkjwEAJ8\nAvhg3wfsKdrDuqPrBmSSsC9O9e99/2bFDe0nu9R678Qa7BN7J7paSAhO2wiDg1WDvaMZJxmadeK+\nW2/t+LieGDTNKPab+1nRZzk9xl49E+4fbkyeZ7coVSUL+83ecT4hI0nYWhcRcqyiqKyvNKaU6Gv2\nwYFbf7q13b4ls5Zg8mpNHDP/NZORYSOZHDvZKH18l/Md4aa260lbbJZ+S3JdEeEfYZTiHHvX2duS\nXO3Edb3t05XvKdpDXFAcJi8Tv173a4BurxT44f4Peeyrx5zuP1x2mJqGmh5E3VZJbQkrb1pJoE9g\nn/yONM2ZQZMk7DdtD+E8ZHu9u/1G+u8b/m3sWzB5ASF+IcZaCo4N37lVrSUJuxV7V7AwbSHQvyWJ\n6MBoACbHTW63786Jd1Lxm7Y3iERzIonmROP9BW9eQFSAmjbdfgO22qztRqa7U7ipdXS5YwnJ7Gs2\nSlGuVN1QzZyUOXx080cAvHD5CwBszN5IXFAcpbWt63zbqyW76on1T/DMt8843Zf0UhJ/3PjHHkbe\nqri6mKiAKEL9Qvt9XXTt9DZoksTk2Mk8Pv3xkx4zLGQY0JoA7D2RbhhzAzGBMZQ/2r4v/JmRZ7It\nfxvQfu6gv275K6BKEv2VJP4y8y98fbvzCf1O7IETbgrnvevf4/D9h9lxzw7MvmajJGGfC8pqs7Z5\nknY3x2VTU4elGq/NfmaXPyU3NTdR31iP2dfMtaOvBWD++PksTl1Ms2wmLiiON2e9aRzvOMtvV5TU\nlACQdiyt3b7Xt6meW56il8OvUT3UIvwjCPQJbPMwo2l9bdAkiQCfAJ6a8dRJj4nwj0A+2Vp9Yb+x\nvzX7rQ6Pz38ony13t9brW2yWNtUftiYbTc1NHK86blRl9bWk8KQ2N87OhJnC8PXyZUTYCMZFjaOy\nvpKSWnXjOlx2GBh4JQl7aW5k2Eij3QcgxNf11U3VDdX4e/u36zI9KkJ1rIsNiuXiMy42HjByqnLa\n/A1kV2afdJ2P6oZqQPU4A8gozeDLI18CGO1AHfWu6g4pJUXVRUQFRBHkG+SWiRC109egSRI9YW+X\ncPYUHRMYQ6BPIOcPUS1EVfVVbdZi9vf2p7yunCNlRwZct8Ols5cCbb+b/ea7+tBqwk3h5FnygIGX\nJOzCTW3bjuwN165UbavucN0PexuT/b9NzU2Aeup3rPIa+sJQFny6wOmCVvbODyU1Jbyw5QWS/57M\nZf+6DFBdp++eeDe/S/tdrxbEKqouwsvDixC/EKMnlqb1l1M6Sdj7vXe2vOimOzdxy7hbKKouwlJv\nwVN4svWnW4nwjyDPkke+Nb/Lg+j6w3lDzmP6UNVR3Fnf/pFhIwd0khhqHsq0xGlttrm64Trfkk96\naXqHDwlY9W3NAAAgAElEQVQpkSmYfc3G34bZTyULIQSfHPoEaE0AS3cuJfq56A4TmL0R/MzIM41x\nN/bPltaWGu1kvfleB0sOMjpiNEIIAn0CuWPVHS5pDNe0rjilk0SQb1Cb6idnPIQHI0NHUmAtwGKz\nMDRkKJPjJhNmCmP8q+OJCohy+QCv3th812bOCFX985wtJJRoTqTAWgAMzCRx9JdH+fOlf26zrScN\n182yud2KfnbjXx1P6rLUDscVJJoT23QCWHfbOjJ/mcmU+CmsP6a65G7K3mTsL64p5pusb9qcQ0qJ\nj6cPn93yWbulbOd/NJ9I/0gWpi4kwj+CV7a+YpRWuqustswY5+Lv7W9M9mgnFgkKrYVOPq1pvXNK\nJ4nuiAmMUUmivrW7qH29CXdMgtgVG+7YwG3jb2uzbVLsJJ679DmuTLrS6EFUXldOqKmL83r3Ew/h\n0a6E15OSxKfpnzLibx1XBRbXFAMn7zZtFxMYw7CQYdwx/g6jdLD60GoWpy7mX3P+xYJJC4ySmZ3F\nZsHLw4srk65kQswEAO6fotYweWfPO4yLHoePpw+NzY08vv5xPkn/pFvfza66obXKzJ5o7Kvp1TXW\ntfmumuZqOkm0sCeJqvoqYwyFvbH6muRr3BmaU9OHTm/X22nbgm08NPUhYgNjW5NEbXmH8xcNND3p\n3WS/aTY2Nzo9pjulqGDfYKNaaUfBDqbET+HWs24lNii2XZLYU7iHlIgUACbGTARax+MA/GTCTwA1\nohzUsq2//e9vu73IkdVmJdBbfQd7z7C0rDRS/pHCXavvAuCl71/qsIeVpvWWThItjJKEzWJM77H0\n2qWsn7+ev1/5dzdH131hpjDKasuob6wn35o/oOclshseMrzb4xTsvYtOXLnPsfTXnSRh9jMb4xDS\nS9MZHTHaiM1xeVmA/cX7jVKK/b+hplCWX7scwFi9cOXNK3nu0ud4e/fbPL3pafYV7aO4upjNOZu7\n9h0dGt/vmXwPoyNG8+WRLzlYcpB397wLwOvbX+eeT+/p8vfUtK7SSaKFY3WT40pqFw2/yOg1NJiE\n+4dTWlvKbR+p6ijHKT4GqtERo6msqzQWX+oK+zKvJ35mZ8FO43W3SxL1ldQ21FJSU2JMEXJ+wvnt\nburldeXGmJRQU6jR/jVv3DzOG3Jem/myRkWMMtbyXrZrGY989QhT35zKr9b+ql0MzbKZOz6+wygd\nVTdUG99hctxkdv9sNwBT4tsu5phemt6t352mdYVOEi2iA6ON6qaBNIVFT9lLEoOJEIKhIUPblAqk\nlCz4ZIHTUcb273hinbxjA3i3ShK+qhvumzveJNI/0nhAGBYyjFxLbpsSSlltmdErypGXhxeb79rc\nps3lipFXMG/sPJ695Fm25m012rv+8cM/2n3+2+xvWbZrmXHDX5+5vk0PLW9Pb+STkteufq3dZyvr\nKlnwyQLd+0lzGZ0kWvh7++Pr5UtOVU6beZ0GK7Ovqjb5YP8HPD/zeXeH02XxQfHcsvIW3tyhRkF/\nl/Mdb2x/g41ZGxGLhJFA7E/Z9ik17COf7cpry42qou4kiQj/CIprivnF578gp6o1Wfl4+uDn5Wck\nq41ZG3n222e7PKeXp4cn717/LnNGzyG7MtuIv7G5kYzSDLblbTOObV17u5jG5ka+yfqG9NL0dueM\nD4oH4KeTfsr3d38PwL7ifbyx/Q2yK7O7/J017WR0knDg7eHN27vfNtokBjPHp9jBlPTig+JJL003\nGmS/z1U3v6vfuxqAWStm8d6e9/D+veqSXFZbRkxgDPmWtqvaldWWcW78uQDtJns8mSDfIGN+sBPH\ncQgEIc+qhuOnNqrR/w3NzkdjdyTBnECeJY8f8n5g1dxVACT/PZmz3zibqUum8uKWF7ll5S2AKh3Z\nk+K8ce2XQLSP6Pb39mdK/BQuSLiA3YWqKqraVt2tuDTNGZ0kHJTWlpJdmT2g5jlyhc4GEw4k9uob\ne33+kh1L+GRea9fRnQU7Wb57ufG+tLaU2aNmG2Mb7MrrygkzhWF9zGp0T+0qe6+2D2/8sM12+yzB\nUkrKa9U8YPb/dpWflx8TYiZwrOJYu/Ebm49v5k/f/cl4f8171/D0pqeZnji9zQJNdvb/r/ZxIoE+\ngcbYGD0JoOYqOkk4WDBJrVbiOCvoYDYpdhLQOt5jMLBPhldeV05jcyNZFVlcOPRCVs1dRVJYEqDq\n6O3KasuYMXwGGaVtex4VWguJCojqcEqOztgHKJ7Y3uDn5QeosQlWm5VHL3iUn5/z826f3z4Ku6Me\nZ47dbOsa63hj+xtcnXy103M9Mf0Jo6ttoE8ghdVqUJ1OEpqr6CTh4LVrVEOg/WYw2G1boOq5B9pA\nupOxV/VE+keSU5lDbWMtgT6BzBo1i0cveBRonQqjuLqY9NJ0JsZMpLC6sE2j8oGSA4wK7/LquG3Y\n2zBO/DsYGTYSUDdgq83Kvefc26OJHx2nvW/4vwYeOr/9WpeOPaPsA/Q68vsZv2dOyhwj7pUHVgL0\nyZTr2sDQ3+uJ9CpJCCFuEELsFUI0CSEmnbDvMSFEhhDigBBipsP2SUKI3UKIdCHEC725fl8ofaS0\n0ynJB5PCXxcyb2z7+uyByt6baHzMeDZkbSDIJ8ioVrlz4p1GY7Sflx8bsjYwJX4KSeFJeHl4tVk0\n6kDJAcZEjulZDE6m9l5761pA3YB7M9XJgkkLuH387YBKGCdOTwKqMdqeFLvafdlxckpXT5SoDQwb\nsjYQ+myoMdK+P/S2JLEHmAO0mdRGCJEC3ASkAFcAL4vWivFXgLuklMlAshDisl7G4FJhprABNU9T\nb0UFRA2qNon7p9zPu9e9y5iIMfyQ90Ob7shCCKNrp9nXzObjmxkRqqbkiAqIMrqM1jbUkmfJ6/HM\nvc66DscHxzMpdpJRkuhpkrgi6QpjaVlobVu4bMRlJAQnUPBQAU/NeKrb1YQLL1xovHY2n5U2eBVX\nF3Ph0gsBNdq/v/QqSUgpD0kpM4AT70KzgRVSykYp5TEgA5gihIgBgqSUP7Qctxy4tjcxaKeWyIBI\n5o2bR1RAFEfKj7RrFxgTOYY5o+cwZ/Qcnt/8vDHRoX18w8GSg5z58pkMNQ/t8ToOieZEY3W/E5l9\nzRRXqzEZJ06J0ltCCLIfzCY6MBohBJeccQnjosZ1+fMpkSnG673Fezs85qYPbuKuVXf1Olat/zl2\n896a1355477Su9VQnIsHHIen5rZsawSOO2w/3rJd09qI8I/gq6NfGd1Y7dbcsgZQCyq9uu1VhocM\nB1TXVYvNwoq9K8isyGzXfbU7ll671OlcUMG+wRyvOt4v3aSvSr6Kq5Kv6tZnll+7HF8vX27+8GYW\nf7OYeybfYyyJC2pNbolkyewlrg5X62OOA0Z3FOzo8JgP9n3A6vTVLr1up0lCCLEOiHbcBEjgcSll\nz6a17IaFCxcar1NTU0lNTe3rS2oDgJ+XH43NjcbcTHb2qhl7VZL9iT/YN5iq+iqOW9QzyMlWk+uM\nj6eP01KC2c/MsYpjA3YuLPuswDd/eDNPpj3JqkOrjA4MoKpTT5Xee6ebt3a2rrB54uDRtLQ00tLS\n2Jq3td2Yod7qNElIKS/twXlzgQSH90Natjnb7pRjktBOH/b1GeaeObfD/R7Cg5U3reT8BLWqYJBP\nEJZ6izEy+cR/RK4S7BNMZkWmMWeTq6y8aWWPG9pPxt6+IhYJdv1sFyF+ITpJDDLT35rOb6f91pjl\nd97Yee2mobE/QD+z6RnKa8vZvmK7y67vyi6wju0Sq4G5QggfIcRwYCTwPyllAVAphJjS0pA9H1jl\nwhi0U8SoiFHIJyWPTnvU6TFzUuYYT/xGSaJlSou+uhGa/cxkVmS6vCQxJ2WOse62Kzx36XOAmkXA\nbkPWhl6vt631v03Zm3hv73tEBUSx5a4t/GLKL7DUWzo8trzW9WvH9LYL7LVCiBzgPOBTIcTnAFLK\n/cD7wH5gDXCvbO3Efh+wBEgHMqSUa3sTg6aBShLFNcWU15bzxjVv8I8r20+c56rrZJa7Pkm42kNT\nH2LrT7e26T6bWZ45qHq6aa2LShVVFxlTxgf5BLXp7u2ovK7cWHPEVXr1WCGl/Bj42Mm+p4GnO9i+\nDeh6lw1N64KLh1/MnavvJMwUxt2T7u6z6wwLGUZhdWGbwW4DVaI5kdyqXGNhpuOW43qZ00Em35KP\nQLA9fzv+3v5Gt2tnJYmKugqXLzCmR1xrp4QZw2eQZ8ljWMiwPr3O9SnXAx1PqTHQRPhHYLVZjaq3\n9/e9T3mdmmtqoC7Jq7VVYC3g7LizqW2sJacqh0CfwE5LEq6ubtIVlNopwdvTm00/2dRmrEBfsI8I\n7+7sr+4ghCAmMKbdvFYA9U31p8z0M6eq4S8OZ2zUWGKDYrE12dhVuItAn0A8hAeWegtSynbVh32x\nVLEuSWinjAsSL+iXJ/xxUeNIHZba59dxhXD/cL7L+c6Y2TghOAGzrxmrzermyLTOHKs4xqfpnzIy\ndCRzRqv5uXw9ffHx9MHTw7PDqTkGXJuEpp2Odv98t7tD6LIwUxiPfPWI8f7KpCtZk7EGq81qrEeh\nDUxxQXGEm8J54NwHiA6Mprax1ig52KucTpzXa8D1btI0bWCzT1b44uUvAmDyMjEsZBhHyo64Myyt\nC5plM2t/vJahIUPx8/LjmUueMfYF+Qa1a7y+e/XdfVKS0ElC005h9inDHzj3AQACfAKYGDORnQU7\n3RmW1gWVdZVtJrh01FHj9ZIdaqoVV4+F0dVNmnYKe/HyF401Ou6fcj/3nXMfH+7/kIMlB90cmXYy\nDU0N2JpsTlfJPLEk0Ze91XSS0LRT2JT4Kcbrv13xN0CN9Vh7RI9hHcgsNguBPoFOBz+eWJLoy/Ul\ndHWTpp1mRkeMZmPWRt7Z/Y67Q9GcqG2oxd/b3+n+cP9wSmtap56xr7PSF3SS0LTTTFJ4EhabhR9/\n9GOklBwuO+zukLQT1DbWnnRFwuiAaGM9c9BJQtO0PvLRwY9IeinJ3WFoJ6htqMXk5TxJhJvCWfTN\nImP1wpqGGmIDY8m4v/3Ayd7SSULTTkMrrl8BwPXvq2lGmmWzO8PRTtBZScLPyw+rzcpHBz8CVJKI\nDoxmZNhIl8eik4SmnYZuOvOmNu/tS7JqA0NnJYlfnfcrfnXur4z1U2oaak7ahtEbOklo2mnoxF4z\nv9/wezdFonWks5KEp4cnZ0WfRUZZBhV1FcYssX1BJwlN0/jHD32z/obWM52VJEB1QMgozeD+z+/n\ngbUP6CShaZprfXvntwT5BHHkgSNGHbc2MNQ11p20JAGQFJZERlkG1Ta1DryrZ3+100lC005TUxOm\nUvVYFWeEnsHIsJF6FPYAUlpb2mlJIiogioamBmNQ3aTYSX0Si04SmqYxLWEam7I3uTsMDdhZsJP7\nP7+fa5KvOelxQgiSwpPIqsgC4Oy4s/skHp0kNE1jaMhQ8i357g5DA17+4WUArkq+qtNjR4SOILsy\nG4AJMRP6JJ5eJQkhxJ+EEAeEEDuFEP8RQgQ77HtMCJHRsn+mw/ZJQojdQoh0IcQLvbm+pmmuEegT\n6HRJTK1/JQQn8MjUR/Dx9On02KHmodQ31bN09tIB23D9JXCmlHICkAE8BiCEGAPcBKQAVwAvi9Y+\nd68Ad0kpk4FkIcRlvYxB07ReCvIJctpwfbjsMGJRxxPNaa5XVV9FuH94l44dGjIUgJjAmD6Lp1dJ\nQkr5lZTGUM0twJCW17OAFVLKRinlMVQCmSKEiAGCpJQ/tBy3HLi2NzFomtZ7Qb7t1yewW/zN4n6O\n5vRmsVkI8gnq0rGJ5kQAogOj+yweV7ZJ3AmsaXkdD+Q47Mtt2RYPHHfYfrxlm6ZpbhTk036lM7u3\nd78N9O2aBVqr0tpSgny7liSGmlVJIiogqs/i6XQ9CSHEOsAxTQlAAo9LKT9pOeZxoEFK+Z6rA1y4\ncKHxOjU1ldTUVFdfQtNOe11pk3joy4cI8Qvhdxf+rp+iOv08vfFpPtz/IbNHze7S8fbqpv3/28/r\nG1/vk5hEb58OhBB3AD8FZkgp61u2/QaQUspnW96vBZ4EsoCvpZQpLdvnAhdKKX/u5NxSP71oWt87\nUHyA2Stmk35/epvtw18czrGKY8b7ybGT2bpgaz9Hd/qwt/3su3cfYyLHdOkze4v2MjZqbNvzCIGU\n0iUNSb3t3XQ58DAwy54gWqwG5gohfIQQw4GRwP+klAVApRBiSktD9nxgVW9i0DSt94aHDienKofG\n5kZjW54lz0gQ3h7eQN91szzdvfT9S3x88GPCTGEcfeBolxME0C5BuFpvly99CfAB1rV0XtoipbxX\nSrlfCPE+sB9oAO51KBLcBywF/IA1Ukq9jqKmuZmflx8xgTEcqzhmTDd9tPxou+NqG2v7O7TTwgNr\nHyAhOAFLvcWoQhooepUkpJROVyuRUj4NPN3B9m3AuN5cV9M014sJjKG4upiRYSOx2qw8v/l5pidO\nZ/74+fxizS+Avl0B7XSXa8klKiAKDzGwxjgPrGg0TXObMFMYZbVlAEx5YwofH/yYH/J+4O5Jdxs3\nrtoGXZLoK82ymWDf4M4P7Gc6SWiaBqhZRO1J4kDJAQC2/lQ1UtuThC5J9C2dJDRNG7AcSxJ2Z0ad\nCbQmCd0m4XpNzU3G60CfQDdG0jGdJDRNAyDCP4LC6kKaZTOewpP99+439t0w5gYSzYm6JNEHrDar\nMcJ6IK41rpOEpmkAjIkcw77ifZTVlmH2M5MSmWLse3P2m6TdnqaTRB+oqq8yRlgPxDYfnSQ0TQMg\nJSKFQyWHyKrIIi4ort3+IF/nU3doPWexWQj2DcbkZWJy7GR3h9NOb8dJaJp2igj3D6eiroJvc75l\n6pCp7fYH+wZTWV+JlJLWSZ213qqqryLYN5jih4v7bLrv3tAlCU3TADD7mqmoq+BYxTGSwtsPgfLx\n9KGxuZHjVcc7+LTWU/YkEeATMCCTr04SmqYBatS1RJJTlUOEf4TT44a9OKz/gjoNVNVXdXlqcHfQ\nSULTNEBNChfiF0JGaQaR/pFOjxuIPXAGM0u9ZUCOj7DTSULTNEOIXwiHyw4TGeA8SQCsPrS6nyI6\n9RVWF5605OZuOklommYYFjKM6oZqpyWJj27+CIDZK7q23oHWuf3F+7s162t/00lC0zRDSoQaG+Gs\nJDEQRwQPdnuL9nJm5JnuDsMpnSQ0TTOkRKTg5+VHgHdAh/udbdd6pqm5iUOlh3RJQtO0wWF0xGgi\n/SOddsW09+PXycI1ci25hPiFdHlNa3fQSULTNMN5Q87j5atedrrf5G0CoEk2OT1G67oCa0GHo9sH\nEp0kNE0z+Hr5cnXy1U73J4cnY33MSmNzI7YmWz9GdmoqsBYQExjj7jBOSicJTdO6JcAngGDfYKrq\nq9wdyqCXb8knOiDa3WGclE4SmqZ1m9nXTGVdpbvDGNQamhr42Wc/Y3z0eHeHclK9ShJCiMVCiF1C\niB1CiLVCiBiHfY8JITKEEAeEEDMdtk8SQuwWQqQLIV7ozfU1TXOPYN9g1meuZ03GGneHMmjZ58C6\nLuU6N0dycr0tSfxJSjleSjkR+Ax4EkAIMQa4CUgBrgBeFq3dJV4B7pJSJgPJQojLehmDpmn9zMvD\niwWfLuCqd69ydyiDVlV9FWdFn0V8cLy7QzmpXiUJKaXV4W0AYJ/UZRawQkrZKKU8BmQAU1pKGkFS\nyh9ajlsOXNubGDRN63+Hyw67O4RBr7K+ckDP2WTX6zYJIcRTQohs4Bbgdy2b44Ech8NyW7bFA47z\nDB9v2aZp2iDy+jWvt3lfXF2sJ/7rpqr6Ksy+ZneH0alOFx0SQqwDHJvfBSCBx6WUn0gpnwCeEEI8\nCtwPLHRlgAsXtp4uNTWV1NRUV55e07QeuOnMm7j5w5sBNSts1HNRrLh+BTePvdnNkQ0elXWuK0mk\npaWRlpbmknOdqNMkIaW8tIvnehfVLrEQVXJIcNg3pGWbs+1OOSYJTdMGnuLqYgAamhvcHMng4srq\nphMfoBctWuSS80LvezeNdHh7LXCw5fVqYK4QwkcIMRwYCfxPSlkAVAohprQ0ZM8HVvUmBk3T3Ovx\n9Y8DUG2rdnMkg8u+on0khye7O4xO9XaN62eEEMmoBuss4GcAUsr9Qoj3gf1AA3CvlFK2fOY+YCng\nB6yRUq7tZQyaprnRkh1LACitLXVzJIPLtvxtzB07191hdKpXSUJKecNJ9j0NPN3B9m3AuN5cV9M0\n99txzw6W71rOX7f8FYCSmhI3R9Sx2Stm86PEH/HQ1IfcHUobWZVZnBF6hrvD6FRvSxKapp2mJsRM\nYELMBGoaanht22sDtiSx+tBqCqwFAypJ1DXWUVZbRmxQrLtD6ZSelkPTtF559epXWT13NaU1pewr\n2ufucDpk8jK5O4Q2jlcdJz4oHg8x8G/BAz9CTdMGvHD/cD7L+Iyxr4xlb9Fed4fTjsnbRENTA03N\nA2OK86yKLIaGDHV3GF2ik4Smab0W4hdivN5fvN+NkbRlnx/J5GVi4msT+fFHP3ZzREpWZRaJ5kR3\nh9ElOklomtZrw0OGG69zq0469Klfrc9cD4DZz8y+4n2s2LvCzRHBsYpjvL7tdSbFTHJ3KF2ik4Sm\nab1m8jax5pY1JIUlGU/vJ7LarP0+dUdRdRHhpnCampuI9I/s12s78+dv/8z+4v3MGzfP3aF0iU4S\nmqa5xBVJV7AwdSF/2fIXvsv5rs2+Hfk7CHo6iKU7l/ZrTMerjpMUnkRNQw0+nj79em1ndhTs4JN5\nnxAVEOXuULpEJwlN01xmSPAQAH7IVRM93/7x7ewu3M2k11XVSn+Pys615JIcntwmSbSO63WPkpqS\nAb9kqSM9TkLTNJexPx1bbBbEIrWEzIjQEcb+QJ/Afo3neNVxrht9HR+WfEh1g0pQFXUVhJpC+zUO\nuwJrAWW1ZYSZwtxy/Z7QJQlN01xmVPgopiZMZWP2RmPbxwc/Nl5bbVb+l/u/fovneNVxfjT0R6SX\nplNWW0ZcUByV9e5ZdjW9NJ3Y52PdmqR6QicJTdNcRgjB1UlXt2mT2FGww3h9qPQQ5/7zXDbnbO6X\neMpqyxgTOYaKugoAQv1CsdqsnXyqb9Q21ALQJJvw8hg8lTg6SWia5lKxQbFYbVZ+O+23LJi0AIA1\nt6zhielPGD2fpr45tc/bBpplM7UNtQT4BDA6YjSNzY0E+gRiqbf06XWdsdjUdWeOmOmW6/eUThKa\nprnUyDC1gkBUQBR51jwALh95ORH+Eaw61LoyQF1jXZ/GUddYh6+XLx7Cg7Tb09hwxwaCfIPcVpKo\nqq9iWuI01t46uCa+HjxlHk3TBgX7GgmhplAmxkwkszwTIQSRAW3HKVTUVWDy7rs5lWoaagjwDgAg\nOjCa6MBoVZKw9X9JQkrJ2sNriQ+KRy2lM3jokoSmaS5l9HCqt7D4osXsvVfN5XRBwgVtjrO3E4Aa\npW1rsrk0jmpbNf7e/m22BfkE8cXhLxCLBGKR4J3d79DY3OjS63ZkZ8FOXvrfSwT5BPX5tVxNJwlN\n01xu1dxV7UYUDw0ZStkjZfzmgt/g7eFNRV0F/z36XxZ/s5ghfx3CHzb8oVvXKK8tP2lPqZqGmnZJ\nItg3mGW7lhnvf/zRj9lTuKdb1+2J+qZ6ADw9PPv8Wq6mq5s0TXO5WaNmdbg91BTK05c8zf6S/eRa\ncnlzx5t8fvhzoG3Joit+9cWvWL5rOfLJjhvAaxpqCPAJaLMtPijeuGHbdfe6PVFao9bauHDohX1+\nLVfTJQlN0/rdjxJ/xFdHv2qznkJ32ydyKnNOur+6oX11U0czr+ZZ8rp13a6y1FuMHlxltWXcMOaG\nQTNfkyOdJDRN63eTYifx2rbX+CzjM2Obr6dvt87hbCJBu4q6Csy+5jbb7D2v7OaPn0+BtaBPJh4M\nfiaYr45+BUBmRSbxQfEuv0Z/cEmSEEI8JIRoFkKEOWx7TAiRIYQ4IISY6bB9khBitxAiXQjxgiuu\nr2na4DI6YnS7bbWNtXx55EvyLfmdfl5KSUZZxkkbgvMsecQFxbXZNjF2IgCHfnEI2xM2ksOSOVZx\nDM/FnhRYC7r5LZyzz1FV01ADwKpDq5gzeo7Lzt+fep0khBBDgEuBLIdtKcBNQApwBfCyaO339Qpw\nl5QyGUgWQlzW2xg0TRtcHNd2jg1Ur6vqq7jsX5fxxPonOv28fazDyaa3yK3Kbff07uPpQ86DOSSH\nJ+Pt6U2YKYy///B3QE0r7iq7CncB6jsBFFcXD5qV6E7kipLEX4GHT9g2G1ghpWyUUh4DMoApQogY\nIEhK+UPLccuBa10Qg6Zpg0zG/RkALEpdxNpb1/Jp+qcAWBs6H+xWUlOCv7f/SWeVzbXktitJQOtM\ntaCWXQU1XUd5bXm34j+ZnQU7AZj/8XxuXXkrpbWlhJvCXXb+/tSrJCGEmAXkSClP7EMWDzi2KuW2\nbIsHHCsSj7ds0zTtNGNvH8i15HLpiEvJteRy3pDz+PLIl+zI33HSz5bUlJBoTjRmdu1IdmV2p0uE\n2tsipiVOo7zOdUnCsVTy773/NqYEGYw67QIrhFgHRDtuAiTwBPBbVFVTn1m4cKHxOjU1ldTU1L68\nnKZp/ejuiXdzZdKVRi+npLAkQvxCmPT6JGofr8XPy6/Dz5XUlJAQnMChkkM0y+Y2vaTscqpySDAn\nnPT6s0fN5rs7v+PVba+y9vBaHvziQTJ/mdnr71VtqyYhOIGcqhyCfIMweZn6dKR1WloaaWlpfXLu\nTpOElLLDJCCEGAsMA3a1tDcMAbYLIaagSg6OKXxIy7ZcIKGD7U45JglN004tb8x6w3i98qaVTI6b\nzH1r7gMgqyKLURGjOvxcSU0JUQFRnBl1Jl9nfs3FZ1zcZv+ewj2kl6aTEHzyJGHyNnF+wvl8uP9D\n/vHDP9qNoeip6oZqEs2J5FTl0NDU0KaKqy+c+AC9aNEil527x9VNUsq9UsoYKeUZUsrhqKqjiVLK\nIvryWz0AAA3RSURBVGA1cLMQwkcIMRwYCfxPSlkAVAohprQklvnAKqcX0TTttDEnZU6b6qHMCudP\n9MU1xUT6RzJj2Ax2F+5ut/+sV88CIMi3a9NgnDfkPJclCGhNEqB6bQ3W9ghw7TgJiaqKQkq5H3gf\n2A+sAe6VrfMC3wcsAdKBDCnl4JoSUdO0PmW/Vdzx8R1szdva4TElNSVE+EcQGxRLvrXzLrOdmT16\ndq/P4chqsxLqp3peNctmo4F8MHJZkmgpUZQ5vH9aSjlSSpkipfzSYfs2KeU4KWWSlPKXrrq+pmmn\nlsLqQv65/Z8d7jOSRKDzJHHx8Is73N4RH08fdv1sV4/i7Ei1rZoxkWOM9xGmCJedu7/pEdeapg1Y\nkf6RHW4vsBYQFRDF5LjJrDuyrs1Kd1JKvDy8+PzWz7t1reTw5G6P+namuqGacdHj2L5gO9A6M+5g\npJOEpmkDyqTYScbrEyfoA7UMaGZFJkOChzA2aiyF1YVMfXOqMc6hpqEGbw9vvD29u3VdX09fbE02\nl0zRYam3EOgTaMwd1VlX3IFMJwlN0waUhakLsT1hY+GFC9sNliupKcH/j/7sLdpLfLAaYvXa1a8B\nGOMcKusrCfYN7vZ1hRCYvE3GWtS9kW/NJy4ojpjAGADCTGGdfGLg0klC07QBxUN44O3p3eFSo+sz\n1wOQEJxAdIAavrVg8gImxkw0ShJltWWY/dpO7NdV/t7+1DTU0NTc1OkEgs7UN9ZTXltOVEAUZj8z\nb81+ixnDZ/ToXAOBThKapg1IQT5BxlKjG7M2kmfJ40jZER6e+jDZD2a3WcAnzBRmlCTGvTKO9NL0\nHl3TvrzpM5ueIeGvCbR2yuy6AmsB0YHRxgC/OybccdI5pgY6nSQ0TRuQzH5mympVh8kfLf0Rt398\nO5kVmQwPGd7u2FBTKGW1ZdQ39m6sQ1xQHLlVuewoUNOCHCk/0u1zWGyWHlV3DVQ6SWiaNiCdFX0W\nOwt2YqlXpYn00nQyKzI5I/SMdsfGBMSwbNcyQp8NJdI/ksrfVPbomonmRLIrs9lfvB+Tl4ni6uJu\nn6OjZVMHM50kNE0bkJLDkympKSH4GfVUnmfJ43DZYYaHti9JJIUnsSZjDbWNtTQ0N/T4ST4hOIEj\n5Uc4Wn6UyXGT+c+B/3T7HDpJaJqm9QMP4cHZcWcD8KdL/oS/tz/HKo4x1Nx+XYbk8GTjdW96JyUE\nJ/DfzP+SaE5kU/Ymnt/8fLfPoZOEpmlaP/nJhJ8AcPnIy6lrrAPA16v9gLeksCTjdUczwnZVojmR\nDVkbODPqTEaEjujROWobanWS0DRN6w/2gXVhpjBsTTanCcC+6tsVI6/gjxf/scfXs1dlpUSksOfn\napmchWkLu3WOmoYaTF6mHscw0HQ6VbimaZq72KfYtg9Gc6xWcuTl4UXxw8VE+PdujqSUiBQAJsdO\nxuStbvQ/5P1wso+0Yam3sP7Yel2S0DRN6w/BvsG8d/17mLxNpN2exrrb1jk9trcJAsDb05t3r3uX\nq5OvBuDVq15tt072ybz4/Yss3bmUxubGXscyUOiShKZpA5YQgrlj5wJw4bAL++Wa88bNM16HmcKM\nsRpdYZ9GxHEG2MFOlyQ0TdOcCDOFUWAtwNZk69LxNQ01vHDZC/x66q/7OLL+o5OEpmmaE5EBkXyb\n8y2XLL+kS8efat1fQScJTdM0p0ZHjAZgY/bGLh1f06iThKZp2mnDx9OH5y59Dm+Prq1NoUsSJxBC\nPCmEOC6E2N7yc7nDvseEEBlCiANCiJkO2ycJIXYLIdKFEC/05vqapml97f+d///w9PBsN215R3SS\n6NhfpJSTWn7WAgghUoCbgBTgCuBlIYRoOf4V4C4pZTKQLIS4zAUxaJqm9QkhBHFBceRbOl5L25FO\nEh0THWybDayQUjZKKY8BGcAUIUQMECSltI9OWQ5c64IYNE3T+kyEf0SXusLWNNQYg/BOFa5IEr8Q\nQuwUQvxTCGFfDioeyHE4JrdlWzzguNzT8ZZtmqZpA1aoX2iXkkRVfdUptZYEdCFJCCHWtbQh2H/2\ntPz3Gvj/7d17jFTlGcfx768C3a1U3P1DNhVEqRIxkhKJSNy2bGxFS63tPzWkFyltmkZNbVbTqk29\n/Cc01KpJbWJbwVuxSmnBlAgaXRNIaLG6gkIFgroIYesF1wjFIjz947y7O5o97GV2dqZ7fp+E5Mx7\n5sy878PsPHMuz3u4B5gaETOB/cDgp0w0M6txAy2q6zrcxYRPDu3WqbWq34rriLh4gK/1O+DxtLwX\nmFyyblJqy2vPddttt/Ust7S00NLSMsDumJkNj4a6hp7bo+ZpWtrEm4feHPL9tcvR1tZGW1tbRV5b\nQ7mHa8/GUlNE7E/LrcD5EfEtSecADwMXkB1OehI4KyJC0ibgWmAz8Dfg7u4T3n28fpTTPzOz4XDz\n0zcz9oSx3DL3lj7XH/zvQcbfPh6AuLX631mSiIi+zhcPWrlzN/1S0kzgGPAa8COAiNgm6VFgG3AE\nuLrk2/4aYDlQB6zNSxBmZrWisb6Rjq6O3PWt61pHsDcjq6wkERFXHmfd7cDtfbT/E5hRzvuamY2k\nhvoG2jvbc9fvemcXd8y7g+vWXzeCvRoZrrg2M+tHY30jB/6TnZO4c9OdHymsO3rsKM/te46FMxfW\nxKGm4eYkYWbWj9JLYFvXtbKxY2PPuln3zuLQkUM9N0YabXw/CTOzfnRfAvveB+8B2ZxOi1YvYsqE\nKbzY+WKVe1dZThJmZv1orG/kwOEDvHrgVQC6Puhieftymic3A7DuO+uq2b2K8uEmM7N+NNQ3sP/9\n/SzZuASAdw+/C8DGPdlhp3mfnZe77f87Jwkzs37UjakDYMVLKwDY07XneE8fVXy4ycxsgBrrG2me\n3MxDWx+qdldGjPckzMwGYOnFS1l1xSoum3YZO97ewcpvrhx18zT1xXsSZmYDcP2F1wP0zOE0/6z5\nbPj+Bo7FsWp2q+KcJMzMBmHySdkcpfVj6zn3lHOr3JvKK2uCv0rzBH9mZoM3nBP8+ZyEmZnlcpIw\nM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkzMwsV9lJQtKPJW2XtFXS4pL2myTtTOvm\nlbSfJ2mLpB2S7iz3/c3MrHLKShKSWoCvATMiYgawNLVPB64ApgNfAe6R1F3991vgBxExDZgm6ZJy\n+lAUbW1t1e5CzXAsejkWvRyLyih3T+IqYHFEfAgQEW+l9q8Dj0TEhxHxGrATmC2pCfh0RGxOz3sA\n+EaZfSgE/wH0cix6ORa9HIvKKDdJTAO+KGmTpGckzUrtpwKld+XYm9pOBd4oaX8jtZmZWQ3qdxZY\nSU8CE0ubgAB+kbZviIg5ks4HHgOmVqKjZmY28sqaBVbSWmBJRDybHu8E5gA/BIiIxan9CeBW4HXg\nmYiYntoXAHMj4qqc1/cUsGZmQzBcs8CWez+JvwIXAc9KmgaMi4i3Ja0BHpZ0B9nhpDOBf0RESOqS\nNBvYDFwJ3J334sM1SDMzG5pyk8Qy4D5JW4EPyL70iYhtkh4FtgFHgKtLbgxxDbAcqAPWRsQTZfbB\nzMwqpKZvOmRmZtVVkxXXki6V9K9UcHdDtftTaZImSXpa0supKPHa1N4gab2kVyStkzShZJs+ixVH\nA0mfkPR8OmxZ2DgASJog6bE0vpclXVDUeEhqlfRSKsZ9WNK4osRC0h8kdUraUtI26LEPqZg5Imrq\nH1ni2gVMAcYC7cDZ1e5XhcfcBMxMy+OBV4CzgSXAz1L7DWQ1KQDnAC+QHS48PcVL1R7HMMajFXgI\nWJMeFzIOaYzLgUVpeQwwoYjxAD4D7CY77wnwJ2BhUWIBfB6YCWwpaRv02IG/A+en5bXAJf29dy3u\nScwGdkbE6xFxBHiErDhv1IqI/RHRnpbfB7YDk8jGfX962v30Fh5eTh/FiiPa6QqRNAmYD/y+pLlw\ncQCQdBLwhYhYBpDG2UVB4wGcAJwoaQxQT1Z/VYhYRMQG4MDHmgc19qEWM9dikvh4IV6hCu4knU72\ni2ETMDEiOiFLJMAp6Wl5xYqjwa+Bn5LV4nQrYhwAzgDekrQsHX67V9KnKGA8ImIf8Cugg2xcXRHx\nFAWMRYlTBjn2IRUz12KSKCxJ44GVwE/SHsXHryoY1VcZSPoq0Jn2qo53+fOojkOJMcB5wG8i4jzg\nIHAjBftcAEg6meyX8xSyQ08nSvo2BYzFcVRk7LWYJPYCp5U8npTaRrW0C70SeDAiVqfmTkkT0/om\n4N+pfS8wuWTz0RKjZuBySbuBFcBFkh4E9hcsDt3eAPZExHPp8Z/JkkbRPhcAXwZ2R8Q7EXEU+Atw\nIcWMRbfBjn1IManFJLEZOFPSFEnjgAXAmir3aSTcB2yLiLtK2tYA30vLC4HVJe0L0tUdZ5CKFUeq\no5USET+PiNMiYirZ//vTEfFd4HEKFIdu6VDCnlSoCvAl4GUK9rlIOoA5kuokiSwW2yhWLMRH97AH\nNfZ0SKpL0uwUwytLtslX7bP2OWfyLyW7wmcncGO1+zMC420GjpJdyfUC8HyKQSPwVIrFeuDkkm1u\nIrtqYTswr9pjqEBM5tJ7dVOR4/A5sh9O7cAqsqubChkPsql9tgNbyE7Uji1KLIA/AvvIipY7gEVA\nw2DHDswCtqbv1rsG8t4upjMzs1y1eLjJzMxqhJOEmZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4n\nCTMzy+UkYWZmuf4HeWNyBvF5vcoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6d2437a790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cointegrated,x,y = GenerateData(1000)\n",
    "plt.plot(x)\n",
    "plt.plot(y)\n",
    "print \"Cointegrated?\",cointegrated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def LinearRegression(x,y):\n",
    "    slope, intercept = polyfit(x, y, 1)\n",
    "    std_eta = std( y - intercept - slope * x , ddof=1 )\n",
    "    return slope, intercept, std_eta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(19.6297413653121, -48.422956376420203, 142.98613260755775)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LinearRegression(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a key result from the paper: calculating the moments and area as derived in the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def CalcLogAreaLog(logf,logF):\n",
    "    lncdf = norm.logcdf([1,-1], loc=exp(logf).real, scale=exp(0.5*logF).real)\n",
    "    logarea = logsumexp([lncdf[0], lncdf[1]],b=[1, -1])-log(2)\n",
    "    return logarea\n",
    "\n",
    "def CalcMomentsLog(logf,logF,logarea):\n",
    "    lnpdf = norm.logpdf([1,-1], loc=exp(logf).real, scale=exp(0.5*logF).real)\n",
    "    logmoment1 = logsumexp([lnpdf[0]+logF-logarea, lnpdf[1]+logF-logarea, logf],b=[-0.5, 0.5, 1])\n",
    "    logmoment2 = logsumexp([logF+logf+lnpdf[0]-logarea, logF+logf+lnpdf[1]-logarea, \n",
    "                            logF+lnpdf[0]-logarea, logF+lnpdf[1]-logarea, 2*logf, logF],\n",
    "                           b=[-0.5, 0.5, \n",
    "                            -0.5, -0.5, 1, 1])\n",
    "    return exp(logmoment1).real, exp(logmoment2).real"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now inference: filtering and the EM update routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Filtering(V,std_eta):\n",
    "    T = V.size\n",
    "        \n",
    "    # DIRECT METHOD\n",
    "    logft = log(sum(V[1:]*V[:T-1]))-log(sum(V[:T-1]**2))\n",
    "    logFt = 2*log(std_eta) - log(sum(V[:T-1]**2))\n",
    "    assert isreal(logFt), \"logFt must be real\"\n",
    "    logarea = CalcLogAreaLog(logft,logFt)\n",
    "    loglik = -0.5*log(sum(V[:T-1]**2))-0.5*(T-2)*log(2*pi*std_eta**2)+logarea  \\\n",
    "                        -(sum(V[1:]**2)-sum(V[1:]*V[:T-1])**2/sum(V[:T-1]**2))/(2*std_eta**2)\n",
    "\n",
    "    # calculate moments\n",
    "    moment1,moment2 = CalcMomentsLog(logft,logFt,logarea)\n",
    "    \n",
    "    return loglik,moment1,moment2\n",
    "\n",
    "def EMUpdate(x,y,moment1,moment2):\n",
    "    T = x.size\n",
    "    xt, xtm1 = x[1:], x[:T-1]\n",
    "    yt, ytm1 = y[1:], y[:T-1]\n",
    "    \n",
    "    # find the coefficients\n",
    "    a = 2 * (T-1) * moment1 - (T-1) * moment2 - (T-1)\n",
    "    b = moment1 * sum(xt+xtm1) - moment2 * sum(xtm1) - sum(xt)\n",
    "    c = moment2 * sum(ytm1) - moment1 * sum(yt + ytm1) + sum(yt)\n",
    "    d = 2 * moment1 * sum(xt * xtm1) - moment2 * sum(xtm1 ** 2) - sum(xt ** 2)\n",
    "    e = moment2 * sum(xtm1 * ytm1) - moment1 * sum(xtm1 * yt + xt * ytm1) + sum(xt * yt)\n",
    "        \n",
    "    # solve simultaneous equations\n",
    "    slope = ((a * e) - (c * b)) / ((b ** 2) - (a * d))\n",
    "    intercept = (-slope * d / b) - (e / b)\n",
    "\n",
    "    # now find optimal sigma\n",
    "    eps = y - intercept - slope * x\n",
    "    ept, eptm1 = eps[1:], eps[:T-1]\n",
    "    std_eta = sqrt( (sum(ept**2) - 2 * moment1 * sum( ept * eptm1)  + moment2 * sum(eptm1 ** 2)) / (T-1) )\n",
    "\n",
    "    assert std_eta>0,\"Standard deviation must be positive\"\n",
    "    assert isreal(std_eta),\"Standard deviation must be real\"\n",
    "    return slope,intercept,std_eta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the real meat of the routine, simple since the inference routines are given above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CointInference(epsilon,std_eta,x,y):\n",
    "    loglik,moment1,moment2 = Filtering(epsilon,std_eta)\n",
    "    slope,intercept,std_eta = EMUpdate(x,y,moment1,moment2)\n",
    "    #std_eta_with_old_regression = sqrt( (sum(epsilon[1:]**2) \\\n",
    "    #                                      - 2 * sum(moment1 * epsilon[1:] * epsilon[:-1]) \\\n",
    "    #                                      + sum(moment2 * epsilon[:-1] ** 2)) / (x.size - 1) )\n",
    "    return loglik,slope,intercept,std_eta#,std_eta_with_old_regression,moment1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, the function we'll expose to check for cointegration using the Bayesian method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BayesianLearningTest(x,y):\n",
    "\n",
    "    slope, intercept, std_eta_coint = LinearRegression(x,y)\n",
    "\n",
    "    # cointegrated case - learn slope, intercept, std by ML\n",
    "    logliks = [-inf]\n",
    "    for i in range(1000):\n",
    "        assert ~isnan(intercept), \"Intercept cannot be nan\"\n",
    "        assert ~isnan(slope), \"Slope cannot be nan\"\n",
    "        assert isreal(std_eta_coint), \"Standard deviation must be real\"\n",
    "        assert std_eta_coint > 0, \"Standard deviation must be greater then 0\"\n",
    "\n",
    "        loglik_coint,slope,intercept,std_eta_coint = CointInference(y-intercept-slope*x,std_eta_coint,x,y)\n",
    "        \n",
    "        if loglik_coint-logliks[-1]<0.00001: break\n",
    "        logliks = logliks + [loglik_coint]\n",
    "    \n",
    "    # non-cointegrated case - use above slope, intercept, use ML std\n",
    "    epsilon = y-intercept-slope*x\n",
    "    std_eta_p1 = sqrt(((epsilon[1:]-epsilon[:-1])**2).mean())\n",
    "    loglik_p1 = sum(norm.logpdf(epsilon[1:], loc=epsilon[:-1], scale=std_eta_p1))\n",
    "\n",
    "    bayes_factor = exp(loglik_p1 - loglik_coint)\n",
    "    cointegrated = loglik_coint > loglik_p1\n",
    "    \n",
    "    #print \"slope,intercept\", slope,intercept\n",
    "    return loglik_p1 - loglik_coint#cointegrated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bringing it all together\n",
    "Here we'll test the routine. First, generate some data to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cointegrated? False\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VGXWB/Dfk9577wHSCKEnCoIGVMSG4irqKiq66i5Y\nVncX276vuvpiW/uKbS2oWBbFpahggSgoEHoSCCmQhPQyyaT3nPePm0lmMiWTKZmS8/188iFz6zOX\nyZnnnqdcQURgjDE2sThYugCMMcbGHwd/xhibgDj4M8bYBMTBnzHGJiAO/owxNgFx8GeMsQnI6OAv\nhIgSQuwSQpwQQuQKIe4bXO4vhPheCFEghNgphPA1vriMMcZMQRjbz18IEQYgjIiOCSG8ABwGcBWA\nVQBkRPS8EOIhAP5E9LDRJWaMMWY0o2v+RFRDRMcGf28DkA8gCtIXwIbBzTYAuNrYczHGGDMNo2v+\nKgcTIg5AFoBpAMqJyF9pXSMRBZjsZIwxxgxmsgbfwZTPlwDuH7wDGPmtwvNIMMaYlXAyxUGEEE6Q\nAv/HRLRlcHGtECKUiGoH2wXqtOzLXwqMMWYAIhKG7muqmv/7AE4S0atKy7YCuG3w91sBbBm5kwIR\n8Q8RHn/8cYuXwVp++FrwteBrofvHWEbX/IUQ5wG4CUCuEOIopPTOowCeA/AfIcTtAMoArDD2XIwx\nxkzD6OBPRL8CcNSy+iJjj88YY8z0eISvFcnMzLR0EawGX4thfC2G8bUwHZN29TSoAEKQpcvAGGO2\nRggBsoIGX8YYYzaEgz9jjE1AHPwZY2wC4uDPGGMTEAd/xhibgDj4M8bYBMTBnzHGJiAO/owxNgFx\n8GeMsQmIgz9jjE1AHPwZY2wC4uDPbNbByoO4efPNli4GYzbJJE/yYswSnt7zNLIrsy1dDMZsEtf8\nmU06UXcCByoOoLGzEZ29nZYuDmM2h4M/s0nP//Y87jvnPkT7RONs81lLF4cxm8PBn9mcMnkZthdu\nx+r01Yj1i0VZc5mli8SYzeHgz2zOS/tewh2z7oCfmx/ifONQKi+1dJEYsznc4MtsSkNHAz7O+Rh5\nq/MAAHF+cSiTc82fsbHimj+zKe8cfgfXpFyDCO8IAECsXyxKm0stWyjGbBAHf2ZTsiuzccnkS4Ze\nc82fMcNw8Gc2pUBWgKSgpKHXsb6xnPNnzAAc/JnN6BvoQ0lTCRICEoaWRfpEoq69Dj39PRYsGWO2\nh4M/sxklTSUI9w6Hu7P70DInBydEeEegvLncgiVjzPZw8Gc2o0BWgKTAJLXlcX7c3ZOxseLgz6xO\nmbwML+97WW15QYPm4M8DvRgbOw7+zOr8Vv4bnvv1ORCRyvKRjb0KPNCLsbHj4M+sTqm8FLXttaho\nqVBZrivtwzV/xsaGgz+zOopAfqDygMryggbNNf9YP+7uydhYcfBnVqdUXop5UfNU5uqXd8nR3tuO\nSO9Ite15oBdjY8fBn1mdUnkpVqSuUAn+BQ0FSAxMhBBCbfsonyhUtVahb6BvPIvJmE3j4M+sChGh\nrLkM1069FoerD6N/oB+A9nw/ALg4uiDUK1StjYAxph0Hf2ZVattr4eXihSifKER4RyC/IR+A9m6e\nCrG+sZz6YWwMOPgzq1ImL0OcXxwAICMyYyj1o62bpwIP9GJsbDj4M6tSKi8dDv4RI4K/jpr/yO6e\nObU5OFF3wqxlZcyWmST4CyHeE0LUCiFylJb5CyG+F0IUCCF2CiF8TXEuZt9K5aWI840DMFzz7x/o\nx+nG00gMTNS6n/Lsno2djbj808tx6cZLUd9ePw6lZsz2mKrm/wGAS0YsexjAj0SUBGAXgEdMdC5m\nx0rlpYj1iwUAzAibgQJZAQpkBQjyCIKni6fW/RQ1fyLCXdvuwrUp1+KmtJtw89c3DzUaM8aGmST4\nE9FeAE0jFl8FYMPg7xsAXG2KczH7VtY8nPN3c3LD1OCp+Cz3M535fmB4oNeHxz5EoawQz1z0DJ5a\n/BS6+7rx9C9Pj0PJGbMt5sz5hxBRLQAQUQ2AEDOei9kJ5Zw/IOX9P8n9RGe+HwBifGNQ3lyOtT+u\nxae/+xRuTm5wcnDC59d+jneOvIPvT39v5pIzZlvG8wHupG3FE088MfR7ZmYmMjMzx6E4zNoQkZT2\n8Y0dWpYRmYH1h9aPGvzdnNwQ6hWKtfPXYlrItKHlYV5h+OCqD3D/jvuRvybfbGVnzNyysrKQlZVl\nsuOJkTMnGnwgIWIBbCOi6YOv8wFkElGtECIMwG4iStGwH5mqDMy21bXXIeWNFMjWyoaW5dfnY+r6\nqdh5804smbxE5/7VrdUI8wpTGwXcN9AHn2d8UP+3ep3tBozZEiEEiEh9yLueTJn2EYM/ClsB3Db4\n+60AtpjwXMwOjUz5AEBSUBKCPYIxNXjqqPuHe4drnP7BycEJiYGJOFl/0lRFZczmmaqr56cAfgOQ\nKIQ4K4RYBeBZABcLIQoAXDj4mjGtlAd4KTgIB5y5/wyifKKMOnZaaBpy63KNOgZj9sQkOX8i+r2W\nVReZ4vhsYlDu46/My8XL6GOnhaQhry7P6OMwZi94hC+zGprSPqaSFsI1f8aUcfBnVqO0eXiAl6ml\nhaYht5aDP2MKHPyZ1dCU8zeVSO9IdPd383QPjA3i4M+sgqY+/qYkhODUD2NKOPgzqyDrlMHF0QW+\nbuab/y8thFM/jClw8GdWQXlCN3OZFjKNa/6MDeLgz6yCOfP9CtzXn7FhHPyZVdDWx9+UpoVMw4m6\nExigAbOehzFbwMGfWYUzTWfMnvbxc/NDgHsASppKhpb9evZXZJVmmfW8jFkjDv7MKvxc9jPmR883\n+3mUUz/tPe248asb8dfv/2r28zJmbTj4M4srby5HbXst5oTPMfu5lHv8rNuzDvOi56G2vXZMvYC6\n+rrw2E+PmauIjI0LDv7M4nYU78CSyUvg6OBo9nMp+voXyYrw9uG38dKSl3DL9FvwwbEP9D5GQUMB\n1u1dh7PNZ81YUsbMi4M/s7jvir/D0slLx+VcirTP/Tvux0PnPYRIn0ismrUKG3M3ore/V69jFDcW\nAwB2lewyZ1EZMysO/syievt7satkFy6Zcsm4nC8pMAnFjcUokZfg/nPvBwBMCZiCpMAkfFP0jV7H\nKG4sRqB7IAd/ZtM4+DOL+q38N0wJmIIQz/F5xLOrkysumnQR3rjsDbg4ugwtXzVzld6pn+LGYtw2\n8zbsKtkFfgods1Uc/JlJnWk6g8d+egxf5H0BWYds1O13FO/ApVMuHYeSDfvupu+wOH6xyrLrUq/D\nL2W/oKatZtT9ixqLsHTKUjgIBxQ1FpmrmIyZFQd/ZhLyLjn+9v3fkP5uOtp72/FJ7ieIfzUeGe9m\n4FTDKa37fVf8HZZOGZ98vy5eLl64OvlqfJLzyajbFjcWIyEgAYvjF3Pqh9ksDv7MaNmV2Uj6VxLk\nXXLk/SkPryx9Bdtu3Ib6v9UjMTARO4t3atyvqrUKZ5vP4pyoc8a5xJqtmrkKH+d8rHObjt4ONHQ0\nIMonioM/s2kc/JnR1u1Zh8cveBzvLnsX4d7hQ8tdnVyRFpKG8pZyjfvtLN6JiyZdBCcHkzxN1Gjz\no+ejuLEYLd0tWrc503QG8f7xcHRwxKK4Rdhdupuni2A2iYM/M0qpvBR7zu7BrTNu1bg+2jdaa/Df\ncXr88/26ODk4YXrodBytPqp1m+LGYkwJmAJAem/+bv48TTSzSRz8mVHeOvQWbp1xKzxdPDWuj/aJ\nRnmzevAnIvx45kcsmbzE3EUck7nhc3G4+rDW9cWNxZjiP2Xo9YXxF3Lqh9kkDv7MYJ29nXjv6HtY\nnb5a6zbRvtGoaKlQW97Y2Yj+gX5E+kSas4hjNidiDg5VHdK6vkhWhITAhKHXi+MXY1cpB39mezj4\nM4N9ceILpEekD6VBNInwjkBNWw36B/pVlpfISxDvH2/uIo7Z3IhRav5NxSrvNzMuE3vK9qBvoG88\niseYyXDwZwYhIrye/TruybhH53Yuji4I9AhEdVu1yvKSphJM8p9kziIaJDkoGZUtlWjuata4Xjnn\nDwDBnsGI9YvV2qOJMWvFwZ8Z5EDlAci75Hr10deU9y+RlyDez/pq/k4OTpgRNgNHa9Qbfbv6ulDT\nVoMY3xiV5S8ueRG3b70d+fX541VMxoxmHX3smNXrG+jD6m9Wo669Du297SiSFeG+c+6Dgxi9/qDo\n8TMP84aWlTSVYFrINHMW2WBzwqW8f2ZcpsrykqYSxPrGqnVNvWjSRXj+oudx2aeXYd8d+xDmFTaO\npWXMMBz8mV5ONZzCztM78erSV+Hh7AEvFy+cG3WuXvtqq/lfmXSlOYpqtDnhc7Dj9A615cWNxSqN\nvcpunXkryprLcPmnl+Pn236Gl4uXuYvJmFE47cP0crT6KOZFzcPVyVdjyeQlmB89X69aPzAY/Fts\nI+0DDDb6Vqk3+o7s5jnS/5z/P5gWMg1/3/V3cxaPMZPg4M/0crTmKGaGzTRo35EDvQZoAGXyMsT5\nxZmodKaVHJSMqtYqtUbfosYinT2bhBBYNXOVzt5CjFkLDv5ML8dqjmFW2CyD9h2Z9qlqrYK/uz/c\nnd1NVTyTcnRwxMywmThSfURl+ciePpokByXrnMiOMWvBwZ+NiohwrOaYUTV/5YFe1trNU5mi0VeZ\nPsE/1DMUPf09ek1nzZglcfBnozrbfBZuTm4I9Qo1aP9wr3A0dDSgp78HgHXn+xVGDvbq6e9BZWvl\nqKkqIQSSApNQICswcwkZMw4HfzYqY/L9gJRGCfMKQ1VrFQCp5m/twX9OxByV4F8qL0W0TzScHZ1H\n3Tc5KBkFDRz8TU0mAz77DHjqKeC224ArrgB+/ll9m+uuA8LCgH/9C+jpsUhRbQIHfzPq7wfWrwdu\nvRU4fdrSpTGcMfl+hSifqKG8v7VO7aAsKTAJNW01KJIVobW7FYWywlFTPsr7TsS8f3s78PrrwCef\nAJWV+u9XVgZs2wZs3gz85z/AV18BpaWA4gmZzc3AE08AiYnAF18AHR3AwoXANdcAN98s/X3V1QHf\nfw/MmAFERwNbtwLbtwOpqdIxu7tVz9nTA3z5JfDyy+rr9NHbC3z8MTB9OjBpEvDnPwNZWUCfDc3y\nwf38zeTAAWD1asDLC1i0CDjnHODee4GHHgLc3CxdurE5WnMUN6fdbNQxlHv8nGk6o3UKaGvh6OCI\na1Kuwfkfno+2nja09bRh7fy1eu2bHJQ86kNh7ElvL/Dee8A//gHMny8F7fvvB4KCgFmzgOBg6XfF\nT2Ag4OsL/PqrFMyLiqS/D1dXwMlJCsZr1ki/n3MOsGcPcOmlQHY2MHmy6rmvu076YkhOBjw9gQ8/\nBC66SFq3Ywfw44/Ak08Cd94JXHABsHQpUFICfPQRMHWqtM/770v7zZmj+33W1wN5eVI51q+Xgv4L\nLwAREcCWLcBf/yp9aV12GbBsGXDJJdJ7qKgAysuBlhbA0XH4RwjpuA4OUvmjokz7/zIaYe4HUAsh\nlgJ4BdJdxntE9NyI9WRPD8EuLZU+jN9/Dzz/PHDTTdJ/8tmzwIMPAkePAq+8It2yKv7zrV3sK7H4\n6Zaf9K75avLX7/+KYI9gPLTgIUS/HI09q/ZYbVdPTQZoAAICQo//tLy6PFy36Trkr9E83cPp01Lw\nkcmkGuMUwy+rxeXkSAE4Ohp47rnhADowABw/Dpw8Kb3PhgYpeMpkwz+zZgHXXw9ceCHgPCKbRiRd\np337gLlzgZQU3eUoKQH8/QE/P83rZTLpb3LnTikldMcdQEKCdJ7PPgMeeAC46y7g4YelLwSF/n7p\n7/WFF6QvpWnTgLQ06W7jHA0PoKuokO5itmwBfvlFug6RkdL18fWVjqf4UbzPvj7pS8XdHZg3T/ri\nuP566YtDFyEEiMjwKEJEZvuBFPCLAcQCcAZwDEDyiG3IHlRWEq1eTRQQQPT3vxPJ5Zq327GDKCmJ\n6NJLiQoKxreMhmhobyDvdd7UP9Bv1HFe2fcKrflmDXX1dpHLUy7U299rohJan87eTnJ9ypV6+npU\nlm/ZQrRgAVFwMNF99xE99hhRYCDRTTcRnThhocIa4aefpPeycSPRwIClS2Ocqiqi668nCg0lev55\norY2orw8onPOIcrMJDp5cuzvsbtb/30GBogKC4k++IBo4UKiKVOINmwg6u0lyskhevZZovPPJ9qz\nZ3ifwdhpeHw2ZudRDw6cC+A7pdcPA3hoxDZ6X0xdiouJnnuO6NZbic6eNckhR9XaSvTpp0TLlhH5\n+hL95S9EdXWj79fdTfTPf0p/+LNmEaWmEiUkEM2dS/Tqq0QymX7n7+0lev11otOn9dt+YICoqIjo\n3XeJHnqI6P77ie6+WwpEe/dq/qD+dOYnWvD+Av1OoMNXJ7+iZZ8to8KGQop/Jd7o41m7Sa9OooKG\n4W/3jz8miogg+vpr6f9fQS4nWreOKCSE6NpriY4ds0BhDbBxo1Tm3bstXRLTyskhuu466b0FBRG9\n9RZRv3H1njEbGCDatYvogguI3N2J4uOJ1qwh+uYboo6O4e2MDf7mzvlHAlAe118BIGO0nQYGgLfe\nArq6gPvu03378+mn0u1mbS1w9dVS/i09XcrhLVWacLKuTmrk0TevVlsLvPSSdEuYmQnccouUu+/u\nBr79VmpE2rkTWLBAukX76CPptk4fLi7AX/4i3TqePSu9dnaW8oLvvw/87/8Cl18O/OEPUp7SQUOz\nfGMjcMMN0r9PPgm8+CKwcqXmVFJNjZSP3bZNut1ctEhqCAsNldofGhuBVaukMtx5p3Q77ucn3UYf\nqjyKmaGG9/RRiPaJRmlTOfaetP4+/qagaPRNDEzE558Da9dK+eepU1W38/UFHnlE+py//baU205P\nl/6/ZswwTVkaGoCqquHc+8gUy0gDA1KD6ebNgIeHVEZvbylnXVsrHauwENi1S/oc2ZO0NOlvu7BQ\naq+LiBj/Mggh/Y0uWiTFreBgM6WIjfnmGO0HwO8AvKP0+mYAr43YRuVbr7ZWSomccw7RhRcSnXuu\ndDs0klxO9PvfEyUnS7effX3D637+mSgykuiRR6Qa9oIFUs08IIBoxQqigwe1f+vW1hLdey+Rv7/0\nbXvkCNHLLxPNnk0UHi4d5+KLpdpzQ4P24xijoYHolVeIpk0jmjxZqhkeOULU1CStz8uTlj/4oFT7\nP35cuntYsYLozJnhmkp3N9ELL0h3GH/7m5Rm0nYbOjBAlJUl3TktXCidOzKSyOPmm+iPb79n9G39\n5u+ryeGhYHJf+BbNefIPKv9f9uiBHQ/Q83ufp02bpFRCTo5++3V0EL32mlTzvOsu6fOoj/Jy9b+D\ngQGiTz6RUjOpqVI5nJykmuSzz6p/fvv7iTZtIpo+nWjmTKI33yRav57omWekv6VnnyV6/32pBlpf\nr1+5mPnABtI+O5Rea0z73HHH43TbbY/T8uWPU2DgbnrsMaKeHunD+Prr0u3XP/4hfZD/+1+iL74g\niosj+uMfidrbNV+Ymhoph3fnnUTffkvU1UXU0kL00ktE0dFEixYR7d+vus+PP0q35vffT1RdrX7M\nggL9/xhNYWCA6MAB6T2kpRF5eRH5+UlfTBs2qG7b0UH05z9L5Xdzk74U4+KMb1uIez6Vpiw8PHS9\njhyRrucHHxB99RXR4cNSENH25dDVRbR2LVFYeD85P+lKt2+6j+Ju+T9atIioosLwco2H2lopF5+X\nR5SbS9TcrP++z/zwFiX+7Q4KDSU6enTs525qInrgAemL++mnNacCBwaIfviBaPly6TMxfbr0f/7c\nc0TPffcJXX1ND6WmEh06NLxPf7/0+rbbpM/SzTdLlag5c6TPV3o60bZttp/Dt0e7d++mxx9/fOjH\n2OBv1t4+QghHAAUALgRQDSAbwI1ElK+0Dc2dS3B1lVrZH35Yut1RVlgotbjL5UBbm5QOWrMGuOoq\nw8rV2yulaf73f6W0zdNPS312//1vYMMG4OKLDXzDZkYkpWj6+qSUjTYdHcCZM0Bnp5RCMFRnbycC\nng9Aw1/k+PA9V7z8snQrHBYGhIRIaYDSUqmnRXe3lB7w85O2aW+X1svlUpe3d94Bzvl0Evzd/fHg\nOX/Fma034tVXpes/d65UzgsusI5usKdPA+vWAV9/Lb1XxS13Y6P0OVyxQv02vK1NukWvrZXSJW/v\n/Bm+yx9Dzp/3wt/f8LIUFEhl2boVuP7GPsgy7sdlzs9jz0+e+OEHKTW3erXUq8zbGzh4EHj9jV58\nHOOB23oP4q0nZsLVVfOx6+qkFIePj9SHPikJRpWVjS9je/uMV1fPVzHc1fPZEeupvr0eQR5BZi2H\nJu3tUl7/6V0vYpb3Uvz3nVSE8XM4hhysPIg/bPsDjv/x+KjbdnRIg3Gam6VA6OkpfRn4+EhfBgBw\nwYcXYO/Zvfj19l9xbtS5KC+XuvIdOiT15e7tlbrIRZrxme4HDgCvvip9SaWnSz/e3lLAP3NG6p73\n7bdS5eLPf1YNhgcOALffLnXNfPZZIDdX2nbHDul9h4ZKPxkZwB/+XIOLNqeh/m/1Jil3dTXwx/Wf\nYKvTSswq+Rir5tyMiy+WAvbIL6JTDaeQ8kYKNl6zEb9P+71Jzs+sj7HB3+yDvIhoB4AkXds8sPMB\nfLx8/AfFeHoC//M/wAbft3DrPC+EhdlZ65WRdp7eiYyIUdvnAUgNgx4eQHi49m2ifKIwQANDUztE\nR0s/K1ZIdzXPPiv1nf7vf6W7gbY26Q7tyy+lfthXjvLsl64u6W5HU32moECaFiAvT2psHxiQGmCf\neUb64po8WRq0M3Mm8NprmvuLn3MOcOSIVBOfP1/qk3355cDjjwNxcapBmCgUvf29kHXIEOgROPoF\nHEVo2ACKw5/BqshVqJz8Me69WfugO8XjJE/WnzT6vMx+WcUI371n92JH8Q69ngdras1dzTjdVIyz\nzaXjfm5rVtlSiVf2v4L9f9hvsmNG+0TDw9kDIZ4hauuEkHq9JCdLPV6WLZPuAs4/X+ppde+9UrB+\n7jn11FBfnzStwJNPDh9rZG04MFAK+l9/Da1pEH24ukrnUZxLGyEEkoKkCd7me8w3/ISDthVsg5uT\nG/512b8Q+VIkqlurEe6t+Zv2VMMpxPjGcPBnOlnF3D5vX/E2/rj9j2jraRv3cyvmbC9rLhv3c1uz\ntT+uxd1z7jZqVO9I0T7RiPeL1zlKdvly4IcfgNhYKX+9ebM0idfRo9J8MfPmSXnqQ4ekUZuHDklp\nlu3bpWH3cjnQ1CTl55V/ioqAP/7RuMA/Vqaa44eIsG7vOjy64FF4OHvg6uSr8VneZ1q3z2/IxzXJ\n13DwZzpZRfBfMnkJMuMyLfL4uyPVRzArbBZK5aXjfm5rtadsD/aU7cGjCx816XGnh07HvKh5o243\nc6bUGB+vNPebvz+waZN0B/D559J4hMmTpXlUHnhAuitITDRpcY1mqtk9d5fuRnNXM5anLAcArJy+\nUufcQfkN+bgq+SqUNZcNTaPN2EhWEfwB4KVLXsKnuZ+ipKlkXM97uPowfpfyOw7+g/oG+nDPd/fg\nhYtfgKeL5+g7jMHC2IV4d9m7Bu8vhNTgunmzdCcgl0sD2LQNbrM0U83r/8zeZ/DwgoeHnpmcGZeJ\nho4G5NXlqW1LRDjVcAozQmcg1jcWRbIio8/P7JPVBP8A9wDMCp817lPhHq4+jCsSr4CsU4auvq5x\nPbc5dfR2YGvB1jHv987hdxDgHoAVqSvMUCrT0zT62VqY4pGOx2qOoaChQKXXjoNwwE1pN+Hj4+q1\n/8rWSng6e8Lf3R9Tg6dy6odpZVV/OjE+MTjbfHbcztfS3YKKlgqkhqQiyidqXM9tbi/8+gJu/OpG\nDNDAmPZ75/A7eGrRU3rNXsl0mxwwGWeazsCY7tT7K/bjksmXwMXRRWX5yukrsTF3I/oH+lWW59fn\nIzkoGQA4+DOdrCr4x/rFjmvD67GaY0gLSYOTgxPi/OJQJrePRt/y5nK8lv0anB2cUdmi/1M12nva\nUdRYhPQII0aGsSFuTm7wcPaAvEtu8DGKZEVIDFRvzEgNSUWIZwiySrNUlp9qOIWUIGn+46nBU3Gy\ngYM/08y6gr+v7uBf0VKBF359AQveX4DTjcY/Gutw1WHMCZcmII/zjbObvP8jPz2CP83905jTaEeq\nj2BayDS4Oo1jlxg7F+wZjPoOwwd6FTYWIiEwQeO661Ovx+b8zSrL8hvykRIsBf+UoBSu+TOtrCv4\n+8VqrH239bRhycdLMP3N6SiQFaCrrwvHa0cfdTqaIzVHMCdiMPj72Ufw31+xH7tLd+PhBQ8jOTB5\nTA2O2ZXZeg/qYvoJ8QxBXXudwfsXygo11vwB4MqkK7GtcJtKWim/YTjtkxSUhOLGYvQNaH+2YHZl\nts71zH5ZVfCP8Y3RWPM/Wn0UDR0NqPpLFf697N+YFzVv6HmwxjhcdRizw2cDkL54Sm18oBcR4c87\n/ox1i9fBy8VrzA2O2VXZyIjk4G9KxgT/voE+lMnLMNl/ssb1KUEpcHF0UakIKad9PJw9EOEdofUu\nuamzCQs/WIiPjn9kUPmYbbOq4B/pHYnatlr09veqLC9qLEJaaBrcnKShncrPgzVUW08byprLkBos\nTelgDzX/r/K/Qt9AH1bOWAlAqvmNKfhXcvA3tWCPYNS3G5b2KZWXIsI7QmsaTgiBKxOvxLaCbQAA\neZccbT1tiPIZfmiFrkbf/5z4D+L84vDs3mfVGo6Z/bOq4O/s6IwwrzBUtqo2UhbKCpEYMHzrG+0T\nbXTPnOM1x5EanApnR+nJFvYQ/PeV78OK1BVD/cHHUvOva6+DvEuuNb/MDGNMzb9Qpj3fr7AsaRm2\nFkpdek81nEJyULJKT62pQVOR36D5WcIf5XyEf178TwR7BuPLk18aVEZmu6wq+AOa8/4j854xvjFG\n1/wPVw+nfAAgwjsCDR0N6O7rNuq4llTVVoVI7+EpMWN8Y9DY2YjW7tZR9z1YeRDpEelDXxzMNIwN\n/sqVHk0WxEidH6paq1S6eSpoq/kXNxajuLEYS6csxaMLHsW6veuM6pLKbI/V/aVryvuPDP7RvtFG\n5/wPVw8CPideAAAgAElEQVT39AEAJwcnRHpHGv2lYkmVLZWI9BkO/g7CAYmBiSiUFY6674HKA5zy\nMYNgD8N7++hq7FVwdnTGJVMuwfbC7VJPn8F8v4K24P9Jzie4cdqNcHZ0xmUJl0FA4JuibwwqJ7NN\nVhf8Y31Va/79A/043XRaZYKxCO8I1LXXqbUNjMWR6iMqNX/A9lM/la2ViPBWfeiovqkfzvebhzE1\n/6LGIr3ScMsSl2Fb4TaVxl6F5CCpx5dyTp+I8NHxj7ByutQ2JITAowsfxf/t+T+u/U8gVhn8lfP5\n5S3lCPIIUplnxsnBCaFeoahqrTLoHJ29nTjdeBrTQqapntsv1maDPxGhqlU17QPoF/yJCNmV2Ty4\nywyMTvuMUvMHgKVTluLn0p9xtOaoWtrH29Ubge6BKnfTv5b/Cndnd5XKz+9SfgdZh0xt0BizX9YX\n/EeM8tX2BxDtY3iPnxP1J5AQmKDWi8KWB3o1dTXB1dFVbTK2pMAknJLpDv6nm07Dy8VL6/zwzHCG\nDvLq7O1EbVstYn1jR93W390fcyPmoqatRuMU3CNTPx8d/wi3TL9FpWHY0cERq2auwrbCbWMuK7NN\nVvEwF2Ujc/7aGr2ifQ3v8ZNbm4u0kDS15XF+cfjhzA8GHdPSKlvUUz6AfjV/TvmYT5BHEBo7G9E/\n0A9HB0e99ytuLMYk/0l673Nl4pWobqse6r2mbFHcIqz8eiUWxizE+bHn46v8rzQ+mjMlOAW/Hf1N\n7zIy22Z1wV+R9iEiCCG0zm0S4xNjcKNvTm0OpodOV1se5xdnsw91qWqtUmnsVUgMTERxY7HO4MPB\n33ycHJzg6+qLxs5GBHsG672fvvl+hZUzViLWT/NdwkMLHsItM27BnrN78EvZL1g1c5XKWAAFfTsH\nMPtgdcHf08UTns6eqGuvQ6hXKAobC7Fk8hK17aJ9o0f9oD6++3HcMO2GoblOFHLrcjUe01obfLcX\nbscnOZ/Aw9kDHs4emB0+G7fPul1lm8rWSrV8PyBdzxDPEJQ1l2GS/ySNx8+uzMa6C9eZpexsOPUz\nluCvTzdPZUEeQbgm5Rqt68O9w7EidYXOqbon+09GmbwMvf29Gu8gmH2xupw/IOX9FSkdbQNd9Bno\n9cWJLzR2X9NW84/0kUYYW9vTj/752z+REJCA86LPQ7RPNB7+8WG1bbSlfQDdqZ+e/h4crz2OuRFz\nTVpmNsyQRl99G3tNydXJFeHe4VZZAWKmZ5XBX5H37+7rRmVLJeL94jVuo6vBl4hwtvksfitXzWHW\nttWib6BPY6B0cnBChHeESeYNMhVZhwxHa47i0YWP4o7Zd2DteWvR0duB5q5mle201fwB3c+S3XBs\nA+ZHz4eXi5fJy84khgT/okbN6U5z49TPxGGVwV/R1/9M0xnE+MZovAUdbaBXfUc9+gb68Fv5byp9\nl3PrcjE9dLrWh5VYW+rn26JvsTh+Mdyd3QFIfbITAhNQ1Kj6eD5tOX9Ae82/u68bT+95Gv/I/Ifp\nC86GGDK/jz5TO5hDYgAH/4nCeoN/c5nOW99gj2C09bSho7dD4/pSeSnSQtPgIBxQIh9+LnBObY7G\nnj4K1tbou7VwK5YlLlNZlhCQoPZsVk0DvBS0Bf9/H/k3poVMw7zo0R+qzgw31pq/vEuO9p52hHuN\nf9dbrvlPHFYZ/BVpH13BXwihs/ZfJi9DnF8c5kfPx77yfUPLFTV/bayp5t/d140fTv+AKxKvUFme\nEJCg9gda2aI97aMY5amss7cT6/au41r/OBhr8Ff0cLPEozQTAxPV7iqZfbLK4K9o8B2t0UvXQK9S\neSlifWMxP3q+St4/pzYHaaG6a/7KdwqWtLt0N6aFTFPrJTLyD7S3vxeyThlCvUI1HifcKxydvZ1o\n7GwcWvbmoTeREZkx9DAbZj5jnd9nrN08TYlr/hOHdQb/wZz/aI1eugZ6lTUP1/x/q5CCf99AH/Lr\n89WmdVCWGJiIggb9n35lTlsLtmJZ0jK15SNz/jVtNQj2CIaTg+aeu0IIZERm4LKNl+GFX1/A0eqj\neP7X5/Fk5pNmKzsbNtaa/1i7eZpSjG8M6jvqtaZTmf2wyuAf5BGErr4uHKs5NnrNX1vap7kMsb6x\nmBU2C4WyQrR2t6K4sRjh3uE6e7Yo8uOWnuCKiLC1YCuuSrpKbZ0i7aMoY2VrpdbGXoXtv9+Oxy94\nHGeazuCKz67AhZMu1Jn+YqYz1uBfICuwSE8fQJrmYZL/JBQ3Flvk/Gz8WGXwF0IgxjcGvQO9Whsx\nAd3dPUvlpYj1i4Wrkytmhc1CdmU2cmt15/sBIMA9AO7O7gZPGmcqR6qPwNPFE0lBSWrrgjyCQESQ\ndcoAQOOEbiO5Obnh0oRL8eYVb6LigQp8svwTs5SbqRvr/D4n6k4gNSTVjCXSTVfqZ4AGsDl/M7Ir\ns8e5VMzUrDL4A1LePyEgQefDRbQN9CKioQZfAEN5/9F6+iikBKVoffrReNlSsEVjrR+QvhwTAxOH\nevzoauzVtr8lGhMnqgD3ALR0t+g1BXlvfy+KGovUZuccT5q6exIRfjj9AzLezcCqLav4ub92wHqD\nv2/sqLe+2p7lK++SQwgBPzc/AFLw31exb9SePgopQSnIr7ds8N9euB1XJl6pdb1y3l9XN09meQ7C\nAYHugWjoaBh12+LGYkT5RMHD2WMcSqZZQqB6b7KbNt+ENd+uwdrz1uKty98yeJpqZj2sNvinR6Tj\nvOjzdG4T4ytN7jYyP6/o6aMwL2oe9lXsw/Ha43rV/JODki1a8+/u60Z+Qz7SI7XPr6/c11+fnD+z\nLH1TPyfqTyA12HIpH0A97XO85jh+KfsFeavzsCJ1BUK9Qjn42wGrm9hN4c45d466jY+rDxwdHNHU\n1YQA94Ch5YqePgqhXqEIcA9AVWuVxvnOR0oJTsGWgi0GldsUCmWFiPOLg5uTm9ZtEgMTh8qoT86f\nWZa+jb55dXk6e6ONh5HB/42Db+DuOXfDxdEFgHEPqGHWw2pr/vrS1ONnZM0fkFI/qcGpes2Pbumc\nvz4BQKXmr2NSN2Yd9A2Y1lDzD/UMRU9/Dxo7GyHvkmPTyU0qlbEQzxCDn0vMrIfNB39NPX7K5GVq\nc5tfPOliLIhZoNcxo3yi0NbTBnmX3GTlHIu8ujxMCx4l+A/m/ImI0z42QN/5fayh5q/coeDDYx9i\n6ZSlCPMKG1of6B4IeZccfQN9FiwlM5bNB39NPX5Gpn0A4JYZt+CVpa/odUwhhM6ZMM0tr370AODn\n5gc3JzcUNxZjgAbg6+o7TqVjhtCn5t/d142SphKL9fFXlhiYiFMNp7D+4Hrck36PyjpHB0f4uflB\n1iGzUOmYKRgV/IUQ1woh8oQQ/UKI2SPWPSKEKBJC5Ash1J+cYiIxvjEok6tOxKYp7TNWKcGW6/Gj\nb+0vMTARWaVZiPCO4K6bVk6f4F8gK8Ak/0lqz5a2hMTARLx56E14OHtgfvR8tfWc97d9xtb8cwEs\nB/Cz8kIhRAqAFQBSAFwKYL0wU3RaELNA7bm7Zc3qaZ+xslTev72nHdWt1ZgcMHnUbRMCEpBVlsWN\nvTZAn/l98uryLDq4S1liYCIOVB7AmvQ1GisWHPxtn1HBn4gKiKgIwMhPx1UAPieiPiIqBVAEwCwP\niV0QswBVrVU43XgaANDa3YrO3k4Ee+j/yDxN9HnwuTmcrD+JpKAkrfP0KEsISEBWaRbn+22APsHy\nRN2JUdt6xktyUDL83Pzw+7Tfa1zPjb62z1w5/0gAyq2wlYPLTM7RwRHLk5fjq/yvAAzX+o290dBU\n8//y5Jdm/0IYS4NfQmACqlqrEOHFPX2snT7BMq/eemr+s8JmIfdPufB08dS4PtgjmGv+Nm7U4C+E\n+EEIkaP0kzv4r/bhp+Ps2qnX4suTXwIY7OljZL4fAKYETEF5czm6+roAAI2djbh9y+145KdHjD62\nLvr09FFQNAxyzd/6BXuOHixP1J2weE8fBSEEonyitK7ntI/tGzW3QEQXG3DcSgDRSq+jBpdp9MQT\nTwz9npmZiczMzDGd7IK4C1AiL0GZvAyl8lK1nj6GcHZ0Rrx/PIobizEtZBpeO/AaLk+8HLtKduFU\nwymzzb2SV5+H++Lv02tbxYA1zvlbP19XX3T2dqK7r1tjg25HbwcqWyv1GoRoDUI8Q3Ck+oilizGh\nZGVlISsry2THM+UIX+U8y1YAG4UQL0NK90wBoHUaQOXgbwgnBydcnXQ1vsr/CjVtNSap+QOD0zzU\n5yPGNwZvHHwD++7Yh5SgFLzw6wt476r3THKOkcaS9vFy8UK4VzgP8LIBQoih1I+mGnV+fT4SAxP1\nauuxBpzzH38jK8ZPPmnc8ziM7ep5tRCiHMC5ALYLIb4DACI6CeA/AE4C+BbAajLzBPmK1I8pevoo\nKPL+/8r+F5ZOWYopAVOwJn0Nvj71NSpbtN7IGKyxsxGt3a2I8Y3Re58PrvoAcyPmmrwszPR0pX7y\n6vIsPrJ3LDjnb/uM7e3zXyKKJiJ3IgonokuV1j1DRFOIKIWIvje+qLotjl+MAlkB9lfsN0naB5CC\n/6GqQ3j1wKt4bOFjAIBAj0CsnL4Sr+zXb8DYWCjmcR9LY/UlUy6xin7hbHQhniGobavVuO5EvfXk\n+/XBOX/bZ/MjfBWcHZ2xLGkZzjafNVnaJyU4BdsKtyEzLlMlx//gvAfx3tH3TD79w1gae5ntSQxI\nxMn6kxrX2VrNn4O/7bOb4A8A16ZcC2cHZ4R7h5vkeEmBSXB2cB6q9SvE+sXi8sTL8ebBN01yHgVb\nq/2xsUmPTMeh6kMa19na/72fm99QAzazTXYV/C+efDHevfJdnU//GgtvV29UPlip8QEwf533V7xx\n8A2TTm5lDZN6MfOZGzEXh6rUg399ez3kXXKTpSvHgxACQR5B3Ohrw+wq+Ls4uuDWmbea9JjBnppH\nCs8Im4F4/3hsOWWaef+JiIO/nUsKTEJtWy2aOptUlu89uxfzo+frNd24NeHUj22zq+A/3takr8G/\nDv7LJMeqba8d6g7I7JOjgyNmhc/C4erDKsv3nN2DhTELLVQqw3Hwt20c/I1wTco1ONVwCifqThh1\nHHmXHCu/XokrE6/k2Tnt3NzwuThYeVBlmS0Hf32eUcCsEwd/I7g4uuDO2Xdi/cH1Bh/jdONpzHtv\nHlKDU/Hule+asHTMGo1s9G3raUN+ve7nNVsr7utv2zj4G+nuOXfjs7zP0NLdMuZ9D1YexHnvn4d7\nM+7FK0tfsbmcLxu7kY2++8r3YVb4LJ3Pa7ZWnPaxbRz8jRTpE4nF8Yvx8fGPx7zvKwdewSMLHsHq\n9NVmKBmzRpP9J6Olu2UoaNpqygcYDP4dHPxtFQd/E1iTvgZvHHwDY53BIqc2BwtjbfMPnxlGCKFS\n+7f14M85f9vFwd8EMuMy0dLdgtNNp/Xep6e/B8WNxUgJSjFjyZg1UjT69vT34GDlQY2PSbQF+kxT\nzawXB38TEEIgIzJDrReHLgUNBYj1jYW7s7sZS8askaLR93DVYSQEJsDXzdfSRTII5/xtGwd/E8mI\nzMDBKv2Df25dLtJC08xYImat5kZINX9bTvkAHPxtHQd/E0mPSB9b8K/NRVoIB/+JKNonGgM0gM/y\nPrPp4O/pLD3isb2n3cIlYYbg4G8icyLm4FjNMb3n+smt4+A/UQkhkB6ZjmM1x2y6wV8IwXl/G8bB\n30T83PwQ4R2B/Pr80TcGp30murnhc5EQkIAwrzBLF8UonPqxXbbxzDgbkRGZgezK7FGDenNXM2Qd\nMkzynzROJWPW5qrkq7ROGmhLOPjbLq75m5C+ef+8ujxMDZ5qsqmnme2ZHT4b92TcY+liGI2f5Wu7\nOPqYkL7BP7cuV+MzAhizNTy/j+3i4G9CM8NmIr8+H119XTq3454+zF5w2sd2cfA3IXdndyQHJeN4\nzXGd2+XU5XBjL7MLHPzNh4jwS9kvZjs+B38TS49IR3Zlttb1RMQ1f2Y3OPibT3FjMS748AIcqzlm\nluNz8Dex9Ejdef+Klgq4ObnZRU8PxoI9gu2ywbdMXmbxL7VCWSEchSOe2fuMWY7Pwd/ERmv05f79\nzJ6Ee4ejsqXS0sUwub/v/jvePvS2RctQKCvEyhkrsatkFwplhSY/Pgd/E0sNSUV5c7nWh7twyofZ\nk3CvcHT1daGxs9HSRTGpQ1WHUNFSYdEyFMoKMTtsNtakr8Hzvz5v8uNz8DcxJwcnzAybicNVhzWu\n52kdmD0RQiA5KFnvke22oKW7BQUNBahotWzwL2osQmJgIu7NuBeb8zeb/MuIg78ZzAmfgyPVRzSu\n47QPszcpwSnIb7Cf4H+k+gh8XH1Q3lxu0XIUygqRGJiIQI9ArJq5Ci/+9qJJj8/B3wxSQ1Jxov6E\n2vLe/l4UyYowNXiqBUrFmHmkBKXgVMMpSxfDZA5WHsTliZdbNO3T0duB+o56xPjGAAAenPcgNhzf\ngIaOBpOdg4O/GUwLmYa8ujy15QWyAkT7RsPD2cMCpWLMPFKC7Kvmf6j6EJZOXorOvk6LTVdd3FiM\nSf6T4OjgCEB6VvjDCx5GdWu1yc7Bwd8MpgZPRX5DPgZoQGU5N/Yye5QSnGJXOf+DlQeRHpmOKJ8o\nVLZapieTIuWjbO15a02aMubgbwZ+bn7wdfXF2eazKst5Th9mjyb5T0JVaxU6ezv13qeztxNEZMZS\nGUbWIYOsU4bEwERE+URZLPVTKCtEYkDi6BsagYO/mWhK/eTU5nDNn9kdJwcnTA6YjAJZgdZtimRF\nWH9wPW79762Y+sZUeK7zxH9P/XccS6mfw9WHMTt8NhyEg+WDfyAHf5uUGpyKE3Wqjb7c04fZK02N\nvhUtFbjvu/sw5bUpyNyQiUNVh7AgegE+/d2nWH/5emzM3Wih0mp3sPIg5obPBQBEedt38OeHuZhJ\nakgqskqzhl7zA1yYPUsJUs/7P7f3Ocg6Zdh8/WakhaRBCDG0LsY3Bg/9+BBau1vh7eo93sXV6lD1\nIdyQegMAIMonSmPHjfHANX8blhqcqvLBya3LRWpIKj/AhdmlkX39iQjbi7bj0YWPYnrodJXADwAB\n7gE4L/o8bCvcNt5F1UnR2AtIwd8SA71kHTL0DvQixDPErOfhSGQmU4On4lTDKfQP9APgnj7Mvo3s\n7nmy/iSICKnBqVr3uT71enxx4ovxKJ5eqlur0dnXiXi/eABAtG+0RdI+ipG9I78wTc2o4C+EeF4I\nkS+EOCaE+EoI4aO07hEhRNHg+iXGF9W2eLt6I8QzBGeazgDgnj7MviUGJqK4sRh9A30AgO2F23FF\n4hU6A9jVyVdjd8luyLvk41VMnQ5VHcLciLlDZbZUg+94pHwA42v+3wNIJaKZAIoAPAIAQoipAFYA\nSAFwKYD1wtxfY1ZIeaQv9/Rh9szTxROhnqEolZcCALYXScFfF183XyyKX4Qtp7aMQwlHd6jq0FBj\nLwAEeQShtbt11CfzafN1/teQdcjGvN94dPMEjAz+RPQj0dBIpv0AogZ/XwbgcyLqI6JSSF8MGcac\nyxYpevwQEff0YXZPMdhL1iFDTm0OMuMyR91ntNRPeXO5wcF3rA5WDef7AcBBOCDCO8KgKauJCHdu\nuxO3b719zOMZbKXmr+x2AN8O/h4JQHlWpMrBZRPKtJBpyKvPw9nms/B09kSQR5Cli8SY2Sjy/juK\nd2BR3CK4ObmNus+ViVfi1/JftdaQr/7ianya+6mpi6rR4erDmBM+R2WZoamf6jZpGoaq1iqsP7h+\nTPuOV/AftaunEOIHAKHKiwAQgMeIaNvgNo8B6CWizwwpxBNPPDH0e2ZmJjIzMw05jNVJDU7FP3/7\nJ9f62YSQEpSC3yp+w9Gao6OmfBS8Xb1x8aSL8fWpr/GH2X9QWXe68TSOVB9BdmU2bp91uzmKPKSu\nvQ49/T2I8olSWR7lE4XylrHP7plTm4MZYTPw1uVvYf7787EwdqFebX4DNICixiIkBCaorcvKykJW\nVtaYy6LNqMGfiC7WtV4IcRuAywAsVlpcCSBa6XXU4DKNlIO/PUkJTkFRYxGOVB/B9BBu7GX2LSU4\nBW8ffhtnms7gxSX6Tz/8+7Tf4+X9L6sF/00nN2F2+GydT8ZT+DT3U6SFpBlcycqry8O0kGlqDdSG\n1vxzanMwPWQ6EgIT8NKSl3DDlzfg0F2HRp3Usaq1Cj6uPvBx9VFbN7Ji/OSTT465XMqM7e2zFMDf\nACwjom6lVVsB3CCEcBFCxAOYAkD7U83tlIezByK9I/H1qa+55s/sXkpQCg5XH8Yk/0mI8I7Qe78r\nE69ESVMJcmpzVJZvOrkJTy96Gvn1+aPOG/Tivhfx45kfRz3XJzmf4FDVIbXleXV5mBY8TW15tI9h\n3T0VNX8AWDljJdIj0zHzrZn40/Y/4bPcz1DTVqNxv/FK+QDG5/xfB+AF4AchxBEhxHoAIKKTAP4D\n4CSkdoDVZI2zOI2D1JBUHKs5xj19mN0L9AhEkEeQ3ikfBWdHZ9w15y6V3PiZpjOoaKnAkslLkBKc\ngmM1x7Tu39HbgeM1x9UmUhypoqUCd227Cxtz1KeVyKvL01hBM6rmr5TmeX/Z+/jsd58hMTARm05u\nwrT107C1YKvKPgM0gC/yvkBKUMqYz2cIY3v7JBBRLBHNHvxZrbTuGSKaQkQpRPS98UW1TanBqXAU\njkgJHp//UMYs6YbUG3B96vVj3u/O2XfiixNfoLmrGQCw6cQmXJN8DRwdHJERkaEz9XO46jD6qX/U\n3PwjPz2CuRFzkV2lnoRQpH1GMiT49/T3oKhR9aFNjg6OmBMxBw/MewCbr9+M7276DndvvxvvHXkP\nANDa3YrlXyxHfkM+/rHoH2M6n6F4hK+ZpQanIjEwUa+eD4zZutcve92gik64dziWTF6CDcc3AJBS\nPtelXgcASI9MR3al9qzx/or9SI9I1xn8syuz8dOZn/D5tZ/jeM3xocFogNQtM68uT+NoZEOCf359\nPuL94nX+zadHpuPn237G/+35Pzz0w0OY//58hHqG4sdbfjT7tA4KHPzN7PLEy/Hq0lctXQzGrN6a\n9DVYf3A9TjeeRnlLOc6PPR8AkBGZoTv4V+7HdVOv0/rMXSLCgzsfxFOLnkKEdwSifaNVZtw923wW\nXi5eCPQIVNs3xDMEjZ2N6O7rVlunzciUjzaJgYn49fZfkV2Vjbvn3I23r3gbLo4uep/HWBz8zczP\nzQ8XT9bZYYoxBmBhzEI4OzpjzbdrsDx5OZwcpM6IKUEpqG6rRlNnk9o+RIR95fuwPGU5Gjoa0NPf\no7bNppOb0N7bjttm3gZA/ctEW8oHkNI14d7hqGqt0vt95NTmYEboDL22DfcOx+5bd+OejHvMPpfP\nSBz8GWNWQQiBNelrsPP0Tlw39bqh5Y4OjpgdPltjL52Klgr0Uz8m+09GmFeYxtG4j+16DC8ueXHo\nebjpEekqbQh5dXk6O2SMNfWTU6dfzd/SOPgzxqzGzdNvxm0zb8MFcReoLM+I0Jz62VexD+dGnQsh\nBGJ8Y9Ty/o2djahtq8WiuEVDy9SCf732mj8w9u6e+qZ9LI2DP2PMani5eOGDqz4YSvkoZERmaOyl\ns79iP+ZFzQMgTcE8Mu9fKCtEUlCSSkplRtgMFDQUDI0d0JX2AcZW869rr0NXX5faSGFrxMGfMWb1\nFHn6kcOF9lfsx7lR5wKQaugj+/prGjTl5uSGqcFTcbTmKPoG+lDQUKDSLXOksQR/Ra3fFiYx5uDP\nGLN6Mb4x6B/oR2XrcE6/u68bx2uPY26ENA1ztE+0WtqnoKFA4/TI6RHpOFh5EMWNxQj3Doeni6fW\nc4/liV5jaey1NA7+jDGrJ4RQ66VzvPY4EgIS4OXiBQAac/6FjVLaZ6T0SCnvP1pjLwBM8p+EAxUH\nUN1aPWo5bSXfD3DwZ4zZiIzIDPxW/tvQa+WUD6A9569prpz0CGng2Gj5fgCYFTYLq9NXY+EHC1HS\nVKJzWw7+jDFmYsuTl2Nj7kas2rIKVa1V6sF/RM5/gAZQJCtCQoD69MhTg6eiuq0ae8/uHTX4CyHw\n6MJH8cC5D+D8D8/HyfqTGrfrG+jDqYZTOp9bbE04+DPGbEJaaBoK7ilAqGcopr85Hd8UfaMS/IM8\ngtDZ14n2nnYA0hgAf3d/eLt6qx3L0cERs8JmYVfJrlGDv8KajDVYt3gdLvzoQtS116mtP1l/EtG+\n0TrbD6wJB3/GmM3wcfXBsxc9i+w7s7F2/lqVlI4QQqXRt6ChQOf0yOkR6XB0cBzTFMorZ6zEuVHn\n4ofTP6it+6XsFyyMWTiGd2NZHPwZYzZnkv8kPHb+Y3AQqiEs2nc49TPag9DTI9ORFJg05vl0Fsct\nxu7S3WrLd5fuVhlMZu04+DPG7Ea0z3Cjr2KAlzZXJV2FDVdvGPM5FsUvUgv+AzSAn0t/VhuZbM04\n+DPG7IZyd88Cme60j7uzO+ZEzNG6XpvU4FS09bShTF42tOxE3Qn4u/vbxMheBQ7+jDG7MbLmb45H\nIgohkBmXqVL7t7WUD8DBnzFmR6J9o3G25Sy6+7pR1VqFeL94s5xncdxi7CrZNfQ6qzQLmXGZZjmX\nuXDwZ4zZjRjfGJQ3l6O4sRixfrFwdnQ2y3kUeX8ikvL9ZT/bXPB3Gn0TxhizDYqunqPl+42VEJAA\nIsLpptNo72lHkEcQIrwjzHY+c+DgzxizG96u3nB2cMaBigNICtTe08dYQggsil+EXSW70NnbiczY\nTLOdy1w47cMYsyvRvtH4seRHs9b8AWBRnJT6ySqzvXw/wMGfMWZnYnxjcLT6qNmD/+J4qdH3l7Jf\nbDL4c9qHMWZXon2iQSCzpn0AIM4vDh7OHnB1dEW4d7hZz2UOHPwZY3Yl2icaXi5eCPMKM/u5Fsct\nNluPInPj4M8YsyvRvtFIDEwcl0cpPrX4KQhY/yMbNeHgzxizKxfGXzj0dC9zs7XuncrEyAcij3sB\nhNHN604AAAUKSURBVCBLl4ExxmyNEAJEZPBtB/f2YYyxCYiDP2OMTUAc/BljbALi4M8YYxMQB3/G\nGJuAOPgzxtgExMGfMcYmIKOCvxDiH0KI40KIo0KIHUKIMKV1jwghioQQ+UKIJcYXlTHGmKkYW/N/\nnohmENEsAN8AeBwAhBBTAawAkALgUgDrxXiMtbZxWVlZli6C1eBrMYyvxTC+FqZjVPAnojall54A\nBgZ/XwbgcyLqI6JSAEUAMow510TAH+xhfC2G8bUYxtfCdIye20cI8TSAWwDIASgeXx8JYJ/SZpWD\nyxhjjFmBUWv+QogfhBA5Sj+5g/9eCQBE9HciigGwEcC95i4wY4wx45lsYjchRDSAb4houhDiYQBE\nRM8NrtsB4HEiOqBhP57VjTHGDGDMxG5GpX2EEFOIqHjw5dUATg3+vhXARiHEy5DSPVMAZGs6hjGF\nZ4wxZhhjc/7PCiESITX0lgH4IwAQ0UkhxH8AnATQC2A1z9vMGGPWw+Lz+TPGGBt/Fh3hK4RYKoQ4\nJYQoFEI8ZMmyjDchRJQQYpcQ4sRgI/p9g8v9hRDfCyEKhBA7hRC+li7reBBCOAghjgghtg6+npDX\nAQCEEL5CiE2DAyRPCCHOmYjXQwjxgBAib7CDyUYhhMtEug5CiPeEELVCiBylZVrf/1gH1los+Ash\nHAD8C8AlAFIB3CiESLZUeSygD8CDRJQKYB6ANYPv/2EAPxJREoBdAB6xYBnH0/2Q0oQKE/U6AMCr\nAL4lohQAMyC1pU2o6yGEiIDUe3A2EU2HlKK+ERPrOnwAKT4q0/j+DRlYa8mafwaAIiIqI6JeAJ8D\nuMqC5RlXRFRDRMcGf28DkA8gCtI12DC42QZIDel2TQgRBeAyAP9WWjzhrgMACCF8ACwkog8AYHCg\nZDMm5vVwBOAphHAC4A5pvNCEuQ5EtBdA04jF2t7/mAfWWjL4RwIoV3pdgQk6EEwIEQdgJoD9AEKJ\nqBaQviAAhFiuZOPmZQB/A6DcADURrwMAxANoEEJ8MJgGe0cI4YEJdj2IqArAiwDOQgr6zUT0IybY\nddAgRMv7HxlPRx1Yy7N6WpgQwgvAlwDuH7wDGNkCb9ct8kKIywHUDt4F6bpNtevroMQJwGwAbxDR\nbADtkG71J9rnwg9SLTcWQASkO4CbMMGugx4Mfv+WDP6VAGKUXkcNLpswBm9nvwTwMRFtGVxcK4QI\nHVwfBqDOUuUbJ+cBWCaEOAPgMwCLhRAfA6iZYNdBoQJAOREdGnz9FaQvg4n2ubgIwBkiaiSifgBf\nA5iPiXcdRtL2/isBRCttN2o8tWTwPwhgihAiVgjhAuAGSIPDJpL3AZwkoleVlm0FcNvg77cC2DJy\nJ3tCRI8SUQwRTYL0GdhFRCsBbMMEug4Kg7f05YPjZwDgQgAnMME+F5DSPecKIdwGGy4vhNQhYKJd\nBwHVO2Jt738rgBsGe0TFQ8fA2qEDW7KfvxBiKaSeDQ4A3iOiZy1WmHEmhDgPwC8AciHduhGARyH9\nh/0H0rd4GYAVRCS3VDnHkxDiAgB/IaJlQogATNzrMANS47czgDMAVkFq/JxQ10MI8TikCkEvgKMA\n/gDAGxPkOgghPgWQCSAQQC2kKfP/C2ATNLx/IcQjAO6AdL3uJ6LvdR6fB3kxxtjEww2+jDE2AXHw\nZ4yxCYiDP2OMTUAc/BljbALi4M8YYxMQB3/GGJuAOPgzxtgExMGfMcYmoP8HvVlwCCkBgAwAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6cfa62ef10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cointegrated,x,y = GenerateData(100)\n",
    "plt.plot(x)\n",
    "plt.plot(y)\n",
    "print \"Cointegrated?\",cointegrated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's try the Bayesian routine to see if the result matches the truth given above when generating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test result: False\n",
      "Congratulations! The result of the routine matches the ground truth\n"
     ]
    }
   ],
   "source": [
    "BF = BayesianLearningTest(x,y)\n",
    "test_result = BF<0\n",
    "print \"Test result:\", test_result\n",
    "if test_result == cointegrated:\n",
    "    print \"Congratulations! The result of the routine matches the ground truth\"\n",
    "else:\n",
    "    print \"Unfortunately the routine disagreed with the ground truth\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing with Dickey-Fuller\n",
    "\n",
    "For comparison purposes we now test for cointegration using the standard test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.tsa.stattools import adfuller\n",
    "slope,intercept,_=LinearRegression(x,y);\n",
    "epsilon=y-intercept-slope*x;\n",
    "adf=adfuller(epsilon, maxlag=1, autolag=None, regression=\"ct\");\n",
    "pvalue = adf[1]\n",
    "cointegrated_adf = pvalue<0.05\n",
    "cointegrated_adf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ROC curve: Bayesian test versus Dickey Fuller\n",
    "\n",
    "Both the Bayesian learning test and the Dickey-Fuller test do the job and provide a test statistic which we compare against a threshold. To compare which test is better, we look at the ROC curve, and in particular, the AUC of the ROC. To do that, we repeatedly generate time series and perform both tests, then plot the resulting ROC curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 250 of 5000\n",
      "Experiment 500 of 5000\n",
      "Experiment 750 of 5000\n",
      "Experiment 1000 of 5000\n",
      "Experiment 1250 of 5000\n",
      "Experiment 1500 of 5000\n",
      "Experiment 1750 of 5000\n",
      "Experiment 2000 of 5000\n",
      "Experiment 2250 of 5000\n",
      "Experiment 2500 of 5000\n",
      "Experiment 2750 of 5000\n",
      "Experiment 3000 of 5000\n",
      "Experiment 3250 of 5000\n",
      "Experiment 3500 of 5000\n",
      "Experiment 3750 of 5000\n",
      "Experiment 4000 of 5000\n",
      "Experiment 4250 of 5000\n",
      "Experiment 4500 of 5000\n",
      "Experiment 4750 of 5000\n",
      "Experiment 5000 of 5000\n"
     ]
    }
   ],
   "source": [
    "from numpy import zeros, mod\n",
    "T=20\n",
    "experiments = 5000\n",
    "\n",
    "cointegratedActual=zeros(experiments, dtype=bool)\n",
    "logBF=zeros(experiments)\n",
    "pvalue=zeros(experiments)\n",
    "    \n",
    "for expr in range(0,experiments):\n",
    "    cointegratedActual[expr],x,y=GenerateData(T);\n",
    "\n",
    "    #classical test\n",
    "    slope,intercept,_=LinearRegression(x,y);\n",
    "    epsilon=y-intercept-slope*x;\n",
    "    adf=adfuller(epsilon, maxlag=1, autolag=None, regression=\"ct\");\n",
    "    pvalue[expr]=adf[1]\n",
    "\n",
    "    #bayesian test\n",
    "    logBF[expr]=BayesianLearningTest(x,y)\n",
    "    if mod(expr+1,experiments/20)==0:\n",
    "        print \"Experiment\",expr+1,\"of\",experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAFwCAYAAAA2ZEWgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VdX9/vH3zsCQkIBIAJkRkKkMBUEUVEAUVNRqrQIW\nxGrFWuWrP0WkqMVq6zxRHClOKE5gBQdkkoCACggGEZA5QICEMSEJkJC7f3+cgEgZMtzcfc/N81or\n607nnvPxrJDHvc8+extrLSIiIn4U5boAERGRklKIiYiIbynERETEtxRiIiLiWwoxERHxLYWYiIj4\n1ilDzBgzzhiTboxZdpJtRhtj1hhjfjDGtA9uiSIiIsdXlJbYG0DvE31ojLkUaGKtbQYMAV4JUm0i\nIiIndcoQs9bOA/acZJOrgLcLt/0OqGqMqRWc8kRERE4sGNfE6gKbj3qdVvieiIhImdLADhER8a2Y\nIOwjDah/1Ot6he/9D2OMJmoUEZH/Ya01JfleUUPMFP4czxTgr8AHxpguwF5rbfqJdqQJh4tv1KhR\njBo1ynUZvqPzVnI6dyXj5/MWsAHmps498jc6eWMy+YF8cvNzWbdnHVUqVDmy7YY9G8jJz2Hrvq1k\nHcziUOAQValHq/wbqWJqYTDE7GvM/AWQubMybOrqfdFGQSCWESPg6quhUyfvbWNKlF9AEULMGDMB\n6A6cbozZBPwdqABYa+1r1tovjDGXGWPWAjnATSWuRkREykRBoICsg1ns3r+bKT9PYeySsWzO2kyt\n+Fqs27PuyHbdG3UHYPf+3fRp0oekuCQqm2o0qdqcHTsgYKF1kmXDymqk/nwasz6uD3lVaPXbajRv\nDk2bFu6oEpx9KfTtCw0bQo0aEFUGF7BOGWLW2gFF2OaO4JQjIiKHWWsJ2ADbsrexP38/W7K2kHkw\nk9S9qVSMqUjABvgx/UcSKiYAsDxjORWiK7Bx70ay87JZt2cd1SpVA2Dvgb1H9hsbFcvvW/2ex3s9\nTqukVhQUwOYfG3IoLwZjDF9/Ddu3w1ujIL2wXy0uDnJzoUMHyM+HzExo1Qqe/wfceWfZBFRRBOOa\nmJSx7t27uy7Bl3TeSk7nrmROdt7yC/LZc8C7WyktK42fd/3M2t1rj3yedTCL79K+I68gj0oxlVi8\ndTHZedlHPm9avSk7cnbQtHpTqlSoQosaLbDWkpOfQ+PTGgPQqU4nWia1JMpEUSu+FknxSZxR5Qzv\n+Plg91dj6lTD+vXwyLXwWRUwBvbt845x8cXe6x07vLC6/Xa4/HI46yxISAjyyQoSE8prVMYYq2ti\nIhIpDhw6wLxN85i3aR7T100nKT6J3Pxcfkz/kaT4JHbk7CA9J5242Dhy83MBSIpLYkfuDs487Uya\nVm9Ko6qNSIpPwlrLocAh2tZqS52EOhTYAlrUaEHdhLrFumaUmwsbNnhhlJYG06bBq69C9i95yNVX\nQ+3aXgsqLg6qV3cbUsaYEg/sUIiJiJxEfkE+GTkZLNm2hFsuuoWMtAzXJflWw4YN2bhx4/+8X5oQ\nU3eiiJRL1lqWZywnYAP/89mO3B2kZaUxePLgX3+QphHWpVGaUYgn3KdaYiISyZalL2P2htmMXjia\n9XvWc3rl0zHGsDN3JwBta7X9n+9k52VTpUIVqlasygfXfkDtKrUxxhxuMYT6PyFinOj8qTtRRMq1\n/IJ89uXt4+Hkh8nNz2VT1ib25+9n1c5V7MjdQb3EenQ4owO3dbyNs+ucfeR7VStVpUJ0hSIfRyFW\nOgoxERHgpUUvsXb3WlLSU1i7ey2bMjcd+ezBCx6kTkIdqlSoQq34WrSv3Z6k+KSgHFchVjoKMRGJ\naNZaduTuOPKHbnv2dlbtXMX3275n/Z71TFo5idioWPID+VzY8EIua3YZ1SpVo3PdzrSr1a5Mrrkc\nTSFWOgoxEYkYBYECZm2YxVspbzHhxwm/+qxmfE0AMnIyaFStEadVOo2ejXvSrlY7Lmx0IdUrV//V\nNEihEu4h1qhRIzIyMoiNjSU6OppWrVoxcOBAbr31Vowx3HTTTUyYMIGKFStircUYw7hx4/jDH/4Q\nkvrKIsQ0OlFEgi4nL4dNmZvIyMlge/Z29h/aD3jz83259ks+WvHRkW27NejGPefew91d7qZapWpU\niqlEdFS0q9J9zRjD559/To8ePdi3bx9z5sxh6NChLFy4kHHjxgEwfPhw/vGPfziuNHgUYiJSIluy\ntrA5czML0xaycudKXl/6OvmB/CPdfQD1E+tTu0ptEiomUD+xPgEbIGADvNDnBa5rfR214muVeRdg\nsOza5bqCojnc0klISKBv377UqlWLc889l3vuucdxZWVDISYiJ2StZUvWFrZnbyc1M5Xt2du5c+qd\nRz6vWrEqtarUolVSK+7reh93dr6T0yqfhsEQGx3rsPLgyM+H8eNhzBhYutR1NSXTqVMn6taty9df\nf+26lDKhEBORI/IL8pm+bjpfb/qanbk7GbfU64KqVqka1StXp35ifa5vfT0v9HmBxIqJVI6t7Lji\nshEIwNNPw/Dh3uuzz4Z586BbN7d1lVSdOnXYvXs3AE899RRjxozBWktsbCwZGf6egUQhJlJOZR3M\n4qr3r2LFjhXUrlKbFTtWcChwCPCuU519xtm8eNmLDG4/mLjYOMfVlq1AADZu9OYaHDYMUlLgwAEY\nOBDeesubh7AogtUzGuyxI2lpaVSvXh2AYcOG6ZqYiPjLgUMHWLx1MVEmirW717IwbSEvLnoRgFcu\nf4XOdTsD0Pi0xkeW7ohkr70G8+d7y4dMmvTLLO5JSVC/Pnz8MbRvD2ecUbz9huPAxUWLFrF161a6\ndevGt99+67qcoFOIiUSIvII8svOy2bpvK7M3zOaV719hU+amXy3ncW69c9lzYA91EuowtPNQRl4w\n8shw9vJgyxbo2hU2bYIBA6BXL2914d69oUkT19UF1+HRiXfddRcDBw6kdevWrksqEwoxEZ/JOpjF\ny4teZk7qHApsAdv2bSM1M5Wsg1kAxMfGE7ABep3Zi2cveZZWSa04Pe50KsdU9s1IwJLKzIQHH4SM\nDK9rLzrae2/9eti503sfYOFCL7wi0RVXXEFMTAxRUVG0atWKe++9lyFDhrguq8zoZmeRMFcQKODm\nKTezLH0ZG/duZM+BPUSbaK5ueTVXNb+K2KhY6ibWpVVSK6pXru663JDKz4esLFiwAJ55BubM8d6/\n+25o1MhbJwu8bsMGDaBxY6hbt+THC/ebncOdbnYWKUe+Tv2aZ799lrmpc9m9fzej+4ymdpXadKrb\niQZVGxBlHK0H71heHjz3HDz0kPf8sMaNYdw4+NOf3NUmoaeWmIgjufm5LEpbRNq+NFL3pvJD+g+k\nZ6czJ3XOkW16Nu7JrR1upXmN5rSv3d5hte5ZC2PHwuGesd694amnoE2b0NWglljpaO5EEZ/JL8gn\nYAO8nfI2eQVes2Fz1mY27t3IBz99AMBva/+Wxqc1xmC4+MyLSYpPoku9LtRJqOOydKeshdRUWLnS\n6yLcsAG+/NLrOhw4EJ54ovgjB4NBIVY6CjGRMFYQKGBZ+jKeWvAUX234ivSc9F99/rsWv6NOlTpY\nLNEmmpZJLbm90+2Oqg1f48fDoEHe8/h46NDBG/beuTP06QPNm7urTSFWOromJhImJq6YyJyNc0hJ\nT8FiWbNrzZHQSqyYyJ/a/4k/tv0jrZJaReysFmXh9dfh5pu94e5Ll0JCguuKJNwpxESKYHv2dvIK\n8tiZu5P7Z97PjPUzOK/+efRt1pfWNVtTIboCDao2oEWNFuV2wEVxbdgA06d7XYQ7dsDo0XDwIJx3\nnjfFU4TfDSBBou5EkROw1jJv0zwuePMCwJuRPSMnA2MMz/V+jiEdh0T8fVfBtG8frFoFX33lhddX\nX0FsLPTv710Da9rUu97VuLHrSk9M3Ymlo+5EkTKUk5fDrA2zGL9sPKt3rWZZ+jIAzqhyBt/c/A0N\nqzV0XKF/ZGbCXXd5cxGuXeu1ug5r2BB69oRPP4W+fd3VKJFBLTEp9/IL8pm5fiaXTbgMgC71utC3\nWV+61OtC1wZdqRRTyXGF/rBpE7z7Lvzzn5CT4703YgR07OjdaNyihf+vcaklVjpqiYkEwZasLSzY\nvIAdOTt4YPYD7D2wF4B2tdqx8M8LqRBdwXGF/pCXB1OmwMMPw/Ll3ntVq8KFF3rXtyJtLkI/aNSo\nERkZGcTExBAbG8t5553HK6+8Qt3STFMS5tQSk4hXECggJT2Fu768i683eQsD1oirQZd6XUismMiw\n84bRrlY7Xd8qhkce8WbMAGjb1mt9XXQRVI7wgZjh3hJr3Lgxr7/+Oj169CAvL4+//OUv7Nmzh48/\n/th1aUDZtMQ0jEoiTuaBTGasm8GASQNoOropMY/E0PG1jizZtoTHLnqM3fftZsewHXza/1PeveZd\n2tdurwArhnff9QLsppu8ARkpKd61rUgPML84HBIVKlTg2muvZcWKFQB88cUXdOjQgapVq9KwYUMe\nfvjhI9/p27cvL7744q/2065dOyZPngzAqlWruOSSSzj99NNp2bIlH3300ZHtvvjiC1q3bk1iYiL1\n69fn2WefLev/xF+z1obsxzucSHAFAgH7/o/v28e/ftxGPRxlGYVlFLbZ6Gb2rql32dkbZtvsg9mu\ny/S1rCxr33zT2qFDrQVru3SxNhBwXVXohfvfsEaNGtlZs2ZZa63NycmxN954ox08eLC11to5c+bY\n5cuXW2ut/fHHH23t2rXt5MmTrbXWfvjhh/acc845sp8ffvjB1qhRwx46dMjm5OTY+vXr27feessG\nAoEjn61cudJaa+0ZZ5xh58+fb621du/evXbp0qUnrO9E56/w/RLliroTxdf2HdxH4uOJAPQ6sxfV\nKlXjud7PUS+xnuPK/Csnx5ujMC0NJk70Vjw+rE8fqF3bm4C3WuSvnfk//NCduGvXLmJiYsjOzqZm\nzZpMmzbtuGuJ3X333URFRfHMM89w8OBB6tSpw8KFC2nSpAnDhg1j//79jBkzhg8//JAXX3yROXN+\nmdPztttuo27dujz44IM0atSIkSNH0q9fPxJOMXJHAztEgIAN8I85/+C7tO/4cu2XAGz8v40aAl9C\n1sJHH3k3GE+Z4s1ZCPD738Nll8Hgwd7UT9HRTsv0BfNwcLql7d9LHpSTJ0+mR48eWGv55JNPuOCC\nC1i5ciUbNmxgxIgRLF++nLy8PPLy8vjDH/4AQMWKFbn++ut55513eOihh3jvvfeOXEdLTU3l22+/\npXrhujbWWgoKChhUODfYpEmTeOSRRxg+fDjt2rXjscceo0uXLqU8A0WnEBNf2Jy5meUZy5mwfAKT\nVkxi/6H9XNvqWj689kMubnIx1SqVw2ZBCRQUwIQJEBMD69bBTz/BJ5/AgQPeKsc33wxXXAHty/eE\n+SVWmvAJWg2FLR1jDFdffTVDhgxh3rx53HfffQwdOpRp06YRGxvL3Xffza5du458b9CgQQwcOJCu\nXbsSHx9P586dAahfvz7du3dn2rRpxz1ex44d+eSTTygoKODf//431113HZs2bSr7/9BCCjEJa9l5\n2fSb2I/P13xO1YpVqV2lNg9c8ACD2g1Sl2ExzJwJixbByJFey6tfP2+Kp0qVvBnhr7vO6yaUyDJ5\n8mT27t1Ly5Ytyc7O5rTTTiM2NpaFCxcyYcIEevfufWTbLl26EBUVxT333MPAgQOPvN+3b19GjBjB\nO++8Q79+/bDWkpKSQpUqVWjSpAkfffQRffv2JTExkYSEBKJD3WQv6cW0kvwQ5hdFJTwcPHTQPpz8\nsO09vveRQRpv//C267J8ZckSa//2N2sHD7a2XTtvMMZll1l7yy3W/vCD6+r8K9z/hjVq1MjGxcXZ\nhIQEm5iYaNu0aWPfe+89a621EydOtA0bNrSJiYn2iiuusHfeeacdOHDgr77/6KOP2qioKLthw4Zf\nvb969Wp7+eWX26SkJFujRg170UUX2ZSUFJuXl2f79Oljq1evbqtWrWo7d+5sFyxYcML6TnT+0MAO\n8TtrLdPXTWfm+pk8/c3TAAxqN4g2Ndtwe6fbiYuNc1yhP/zf/3k3GoO33la/flCjBnTp4k31JKUT\n7gM7Smv8+PGMHTuWuXPnlsn+NbBDIs5l717GvE3z2Je3D4C2tdoy8vyRPNLjEd27VQwHD8Ktt8Lb\nb3vXtN5/H+KU+1IMubm5vPTSS9xxxx2uSykWhZg489nqz5i6dirT/jiNZtWbUb9qfWKi9Ct5KlOm\neCG1a5e35taOHb989sILMHSou9rEn6ZPn84111zDJZdcQv/+/V2XUyzqTpSQ2p+/n/tn3s+YRWMI\n2ACtk1qz/PblrssKe4EALF7szYyxYwdccIHXVRgXB+ec43UdVq3qusrIF+ndiWVN3Yniay8vepnb\nv7gdgDs63cHfu/+dGnE1HFcV3tau9YLr55+91/Hx3jRPbdu6rUskXCjEpMykZ6fzxZoveGL+E/y8\ny/srPKDNAMZfPV6rH5/C2rVw++0wY4Y3MOOrr6BHD9dViYQfhZgE3bZ927j1s1v5bPVnVIyuSNcG\nXXm81+P0aNSDqpXU53Uiycnw+eeweTN88IG3eOTMmd7s8CJyfLomJkGzYc8GRn83mue/e57KMZV5\n+fKXubH9ja7LCmsFBd5che3be6sfN20KN9wAjRp50z1JeNE1sdLRNTEJO7v37+b7rd8zdslYPlrx\nERWiK3Bn5zt57KLHiK8Q77q8sJWeDn/+M3z66S/v/fgj/OY37mqSU2vYsKFu/SiFhg2DP7+pWmJS\nIpkHMqn2hDdfYVxsHO1qteOOzncwoM0Ax5WFtz174PLL4ZtvvNcTJsDVV3vTP4mUV2qJSchYa7l/\n5v08ueBJAPIeyCM2OtZxVf4wZw507+49f+klGDTIG20oIiWnEJMisdZy66e3MnHlRPYe2MvwrsN5\nvNfjrsvyhW3b4D//8VZDjonxWmNVqriuSiQyKMTkuAoCBezI3cEHyz9gyuopLM9YTkZOBsO7DueW\nDrfQtHpT1yWGvTlzvDW5Dq92cffdMGqUAkwkmBRiAsCOnB08veBp9h/az48ZP5K8MfnIZx3O6MDo\nPqM5p945NKrWyFmNfpGX563V1b07/Pa33hIojRu7rkokMmlgRzmWX5BPSnoKncZ2AqBCdAX+1u1v\nFNgC2tVqx9Utr9ZNycW0fbs3BRRA/fqwcSNE6RSKnFRpBnYoxMqh9Ox0pq6dyk2TbwLgrNPP4vtb\nvyc+Nl7Dh0vh8CKTh59XqOC2HhG/0OhEKZI5G+fQ8+2eBGwAg+GSJpcwud9kKsVofHdp7NkD1av/\n8nrjRgWYSKgoxCLUgUMHGLNwDJ+v+ZzEionsyt3F/M3zaZXUilmDZlG7itaiL60ffvB+bvIatGRk\nePMcqjErEjrqTowwW7K2UP+5+kdeX9PyGq4860qqVKhCzfianN/wfIfVRYann4Zhw7znrVpBzZre\nHIfR0W7rEvErdScKH6/8mCfmP8HCtIVEm2i23rOVmvE1XZcVcQYPhrfegksvhZdf9ibpFRF3NG7K\n5xamLaTJ6Cb8/sPfc1ql05gxcAaHHjqkAAuyvXu9xSffestriX3xhQJMJByoJeZTW7K28OriV3n0\n60epm1CXpUOW0r52e9dlRZw9e+D+++G117zXX34JvXu7rUlEflGka2LGmD7A83gtt3HW2ieO+TwR\neAdoAEQDz1hr3zzOfnRNrJTGLBzDmIVj+HnXzyRWTGRwu8G8cOkLrsuKWJde6gXXk0/Cvfdq0IZI\nWSjT+8SMMVHAauAiYCuwCOhnrV111DYjgERr7QhjTA3gZ6CWtfbQMftSiJXCrtxd1HiqBjf/9maG\nnjOUtrW0Rn1ZefNNuOMOb62vl16Cv/zFdUUikas0IVaUa2KdgTXW2lRrbT7wPnDVMdtYIKHweQKw\n69gAk9L5bPVn1HiqBgBjrxirACtDb7zhDZs//3zYskUBJhLOihJidYHNR73eUvje0cYArYwxW4EU\n4P+CU54AfLLqE6547wp6Nu7JwQcOalaNMjRhAvzpT3DxxTB1KtQ99jddRMJKsAZ29AaWWmt7GmOa\nADOMMW2ttdlB2n+5UxAooN+kfkxcMRGAc+udy6xBsxxXFZkyM+GVV7wBHOBN3Pvf/zotSUSKqCgh\nloY3YOOweoXvHe0m4DEAa+06Y8wGoAWw+NidjRo16sjz7t270/3wKoHyK53/05kl25Yw7spxDGo3\niJgoDSQNlsxMePZZ2LTJu/Z12BVXeC0xLZUiUraSk5NJTk4Oyr6KMrAjGm+gxkXANmAh0N9au/Ko\nbV4EMqy1DxtjauGFVztr7e5j9qWBHadQECjgiflPMPKrkUy9YSp9mvZxXVLEmD8fHnkEpk3zXg8b\n5k0TNWQIVK3qtjaR8qxMZ+yw1hYYY+4ApvPLEPuVxpgh3sf2NeBR4E1jzLLCr913bIDJie3M3cmm\nzE0s3baUkV+NJD0nnZva36QACwJrYcECuOYab27Dli29gRuDB7uuTESCQXMnOrR+z3oen/c4Y5eM\nJT42nmqVqtGiRgse7v4wXRt0dV2e7337rXdjclYWNGsGkyZBmzauqxKRY2nuRB/5fPXnLNm2hDdT\n3mT9nvUkxSXxRK8nuK/rfa5L872ffoLp0+GDD+C777z32rf3Xp91ltvaRKRsqCUWQi8teom/fvFX\nep3Zi5rxNfnL2X+hW4Nursvyvf37oVs3WLLEC6vGjeGvf4VevaByZdfVicipqCUWxgI2QKsXW7E9\nezuZBzO5reNtvNz3ZddlRYx//QueeMLrMpwxwwsuESk/FGJlKL8gn8YvNCZtXxpf3vAlbWq1oU5C\nHddlRYzhw705Da+80hsy36SJ64pEJNTUnViGUran0P7V9iz+82I61unoupyIkp4OtWt7Q+YfeMB1\nNSJSGmU9d6IUk7WWsd+Ppf2r7WlyWhMFWBAtXerNJn94La+RI93WIyJuqTsxyGaun8nF4y8GoE/T\nPkzuN9lxRZFj/Hh46ilYv967Ufn3v9fSKCLlnboTgyhgA0T/I5qLz7yYaX+cpol6gyA/H0aMgHXr\n4JNPoF07mDgRmjZ1XZmIBEuZricWTJEcYhk5GdR6uhYAWfdnkVAx4RTfkOM5eBC++cabgHf1am9B\nSoCbb4Zzz/UeRSSyKMQc25y5mQbPe3MkH3zgIBWiKziuyJ9efhluv9173rAhDBzoTRM1YIDbukSk\nbOk+MUestbR6qRWrdnqLXGfcm6EAK6G+feHzz73Hd97RhLwiUjQKsVJ4cv6TrNq5iu9u+Y7OdTu7\nLseXCgrgllu8AHvrLRg0yHVFIuInCrESenL+k9w/635Gnj9SAVYK554LixbBQw8pwESk+BRixZCb\nn0vK9hSWpS9j+Mzh3HPuPTza81HXZfnSoUPelFGLFnnre11yieuKRMSPNLCjiDZlbqLh894dti1r\ntKRRtUZ8NuAzoozuFy+JCy+EuXPhhhu8a2AiUn5pYEcZ+2bzN5z3+nkApN6VSoOqDRxX5F95eXDm\nmZCW5o1GvO021xWJiJ8pxE4iJy+HKo9VAaBxtcYs/PNCasTVcFyVv51/vhdgS5bAb3/ruhoR8TuF\n2Em8tOglADLvzySxYqLjavxv505YvBjef18BJiLBoRA7ges+uo6PVnzEDW1uUIAFSe3aEAh4C1iK\niASDRiUcxwNfPcBHKz5i/NXjGX/1eNfl+N706d6aXwUF3lRSdeu6rkhEIoVaYkf5ZvM3XPrupWQe\nzGTk+SP5Y9s/ui7J11asgD/8wXvs0cO7mblZM9dViUgk0RD7QtPWTqPPu32oHFOZGQNn0LVBV9cl\n+Zq1EB/vTej75pvePIgiIsejRTFL6ZNVn9Dn3T50rd+VXfftUoCVQn4+PPEEREXB/v3w888KMBEp\nO2qJAXH/jKNeYj1W37nadSm+lZnp3bx85ZXe6/PPh0mTICnJbV0iEv7UEiuh1L2pxD4Sy/5D+5l6\nw1TX5fjS3LnezcvVqnkB1rIl/PCD974CTETKWrkNsYJAAYMnD8ZgyLg3gybVm7guyXfeecebPsoY\n+Owz7zrYihXe6ssiIqFQbkcnPr3gaZI3JjPuynEkxavJUByzZ8Prr3shdsEF3uuocvu/QyLiUrm8\nJnYocIg6z9ThkiaX8M41mn22qAIBePttuOkmSEiAsWPh+utdVyUifqcJgIvBWkvFRysSsAHuPe9e\n1+X4grXe/V6TJnmvmzWDVavU+hIR98pdiO3ev5uADbD57s3US6znupywFwhAdLT3/O9/hxEjoGJF\ntzWJiBxW7kLswdkPAijAiujSS73HvXuhalW3tYiIHKtcdQjt3r+blxe/zN+6/c11KWFv3jxITPTm\nPfzgAwWYiISncjOwIy0rjXrPea2vAyMPUDFGfWLHs2IFdOjgTRd19tnesilNdPeBiJQh3ex8CvkF\n+QqwIhg7Flq3hipVvJuVFy1SgIlIeCsXITZ93XQAfr7jZwXYCbzxBtx6qzfP4c6d3rRRIiLhrlx0\nJ/Z6uxf78vbx3S3fhfzYfvD663Dzzd5MG0uXejNwiIiEiroTT8BaS8sXWzJrwywGttVU6sezcaMX\nYNdf7815qAATET+J2JbYvoP7SHw8EYAZA2fQ68xeITmunzzyCDz0kPd861Y44wy39YhI+aQZO46j\nx1s9AMh/MJ+YqIj9zyy29HT48ENv9o05c+Cii+Djj73h9CIifhORf93v/vJuvt/2PQv+tEABdpTc\nXOjYEdLS4Lrr4L33oF8/11WJiJRcxP2Fz87L5vnvnueRHo9wbv1zXZcTNqZOhcsu857PmePNPi8i\n4ncRd02s/SvtSUlPIfP+TBIrqo8MYPFi6NQJevaEWbNcVyMi8msanVjokvGXkJKewpzBcxRghVas\n8AKsRg2YMcN1NSIiwRUxIfbCty8wY/0MZg2axQUN1VcGMHSoNwNHfDxs26alU0Qk8kTEn7Xd+3dz\n17S7uO+8++jZuKfrcsLC6tXw7397y6dkZ0NMxF39FBGJkGtiF7xxAV9v+prAQwFMOb9b11r45ht4\n8EH4/nu3iYTeAAAYH0lEQVRvCRURkXBWrq+JHTx0kK83fc0bV71R7gNs/Hivy7BrV9i0CV580XVF\nIiJly/edTHdPuxuAG9vd6LgSt9LSYNAg7+blL79U96GIlA++/lO3bd82Xl78Mo/2eLTct8KefBJi\nY2HmTNeViIiEjq+7E99Z9g6VYyoz8oKRrktxbvRobxJfEZHyxNchNvKrkVzR/ArXZThlrbcCM+ga\nmIiUP74NsbW715IfyOefPf/puhRndu2C3/zGG4X4n/9oEl8RKX98eU3MWstDsx+iWqVqNK3e1HU5\nTgQC3kCO1ath5Upo0cJ1RSIioefLELt5ys28t/w9nuz1pOtSnLnwQpg3DyZOVICJSPnluxB7aPZD\nvPHDG7x7zbsMaDPAdTkhl5PjhdaWLd5ciL201qeIlGO+uiZ28NBBHpn7CHd0uqNcBhjAtdd6Afbf\n/yrARER8FWJPzve6Dx/v9bjjStzIz/duZH7gAfjd71xXIyLinq9CbPTC0dzY7kbiK8S7LiXkkpOh\nQgXv+bBhTksREQkbRQoxY0wfY8wqY8xqY8zwE2zT3Riz1Biz3BgzO5hFBmyAx75+jJ25OxnRbUQw\nd+0Ly5ZBjx7enIjWaii9iMhhp5zF3hgTBawGLgK2AouAftbaVUdtUxVYAFxirU0zxtSw1u48zr5K\nNIv9Mwue4d4Z9zK4/WDeuOqNYn/f7w7PqLVqFTRv7rYWEZFgK80s9kUZndgZWGOtTS082PvAVcCq\no7YZAEyy1qYBHC/ASmrm+pncO+NeBrQZUC4D7LbbvMeCAi1qKSJyrKL8WawLbD7q9ZbC9452FlDd\nGDPbGLPIGDMwWAUmb0ymeuXqvHP1O8HapW8MGwavvgr33qsAExE5nmDdJxYDdAB6AvHAN8aYb6y1\na0uzU2stYxaOYUjHIeVulvp9++Dpp+GFF2DoUNfViIiEp6KEWBrQ4KjX9QrfO9oWYKe19gBwwBgz\nF2gH/E+IjRo16sjz7t2707179xMeeNiMYWQezGRw+8FFKDOyHD4tt9/utAwRkaBLTk4mOTk5KPsq\nysCOaOBnvIEd24CFQH9r7cqjtmkB/BvoA1QEvgOut9auOGZfxRrY0f6V9nSq04mxV44t8nf8bv9+\niIvznn/6KfTt67YeEZGyVqYDO6y1BcaYO4DpeNfQxllrVxpjhngf29estauMMdOAZUAB8NqxAVZc\n2XnZpKSn8Fzv50qzG9/59FPvMSsLEhLc1iIiEu5O2RIL6sGK0RJ7f/n79J/UH/v30NUXDjp0gGrV\n4KuvXFciIhIaZT3EPuQOBQ7Rf1J/Lmt2metSQuqDD2DpUpg1y3UlIiL+EJYDt19Z/AoAH1/3seNK\nQmfFCujXz1uluWdP19WIiPhDWIbY1LVT6dO0DxVjKrouJWQGD/a6ERcscF2JiIh/hGWIbc7czPWt\nr3ddRshcdBEsWgRvvw2xsa6rERHxj7C7Jmat5ceMH2lUrZHrUkJi9WpvEMe8ed4EvyIiUnRh1xIb\n+F9vxqpuDbo5rqTs7d7tTehbu7YCTESkJMIuxHLzc3m176vERIVdIzGoJk6EpCTvedqx85+IiEiR\nhFWIbcnawn9X/ZcacTVcl1Km3nsP/vAHuPxy+PZbTe4rIlJSYdXcGfzJYACuaXmN20LK0H33wVNP\nwSWXwJQprqsREfG3sGoDzNowiwnXTHBdRpnZudMLsDvvhC+/dF2NiIj/hU2I3TPtHgCuanGV40rK\nRlraL9fARo/+ZbVmEREpubAIsVU7V/Hst88yotsI4mLjXJcTdCNGQL163j1gGsQhIhI8YRFir33/\nGnUS6vCvi/7lupSgWrvWa3E9/jiMHAl5eVCnjuuqREQiR1gM7Ji/eT79f9PfdRlBZS1ceaX3fOdO\nOP10t/WIiESisGiJ7crdxbn1znVdRtCkpnrzIK5cCVOnKsBERMpKWLTE1u9ZT8NqDV2XERR79kCj\nRt7zlBRo29ZpOSIiEc15iG3btw2LpfnpzV2XUmpZWXDOOd7zQEAjEEVEyprzEJu4YiIGQ0LFBNel\nlMq6ddC0qfd8wQIFmIhIKDgPsTW713D9b/y97MrkyfC733nPs7MhPt5tPSIi5YXzgR3jl42nThX/\njjvfs8cLsJo1YeNGBZiISCg5bYllHshk74G9/KXTX1yWUWLWQvv23vMffoAzznBbj4hIeeO0JZZX\nkEeUiaJp9aYuyyiR9evhvPNg0ybYsEEBJiLigtOW2JPznyRgAy5LKJHsbGjSxHs+ceIvQ+pFRCS0\nnIbYa0te4+oWV7ssoUT+/GfvUcPoRUTcctqdWDG6ou/mS0xPh/ff9+ZCVICJiLjlNMQSKiYQGxXr\nsoRi2bfPm40e4JFH3NYiIiKOQyy/IN/l4YslLQ0SE+HQIVi0SK0wEZFw4OyaWG5+LpuzNpNYMdFV\nCcXy3HNecOXmQqVKrqsRERFw2BJbs2sNFaMrkhSf5KqEYhk/Hm67TQEmIhJOnIWYMYazTj/L1eGL\nJTMTMjLg3ntdVyIiIkdzFmKL0hax58AeV4cvlj/+0Xts3NhtHSIi8mvOQuyjFR9RM76mq8MXmbXw\n2WfwzDMazCEiEm6cDeyoVqkag9oNcnX4IrvnHu/x7rvd1iEiIv/L6TUxQ3g3bQIBb1Ti/ferFSYi\nEo6cL8USzm67zXv8xz/c1iEiIsfnLMS+WPMFeQV5rg5fJJ995nUnxvpnUhERkXLFWYhlHcyiZ+Oe\nrg5/SsOHw7Zt8Bd/LnUmIlIuOJ3FPhxHJ1oL558P8+d7gzkOL7kiIiLhx0mIbd231Tt4lNMMPa4O\nHbxVml99FW66yXU1IiJyMk5SZPq66QBER0W7OPxxWQvXXOMF2Ny5XmtMRETCm5NrYrM3zuaq5le5\nOPQJff45fPIJfPqpAkxExC+chFjK9hQub3a5i0Of0MMPQ/360Lev60pERKSoQt6dmJOXQ0p6Cg2q\nNgj1oU/ozTdh8WKYPt11JSIiUhwhD7F3f3wXgIubXBzqQx9Xq1awciXcdRdcHB4liYhIETkZ2DGg\nzQCijPvJQlJSvABbtQqaN3ddjYiIFJeTJImPjXdx2F/Jzob27aFyZQWYiIhfuW8OOTJqlPf4009O\nyxARkVIotyH2zDPezcxa6FJExL9CHmIFgQKstaE+7BGHDv3SffjSS87KEBGRIAh5iD04+0HyA/mh\nPuwRDz0Eq1fDlClQqZKzMkREJAhCHmKHAocYes7QUB8WgF274LHHvG7EK65wUoKIiARRyEOsWqVq\nVK9cPdSHBeDSS73HsWOdHF5ERIKs3AzseO01WLQI/vMfiA6feYdFRKQUyk2IDRkC3brBzTe7rkRE\nRIKlXITYe+95j1Onuq1DRESCq1yE2PLl0L8/VKniuhIREQmmkIfYlqwtIT1eXh78619Qp05IDysi\nIiEQ+pudbQG14muF7HirVnmPh6eZEhGRyFGkEDPG9DHGrDLGrDbGDD/Jdp2MMfnGmGtOtr/KsZWL\nW2eJBQLQrJm6EkVEItEpQ8wYEwWMAXoDrYH+xpgWJ9jucWDayfbXvnb7klVaQpMmwf79IT2kiIiE\nSFFaYp2BNdbaVGttPvA+cNVxtrsTmAhkBLG+Unv2WTjvPNdViIhIWShKiNUFNh/1ekvhe0cYY+oA\nv7PWvgyY4JVXOjNnQm6uN9WUiIhEnmAN7HgeOPpamfMgy86Giy+GunXhzDNdVyMiImUhpgjbpAEN\njnpdr/C9o50NvG+MMUAN4FJjTL61dsqxO0v9JJVR20YB0L17d7p3716Csk/t97/3HlevLpPdi4hI\nCSUnJ5OcnByUfZlTre1ljIkGfgYuArYBC4H+1tqVJ9j+DeBTa+3Hx/nMnv/6+cy9aW6pCz+ZuXPh\nwgu9ZVcefrhMDyUiIqVkjMFaW6IevFO2xKy1BcaYO4DpeN2P46y1K40xQ7yP7WvHfuVk+6tdpXZJ\n6iyW11+HFi0UYCIika4o3YlYa78Emh/z3qsn2PZPJ9tXhegKRS6upFJT4cYby/wwIiLiWMhn7Gha\nvWmZ7n//fkhOhrZty/QwIiISBiIuxN55x3s8vACmiIhErpCHWNtaZdtE2rYN+vYF43yQv4iIlLWQ\nh5gpw1vI9u6Fv/8devQos0OIiEgYiaj1xHJyIC4O/t//c12JiIiEQsSEWEYGtG4NBw64rkREREKl\nSEPs/aBfP8jMhO+/d12JiIiESkS0xFJTYfZsmDABOnRwXY2IiIRKRITY/v3QvDn07++6EhERCaWQ\nh1idhDpB3+fu3ZCVFfTdiohImDvlBMBBPZgxNr8gn5io4F6K69gR9u3TjPUiIn5UphMAh7vcXFiy\nxLseJiIi5Yvvr4ktXuw96nqYiEj54+vuxIICuOgi75rYsmVB2aWIiISYr7oTo0100PY1dCjMmfPL\npL8iIlK+hH7uxCDOzLtiBfzzn3DDDUHbpYiI+Ihvr4nt2+etG3buua4rERERV3wbYqmpEB2tGetF\nRMozX4bYnj1eF2JcnOtKRETEJV/eJ9a+PWzapAEdIiLlXciH2Jf2eHv3wmmnwfLl3tIrIiLib6UZ\nYu+rEEtPh9/8BnbuhBCWLSIiZag0Ieaba2Jz50Lt2l6ALVzouhoREQkHvrkm9tJLcNZZsGoVBPFW\nMxER8THftMR++gnuuEMBJiIiv/BFiL37rjeQ45xzXFciIiLhxBcDO5o1gwYNYOZMtcRERCKNryYA\nLq7du2HtWm+9MAWYiIgcLexbYjVrwo4d3rIrUb7o/BQRkeKI2CH2zzzjBdiiRQowERH5X2EdDUuX\nwrXXwtlnu65ERETCUViH2LvvQu/erqsQEZFwFbYhlpLiPfbt67YOEREJX2EbYs8/D/XqeVNNiYiI\nHE/Yhtjbb8Of/+y6ChERCWdhO8S+Rg1YssS7yVlERCJXxA2xX7cOdu2CKlVcVyIiIuEsLENs8mSo\nXx+qV3ddiYiIhLOwDLHZs6FFC9dViIhIuAvLEFuzBv74R9dViIhIuAvLEMvJgaZNXVchIiLhLuxC\nbPt22LIFTj/ddSUiIhLuwi7E5s2D2Fg46yzXlYiISLgLuxCLiYHLLtPaYSIicmphF2JbtkB+vusq\nRETED8IuxNavh7g411WIiIgfhF2IvfUWNGzougoREfGDsJs70RhYvhxatw5RUSIi4lTEzJ14ON/q\n13dbh4iI+ENYhdioUd5jfLzTMkRExCfCKsS+/hpuvx2io11XIiIifhBWIbZ8OVxwgesqRETEL8Iq\nxLKy4Le/dV2FiIj4RdiEmLVw8KAWwhQRkaILmxDbvNl7rFnTbR0iIuIfYRNi8+ZBpUre3IkiIiJF\nETYh9tNP0LOn6ypERMRPwibE5s6F2rVdVyEiIn5SpBAzxvQxxqwyxqw2xgw/zucDjDEphT/zjDFt\nilvIzp1wxRXF/ZaIiJRnpwwxY0wUMAboDbQG+htjWhyz2XrgAmttO+BRYGxxirAWVq2CBg2K8y0R\nESnvitIS6wyssdamWmvzgfeBq47ewFr7rbU2s/Dlt0Dd4hRxeM5E3SMmIiLFUZQQqwtsPur1Fk4e\nUrcAU4tbiDFazVlERIonqAPajTE9gJuAbsX5XkrKL60xERGRoipKiKUBR1+tqlf43q8YY9oCrwF9\nrLV7TrSzUYenqge6d+9O9+7dmTULOnQoaskiIuJnycnJJCcnB2Vfp1wU0xgTDfwMXARsAxYC/a21\nK4/apgEwCxhorf32JPs67qKYDzwAqakwfnyJ/htERMTHynRRTGttAXAHMB34CXjfWrvSGDPEGHNr\n4WYPAtWBl4wxS40xC4tTxJQpcMYZxaxcRETKvVO2xIJ6sBO0xFq0gKefhr59Q1aKiIiEiTJtiYVC\nfLxaYiIiUnxhEWIiIiIloRATERHfUoiJiIhvOQ+xQACWLNHNziIiUnzOQ2zLFu+xxbFTCouIiJyC\n8xBbsQIqVIAqVVxXIiIifuM8xNasgU6dXFchIiJ+5DzEAgFo1Mh1FSIi4kfOQ+zHH9WVKCIiJeM8\nxEAtMRERKRnnIfbhh2qJiYhIyTgPsdxc6NHDdRUiIuJHzkPs9NOhRg3XVYiIiB85DbHMTMjIcFmB\niIj4mdMQe+YZ71EtMRERKQmnIfbllzBgAERHu6xCRET8ylmIWQubN0P//q4qEBERvzM2hNPHG2Ps\n4ePl53tzJm7dqlWdRUTKM2MM1lpTku86a4kdHlavABMRkZJyFmJpaTBpkquji4hIJHAWYlu3wpln\nujq6iIhEAifXxA4dgthYyMqChISQHV5ERMKQ766JrVvnPWrORBERKQ0nIZaSAklJYEqUuyIiIh4n\nIZacDG3auDiyiIhEEichduAAdOzo4sgiIhJJnISYMdCsmYsji4hIJHESYvPmeTN2iIiIlEbIQ8xa\nWL0afvObUB9ZREQiTchDbNcu71EDO0REpLSctMRq1IDTTgv1kUVEJNI4XU9MRESkNEIeYllZv3Qp\nioiIlEbIQ2zcOIhS+09ERIIg5HGyfz/cfHOojyoiIpEo5CEWCEDz5qE+qoiIRKKQh9jo0V5rTERE\npLRCvp4YWHbv1hB7ERHxlGY9MSchFsJDiohImPPVophxcaE+ooiIRKqQh9hZZ4X6iCIiEqlCHmI5\nOaE+ooiIRKqQh1hubqiPKCIikSrkIda/f6iPKCIikSrkIdaxY6iPKCIikSrkIaZ5E0VEJFgUKSIi\n4lsKMRER8S2FmIiI+JZCTEREfEshJiIivhXyEKtcOdRHFBGRSBXyEGvUKNRHFBGRSKXuRBER8S2F\nmIiI+JbWExMREd8KeYhVrx7qI4qISKQqUogZY/oYY1YZY1YbY4afYJvRxpg1xpgfjDHtT7SvqlVL\nWqqIiMivnTLEjDFRwBigN9Aa6G+MaXHMNpcCTay1zYAhwCsnPKCuwhVbcnKy6xJ8Seet5HTuSkbn\nLfSKEimdgTXW2lRrbT7wPnDVMdtcBbwNYK39DqhqjKkV1ErLMf3DKBmdt5LTuSsZnbfQK0qI1QU2\nH/V6S+F7J9sm7TjbiIiIBJU690RExLeMtfbkGxjTBRhlre1T+Pp+wFprnzhqm1eA2dbaDwpfrwIu\ntNamH7Ovkx9MRETKJWutKcn3YoqwzSKgqTGmIbAN6Af0P2abKcBfgQ8KQ2/vsQFWmiJFRESO55Qh\nZq0tMMbcAUzH634cZ61daYwZ4n1sX7PWfmGMucwYsxbIAW4q27JFRESK0J0oIiISrspkYEcwb44u\nT0513owxA4wxKYU/84wxbVzUGW6K8vtWuF0nY0y+MeaaUNYXror477S7MWapMWa5MWZ2qGsMR0X4\nd5pojJlS+LftR2PMYAdlhh1jzDhjTLoxZtlJtil+Llhrg/qDF4xrgYZALPAD0OKYbS4FPi98fg7w\nbbDr8NtPEc9bF6Bq4fM+Om9FO29HbTcL+Ay4xnXdrn+K+PtWFfgJqFv4uobrul3/FPG8jQAeO3zO\ngF1AjOvaXf8A3YD2wLITfF6iXCiLlphuji6ZU543a+231trMwpffonvxoGi/bwB3AhOBjFAWF8aK\nct4GAJOstWkA1tqdIa4xHBXlvFkgofB5ArDLWnsohDWGJWvtPGDPSTYpUS6URYjp5uiSKcp5O9ot\nwNQyrcgfTnnejDF1gN9Za18GNELWU5Tft7OA6saY2caYRcaYgSGrLnwV5byNAVoZY7YCKcD/hag2\nvytRLhRliL2EGWNMD7wRoN1c1+ITzwNHX7tQkBVNDNAB6AnEA98YY76x1q51W1bY6w0stdb2NMY0\nAWYYY9paa7NdFxaJyiLE0oAGR72uV/jesdvUP8U25U1RzhvGmLbAa0Afa+3JmublRVHO29nA+8YY\ng3eN4lJjTL61dkqIagxHRTlvW4Cd1toDwAFjzFygHd41ofKqKOftJuAxAGvtOmPMBqAFsDgkFfpX\niXKhLLoTj9wcbYypgHdz9LF/LKYAg+DIjCDHvTm6nDnleTPGNAAmAQOttesc1BiOTnnerLVnFv40\nxrsudns5DzAo2r/TyUA3Y0y0MSYO72L7yhDXGW6Kct5SgV4Ahdd0zgLWh7TK8GU4cU9IiXIh6C0x\nq5ujS6Qo5w14EKgOvFTYqsi31nZ2V7V7RTxvv/pKyIsMQ0X8d7rKGDMNWAYUAK9Za1c4LNu5Iv6+\nPQq8edRQ8vustbsdlRw2jDETgO7A6caYTcDfgQqUMhd0s7OIiPiWZrEXERHfUoiJiIhvKcRERMS3\nFGIiIuJbCjEREfEthZiIiPiWQkxERHxLISYiIr71/wGRJt27mKFEogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6cfa0f9250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "roc_auc_DF 0.777434408799\n",
      "roc_auc_Bayes 0.884706972838\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr_DF, tpr_DF, _ = roc_curve(cointegratedActual*1.0, -pvalue)\n",
    "roc_auc_DF = auc(fpr_DF, tpr_DF)\n",
    "fpr_Bayes, tpr_Bayes, _ = roc_curve(cointegratedActual*1.0, -logBF)\n",
    "roc_auc_Bayes = auc(fpr_Bayes, tpr_Bayes)\n",
    "\n",
    "plt.figure(figsize=(7,6))\n",
    "#plt.axis('equal')\n",
    "plt.plot(fpr_DF,tpr_DF)\n",
    "plt.plot(fpr_Bayes,tpr_Bayes)\n",
    "plt.legend(['DF','Bayes'])\n",
    "plt.show()\n",
    "\n",
    "print \"roc_auc_DF\",roc_auc_DF\n",
    "print \"roc_auc_Bayes\",roc_auc_Bayes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The above curve shows the efficacy of the classification of the test between cointegrated and non-cointegrated. Perfect classification occurs in the top left of the chart.\n",
    "\n",
    "<a class=\"typeform-share button\" href=\"https://bayesiancointegration.typeform.com/to/gRS0pq\" data-mode=\"popup\" style=\"display:inline-block;text-decoration:none;background-color:#267DDD;color:white;cursor:pointer;font-family:Helvetica,Arial,sans-serif;font-size:20px;line-height:50px;text-align:center;margin:0;height:50px;padding:0px 33px;border-radius:25px;max-width:100%;white-space:nowrap;overflow:hidden;text-overflow:ellipsis;font-weight:bold;-webkit-font-smoothing:antialiased;-moz-osx-font-smoothing:grayscale;\" target=\"_blank\">Contact me</a><script>(function(){var qs,js,q,s,d=document,gi=d.getElementById,ce=d.createElement,gt=d.getElementsByTagName,id=\"typef_orm_share\",b=\"https://s3-eu-west-1.amazonaws.com/share.typeform.com/\";if(!gi.call(d,id)){js=ce.call(d,\"script\");js.id=id;js.src=b+\"share.js\";q=gt.call(d,\"script\")[0];q.parentNode.insertBefore(js,q)}})()</script>"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}