MeanReversionTutorial_2_j.ipynb

meanreversiontutorial_2_j.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "MeanReversionTutorial_2_j.ipynb",
      "version": "0.3.2",
      "provenance": []
    },
    "kernelspec": {
      "name": "python2",
      "display_name": "Python 2"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "[View in Colaboratory](https://colab.research.google.com/gist/gjlr2000/412f445c83afc5de0146da1837283bd0/meanreversiontutorial_2_j.ipynb)"
      ]
    },
    {
      "metadata": {
        "id": "EVVS23TJyeDW",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# Mean Reversion: How to find it?\n",
        "\n",
        "Paraphrasing Peter Drucker's *\"you can't manage what you can't measure\"*:\n",
        "## If you can't test, you can't find it ... \n",
        "\n",
        "Before finding out how to find the signals in the whole financial spectrum, it would help to know what available tests are out there to identify whether a signal is mean reverting or not.\n",
        "\n",
        "In the previous [post](https://www.markov.finance/single-post/2018/03/20/Educational-Series-Mean-Reversion) I introduced the OU and the AR(1) processes. Remember, AR(1) uses assumptions that simplify the math required to solve the equations. While we have not used AR(1) to calibrate the mean reverting parameters of a (suspected) mean reverting series, most of the tests available in the open source universe deal with AR(1) processes, like:\n",
        "* Hurst exponent,\n",
        "* Augmented Dickey-Fuller test,\n",
        "I do not want to rewrite the explanation, as you can find it in [Basics of Statistical Mean Reversion Testing](https://www.quantstart.com/articles/Basics-of-Statistical-Mean-Reversion-Testing) [BSMRT], but just for fun I will apply it to our USD 2-5-10 example instead of the typical Google/Apple pairs.\n",
        "\n",
        "**Warning:** Blindly using statistics can be harmful for your health. Most statistics come with a list of side effects that appear when the 'small print' in the label is ignored. My goal in the following lines is to use them to help us find mean reverting signals, but only as a part of more fundamental analysis. If you use them without an existing background you might end up in the land of [spurious correlations](http://www.tylervigen.com/spurious-correlations) where per capita consumption of chicken correlates with total US crude imports. Statistics tests have strong assumption on the input data which you should be aware (if you were looking for a fully non-human intervention method). For starters, all the tests below assume AR(1) time series - meaning the $\\Delta t$ between data points is constant, but we are using time series with business days where some data points have a $\\Delta t$ of 3 days instead of 1 day.\n",
        "\n",
        "Below, I load the whole term structure os US rates, and then compute the 2-5-10 fly with 50/50 weights: \n",
        "\n"
      ]
    },
    {
      "metadata": {
        "id": "YRS7HOt8ydPG",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "#@title Default title text\n",
        "# -*- coding: utf-8 -*-\n",
        "\n",
        "# In this section I add the initial libraries required for this tutorial\n",
        "# 'requests': allow HTTP-post API calls\n",
        "import requests\n",
        "# 'json':     to format as json calls to API\n",
        "import json \n",
        "# 'pandas':   to plot the timeseries\n",
        "import pandas as pd\n",
        "# 'datetime': to convert different date formats\n",
        "import datetime as dt\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "import time\n",
        "\n",
        "# defining the api-endpoint \n",
        "#API_ENDPOINT = \"https://markovsimulator.azurewebsites.net/api/MarkovCalibrator_01?code=UFYL9CcFRERnVbZJ3ELdRRY7u59vD3BgONzaGLNEOekTfKf0Ug9Bfw==\"\n",
        "API_ENDPOINT = \"https://gmaatest.azurewebsites.net/api/HttpTriggerCSharp1?code=COXmk9ou6Mr8/dUCYj045SaeNm6jCkPzJ3ilU/dCiTvNcFNQax4BMQ==\"\n",
        "\n",
        "# Treasury rates: we can download them from Quandl\n",
        "QUAND_TS = \"https://www.quandl.com/api/v3/datasets/USTREASURY/YIELD.json?api_key=o6YKGjyDcseE3LFKrSK3\"\n",
        "\n",
        "\n",
        "# Call Quandl to get all the data\n",
        "q = requests.get(url = QUAND_TS)\n",
        "q_json = q.json() \n",
        "q_data = q_json['dataset']['data']\n",
        "\n",
        "# Compute derived series from Quandl Data\n",
        "dates = []\n",
        "values = []\n",
        "q_len = len(q_data)\n",
        "column_len = len(q_data[0])\n",
        "column = column_len-1\n",
        "columns = q_json['dataset']['column_names']\n",
        "#print (columns)\n",
        "\n",
        "\n",
        "# The following could be easily done using Quandl's\n",
        "# python API - but I keep it the 'hard way' to allow\n",
        "# compatibility with Google Colaboratory\n",
        "df_rates = pd.DataFrame.from_records(q_data, columns = columns)\n",
        "df_rates['Date'] = pd.to_datetime(df_rates['Date'])\n",
        "df_rates.set_index('Date', inplace=True)\n",
        "#df_rates.plot()\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "ByC83D4E3HwP",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 353
        },
        "outputId": "fc113e3f-ea09-4da9-8891-cca12856d118"
      },
      "cell_type": "code",
      "source": [
        "\n",
        "\n",
        "\n",
        "# columns names: [\"Date\",\"1 MO\",\"3 MO\",\"6 MO\",\"1 YR\",\"2 YR\",\"3 YR\",\"5 YR\",\"7 YR\",\"10 YR\",\"20 YR\",\"30 YR\"]\n",
        "# as the index start from 0, columns 2,5 and 10 are:\n",
        "cols = [5, 7, 9]\n",
        "\n",
        "columns_fly = [columns[i] for i in cols]\n",
        "TimeSeries = df_rates[columns_fly]\n",
        "\n",
        "# I use this line to keep each weight linked to the name of the maturity\n",
        "Weights = pd.DataFrame([-1, 2, -1], index = TimeSeries.columns)\n",
        "print (Weights)\n",
        "\n",
        "# Now all I have to is the 'dot' product of Weights and TimeSeries\n",
        "butterfly_ts = (TimeSeries.dot(Weights)).dropna()\n",
        "butterfly_ts.plot()\n"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "       0\n",
            "2 YR  -1\n",
            "5 YR   2\n",
            "10 YR -1\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7fda2f273590>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 2
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXecFOX9+N+7VzjuuDsOOHpvDwgI\nCCIIKog1scaaRBMLsURNTIyGr/5SjBo1iTUaYzcmtmjUxI4KKEqRJp2H3ttxtDuOA+5uf3/M7N7s\n7OzuzOxsub3n/XrxYm/qZ2dnPvN5Ps+n+AKBAAqFQqHIPvzpFkChUCgUyUEpeIVCochSlIJXKBSK\nLEUpeIVCochSlIJXKBSKLEUpeIVCochSctMtQJCKiqqIeM2yskL27q1JhziOaCpygpI1WShZk4OS\nNT7l5cW+aOsy2oLPzc1Jtwi2aCpygpI1WShZk4OSNTEyWsErFAqFwj1KwSsUCkWWohS8QqFQZClK\nwSsUCkWWohS8QqFQZCkJhUkKIR4BRgMB4OdSyrmGdTcBVwD1wDwp5a2JnEuhUCgUznBtwQshTgH6\nSSnHANcCjxvWlQC3AydJKccBxwghRicqbLp4/fPV/O2dJZbrVm3ex0/++Ck79jSNWF2FQtF8SMRF\nMxF4F0BKuQIo0xU7wBH9XyshRC5QCOxJRNB0MmXuZubJCh7+97cR6557fzk7Kmt4f+aG1AumUCiy\niscff4jrr7+aG264hhUrliV8vEQUfEegwvB3hb4MKWUtcDewDtgIzJFSrkrgXBnB0nWR7yi/T0si\na2hQjVMUCoV7Fi6cz5Ytm3n66ReZPPk3PProXxI+ppelCkLpsrolfyfQHzgATBVCDJVSLoq2c1lZ\noWUmWHl5sYciJo5Znrw87R2Zl5+bcbJGo6nICUrWZKFkTQ6JyLpixSK+852zKC8vprz8WGpqqmnZ\n0kerVq1cHzMRBb8N3WLX6Qxs1z8PBNZJKXcDCCFmACOAqAreqoZDeXkxFRVVCYjoPWZ5gpb7oUNH\nMk5WKzLxmkZDyZoclKzJIVFZN2/eTrduvUPHKC4uRcoNdO/eI+55o5GIgp+C5oZ5WghxHLBNShn8\ndhuAgUKIllLKQ8BI4MMEzpWx+P26i0Z5aBSKrODfU9cwd+Uux/vl5Pior7dWBMcPaM+lp/Z1dDwv\n+mW7VvBSyplCiPlCiJlAA3CTEOIqYL+U8h0hxJ+BaUKIOmCmlHJGwtJmIMHfQPngFQpFIrRr147K\nysrQ37t376Zdu3YJHTMhH7yUcrJp0SLDuqeBpxM5flNg2+6DAKzdtj/NkigUCi+49NS+jq1tSNxF\nM2rUaJ5//mkuuOAipFxJu3btKCwscn08yKB68JnK0nWV8TcCqmqOJlkSRSwCgQANgQA5fpWcrWia\nDBkyFCEGcsMN1+Dz+fjlL3+d8DGVgo/Dw/8Onxdeu20/fTqXAvD1ku1WuyjSwGNvLWbx2kqe//UE\nfL6o/Q8Uiozmxhtv8fR4ytxxyOxlO0Of//f1+tDnEaI8HeIodBav1UZaR442pFkShSJzUAreITl+\nn+Fz4+XLz1WXMhN4Y9qadIugUGQMSivFwCpMyW8Y/ufmNH6evXxnxLaK1DN94VaueWAqMxZtS7co\nCkXaUQo+BlWHIidOfYYr5jdY8x6ErCo85MWPVqZbBIUi7TRrBX/kaD2ffLOJAzVHLNfvPXA45v6b\ndlaH/b17/yHPZFPEZ+feGj6fv8WThBCFIhtp1lE0H8/ZxLtfreeNqWv4/dXH071DeMrv/a/Md3S8\nepXslFJ++/w3HK1roGt5YrHCCkW20qwt+J17Gy3u3784N2K9VUTGQYPbpm1Ji7B1KjwvtRyt034f\nlYOgUFjTrBW8m6H9l4saY9/7dWsdti4vp1lfzrRR16BCIxUKK5q1RorlUrFTW8YYEw/eFAdSOKf2\ncH3Esr5dStMgiUKRWTRrBR9Lif/m+TnOj6cUfFp4+RMZscyYr6BQNFeat4KPoZC3VzrvsarmWDMH\nv1LwCoX7KBohxCPAaCAA/FxKOdewrhvwGpAPLJBS3pCooMnAa4NbuWgyB2XB22fG4m18Pn8Ld105\ngjyLrmqKposrC14IcQrQT0o5BrgWeNy0yUPAQ1LKUUC9EKJ7YmImBydBL+ZaM5t2RpYFVTXhMwdl\nwdvnxQ9XsmlnNWu2Hki3KAqPceuimQi8CyClXAGU6X1YEUL4gZOA/+nrb5JSbvJAVs9ZuHq37W0v\nP7Vf2N+vfbY6YhtlwGcOfhchq598s4lZy3YkQZqmgXolZh9uFXxHoMLwdwWN/VnLgSrgESHEV0KI\n+xOQL2NoW1oQ9neFIWt1zJBOgHLRZBJuUhLemLqGZ99b7r0wTQSVxpF9eJXJ6jN97gI8htab9QMh\nxHellB/EOkBZWSG5Fv6/VHZUb9u2VdShfXl5MV3KW1FVc4Si4gL2GMoYBK3FOp+/SXSAbwoyBnEr\na35+rut9U71fOrCStaysKCO/QybKFI1Mk9Wtgt9Go8UO0BkIZgDtBjZKKdcCCCE+BwYBMRX83r2R\nUSup7qi+q+JA1I5AFRVV7K8+TPWhozz++oKwdcvXa7XI735uNi9MPtXWufZWHWb3/kP069o6/sYe\n0ly61NfWHnW9r5v9suG67ttXQ0VFfhokik42XNdUnDcabl00U4CLAYQQxwHbpJRVAFLKOmCdECLo\ntB4BRAYqZyB1dbFdLNV6mYKl6/eELQ+mzDvhV09+zf3/WkBNrUqzTwb1Dt1llftrkyRJ08HNvIUi\ns3Gl4KWUM4H5QoiZaBE0NwkhrhJCXKhvcivwor5+P/CeJ9ImmRc+XGG5vLBF+EDHXKOmtFV4TRo7\nBNVPXb3y23uNDzh8JDK7NRa3PzUzOcJkODW1dY1/KP2edbj2wUspJ5sWLTKsWwOMc3vsdDF35S5u\ntFj+h2tHhf196HBd2N8t8tzHDm/YUcWxfdq63r+5EmzRF2T8sM5M/1Zr8lHQIpdahwq+uXJQjSCz\nmmabyRqMeOnfLb4PvE2JFkFT3rrAcv26bftdy7Fi4574GykimCd3hf1tHEXl5fioq1cFyOxgHD8q\nD0320WwVfLDQmJOMx8G9vLe0ld/THbWmUZTxZ8zJ8VPvwPWlaghp+JSPJutQCj7H/k1t9TK46cIh\nYX9/MGtD3OMY4+VVxqU7VmzcG77AZ2yG7nNUQlhlIGsEyP7rsHB1BQ+9vpAjR5uHC6/5KnjdwsuN\nEhZpRV5e5LY9OrYiP7dx+X++WBf3OMaJVdUkxB2HTQ+o8Srm+H1heQrxUApe47n3lme9a+uv/1nC\nsg17WbCqIv7GWUDzVfC6hWe2yldt3hd1n5b5kXPSOX4/w/q3d3Tu8IdIKRenvPf1+ojoI+Nfu/WQ\nR7tWmnLRaOzce4h5K3fF3zALONRMJuGbrYIPRlnUmHy5D7yywGpzAPItomX8fp8jNw/AUYOCf3/m\nRkf7Nnd27KnhnRnrI5b36qQle0wY3iVUVsJuj1w73pz9B4/w6dzNWW/hbnNRJrsp0lxGbc1WwX80\nR6t/FuHLNTG4V5vQZytvSo7fR66pVV9VzZGYx3QyAagIJ5rSbl3Ugqd/NZ4rzxTs0nvtvjdzA3sO\nxE9gsmPB3/uPubz2+WqmLdzqTOAmxvszN6RbhJRg9+Xf1Gm2Cr5iX2OxsNKi6OnZRh+5lZ3eIs8f\n4ea585nZMc99NMutwGSSG2VSuiEQIC83/Hb+eM4m7nw29m8B9qy5St2nLzdFd+E1SZqBe6q+oYEp\n32yi9khd2LLmQLNV8HWG8gK/+fHIqNsZrXarCdHcHH9EJMzB2rqI7aKdW+EMp408zFnHVpituVgK\n32gYZAPNQL/z4CsLeX3qGn768JehZeZkxWyl2Sp4aZhMLSyIntBrjFOfuXR7xHqfz+c4lj3b/bjJ\nJFpYaUEL94VRzWWeYw3fN++qdn2eTMTKPZVtZa/XbI1MRDx4SCn4ZoPZMjdacEaFsn67daU4p7Hs\nxgiQY3qWOdpXYU371i1tb2tWYGaLvTlF1Vh91cfeWpx6QZJEsECgmeYSndzsFXyfLiWY9XO1oT6H\nnRsh1gjACqMFv3xD7EleRTiJ6t4/v7aQu1+cG7bMXHkynk8+myxcq+9irvPTlNleedByedfyVimW\nJD00ewW/v/oI5ulTo4/caN3/8PT+lse49DTr5W9NX8vjby2OeIjUJKt7Es22XLFxL5tMbhazQo+n\nv7PJws/2YJJoiYQqiqaZ4Pf7yDXFsRt//L1VjWF2Hcqs3QDFhdZROB/O3si3a3ZHKITmEoObLob3\na+doe/Pv8faXa/l4TvQ2wtkUgJHtk43RRuAbd1Tx2JuL2H8wdkhzU8e1ghdCPCKEmCWEmCmEOD7K\nNvcLIaa7li4F5Ob48fl8jB/WObTMqJDXGjrNO4ngMMbCmxVCME47yF//s5g/vDQ3qyzDpGHjEkUb\naUXD/L6dumAr/562hgdeWWA5xM+m3ynb4/qnfLPZcvlXS7azaG0lb3y+OsUSpRZXCl4IcQrQT0o5\nBrgWremHeZtjgJMTEy/5BK13Y5aq0aJraYjOyDElNJ0+slvU405d0PjgmBXC66abauHq3WzYUZV9\nMdZJwI5qbVNSwHljewLQtiR+M5ZoI6pVm/dZTjhm0wgsWsjusg3ZUcZ6bpzSC7OX70yRJOnBrQU/\nEXgXQEq5AigTQpSYtnkIuCsB2ZKG0Sd+bB9tOJ9vKCRmdNF0a984GWO24M/VlYgVWyoa/bxGhVBX\n3xDV/6fa98XH6sqdPLRTxLILTupNi7wcigry4h4zlkV+wGIIn1WTrFGWP/T6tymVQ5Ec3AYPdwTm\nG/6u0JcdABBCXAV8AWywe8CyskJycyNrvSSjS/mbn68KfR4+oAPl5cW0Lm30r7duXRj6nJeXE5Jh\nnymBqby8mFYtNQXyi+8P55HXFobWzZeN1eratG0V2u7tadGHhK98tpqzxvVx85UckWmd32NhlrXO\nF2mTjBrc2fI75ef5we+zXLe7+igD9TIUlTXRX6y5Of6I/cvatKLEIvu5KV7X/Pzo3cj2HqqjrLiA\n8ihzT6nC7XVdvdlehJqXv1um3QPus0PCCZm2Qog2wNXAaUAXuwfYuzeyyFGyupS/bOi9mkuAiooq\n2hk6Au2ubLS+jx6tD8mwb1+4jHsqqznUIpfy8mKqq6KXp50xbxMjB2gVJ1fHqH2zr+pw0ruyN/Uu\n9XtM98ktFw1hYNcSy+/k8/k4cqTect0dT8zghcmnAlAZJZQOtEYiO3cdCFtWUVHF4ZpwBd9Ur2tZ\nq+hlOm57TMv8DF6ndJDIdd2wxZ7L06vfLV33QKyXilsXzTY0iz1IZyCY5nkqUA7MAN4BjhNCPOLy\nPEkn6IIZ0luz5kpb5Ye5ULqWF4U+m10rxgzW1sXRfb0V+xsnVb9aEpkNa0T1yIyN2aXQrX2rqKFw\nOX6frZojsVwuNYfrmbsi3I+bTSF2A7pnb6Jdlc0ImUAgkLXuUbcKfgpwMYAQ4jhgm5SyCkBK+ZaU\n8hgp5WjgQmCBlPIXnkibRHw+H21LWpDr94d85qVF+Vwyvm9oG3MVSKNeGdC9NT855xjLY7d0kEZ/\ny6MzWL/9QPwNmysm3ZoTo2GLpuDjK+NY29TVN/D0/5aFLTtalz21xLNoOiGMPQdqefGjlXG3O6Zn\nGS99tJKbH53BTgsvQlPHlYKXUs4E5gshZqJF0NwkhLhKCHGhp9KlmBy/n/qGhpCCP3loZ1oYfJTm\nGjLGEgU+n48xgztiRQc9jd4cPRONe/4xj1nLdjiSvblg1kexykT4bSp4p2GPtVnVLCL+dz/aBIvj\n2a0ZlJfjZ8ZibVS9bmv2GVauffBSysmmRYssttkAjHd7jmQzYXj4FEFOjo99+47w1/8s0f42KQ9z\nFyErz8DPLjqWx/8THlrXENAiaabMtY7JteLZ95YzZpD1C6O50RAIcPeLcxnaty1D+4YnMcXKTcjx\n+2LW3v/bu0sBGDckMgonFmu3HaB7h8yaTHOLHW9TTe1RSlvFDzfNFN7+cp3tuvZGA+HZ95fTq3MJ\nHdsUxtijadEsM1mDPVQvPLl32PLgcD/Y5clnUh71ZgveQsOL7q0jltU3BCybgJx1QncHUjdfag/X\ns3lXNe/P3Mi3q3eHrYut4P1hIaovfyLD1s9buYt5K3c5tuDfnLbG0faZjK2QzyZWmctJ0xLzMzzl\nm+gZzE2RZqngj+hDzgJTiJi59d5BUyU683Dfuj585LKGhkDEvoUtchnr0HLMFvYcqOXZ95bZ6rYE\nYHSzD+hRZloXx4I3XPfp0bI2ber3NnrS1GkxEtyaHDa+ezYldpnJ9tLdzVLBBzErh407wkOczC4V\nO5ae1aRffUOAfdXhFvwdPxjuuHlFtvDKp6uYtWwn/zRZ1HYwJ+DEtOBzbEbR2NTwQV+03RdTU8CO\n7t60s4pJD05jrUVd9UzDaRTaIg8qZzYEArzz5To2ZmBwRPNW8A6Hnn06l8Y/poXC+XLRNt6aHj6s\n12rgODp91rBQd7PY7Wwf671qd5I1pivCpoEaLCo3c2n2TIDbcdE89tZiGgIB7vvn/Ljbppu/xqll\nf/64Xp6fc+m6Pbw3cwM3/2Wa58dOlGat4J3StrSAv992iuP9lqyrjJigzc1x3gkq2zh8NPFolFjX\nMNfvIxDQLCwP9HvUaqKKzGHVlshRxmA9xwXiFwzcVuk8VNKL+zhZNFsFb+dhbWGRxm0sSuYEcyx8\nbo4/q6oSusFuFya3lyn4MDc0BGJe68/m2Ytu6tnJXG6p6dMc7sHWRY0RQOZ5NjOrNjsv+Gd8Z+ze\nf4jK/Znjwmu2Cr7YopaIGdEtMiIG4JFbxvHgDWMcnW/JunBfX06Onzy9OqVxYnaAIQpnz4HarCps\nBeHxyfa/m7tr4NfnQ+rrA2yPYZlt3GkvvTzYXrF35+xR9Fl2e1lifL5ybIyaDzvMczC6Ce94aha3\nPzXT0f7JpFkp+EAgwBp9CLfGYihnJhhOaaa0KJ9yBz1Arcjx+2hTUsCtlxzLA9c3viyMVuKv/jaT\nr5dkj78XYOeeRkVrV7ms2+Zu8ipowdc3NETU4DeSrxe569s19hxLrv7CcCtPJuLUgNi9P/p1PFrX\nkJEGiTGu3Vzy24r9FiHNschkV2uzUvBP/28Zf/yX/Ymivl2tLfh4/OziY+NuE7Qqju3TjjYlBaFy\nxQ0NAdqVFoS2i1e7pqlhtHbsugdW2SwaZSY4HI+XzRpcW9wydmnheMP7pohTdbwzyovyaF0D1/9l\nOo++mXkNu08Z1pjQWGkjAsqc7xKPWBP96aZZKfhvVsQu/n+BaYbdbRjjsL7xW8aZj90yX/PRFxfm\nhUXXZJuP1Gjt2P1quRahpyMHtOfqswfE3C94/Nc/XxPztwyuMVp3owa2j9jOeIw5WdIowqnFHW1U\nG0wONLsiM4Hc3MbfbWWMaq5BnIb9WxnwmfLcNisFH49YFSG9xhwvf9tlwzh5aCfOOL5bWAKVHVdS\nU6LUUJ722zW7Y2zZyMCekRUPR4pyThra2WLrRoLx6rOW7aBDm+gutUF6XXijAr98Yj/OPqE7JYWN\nVr3RUjMXIGuqONVD+RY9G8C9MZQKjEaFHW9KfX0DD73xLV8u2mbr+FapFpmSHKYUfAyS6VozH7tr\n+1ZcdfZA8nJzYvqLmyNWPk47RcSMCrkqRlOPoJIznqV1qxZcMqEvow31gOxM0DU1nCp4q0ztTMfn\n83HWKK0syE0XDom7/fbKGpat38NLNqpRQuZY61Y0WwX/3TE9IpYN6xfuWknkVj7j+G50aVcUdX20\nGuZesOdALf+aIqk+lHk1rt08DFZ7DOndNu5+RgU/O0Z1zmC2a8+OxbQpacGlExpLRBuPkcm+VreY\nXTQ3nD8o5vbR7lvj7/rc+8s5dLjOcrtUExwxXnpqX16YfCptSgri7OG8eqaVtR6ryF0qca3ghRCP\nCCFmCSFmCiGON62bIISYLYT4WgjxghAi414kVgqipDCfFi7j3M1cPrEf90w6wZNjOeXZ95YzdcFW\n3vlyXVrOHwuzfv9ozsb4St9idas4E6IQbvn36Bi9+mPwAc3Py+EvPx0bVgQufHifhQre9PeogR1i\nbh9NcRt/wplLd/DxnPQU7TLnm9xzrfNn8KhhktVOxJSVgo9X7nvpukrmxWkI7gWuFK8Q4hSgn5Ry\nDHAtWk14I88AF0spxwLFwFkJSZkEotWsMFpp/VxG0cTCSfMPt+yt1toH7tgTOytvzvKd7K06zOxl\nO9hfHb3loJeYLcY3p62NqBBpxu0Q2GhwW/VQbZRJ+99KfxunSrIyisbhtX07itFgPs5iD2q8uMEs\nh5VLqVen2KWejRb8vS/Pi3vOTbsi8yjMlUvNPPzvRaFy1cnErWU9EXgXQEq5AigTQhizP0ZIKbfo\nnyuA+OPpVBPlvjZWl+uqt/NLBHPI5KO3jOPJX5yc8HGtCAQCrNiwJ+TDXxEjYuD5D5bz9P+WcduT\nX/PMe8v5s6mIV7Kw0ifxhvMzFtub7DJjfFnH0mPBF4iVhW604DN5ItEtxuvSvUP8+31fFEPAfH2T\n7R6sPVLHio17IxS62Riwmr85rn95xLIHrh8d+mzu2BWvuNz7MzfGlddIKttyulXwHdEUd5AKDD1a\npZQHAIQQnYAzgA/dCphqvO5eU2iy2PNy/Y6s+LwoYWlWzFmx05aiXr1lX0QC1bbd0RtPe4mVxRjP\nMt65x92ks/GosSzVmMXMjC4afKGIm0ynrr4hLKksGsHrcv64XvzfD0cAsUeZeVEShczX104Vz0T4\n2ztL+fNrCyPCMs2ntXppnz26B3dfMyr09xVn9A8b4Zl1wMI4I0wrcmMkVD34yoLQ54O1Ry17RXiF\nV/6CiKsohGgPvAf8VEoZd7xWVlZIrkUIVqyO4YnQvUtry2P36FgcKhvs5NzRtt1kSJH/zTUnOP4+\nR+sabO+za3+kJVF9tIFepiqYNUetH75kXWsj2y3qdLQuLYx67vLyYnIsXnJ2ZM3Lb7y9i4ujh0nm\n69uVFBdEHLeVIXS2fftifjtpNJfd9SH9ukXeP6m4fna578U5zF66g0duPYW+FiU3grIW6XVaRK+2\ndO2ibffE7RO49t5PLY/bsmWe5fcM5IQ/uzW1dZ5dD6vjLF2/B4D9h8znCX/RtG9fbKlsO7RvdDh0\n71xKty5lnDCoI3OW7SAnz6QWc/xRv8uOSmvD6LyTekfdZ0tF4z63PDoDgPceOt9y20Rxq+C3YbDY\ngc5AKOVSd9d8BNwlpZxi54B7LRrelpcXU1Fhr05INBoCActhWtuiPMtjT/7Bcdz48BcAts8dS86K\n3Y21V3q1L3L1fTZu3kthQfyf6tM5kQp+wfIdFOX6QpZMeXkx1VGG2YleayvM19/qd/7TP+cxoEtk\nfZfgdS2wUPB2ZD1ypNH1s29/dGu29rA2ZK6qro047iyDe6hyd3WoAN3qzfvCtvXiXvWS2XpJ40Vy\nJ6UF4crXKOv+Ku2FW13V+N1jZXIW5udYfs/d+8JHWUfqGjy5HvGua3X14bD1xvDZs07ozt49Nkam\n9ZqsPdq3Ys4y2Gsqx9Amiq4AkBv2WMt18LCj75/ItYr1InXropkCXAwghDgO2CalNEr4EPCIlPJj\nl8f3hLe/XMekB6dx4OAR25NJVhUkE2Fwr8SnH5autzdhdbA20pf90kcrud8wJAR4/oMVlvs7LbIU\nj9c/X82kB6dRa1C0buZLozUzd0LMcsGhOPhIQ2CtoRGzRUJtkyegK0S74aDH9rG+nzOlBo1RDGO4\nqxW/v/p4Lp3Qlz568bjg9z58xPQcxfhq0VwxTvovg/2Kpk5xZcFLKWcKIeYLIWYCDcBNQoirgP3A\nJ8CPgH5CiEn6Lq9KKZ/xQmAnBHszrt6yn+H94pcPSAZeRF78/b/LKG6Zx8Ce7vy/drNhb3z4C354\nen8mjujq6jxBvlmxk7emr2W37o7ZtLOa/rqbwI0iMPpR+3ct5bKJ/RwfI9Z5g37cmrix201vkjWe\nxPUBCwXvIhw0Q/S7I7p3KA5rnh6cRD9scmHGeoa9iqx69bPVSWkF6doHL6WcbFq0yPA5o1qw+3zp\nyzbzKjnmz69/ywuTT/XkWLF45dNVCSv4v/83PI3fGFc8LUpf1MVrK5m6YAs3XTiYPNNcjDHO+Moz\nBV3K7UU3GRsxbLUxibxwVQVnHB/9ITPrvZ17auhgqFT4rymSqQu2ctOFgxkhImvZpAWf9t1f+2wV\nV509gHal4XMRwWtrV6n/42PJyUM7R0xeZko2Z3AO7ZaL4mesmmlU8OEj2VhZ0wX5zlSouS1ossnC\nQac16RpCZkopUa+jg5zgA96ctoaFqyuiRiQ8+uYiFq+tZNGaSHeUUXk4STa68gwR+vzBrPihbE5f\nxv/3zOywv6cu0F5eT76T/Phmu+zYU8NvnpvD8g17efXT1RHrg7rLifvJqqFFhuh3fGgF0Yb3iwyF\njEejiyZcwccKCKpz+Fw99tai+Bt5SPKzblJM5f5apn+7lXNO7Bla5iN82BXPEvbaUi4pyqdfl/j9\nXJPJ1AVbOHNU9/gbAvurD1PayrtBWFXNUT6aswnmxN/2SF3kPMArU1aFPjtRwt07FDOoVxuWrbee\nCDMTrVJiEHMcfIHH8zXJ4KPZjRmlVsXdgha8kzo7hy2UWroMqODLf2/VYaYu2KL1+XVpUwWNoDWm\n5uKxRidzVjirKrqv2jok8qRjOzk6jl2yTsE/9tZitlRUh8WfVx86Sm0aa2M8cvPYtKe5vzF1DafH\ncD8Y8boQnpOsxhqLiWKjb3zHnpqwBg7xcJKcFC/nwPwb1no8KZ0OGiwmWePx9hdrueWi8AS+dFnw\nwfM+879lSBft9ozMXGrdeyFWTH83QzKk6F6G3NSYXLh5V3XY+lgkK8M961w0Wyq0sERjRMmLH61k\nZpzaEMkkVcrdWIrXillLd/DEm/ETobwudbpgdUX8jXTMHXc27AivBeK0k5YTF5mVgu9n0eUpXnRG\nU6LBYpI1HmY325Gj9VEjs5JNcORQEaPTlF36dLYeZQdieGHa6HkS557Ykz/dchJ/+emJoXW/e+Eb\nFqyqiMiMtSJZWdJZp+CDmGuJTl/wAAAgAElEQVRQvDtjfZok8Y6KfYeoq29g/fYDlkPi/VGGf0Ge\n/2AFn8yO74uuc9jRJh6JhF/+4aXwWiAdY9R1t8KJ4rJS8FaTrk25Jk1NbR0bDS9Np5OsVvzqbzNt\n97X1muBjsOdA4rWUohVaizXJWhd0ceX48OttOI088fYS3v5yHQ0NAdZvj164LFbmayJkrYKPNtw2\nNrXOVKL92L/++yzemr6We/4xj68WJ97K75FbxvH7q4+PWF6dwloZZoyp8FYTw04VkRPDqIU5gxHr\n0Ve5IRLF6J8d3FsLY22fYL/eRKk1x3Eb+M3zc7j5z9NCtWLcWPBm0lmW2svonaNRDJtY5wjNYcS4\nfsvW72HK3M3c849GY2X8sPBmNZ1jlBZPhKxR8GYrMVrNjJWbEvPTpYKu5dqPXWSRvbpQd3es2BS/\n9Vg8Sovyw+KAg9z3snXf2kAgENe6T3Sy7d/T1oQ+W53LqbvLieKy6vpktffQvo3JPs+9vzz0OXjP\n5TqoH5QMYs0N7K3SLN09B2oJBAKhCJFoL84TjoldPjjdNAS8axNorhvVUy8xHctlGaz7bu7QZqSo\nIC8iWfH8cb34041jmHTOQG44f5Bli0gvyIpJ1tc/X82UuZspM9QN+dCGKyJTueWiY/lmxU78fh+v\nfRYe2ha07r1sKDB6UAdmL4sfDfDUf5cxb+UunrrtlKh183ftS8wX6rU1aFbww/u1Y2vFQUs5T7Zo\nAWgVKWN8ycxetpPrztWaZATfbenO6rQzyvn9i3MZ3q9d6JkxW6D3XDuKtdsOcOLgjhnefzbAfOlN\nXXVzz4DQsxZLwetvyJg9f32NL9YgwSg1c16C12SFBR9MCzZexAMxWrRlOmXFLThzVHcmDO8SsS54\nI8WzpC8Z38f2+Ypbxp6c3bijimsfmBpqULCvKrq/8yMPX6xeKErzi+v0kd04cYh16QMr15joUcZ3\nx/Tgt1eNjHuuYBji9sr4VRyTis1By8LVuxtdDKZ5hS7lrTh5aGdyc/xcOqFvaFRpF7M17BXrtx/g\n2gemhv5evKaSAweT86wH5/FiuWjqbUQhrdy0L233RFYo+GzFSuEEh4ILV++O6Yc3diWywlhzJ2BR\nbKOhIcDho/X84omvuPulueFbxFAgXy6KLpOxgfUpw2I3zA7K4DV+v4/TR3aLGBK3jhKB5Pf5uOiU\nPvTsGFkMLch9/5zHgYPhE9zpTCz7u4NGEnU2FNRZJ3SPWSZZWFSrTFZm6wezNobdi5t2VTsqqe2E\n4PMX00UT5QWZKSgF38QwPogvfNgYmma0di8e3yfMjWBVICpeZMukP01j2oKtlpE5z72/3FboVwQG\nmX54en8uPzV2uGEyGtP7/T5atsjlhvMHhw2rEwlTW7v1AJ/M3hC2zCphy8iND3/BL5/4yvU5oxEI\nBBzNM9mZJARCTautaF0cmRRXe6Q+KUreSspkKfjgSG9rxUGefHsJ8+Uunnp3adi97yZRLJUoBZ/h\nmDPczKFW82UFb3+5Lsw9FbROTx7amW7tWzHewtVjnK+IVi3POOFpZO3WA3zxbXiXpZraOl79dJXl\n9kGMVm5ujp/j4/T/9Lq6JYT7py86pdGNVZlgmJ1ZQcby24L23aJlNSaCU50askDj1CqIldkcTbUl\npe6KxcnMFWC9CjlsqdeZ+WrJduavquDJd5Yyd+UuPpu3JbSNVaLYHd8f7sn5vaDZKfimECZp5Orv\nDOTiGP70J99ZwvszN4SlV7fWH8arzh7A3deMsgzbM/pIRw/SLBUnVuyrpsnfD2dv5LP5W6JsbU2Z\nheVn5I3PYve1tIN5wtGYH9HXIonJ9Xn8vrBSB4FkDD9s4NRqrrdpwWcKZik7lLWMyBZNJI9jhGis\nYWMuWRDkzelrQ5+DV9t4mw3oUeb6/F7jWsELIR4RQswSQswUQhxvWneaEOIbff1vEhfTO+L5pjMR\nO1ER/4pRr8VcHQ/CH4LenUt45vbx/PySYyO2s8s3DmtyxCLobrKTlBUPc3xxUUHjPICXKu3I0YYw\nqz2eBe+GQ4freHfGOhauip4Z7PS8weYeyfAhJ8UNb3oWRg5oH7K0vcA412BHUQfvVateAtGwmrNI\nFq4UvBDiFKCflHIMcC3wuGmTx4GLgLHAGUKIYxKS0kOS3CoyKdhx7xnDC+28EIabGg/n5vgZ3Kut\nq4zGhkAgVPvdC7xUjgcOhrteClo0DuftlBC2yz8/WhEmdzImiP81RfK/rzfw17eXWNbscXPeDbob\nxY5bI9qIJ3oj7iRMkptuz/qGgKfnWWloVB9rtP/RnI288umqxnkiB4/N7Sl04bh99U0E3gWQUq4Q\nQpQJIUqklAeEEL2BPVLKzQBCiA/17ZdHP1zyyMv1h0U0ZErdaickWsumpLAxQuTUkd04aXBHura3\nDntzc328fpAbGgJgEWb/5xtPjFwYh75dW7PAYPEaLa1ECtDl5/k5EqW3LTQ20rDCrQthtaFxS/Wh\nI5ZtHJ3+fsG5m2iJgUZyDdp1975DtNNdf9EmdVPxrDU0BDwdKbQvayxkl+P3k+P3WRocb07T3DTt\nSrXSBE6eUK96RNg6l8v9OgLGcWIFjT1azet2AcmphWkDc7haMiyrZJPoBH3b0gJuu2wYV509gJ9e\nPJQeHYvjTqrFw5gAlMgDZlXMK5oF37a0wHJ5LM42ueSM17KsxPnxgsTzWZvvs0AgEOpb6jaE0o4L\nyO397VTp2GlJ51WGqZHNu6rD/m5oCHj6IrlgXC/OOL4bD94wBiCsoYsVoZGr6fJZuYJPG9GVu68Z\n5YmcdvHKeRXr7rB155SVFZKbG2m2edGZ/b4bT+Sup2YC0LZNkWfd3o0k45hBSortK6KObQstZRlv\nWOaFrGeO7hk6jtOQybt/Mia07w/OGsjdz4U3zihrU0RxYWRcuhu5q01Wdrt2rSjU/fCnlxXylCFm\n3MnxNZdG9O+9oeIgxw7oSCAQwOfzce5t/wW0OYEbDaV2nZzTGClVUmr9O//9eRtF9y2wI0fn9sUh\na33J+j1x99m066An95rxGOaEoRYFebQyRfj88caxCZ33lsuPC33Oj5Kxbaa0RBvNBM9bYJjrCXJM\nn3YcNyi1tq5bBb+NRosdoDOwPcq6LvqymOzdG5np5bZT/clDO4Ul3HQyWH5Hao940u3diFs57XLw\noP0Qvl4dY8vilay1hut4xGISNxZlhbmhfWssvtuuXVXUFoUreJ/PXef5/aYyspWV1RzUJ+XMlp+T\n45cU5VMVI1v62XeXIjqX8Ku/zQwrL7xt90GmG6zf6d9sjJlEFI2K3VUU5UbaTt8sd1cW2853v2Bs\nT6bqzaF37qmJu8/iNbsTvtfi3a81h46w/0D4/E+grs6z53FDjAqQRqqqNBmC5/VZTPZt3XkgKXoi\n1svM7Th9CnAxgBDiOGCblLIKQEq5ASgRQvQUQuQC5+jbJwUr6/Gqswdymqmv6K9/MJyJx3WlXwpn\nsL3CSQhbnYc1aoL0t8pUNNy/dobI91zbODQ1hiqayzpDo/vhIoNi/O2PI6te2iGyC1OjTZNIidyb\nL4zf8zPY6MScT1BnuHgfzNrg6vzrt0cqij+9usDVsezSqmWkVZp2ApFzQE4awniFeZ7sVF3/nHNi\nj9Cy5RsSLxDoFFcKXko5E5gvhJiJFjFzkxDiKiHEhfomNwKvATOAN6SUsTNgEmBDlGSKC07qzZDe\nbfm/K7Thluhexg/P6J8xPVKdEIxTt0Oifs87vj+cE47pQK9Ojan5l0/sS1dTo+s1W/fx8L+/pab2\nqGVk0vB+7cL+7tS2cVLX6P83nifIbU9+zbSFWynSFcovLx0aUQjKLkbf8pDekRm9bunQppBrvjMw\n5jYvf2Idx2/s47ly0z5mLbVndY8Z1JgY9k+LYzeFSqnJwGxfpKN7mvmMRQV5vDD5VL53cmMOSzpU\nj2sfvJRysmnRIsO6L4Exbo/thDVbrJMRCgty+cWlQ1MhQtIxV24Mdo63ItE2cgN6lDGgRxnXGAo6\nrdmynz9cO4r5soIn31kCNFqQny/YGlEU7acXDKZ35xI27qxiz4HD/PgsEaZojTHX0UYn//xEcsXZ\nA7Q/Engwcj0qR2DFuGM7hZWLsEvPTiV8bVDqz76/nDGD47/Ejd2uzNe8KUaHeUW6q3faJVlNPWLR\n5DNZjVll2UwXQzW/VEcCBWt9jBDl3DvphLB1DQ2BCHmG929Hm5IC/vLTsbww+VROGRaujIyjqJjW\nln5YJ0kkZvxJVPBuidfcOxrG8Mode2q446mZfLlIm95qitFhXhAgOTWLnGLHOjdmyaaKJq/gmwt3\nXjEi9LmkKHZ5X68xRhKYMx4DgQCzTfXCvXKDhZ7bBA4XbeTgFd072GuqbMQqDt7K5RKxn8G1s2Lj\nXnbvr+Wlj1YCTceKTQaZ8N23VsRPmmvjIBrOK5SCbyK0bJFLeWvtBsnP9fPUL0+xjCFPBj0MXZ/M\nFndVzVFe/zy8Lo1XPtBAyIJ3j9FFk4zuROaOWCNtWGl7LOrpT1u4NWp2apBYE+ipytDu20W75zq1\ntZ7I7N+1lO+d3Ds1wgAEwhOd0jXHVhMjaS6YERt8flNJ1in4LknqbZgJBCNAAgGtgt6vf3gct102\nLOnnNdZzMbcRNJfFffjmsVGP88StJ/PILeNsn/eQ/tAk8sgaLfjh/bwfIh9nKvkQVICx+GCWdY2d\neBmu0dYfrWuI6oP3Wtleppd47hll0rtFfi7jjk1drPehI/UhC/5HZwr+eutJKTu3kaIYEUa/vGwY\nf7phTFiWbKrIKgV/55UjQlEz2UjQOgk1Svb5ImKok+1nLjIlcJgtptYxysoWFuRSGsO9ZMz+G9ij\njHema+GF+w66L6ub7OthDPO86cLBjE1Auc2PUUQMoiv46/8y3dJN8fSvxtOqMH5o43lje9qSDxon\nCmdFafHY0NAQugeS1UjayJzlO0M++NatWtAySZ2k4tE6xn2dm+MPlXVINU2+J+uoge35ZsUuTjim\ngy3rqSkTjC6MFTHx3TE9oq5zgznV34wXI+L7rx/Nzj01YTXZjf5yp4lURhItyRAP4wtvhGif0GTn\nPz+RYdExQYUeVKpHY1j45tOOHdyRvFx/3HaMEPulbMaqpEHLFjkcOqz9Rm31HqM+oGV+TiiTN5kE\nO5KlwjtzxvHdOKZnG2pqj1KQn8vj/1kMpLa+jBOavAX/47MGcMn4Pvzw9P7pFiXpBK3lWHNKVin+\niVAQxyLyIoKhQ1khx/ZpR0tDfZtcv59uun97aN920XaNS7IferOrwsn5rCznWcu08MlDh+u47s/T\nue7P00NNpddujZ5Vabbuy0o0pT28fzsuO7Uv+XnRH3UnoxyrTQMBrejWJeP7cOkELe7b5/OxdtsB\n7n15nu1j28EqMS40V5MCHXvi4I4c26ctowd1ZFi/dqFzRmtCn26avAXfskUuZ4/21mrNVILDz1hh\nduYOUIkS79n3srHFyAHteeY9rehobo6Pzu2K2LyzKqGWbD6fj8tO7Ru3aJRbA8Hn8/GjMwUd9XRx\nJ9aqOXkM4Nn3lvPse8vD2iy+/vkaRoj2EdsaMY9yTh/ZDdCMgjNHdadT2yIefXOR1a6OoousLNUA\nmvvN+Bz6fNoKq4zbRLCaaP5cbzSTaA6IFaWt8sPaVpqbi/z+6lF8vWQ7w/u7N0KSSZO34JsTPzpT\nMEKU8/3T+oUtv+CkXqHPdosj2SXe0NPLBBtjIkhujj9sriERzhzVnWFRRgF3fH84o4/pYKsJeDTG\nD+/CSRZtEc2ce2LPsL+H9WsXVa5gmQOAygPxa+2bq0uaR3LFMXzxTtxYltZ+IDJXIZXBLMEibEnJ\n5NUva36unzOO7xbxAu/WvhWXT+yXdFegWzJTKoUl7Vq35KYLh9CuNHzC5ryxvbhkQh9GD/IuDPC2\ny4fRt0spJw+NrfiSlWSSk9NYhzuZ/s0BPcq47rxBKckyPKZneIeg3Bw/P71wcKimeCxuefRLANpG\nKXEcr0lKLIXryEVjacEHIkKdCjzssmSXZCR73fy9IfTpUsL914/h8on94u+QYSgFnyWcfUIPrjt3\nkGfHG9SzDXdeOSIiasZMspJMcvz+xobGTax+0FXBEgsWmF1BuTl+HrghflWPg3qMfDRrfq9FbL2R\nWNnATn5BS0s1Ur/HDBtMhA5l0aNRkmEI9OlSyl1XjozbPzhTUQpekRDGrkZeZNgG4+wr9x/iWz1s\nMFNKDNhlYJRenn6/j4kjutK+dcuwyoxevMAe+be1fz1IrFMsWWu/QF0UD03E8Q8eil5KOREKYxgc\nowbEnqdojigFr0iIb9fsDn2eaCrR7IY7r9RKMiwzlFZtYgY85a1bcu13tUqT7Q3xz8FElz9eP5pH\nbomeEJZqDjlpXWj4MQ7o+QnaIC78R6pOkoK3Gm8M1nNBenWOrEza3HHlKBNC5AEvAT3Q2tpcLaVc\nZ9rmMuA2oAH4XEp5V2KiKjKdczyIwbeKQklH+ddEGTukE4N7taGkKJ/VW/ZT2io/lOTl9/kSemvd\nc+0ofvP8N472iTVK6O6kFLPBJffgqwu47yejgdS9hK08grv2ak1drEIomztuLfgfAPuklOOA+4D7\njSuFEIXAg2jNtscApwkhjklEUEXmcNEp1unvXijibHpES1u1wOfz0b9bazp4mKbuqvFGjAs7pLf9\njlLGKK3G9nnxvfhVNUeoqU3cqrc60y69121Tm6tJBW6nuicCL+ufPwNeMK6UUtYIIYYEuzwJISoB\n77otKNJKmwSaVcelmT6j9046gepDRzlYe5S//mdJaHmLvBwOG2LcR4pyVy/S2E2T7R/PqhRAXX0g\n7hF+/vhXALww+VTb53JCbo6vSY70ko1bC74jUAEgpWwAAkKIsBk2g3IfAvQEZqPICoyWUns9qmG8\njThwOzTXR7RzuyL6d2vN8H7l3HJRYzvAX30/vJhcTo7f3WR2kpRfsPHM2m32epcmTJTBQjJaVWYD\ncS14IcQkYJJp8Qmmvy3vHiFEP+BV4AdSypjjs7KyQnJzI5N0vOjKngqaipyQuKytWzc+zMGwuZ6d\nSz25BvUWYXhN5dp6JecZ5cV0bF/M5p3VjB7aFZgfWpebl2PrPOZtamMUqmxdVuha9rmGAmnRjmFc\n7uY8xn1yYmQ1Z8J9kgkyGImr4KWUzwHPGZcJIV5Cs+IX6ROuPinlEdM2XYF3gSullN/GO8/evTUR\ny+J1VM8Umoqc4I2sVYZY7O2VWqODgzWHPbkGuy3ug6Zwbb2+BzqVFtCptCDimIcP19k6j3mbvXuj\nN6TYt7eGihbuMqBrDNEy0eQyLnd6jczXtU53V504uCMzTb1s032fpEsPxHqpuHXRTAEu0T+fC0yz\n2OZ54EYpZXJbvStSjlVCSSJt9YzsORA7Yac5YqwQOm+lVnjsvLE9HVUOTZZ/Oh2tAlu2yGHSOSpm\nww5uJ1nfAE4XQnwFHAauAhBCTAa+ACqBk4A/CCGC+zwspfxfQtIqMgIrXeGV/shLQ2PiTOeiU/qE\nmoQM0ssdXHCSFslk1TzkpxcMjlgWu/WteyVtJ5M5kXLPEefz7EjNA1cKXkpZD1xtsfwBw5+pb1+i\nSAlW1nq8dnN2MTYXVzRy04WDefKdpZw7tlfM7R6+eax1ffckaUY7lSq8THoyJlXddvkwHno9rve3\nWdPkywUrUs+S9fZT252Sro48mc4I0Z7nfj0hbqy3uaVikLokuVLsWPDPf7DCyzOGzIv+XVt7eNzs\nRI2HFY45ZGGte1k2WGGNnUSeaL72eotuUMEqlmXFCeQ12HDNbdsdfYLX1Sn1cybSJ6C5oMwlhWOs\nVPnUBVu46JQ+KZdFEU60iopWceL3TjqB/QePxOyTGw878efxShk7YceeGltuIYWGegUqHGM1LA/2\n5PSCP1432rNjNTeiWfmFFq6b/LwcyhNsBm1nVNHSZQimmepDR6mrD3j6wsh2lIJXOKbE476vZspb\nJ7EUQpZh1/ruaGpZ+Lurjvfm/K3in79iX/yOVHaoqjkSfyNFGErBKxxz4pCOST1+prY/y0Tu+MFw\nV/v1cFJBMhYpNKZVrRnnqCdJ4RgrH+gEj2rRKJxhp91fMlm3PUU1aBSuUApe4ZgeHSKtv6KWar4+\nHeTl5jDu2E5pO/+qzUlodB0FZb87Ryl4hWP8fh8nDg5305x6XOLdnIwM6a2qS9vlmu8M5LyxPRk3\nJH2K3i6OukeZsPIG/Vp3UXUtb+X6uNmMUvAKV5wyrHPoc+tW+dbZkwnwozNF/I0UIS44qTfX6G0C\nM4ER/cstl09dsMX1Ma2it0T3Ml6YfCp/uHaU6+NmM0rBK1xhLDL1y0uHxdhS0Rz5TpRCaJt3VQNa\nYtyn8zZTud9+hE06Cps1dZSCVySMyihsGpw8VHPhPPTzkz0/dlubXb7mrazQ/9/Fa5+t5sFX7Reb\nrfWwaFlzQc2MKRKmyE2P0DgEE3MKVW0az/jRWQM4f1xv+ncv87xuud0IxoZAgO2VB9lXrcW073Zg\nwf93xno3ojVrXD09epOPl4AeQD1wtZRyXZRtXwMOSymvcimjIsOxypJMlJYtcnni9gkEjnpTpVKh\nZZ2WFXs7V+KGu56dw6UT+jre77Cy4B3jdmz9A2CflHIccB9wv9VGQojTAVWgJMtJVjf7Hh1LKCrw\nfnSg8B6n94CbGvSjj+kAwNXfGeB43+aKWwU/EXhH//wZMNa8gRCiBfD/gHtdnkOhUGQgndpatHpw\n+I6vd9EkO1iDpmW+ctvZxa2C7whUAEgpG4CAEMJclOL/gKcAleqWheSozkvNllsvGcpZo7qHLXM6\nhjtS59zdEoyiyYlSMVMRSdxXoRBiEjDJtPgE099hV1wI0Q8YKaX8vRBivB1BysoKyc2NrDqXaV3K\no9FU5ARvZG3bthVnj6lk7NDOSf3uze26popEZC0vL2Zg33I+/mZTaFlOjj/smPvidPgqKGi0B+PJ\nElzfUi9yV1ZWmLHXOtPkiqvgpZTPAc8ZlwkhXkKz4hfpE64+KaWx1Nt3ge5CiNlACVAuhLhDSvmn\naOfZu7cmYlm6upQ7panICd7KeskpWl/QZH335npdk41Xsub4fSG3SUNDIOyYVs+zkUOHGtXFY6/N\np3PbIsZb1DMyyrr/gBZxU11Vm5HXOl33QKyXiltn1hTgEuAT4FxgmnGllPJR4FEA3YK/KpZyVygU\nTQ+/QcFfeHLvsHVH6yI7SBkx+uA/m6dlt1opeCNBF020piaKSNw6Ut8AcoQQXwE3ofnbEUJMFkKM\n8Uo4hUKRuRgVbY8O4bVg4jVPd9PisV754B3jyoKXUtYDV1ssf8Bi2XRgupvzKBSKzCXHEBpprtVe\nVJDHoz8bx62Pf2W5r5vI2uBLQfULsI+6UgqFwhU5OUYFb7HeoaW9c08NW/RaNVbsP6j57ZWLxj5K\nwSsUClfUGCJlrBKdfDGCJ+eu3BWx7P+emc1vX/gm6j5fLd4OwPbKg07EbNYoBa9QKFwRr/l1LDfM\nngOHXZ/Xhfu+2aIUvEKhSBirfqlel7AY1LMMANG9tafHzWaUglcoFAlj5Rb3ukRRSVEL/VzKB28X\npeAVCkXCWFnwVsucMuPbrfzssRnsOVDLrGU79OMmfNhmg6rao1AoEsZK6bpVxMHWfM++v5zZy3YC\nMHPpjtB6FUVjH6XgFQpFwnjpg6+rbwB8IeUO4Q23vRgZNBeUglcoFAnjpQ/+xQ9X0q19eGbsO182\n9hNSBrx9lIJXKBQJ46UPfvbyncxevjPq+paqjaNt1CSrQqFIOv26lvL0r07hpxcMTvhYuaoXgW3U\nq1ChUCRMPH/77d8fTm6On5ED2qdIIgUoC16hUHhAPG+M0epWceypw5UFrzf5eAnoAdQDV0sp15m2\nGQo8r//5XynlPQnIqVAoMgxjww8nOjuQQK2BP1432vW+zRG3FvwPgH1SynHAfcD9Fts8A1wHjAKO\nEUJYdOpVKBRNFWM8erQJ1TuvGMHN3xsStsxKvY87thM9OsZvd9exjVIjTnDrg58IvKx//gx4wbhS\nCNEBaCWlXKAv+r7L8ygUigwlXMFbb9O3a6mtY405pgNfL9nuhVgKA24t+I5ABYCUsgEICCHyDet7\nAnuEEC8JIb4WQtyamJgKhSLTiNXww/GxcvyqSmQSiGvBCyEmAZNMi08w/W3+dX1AL+AC4BAwSwjx\nqZRyWbTzlJUVkpubE7E807qUR6OpyAlK1mTR3GTNzfWDXvW3fXlxQkq+TVnsFn9BMv0aZ5p8cRW8\nlPI54DnjMiHES2hW/CJ9wtUnpTxi2GQnsExKWalv/xUwCIiq4K26sDeVTvVNRU5QsiaL5iirUZ3v\n3h29E5MdWuX7mHBcF6Yt2Bpzu0y+xum6B2K9VNy6aKYAl+ifzwWmGVdKKdcDxUKINkIIPzAMkC7P\npVAoMpBEi37ddOEQOrUtpFPbQgryc+ndqcQjyRRB3E6yvgGcrlvmh4GrAIQQk4EvpJSzgF8AH6FN\nmn8spVyUuLgKhSJb2L3/EPdMOiE0EojWIeq5OybQEAg47vGqcKngpZT1wNUWyx8wfJ5DpK9eoVAo\nAFi6rpIzR3UP/d2l3NoP7/f78Mfo76qIjspkVSgUrhjQvSyh/VsV5of93aNDpC/5vJN6J3SO5o5S\n8AqFwhWXTeyb0P4jRXnY38YgnJYtcjl/XC+uOmdQQudo7igFr1AoXNEyP7FahYWmsr8+gxumfVlL\nzh/Xi7xcpaISQVWTVCgUrvC71L33TjqBb9fsZkCPcBeP0YIva9UiAckUQZSCVygUKaVzuyI6t4uc\nUDUmSl15pkilSFmLGv8oFApX+JIY2VJWrCx4L1AKXqFQKLIU5aJRKBTuSIIBX1qUTz+bFSgV8VEK\nXqFQZAwP3zw24cqUikaUi0ahULgiGWpYKXdvUQpeoVAoshSl4BUKhSuUtZ35KAWvUCgUWYqrSVa9\nycdLQA+gHrhaSrnOtNF50o4AAAxjSURBVM19wHi0l8g7Uso/JSSpQqFQKBzh1oL/AbBPSjkOuA+4\n37hSCDEYmCClHAuMBa4WQnRMSFKFQqFQOMKtgp8IvKN//gxNiRvZDxQIIVoABUADENmTT6FQKBRJ\nw62C7whUAEgpG4CAECJU3FlKuRl4E9io//u7lPJAgrIqFIoM48TBHbl8Yr90i6GIgi8QsG6TFUQI\nMQmYZFp8AjAs2IZPCLEF6B1svC2E6A28juaDzwNmorlsdkU7T11dfSA3N8fl11AoFIpmS9RwpriT\nrFLK54DnjMuEEC+hWfGL9AlXX1C56xwPzJFS1ujbLwYGA1OjnWfv3kgPTlPpVN9U5AQla7JQsiYH\nJau980bDrYtmCnCJ/vlcYJpp/RpgpBDCr78AhgDrUCgUCkXKcFuL5g3gdCHEV8Bh4CoAIcRk4Asp\n5SwhxBTgK33756SUGxKUVaFQKBQOcKXgpZT1wNUWyx8wfP4d8Dv3oikUCoUiEVQmq0KhUGQpSsEr\nFApFlqIUvEKhUGQpcePgFQqFQtE0URa8QqFQZClKwSsUCkWWohS8QqFQZClKwSsUCkWWohS8QqFQ\nZClKwSsUCkWWkhEKXgihuvd6jBCis/5/RvzG2YAQoiTdMjhBPVfe09Seq7TFwQshWgOXAf8F9kop\nDwshfFLKjAvM12W9FZgDzJVS7s5gWUuBO4AfA6OllFvSLJIl+jUtl1KuTrcs8dBl/TVwCHjAVBo7\noxBClAHXAu8BW6WU1Rl8r6rnKsmk5S0khDgX+AA4FrgB7eEhQ3/Y84B3gULgNOBByFhZrwP+p//5\nLFCfiVacECIHreT0/wkheqRbnlgIIW4EPkVrQ/lghiv3iWgGUwe0ct4PQcbeq+einqukk1IFr9eG\nB+gFvCylvAl4GK308Gn6Nhkx9DHI0QP4h5TyDrQbcqVhm4z5kfVG552BK6SUdwGjgBaZ9MAYrmlv\ntDLTdcBxxnaPmYQQoh0wBq0E9gP6KLO1YX2m3KvBVmhd0GS9XUp5L3CSEOJ7+jYZIauBbjSN52oE\n2gszY5+rWKTERSOEOAa4Ea2P6xPAncAq4AUpZZ0Q4h9ADynl+KQLEwddUV4BrABeRrOEvgJqgenA\nfGCelPLJdMkYRJf1SmAp8IreHze47v8B86WUH6VLPoMsg4FrgE3A00ALoA1wEjAS+JuUckX6JGzE\nIOsG4BlgAjAOWIbWbD4X2Cel/Hm6ZAxiuq5/1z8XAa9LKTcLIZ4GjpNSHp9GMQEQQvQBxkgp/6X/\nfQVao6Aa4Asy67nqA5wopfynxbqMea7skLS3evAtLIQ4A3gSmAX0Af4f8CpwNvBrIcR9wLfAAf1H\nTzkGWY8F/ob2MJ8A3AXMllJuA1oDjwL3A5cIIX6n75PqUZBZ1qVoVsXdQoie+rpcoAyoToeMJjn7\no/3+i9A6ez0I9JRSrkNrzF4AjNN9x0ZrNN2yDgPuAY4CB9B8xV8CvwCGZcDvb5R1KNpzVQ10An4j\nhLgf2AZUCSF+btw31bLqPAj8TAhxuv73K1LKrWTYc2WQ9RYhxKlBWYQQPv3eTOtz5ZRkChh0x/QE\nlkkpXwVuB1pLKRcAv0G7UPXAU8C/0X7sdBBsfDIY2K6/uSfrf39PCFEspVwnpXxBSrkKbd7gYiFE\ngdFqTqOsdwH9gPOFEK2llHVoFuitAGmQESDodjkGqJBS/gNNOVYBE4UQHfWeve8Do4F2+vbpGJ5b\nyXorcAQYiGaA/BbNMt4DXA9clqbf30rWn6Ndt2LgNfQ5Aynl3cDvgR5CCH8a3Ap5EHoZHUUbEV9h\nnEiVUq7PkOfKStYf67I2AH690dF60vtcOcJzBS+EmCCEeAX4nRBiJNoM+UlCiF8CC4ByIcSTwEop\n5WPA3VLKWjTrbp/X8tiU9W4hxDBd1kIhxDAp5QFgK5oiHS40Buq79gY+0+XOJFmH6P8AXgJ8QojR\n+v4pUZy6nP8BHhJCnKjLmSOEGKDL+SnaXMEIACnlf4HNwGQhxEzg8lTIaVPWz9B8xblSyo+B4Oii\nJ/BRGn7/eLL2AUqklG+iGU2gjUTXplIZmWQdpyvv36JNVB5G7wYnhMgVQgzMgOcqrqw0Gh7/IA3P\nlVs8VfBCiN7AvcA/0ZT5xWgPw3fQImZuk1KeC7QFrhNCDAfuFULMQBsSL/NSHgeyzgMuBU5Gi0J4\nRgjxIlCOFhrXUf8ek4UQH6FZoinzwTmQtQbtIQfN0lsCnAipiU4QQnQC7gOeA74GfoT2gHwEnKfL\nMR3N5dFH36cAzYIfCtwX9NFmiKzTgINAL/2lercQ4mO03//jVMjpUNb9wEAhREfgNiHE12iu0AVp\nknUmmhV8tR4Ouxv4EPiOEKKzPtIcossavK6pfK6cygqaol9GCp+rRHDbdDuE7pc6EZiLFh0zT0r5\nsdAiIzoCP0GbnCxFu2igWRd/QJvEWg1M1C25pBJH1k7ARWjKcwbQX0r5nhBiAvBrKeVZ+ovoZN2a\ny0RZT0Ebkr8kpawUQjwnpdycAjl/g/bbBoB3pJQf6Yp7C/AA8CLQRggxXlfwM4G7gcf17/hvKeUL\nyZQzAVlnAXdJKZ8QQmxC+/3fzVBZZwK/lVI+AfxWCDFDSvlpmmXdBtwnhHhPj3NfiPZCvwR4DG0O\n4UNgrJTykwyV9WLgcSFEHynlWiHE080pDv4x4I9oQ24JnCaEGK7HC/vQhjm/RBvyXKLvU4YWRVMk\npaxOhXK3KesR4GdSSolmuQH0B6brvriaVCj3BGQdAIQekhQo985ocyelaCOdv6L5WAv1YfZstInJ\nsWgvqt8JLVS2LfC1/rBNS5FydytrG2COEKKllHJPipR7Itd1lhCiJUCKlHs8Wb/SZbxdl2mjvv11\nQohlwEhdB6RCubuV9Xpd1lH68iah3CFBBS+EKAIEmpUzXv/iT6HNQH+B5nZ5Fi1SYjEQEEK8jzZJ\n9ZiUcn8i50+CrM+guToAzhJC/BdtOPy/VA7FEpU1VXLq528jpfyFlPJ54BVd7of09Q3AW2j32Wdo\n3+d5tN//JSllfQqva6KyHkqRnNkmK8C/gE5CiM66W+QRNAPvR1LKV5qIrD+WUr6WQlk9IeE4eN2P\nnoMWjz1NSvmubk0MlFIuFloM/C+klD/R3Qs9ZJrS023IOhC4Q0p5tb68XGohkkpWaxk7AoPQ4pn9\naKF6X6BFIFwgpZwvhOiHlql8vb5bsZQypZPpStaMkPVXwE1o1vPJUsp3lKzJx7NEJyHELWgZX/+R\nUi4UWtr0fuAMtAnKW6SUhz05WYLYkPVmmSEp6U1FVqHF3k8Fzgd+CJwD/Blt8vQktEzAg9GPkDqU\nrMnBhqxXSimr0ydhI01J1kTwYpI1GNM6BZiE5gdeqP9/LJq/PSOUuwNZM0FhNhlZdYYASCn3Ak8I\nIQ4Cp6JNtN+YKUpIR8maHOLJmkkKsynJ6hpPSxUIIc5Gcyn0Qgsje1CmMKbVCUpWbxFCnIMWM/4u\nWtr8N8AfUzl3YRcla3JQsmYeCVvwJi5Esy7/JKV82eNje42S1VvaoqWcX4g20fdqmuWJhZI1OShZ\nMwzPFLwegvQNGeKOiYWSNSlsQJv0+1sGuY2isQElazLYgJI1o0hbww9FdiEytFGDFUrW5KBkzTyU\nglcoFIosJePLXSoUCoXCHUrBKxQKRZaiFLxCoVBkKUrBKxQKRZbidRy8QtFkEFqLQ4lWrAu0rj4z\ngD9IrdtUtP2ukCmqW69QJIKy4BXNnQop5XipNXyfiNb2LmrSi17i+Lcpkk2hSAhlwSsUOlLKWiHE\nrcBqIcQgtKY0bdCU/ptSygeBF9B6nE6RUp4hhLgUuAWtRn8FMElKWZmmr6BQhKEseIXCgJTyKFpb\nxHOAd6WUE9Aaa9wphCgBfodm9Z8hhOiG1vD8NCnlOLTOZXemR3KFIhJlwSsUkZQCO9Caxd+I1j2r\nAM2aNzIGrX3iJ0IIgBbA+hTKqVDERCl4hcKAEKIQrWPWdDSFPVZKGRBC7LbY/DDwjZTynBSKqFDY\nRrloFAodvTPW48CnaE1WluvK/TygEE3hN6BF24DWv3OU3ikIIcQlQojzUy+5QmGNqkWjaLaYwiRz\n0JqoTEHzow8AXgO2A/8FBgPDgdHAfKAOOBk4D7gNqNH//VhKuTOV30OhiIZS8AqFQpGlKBeNQqFQ\nZClKwSsUCkWWohS8QqFQZClKwSsUCkWWohS8QqFQZClKwSsUCkWWohS8QqFQZClKwSsUCkWW8v8B\neknPRvQEKCYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fda2f273110>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "ZGc4wk6A32qA",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Hurst exponent in practice\n",
        "I'll quote from the [[BSMRT]](https://www.quantstart.com/articles/Basics-of-Statistical-Mean-Reversion-Testing) blog I mentioned:\n",
        "\n",
        "\n",
        "> A time series can then be characterised in the following manner:\n",
        "* $H<0.5H<0.5$ - The time series is mean reverting\n",
        "* $H=0.5H=0.5$ - The time series is a Geometric Brownian Motion\n",
        "* $H>0.5H>0.5$ - The time series is trending\n",
        "\n",
        ">In addition to characterisation of the time series the Hurst Exponent also describes the extent to which a series behaves in the manner categorised. For instance, a value of HH near 0 is a highly mean reverting series, while for HH near 1 the series is strongly trending.\n",
        "\n",
        "\n",
        "Lets compute the Hurst parameters with python:\n"
      ]
    },
    {
      "metadata": {
        "id": "23abicIZ2fsV",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "dadb2550-83a3-41ac-b862-31212c60740e"
      },
      "cell_type": "code",
      "source": [
        "from numpy import cumsum, log, polyfit, sqrt, std, subtract\n",
        "from numpy.random import randn\n",
        "\n",
        "def hurst(ts):\n",
        "\t\"\"\"Returns the Hurst Exponent of the time series vector ts\"\"\"\n",
        "\t# Create the range of lag values\n",
        "\tlags = range(2, 100)\n",
        "\n",
        "\t# Calculate the array of the variances of the lagged differences\n",
        "\ttau = [sqrt(std(subtract(ts[lag:], ts[:-lag]))) for lag in lags]\n",
        "\n",
        "\t# Use a linear fit to estimate the Hurst Exponent\n",
        "\tpoly = polyfit(log(lags), log(tau), 1)\n",
        "\n",
        "\t# Return the Hurst exponent from the polyfit output\n",
        "\treturn poly[0]*2.0\n",
        "\n",
        "print (hurst(butterfly_ts))"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[0.43670626]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "WWnQCtWxIDOZ",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "when I computed the result it was 0.43688, less than 0.5, and indicating the signal is mean reverting. However, it does not seem to be _very_ mean reverting as we would have preferred for it to be closer to 0. \n",
        "\n",
        "### Augmented Dickey-Fuller in practice\n"
      ]
    },
    {
      "metadata": {
        "id": "9sw_8rXO4uBV",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 88
        },
        "outputId": "fbf23f51-f880-434a-bce8-c5a890c5fdee"
      },
      "cell_type": "code",
      "source": [
        "# Import the Time Series library\n",
        "import statsmodels.tsa.stattools as ts\n",
        "\n",
        "# Import Datetime\n",
        "from datetime import datetime\n",
        "\n",
        "# Output the results of the Augmented Dickey-Fuller test for Google\n",
        "# with a lag order value of 1\n",
        "print (ts.adfuller(butterfly_ts[0], 1))\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(-4.2472975794099765, 0.000547416540345047, 1, 7096, {'5%': -2.8619473981777572, '1%': -3.4312718809379334, '10%': -2.5669868539820775}, -28994.736676103945)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python2.7/dist-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
            "  from pandas.core import datetools\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "metadata": {
        "id": "g8iT57ro5Pya",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "The first value is the calculated test-statistic, while the second value is the p-value. The fourth is the number of data points in the sample. The fifth value, the dictionary, contains the critical values of the test-statistic at the 1, 5 and 10 percent values respectively. At the time when I computed the test, the calculated test-statistic was -4.2298 which was less that the 1% critical value of -3.4313 (and less that the critical values at 5% and 10%), indicating the the butterfly 2-5-10 is mean reverting (although students of Statistics might argue that the test rejected the hypothesis that the Time Series is a random walk). The p-value corresponding to the test was 0.058% - so well below the 5% probability threshold usually considered as signal that the time series could be non mean reverting. Read [What is the meaning of p values and t values](https://stats.stackexchange.com/questions/31/what-is-the-meaning-of-p-values-and-t-values-in-statistical-tests) if you need refreshing your memory.\n",
        "\n",
        "While Statisticians might debate the merits of the tests, they provide some basis for us to go beyond the accumulated market memory of mean reverting signals.\n"
      ]
    },
    {
      "metadata": {
        "id": "OmDHx1r4tHjf",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## How can I find Mean Reverting tradeable signals ?\n",
        "\n",
        "We could find some asset classes that are mean reverting by nature (some FX crosses like AUD/NZD for some periods of time), but the opportunity set increases when we can enter in linear combinations of asset classes, or portfolios of long/short assets that we can enter and exit over a period of time (without rebalancing - to minimize transaction costs).\n",
        "\n",
        "### Pairs trading\n",
        "The internet is already littered with examples of pairs trading (taking a couple of assets and finding if buying one versus selling the other weighted by a certain hedge ratio results in a mean reverting strategy), so instead of writing another example please read:\n",
        "* [Identifying Asset Pairs for Pairs Trading](http://hkopp.github.io/2017/07/identifying-asset-pairs-for-pairs-trading) - were the basic recipe of magically finding two stocks, do a linear regression on them to find the hedge ratio and test for mean reversion (as above) is explained, counterbalanced by \n",
        "* [Why I won't teach pair trading to my students](https://www.marketwatch.com/story/why-i-wont-teach-pair-trading-to-my-students-2012-10-01) (with real cases of blow-ups and a reminder never to forget the fundamental drivers which should drive the mean reversion),\n",
        "\n",
        "The basic recipe of using linear regressions to identify the hedge ratio suffers from the fact that the ratio will differ depending on which asset class is used as the explanatory variable (see [Effect of switching response and explanatory variable in simple linear regression](https://stats.stackexchange.com/questions/20553/effect-of-switching-response-and-explanatory-variable-in-simple-linear-regressio) - it turns our the formulas give different results unless all the points are exactly on the line~) annoys all the practictioners trying to find the correct ratio."
      ]
    },
    {
      "metadata": {
        "id": "N6MVN4VsY-4j",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 488
        },
        "outputId": "aeab9fb0-d7e8-4a4b-c185-eb1b9fb4e22f"
      },
      "cell_type": "code",
      "source": [
        "\n",
        "wings_body = butterfly_ts.copy()\n",
        "\n",
        "\n",
        "# A common practicioner trick is to express the mean reverting signal as\n",
        "# a scatter of two spreads.\n",
        "\n",
        "# We start with the 'wings', which is just the average of the outlying \n",
        "# maturities (2y and 10y)\n",
        "wings_body['wings'] = (TimeSeries.dot([0.5, 0, 0.5]))\n",
        "# the fine line: the cointegration test requires the underlying time series\n",
        "# to be non-meanreverting, so we test them using adf.\n",
        "print ('wings adfuller p-value:',ts.adfuller(wings_body['wings'], 1)[1])\n",
        "# versus the body (5y)\n",
        "wings_body['body'] = (TimeSeries.dot([0, 1, 0]))\n",
        "print ('body adfuller p-value:',ts.adfuller(wings_body['body'], 1)[1])\n",
        "print ('if both the above p-values are > than 5% the series are not mean reverting')\n",
        "print ('and can be used in the cointegration test')\n",
        "# To plot the scatter with a line:\n",
        "import seaborn as sns\n",
        "sns.lmplot(x='wings',y='body',data=wings_body,fit_reg=True) \n",
        "#wings_body.plot(x = 'wings', y ='body', kind = 'scatter')\n",
        "\n",
        "import statsmodels.api as sm # import statsmodels \n",
        "\n",
        "X = wings_body[\"wings\"] ## X usually means our input variables (or independent variables)\n",
        "y = wings_body[\"body\"] ## Y usually means our output/dependent variable\n",
        "X = sm.add_constant(X) ## Intercept the model\n",
        "\n",
        "# body = beta_w * wings + constant_w\n",
        "model_yX = sm.OLS(y, X).fit() ## sm.OLS(output, input)\n",
        "\n",
        "# Print out the statistics\n",
        "print (\"Wings as explanatory variable have a beta_w of \", model_yX.params[1])\n",
        "\n",
        "# wings = beta_b * body + constant_b\n",
        "# Note the difference in argument order\n",
        "model_Xy = sm.OLS(X, y).fit() ## sm.OLS(output, input)\n",
        "\n",
        "# Print out the statistics\n",
        "print (\"Body as explanatory variable have a beta_b of \", model_Xy.params[1][0])\n",
        "print (\"and beta_w is not the inverse of the beta_b:  \", 1 / model_Xy.params[1][0])"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "('wings adfuller p-value:', 0.9091249800366502)\n",
            "('body adfuller p-value:', 0.8583858294638854)\n",
            "if both the above p-values are > than 5% the series are not mean reverting\n",
            "and can be used in the cointegration test\n",
            "('Wings as explanatory variable have a beta_w of ', 1.0217157779654704)\n",
            "('Body as explanatory variable have a beta_b of ', 0.9864733462355418)\n",
            "('and beta_w is not the inverse of the beta_b:  ', 1.0137121330404688)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt81NWd//HXXELuQJQBVITIuhxE\nBbWtyEUpai3beCu20F1dd1tt1d9q96ft1t3fr7uldn/ttlVpa7sFLdVad4usomBSqWJbFBPFIoJD\n5CgyAVEuYeQSSAi5zO+PSWjAEGbI9zKX9/OfJl9mzvfz+JbHm+P5nksgkUggIiLeC/pdgIhIvlIA\ni4j4RAEsIuITBbCIiE8UwCIiPgn7XUBfGhubcnqKRkVFCbt3N/tdhu/0HJL0HJJy8TlEIuWB3q6r\nB+yjcDjkdwkZQc8hSc8hKZ+egwJYRMQnCmAREZ+4OgZsjDkHWALMtdb+1BhzOvBrIARsA/7WWtvq\nZg0iIpnKtR6wMaYUeAB4ocfle4CfWWsvBjYCX3Lr/iIimc7NIYhW4DPABz2ufRJY2vXzM8DlLt5f\nRCSjuTYEYa1tB9qNMT0vl/YYctgJnNJXGxUVJTn/RjQSKfe7hIyg55Ck55CUL8/Bz3nAvc6L6ynX\n5gIeLRIpp7Gxye8yfKfnkKTnkJSLz+FY/6B4PQtivzGmuOvn0zhyeEJEJK94HcDLgeu6fr4OWObx\n/UVEMoZrQxDGmI8B9wGVQJsx5nPA9cAjxphbgM3Ar9y6v4hIpnPzJdxqkrMejvYpt+4pIpJNtBJO\nRMQlm7c38cizG4755xm9G5qISLZa926c/3z6TQ61dfJPx/iMAlhExGF167fzy5q36Ojse0ddBbCI\niIOee+09Fr7wzuHfr5s2+pifVQCLiDggkUjw5IpN/PaVzQAEAvB3M8ZyyYRTj/kdBbCISD91dHby\n6DLLS+u2ARAOBbntmrM5f0ykz+8pgEVE+uFQWwfzlqznjY27ACguDPPV687FjKw47ncVwCIiJ6j5\nYBs/eWIdb2/dC8CgsgHcNes8Th9altL3FcAiIidgz/5W7n/8DbY2HgBgaEUxX5t9HpHBxcf55p8p\ngEVE0rT9w2buf/wNdu09CMCoYeXcOWsCA0sHpNWOAlhEJA0N2/cxd9FamprbADhrVAW3zzyX4sL0\n41QBLCKSovqGD3lg8Zu0HuoA4OMmwpevOpuC8Int6qAAFhFJwaq3dvDQM/WHV7dNv+A0rr98DMHg\ncc+WOCYFsIjIcbyweiv//fzbdC8svmbqGVw9pZJA4MTDFxTAIiLHlEgkWLIyxtKXG4DkOWo3XDGG\n6ReMcKR9BbCISC86OxM89pzlj28kT04LhwJ85aqz+fjYoY7dQwEsInKUtvYOHnymntW2EYDCASG+\nOvNczqo8ydH7KIBFRHpoaW3ngSfXsWHLHgDKSwq4c9YEKocPdPxeOhFDRKTL3gOH+P5/v344fAGC\nAYjGPnTlfuoBi4gAO/e0cP/CN9i5p+WI6/ua21i6MgbAlZMqHb2nesAikve27Gjie79e/ZHwBUgk\noL0jwfOvvef4fdUDFpG8Zrfs5idPrqOltaPPzx1oaXP83gpgEclbq20j85eup72j87if7e+ii94o\ngEUkL7249gN+tWwDib7PzTzslJNLHK9BASwieSWRSFBdt5mnXtyU8ndKi8LMuvRMx2tRAItI3uhM\nJFi4/B2Wr956+FogQJ+94CGDCrlxxljOOeNkx+tRAItIXmjv6GRBzVu8Wr/jiOt9h28RP7htsms1\nKYBFJOcdPNTOz56Ksj7FBRXBABSEQ9w4w7halwJYRHJaU/MhvvfY62z/sDmlzwcDUFQYZsbEka4M\nO/SkABaRnLVrbwv3PPIn9qc4hzcQgL8cMYiqyZWuhy8ogEUkR73fuJ/v/OpPHGo//hzfbmNOH8zd\nf3OBi1UdSQEsIjnnude2sPCFjWl9Z3DZAKomjXKpot4pgEUkpzxUXU9ddHta3xkRKWXWpWd6MuzQ\nkwJYRHLG479/J63wDYcCfPVz4z0P3m7aDU1EcsKyV7fwu1Xp7Vh29dQzfAtfUA9YRLJcZyLBz5+O\nHj4+KFVTzh3u+P6+6VIAi0jWau/o5N8f/RNbduxP63vjKiu4qWqcS1WlTgEsIlmpta2Dbz/8WsoL\nLLrNnDba955vNwWwiGSd/S1t/Mdjq9MO30wYduhJASwiWWV3UyvffngV+5rTO6EiABkx7NCTAlhE\nssa2+AG++YtXU95EvafPThvtfEH9pAAWkayw6YN9/Pujf0r7e8WFIf7qolEZNfTQTQEsIhkvGotz\n/+Nr0/5eJr1w640CWEQy2iv123lwaX1a3ykIBbjDxxVuqVIAi0jGev5P7/Gb5e+k9Z1BpQXMveNi\nlypylgJYRDJOIpHgwWfW82r9zrS+FxlcxPdvde8IIacpgEUko3R0dDL3f9YS3ZTa8UHdsi18QQEs\nIhmkrb2DO+77A++lubR45LAy5nzxQpeqco8CWEQyQvPBNuY8/Bq79h5M63uZPtOhLwpgEfHdnv2t\n3P/42rTD967ZEzJ+pkNfPA1gY0wZ8ChQARQC37bW/s7LGkTEO9V1DTxbt5mWQx2Hrx29MGLH7mb+\n47HX2XvgUMrtBgIw+ZzhWR2+4H0P+O8Ba639F2PMqcDvgbEe1yAiHrjzgZfYe+Cj+zW0tHaw5KUY\nAOeecTLffWw1bWkcnAnw2Uuyd9ihJ68DeBcwvuvniq7fRSTH3D53Bc2tHcf8847OBM++spmnXtyU\n1r4OgUDuhC9AIHEiu1r0gzFmGXAmyQCusta+cqzPtrd3JMLhkGe1iUj/vG538u1fvEJnp7O5EgjA\nwNIBXH3JXzDrsjGOtu2RQG8XvR4DvgHYYq2dYYyZACwAPn6sz+/end5en9kmEimnsbHJ7zJ8p+eQ\nlO3PIRqL85Mn1jkevsNPKuG7X7no8O/Z+IwikfJer3t9KOcU4HcA1tq1wKnGGHVxRXJATW0D7R3O\nhm9pcZiRw8ocbTOTeD0GvBGYCDxpjBkF7LfWHnugSEQyXnVdA8+veo+mlvQ2SE/FkEHFTB1/iuPt\nZgqvA3g+8EtjzIque9/q8f1FxEHVdQ1Uv9zAoTRnMaTqummjs36qWV88DWBr7X5glpf3FBH3PP/a\ne66F77jKipwOX9BKOBE5AdFYnJraBprSPJctVeMqK/j6F853pe1MogAWkbRU1zVQU9tAa5uzPV8z\ncjBVk0Yx/cLKrJzpcCIUwCKSsuq6Bp5+cRNOzTQLAKFggKsvPiNnFlekQwEsIimJxuIse3WLY+Eb\nDAQYXD6AT55/Wl6GLyiARSRFK9dto73DmWGHUDDAQ9+Y7khb2czrhRgikoWisTjr3o1zyIFx30AA\nPmYiDlSV/dQDFpE+Vdc1sOSlGB0OjT0MHZzbiyvSoQAWkWOKxuKOhu+ISCmzLj0z5+f3pkoBLCK9\nisbiPLS03pHwDYUCXDM1P2c69EUBLCIfEY3FeXLFJvb3c3+HogEhRg0vp2rSKPV6e6EAFpEjOLG/\nw5BBRdw4wyh0j0MBLCKHzXl4FVvSPBK+pyGDCrlxxlgFb4oUwCICwL0L1/QrfAF+cNsUh6rJD5oH\nLCJEY3HqG3b3q41xlRUOVZM/FMAiws8Wv9mv70cGF+XF7mVO0xCESJ67e15tv3Y2m3LucG6qGudg\nRflDASySxxbU1NO45+AJfVeLKvpPASySB6KxOIt+/w7b4i0kEglKiwqoGFh4wi/d1Ot1hgJYJMdF\nY3HmL1nPgYPth681tbSd8CGaIyKlCl+H6CWcSI57dNmGI8K3P0qKwsy69ExH2hIFsEhOu3fhGnbt\nbXWsvRkTR2rM10EaghDJQdFYnEeXbXAsfIMBuPaS0dpMx2EKYJEc0v2ybWtjs2NtaraDexTAIjkg\nGbwb2dp4wNF2w6EA99w00dE25c8UwCJZrr8b6PRl4rhhrrQrSQpgkSzmVviGggEuOnuYppu5TAEs\nkqWisbgr4atFFt5RAItkmWgszsp123htw07H2w4GUPh6SAEskkWisTi/rHmLPfsPudL+tZeMdqVd\n6Z0WYohkkUeXbXAtfKecO1zzfD2mHrBIllhQU+/oqrZuA0sLuPnKcZrn6wMFsEiGq65roLq2gUP9\n2LP3WPTCzV8KYJEMtqCmnpff3O5K2yWFIYWvzzQGLJKhorG4a+FbEArw0zunudK2pE49YJEMNe/p\nqCvtlhaFueWas11pW9KjABbJQNFYnObWDsfb1cY6mUUBLJJBuhdZrHs37njbM6dpO8lMowAWyRCv\n2508uWITOz5s5uAhZ3u/Ct/MpAAWyRDLV21h6879dHQmHG13XGWFwjdDKYBFMsSr67c5Hr6a55vZ\nFMAiGWDOw6scX2hx1+wJetmW4RTAIj5y+uy2bjOnjVb4ZgEFsIhP3NrZTJvqZA8FsIjHFtTU89pb\nOznU7vzeDuMqKzTmm0UUwCIecvP8Nr1wyz4KYBGPVNc1uBa+muebnbQZj4gHorE4S16KudK2xnyz\nl3rAIi6LxuL8aNFaHJ7iC2jYIdspgEVc4tYUs24adsh+CmARhy2oqacuut2VHm83LS/ODZ4HsDHm\neuAbQDvwb9baGq9rEHGLmydYdBtUWsDXv3C+q/cQb3j6Es4YczLwLWAqcCVwjZf3F3HbK+t3uNp+\nYUGQuXdc7Oo9xDte94AvB5Zba5uAJuArHt9fxDULauod30ynp8jgIr5/62TX2hfvBRIJFweqjmKM\nuRs4CzgJqADmWGtfONbn29s7EuFwyKvyRE7Yjxa+zguvveda+3/7mbOYddkY19oX1wV6u+h1DzgA\nnAx8FhgF/MEYM8pa2+u/Art3N3tZm+cikXIaG5v8LsN32f4c3B737d7VLJufUTqy/e9DbyKR8l6v\nex3AO4Baa2078K4xpgmIADs9rkPEEW6GbwBY8M+XutK2ZAavV8I9B1xqjAl2vZArA3Z5XIOII9zu\n+d45e4JrbUtm8LQHbK193xjzBPBK16U7rLXObwkl4qJoLE5NbQP2vb2u3UObqecHz+cBW2vnA/O9\nvq+IE6KxOE+u2MT2uHvvJ36pYYe8oZVwImlY9PuNbIs3uzbdrLykwJV2JTMpgEVSVF3XwPu7DuDW\nzM2CUJBPfeJ0dxqXjKQAFknRs69sdi18gwH4wqcN08ef4s4NJCNpP2CRFCyoqaeltcO19q+9ZLQW\nWuQhBbDIcVTXNbg63Uw7m+UvDUGI9MHtub7aUD2/qQcscgxe9HwVvvlNASzSCzfPcAMYOaxMe/qK\nhiBEjragpp7a6HZXZjwEgMkadpAuCmCRHtwadggGoKgwzIyJI/XCTQ5TAIv0sHSl88MO4VCAvzht\nEFWTRml/BzmCAliE5Jjvzxa/SXuHs+MO4WCAr35uvIJXeqUAlrwXjcW5//G1jrcbDMDVF5+h8JVj\nUgBLXnNrnm8oGOCis4dpvFf6lFIAG2MWAr+w1i53uR4Rz8x5eBVbdux3vN0LzxrK1PGnqOcrx5Vq\nD3gxcKsx5sfAQuBha+1W98oScZdb4VteUsCt15zjeLuSm1JaiGGtXWSt/RxwIfAO8BtjzG+NMZe7\nWp2IC+6eV+tK+AYDaDtJSUvKK+GMMSXAdcDNXd97BviqMebfXapNxHH3LlxD456DjrcbDgW49pLR\nGvOVtKQ6BvwwcAXwNPA1a+3aruvzgD8B33StQhEHVNc1UF3bwKE2548gHFw2gC9VnaUxX0lbqmPA\nLwO3W2sP9LxorU0YY252viwR51TXNbB4xSZX2i4uDCl85YT1GcBdPd/umemTjTFH/Lm19kvW2jUu\n1SbiCDdWtwGUFIW59ZqzFb5ywo43BrySZO+3EzgJWAtEgWGAe8fCijgkGos7vroNkvN8Z0wcqfCV\nfumzB2ytXQBgjJlpra3qvm6MmQs85XJtImmJxuKsXLeNxj0tRAYXM3X8KSz6/UbH71NYEKJq8ii9\ncJN+S3UMeKQxZrC1dk/X7+XAaJdqEklbdV0Dz6yM0dbV2928fT/R2Ic0H2x39D7m9EFUTa5Uz1cc\nkWoA/xzYaIyJkRwTPgP4f65VJZKGaCzO0y9uorPHSENnIuFo+JYUhTnnjJO0yEIclVIAW2v/0xjz\nGHAmyT2l3+3RGxbxzb0L11DfsNv1+5QXFzBVR8aLw1KdB1wG/CPwCZI94FeMMT+y1ra4WZxIX9xa\nTny0IYOKuP6KMRp2EMelOgTxELAVmE+yB3x517UbXKpL5Jiq6xp49pXNtLR2uH6vIYMK+cFtk12/\nj+SnVAN4mLX2r3v8Xm2M+aML9Yj0qbqugaUvxWjvdOHAth4CQKSimBuuGOPqfSS/pboXRGnXXhAA\nGGNKgSJ3ShI5tmfrNrsevuFQgDEjB3ODhh3EZan2gOcDbxljVnf9/jHgX90pSaR3d8+rpeWQu8MO\nQwYVcuOMsQpe8USqAfw4cBrw1yR7vg8CT7hVlEhP1XUNLHkpRofLPd+K8kJ+cNsUV+8h0lOqAfwE\nsBP4CcnhsSkkQ/kql+oSAZLh+9SKTbgZvcEgDBlUTOXwchfvIvJRqQbwQGvtX/X4/efGmBfdKEik\npyUvxVwN34JQkIqBhRQXhjXPVzyXagC/Y4w5xVq7DcAYM5zkyRgijquua2Dpypgrm+gcbcTQMiKD\ni3SGm/jieNtRvkRy4UUR8K4xZgPJndHGAq+7X57kG68WVwCYkYO5+28u8OReIr05Xg9YJ12IZ+5d\nuMaz8C0tClM1aZQn9xI5luNtR7nCq0Ikv1XXNXiyp0O3W7SRumSAlA/lFHHLohfedu3IoN7MnDZa\n4SsZIdWXcCKOi8biPLpsA7v2tnp2z3GVFdpIXTKGAlh84dXiip5mTtOx8ZJZFMDiOTdPKe5NMAD/\ne9YEDTtIxlEAi6eisbin4TuwpICbrxqn8JWMpAAWTz3wxDpP7hMIwJ3q9UqG0ywI8Uw0Fj98aKbb\nFL6SDRTA4oloLM5DS+s9udeUc4crfCUraAhCXBeNxZm/ZD0HHD4ivjfjKiu4qWqc6/cRcYICWFz3\n6LINrodvKBjgorOHKXwlqyiAxVXRWNz1hRYDCoLcPvNcDTtI1lEAi2sW1NTz8pvbXb1HKBhQ+ErW\n8iWAjTHFQBT4jrX2ET9qEHfdu3CNq5vrDCgIEg4FmTFxpMJXspZfPeBvAh/6dG9xUTQWZ97TUZpb\n3Ts8syAc5Lwzh2gTdcl6ngewMWYsMA6o8fre4q5oLM7PFr9Ja1una/coLQ5zy9XaSlJygx894PuA\n24G/O94HKypKCIdD7lfko0gk+w+CfN3u5H+Wv019Q5xOl7I3HAoytrKCz182hgvMUHdukgFy4e+D\nE/LlOXgawMaYG4E6a23MGHPcz+/e3ex+UT6KRMppbGzyu4x+icbiPPbc2+za2+Ja+I4cVsacL154\n+Pdsf2bHkgt/H5yQi8/hWP+geN0DrgJGG2OuBEYArcaYrdba5R7XIQ6pqd3Mrj0HcWtXySGDio4I\nX5Fc4mkAW2tnd/9sjJkDNCh8s1c0Fie2fR+dCefTt7AgxOCyAVx/xRjH2xbJFJoHLGmJxuKsXLeN\nzTua2L2vlUPtzo87FA0IMf4vTtYsB8l5vgWwtXaOX/eWE1Nd18CyV7dwqK2Djs4ELnR8GVdZwde/\ncL7zDYtkIPWAJSXRWJxlr26hvb3TlfAtLAjy86990tlGRTKcAlhSUlPbQMvBdtx41zb85BK+++WL\nXGhZJLMpgOW4FtTUY9/b60rbd82ewPQLK3Nu2pFIKrQhu/Spuq6Buqg7G+rcNVunVkh+UwBLn55/\n7T1X5vjq1AoRDUFIL6KxODW1DWx8fx8dLqRvRXmhNk4XQQEsR+leWty4u8WVF27DKooZNTw/1vmL\nHI+GIOQIK9dtY/e+VlfCtyAcpKgwzNTxp7jQukj2UQ9YjvD21j20dbizq87oUwdSNWmUxn5FuiiA\n5bAFNfXsaTrkeLuBAHz2ktFcOanS8bZFspkCOM917+3wzta97G5y/vDMUDDANRefofAV6YUCOI9F\nY3GeXLGJg63troRvcWGI2649R0MOIsegAM5T0VicR57dwJ6mQ65tJ6nwFembAjgPdfd8P9znfK8X\nkmO+VZP1sk3keBTAeaimtoH3driz90JBKMBVUzXmK5IKBXCeicbivPP+XseXFwcDUFpcwJevGqee\nr0iKFMB5pqZ2s+OHZ5YUhYkMLua6aaMVviJpUADnie79HdzYVvL0SClVkysVviJpUgDngWgszs8W\nv0lrm/Mr3KacO1wb64icIAVwHpj3dNSV8J05TavbRPpDm/HkuHsXrqG5tcPxdsdVVih8RfpJAZzD\nqusaqG/Y7Xi7OrlYxBkagsghyRdtm9nauJ/2jk5Hhx1CwQAjhpZppoOIgxTAOSIai/Nfz73N7qZW\nOjsTtDs00TcQABJwysklCl8RhymAs1z39LJ33t9Hp0vHB33y/NM03iviAgVwFquua6CmtsGVGQ4D\nwkGunFKp4BVxkV7CZaloLM6yV7e4Er7BYEDhK+IB9YCzTPcG6ms37nIlfMOhIBeMGaLwFfGAAjiL\ndL9o27O/1Z2ebwAGFAR1aKaIRxTAWSIai/PQ0nqaWtpcu0dBOMSMiSM100HEIwrgLNC9gbqb4Tsi\nUsqsS89U+Ip4SAGcBVau28bB1nZX2g4GA1yrQzNFfKEAznDRWJw17+yird35Md+SwjC3Xnu2er0i\nPlEAZ7BoLM4va95yJXy1k5mI/zQPOIPV1DawZ/8hx9sNBQMKX5EMoB5whprz8Cq27NjveLsBYFDZ\nAMfbFZH0KYAz0N3zamncc9DxdoOB5EKLT55/muNti0j6FMAZZkFNvePhW1wYorMzQWlxgTbWEckg\nCuAMEo3FefnN7Y62OSJSyj03TXS0TRFxhgI4g/xo0VrH2gqFAoyIJDdQF5HMpADOEHfPq8Wp7XwD\nwJmnDtRR8SIZTgGcAeY8vMrRcd87Z09Q8IpkAQWwz+5duMax6WbBIFx7sY4NEskWWojho3+dX+vo\nqcXXXqzVbSLZRAHsk+q6Bt54u9Gx9gaWFCh8RbKMAtgnT63Y5FhbBaEAN181zrH2RMQbGgP2WDQW\n5xfP1OPEhIfSojAjhpZRNWmUxn1FspAC2EMLauodWWgRCCSPi7/3f01xoCoR8YsC2APRWJx5T0dp\nbu3od1vBIISD2s9BJBd4HsDGmB8AF3fd+3vW2sVe1+ClaCzO/CXrHQnfglCA8tIB2s9BJEd4GsDG\nmOnAOdbaScaYk4E1QE4H8Mp12zhwsH/HCRWEg5w0sIjvfeUih6oSkUzg9SyIF4HPd/28Byg1xoQ8\nrsFTq97a2a/vF4QCnDqklFHDyhyqSEQyhac9YGttB3Cg69ebgN92XctJt/zwD/1uo2JgEQBTx5/S\n77ZEJLP48hLOGHMNyQC+oq/PVVSUEA5nZwf58/9STVtH/yabDSgIMnbUSVx24UguMEMdqiwzRSLl\nfpeQEfQckvLlOQQSCYe24EqRMebTwHeAGdbaD/v6bGNjk7fFnaBoLM7Kddto3NNCZHAxa95u7Ff4\nFoSDnFReyPVXjMmL+b2RSDmNjU1+l+E7PYekXHwOkUh5oLfrXr+EGwT8ELj8eOGbDarrGnhqxaYj\nFlXEtvXvL04wGOD8vxzC1PGn5EX4iuQzr4cgZgNDgEXGmO5rN1prt3hcR785taiip9LiMLdcfbaC\nVyRPeP0S7kHgQS/v6QY3TiweWFrAzVeOU/iK5BGthEvT7XNXOLKooqdBpQXMveNiR9sUkcynAE7D\nbff9kda2TkfbHFdZwde/cL6jbYpIdlAAp+iWH/6h39PKjlY4IKTwFclj2g84BbfPXeF4+ALMunyM\n422KSPZQAB/H3fNqHR/zheTQw6zLFMAi+UwB3AenTyvupnFfEQEF8DEtqKl3fKoZJKebKXxFBBTA\nvYrG4o4vsuh285U6u01EkhTAvfjFM/WutDtyWJkWWojIYQrgo1TXNbCvuc3xdgtCAeZ88ULH2xWR\n7KV5wD24sb9Dtzs+N96VdkUke6kH3KW6rsG18NXQg4j0RgFM8qXb4hWbXGm7pDCkoQcR6ZUCGJj3\ndNSVdksKQ/z0zmmutC0i2S/vA3hBTb0rK93CoYDCV0T6lNcB7NZLt0AALhgTcbxdEckteRvA0Vic\nuugOV9oOBQM6xVhEjitvp6HV1DbQ6cKBpIUDQlQOL9esBxE5rrztAW9tPOB4m8EADD+phKpJoxxv\nW0RyT94GsBtOHVLKddNGq/crIinJ2yEIp0cfAgG456aJzjYqIjktL3vAyaln7Y62OSCcl49SRPoh\n73rA9y5cQ33DbkfbDACVpwx0tE0RyX151W1zI3wBBpUN0Is3EUlb3gRwNBbnLRfCt7ykgC9VnaUX\nbyKStrwZgli5bhv9fe82uGwAg8oKj7imWQ8icqLyogecSCSwW/b0q42Rw8r4UtVZDKsoJhgIMKyi\nWOErIv2S8z3gjs5OfvWsZe+BQyfcRmRw0eEtJRW4IuKUnA7g1rYO5i9Zzxsbd/Wrne/fOtmhikRE\n/ixnA/jAwTZ+/MQ6Nm7dC0BBOEhbe2fa7YwcVuZ0aSIiQI4G8O6mVu5f9Abvd+33MLhsAHv2pz8E\nMai0QKdZiIhrci6At3/YzH0L3yC+7yAAwyqKaWpJP3wLQgHm3nGx0+WJiByWU7MgYtv28d1frz4c\nvmeNqqD5YBvNB9M78SIUhPn/NN2NEkVEDsuZHvD6hg/56eI3aT2UDNtPjB1KU8shmlrS3/PhoW9c\n6nR5IiIfkRM94FVv7eBHi9YeDt9LLziNKecOZ8Pm9Of+3jV7gtPliYj0KusD+IXVW5m/ZD0dncl1\nbtdOPYPrPzWGx56zabc1UwsrRMRDWTsEkUgkePqlGM/UNgDJHclu+LRh+vmnARDf25pWezOnjebK\nSZXOFiki0oesDODOzgS/fs6y4o0PgOQR8F+56mw+Pnbo4c+ks+/DXbMnqOcrIp7LugBua+/gwaX1\nrH67EYCiASHuuG48Z42qOOJz4VCA9o6+Y3hQaYGmmomIb7IqgJsPtvPTxevY0LWxzsCSAu6cdR6j\nhpd/5LMTxw3j5Te399qOOX3BAOLmAAAFwUlEQVQQVZMr1esVEV9lTQDv3d/K3EVr2bJzPwBDBhXx\ntS+cx7CKkl4/f1PVOABee2snbR2dFISCfOKsoYevi4j4LSsCeOfuZu57/A0a9yQXWJw+tIy7Zk34\nyN68R7upapwCV0QyVsYH8JYdTdy/aC37uraTNKcP5o7rxlNSlPGli4j0KaNTbMPm3TyweB0trckF\nFheMiXDL1eMoCId8rkxEpP8yOoDvX7SW9o7kFpKXTDiVGz9tCAYDPlclIuKMjA7g7vC9cnIln734\nDAIBha+I5I6MDuARkTIu/dhpfPK80/wuRUTEcRkdwPfcpM3QRSR3Zf1mPCIi2UoBLCLiE8+HIIwx\nc4GLSO6X84/W2te8rkFEJBN42gM2xkwD/tJaOwm4CfiJl/cXEckkXg9BXAY8DWCtfQuoMMYM9LgG\nEZGM4PUQxHBgdY/fG7uu7evtwxUVJYRzfNVbJPLRndzykZ5Dkp5DUr48B7+nofW5smL37mav6vBF\nJFJOY2OT32X4Ts8hSc8hKRefw7H+QfF6COIDkj3ebqcC2zyuQUQkI3gdwM8BnwMwxlwAfGCtza1/\n6kREUuRpAFtra4HVxphakjMg/sHL+4uIZBLPx4Cttf/s9T1FRDKRVsKJiPgkkEikc4C7iIg4RT1g\nERGfKIBFRHyiABYR8YkCWETEJwpgERGfKIBFRHyiABYR8Ynfu6HlLWPMD4CLSf5/8D1r7WKfS/KN\nMaYYiALfsdY+4nM5vjDGXA98A2gH/s1aW+NzSZ4zxpQBjwIVQCHwbWvt7/ytyl3qAfvAGDMdOKfr\nZJAZwI98Lslv3wQ+9LsIvxhjTga+BUwFrgSu8bci3/w9YK2100lu2vVjf8txn3rA/ngRWNX18x6g\n1BgTstZ2+FiTL4wxY4FxQN71+Hq4HFjetTNgE/AVn+vxyy5gfNfPFV2/5zT1gH1gre2w1h7o+vUm\n4Lf5GL5d7gPu8rsIn1UCJcaYpcaYl4wxl/ldkB+stQuBkcaYjSQ7KV/3uSTXKYB9ZIy5hmQA3+53\nLX4wxtwI1FlrY37X4rMAcDIwk+R/hj9sjOnztJhcZIy5AdhirT0TuBT4qc8luU4B7BNjzKeB/wv8\nlbV2r9/1+KQKuMYY8wpwM/CvxpjLfa7JDzuAWmttu7X2XZLDEBGfa/LDFOB3ANbatcCpxpicPhRS\nY8A+MMYMAn4IXG6tzduXT9ba2d0/G2PmAA3W2uX+VeSb54BHjDHfJzn2WUYejH/2YiMwEXjSGDMK\n2J/rQ3MKYH/MBoYAi4wx3ddutNZu8a8k8Yu19n1jzBPAK12X7rDWdvpZk0/mA780xqwgmU23+lyP\n67QfsIiITzQGLCLiEwWwiIhPFMAiIj5RAIuI+EQBLCLiEwWw5AVjzHnGmAf8rkOkJ01DExHxiRZi\nSM4wxsSA8621e4wxi4AD1tovGmOGA8uBPdbaqcaYP3b9PhkYA3zLWvtfxpjRwK+BBMnd6j5DcnvI\nIuBBoBUoAe7Jx/16xXkagpBc8gIwtWsjm2HA6K7r00mGZ09l1trPkNwM6Rtd1+4BHrfWTiW5PHhM\n1/UvA0u69qm9iuTGOSL9pgCWXPI8cAlwLrAB2G6MOZ1kAP+foz77x67/3Qyc1PXzed3XrbXLgP1d\n158EbjHG/CfwcZK9ZJF+UwBLLllOcket6cAK4CVgGnARH/273t7j5+6tH4NAzz0YOgGstS8C5wDP\nktwu8r8crlvylAJYcoa1Nk7y7/RnSPZkXyS58dE2oCWFJjaQHBfGGPMpoLzr5zuAEdbaZ0gOWUx0\nunbJTwpgyTV/BCqttR8AbwKTSI7npuJbwD8YY/5Ashe9lWRPeQPwm67rNcA/O1205CdNQxPpYoz5\nOFBkrV1pjBlGMniHWmvbfC5NcpSmoYn82X7gx117NA8AblH4ipvUAxYR8YnGgEVEfKIAFhHxiQJY\nRMQnCmAREZ8ogEVEfPL/ATmT0Rrt4ZoVAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fda2b341dd0>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "xe4ewK6vYel2",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Here is where Statisticians introduced the concept of [cointegration](https://en.wikipedia.org/wiki/Cointegration) which a soothingly complicated name to say that if you can get a mean reverting series by linearly combining other series then you can call them cointegrated. The small print: the series have to be [Integrated of Order 1 I(1)](https://en.wikipedia.org/wiki/Order_of_integration) but luckily the ADF test checks for that. Some practicioners might skip this test (as it is not readily available in Excel or Bloomberg) and will 'visually' inspect the timeseries - if they look 'trending' then they will assume they are I(1).\n",
        "\n",
        "When you compare only two series, you can combine the linear regression plot + ADF test and call it [Engle Granger two step](https://en.wikipedia.org/wiki/Cointegration#Engle%E2%80%93Granger_two-step_method).\n",
        "\n",
        "\n",
        "\n",
        "### Using PCA\n",
        "\n",
        "As we saw above, linear regression can give slightly (or in some cases, big) hedge ratios (in cointegration you do not define one explanatory variable but use the set as a cointegrated set).\n",
        "\n",
        "That is why some practitioners have chosen the Principal Component Analysis (PCA) methods (as it is symmetrical) - the hedge ratio results are more 'intuitive' (as they do not depend on which asset you put first). \n",
        "\n",
        "Once PCA was found for pairs trading hedge ratios, it was normal that people would apply it to multiple simultaneous assets (like the term structure of rates!). \n",
        "\n",
        "I will briefly mention how PCA is used by some traders, but keep in mind that there is a 'proper' method (the [Johansen test](https://en.wikipedia.org/wiki/Johansen_test)) finely tuned to test for mean reversion of linear combinations of Time Series - I will explain it later on.\n",
        "\n",
        "For two assets, the PCA method which gives a linear fit solution that is symmetrical is"
      ]
    },
    {
      "metadata": {
        "id": "1tD9pfJdlp6N",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "6b4b889d-16d5-4feb-9453-13fd5ad7be26"
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "wings_body_pca = wings_body[['body', 'wings']]\n",
        "\n",
        "# Demean the time series (we only want the 'hedge ratio')\n",
        "TS = wings_body_pca - wings_body_pca.mean()\n",
        "\n",
        "# Trick: compute the covariance matrix on the levels\n",
        "# by doing this we are finding the linear combination which \n",
        "# return the signal with the lowest variance.\n",
        "# Many Statisticians would be horrified, as we are not doing \n",
        "# it on vectors of Normal Distribution. More on this later,\n",
        "# but it turns out this method fits a line with \n",
        "# symmetric errors.\n",
        "TScov = TS.cov()\n",
        "eig_val_sc, eig_vec_sc = np.linalg.eig(TScov)\n",
        "# Make a list of (eigenvalue, eigenvector) tuples\n",
        "eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))]\n",
        "\n",
        "# Sort the (eigenvalue, eigenvector) tuples from high to low\n",
        "eig_pairs.sort(key=lambda x: x[0], reverse=True)\n",
        "\n",
        "# Select the eigenvector with the lowest eigenvalue\n",
        "# which corresponds to a signal with low variance\n",
        "eig_2 = eig_pairs[1][1]\n",
        "\n",
        "print (\"beta_pca :\",eig_2/np.min(eig_2))"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "('beta_pca :', array([ 1.        , -1.02368157]))\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "owSB4sn8oi42",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "The 'symmetric' hedge ratio of 1.0238 is quite similar to the ones above, but at least does not bring uncertainty. \n",
        "\n",
        "Think of PCA as a method to find the axis of an ellipsoid that can fit the data (read [Wikipedia PCA](https://en.wikipedia.org/wiki/Principal_component_analysis#Intuition)) - it will work even if the data does not have a linear relationship and we can choose the axis (or 'eigenvector') with the lowest associated variance as a mean reverting candidate.\n",
        "\n",
        "Once the method was derived, it was not long before someone thought that the PCA could be extended to many assets. Instead of arbitrarily using the average of the rates they decided to use it for 3 (or even more) assets..."
      ]
    },
    {
      "metadata": {
        "id": "Kj0wLmhQpaco",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 302
        },
        "outputId": "f4c8218f-333f-4b57-e77a-58008ab56d1b"
      },
      "cell_type": "code",
      "source": [
        "\n",
        "\n",
        "# Demean the time series (we only want the 'hedge ratio')\n",
        "TS = TimeSeries - TimeSeries.mean()\n",
        "TScov = TS.cov()\n",
        "eig_val_sc, eig_vec_sc = np.linalg.eig(TScov)\n",
        "eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))]\n",
        "eig_pairs.sort(key=lambda x: x[0], reverse=True)\n",
        "\n",
        "# Select the eigenvector with the lowest eigenvalue\n",
        "# which corresponds to a signal with low variance\n",
        "eig_3 = eig_pairs[2][1]\n",
        "\n",
        "# I normalize by the middle weights (5yr) and \n",
        "# multiply by 2\n",
        "pca_weights = 2 * eig_3/eig_3[1]\n",
        "print (\"PCA Weights :\",pca_weights)\n",
        "butterfly_ts['pca'] = (TS.dot(pca_weights)).dropna()\n",
        "butterfly_ts.plot()\n",
        "print (ts.adfuller(butterfly_ts['pca'], 1))"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "('PCA Weights :', array([-0.83518678,  2.        , -1.25708147]))\n",
            "(-5.08210730611698, 1.5228083481505643e-05, 1, 7096, {'5%': -2.8619473981777572, '1%': -3.4312718809379334, '10%': -2.5669868539820775}, -30116.630501359672)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD7CAYAAACG50QgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXd4HNXVxt+Z7eptZVu2JbmOe7dx\nAWxjSkJLCCT5AoEQICSE0FLAQBoQakIChA9CDaE7dPhoptiAey+y5XG3ZfVmSStp68z3x+zs9NnZ\nppW08+Pxw+7Uq9mZM/eee857CJZlYWJiYmKSGZDpboCJiYmJSd9hGn0TExOTDMI0+iYmJiYZhGn0\nTUxMTDII0+ibmJiYZBCm0TcxMTHJIKzpbkA0mpu7JDGlhYVZaG/vSVdzYsJsa2ow25oazLamhnS1\n1e3OJdSWD7ievtVqSXcTDGO2NTWYbU0NZltTQ39r64Az+iYmJiYm8WMafRMTE5MMwjT6JiYmJhmE\nafRNTExMMgjT6JuYmJhkEKbRNzExMckgTKMfJz3eIP7w7EZs3tekuv7Vz/fj3n9v7ONWmZiYmOhj\nGv042bq/CbUt3Xjy3Socb+xSrP98ywlsqGpIQ8tMTExMtDGNfhL48783p7sJJiYmJoYwjb6JiYlJ\nBpGQ9g5FUf8AMB8AC+AmmqY3h5cPB/CKaNPRAJYDsAO4B8Ch8PLPaJq+N5E2pAsCqrIWCliWBUEY\n29bExMQk1cRt9CmKWgxgHE3TCyiKmgjgeQALAICm6VoAS8LbWQGsBvA+gEsArKBp+reJNTv9GLXj\nDMvCYhp9ExOTfkIi7p1lAN4FAJqmqwEUUhSVp7LdlQDeomnak8C5BiwMYxaeNzEx6T8k4t4ZCmCr\n6HtzeFmnbLtrAJwt+r6YoqhPANgA/Jam6e0JtCFt1DQZe4cFQyxs/V7A2sTEJFNIpjlS+DAoiloA\nYB9N0/yLYAOAZpqmPwyvexHAVL2DFhZmKaRJ3e7c5LQ4ThiGxcrNNZJl4jZ9vul45HNRUTZysux9\n1rZESPd1jQWjbf180zFMHFWM4e6cFLdIm8F4XfsDZlvjIxGjXweuZ89TBqBets35AD7nv9A0vQ/A\nvvDn9RRFuSmKstA0HdI6ibz4gNudi+ZmZVx8X7JFJSFL3KZHVwiDl8bmLvQOAKPfH66rUYy2tbal\nG4+u2AEAeH75GaluliqD8br2B8y2GjuvGon49FeCm5gFRVGzANTRNC3/y+YC2Ml/oSjqVoqifhT+\nPAVcr1/T4PdXunr8hrcNhUyffroIBpnI515fMI0tMTHpP8Rt9GmaXgdgK0VR6wA8BuB6iqKupCjq\nItFmwwCIu8WvAriWoqivADwF4Op4z59OYpmbDQQH3Dtt0GC3Cbf33f/ZgvfXHgHLmi9hk8wmIZ8+\nTdPLZYt2ytZPlX0/AWBpIufsD7y5+lD0jcK8/Nl+/PoHM1LYGhMjNLb14N1vjoAaWQCqvDDdzTEx\nSRtmRq4Kvb4gGnUKGfsCxnvvVYfbktEkkxho6eiFpzcAtU5980lv3zfIxKQfYQYTqvDH5zahtdOL\nm78/DdPGlKS7OSYxcuuT6wEAd189T7HOae9fRapNTPoas6evQmsn1xt85I1dWL29Ns2tMYkXtZ5+\nXnb/j6QyMUklptGPwouf0pLv7V0+1e2qjrT2RXNMYkBt0tacyDXJdEyjHyO1zeqZuLzvnjGNSr9B\nrcCNKYthkumYRj8GOjw+/P2/O1XX8bZeHpdfnOdMdbNMNPhw/THFMkZlOxOTTMI0+jFwsFYqKzSx\nQgj96+jm3D6e3oBkm+I8R+obZmIY071jkumYRj8K5y2o0Fz3rVPKI5+rj7UDAA7Vdki2MXuW/QvT\n5hvnjqc3KOa0TAY+ptGPgsOmHeIXEKX5d/VwPXz5Q8KaPuR+hdnTNwbDsmho68mY6LV270n0BLRz\ncwYTptGXEQxJ++ZiE2GxSIVED9VJe/WA0r1j2vz+BRPj0Kuloxf3vrgFR+rliuGDm0yb8P79uvvw\nu2/+nO5m9Amm0Zex82CL5Pv2/c2Rz1ZSavQZhsXV500EAMwYq0zictgtZjRPH+OPki0da0//3W+O\n4FBdJ558tyqRZg04MsnoZ9rozzT6Mnp9UqNxtEEQDpXXumVZYFhxNgBgaFEWvt5ZJ1nv84dwrKEr\nqiEySR5fyX4DOfHasgyzCxnVWWHYzJp5M42+DL23vrzUbV62HWT4Cm7b34wXPt6nut+GvY2Gz/36\nFwcUow0T4wwpzNJdH2uvLlOrG2dSTz/EZlanzDT6MvxB5VufNxTynj4BwBOewG062at5TK/f2E3V\n3OHFys01ePTNXQZbayLH5dDX1gmEjPfqunr8WFvVEP6WOUYQyCxhuiAjPJ+Z4Ooxjb4MNUGuQ3Xc\nJJ7MpY/pY0skETxaGNXU95qFPhLmYK1ych0ARg3LA2D8BQwAb399OPJ58JsCKf/7zu50N6HPEPf0\nM6HXb6psyuA783YrGen18z55+Yi3rCQbLkf0S2i0elaP1zT6ifLGKvVaB/zL3Os3fo3FkViZ5uZp\n61TXmBqMtHtPRj77Q35YycFtFs2evgzeQM8a744s+2QjV+hcbehXmKuecZvtFN04Bi2GuGdpklys\nFu5Wj6V8pfjnziAXN4DMmsi1W2yRz6/sewuN3UrNpsGEafRlBMNP9/iRBZFlVUeUYmrjR+TrHuf+\nny+IfF6zq95QL76hLTOSQ9KBNZxjEYrTemupq5oMfPa3C52tHc27cffGv6WxNakn7nEMRVH/ADAf\nnLvzJpqmN4vWHQVQA4B3kF1G03St3j79hVB4oo/vGYoRd370Su7ZrSRyXELvoaXDi/98sg/XfXeK\n7rljmWQ00YbIPgkEbWB92ZFllvCETIvOhLucTJjUM0JtswfD3TnpbkbKqO2uVywLhAKwiUYAg4m4\nevoURS0GMI6m6QXgips/prLZt2maXhL+V2twn7QTDA//bVbppWEYVhLG1tHt1zyGWgTQYZXsXTmx\nuB5MpPj8IVQd5moaOCdvgHP6NwAAMqcdsPojv2dnT0DzGCbq/OG5TeluQkrJsroUy+q7jYVZD0Ti\nde8sA/AuANA0XQ2gkKKovBTs0+d09nDGXG70v9pZJ+nplw+JrefTamBibLhb6JnWt3bHdPxM5z+f\n7FPIXhOOHjgmbYRj8jpMGV0MILaJXPHvrTV309LRiyadesom/Z+xBaMUy9hBHK8Vr9EfCqBZ9L05\nvEzMvyiKWkNR1AMURREG90k7/KStTxbad6yhS+LT541IMhFPHp/0aI8kTJTw8y5iCDvnyiEdXkwN\n/14HTnSgpsljyHUj3kauqcRz65PrsfypDfE0eUCx96jy+g4WxHH6PC7r4K2DkazYJHl8yh8BfAKg\nDVzv/mID+6hSWJgFq1UaO+9258bRxNiw2KWXxuWyITdXuBHcJTlw62R/qrWxN8SifKj24MbhEHyI\nubnOPvk7xfT1+RJB3la7jQTk7nrRHXaSFLKc//T8Jtz503mYP2WY7jnsonBcp92qOKdYnE/v2g3k\n68rT3hPsd39Hou3x9AZQ3+JBdq7Sd2/JYpL69/anaxev0a+DtJdeBiAyG0LT9Iv8Z4qiPgIwNdo+\nWrTLhs5udy6am7s0tk4e3R5pRqLPG8DJDsGqnGzvAaGTdKXWxhsfXo1nbl2quU9Hl3D8+qYujChS\n+hpTRV9d12Sg1lZCrQ8h6szfv/YfsFWUI3BsEgBg654GjIniovOJkuVyXFbFOcWjQa1rN1Cva3Ge\nQ+KS9Hi8/ervSPS67myuwlPbXodvzwJcfG6BYn1TaweKkZy/N133gNaLJl73zkoAlwAARVGzANTR\nNN0V/p5PUdSnFEXZw9suBlClt09/odsrDOEZFhhTJvTKCYKQDPct8vRcEQ6VrF4gerhgMMiKPpuR\nPLHQ2qmUDWCDdsl365DjXGQPEHO2VZFK2Uuxu09ePGegM3fiEMn3FV8eTFNLUsMzu18CYfPDOuIA\ndh9Ral0F2cGbKBmX0adpeh2ArRRFrQMXhXM9RVFXUhR1EU3THQA+ArCBoqi14Hz3b6rtk5w/IXmI\nY7GzHFapYSCAjSLhNFLH6MvnA4widhfkZA3OcLG+5NsLRiqWkU5ugjwQkL5U5XUUAOlErtocgHjZ\nvS9tjbeZ/RK1u7vqSGuftyNV8BO11pI6HHesVaxX8/MPFuL26dM0vVy2aKdo3aMAHjWwT7+CFAmq\nzRhXgi+3nYh8JwBsPyD0CPR6+tFgWFZyLh6xjs97a45EJh9NYkEwxCPc2YBHttbHzcN0+4RRXUNb\nD+54egMuOm0ULlgkRHKIjbqa6uRgztJVy8jdsKcRU0Zlxj25oX4zZpZOTXczUoKZkStCXB3JbiOR\nLUqwIkBg2ewRke9qRlvO9HHSwirfOXUUOjw+3PzYGny88Zhie3Fv83BdZlVqShqkcA31xLMsot9v\n1yGuB/vON0ck24jNnlqwz2CWH1b9ezMoWa2qVV0mfTBgGn0Rn28RevYWksQV51CR7wQh1dMhda7c\nxYtHAwBmjC+VLM9yWnGwtgOe3oCqMFgmPVTJZjYVDnclBKPvV5toJ5Qy2Zrvb9Hv0djeg5se+wZ7\nRKGLjYM4Pn/l5hoQri7AJsyVDOaXnBr+UAChQejmMY2+CLmRkEzeEfKJWO2efkk+F3Ujl2l+7fMD\nuuffVC0VerrqgS/xzAd7dfcx4cjiwytJ4Tfs9KtNrnK/4ahhQmSDkVHbSY8fXT0BPPz6Dix/aj0C\nQQZf79Cv0jXQcU5dC9fM1ZHvg8XoN7b3gOmJnlx5y1d34k/rH+yDFvUtptEXccasEZrrQiGpDIOe\nnZg1nnPrLJtbrrI2trmA9Xsaom9kEoEQ9fRXNr2vWJ8VHq2J7ZfW9IyWiWtq78WJZo+iqM7gQvnX\nb6GbcayhXwXcxcXbm3aCzPJE3xBAu+9k9I00ON51Aq/RbyMYMhYJdLzzBFp6U58EZxp9EftrtH9g\nrz8k6enrPe82qwVBJgjSoowIUXU5APhsS43m8Uzhr+hEJsFJ/VDXs+ZyL/aQOFpH48fUu+wkQUhC\nfItVQjoHNup//L/e39PH7Ug+B3yxRVr93+FP4zrPg5sfw5raDfj6mDHtoge3PIY/rX8grnPFgmn0\nRWzep62j7Q+EpD39KD32O9fei8vfukmxfBvdrLK1vutHLQbdREpEoZTQN/p8r17rBU4fb4981nvZ\nWiwE7DbBfSfWTRoUaNzerR29+Odbu9DSYVyttL/Rw8YWJPHx0S+wpWF73Ofb1bAX1395KzbrHKMv\ni7ObRj+Mlr9yQjmXrVc+NBchsRGIMrL3BNQF01riMOCrtw9u33EyyMr3wj5hEwi7vrAdoWL0xT79\nB1819nATkFY6G3yT8Op/TzDEYvuBFry8cn8ftyd5WHLVR/RMTw58e09RXffvva/Ffb51NdzI4gWd\nY1S3CdczxITw9K7/oCFFxVxMox9GK1v2vIWVALhen7Snb4wfLRsn+R6PT/TzLTVYsyuqYkVGswVv\nwJLXBgelP3Qnwnd8vMVUeP7w3CbsPiwkKw22DOrZE0p01/cM0HrOtc3qvvxA7Rj4qk4F09P3Gjks\ny+KJnc9Hvt+4+nbsbNmDe1JUzMU0+mHEBn3pzOGRz9awP0A8kXv/z+cbnsQ7a+5IPP27JbrbfLH1\nhO56f5DB8x9V4+lB4E9NN0JPXzDSBiI2o+ILDK7QPjaKu+HgiYEpO3H/y9tUlzOdRdwHtu9NYptX\ney6xXWddvJhGP4zYCJy3oCLy2RKuoPXBuqNoDlddssQQtcGwjGoVLp7Obj9e+czYUHnD3saI3r8J\n0NXjxwfrjqI3hl4nn18hKVgj+zn/b91RtHZ4Y1JUP1I/8KNaxDAYXCMXAPhsc43mCIXxhEXXWP1n\n++DJI7hz7b2o88QXVdcTUOZ2vH/4Y83t/Uzyi/6YRj+MeLgvjs/na6sCwL7j3FtXT3cHANq8wmQg\nP0FzyZIxqts++2Fscfirt9XGtP1g5sVPabzz9WG8/VX0gvJFTq68pdyn7/OH8NGG45Jt3/76MP7x\nxk5JdE6m0UMqRcgGMscbu/DaFzp5MpEevvazvbVxB/6x7Umc9HXgpeoVcbWjoUcZyBHSGVWlYoLX\nNPpheJG0iRXS2rdWldRbXndnf/sh3LTqdkG5MUxvUJis5X9Ql4byZk2j1MdYOTRXs0oTABxvMhZf\nPBA50eTBzoPGjU1LB3edmw1EkvDRVhGjH+7pv7f2CBpVCtLXtXQjL8uuWK7HYJrMrclele4mJJUP\n1h2NfGZ8auG10UfvK+h3I5+Pd8XX+cqxKaO8GnRKM5a4kq91ZBr9MJ9u4uLkq4+1S5ZbLMqbge/p\nP7r9KQTZEOwV1ZL1/pDgguHTuC0aLh5xrd28bDvuuHy2JBRQzrb96iGfg4E/Pr8Jj765y7Dx5F++\nvHaOHvwcDBGWYeDdeWoGn8flMKZHOCIcrhmvumq/xIAHc9ehVry35siAyNSVJlbG197uYOKyG4yK\nHpSe1IONTFadKwHT6IfR6i2qFSuXK2ySOdJJLbEsKy/6pWYQ5LHO2U4rrBYS7AB4iJKNeE7FaDKa\n2M0WbNWvvDkypwwAwCAUPp/+OYgY2nGimQvPXbt7kEZYWdTdXI+8sRPvrTmCNQPg7+YVcsncVtWw\n3rkTBJ2sQK26KzYWAiH1a1Yn69WzLIum3r51pZlGP4xmco1Kj4cgCElvHgAWThGMjljdkffJqRmZ\nW59cL/nO90b4nq64iMtg5+kP9sA6/ADI/CbDUTOSrNqQeo/ob6ffjRtm/Ay5dk5r5aMTH3Kb89da\nw/jzSwlHDxxT1oDMaVfdLj9bcAHVDFLXm2v2F7rr11f1b6kQsXqtY+Jm1W2GFQulT4O141S3iYXl\na+5RXf5K9ZuS7//a9ULkM1U4NuHzGsE0+mFGlHBG4fuyCVebilvGZiUV2XVLZnBhnqOG5Ul6iPwL\ngDUQC9LYzvX8eYOkF/Uz2Nh6Yh9sww/BQW1DMGTMTVIkqlksFloT47I6MaFoHE54uAQ3P8O9rPlr\nrFY8hYdlAWvZQZBZHtjG7FTdpjhfaIPeXMxAh3BouzaaTvbv7Fy10ToA/G7GryOfk62j5A2pJ2Hy\n9x9PVavgGh5vGv2+hU/jz3FJK1ap+XUDrA+v0m9Jlo0dkY9bfzQTt/xgOh7f+WxkeWM354MfWqRd\nRF1O5VAuQWTM8PzIMv6ezHZa0eHRzzodiBA24WH48MhKQ/tIfP8h/UpjZ1dwtYmnFk3hNg//3nuO\nqvfgjUIQwLnzuRDfkaXRlRsHKo5JGzTXiTOT+yNac0Ru0SSp3JUnV+FkGWMvhQAThDeo/XzqReNs\naRQ6kguHzcVtc240dM5YiXuWgKKofwCYD24kfBNN05tF65YCuB9ACAAN4BoApwN4AwCfYbSbpukb\n4j1/MvH0BvDCx1zRBHmvJS9bGcGxSUNDY4Is8gcAHt/5LP73jIdQkGO8F3jVeROx82AL5k8aio82\ncMVWrvvOFDzxbhW6vUHc8vhaPHfb0sGl8ijSzNnavB0XU+dF3UVcyYwN6hv9LCv30vUz3ANpJCOX\nYdmoE5okQaAprKv/wsf7MJsq1d9hgELY/CAc3WB9SjeoXmIay7Lo7AlI3GB9Df9bEy5pLoXVQsJb\ntQAgWLDTpfv4qhbBVrkH1lL9xEkxnf4u3K7h1hFz06rbYSWteHixdNsCRz6spBVFzkJcNvH7hs8b\nK3EZfYqiFgMYR9P0AoqiJgJ4HsAC0SZPA1hK0/QJiqLeAPAtAD0AvqJp+pJEG51MgiEGNz76TeT7\nCQN+WSuhHV2jxahhxv3z2U4bFk4ZJlmWK6uZ6wuE4LQnf2Y/bYjUMQNM8nuO/AuSPnkAwBh0x9w7\nVbf+BEHgQLgoeuzHHFhYyw4hcGRaTPus3FyDFV8exC+/OwVzJqTnhcgwLGDzwjlVWguXJAG2hxtN\nK0cDBJiuQoA3+gSQZXWhJ6h0ZTV0N8FGWvH6/ncMtSfIhhAMhRS9/iFZblw/7uqUd+bide8sA/Au\nANA0XQ2gkKIosVWbTdM0/4psBtBvC2sGZJopbeHi6HqRGyQR32W7/Ozxce0HKH2OiWrH9DtEYXQ9\nwR50+hPPcD2n4ozI52HZQyKf7WO3YX9dE/Ye1dcu77W0wFoSFrsLN2/B5CGSbeR5ep7eQZzQFYct\n+nond/02VWvHoqeaEMOCdCoFEMXRX2ov7FBrGXzVc8F05wFgMTRb+tuzLAuWZXHPxr/hj+sfwN5W\nOqZ2PVf1suS73WLvk9F7vEZ/KDhjztMcXgYAoGm6EwAoihoG4GwAH4VXTaIo6n2KotZQFHVWnOdO\nKvJrnJdtx1sHPsDvvvkTgho9znjNbSx+fZ4bLp6Km78/TVHdqbVjcMkt20ZKa5L+c/szUfex2/Rv\n3zEFlcK2pDBSshQ1wTr0GLbsa4K7QFsH/6RTkMcgnb3Iz7bjR2eOx/kLheMSBAGvT3Bv1Leqq6sO\nBsic2HVg+GCEoMZkal/Asiwspcp6FeJnimHUfO0EmK5isCELCAJo90gnsx/a8k8sX3N33O3a0Vwl\n+X5m+eK4jxULyfIPKF5PFEWVAvgAwC9pmm6lKOoAgLsA/BfAaACrKIoaS9O0rphMYWEWrFapO8Xt\nTp4SXo8s1f6K8ybhz5teBQCwLh/ceUo/fXtIGVer1yZ+XXtvtOE/ixWH38KsYVOwsHwOAODs8L57\nDksTkGpaezF7SlmU48VGMq9rrBCygjN13Q1wu3Nxz+pHMDx3GK6a/UPJerc7F2UlOTgaKWavNCo5\nuXbtv4lg4XLZ0XxS++Upz8f4zuIxGFVehJ+XF+GLrTXo9QVxdMhrQIET2LkEAODKcijOmc7rGit8\nW7O9Feh2HpOsI53aETxaf6MrXKmMtJBJvw5Gj8da1d2xpaV5GDsiH+7CLLhcOnMOYT2e9pBU6vh4\nV3R/v//QVDDd+YAlBOfk9brbVpYN0V2fLOI1+nUQ9ewBlAGIZGiEXT0fA7iTpumVAEDTdC0AXrDi\nEEVRDQCGAziid6J2WfFptzsXzc3JE7eSG/2SbKFHePPHd+EvC+/AH34yB/f8Z0tk+UcHlCnqem3i\n1xU4LThzzghUHW5Dg0omKOHsxtdHN+LroxsxzkVJ1tU3Sgs/dHb2xnwdnnp/D2qbPbj7aqVmeLKv\nazJobu7C7kYauxtpXFB+bmQ531Yv2QrHpC3wH1L3M1v9Lt2/yevTd8XIXWjzJ7gjx2NZRMJESYcX\n3EuHQENTF5qLXIq2DgTEbSUDWUAMxcC0/sa28Ih0G92EW/6+GndcPjvhdgKxXdfWjl6QWcptm5u7\ncMePufa8rqPLY8mXugHz2GHoJIwlpIU63EDQDiP+gWTfJ1ovxXjdOysBXAIAFEXNAlBH07S4xQ8D\n+AdN05/wCyiKuoyiqN+GPw8FMARA2tXD5K5xuU/t9+vuk03CxjZMFVfYIggCl545HotnaPTQNWLN\n1doVjMOnv3FvYyR7dDDgKdkKMqcDtvJ9irHmrXNuwPCcYeo7AiDsXmEXqw+O6atBFgpJRjYrqfip\nxRPncregfRwX0dU1SFRQk+WMaRG5IQ/WpkeOORAKgHTp3/dil1002juj/8aR5z6i2qnvq7eTNnT1\n+PG7J9Zi497Uzn/EZfRpml4HYCtFUesAPAbgeoqirqQo6iKKorIAXAHgGoqiVof/XQvgfQCLKYr6\nBsB7AK6L5trpC2KqP2v1wzVPv16mfEa+Mk9ZHP306UqjP3NcCS47Wzs5QyHsmYC4l3zyWkzV4Vbc\n+9IWrKuqx30vb4XXn76IFD1NEgAIBMLXQEVLpSJvpP7Bg1YQBIG8LBssxfUgHV44xu2IrLaQRHTD\nJzqvpZAb+vPujIFGVUs1TvaKjbLGPaIhyaAFX3kunYQMPCvy/BxdopTkBETJmCKpZjagfY7lc2/C\n5n1NaO304akU182I+w6laXq5bJE4ZVErKP2CeM+XKuQdZt2ondzoiTxymdQ2rzJCRC3h64aLp2F/\n+yHOcaaCvKc/3G08EWjLvias3yP0Yg/WdijURAFg16EWPPLGLgDAoVrOnbTzYCtOmdQ3vkY5vpB+\nn4Dl5XAJBmSufiSOYt+QDQQAQkMmW34XjMqrkC0hVF82Wtmf/Zk6TwOe3PVv5NO5uG/RHwBo6+nb\nKvYicHi66jo1spwxGNMk8ObqQ7BbSVx46qjIMlZF5OzXs36pWPY/Z4zFwdoObAnXsb7j8tl4/sNq\nyMcnlrwYEvpERVm8O5bANfczxSY3zrgWQ7JLcfvKL40fNwEyPiNXbuSPdUln+W1kbDet/Hgd/i7U\nRJFh/W74BhX3bOWCTXLbFEvloiferZIkMr3ztbr+PG/wxVhVVEZTgdOrdMX0RFM1DOvtENaAZt1T\nMRZRfgXh6MWhug7FgIlwhN0ALKDZ2wU3WCesypdSf1dXfufrw/hsi/Qe58NjO3yCh7YnV93HbcmP\nrmiaTj7acAzvrpFOE6rp1Ysju3jOnleOn10wOfJ97PB8XH3+xMQaJC7KolGVa3yhUuDtmQ/24kh9\nbAXcjWIafdlD2h2QJl+cE07f5zaOfryASqWbT48Kb3CGZfDmgfdB5goPD98rCbKCK8Un0+iQ9/Q/\n2SQt/BELsfhWtSShk43Vqny5vL1fqCikNgJjQ5wRJ7ONTYD9af6tkc+MJ5+rdsWykgfTOZ1L1GPD\n/4nOhhATipSv6/EFYR16VHGOVD2oyeKDdUfx2udSgx5L3olYLkOM1ghZLhnS1tn3ocbq4ZjqyCO2\nxpTlw9U8I4GzExqfRUsJQnH91u9pwH0v6dd7jpeMNPqfbDyOLfs4H2x7l3BTjh+Rj6DMaBPhB2LB\nZKV076UTLo583tXM+eGOdBxTbLe9eTdqPdxs//72Q1hVs0ZV7U9cfEXe0092zsaLn9KGhLL2HUtM\nm0aN99YcUWrgq7hKDjYLvq4Of+LGtNhViGun/kRyvqHFGuqqkIrksQCe2Pk8fr/uPrT2cteEzFeG\n7q7anvbYhJiJxejLNWl4tEY4ciN657MbDZ8rWQRjqD5FEEBJvlPyvLs80p64b/+suNsyI18ZOQeo\n6wOlKgEzI43+f1cdxBPvcomoIlZ0AAAgAElEQVQRr3wmZNE1tPUoJACOdXJDYc4nLr2BxT7n/ScP\nAZDKKot5vuqV8D7aYkz/2ft65LNculktU++hV9WLPBth9fZa/O/buyPftTJJV26uwRurDsKfhMLf\na3bVY31VA95bcwSPvCFVrWRVXCniB8HjV0ZfWLNiL2phJbnRARGWfbBZSaj2wFjA0yv8Bv6QH/va\nuR5yQ08j7OO3hkM1OYhg4gqbLMuiobspJSXy9Iipp29Vv0+MFr5JR6GZWII1CILAg79YgJ9dMCmy\nTF49jw3E/1sf3lqGh0+/G3l2aTil1lxQSwoUTDPS6IsRG7vOnoCih72rRXsmvVFU75K/sbRu/oYe\nPrHDWJddLkMg7zEBXM3eQDD+h0ii/67zYHy88Ti+2GpceEqL5z+qxjP/p14TuKNdX89o9Ym1yoWu\n2Hv/Eb8+wWDmuBJUH21XuO3IvBawZFCiB1TXLUyEP7HzeVgKpBXMSp3aRVxuf3o9rnrg86jGZ339\nFtyz8W94Y/97Bv+axKg+2oZ1VfWwqBh9W5d69JNaARKgf7u15M8kEeUZlHewxNXz2IBNt3j68rk3\n6x675aQXTqsTP5t6hWT58x9Vq27v04m0i5eMN/riVOzh7mxVnzzDMjhvgTx6Q5raz7tv9HppRzuV\nfnheRvmbWmm2XlAWrqilyREIJmkIGMV/lIrYYU9vACs312iqNIYI4bdYX69e/CJWBKPPoqGtBwzL\ngJBlmjombAFZucNQaB5PY1Bw68lHRZ1j34Fr3kocbdA3jK+F5bq/rtXP3EwEMq8lUtP5r6/vwLP/\nV61qBI3UfxBz/8vqo87+MK8d698iR6zRQ9gCIOzqve97Ft6Okbn6WfK8y2Z0fgUWlc3Dr6ZfAwDY\nVN2kun1MoaQGyRijf7C2Q1H/FpDelHdfNU9VC/vQyaO4ePEY/PTb0pl8sYE/cPIwOnxdWH1ijWYb\n2r0dCg0d3nf4Oi1V6BNP6gLaNjkYwyRVIugVZG/sbkJPQH8YyrIsLO4akPlCD/mxt3bh9S8O4H1R\ntEWwaYTwmdA+ZrzupsjLk2BR39oDS0ktbMOOKrYjC5piMvoc3N2kpbZ54KR61BRPX7h1HBO2wDl5\ng3Q+R827pTMqsY3ZoWn4+htH6juxv0Ya2SXvZUejSZQ9z/TkRFyDcoqcyjBoPS6dcAkmFo/XHa3n\nZZlGP27ue2kr/vradsnELSA1HgRBwGVzyXdFm5d7Wcg18bOs0m3vWHsPDocncstzR0CO3OD/7ZcL\nceYc5XaAsqdflKeeE59oXDj/QLTHGVXRG/Ti7o1/w90b/6q7Xae3G/ZRe+CghIgEPuz0i20nwBvM\nQK2QoMYQ2olAfJ2BWBH819z5LAXqPSyWhebDrUl4cnhdlZCiLy7H+MUx7Q5BX9PVLcxXhFRHi9r3\nlbW4AbYKdTddtKPYrX1rcu75zxZ8sFboVJxdsRTT3ZN19lDS2RMA4+We9WD9qLDqpjZ3LxBSmL4j\nyhfgkSdHfrxROxIvFaqbGWH0xb0Wv+it+sbqgzjpkU6Y8sqa4nqVL1avCCviSX8AqmgcrpjICYHN\nLJXqv6iJMVW3HUC9R3CTFOU5NX/UZ6teknzXKkKRaMbsA69sA328HW+t0tYe0YOfe+jy69ch6OzV\n67UzqtE7WjAMiy+2ngDTqx15o0XElWEkOz6GNgEAYeM6FAdEORRif7LXElsCWbIRPwc9PuG++Xy7\n8gUa1SViUd53LR3S3/hQXYcin0QvGzxliH7jeGXRfdWnIHCcQqi1DKxPqZYrdvUWu4oinxdMGYqn\nbz8T+TnC8/vySloSNt3QGntAQiJkhNHfLVKoFE+IfrxB+oYNMSG8eeB9AMApQ6XCUHtb9ynsw9iC\nUSjNcgMAig0M7b6pXY/3Dgux5/IIHTnr6jZhX5u+MX7ote3w+oNo7fCqyvrKBeXUePDV7fhis1J6\n1ghbG3dE3wjAlzuE49tG7YKlRAhtlFxXjUmyGe6pkc+frD8CXyEdVU9FDWG0xcr+r0ZsRt86lOtR\nisNR/UHh+o/JmqDYJ8gE0eHrm0lQscE92iAECtBe5XxJLBIUPB3d0vv53heVceap9PHLXVLCd2F5\ntElcTQJOBBtGIZyWp1j9+1N+o7qblSQwrCQb3z5FmBP8Zlc97ntpK4IhBh0eH7YdaFbdN1VkhNFv\naBN6IFqhr7/54QxFNq4UdS0WPvIhHn/sLV/9Hv6QHx8c+kR1/Sv73sQ/dzyj61/t8Pjx2Ju78Lsn\n1+HOZ5Qx0N44QuTuvHw2Hr3xVMVytSLievVAxeztEMJDre462EfvBmHvha1iL+ZNKRRFhag/lIUO\noV7wM5vfhG2k8mVopKIZP7IieKOl4cLhNovRQKgYwh6/cH3ybcqOwU2r78Ada/8CT0B4geXYYh/B\nGOH2p4U6t+Ks7NBJZUWraJFGhMrf+mUSIrwSYaOsUMu2/co8imS4S9Sy1MW9ezF8J1Ot/oU/wOCW\nx9dyI90wl545LuH2RWNQG/3Obj+ONXQhJJrs1CoqbpP5GuVl0QgArIph52+iL2u+UayTU5GrDIO7\n5avf45Nj+pobW5t26q7fdzz24hZanDVnJMYMz8eG5nVwzftEojypJj8rVrJcVSP1WTf1NOPxHc+i\nuacVvS7lC9U2ZhesQ46j1r4Flii6Rh0i95BtmK4aty5E+JbntXqSISvAV+WyDhH+Rj5SxxsQ7jdx\nUZ6N9Vuxvk7oYZ/0dqDQwYmTiV8AyaS9S2R4bMJnj0cw4E3tPWGDLzXqdy2QSW2pGP29BhL5UjEx\nyXO0Xhrm/L/vhDsaYvdOHCZP7lqlyoWX9z8W/wWPL31QsU9l17fh3XlaJKO9IEfpnvWrTOCeOWck\nbrt0Jm743lTcfdW8mNtqhEFt9G/+5xrc9cJmvLHqUGSZVmhZiGEVMshiCIJU9XOe9BmXNBhdoAz7\nNAI/Ofz7K+bgzih65Mohbmzn6gonJL17iCt25hi3A2ROO2DzYv0eoSf19Yn1eHb3y3ixekVk2ZsH\n3pecfwX9Lqrb9uN1+m34PcoJct7Qt5KiqBaN9vb6o7upJpdw0VWnlqlnPQLCHASZ3RU9OsdA9E6J\nbRi6A2KfLPcH8DkQm+q2R9Zs61kNgMvveLF6BV7e90ZkXYgNwUYKQnypjuQhHerRO8uf2oCN1Y2K\ne91pkQYx2J1Kg2WkZvN8lcz2ZBAMMejWdGUKf8vM0ikxH/t3P5op+S52EWuVOHQEisH6siPuxGlj\nSxTbyBPVzpg1HAD3Upk53o0RpcZFFWNhUBv9mGBZeHWyZQmoP4jFTvVh3Qy38uaS94SNwmf5ji7L\nw5jh+Th3vvbLQzxBByhDPedHUcwcqXKjOSZthGPSBvSGj72hfgtW7H8H25uVAm3ia8i3uyvgAdOr\ncwNL6gioD7+P1Ed/uVbmjcS9i+7ED6mLNLcRj9YspVH0i1QmK+XkuuySFvNzFSzLGaI1R7Yr9rl7\ngzLS6f1Dn0juLzWRsIQxODFddbhN0XmwyLJSg5ZuRc3nQpXerJxUCdL96flNWLu7AWRBE8j8ZvSo\nhM0Gm0Yo6twaoaxE6m4jCQLe3Ytw41SlUmfkXGHvAp/YpXZXyyWU4ymnGg+m0Q/DQqp9o5zwIVSz\nbYdmK/2hAJBrT15puHV1m9Dh68ITO5/Hia46lBZpp4Hf8Mg3knwEca/khounYl4Uoz93gvrfw0sO\nBEMMXqr+r+b+XtE15CMlaj31sA0/pLWLhJ9doK5q6BH14hx+t3obCRIFjnzdCA1xryzaRDDrje5b\nJwkSHaLsaVsFl1n5wsf7cMOj3+Bkj/Qct6+5R/U4+9oPwCMaMWjVZ9aipqtW1y20qbrRsNEnoDbg\nUpot+5BaLp8hTMXQ6Pd8KEV5JfXhCBjH+G1wUFvR0a3swBkZiRiBJAmwvbkYll2GYIhBiGEUtoEP\npeafP7n7GJBOpgNAY1vf5D4MzIoPBuiMsYIRC+kQdrp7Mr46sQ6NYfkEglAPY9MyMLNLpymybI0w\nu3S6wofPsAw+OvoZ9rTuw55WvoD4tzSPsZVuwoTyAhAEgQ/XC+F4M8e5UX1UJ2yQDOHjmo+wYau6\nKBZZ2IBfvv4v2LQLUkl6q0YjJUInS2Ap4Cbdpo0pAXvcCsIqM3oig+Wzq0c7WAxM5IpDcQlrlElo\nDSlcMfLcC8IijFp8/hDcxGh4WKG9cnkNMd6Q8MLs9HXCZTVWr/Chzf/Esa4aZFuz8NDpf1bd5qsd\ndZCYcoIF4fSADUlNgHUEjQBbJNn2/FFnq/6Wr9JvwTEe6N3E3YtG7PmX22rx47Op6BsmiDg/gr93\nJparj8pjhc/Q/eOzG9EbdtEU5zlx37XzFefnjb67QOnelNNXlcUGbU+/RieDVA0rSUiMeoEjX+Jj\nJUDAI4pFn1KsDL8TM05FI9sIV025THX5mtoNku/nnDJc8xhfbqvFNQ+twr/eq8LnW0+AyOqMhBNS\nFYU4d34F/njlHMV+9rHbsaFJWwXRMW6HavaqmHgyhNmg4BYgAFxL/UJtq6jHMRKZYRX/ptnRQiWV\n5zyrfInku4WwYEx+ZeR76KTUdzu8JL4Rn19FDkQLPuqsW6f+AMuyIPOEF75t9G44p62Ba+ZqyXa2\nsiNotApuh/NGnYWzK5bqXlvCxb3IjIqu9QVqUXoOW/ROgRF4hd7OngACQQaBIIOGth6JS4xhWS64\nU3Tdxo/UryKWk8JJbjFxG32Kov5BUdR6iqLWURQ1V7buTIqiNoXX/8HIPskm1j9s3MgCPLHzeekx\nRD8YAQIr9r8b+T5vaHR51TsX34CLxp6H8YXaZRDViFrqD8B5p6pn8vKwrKDn4ZyyDrZyGpsatoEk\nCFyyZAwqhyqzCvnediJ4vNwIq6W3NaJKGRVZbP6MCmWJSSPFXNSEw/TPS4LxxVD9G8Bpwxfg4dPv\ngdPC7UcQBAqdwsMsv4Zrmj+PrU1htCZyGZZBTVdd1O3UsI8V5hckE7kygvBGesdnVyyFhbTojtqs\nbi5UM6bSoynGbuv7/qw4T4EFFB6xGy+eCj0qDbjHkkFcV4aiqMUAxtE0vQDA1eDq5Ip5DMDFABYB\nOJuiqEkG9kkqscbjyofpgKy+LSH93qUi9Stn+tBJOLN8MUboFOhWw4hbZE/LvqjbyPnP3tdxz8aH\nI98lfkYyObVw/QHuOM9XvRrX/g6L+nxFvoFJwjy7fnq8ApYAYWCyVoyVtMBpdcAfLnJDEiS+P+47\nWieIrT0itCS6Pzy8Eg9sfgSbGrgoNLkUuFZdYU5Wwlh7JJ4RaE9E8liHHlPsF9nf1QXXvE8io4G+\nQuLeAa+Am9pzijOxwSqfY3uUkca3TlF2dlJBvK/DZQDeBQCapqsBFFIUlQcAFEWNBtBG03QNTdMM\ngI/C22vukwqSIVnx3bHnRT6HGAbTSwTNjglFxnvv51ScEdf51fR7ePR8w3o0dDeiqYfrjVpFVbG+\nsyS+cFI5fNSJJxCbe40NWZBHFsMS1ruXdxoZNqRbKP3icRdgaklspe1IV7csckjaHr63K3YJOSz2\ncHu4vzPIBJFjz5a4eMCXUdTQnjeCuHShmM3hDGh+bkeuCqvl4onF3nV7/ZE9BIG66A/UalEBGT6R\nzzl1beT/4mpxqYYLwQYsQ47COXUdAGBbW+JKrSX5Tlx9nvp99uS7VbjqgS9x1QNfggWruGSkRj1m\nnuw+qiccr9EfCkA8m9YcXqa2rgnAsCj7JJ1oPf1o/jVAeMABoMvfJXnAYtHwyLHHlmE5rYQr4DBd\nJeyT5/PjX8V0TDF3bXgIAPDrH3AFrsePLMCFiyvjPp4Y3jDHku7O+h3hfYRrGqyRhgN22U7gxtW3\nax7jjJGnxaerYqCZ9y68M/JZPhLhX77iDoKliBNcs5UZi1hS49mql3D9l7fileo3Ud26H71BL1iW\njdyDWxp3YGfzHkX9B628kfaQcWnsXl8ocl3435GMMTt550Glq1CtWlyqYBgWTocV9orYR8RqLJnJ\nzaH9/oo5GKZTbY3H52eURl/DJo0uy8OtslyAVJKs6B1d2ao49olQWJgFq1U6LHK7o/u+jrXoixid\ns6ASPb6DOBGe8BUf8+yxp0e+jyociSPtNZhcNhadPsGlU1KcC3eOdjv4/Y20Vb7fpcUX4NSxszAi\nbxg+OKwu0RA1a5NguF5sSL334MgF5s8YgQ9mcKOJus4G1e2MEmobAktRI7JyHHC7c2GxGp80s5Zy\nE5EFuc7I9SIL1dUv84ucYFlCIQMQ63UWGm4B5JFCEbhzlA9z4+/f+iO6/B6UuqWD0/ruRrjduXDk\nVQJhqRl7ZTV624ZG3B6JsK5+E9bVb1Jd9/Tu/+CPS26SLNvetgOzRyt7oh25O6H9i2irbJaWcn+v\n3ihLDP87bH6fU+EMtbthKWxW3SZRlMcR/o78/CyQJKGoxxbvuX/z4zn4zY+54Ae2JnoG/IlmD2xW\n0pAd+PvNi/usFjUQv9Gvg7SXXgagXmPd8PAyv84+mrS3S423252L5uboro0P1yi1y3NctkilLBvB\n4s9XzsU1D60CAMkxlww5PfL9WyPPxJPt/4bH45MY/fa2XpC96u0YmVOG5uYuSVuHZJVGwj8BLupD\nzW/Lb5+NArRHi9slQgArf5QZgGTgnLkKhCWE3k3nqO769s6VOG/02ZHvm1uUiVZGYf0OsAFuVFTf\n2oLm/C4wKjo9WhC2ANgQCZZhI3+/mrYLAByurQPrzQIcPQieGA9bOVfu0sg9oWh30AqmNweWXI2H\nONwtaWnxwEHkwIEc1fOoLXPNWmW4Hd6qBXBOia9wyqf7pAl/43PGo6mpUzLSZVk2Zhkhbj/hb4s2\nYUzmtoHpKoqc2x9Rf1WeuLGxM6qrIxpyO0DmtMM6/GDke1t7t2omWDz3iZyTJ42rYsrtgBqtrZ6U\nSChrvWjifb2sBHAJAFAUNQtAHU3TXQBA0/RRAHkURVVSFGUFcH54e819UkGBTC/jR2eOw2M3nRb5\nToAASRIoyXdi0RThXVToKEC+Q+jN2cMuHl/IJ7nx5RmKYs4oP12xTCwYBkiLqvNMKY7NJ+2a+5lC\nKsAxdS1ccz4X4sWt6vkKHx39HM9WvQwA+PL413hl1zuKbfINToyyDBnRnVlx9FXsad2nehOHOoo1\nj0FYGEM3/jFesjpki6prHg3CGtQ2+GBhCYc4xq3MaBC2Jx9zijjpCP/h2GQCNjZIlSwf3/ksfrXq\ntoj6aUe3H/e9vBUxTyoTLMQGO9o1cEzkRiMtcmExFdVUtWJGieKYtDGqjtLwGAMqtGg1WHtC7YpR\nYbfyt+cLk7apMPh6xGX0aZpeB2ArRVHrwEXhXE9R1JUUFcl/vw7AawC+AbCCpun9avsk3nxtPpcp\n/vEvfX6GnM8efPAXC3D1+ZMik5vtPqkR4CfxgkxIkhqvlwik5lu+dtpPJN/nD5PGyd808+f42dTL\nNY+pha1yjxAZYfUpskwJm3aS2vamXbj+y1vx1sH/U6zLsWVjYvF4lb04GJ8TF5VcjWDrUPj2LpCs\ne2Ln8xL/PAD0bl+K00uNT2gX9E5SXe52lcBqIUCAwM0/jF1HxTCi59DoQ/mTSf+juY6ftwCAYIuy\npN54y0L0bjoHTGdyEoj+s5fTRHr1s/04VNsZc20AANJcLoLAorLoAmBWuZtCcV4GbV3xFezRQk0+\nnGW5hCkx3xlzblLOl0hNgNsum4XnbluKWePVs8r7grh9+jRNy2T3sFO07msAC2Tr1fbpc36wdCwu\nWTwmMrzkH+j97QdVt7eGo0mCbFCi2+I0mC3J47DY8fjSB/Hhkc8wOZzY9SPqe3iNfhsAMF4jmeu+\nRX9AiA1iQ/0WfHjkM2X73LWwumvRu30JHJM2KNYT9tgesPNGnYV5Q2cjz56LFfuVvX8e0uHFmdMo\nnDFlPK55aBUXEig2lPIdAg58e/oUrN30hnxNhGOdgkqlA+qTZftFiqJBjbDG5BC7kSzLVo9L8FYt\ngG34QVjCmcSBw1NhLamTbPPch9UACLD+5OivhNgQPIFudPGZ6TrFvNWijNTca2PyR2Ftnfr8giay\nkah12FF4epIrHxzUiMW0yqQPfDraWrEwqszYCNOv8XIgCEIij9LXDNqMXAUi/56aP5E3vnKsBPde\nDMl6+nrRDG6XuhuDIAicP/psjMrnRhsua/TU7HxHLoqchRhboCy7JsY1c3VEH0d60th6JbNLp6PE\nVQS7xWbIrUGSBJ69bSkq/Isky/kSk2IKc6R/rzgEVo7Xq/4guxxCKGUqlSiNxrRPLBJGQ3y4qfJY\nDApyw6GenjwABLx75oMN2uDdsTiyXfmQ5KoqNveIImh0JCWUES4syJwOxTUwEjHGRAmGtw47gjdW\nxx/VpEZQxbiG2BCCrHSUW9NVq9guHopytbWvjNJ8MrmjnVjIGKMfb2JGpKfPhCRp5npDfiMZtYC6\nPr8WFiK+QZlRoTMecUz67FIupPM7Y76tuw9JEDhwXDq5pSUj8M+lD0Q+j8hVujl4bCHO/RZskCas\n2CKZlkSfFBKPxiSx0dcIGWW82SLBOO6+YbsL4N22DKxfeBFGU7DgS3OeWb5Yf8MwNtIWue/5Uo6G\n0HAFleXoR1iTOe0IRSu+kkDughZqPf12fyt6c6TBHHYyeoKfERKdhAaA4SWpKZRjhIwx+kaZLNPU\nIcO++3X1m8Ak2Z3A6/OcN+qsqNtqxfhG3S+qvoz8PMItMbF4PP52+t04u2Ipls+9WbMkHIf2w+4/\nNgETygsix78oHNM+3T0Fp5UuUd1n6qhS9G46B4HjUt8+y7KRQVt/MPriv1o+jxEhaEegm3MH6tX1\n5WUMvFUL4TswQ7H+lGGz8fjSByPXLxr/3vta5JhklvFkOWupegUsa5SOh23Mrqg9/ZSg9qJhAU+h\nNCJtks4cVSzE+yyK4eWa0+HmGbQqm3KMDp3lxSLEdWyZ8CM+sWh83AWWxeQ78vDPpQ8YOlYyzmcE\n+QiGV3ocqdMrB7iKSOK+JNlbCMbVjlBnIUKNlbjlt4IRO7N8MZaOOBUW0gKnU6k3DwCnTByCTzYq\n9e5ZlsuOJkg2rVov/0N9D6/Tb0tq97qztKOTAjXjwXizEGrRFsqrbeEm4dmePO5f0BbpGd+36PcA\nYov0aOhuhK+2E8mqTDuzdCo2NGzR3oBFRMd+RGk2dhxsUR01TBmdnMlqHrWclaBKh2Bkrva1TwcP\nX78IjjRoBA36nv4jN56KP105F5Mqjd1ovKYKT5ZNNPwO30hXTv5R0tpn1JgPi6P4Qzxo9lajsHCK\ndOhv8XEl5QLHuJ66XE+c93+Le+tiH39h2G+am2WTRLsEw0lCLEsktdDIJeMuxD+XPoBHFt9raPvT\nhs/H40sfRLFLWfdWzl+vWwgwVoQaKzWT5dQI1AqT++IwYgC4bMIlku8so/4yYFgWhDM55RdHi6Um\nVE9mwTMfcAqdFUPyYCvfqxpGOZdSr9kQLyFWmVy3bZ9SfruvOk5GKcx1IKuPpBfE9K+rkAJcdmvU\n4g5dIslkp0U62ZhrE0YIO5urAMSh5pgE7Jbk+COjYbQjKfbNc/vJy0uGe3h6USMAipyC0TxX5ObK\ny7bj91fMwT3XnIJQu2AkQgwL0tkL0uGNaU5ETL5KgZulI08FSZCwWYw/hGq9bn4eRExxvhPZzjgG\n1SHlPsEQg6ojrRjtmIzLJ/5AWN5Qqd3OJPnRo40yyCwPGtv5hEIW1qHqlcmONnYltZiK2rFqW/qm\nIAnPr384HfdfOx93Xj4bV5yT+noBiTBojf5D1y3A8stmqVaskdPYI/QK5gyRPrRqN3q8veF0YyOt\n+OMpvwUAjCsYrTqXYNRPLu81yfcTIj/0DYXTKrjT5BO7o8vykJdlB9MhaNQHQkKvjo9oWjJCGjkU\njXsW3mF42yJH9J68mKumXIZfTr9asXx0Wb7K1lyWuJzll3Gy3aHWYQi2DINv7yk4XNeJQ3Ud+HD9\nMfx9xU7c8cxGzBsi1EsOntDxV8cTow+gRFYKNBbvs57nbdW2Wry3Jv7i9nIYFfdVqip0aTFlVDGG\nFGVhzPD8iE5Pf2XQ+vRL8l0oyY8eEgkAWaLQySkypUa1IaFWaF5f4bQ4dOv5anH3wtuRZ8/Fn+ff\nhgInVyRmhnsqXqPfxuGOowDizz6VP3b+YHjSO0pP38iQ+xcXTMezWw/BVnYEHzYIcf5Dsks5bXtr\nbCF0sfx+bb7Ys0cnFyt7elrzD7wsiJhRw8KuHNaCwGGuE/KXF5W+dGnREr3rLD33tVOvwNO7X9TZ\nnkMehaZWOS5ettLN+N7p8RUakqN2bZMddBEryy+blbSiLclmYHZZk8z7GqJmWogravUlt8y6DjfP\n/Dn+evpdcRV4zgu7NdxZxZG/oSxnKK6afGlkG7nvWI7/yCT4DyuLQQy3S/MIgkR4eB3FTvCGRc0t\nEml3luDa8jHcy45kufbHavB5RuboT0wnypwh0ugbd6F2B+TGi6dJvhsZnQLA3qNt+F7+r9C7RT/6\nSxwWCuirt4qRl+20kon7n0NdXBQXw3AT8R4DdSmio9LT76PIrmvOn4gHfqHIQ8X4kQWGagang4wx\n+o/veBaPbX9asTwQCmB3y17Dx0m1DoseYwtGYVzhGJAEiTNHKvV99KjI1c4dKHQW4M9Lb8F9i/6g\nuQ3P9YvPx03LlOnsNlkMtFBBSv96Dcly495Fd+pKGARVxNscjLq7xCg1njrNdctE19Ya5wteHhN+\nyWLtXu2McSWKZQ9fH91l9cgbu/DKZwcBhutRhjqUwQqEvQeOqd9Evv9mNqd+8sCpf4x6fDk20oq7\nFyzHwmHRi95pvet5raMQw+KNA+/jtjV3STKx40HtXIGgtKevNo+TCL/9nxm44XtTsXDKMJQaqH/b\nn8gYo1/dth+0itRCrD2CHJt+UkVplvIBTgV6ujhqXDBGXW2TZ1LpeOQ7oj8YM8aWYMooZWiiVn1U\nliVwx+WzVdfxFDjydY/R6qkAACAASURBVF0uBElA/vJI5bxKfY+gPW+Ps3fLzzecPYZ7gbgc+i8P\n3rfvtHPXoTCOrE+mQ6nnYquslmTWDs3itsm1x5f9W+wqMpRJHi2clmFZfHWCK7By4KRSETcWWJXc\ngF5/UKJ3lGwmVRZhZhr1cxIhI4y+3g0Ya4KPXmJGtjULN8+8LqbjxUuBQ7unm2vLURirVMe0ax+e\nQGFOYg/fxArlZCqhIyuQKDNFsfffH69VClGfuUNn4jezr8eVM78fWXbXVdqCZX//1SKcM28kfn+F\nsmC9Udgg95tL6v7KtW9iGLlo3WNnVSzRa4WhY4uTuBIdPWs9wYxX0DEi+lm4ZjoZ1FciwAQRYkK6\nE1Bioz+hMLoQlN7Nc/7ocwz1llPJspGn408LfoeFMkXEVGevaj62LIGivMSMPkkQGD1c2jNlNOLS\nk0G+yNjFWoKRhyRIjM6vgNUiGNmRpTm46ZJpuGBhJS46jRsJXLioEgCnTvnDM8ZFMjWB2KJlclw2\nhFrLEKgZBz8tuF/I3DbJdrEYfXlQA494lDA0qxS/m/MrYSXBorPbrxu9A0iNfocvtqxxBQbubXsM\nobiDnUFt9G9efQdu/ebPUuEpGV2iWq5qoV9y1MrR3TDjZxhfOBbzhs6Kr6FJ5HvjzofL6lJG02ho\n4SSL2ZT6UPfcUyqSohfezkp98J6e5BRyBwCnRaqYKs3DSO7LZfrYElx0+mhcsGgUnl9+Br572mjN\nbb+/1Hgd5oohOQBLIlg/BqxPcL/IRdNiSVDSc23xE/CFzgJU5on0kQgujyBapI/YI/NFzdeG26SG\n9iBTWHPhaH39qExiUBt9APCGfDgqmiiSuzlW0IJ8sDcYX0LHhKJxuGnmtXFHksTLGSNPk3y/cca1\nkc/yv3NH0+6UtsWmUR4xWeHSHmi/uBNlyUjppGkyQxMTYems4Zg5rsSQ9vr3xBPFUcJkjRJgtF+s\nZ5UvAQAsHrFQuoJkuF5+lEsYSqJGj97vRbAk/rn0AcwsVUacZSqD3ugDwIvVKyKf19RtlKwTJ2Yd\nT5L0al9x8bgL8MiS+yLfeclmDumD0FcZvXIWTNJXZoybOBOO1FAkqYkOnc7UfYfNghsunoZffW+q\naligmOI8JxZNje1aR1PrrGqp1lw3s3QqHltyP6aWyIrdEAw+21IDuqZNfUcAsPmSK8ymcihLYWM4\nK5zod/IL6SbjrsbrMt38oE5vhieemPi+QpwzIK7mJXc1/SDOCclEyXMlpyiInFjVQ+XwPdVThs5W\nGAWrKJLIkaaXpZzSAhee+u0SjB+Rj6vPm4ihRdLrarWQqBzK51gkp6d/29wbdderRVwRBIvjjR6s\nqdUutkLmtGtGe8WDmlvWUlwPgE1riHV/Ja4gZIqibABeAFABIATgpzRNH5Zt80MAvwE3uf4FTdN3\nUhR1JYB7APAi75/RNG1M4SqJVLVURyapekQuHbmeDM8tM3+B29bc1SdtSwTxQzgqvwIkQUYmcNPV\n0yeSJB1b6B+Pdvt+4bgJ6smcN+osjMwtw1SVQi5jC7T97OnEZiWx/Mdc+GtdSzc+FqmQ2qyEaj6D\nHtHCUWMJ6xyXS+FAFw3HtG/g3b4UpGbtYYB09ILpEEXvtBurPxELBMmAJfq+/uxAIN6e/qUATtI0\nfSqAewHcL15JUVQWgAcBLANXNvFMiqL4ceAKmqaXhP+lzOCrTbjyPLnr33h135uK5VrDwBx7+goe\nGOGxJffj4dPvVizv68xh/yGl3zRZj1zdjlHw7pmfpKMBNosNs4fMUI3qGAiG4tRp0iLfFguJOVQp\nLCSBq841FnGUl8RIM1v4BUKQDCwFTbo+fVs5Lenpq4VB13TV4Y41f8GRjmNRz60qvEcw4BrR/3/L\nviZeo78MAD8D+jkAyUwYTdM9AKbSNN1F0zQLoBWAtth4CmjsVkqripHX+vxukoompwMLaVGt2Xuu\ngeIsySTUqhSaSp4BJcB2FyTpWNH54fiL8POpP4m+YZoYVpyN525bGvlOEgSK8514+ndLFC8Engsr\nz5ctSZ5B3NtRJT1slEOLvTtqeXnvH/4YHf5OvLH//ajnVvMUEdYgQJjuHTXiNfpDATQDAE3TDACW\noiiJ/4Cm6S4AoChqKoBKAHzV7sUURX1CUdQXFEXNjPP8UTFS3UYsqVyRNyJVTUkbhToJXH0By5CG\npZqjMW4E97cw3X2TB3H6iAWY5tau4dsfIAgCI9zZimVajMiXzk1FK4wTPyzIrC7ttUHpCJRQqeMc\niT4zcP90BrTnd0yjryTq+J+iqGsAXCNbfIrsu+qVpShqHIBXAVxK03SAoqgNAJppmv6QoqgFAF4E\noBtLVViYBassHNDtjv7gP1e9Nuo23zQJ2xQWZBs6rtHzx7NtsinwCgbBSDuS3VY/PRsl381VlQ+O\nlZ9eOAV3PLEWbFhnJis4NK3XNhZS2c7/vXUZGJaF1RK9/5af75S0xe2eDOgUwkqk3WxAOofksjph\nIS3w+LsRapNGGfV4GcW57GE5CrvVotkOfvnbX65QXU9meRBAep9Bnv7QBp6oRp+m6WcBPCteRlHU\nC+B6+zvDk7oETdN+2TYjALwL4HKapneEj7UPwL7w5/UURbkpirLQNK2pg9reLi247XbnorlZuxfB\ns62+Kuo2Xd3CsdtP9qCZiH5cAIbODxhva6rwdAk/SbR2pKKtjKcArS0e9MZTQERGVyc34U6Gq0D5\nWE9ar61R0n0PiDl5sgfNVuNtibvdBAuEpB21P8+/Dbtbq/Fy9X/BeGRuupBFcS6fj4uqCwYZ1XbE\ncl3Tff3TdQ9ovWjide+sBMCLilwAYJXKNs8BuI6m6W38AoqibqUo6kfhz1PA9frTJnz9de36yOd0\n1ltNFemo8CWBtSTNveMOKxkSNi5qJ2QzXug7Ewm1pS/M2D5qj2LsTxIESH6hPMdCJeeCT7gaCJPq\nA414u2ArAJxFUdQaAD4AVwIARVHLAXwFbuL2NAB3U1SkoMTfwbl6XqIo6hfhcytLDKWJ/pKFmUzG\nFIyCy+rCuZXL0t2UhClIULQt49AxpH0B482CuK9PEKSmf91SUouarjrVOYbDHcfQG+w1pOxpYoy4\njH64d/5TleXiQHetrJylGstTzoTCcdjXfkB1nZ5qJQB8u3IZPj76RSqalTJcVif+dnp68wuMTKgb\npTjPge6AHYTNj7xQ/y5Jl24shU1pPb/chUMSpOhekL58CJLFA5sfwWUTLlEIBQLA07tfwk0zr1Us\nN4Jl8BYHjJuMyshdVq5deGRodqnmOiC96fgDGYsleUY/J8semci1ssoQVRN9orkwZ2WdhbsX3A5A\nKUIXjckF+to2JAg4wtpUlpI6gFB6dV8R5c6I27pfpQ6GUVw+s3MgJ6MsmUQNMEashNljiIZcGgCA\noagSo1SKys8NwimYpFLEjFIsi3bJhlgrUOwqxF0LbsNfFt0e0/lmF+lX0yIIAtnhAkSW3JOwle/T\n3V5cyCYW/nDKbyTfvd6+LZA+EBi0Rr8sWxoW9q3KZciyxe8XTJeMwUDiV99LrZLhsOJsozU6Mp5f\nL7xcsSyapjyvYVPiKo7Zh24hpaaEkM0pECAkgQVkno4gG7gyprHiPzRNoZOVn23OBcgZtN3XO+bd\ngiAbQkN3I4ZlD4m71inPgrK52N9+EMuiKBNmMuICIKngjFnD8c6nKT3FoKEwKwcXj70Abx38ILKM\nKtTX57eS8fcBFXM3cqNPEBJBQNKlXxA9yErdP9ubdkeVRw61cpnI2bYsdAe4cOzJFX0qBDAgGLQ9\nfYIgYCOtGJk7PGGDD3Bqi9dO+wnGFFQm3rgM4ZEbTk3q8awWEhZSfTLQRMmIXEGOgTAgMVxZGn/N\nV3lFOdvoXdL1IHTrIIvp9vkQYqRG/9mql4y0AgBw4ehvCe2Is8bxYGbQ9vRjoa+KmWcaedkpdImZ\n4dtRGVcwBjfP/Dly7DnIt+dF3f7U6SPQ1hpf/oPY5LNBK6d9E+ZnUy4HQRjTtb/783+jkaxGZW4l\njnYdjakNWeHi89Pck/FaWEI9GR2+wcag7enHwp3zfp3uJgwamB5OjpcNpKaHld/FDfHLEF/t2kyC\nIAiMKxyDYdlDNOezcsKTq+WOsaJRVBzn0ti3xFmMGWG3jNi9o0UjyRVu6eiN/eVz7YWckC8pMmtW\ng6OLTMJ8DcLsDSQTX/U8WAqaEWotA85J/vFvPutcfLF9Gs4/rTL5B89AeEM8oqgwweMIhlbcyxf3\n7mPJEG8PapfH9AS6UcIqtf6L8rgwU3ESWm/Qa/icmYLZ0zdJLiE7Qq3DMX5EamSQ3QUu3PCDWUkR\ncTMRZA4SnSEpzFWP65dE7Bjo6UdjZ3MVbvvmLqw6sk6xbng4kGAwZtcnE9Pom6SECRWJ9RxN+gZe\nGiFRQ6nlr5cuT9wYr6nlalyvPPg1AGAFzZX1YL3ZkReYOLHMHMUryXij/9fT+n8ZRBOTlJOgPdYS\nRhOHcib6YjnedQJ722gA3AiCZdmIaCLBCqZMnI9gZB4h08hYoz/TPRV3L1ieUMKWiZLf/s8MjCzN\nwdJZg68ozWAkyHL+9wCTWM1hLTE1cShnnj12TfkipzBifHDzY5HP4vrPAACXunRxkAmqLs9kMtbo\nLxp+CopdReluxqBjUmUR7rpqHvJTGa5pkjT46nFbm3YmdBwtYb1E5b1vnHSz6vIQyygSuNTwBPST\nwDKRjDP6FXkjAXChZCYmJsnBSE8/VhhPPv79UbXquoNtRw0Jsc1wT4n7/IOVjJvluHnmz9HmbYc7\nyzT6JibJwhvyqS5PRFqbZUn4Atq9+X/teiHqMcyJXCUZ19O3W+wKUSYTE5PE0HLjiMM0Yy5SzhIA\nCEwqoqJuqoXLakpwyzFfgyYmJgmjZdALHIL8Q8ylD1kCXT1+dLbXR900UCeVkr5tzo3Y20ajPNcM\nKJBjGn0TExPMGTIjsQNoGPSLx14Q/zEZEi0dXrjYzqibBk9IRwPleSNQnmcafDXiMvoURdkAvACg\nAkAIwE9pmj4s2yYAYK1o0TJw7iTd/UxMTPoOG2lFgAlGLRcaDa0+fI5dX26b6c0C6epRX2kxwy1T\nQbw9/UsBnKRp+jKKos4GcD+AH8q26aBpeol4AUVRPzawn4mJSR/BlyWM2d+uIL79g3VjEWovBeHq\nhnPyesk6S1674eOUlyq1eEzUiXcidxmAd8KfPwewKMX7mZiYpIBxhWMARK8RHY24a0iTIYCx6mYE\nB1vKoh7GZsu4mJS4ifdKDQXQDAA0TTMAWIqi5Nk4ToqiXqUoai1FUb+OYT8TE5M+4qrJl+HqKT/G\nvKGzEjpOsbMQi8pOwUK3fj/OJguhtJTUhj9pjxTYXv1efKBmPA7VRvf7m3BEde9QFHUNgGtki0+R\nfVf7xX4L4GVw7/CvKYr6WmWbqGPCwsIsWK1S/Qy3O/Z07nRhtjU1mG1NFrmoKBN6+Ym09abSK7H+\nUDXWNQtTefLjWUgLAmJphICD+z+rZwqEYUCOJR+eUIdkbbB+tOq5+hP9qW1RjT5N088CeFa8jKKo\nF8D12neGJ3UJmqb9sv3+Jdr+CwBTAdRF209Oe7t0ksftzkVzs7rORn/DbGtqMNuaGpLR1u5OaZKW\n/HgMK/XjBBsquA8EAzmsP/xCENXbza9fgtbWOjgmblZs31+vc7ruAa0XTbwTuSsBfB/ApwAuALBK\nvJKiKArAnwBcBsACznf/JgCf3n4mJiYDHcFjHOqKXlOB8XD6VwSpNPq+fXO5I+a3RpYdrPEAoWJ4\nd50Ka9khBI6ZFdRiJV6jvwLAWRRFrQFnyK8EAIqilgP4iqbp9RRF1QDYBIAB8D5N05soitqqtp+J\nickgQdyR13XZcHzrlHIQAD7d06Y8lJfz5VtyRVE8DBlZFzg8XTjOvPK4mpuJxGX0aZoOAfipyvIH\nRJ9vM7qfiYnJ4GBocQ5wgP+mNPryJRedNho2Kwk2uwVfeTZFlofahXkGpicHZFa4Zi6rro8/rCQr\ngVZnFmZGromJSdKwEiKTYqBmis3K9dyHFLmAsF3v3XwWICqKwvqygCwPmF7tRK/CHEdc7c1EzOBW\nExOTpFHsEoqekNnKMEqxeFpkohZAoUNUXpO1QDWwj9E2V5MqzdoYRjGNvomJSUogrEoZhYvHCVo8\nowoEbZwCeyECxyn4D0+VbL9g8hAIQwbtOQKSTDSjOHMwjb6JiUmfUegsQHnucABAlktwBREEgWDD\nKIRahku2/9kFk+NVeDDRwDT6JiYmfQrBmx2Rz19fdZlVbG8SP6bRNzEx6VvCBp4VWXE9m890cf76\nUKdZ7S4ZmNE7JiYmfcrEovE41lmDUXmi2Hqdrn6wvhKMJx+Mp1BzGxPjmD19ExOTPuWs8iX48cQf\n4MyKJZFl4nnYB36xAHYbiWsvnMSvBdNVLAnjFGMmZsWG2dM3MTFJDRo+eKfVgQXD5kiW9fqEAuil\nBS786zdLdA/99oMX4GR7d6ItzEhMox8njz32MPbsqQJBELjppt9g4sTJ6W6SiUn/woAMA4/PH1uV\nrFjL7ZoImO6dONi+fStOnKjBU0/9G8uX/wGPPPK3dDfJxKQfYtwyl5VoZ9vyWbtiSNPqx43Z04+D\nrVs347TTlgAAKitHoaurE93dHmRnmyXbTEwixNDTz9eRUWAYwU/0wM/nIy/bbiZjJcCAN/r//fIg\nNu9rSuox504oxQ/OGKu5vrW1FRQ1IfK9oKAQra2tptE3MZEQu2FW68GLNfjzsu1w2ge82Uor5tVL\nAixrZo2YmCiIoacPAP/41SLYbUoVzdLCLDS2ccWU7FZ1lU0T4wx4o/+DM8bq9spTQUlJCVpbhcIO\nLS0tKCkp6dM2mJj0e2I0+louniGFLjS29cBCEqZbJwmYE7lxMG/efKxe/QUAgKb3oaSkBFlZ2hNR\nJiaZSXIMND+QHllquk+TwYDv6aeDqVOng6Im4he/uAoEQeDXv1bUizExyViGMBPQSO5L2vHKh+Rg\n9+FWjC7LS9oxM5m4jH64qPkLACoAhAD8lKbpw6L1swE8LNplEoDvAjgbXN3c2vDyl2iafi6eNqSb\n6667Id1NMDHpl7BJVka7cFElRpbmYPpY04WaDOLt6V8K4CRN05dRFHU2gPsB/JBfSdP0VgBLAICi\nqAIA7wHYAM7oP0rT9OOJNNrExCRzsFktmDdxSLqbMWiI16e/DMA74c+fA1iks+1vATxC07Sy3P3/\nt3fvUVKUZx7Hv3OJAgoM4LCIyREF/UkiBtQYBXFRkU08IAkR4rJ420Vco0aUlYsILkQTcA/ZgBcW\nQ7gY0d242SAaMOAJEhWjEV0vrD4xisZgomMiKBeBgdk/3newIQMzNHRVNf18zuEwPd09/evqqqfe\neuvtep1zB599PIHrkpVv0e8A1ADEYl4n6ZDdHySpOfB3hJZ+vcGSlkl6RNIxeb6+cy6jdl4vv9zb\neVnUaPeOpOHA8N1+/eXdbu9p1/414Oc5rfzFwC/N7FeSLgLuAPrv7fXbtGlB5W5jc6urWzYWOzM8\na2F41sI4EFlzJzcs5HsvteV6oDRa9M1sNjA793eS5hFa+y/Gk7plZra1gaf3B2bm/K1nc+5bBExt\n7PU//HDTLrerq1tSU/NxY0/LBM9aGJ61MA5U1q3baqECoKxg770Ul2s+r9uQfLt3lgKD488DgOV7\neNyXgBfrb0iaLql3vNkHeCXP13fOZdSBHr3jDqx8R+/8F3CepCeBLcBlAJLGAivM7On4uCozy93F\nzQZmSdoG7ACuyPP1nXMZV+YzmmdSXkXfzLYDlzfw+ym73W6/2+2XgZ75vKZzrjh0+WxrVn0AzQ/1\n735mkX8qeVi8+GGeeWYlGzdupKbmfYYMGUr79u2ZNetuysvL6du3H0OGDGXRokXMnTufiopyOnXq\nzJgx49OO7lzBNTsk9BpXlPtVXrKo6Iv+//zuEV54/+UD+jd7tO/GoC57HVTEmjVvMmfOAjZs2MBl\nl/09zZo1Y+bMObRq1Ypx40YxcOAgNm/ezLRpd9CyZUuuvvoK3njjd3TunOzF4ZxL2qZtmwFoXtk8\n5SSuIUVf9NPSvfvJVFZWUlVVxeGHH05dXR1t2rQB4PbbfwBA69atGTduFABvv72G9evXpZbXuaTU\n1oX5bivL/TLIWVT0RX9Ql/6NtsoLIXc2n7q6OsrKdj2U3bZtG5MnT2bOnAW0a3cEo0ePTDqic6mo\nn1+i3E/kZpJ3uuVp9eqX2L59O+vWrWPz5s3U1tZSU/M+dXV1jB49kk2bNlJRUUG7dkfw3nt/4rXX\nXqW2dt8mf3auGNUP2SzzeWwzqehb+mnp0KEjEyaMZe3adxgx4ltUV7fn5pvDJZbPOacvrVtX0atX\nL4YPv4QuXY5j6NCLmTHj+8ybdz+Vlb7Y3cGrvqVf5m3KTPLqk6ejjvos11yza5fNrFlzd7k9ZcqU\nXb6Jd9FFwxLJ5lyavKWfbb4rds4VREOTnLv0eUs/D+efPyDtCM5l1pDjB3Lfqw8y9IRvpB3FNcCL\nvnPugGrfopobTvlW2jHcHnj3jnPOlRAv+s45V0K86DvnXAnxou+ccyXEi75zzpUQL/rOOVdCvOg7\n51wJKau/ToZzzrmDn7f0nXOuhHjRd865EuJF3znnSogXfeecKyFe9J1zroR40XfOuRKS6aIvyWdh\nOIAkdYz/Z/pzLzaSWqWdYV/4dnXgFdO2lblx+pKqgG8CDwEfmtkWSWVmlqmgMedI4BngN2b2QRZz\nAkhqDYwGLgVON7M/pBxpj+JyrTaz19PO0piYdQywGZhiZltTjrRHktoA/wQ8DKw1sw0ZXl992yqg\nTO2VJA0Afg6cBPwzYYMiax+2pAuAhUALoC8wFbKXE0DSCGBRvPlDYHtWW3qSKoClwDhJR6edZ28k\nXQUsA9YDUzNe8M8lNKL+BhgMTIPMrq8D8G2roDJR9CV9Jv54DHCvmV0NfB84T1Lf+JjUs+ZkOBqY\nb2ajCSvoazmPycyHLulEoCMwzMzGA6cBh2ZtA8pZrscCW4Ba4GRJh6SXas8kHQGcAawwsynxaLQq\n5/7U11XYuRMFOIqQ9UYzuxXoLWlQfEwmsub4HMWxbZ1C2IlmettqSKrdO5I+D1wF1AB3AjcBvwXm\nmFmtpPnA0WbWJ7WQ7Cyew4BXgXsJraUngU+Ax4FVwHNmdldaGevFrBcDrwALzGxHzn03A6vMbEla\n+XLFrP8I/B6YBRwKtAV6A6cCd5vZq+kl/FRO1reAe4CzgTOB1cC5hKlH15nZdWllrLfbcv2P+PNh\nwH+a2TuSZgEnm9mXUowJgKTOwBlmdl+8PQxYDmwCVpCtbasz0NPMftzAfZnatvYm8b18/d5aUj/g\nLuBpoDNwM3A/8FVgjKTbgP8FPoorQlo5TwLuJmzcXwbGA782s3eBKuAHwPeAwZJuic9JdLk2kPUV\nQstjkqRO8b5KoA2wIY2MDWQ9nvD5vwh0IxzGdzKzN4EHgWbAmbEvOrfVmnbW7sB3gG3AR4S+518B\n1wPdM7AO5Gb9ImG72gAcCUyQ9D3gXeBjSdflPjfprNFU4NuSzou3F5jZWjK2beVkvVbSOfVZJJXF\ndTP1baup0ghX35XTCVhtZvcDNwJVZvY8MIGw8LYDM4GfEFaApNVPGn8i8Me4dx8bbw+S1NLM3jSz\nOWb2W8I5iAslNcttXaeYdTxwHDBQUpWZ1RJaqSMBUshYr77L5vNAjZnNJxTMj4FzJXUws03AI8Dp\nwBHx8Wkc2jeUdSSwFehKaJRMJLSg/wJcCXwzpXWgoazXEZZbS+AB4jkIM5sE/CtwtKTyFLokPgM7\nd1DbCEfPw3JP1prZmoxsWw1lvTRm3QGUm9l2YA3pb1tNkljRl3S2pAXALZJOJZyZ7y3pBuB5oFrS\nXcBrZjYdmGRmnxBagetSyDlJUveYs4Wk7mb2EbCWUFx7KOgan3os8FjMnKWs3eI/gHlAmaTT4/MT\nK6Qx60+BaZJ6xqwVkk6IWZcRzj+cAmBmDwHvAGMlrQQuylDWxwh9z5Vm9ihQfxTSCViSwjrQWNbO\nQCsze5DQkIJw1PpGkgVqt6xnxoI+kXAydAtweXxcpaSuGdi2Gs3Kp42R+aS0be2rRIq+pGOBW4Ef\nEwr8hYQN5HzCSJ1RZjYAaAeMkNQDuFXSE4TD6dUp5HwOGAKcRRj5cI+kuUA1YYheh/gexkpaQmit\nJtaftw9ZNxE2egitwZeBnpDciAhJRwK3AbOBp4BLCBvNEuCCmOVxQndJ5/icZoSW/heB2+r7fDOS\ndTmwETgm7mwnSXqUsA48mkTOfcy6HugqqQMwStJThG7U51PKupLQWr48Ds39AFgMnC+pYzwq7Raz\n1i/XJLetfc0KofivJuFtKx+VjT8kP7GfqyfwG8KonOfM7FGFERkdgCsIJ0FbExYkhFbIZMKJsteB\nc2OLr2AayXkk8A1CQX0CON7MHpZ0NjDGzL4Sd0xnxRZfQeWZ9W8Jh/LzzOzPkmab2TsJZZ1A+Gzr\ngJ+Z2ZJYzP8ATAHmAm0l9YlFfyUwCZgR3+dPzGxORrM+DYw3szsl/Z6wDizMaNaVwEQzuxOYKOkJ\nM1uWctZ3gdskPRzH4b9A2MkPBqYTzkksBnqZ2S8ymvVCYIakzmb2hqRZpT5OfzrwXcLhugF9JfWI\n45nLCIdINxAOlwbH57QhjN45zMw2FLrgNzHnVuDbZmaE1h3A8cDjsV9vUxIFfz+yngDs3GgSKvgd\nCediWhOOiu4g9Nm2iIfovyac/OxF2IHdojBstx3wVNwAlydU8PPN2hZ4RlJzM/tLQgV/f5br05Ka\nAyRU8BvL+mTMeGPM9HZ8/AhJq4FTYw1IouDnm/XKmPW0+PvMF3woUNGXdBggQmuoT1wYMwlnvlcQ\numx+SBih8RJQJ+kRwomw6Wa2vhC58sx5D6GbBOArkh4iHEYvSvIQbn+zJpUzqgbamtn1ZvYjYEHM\nPi3evwP4b8L6CRL8NwAAA3JJREFU9xjhPf2I8PnPM7PtCS7b/c26OaGcB1tWgPuAIyV1jF0q/05o\n9F1iZguKJOulZvZAgln3W8HG6cd++QrCmPHlZrYwtjq6mtlLCmP0rzezK2L3xNGWwlfvm5CzKzDa\nzC6Pv6+2MFwzccWSNfYdf4Ew3rqcMGxwBWHkw9fMbJWk4wjfuL4yPq2lmSV2wt6zZibrvwBXE1rZ\nZ5nZzzxrYRX8y1mSriV8c+2nZvaCwlfC1wP9CCdCrzWzLQUN0QRNyHmNZeSr9kWWtRL4JTAQ+Aeg\nP/BvhBO0vQnfaNy457+QHM9aGE3IerGZbUgv4aeKKWu+Cnkit37M7VJgOKFv+YX4/0mE/vvUC/4+\n5Ey9iBZT1hzdAMzsQ+BOSRuBcwgn86/KSmGKPGthNJY1S0W0mLLmJZHLMEj6KqFL4hjCkLapluCY\n26YqlpxQPFkl9SeMaV9IuCTAs8B3kzwf0lSetTA8a7YUrKW/m68TWqK3m9m9Cb1mPoolJxRP1naE\nr9N/nXAy8f6U8+yNZy0Mz5ohBS/6cTjUs2SgK2dviiUnFFdWwqUfxhAunpalbqeGvIVnLYS38KyZ\nkblJVNzBRRmd/KIhnrUwPGu2eNF3zrkSkulLgDrnnDuwvOg751wJ8aLvnHMlxIu+c86VkKTG6TtX\nFBSmlzTCxcogzJz0BDDZwoxee3reMEvomv/O7Q9v6Tv312rMrI+Z9SFMet6SMH9zg+KloCcmlM25\n/eItfef2wsw+kTQSeF3SFwiT/LQl7AgeNLOpwBzCfLNLzayfpCHAtYQ5DmqA4Wb255TegnO78Ja+\nc40ws22EKSn7AwvN7GzCRCU3SWoF3EI4Ougn6XOESen7mtmZhNnhbkonuXN/zVv6zjVNa+BPQG9J\nVxFmKWtGaPXnOoMwdeUvJAEcCqxJMKdze+VF37lGSGpBmJnscUIR72VmdZI+aODhW4Bnzax/ghGd\nazLv3nFuL+IMZDOAZYSJa/4vFvwLgBaEncAOwigfCHOpnhZnY0LSYEkDk0/uXMP82jvO5dhtyGYF\nYWKapYR++ROAB4A/Ag8BJwI9gNOBVUAtcBZwATAK2BT/XWpm7yX5PpzbEy/6zjlXQrx7xznnSogX\nfeecKyFe9J1zroR40XfOuRLiRd8550qIF33nnCshXvSdc66EeNF3zrkS8v+EPZU82i966gAAAABJ\nRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fda2ab2f1d0>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "H_zBxRZvriGh",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "... which does look even more mean reverting than the 50-50 butterfly (look how the t-statistic of -5.075 is even more negative than the one we found before). When I ran the code the weight results were assymetric: -0.83 for the 2yr weight and -1.26 for the 10yr, which is not standard for discretionary traders. In practice, traders would prefer simple and liquid strategies instead of adding complication for the sake of small improvements.\n",
        "\n",
        "But a note of caution: while PCA seems to give nice mean reversion signals as result of the implicit variance minimization, it does not *assure* it will find them. Of course, we could marry PCA + the mean reversion test, but as it happens Econometricians derived the [Johansen Test](https://en.wikipedia.org/wiki/Johansen_test) as an all-in-one method that finds the mean reversion signals in multiple time series.\n",
        "\n",
        "### The 'right way': Johansen Test\n",
        "\n",
        "It is well beyond the purpose of this blog to explain how the Johansen Test works, please read [Johansen Test for Cointegatring Time Series Analysis](https://www.quantstart.com/articles/Johansen-Test-for-Cointegrating-Time-Series-Analysis-in-R) (or another academica literature). Furthermore, python does not include a common library (although you can use \"pip install johansen\" for one) which means the Google colaboratory does not support it. The below code might not work, so please take it with a pinch of salt, but when I ran it the test confirmed the presence of 1 linear relationship corresponding to weights [-0.8095,  2,         -1.2880] which were similar but not the same as the PCA weights found at the same time [-0.8323, 2, -1.2613] (remember, if you run this code you might get slightly different results as Quandl would have updated the time series).\n",
        "\n",
        "The big difference is that the Johansen test was designed to find out and test all the possible linear combinations that yield mean reverting processes, while PCA (as used above) tries to find signals with low variance which may turn out to be mean reverting. \n",
        "\n"
      ]
    },
    {
      "metadata": {
        "id": "SgvQgjSA9Ap_",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 306
        },
        "outputId": "2d338c3f-1eeb-4e82-a8c0-1e6a38d3b912"
      },
      "cell_type": "code",
      "source": [
        "!pip install johansen\n",
        "from johansen import Johansen\n",
        "#\n",
        "x = TimeSeries.dropna().as_matrix()\n",
        "x_centered = x - np.mean(x, axis=0)\n",
        "\n",
        "# Because this library is not standard, you need to read\n",
        "# the code itself to find out that\n",
        "# significance_level = 2 means the test will use 99% confidence level,\n",
        "# model = 2 means \"This case should be used if\n",
        "#            the input time series have linear trends but the cointegrating\n",
        "#            relations are not expected to have linear trends.\"\n",
        "johansen = Johansen(x_centered, model=2, significance_level=2)\n",
        "eigenvectors, r = johansen.johansen()\n",
        "\n",
        "# r will have the set of linear combinations that yield \n",
        "# mean reverting processes at confidence 99%\n",
        "for i in r:\n",
        "  vec = eigenvectors[:, i]\n",
        "  vec_min = np.min(np.abs(vec))\n",
        "  vec = 2 * vec / vec[1]\n",
        "\n",
        "  print (\"The\",i,\"th cointegrating relation: {}\".format(vec))\n",
        "\n",
        "  "
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting johansen\n",
            "  Downloading https://files.pythonhosted.org/packages/c4/e4/f5cc4a3657ed82a1068fc33cff2e465c24d291e6afe270e56b97f5dfed01/johansen-0.0.4.tar.gz\n",
            "Requirement already satisfied: scipy>=0.18.0 in /usr/local/lib/python2.7/dist-packages (from johansen) (0.19.1)\n",
            "Requirement already satisfied: statsmodels>=0.6.1 in /usr/local/lib/python2.7/dist-packages (from johansen) (0.8.0)\n",
            "Requirement already satisfied: numpy>=1.11.1 in /usr/local/lib/python2.7/dist-packages (from johansen) (1.14.3)\n",
            "Requirement already satisfied: pandas>=0.18.1 in /usr/local/lib/python2.7/dist-packages (from johansen) (0.22.0)\n",
            "Requirement already satisfied: patsy in /usr/local/lib/python2.7/dist-packages (from statsmodels>=0.6.1->johansen) (0.5.0)\n",
            "Requirement already satisfied: pytz>=2011k in /usr/local/lib/python2.7/dist-packages (from pandas>=0.18.1->johansen) (2018.4)\n",
            "Requirement already satisfied: python-dateutil in /usr/local/lib/python2.7/dist-packages (from pandas>=0.18.1->johansen) (2.5.3)\n",
            "Requirement already satisfied: six in /usr/local/lib/python2.7/dist-packages (from patsy->statsmodels>=0.6.1->johansen) (1.11.0)\n",
            "Building wheels for collected packages: johansen\n",
            "  Running setup.py bdist_wheel for johansen ... \u001b[?25l-\b \bdone\n",
            "\u001b[?25h  Stored in directory: /content/.cache/pip/wheels/29/4a/ea/9a38d9ff2246e7a710488e8eaa22da457614a6f639979a59b4\n",
            "Successfully built johansen\n",
            "Installing collected packages: johansen\n",
            "Successfully installed johansen-0.0.4\n",
            "('The', 0, 'th cointegrating relation: [-0.80934929  2.         -1.28806607]')\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "pMUZxayxMRMF",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Bringing it all together\n",
        "\n",
        "If you have read this so far, I recommend you read the following 'drain the lake to fish' approach to pair trading: [Efficient Pair Selection for Pair-Trading Strategies](https://www.maths.ox.ac.uk/system/files/attachments/593233.pdf) (but remember that you might get nasty surprises if you neglect the fundamental story).\n",
        "\n",
        "In the [next post](https://colab.research.google.com/drive/1QNtBrFLsn7sqayEKgJVXg0UJz1JNYgVI) I will flesh out a ful example on the Term Structure of US rates to highlight the differences between PCA and the Johansen test in a case with multiple time series."
      ]
    }
  ]
}