docs/examples/academy/35. Hypothesis Testing II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 8\n",
    "\n",
    "## Video 35: Hypothesis Testing II\n",
    "**Python for the Energy Industry**\n",
    "\n",
    "In this lesson we will look at the statistical tool of Granger causality. Working with the same time series data from the previous lesson, Granger causality allows us to test whether one time series is indeed predictive of the other with a certain time lag.\n",
    "\n",
    "*Remember: when you run this notebook, you will be using more recent data that was used originally. So, you may find different or even contradictory results! This is the nature of working with real data.*\n",
    "\n",
    "Start by getting all the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from datetime import datetime\n",
    "from dateutil.relativedelta import relativedelta\n",
    "import vortexasdk as v\n",
    "\n",
    "# The cargo unit for the time series (barrels)\n",
    "TS_UNIT = 'b'\n",
    "\n",
    "# The granularity of the time series\n",
    "TS_FREQ = 'day'\n",
    "\n",
    "# datetimes to access last 7 weeks of data\n",
    "now = datetime.utcnow()\n",
    "seven_weeks_ago = now - relativedelta(weeks=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "crude = [p.id for p in v.Products().search('crude').to_list() if p.name=='Crude']\n",
    "assert len(crude) == 1\n",
    "\n",
    "china = v.Geographies().search('China',exact_term_match=True)[0]['id']\n",
    "SEA = v.Geographies().search('Southeast Asia',exact_term_match=True)[0]['id']\n",
    "\n",
    "SEA_exports = v.CargoTimeSeries().search(\n",
    "    timeseries_frequency=TS_FREQ,\n",
    "    timeseries_unit=TS_UNIT,\n",
    "    filter_time_min=seven_weeks_ago,\n",
    "    filter_time_max=now,\n",
    "    filter_activity=\"loading_end\",\n",
    "    filter_origins=SEA,\n",
    "    filter_destinations=china,\n",
    ").to_df()\n",
    "\n",
    "SEA_exports = SEA_exports.rename(columns={'key':'date','value':'SEA_exp'})[['date','SEA_exp']]\n",
    "\n",
    "china_imports = v.CargoTimeSeries().search(\n",
    "    timeseries_frequency=TS_FREQ,\n",
    "    timeseries_unit=TS_UNIT,\n",
    "    filter_time_min=seven_weeks_ago,\n",
    "    filter_time_max=now,\n",
    "    filter_activity=\"unloading_start\",\n",
    "    filter_origins=SEA,\n",
    "    filter_destinations=china,\n",
    ").to_df()\n",
    "\n",
    "china_imports = china_imports.rename(columns={'key':'date','value':'china_imp'})[['date','china_imp']]\n",
    "\n",
    "combined_df = SEA_exports\n",
    "combined_df['china_imp'] = china_imports['china_imp']\n",
    "\n",
    "# dropna in case NaN values returned by search\n",
    "combined_df = combined_df.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input to the granger causality test is a numpy array of values, with 2 columns. The test will determine if the second column Granger causes the first. We create two arrays to test both directions of causality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test if china imports causes SEA exports\n",
    "x1 = combined_df[['SEA_exp','china_imp']].values\n",
    "\n",
    "# test if SEA exports causes china imports\n",
    "x2 = combined_df[['china_imp','SEA_exp']].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now do the Granger causality test on x1. We specify the maximum number of lags - this must be less than 1/3 the size of our data. So we go with 15 lags here. The output is a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.tsa.stattools import grangercausalitytests\n",
    "\n",
    "gc1 = grangercausalitytests(x1,maxlag=15,verbose=False)\n",
    "print(gc1.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({'ssr_ftest': (0.11191210577788067, 0.7394982261424798, 46.0, 1), 'ssr_chi2test': (0.11921072137209027, 0.7298921072736446, 1), 'lrtest': (0.11906594393622072, 0.7300497647785229, 1), 'params_ftest': (0.11191210577793151, 0.7394982261424223, 46.0, 1.0)}, [<statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x000002B329B7BF40>, <statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x000002B326FF0BE0>, array([[0., 1., 0.]])])\n"
     ]
    }
   ],
   "source": [
    "print(gc1[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the test, for each lag, we get the results of a few different types of signficance test. The second value in each of these results tuples is the p-value. So, we pick a test (ssr_ftest), and for each lag value we get the p-value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.7394982261424798, 0.34406244046901324, 0.4379876250464344, 0.5978693386431604, 0.4297746020366883, 0.5913047894842787, 0.6333656452613916, 0.39251634523943296, 0.3460044222271352, 0.42328369935264354, 0.2947044888072325, 0.4534537474070228, 0.45562416494890623, 0.6404053509496745, 0.6127531951186226]\n"
     ]
    }
   ],
   "source": [
    "test = 'ssr_ftest'\n",
    "lags = range(1,16)\n",
    "p1 = [gc1[lag][0][test][1] for lag in lags]\n",
    "print(p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, none of these p-values are less than 0.05, so there is no statistically significant evidence of China imports Granger causing SEA exports. We now repeat the test for the other direction of causality, and we plot the p-values for both tests:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABNdUlEQVR4nO2dd1hVR/rHP0MXpIqK0hEbihrF3ms0GjXFNBNb3MT0stkku5u42dRNNs1NzPpLNZuYrmnGJPauUewFVERARAVFEETqnd8fA4hKucDtzOd5fK6cM+fMC1y+d847bxFSSjQajUZj/zhZ2wCNRqPRmAYt6BqNRuMgaEHXaDQaB0ELukaj0TgIWtA1Go3GQXCx1sSBgYEyIiLCWtNrNBqNXbJjx44zUsqW1Z2zmqBHREQQHx9vrek1Go3GLhFCpNZ0TrtcNBqNxkHQgq7RaDQOghZ0jUajcRCs5kPXaDS1U1JSQnp6OoWFhdY2RWMFPDw8CAkJwdXV1ehrtKBrNDZKeno63t7eREREIISwtjkaCyKl5OzZs6SnpxMZGWn0ddrlotHYKIWFhbRo0UKLeRNECEGLFi3q/XSmBV2jsWG0mDddGvK7tz9BP30AVj4HF3OsbYlGo9HYFPYn6OdSYONbkJ1sbUs0GodmwIABFp0vJSWFL774wqJzWpKFCxeSkZFh1jnsT9B9Q9Vr7nHr2qHRODibN2+22FylpaWNFvQLFy5QXFxsQqtMR1lZmRb0avErF/ScNOvaodE4OM2bNwdg7dq1DB06lFtuuYUOHTrw9NNPs2jRIvr06UNsbCxHjx4FYMaMGcyZM4fBgwfToUMHli5dCqjN3ZkzZxIbG8s111zDmjVrALVinTJlCtdffz1jxozh6aefZsOGDfTo0YO33nqLAwcO0KdPH3r06EG3bt04cuRIrfYePnyYjh078uc//5mEhIQ6v7/ly5fTv39/evbsyZQpU8jPzyc1NZX27dtz5swZDAYDgwcPZvny5aSkpNCpUyemT59Ot27duPnmmykoKABg1apVXHPNNcTGxjJr1iyKiooAVd7k+eefZ9CgQXz55ZfEx8czdepUevTowcWLF3n66aeJiYmhW7duPPHEEw37JV2B/YUteviBuw/k6BW6punwz58PcDDjvEnvGdPWh39c38WosXv27CEhIYGAgACioqKYPXs227ZtY968ebzzzju8/fbbgHKbrFu3jqNHjzJ8+HCSkpKYP38+APv27SMxMZExY8Zw+PBhALZs2cLevXsJCAhg7dq1vP7665UfBA899BCPPPIIU6dOpbi4mLKyslptvOaaa9i7dy9ff/01s2fPRgjB3XffzS233IKXl9dlY8+cOcOLL77IypUr8fLy4tVXX+XNN99k7ty5PPXUU8yZM4e+ffsSExPDmDFjSElJ4dChQ3z00UcMHDiQWbNm8d577/Hggw8yY8YMVq1aRYcOHZg2bRr//e9/efTRRwEVS75x40YAPvzwQ15//XXi4uLIzs7m+++/JzExESEEOTk5Rv0e6sL+VuhCKLeLdrloNBajd+/etGnTBnd3d9q1a8eYMWMAiI2NJSUlpXLcLbfcgpOTE+3btycqKorExEQ2btzIXXfdBUCnTp0IDw+vFPTRo0cTEBBQ7Zz9+/fn5Zdf5tVXXyU1NZVmzZrVaae3tzezZ89m06ZNvP/++3zwwQe0adPmqnFbt27l4MGDDBw4kB49evDpp5+SmqpqXs2ePZu8vDwWLFjA66+/XnlNaGgoAwcOBODOO+9k48aNHDp0iMjISDp06ADA9OnTWb9+feU1t956a7V2+vj44OHhwezZs1myZAmenp51fm/GYH8rdFBuF+1y0TQhjF1Jmwt3d/fK/zs5OVV+7eTkRGlpaeW5K0PthBDU1oj+ypVzVe644w769u3LL7/8wrXXXsuHH37IiBEjKs9///33/POf/wTU6jcuLg6A1NRUFi5cyJdffkn37t157rnnrrq3lJLRo0fz5ZdfXnWuoKCA9PR0APLz8/H29m7Q91bb9+fi4sK2bdtYtWoVX331Fe+++y6rV6+u9V7GYH8rdFArdO1y0Whsjm+//RaDwcDRo0dJTk6mY8eODBkyhEWLFgHKz52WlkbHjh2vutbb25u8vLzKr5OTk4mKiuLhhx9m4sSJ7N2797LxN9xwA7t372b37t3ExcWRkpLCqFGjmDRpEn5+fmzatImvv/668mmiKv369WPTpk0kJSUBSsQrnhqeeuoppk6dyvPPP8+f/vSnymvS0tLYsmULAF9++SWDBg2iU6dOpKSkVN7ns88+Y+jQodX+bKp+f/n5+eTm5nLdddfx9ttvs3v3bqN+vnVhpyv0MCjKhcJc8PC1tjUajaacjh07MnToUE6fPs2CBQvw8PDg/vvvZ86cOcTGxuLi4sLChQsvW/FX0K1bN1xcXOjevTszZsygsLCQzz//HFdXV4KCgpg7d26tczs7O/Pyyy/Tp0+fOu1s2bIlCxcu5Pbbb6/cxHzxxRc5efIk27dvZ9OmTTg7O7N48WI++eQThg8fTufOnfn000+59957ad++Pffddx8eHh588sknTJkyhdLSUnr37s2cOXOqnbNi07hZs2b8+uuvTJo0icLCQqSUvPXWW0b8dOtG1PXIYC7i4uJkgxtcHPgevp0BczZBUFeT2qXR2AoJCQl07tzZ2mYYzYwZM5gwYQI333yztU0xOSkpKUyYMIH9+/dbdN7q3gNCiB1SyrjqxtupyyVMvWo/ukaj0VRivy4X0JEuGo0NsXDhQmubYDYiIiIsvjpvCPa5QvcKBJdmeoWu0Wg0VbBPQRcCfEO0oGs0Gk0V7FPQQcWia5eLRqPRVGLHgh6mY9E1Go2mCvYr6L6hUHAGigusbYlGozEBb7/9dmXBK3OzYMEC/ve//1lkrgp0tcXa0JEuGo3DUFZW1mhBP3funNFj58yZw7Rp0xo8V33R5XProkLQtdtFozELn3/+eWX52nvvvZeysjK2b99Ot27dKCws5MKFC3Tp0oX9+/ezdu1ahgwZwg033EBMTAxz5szBYDAAKk0+NjaWrl278tRTT1Xev3nz5sydO5e+ffvy0ksvkZGRwfDhwxk+fDhlZWXMmDGDrl27Ehsba1Qm5eTJk5k4cSI//fTTZfVlquO5556rLLw1bNgwHnvsMYYMGULnzp3Zvn07N954I+3bt+eZZ54B0OVzzU5lowsd6aJpAvz6NJzaZ9p7BsXCuH9VeyohIYGvv/6aTZs24erqyv3338+iRYuYNm0aEydO5JlnnuHixYvceeeddO3albVr17Jt2zYOHjxIeHg4Y8eOZcmSJQwYMICnnnqKHTt24O/vz5gxY/jhhx+YPHkyFy5coGvXrjz//PMAfPzxx6xZs4bAwEB27NjBiRMnKmO/jSkvu3btWtavX8/HH3/M448/zpQpU7j77ruJjo6u81o3NzfWr1/PvHnzmDRpEjt27CAgIIB27drx2GOPAejyuWbFOwicXHTookZjBlatWsWOHTvo3bs3PXr0YNWqVSQnq7aPc+fOZcWKFcTHx/Pkk09WXtOnTx+ioqJwdnbm9ttvZ+PGjWzfvp1hw4bRsmVLXFxcmDp1amV5WWdnZ2666aZq54+KiiI5OZmHHnqI3377DR8fnzptFkIwdOhQPv30U3bu3ImTkxOdOnVi8eLFdV47ceJEQJUD7tKlS2Wp4KioKI4fV14AXT7XnDg5g0+wdrlomgY1rKTNhZSS6dOn88orr1x1Ljs7m/z8fEpKSigsLKwsEVvf8rIeHh44OztXe87f3589e/bw+++/M3/+fL755hs+/vjjyvNlZWX06tULUGJcscq/ePEi33//PR9//DE5OTnMmzeP0aNH1/n9Vi0HfGWp4Ar3jS6fa278wvSmqEZjBkaOHMl3331HZmYmoES8ogHEPffcwwsvvMDUqVMv84lv27aNY8eOYTAY+Prrrxk0aBB9+/Zl3bp1nDlzhrKyMr788kujystWtIC76aabeOGFF9i5c+dlY52dnStL51aI+ZNPPklMTAybNm3i3//+N/Hx8TzwwANGre6NQZfPNTd+YXC08Z9qGo3mcmJiYnjxxRcZM2YMBoMBV1dX5s+fz7p163BxceGOO+6grKyMAQMGsHr1apycnOjfvz9PP/00+/btq9wgdXJy4pVXXmH48OFIKbnuuuuYNGlStXPec889jBs3jjZt2vD2228zc+bMyo3V6p4UrmTYsGE8//zzeHh4mPRnUYEun1sLjSqfW8GaV2Ddq/BMJri4mcYwjcZGsKfyuVf2A3U0dPlcS+AXBkg4n25tSzQajcbq2Lmgl4cu6o1RjcaqDBs2zGFX56DL51qGilh0HbqocVCs5RLVWJ+G/O7tW9B9gkE46UgXjUPi4eHB2bNntag3QaSUnD17tt4bvPYd5eLiBt5ttMtF45CEhISQnp5OVlaWtU3RWAEPDw9CQkLqdY19Czoot4t2uWgcEFdXVyIjI61thsaOMMrlIoQYK4Q4JIRIEkI8Xc15XyHEz0KIPUKIA0KImaY3tQb8QnU9F41Go8EIQRdCOAPzgXFADHC7ECLmimEPAAellN2BYcAbQgjLBIb7hcH5DDCUWWQ6jUajsVWMWaH3AZKklMlSymLgK+DKVC8JeAtV7KA5kA3UXr/SVPiGgqEU8k5aZDqNRqOxVYwR9GCg6q5jevmxqrwLdAYygH3AI1JKg0ksrAs/Hbqo0Wg0YJygi2qOXRlHdS2wG2gL9ADeFUJcVRFHCHGPECJeCBFvsp17v3D1qiNdNBpNE8cYQU8HQqt8HYJaiVdlJrBEKpKAY0CnK28kpXxfShknpYxr2bJlQ22+HN/ysB69MarRaJo4xgj6dqC9ECKyfKPzNuCnK8akASMBhBCtgY5AsikNrRHXZuDVUrtcNBpNk6fOOHQpZakQ4kHgd8AZ+FhKeUAIMaf8/ALgBWChEGIfykXzlJTyjBntvhzfUO1y0Wg0TR6jEouklMuAZVccW1Dl/xnAGNOaVg/8wuC07RfO0Wg0GnNi37VcKvArX6EbLBNYo9FoNLaIYwi6bxiUFcEFXfNCo9E0XRxD0P3C1KuuuqjRaJowDiLoOrlIo9FoHEPQdaMLjUajcRBB9/ABD1/tctFoNE0axxB0UH50HYuu0WiaMI4j6L5h2uWi0WiaNHYn6JuSzjB5/iayLxRffsIvVLlcdP9FjUbTRLE7QXd1dmL38Rx2pp67/IRfGBTnw8Vz1V+o0Wg0Do7dCXq3EF9cnQXxVwp6RaSL3hjVaDRNFLsTdA9XZ7oG+xKfkn35CR2LrtFomjh2J+gAceH+7D2RS1FplT6iutGFRqNp4tiloPcKD6C41MD+E7mXDjbzB1cv7XLRaDRNFrsU9LgIfwDiU6r40YUor7qoXS4ajaZpYpeCHtjcnchAr+o3RrWgazSaJopdCjpAr3B/dqSeQ1aNO/cL0y4XjUbTZLFbQY8L9yf7QjHJZy5cOugXquLQi/KsZ5hGo9FYCfsV9IgAAHZU9aNXVl3Uq3SNRtP0sFtBb9fSC39PV+JTq8SjV4QuareLRqNpgtitoAsh6BXuf3mki04u0mg0TRi7FXRQ8ejJZy5wNr9IHfBqBc5uWtA1Gk2TxK4FvSIefUdF+KKTE/iGaJeLRqNpkti1oMcG++Lm7HRJ0EE3utBoNE0WuxZ0D1dnYkN8L08w0slFGo2miWLXgg4qHn1fei6FJeWFuvzC4EImlBRa1zCNRqOxMHYv6L3C/SkuM7CvolCXX5h6zU23nlEajUZjBRxC0KFKoa7KRhfa7aLRaJoWdi/oLZq7E9XSix0VCUY6Fl2j0TRR7F7QQfnRd6Sew2CQ4N0WhLNjR7pICYnLoCjf2pZoNBobwkEEPYBzBSUkn8kHZxfwCXbsWPR1r8FXt0P8x9a2RKPR2BAOIei9rmx44ciNLvZ8DWtfVv8/utq6tmg0DsyFotJLWeh2gou1DTAFUYFeBHi5EZ96jtv6hKmN0ZSN1jbL9KRsgh8fgIjBENgBdi9S4ZmuHta2TKNxKKSUzFy4nW3Hsuke6sfozq0Y2bk1nYK8EUJY27wacYgVekWhrsqMUb8wyMuAshLrGmZKzhyBr+6AgEi49TPocC2UFsLxP6xtmVW4UFTK19vTOHRK177XmJ7NR8+y7Vg247oGIYDXlx9m3LwNDHp1Df/4cT8bjmRRXGqwtplX4RArdFAboysOnuZMfhGBfqEgDXD+BPhHWNu0xnPhDCy6GZxc4I5vVEPs8AHq6+Q1EDXU2hZajDP5RXy6OYX/bUkl92IJceH+fHffAGubpXEgpJTMW3mEIB8P3r6tB+4uzmTmFbImMZMVBzP5Ov44n25Jpbm7C0M6BDKqc2uGd2yFv5ebtU13IEGv4kcfW7XRhb0LekmhWpnnnYLpS9UKHcDdG0L6QPJaq5pnKVLOXOCDDcl8tyOd4jIDozu3xt/Tja/jj5OclU9Uy+bWNlHjIGxJPsu2lGz+ObEL7i7OALTy9uDW3mHc2juMwpIyNh89w4qDmaxKOM2yfadwEionZlTn1ozs3Jp2Lb2s4ppxGEHvGuyLm4sTO1KzGduvIlvUziNdDAb44T7lVpmyEEJ7X34+ahisfQUKssEzwBoWmp1daed4f30yvx04hauTEzf1Cmb24CjatWxO5vlCvtuZznc70nlybCdrm6pxEP6z6gitvN25tXdotec9XJ0Z0ak1Izq1xmDoyv6MXFYmZLLy4Gle+TWRV35NJDLQi5GdlN+9d4Q/Ls6W8W47jKC7uzjTvaJQ17Vx6qC9x6KveREOLIFR/4QuN1x9PmqYing5th66TLa0dWbDYJCsPZzJgnXJbDuWjY+HC/cPa8f0ARG08r60AdzKx4OhHVqyeGc6fx7TEWcn292s0tgHfySfZWtyNnMnxODh6lzneCcnQbcQP7qF+PH46A5k5FxkVaIS9/9tSeXDjcfwbebKsI4tGdm5NUM7tMS3mavZ7HcYQQfV8OKjjckUShc8mgfZd+jizs9gwxvQczoMfKT6McG9wM1b+dEdQNCLSw38uPsEH2xI5vDpfNr6evDM+M7c1ieM5u7Vv1Wn9AphdWIm649kMbxjKwtbrHE05q06QmBzd+7oG9ag69v6NeOufuHc1S+cC0WlbDhyhpUJp1mTmMmPuzNwcRL0iQxgWv9wxnZtY2LrjRR0IcRYYB7gDHwopfxXNWOGAW8DrsAZKaXFd+riwv1ZsE6y53gOff1C7beey9E1sPRRaDcCxr8BNfninF0gcrDd+9HzCkv4clsaH29M4dT5QjoFefPWrd2Z0K0trnU8qo7s3Bp/T1e+i0/Xgq5pFNtTstl89CzPjO9s1Oq8LrzcXRjbNYixXYMoM0h2Hz/HygTld0/LLjCBxVdTp6ALIZyB+cBoIB3YLoT4SUp5sMoYP+A9YKyUMk0IYZW/rMpCXann6OsXBid2WsOMxpGZAN9MU3HmUxaCcx2PZ1HD4dAyyD52acPUTjh9vpCPNx3ji61p5BWVMqBdC169uRtD2gcavaHk5uLEpB7BfPFHGjkFxfh5Wj/SQGOf/GfVEQKbuzG1b7jJ7+3sJOgVHkCv8ACeGtuJMoM0+Rxg3Aq9D5AkpUwGEEJ8BUwCDlYZcwewREqZBiClzDS1ocbg7+VGdKvmKh49OBQO/qQ2Fp3sJNw+PxMW3QKuzVR4oodvjUPLDFL5jKOGqQPJa+1G0JMy83h/fTLf7zpBmUEyLrYN9w6JoluIX4Pud0tcKAs3p/Dj7gymD4gwqa2apsGO1HNsOHKGv13XiWZujV+d14W59nuMUbpgoOruYnr5sap0APyFEGuFEDuEENOqu5EQ4h4hRLwQIj4rK6thFtdBZaEu31AwlED+KbPMY3KKC+CLW6HgDNzx9aWqkdVwNCufPi+t5Ptd6RDYXtWuSV5jQWPrj5SS7SnZzP50O6PeXM9PezK4vU8Ya58Yzvw7ejZYzAFi2vrQpa0P38Tb+Sa4xmrMW3WEAC837uxn+tW5JTFG0Kv7KLnyecEF6AWMB64FnhVCdLjqIinfl1LGSSnjWrZsWW9jjaFXuD+5F0s4Sfn97SHSxWCAJX+CjF1w04fQ9ppah7/0SwJnLxTz8rJECkrK1Cr92HowlFnG3nqy53gON/53M1MWbGFH6jkeHdWezU+P5PlJXQlr4WmSOab0CuFAxnkOZpw3yf00TYddaedYfziLPw2OwtPNvuNEjBH0dKDqcjEEyKhmzG9SygtSyjPAeqC7aUysH3ERKh5713kfdcAeYtFXzoXEpXDty9BpfK1D1x3OYnViJtd3b0tWXhEfbjim/OgXz8GpvRYy2HgKS8p44IudnDh3kecndWHz0yN5dFQHAkycVTepRzBuzk58u8MOft8am+I/q47g7+nKtP72vToH4wR9O9BeCBEphHADbgN+umLMj8BgIYSLEMIT6AskmNZU44ho4UkLLzc2nimPV85JtYYZxrP9I9j8DvS5B/rdV+vQ0jIDLy49SEQLT96Y0p1ru7Tm/9Yd5WzrfmrAUdtzu3yyKYX0cxd569YeTOsfYTb/pL+XG6NjWvPDrhM2WWNDY5vsOZ7DmkNZzB4chVcNobH2RJ2CLqUsBR4EfkeJ9DdSygNCiDlCiDnlYxKA34C9wDZUaON+85ldM0II4iL82XK8EJoF2LbL5cgKWPYEtL8Wrn2l5vDEcr7YlsaRzHz+dl1n3FyceHJsJwpLDfxnay607mpz4YtZeUXMX5PEqM6tGBgdaPb5bo4L4VxBCasSTpt9Lo1j8J9VR/DzdHWYzXSjwj+klMuklB2klO2klC+VH1sgpVxQZcy/pZQxUsquUsq3zWSvUcSFB5B6toASn1Dbdbmc2gffzoDWXeDmj1VMeS3kFBTz5orDDGjXgtExrQFo17I5t/UOZdEfaeS2GQhpW6HkogWMN443VxymsKSMv13X2SLzDWnfktY+7ny7QzcI19TN/hO5rErM5O6BkTUmrtkbdhLPVz8qGl6cdWllm9mi5zNUeKK7jwpPdK+7sNS8VUc4f7GEZyfEXBaj/cio9ri5OPF5ZiSUFUHaFnNabjSJp87z9fY07uofbrHCWc5Oght7hrD2UCaZ5wstMqfGfpm36gg+Hi5MHxhhbVNMhkMKete2vri7OJFS2kK5XKR5gvgbRFG+Ck8sOg9TvwGftnVekpSZz2dbUrmtTxid2/hcdq6VtwezB0fxbnIrDE6uNuFHl1Ly4tIEvD1ceWRke4vOPaVXCAYJS3adsOi8GvviQEYuKw6e5u5BUfh4mK+2iqVxSEF3c3Gie4gf+/J9ofQiFJy1tkkKQxksvhtO71dZoEGxRl320i8HaebqzJ9HXxUJCsA9Q6Lwau5DgksnpA340VcnZrIx6QyPjmpv8czNqJbNiQv359v440hb+iDX2BT/WXUEbw8XZjjQ6hwcVNBB1UffkVP+qG8rkS6//RUO/wbX/RvajzbqkrWHMllzKIuHR7anRXP3asc0d3fhkZHtWXahM+LUXtUQw0qUlBl4aVkCUS29rJakMSUuhKNZF9iZlmOV+TW2TcLJ8/x+4DQzB0aatfKhNXBoQU8ztFBf2EKky9YFsO3/oP+D0Hu2UZeUlBl48ZcEIlp41rkLf1ufMI75qLLBhuR1jbW2wXy+NZXkrAv8/brOdRbWMhfju7Wlmasz3+mYdE01vLP6CN7uLtw90D5KZdQHhxX0nmH+pMvyUDlrR7okLoPfnoZOE2D0C0Zf9sUfaSRl5vP38TG4udT+q3J1duL6ceM5Lz1J2ba0sRY3iJyCYt5eeYRB0YGM6GS9yofN3V0YFxvEz3tOcrHYNrNnNdbh0Kk8lu07xYyBEfh6OtbqHBxY0P083WjdsjUFwtO6K/SMXcpv3vYauPEDowuF5RQU89bKwwyKDmRUZ+PEcWxsMAc9etDs+HouFpU2xuoGMW/VEfIKS3hmQmerd0a/JS6U/KJSfjtw0qp2aGyLd1YfwcvNmVkOuDoHBxZ0gLjIFqTLQKQ1fejfzwHPFnD7V+BmfN2St1eqMMX6iKMQgqAe42jDGRav2tBQixvE0SwViXNr7zA6BfnUfYGZ6RsZQFiAJ9/G65h0jeLI6Tx+2XeS6QMibKKhszlwbEEP9yetrAVFZ60k6DnHISsR+t0P3q2NviwpM4/PtqZye5/6i2NEH1UL5tgfS8m+UFyvaxvDK8sS8HB15vEaInEsjRCCm3uFsPnoWY6bqZmAxr54Z3USzVydmT04ytqmmA3HFvQIf9JlS5xyrbRKO1a+ORlVv+ZNL/6SgKdbA8UxIIoS7xB6G/bw7uqk+l/fADYlnWFlQiYPDI+mpXf1kTjW4KZeIQgB3+nM0fpTctFmq3c2hKTMfH7em8G0/hEmLwxnSzi0oIcFeJLjFoRbaR5czLG8AcfWg2cgtDQ+9X3NoUzWHsrikVrCFGtFCFyjhzPUNYFFW5NJO2ve1WmZQfLC0oOE+Ddjpo3F9Ab7NWNgu0C+25GOwUwdYhyS4gJ4pxf8+qS1LTEZ89ck4eHizJ8GO6bvvAKHFnQhBJ4tI9QXlo50kRKS10HkEKM3QkvKqylGBnoxrX9Ew+eOGkYzQz7dnI7x+vJDDb+PEXwTf5zEU3n8dZxp+jCamilxIZzIucjWZBtJLrMHdnwC509A/Mdw+oC1rWk0yVn5/Lj7BHf1D2/YIsmOcGhBB2gVptwWOSePWnbiM0dUt6TIIUZf8vnWVI6Wx3DXFaZYK+Vt6R6OSOenPRnsS89t+L1qIa+whDeWH6J3hD/XxQaZZY7Gcm2XILw9XHQ3I2MpuQib5kFIb3D3hhVzrW1Ro3l3TRJuLk78yYF95xU4vKBHt1fujhMphy07cT395+cuqBjuwe0DGWlkmGKNeAVCUCz9xT4CvNx45dcEs6TBv7f2KGfyi3lmfIzVwxRrwsPVmYnd2/Lr/lOcLyyxtjm2z87/Qf5pGPUcDHkSklZC0iprW9VgUs5c4MfdGdzZN9ym9nfMhcMLeseoSC5KN86fTLbsxMlrwTcM/I3z2VXGcJtKHKOG45K+jceGBrP56FnWHTZtD9fj2QV8tPEYN14TTPdQP5Pe29RMiQulqNTA0j06Jr1WSotg49sQNgAiBkGfP4F/BCx/1m43SOevScLFSXDPUMdfnUMTEHRXF2eyXVphsGQZXUMZpGxU7hYjxPnIaRWmeEffMDoGeZvGhqhhYCjhttbphAV48q9fEykz4cbgv35LxEnAX8Z2NNk9zUX3EF/at2qu29PVxa7PIS8DhpZvhrq4q5V65gHYvciqpjWEtLMFLNl1gjv6htHK28Pa5lgEhxd0gOLmIXgXnqSg2ELZk6f2QmGO0e6WS2GKJhTH8AHg7I5ryjr+cm1HEk/l8b2JSsruSM3ml70nuWdIO9r4NjPJPc2JEIJb4kLZlZZDUmaetc2xTUqLYeNbENKncg8GgJjJ6tjql1TpZzti/poknJ0Ec4a2s7YpFqNJCLpbYDjBIovdx3MsM2FFcSwjNkTXJGay7rAKUzRpfKxrMwjrC8lrGR/bhm4hvry5/BCFJY17dDYYJM8vTaC1jztz7OgxdvI1wTg7CZ05WhN7vlSRYEOfuvypUgi49iW1wb/5HevZV0+OZxeweGc6t/cOpbVP01idQxMR9BbB0bQQeew5mmGZCY+th8CO4F175EdJmYEXfjlIVGPDFGsiajic3o9TQRZPj+tERm4hn25OadQtf9qTwZ7jOfzl2k54utlP266W3u4M79iKJbtOUFqmm0hfRlkJbHhD1RuKHnn1+dA+0OUG2PwfOG8f+xDvrT2KkxDMGdZ0VufQRATdI1BtTKYeM29MNqAeXdO2GOVu+WxLeanZ8Y0MU6yJikfn5HUMaBfI8I4tmb8miZyChpUEuFhcxqu/JRIb7MuN1wSbzk4LMSUuhKy8IpNvENs9+75VPQOuXJ1XZeQ/lPCvedGytjWAEzkX+W7HcW7tHWoXLkFT0iQEHb9QAHIzjpp0Y7Ba0rdDSQFE1i7oKkzxMIPbm7HUbJvu0MxfRdwAT43rRF5RKfPXNKwkwAcbkjmZW8izE2JwcrLNMMXaGNGpFYHN3bTbpSplpbD+ddU9q8PYmscFRELfe2HXIji133L2GcPZo7Dtg8pWk++Vv7/va2Krc2gqgu6rBD2g9DSHT5t5U+zYehBOEDGw1mFvrTxMflHpVU2fTYqTs/LjJ68BKekU5MNNPUP4dHMq6efqVxLg9PlC/rv2KOO6BtEnMsA89poZV2cnJvcIZlXiac7mF1nbHNvgwBLIPlr76ryCIU+Ahy+seNYythlDYS4suhmWPQEpG8jIucg38ceZEhdKW7+mtTqHpiLo3kFIJxeCxRniU8+Zd65j6y6tjGvg8Ok8Fv2RxtS+4XRobaIwxZqIGqbSuM+qVcvjozsgBLy5vH6JVv/+/RBlBslfxxlfl8YWmRIXSkmZ5IfdFtpPsWUMZbD+39AqBjqOr3t8M38l/EdXw5GV5revLqSEHx+Ec6ng1hy2f8SCdSoj/P4muDqHpiLoTs7gG0I712x2pGSbb57iC8rlUou7RUpVzMrLzZnHLFFqNmq4ej26BoC2fs2YOTCS73ef4ECGcSUB9p/IZfHOdGYOjCCshfE13W2RjkHedAvx1U2kAQ7+CGcOw5C/GF1viN6zVbLc8meUu8aa/LEAEn6CkXOh1wxk4lJWbtvHzb1CCPG37/dpQ2kagg4I31Dau2ebd4WeugUMpbWGK645lMmGI2d4ZFQHy5TxDIgEv/BKPzoo36JvM1f+9WtinZdLKXl+6UECPN14YES0GQ21HFN6hZB4Ko8DGeetbYr1MBjU6jywI8RMMv46FzcY/U/ISoDdn5vPvro4vl19qHQYBwMehrhZCEMpN4nV3D/MMd6nDaHJCDp+YbSWWaSfu8ip3ELzzHFsLTi7QVj/ak+raooJRLX0Ylr/cPPYUB3thkPKhsoVlW8zVx4cHs2GI2fYcKT2iI/fD5xi27FsHhvdAR8Px+jBOLF7MG4uTnzblAt2JS6FzIPlq/N6VsnsPBFC+1kv2aggG76dAT5t4Yb/gpMTma7BbDLEMqvZOkL9HL9mS000KUH3LDqDGyXEp5rJ7XJsvcqqq6HV3P+2pJJ85gLPjO+Mq7MFf/RRw6DoPGTsrDx0V/9wQvyb8a9fE2usFV5UWsbLyxLp0Lo5t/UOtZCx5sfX05VruwTxw+6MRida2SVSwvrXoEU0dL2x/tdXJBtdyFSVGS2JwQBL7lFzT/m0cq9qwbpkFpWNxL/kNBxZblmbbIimI+i+oQgkka45xKeYwe1SkA0n99bobsm+UMy8lYcZ0qElwzuaKUyxJiKHAqLSjw7g7uLME2M6ciDjPD/tqX6D8NPNKaRlF/DM+BhcLPkBZAGm9Aoh92IJKxNOW9sUy3P4Nzi1Dwb/uf6r8wpC4qDLjSp79LzpN5illJSWGSgoLiW3oISsvCIyci5ybvmrkLSCk/3msscQRXxKNmsOZbLoj1Sad58IzYNg+0cmt8desJ9Uv8ZSHos+uNVF/jCHHz1lAyBrTCh6a8VhLhSX8cx445s+mwzPABV5k7wWhj1VeXhi97Z8sCGZ15cfYlxsEO4ul/64z+YX8c6qJIZ3bMmQDi0ta68FGBgdSBtfD76NT2dCt7bWNsdySAnrXlVVFGOnNO5eo/6hXDerX4TJ79U5PPtCMX//fh+nzhdSXGqgpMxASZms8v/yr8v/f+WedX+nA3zu+io/Gfrz8KpIWLWp8pyLk+D+EZ1g33RY9xqcS1HfYxOj6Qh6eSx6b798Pkk4z4WiUrzcTfjtH1sPrl4Q3OuqU4dO5bHoj1Tu7GeBMMWaaDdcraaK8sG9OQBOToK/juvMnR/9wWdbUi9rnvvWysMUlJTx9/H2HaZYE85Oqon0/DVJnMotJMi3idT7SFoJGbtg4jvg3Mg9Ef8IlWy0+V3oOwfadKt1+L9/P8SKg6fp364F7l5OuDpf+ufmIi7/2lmUH1df+5SeZfyWhyhwjcBl2Dss8PC+7JrWPh5EBHpBz+kqUWrHQlUpsonRdATdJxiEEx2b5VBmkOw+nsPA6EDT3T95XXmFw8v/SKSUvPjLQZq7u/DYKAuEKdZE1DBVTS91E3S4tvLwoPaBDG4fyLtrkpgSF4pvM1cOn87jiz/SuKtfONGtrPQBZAFu7hXCO6uTWLwznQeGN4HICClh7b9Unf5ut5nmnoOfUGV3lz8D036sMTnpQEYuX21PY+aASOZeH1O/OcpK4X8PQlkBzFrKda3b1zzWNxg6joOdn8Gwv6oSwE0Ix3KM1oaLG3i3oS1ZCIFp/ejnM+DskWrdLasTVZjio6M64G/NbuOh/cDF47LwxQqeHteJ3Isl/HetSsp48ZcEmru78Kg1P4AsQHgLL/pEBjSdmPTkNXAiHgY9qv4eTEEzPxj6tEqoO7Ki2iFSSp7/+SD+nm48MrIWMa6JNS9B6kaY8Ba0NuLDIG4mFJyBhJ/rP5ed03QEHcA3FLe8E3Rs7W3aSJdj69XrFRuiJ3IuMvfHA0S19OIuS4YpVoerhwqnrLIxWkGXtr7c0COYTzYd48ttaaw/nMXDI9tb9wPIGAqyVQ/MRjClVwgpZwvMn0FsbaRUvmXvtnDNnaa9d9wsCGinSgJUk2z06/5T/HEsm8dHd8DXs55unsO/w8Y3oec06HG7cddEjVDJT01wc7RpCbpfGOSmERfhz660HNMV6kpeB80CoHVs5aFTuYXc8cFWzheW8J/brrFsmGJNRA1TCSF5p6469fiYDkgJf12yj0hzlfM1Jfu+gzdj4Le/Nuo218W2wdPN2fFj0lM2qiqggx4zvRuiMtkoEXb977JThSVlvLwsgU5B3tzeJ6x+981JUyGKQbEw7jXjr3NyUqv0tM2QmVC/Oe0cG1AZC+IXCrkn6B3qQ35RKYmnTJApKKVaoUcOrkyfzswr5I4Pt3Imr4hPZ/Wha7Bv4+cxBe3KywBU43YJ8fdk+gD1FPHXcZ3MU87XFJSVwu9/h8V3q6zcxF9UbHID8XJ3YUK3Nvyy14IdrazBuldVSF/Paea5f6cJqhfpmpeh6FIBvA83JJN+7iJzr4/BuT4VOkuLVfKQNKh4c9d6FtrqcSc4u0P8x/W7zs6x0b9aM+EbCrKMPoGq0t4OUzxmZyfD+fRKd8vZ/CLu/PAPTuYUsnBWH3qG1Vyky+K0jgXPFtUKOsCTYzux5P4BjI5pbVm7jKUgGxbdBFvehd5/gvFvqAST0/saddspcaFcKC5j2b6rn1xsnTWJmTXmEVSSukWF1Q58WLnezIEQMOZFuJClGk2jnlLfW3uUsV2CGNCungEIy5+BEztg0nxo0YBCW14toMtk2POV3bXOawxNS9D91CNfkMwiyMfDNBujFeIYOYycgmLu/GgbqWcL+GhGHL0jbKzMrJOTSjI6uoargnxR5WV7hvlbPk7eGE7tg/eHQupmmPgujH/9Uv3upMZV/osL9ycy0Itv7MjtIqXknVVHmLlwO498tYstR8/WPHj9a+DVEnrNNK9RIb2g683qAzf3BK/9lkhpmeRv19Uz9PXA97Dt/6Df/RAzseH2xN2tMqT3f9fwe9gZTVLQRe5xekX4m2aFfmw9+AST6xnGXR9t42hmPh9Mi6v/isRSRA1T/SGzLNC9yVTs+w4+HK3cLTN/g553qePerSGoGyStatTthVAx6duOZZN69oIJDDYvRaVlPP7NHt5YcZjJPdoS0cKLx7/ZTW5BydWDj29X5W4HPFRjSQqTMnIuSMnZn59hya4TzB4cWb8KnWeS4MeHIKQ3jPpn42wJ7QOtuqjN0aYQxURTE3TfEPWac5y4cH9O5FwkI6cRURIGAxxbT0nYIGYs3E7iqfMsuKunbWdWVvrRr452sTkMZbD8WeUvb9sD7lmrVoFViR4Fx/+Awsbth9zYMxgnAd/tsO1uRmfzi5j6wR98v+sET4zpwFu39uDtW3uQlVfE377fd3X45frX1IZ93N2WMdA/HNl3Dv5J3zOo+Qnur098f8lF+Ha6yuWYsrDxoZVCQO9ZcGovnNhZ93gHwChBF0KMFUIcEkIkCSGermVcbyFEmRDiZtOZaEJcm6lHz9y0SndIo8LVTu+Hi9ksOB7CvvRc3r2jJyM62aj/uQK/MAiIqtGPbjMUZMPnN6nGxHF3w7Sf1Ir8SqJHqc3RitDRBtLGtxmD27dk8Y5087cpbCBHTucx+b1N7DuRy/w7evLgiPYIIege6sdjozvwy76TLN554tIFJ3aqQlUDHqzMDrYES/1uJ0d68abvtzR3q0etmGVPwOkDcOMHlxZfjaXbrar5RXzTCGGsU9CFEM7AfGAcEAPcLoS4Krq/fNyrwO+mNtKk+IVBznE6BXnj6ebcqIYXJUlrAfgyM4J5t13DtV2CTGSkmYkarsLYyqp5RLcFTu2H94eprNaJ78CEN2terYX2ATfvRvvRQTWRzsgtZFPSmUbfy9SsO5zFje9t5mKxga/v7c/4bm0uOz9naDv6Rgbwjx/3X3Ibrf83ePipDWQLcaGolBdXnuDb5nfS6uw2FUduDLsWqYzTIU9A+1GmM8jdG7rdAvsXq0WCg2PMCr0PkCSlTJZSFgNfAdVVxH8IWAxkmtA+0+MbCjlpuDg7cU2YX4NX6IUlZRzY9DPJhjY8ecvIq/7AbJqoYVCcD+nx1rbkavYvgY9GQ1kxzPy17jA7Z1eVoZu0qtF+0lGdW+PbzJVvbczt8tmWFGYt3E6wfzN+fHAgPUL9rhrj7CR469YeODsJHvlqNyUn9sChZWpj0cPHYrYuWHeU0+eLiLv5z6o8bw3JRpdx+gD88mcVKTascXkF1RI3C0oLYc+Xpr+3jWGMoAcDVbf/08uPVSKECAZuABbUdiMhxD1CiHghRHxWVu2NFcyGXyjkpoPBQK/wABJOnie/qH7xx8WlBh76fBvRF/diiBzC5GuC677IlogcrBpZ25LbxVAGK+bCdzPVRuc961SJVmOIHgm5aXDmSKNM8HB1ZlKPtvx+4BQ5BcWNupcpKC0z8NxPB3j2xwMM69CS7+4bQHAtjY/b+jXj5Rtj2X08h2NLngN3H1U8y0Iczy7g/9YnM6lHW3pFtoLRz6sWdzsX1nxR4Xn4ZppqPn3TRw0v51sbQbGqT0H8xw6/OWqMoFcXw3blT+Vt4CkpZa3dAqSU70sp46SUcS1bWmnj0C8cyorgQhZx4f4YJOxOyzH68pIyAw9+sZOzh7fSXBQS3deI5rq2RjN/aHuN7WyMFmSrzu2b5qnV1PSfq/eX10S7kerVBG6XW3uHUlJm4Lp5G/jijzSKSxuetNQY8gpLmP2/eBZuTmH2oEjenxZHcyOqg07o1pYHYorpcHY1GZ2mq1orFuJfvybiLARPj+ukDnS8DsIHwppXqt+0lhJ+fljlctz8MTQ3Y5+A3nerRumN3GuxdYwR9HSgaruaEODKTIY44CshRApwM/CeEGKyKQw0OeVldMlJ45owP5wEbDfSj15aZuDRr3ez/OBpnonJBAREDDafreYkaphyuTQyOqTRnD4AH5T79K//jyrAVN/oBv9wCOxgEkHv0taXRXf3pbWvB3/7fh8j3ljLN9uPU1JmOWE/nl3ATf/dzMYjZ3j5hliemVC/LMtH3X+gAA9mJfbmfKFl9km2Jp/ll30nmTO0HW18y58iKpKNCs6oSp9Xsu0DFXM+4lmIGGheA2Mmq2gfB98cNUbQtwPthRCRQgg34Dbgp6oDpJSRUsoIKWUE8B1wv5TyB1MbaxLKG12Qm4a3hyudgnyMikcvM0j+8t1eftl7kr9f15meZfvUo5ynjSUPGUvUcJBlSkitxf4l8OEoKCmEGcug1/SG3yt6lNpEbWSxLoAB0YEsuW8AC2f2poWXG08u3suoN9exZGc6pWYW9h2p55g8fxOncgv5dFYf7uhbz/onWYdxTfiR3NiZHMlzZe4P+81jaBXKDKqaYrBfM+4ZEnX5yeCeEHsLbH0Pcqp4bk/sgN//ppLDBj5qdhtx9YBrpqpSEdXUMnIU6hR0KWUp8CAqeiUB+EZKeUAIMUcIMcfcBpqcyhW6enOpQl3nav1DNRgkf12yl+93neAv13bkT/2CVOxzDd2J7ILQPuDqaR0/uqEMVvyj3F8eC/eug9Dejbtn9Ei18ZWyqe6xRiCEYFjHVvzwwEA+mq7cHY9/s4cxb6/nx90nzBLa+OPuE9z+wVaae7iw5P6BDavXv+F1cG1Gm7FP8PCI9vywO4Mfdp2o+7pG8E38cQ6ePM/T4zrRrLowxfJkI1a/oL4uyIZvZoB3G5j838oaSGan10wV4rrzf3WPtVOM+klKKZdJKTtIKdtJKV8qP7ZASnnVJqiUcoaU0nZzbT18VChXrhL0XuH+XCguI/FUXrXDpZQ8++N+volP55GR7VUjhONbVRRGpB0Luou7ashhaT96QTYsmgKb3lZ/YNOXgrcJwj3DB6p67yZwu1RFCMHIzq1Z+tAgFtzZCzdnJx75ajdj317PL3tP1thguz5IKXlrxWEe+Wo3PUL9+OH+gUS3akDc+NmjsO9b5S/2CuSB4e2IC/fn2R/2czy7oNF2Vsf5whJe//0QvSP8mVBTpJdfKPS/H/Z+rVbmP9wHeSdV8pAln3BbtFNPpjsW1h15Y6c0rUzRCvxU6CJAXEWCUTV+dCkl//z5IIv+SOO+Ye14dFR5cf5j68HJRdUXt2eihqkohFzzruAqqfCXH1sP18+D6982XaMF12YQMcjkgl6BEIKxXYNY9vBg5t/REwk88MVOrvvPBn7bf6rBDTIKS8p4+KvdzFt1hJt7hfD53X0bXod+wxvg7Ab9HwLAxdmJt27tAcBjX+82i7vonVVHyC4o5h/Xd6m9BtCgx8EzUCWLHf4Nrn356qxfS9D7bjh/QiVcOSBNU9B9wypdLsF+zWjr63FVPLqUkld+TWTh5hTuHhTJk9d2vPSGTV4HwXEWzb4zC1E1l9M1OQd+UPVYSgph5jLoNcP0c0SPUp2jzqWa/t7lODkJxndrw++PDmHebT0oLjUw5/MdTHhnI6sSTtdL2LPyirj9g638vCeDp8Z24t83d2t42eJzKaqyYK+Zl0UIhQZ48sLkrsSnnqvsSGUqkrPyWbg5hVt6hdZdItrDB4Y9DRfPQZcboI/lkp0uo8M41eTDQTdHm6ag+4Upl0v5H1+viICrNkbfWH6Y99cnM61/OM+M73xJzC/mwMnd9u0/r6BVjCqFYE5Bv3gOlv1F1eho3UXVYwntY565osszDI82rliXMTg7CSb1CGb5Y0N4Y0p38gpLufvTeCbP38TaQ5l1CnviqfNMnr+JhJPnWXBnT+4b1q5xVS43vKmeGgc+ctWpydcEM6lHW95edYRdaabrzPTSLwm4uzjzxLUdjbsgbhbc+rkqiWutip7OLmrzPWkVZB+zjg1mpIkKeqjKlLyo3txx4f6czC3kRHmhrv+sOsK7a5K4vU8oz135KJm6SRXdt2f/eQVOTsrtkrzW9AkXhjKI/wTe6QXbP1Rd4WcsBR8zZtS2iFYf1o2svlgfXJyduKlXCKv+PJTXburGmfxiZnyyvTLssDphX5OYyU3vbabUYODbewcwtmsjfyY5abD7C5VVW8PP9/lJXQny8eDRr3fXO5GuOtYdzmJVYiYPjYimpbeRHZCcnKHz9eDm1ej5G0XPaSqxbscn1rXDDDRNQa8Siw5qYxSUH/2/a4/y5orD3NQzhJcmx+J0Zfxv8jpwaWZ8FqOtEzVMNYnIPGi6e6ZtVbVYlj4KLTvBveth3Kvm78AuhFqlJ69THW8siKuzE7f0DmXNE8N46YaunMwt5M6P/uDW97eyNVnVKpdS8vHGY9z96XYiAr348YFBxIaYoJtVeUMJBj1a4xDfZq68dWsPjmcX8M+fDjRqupIyAy8sPUh4C09mDIxo1L2sgk9b6HSdqh1TWmRta0xK0xT0ylh05UfvFORNc3cX3lxxmFd/S2RSj7a8dnO3q8UcVHfz8P7mFydLETVMvZrC7XI+AxbPho+vhYKzKvtvxi8qNNFSRI+C4jxI32a5Oavg5uLE1L7hrP3LMJ6f1IWUMxe47f2t3PHBVp74di/PLz3IqM6t+XZOf4J8TdA96HwG7PpMxVjXUaGwT2QADwyP5tsd6fyy92SDp1y0NZWkzHyeGR+Du4sZUvUtQdws9R49+FPdY+2IJiroqndmxQq9olBX6tkCrosN4o0p3avPzMs7rRrhOoK7pQLfEGjRXnUxaiglhbD+dXgnTv2BDHkSHtwOXW+yvK80YrDyJZsp2sVY3F2cmdY/gvVPDufZCTEcPp3H4p3pzBnajgV39sLTre40/jopLlD7E9KgokiM4OGR7ekR6sdfl+xtUC+AcxeKeWvlEQZFBzKqsxlT9c1N5DBVRtrBNkebpqA38wdXr8sy1+4d0o57h0Yx77ZrcHGu4cdSUQeivH+owxA1TO0N1NdNISUkLoP3+qqkkegR8OA2GPF36/lJPXwgtJ/VBb0CD1dn7h4Uyfonh/PLw4N4elyn6p/86suZI/DhSEhcqlLn/cONuszV2Yl5t/WgzCB57Ovd9U6QemvlYfKLSnl2Qoxttio0FicntUpP26LCaR2EpinoQpRXXbwk6IPaB/LXcZ1xrUnMQblbPHyhTXcLGGlB2g2HkoL6uSmyDsHnN8JXt6uEnrt+UBEM/hHmstJ4okeqHqQ2lOLt6eZCl7Ym8JeDasn3/jD1/U1dXKvvvDrCW3jx3MQu/HEsm/fXJxt93aFTeXy+NZWpfcPoGORdP5ttkR5TwdldVWF0EJqmoEN5o4u0+l1zbF35I72d+g1rImIQCGfj/OiFufDb3+C/AyB9B4x9FeZsvNTazhaoDF9cbV07TE1JISx9XLXka90F5mxocDOIm3uFMD62DW8sP8S+9Nw6x0speX7pAbw9XHlsVIcGzWlzeAZA1xthz9dQlG9ta0xC0xV039D6CXr2MTXekfznFXj4QnCv2v3oBgPs/EyFIW59D665Ex7eCf3mqCYTtkRQLDRvbTNuF5OQfQw+HqN8vgMeUpvNjWjTJoTgpRu60tLbnUe+2kVBce2hjCsOnmZT0lkeG9W+4ZmstkjcLLWJvu8ba1tiEpquoPuFQmEOFFVfw+UqKvznjpBQVB1RwyBjp0qcupLj2+DDEfDTg2oj6Z61KnXfqwHFoyyBEKpG+tHVKh7e3kn4Gf5vqMoGve0LVZLWBB+ifp5uvHFLd46dvcALSxNqHFdUWsZLyxJo36o5U/sZ56u3G0J6Q+tY2O4YzS+asKCXlyWtWtKzNo6tg+ZBqu62I9JuuIqWSNlw6dj5k7DkXtUSLu8U3PghzPod2vawmplGEz1SJY5l7LK2JQ2ntFi5t76+E1pEqXj+TqZtqDKgXSD3DmnHl9vS+P1A9XsOn2xKIfVsAc9OiKl9j8keEQJ6z4LT+2yzJWM9cbDfTj3wrRB0I9wuUqoVeuQQ66Usm5vgOBX5k7xWJVtsfAvejYMDS1RI3IPx0G2K/Xz/7UYAwqJZoyYl5zgsvA62zoc+96oPUjNtOD8+ugOxwb48vXgvp88XXnYuM6+Qd1cnMapzK4Z0sFKXMXMTe4tqNO4AIYxNV9CvSC6qlcwEuJDluO4WUFUPIwZBwlJ4rx+sfE59gD3wB4z6h/0VIvMMUPsC9uhHP7wc/m8wZCaqErPXvWbWRDY3Fyfevq0HhSUGnvh2z2UlgV///RBFpWX8fXyM2ea3Ou7NofutquFKgXHdy2yVpivoXq1UqVFjVujH1qlXR4s/v5J2IyD/lIp4uXMx3P6l8pnbK9Gj4ES8/fyRlpWqD9IvpoBPiGr80eUGi0zdrmVz5l4fw4YjZ/h4kypatS89l293pDNzYCSRgVauv2Ju4mapXsO7F1nbkkbRdAXdyUlFuhizQk9eB/6Rl/zujkrcTLjjG7hv86XQP3smepTaF7BGV6b6cv4k/G+icnX1nAazV6iGDBbktt6hjIlpzWu/HeJARi7//PkALbzceHBEtEXtsAqtu6iEtPhPVESXndJ0BR0ua3RRI2WlKovSkd0tFbi4Q4drTdd0wtoE91TdqWzdj350DSwYpDZwJy+Aie+ohh0WRgjBv27qhp+nK1M//IP41HM8MaYjPh42FpZqLnrfDdlHLz2Rm4viguqjyUxA0xZ039C6o1xO7oai844Zf+7oODmr6J2klbYZkmYog7X/gs9uAM8W8KfV0ON2q5oU4KVCGXMKSujS1ocpcaFWtceixExSvwdzbI6WFsOh31Txun9Hw9b/mn4OwAQVguwYv3BVOrbkYs0roqbiP3dUokfBge9VvY6grta25hL5WbDkT6qna7dbYfybNrPxPLh9Sz6d1YfoVs2rL1LnqLi4q4S5ze+qKpY+bRt3P0MZpGyE/Yvh4I8q78XDD2JvhvajTWHxVTRxQa+IdEmHwPbVj0leB6272m4SjaZ22o1Ur0krbUfQUzfDd7PUZu31/ylvuGBbwjnUUUMU66LXTNg0D3b+T7XMqy9SqkbY+75TC4n8UyocuNN4JeRRw83q0mzagl610UV1gl5SCMf/UDvgGvvEp436QE5aWe8iVibHYIDN82DVC6o64uyV0KabdW3SXE5ApFoE7PgUBj+hWtYZw+kDaiW+f7HK6HV2g/ZjVAnpDmPBzdOsZlfQtAW9rlj09G1QWqj95/ZO9EjY8p4q8+BupSqBhbmw+E9w5HeImaw2Pj18rGOLpnZ63w1f3QGHf4POE2oel52sBHzfYshKUOG+UUNVP4DOE1SNJAvTtAXdu636JdS0MZq8Tp0PH2BZuzSmJXqUeow+tkG1HrMGK/6hnhLG/Vt1vLcxF4umCu2vBZ9gtTl6paCfP6lcKfu/U64VUOGO172uPqibW9dV1bQF3dlF/eJqCl08tq489E2vpOya0H7Kj3l0lXUEPeuQ8sn2ng1977H8/Jr64ewCvWbAmpfg7FHVEOfgj2o1nrIRkBDUDUY/D11uvPSkbwM0bUGHqxpdVFJ4Hk7shEGPWd4mjWlxcVOPwkdWqE0rS6+OV/xDdXAa+qRl59U0nJ7TVEjp5zeqoAlDKbSIhqFPqc3NmoIorIwWdL8w9Sh+JambQZY1jYSipkD0SDi0TPk9LZmBmbIRDv8KI+fqSCl7wjsIet6lSjD3u1+JeFA3m3eVaUH3DYW8DCgrubzG9LH1qrVaSB/r2aYxHVXDFy0l6AYDLH9WufX63W+ZOTWm4/p51rag3jTtTFFQLhdpgPMnLj9+bB2E9gVXD+vYpTEtAZHqkdmS1RcPLFFNQ0Y8Y5VUfk3TQwt6ZSx6FT96fhac3q/dLY5G9CjlXisprHtsYyktglX/VDHw3W41/3waDVrQq3QuqhLpUtG1R8efOxbtRkLpRUjbbP65tn+o3lOjn3e8puIam0ULekWj3aqRLsfWgbsPtOlhFZM0ZiJiIDi7m7/64sVzsO41VV8+eqR559JoqqAF3cVd9Qqt6nJJXqe69xib9quxD9y8VJKYuf3oG95QmaGjnzfvPBrNFWhBB+V2yS13ueSkwbljurqioxI9CrISjW8OXl/OpcIf/wfdb4egWPPModHUgBZ0uLzRxbH16lX7zx2Tik5MR83kdln9IggnFdmi0VgYLehQ3oruhIobTl4HXi2hVWdrW6UxBy07qn6d5vCjZ+yCfd+omHPfYNPfX6OpAy3ooFwuhhLIO6lW6JFDbD4jTNNAhFAblclrVTKZqZBSJRF5trB+mV5Nk0ULOlwKXTy6ShWk1+4WxyZ6lGormB5vunseWaHCXYc+ZZWyqRoNGCnoQoixQohDQogkIcRVbTyEEFOFEHvL/20WQnQ3valmpCK5aNfn6lVviDo2UUNVWWRTRbuUlcKKuRAQpTreaDRWok5BF0I4A/OBcUAMcLsQIuaKYceAoVLKbsALwPumNtSsVJS/PP6HWq0HRFrXHo158fBVZR1MJei7F6kGByP/Ydb2YhpNXRizQu8DJEkpk6WUxcBXwKSqA6SUm6WU58q/3AqEmNZMM+PmpXyfoN0tTYXokXBytyrz0BiKL8Cal1URt5hJdY/XaMyIMYIeDFQN2k0vP1YTdwO/VndCCHGPECJeCBGfldXIPyRTU+F20YLeNKgMX1zduPtsma/2Xca8qDfSNVbHGEGv7l0qqx0oxHCUoD9V3Xkp5ftSyjgpZVzLljbWVbzC7aL9502DoG7gGdg4t0t+pmpt1/l6COtrOts0mgZiTG57OlC1x1IIkHHlICFEN+BDYJyU8qxpzLMgMZOhWQB4t7a2JRpL4OSk3C5JK1X+gVMDAr7WvqKaiI98zuTmaTQNwZh38XagvRAiUgjhBtwG/FR1gBAiDFgC3CWlPGx6My1A7M0w8T/WtkJjSaJHQcFZ5UuvL1mHYcenKqolMNrkpmk0DaFOQZdSlgIPAr8DCcA3UsoDQog5Qog55cPmAi2A94QQu4UQJgzw1WjMRLsRgGhY1ujK58DVU8WdazQ2glHlBKWUy4BlVxxbUOX/s4HZpjVNozEzXoHQtodKKBv6F+OvS90Mh35R9Vqa29hekKZJozNFNU2b6FFwfBtczDFuvJSw/Bnwbgv9HjCraRpNfdGCrmnaRI8CWaaamhjDge/hxA4Y8Xdw8zSvbRpNPdGCrmnaBMeBu69x4YsVfUJbdVH1zjUaG0O35NE0bZxdoN0wtTEqZe3JQds/gnMpMHWx7hOqsUn0Cl2jiR4F50+oTkY1cTEH1r8GUcN0n1CNzaIFXaNpVy7QtbldNr6pRH308zrFX2OzaEHXaHyDoVVMzYKekwZbF0C3W6GNfVWG1jQttKBrNKCSjFI3q+qJV7L6JfWq+4RqbBwt6BoNKD96WTGkbLz8+Mk9sPdr6HffpQJuGo2NogVdowEI669S+au6XSr6hDbzh0GPWc82jcZItKBrNACuHhAx+HJBT1qpEo6GPgnN/KxmmkZjLFrQNZoKokdBdrL6ZyhTfUL9IyHubmtbptEYhU4s0mgqqIgvT1oFLh6QeRCmLNR9QjV2gxZ0jaaCFu3UijzhZzhzWJUFiJlsbas0GqPRLheNpirRo5TfPO+k7hOqsTu0oGs0ValoHt1pAoT3t64tGk090S4XjaYq7YarOuf95tQ9VqOxMbSgazRVcXGHsS9b2wqNpkFol4tGo9E4CFrQNRqNxkHQgq7RaDQOghZ0jUajcRC0oGs0Go2DoAVdo9FoHAQt6BqNRuMgaEHXaDQaB0FIKa0zsRBZQKpVJq+ZQOCMtY2oB/Zkrz3ZCvZlrz3ZCvZlry3aGi6lbFndCasJui0ihIiXUsZZ2w5jsSd77clWsC977clWsC977clW0C4XjUajcRi0oGs0Go2DoAX9ct63tgH1xJ7stSdbwb7stSdbwb7stSdbtQ9do9FoHAW9QtdoNBoHQQu6RqPROAha0AEhRKgQYo0QIkEIcUAI8Yi1baoLIYSzEGKXEGKptW2pCyGEnxDiOyFEYvnP2GZ7uwkhHit/D+wXQnwphPCwtk1VEUJ8LITIFELsr3IsQAixQghxpPzV35o2VlCDrf8ufx/sFUJ8L4Tws6KJl1GdvVXOPSGEkEKIQGvYZixa0BWlwJ+llJ2BfsADQogYK9tUF48ACdY2wkjmAb9JKTsB3bFRu4UQwcDDQJyUsivgDNxmXauuYiEw9opjTwOrpJTtgVXlX9sCC7na1hVAVyllN+Aw8FdLG1ULC7naXoQQocBoIM3SBtUXLeiAlPKklHJn+f/zUIITbF2rakYIEQKMBz60ti11IYTwAYYAHwFIKYullDlWNap2XIBmQggXwBPIsLI9lyGlXA9kX3F4EvBp+f8/BSZb0qaaqM5WKeVyKWVp+ZdbgRCLG1YDNfxsAd4CngRsPoJEC/oVCCEigGuAP6xsSm28jXqDGaxshzFEAVnAJ+Uuog+FEF7WNqo6pJQngNdRK7GTQK6Ucrl1rTKK1lLKk6AWJ0ArK9tjLLOAX61tRG0IISYCJ6SUe6xtizFoQa+CEKI5sBh4VEp53tr2VIcQYgKQKaXcYW1bjMQF6An8V0p5DXAB23EJXEa573kSEAm0BbyEEHda1yrHRAjxd5Src5G1bakJIYQn8HdgrrVtMRYt6OUIIVxRYr5ISrnE2vbUwkBgohAiBfgKGCGE+Ny6JtVKOpAupax44vkOJfC2yCjgmJQyS0pZAiwBBljZJmM4LYRoA1D+mmlle2pFCDEdmABMlbadCNMO9eG+p/zvLQTYKYQIsqpVtaAFHRBCCJSPN0FK+aa17akNKeVfpZQhUsoI1Ibdaimlza4ipZSngONCiI7lh0YCB61oUm2kAf2EEJ7l74mR2OgG7hX8BEwv//904Ecr2lIrQoixwFPARCllgbXtqQ0p5T4pZSspZUT531s60LP8PW2TaEFXDATuQq12d5f/u87aRjkQDwGLhBB7gR7Ay9Y1p3rKnyK+A3YC+1B/HzaV+i2E+BLYAnQUQqQLIe4G/gWMFkIcQUVj/MuaNlZQg63vAt7AivK/swVWNbIKNdhrV+jUf41Go3EQ9Apdo9FoHAQt6BqNRuMgaEHXaDQaB0ELukaj0TgIWtA1Go3GQdCCrtFoNA6CFnSNRqNxEP4fQPD9Vdpf3AsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gc2 = grangercausalitytests(x2,maxlag=15,verbose=False)\n",
    "p2 = [gc2[lag][0][test][1] for lag in lags]\n",
    "\n",
    "p_df = pd.DataFrame({'imports -> exports':p1,'exports -> imports':p2},index=lags)\n",
    "\n",
    "p_df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this, we can conclude that there is no Granger causation between these two time series."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Using the two time series from the exercise in the previous lesson - confirm that there is no Granger causation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}