spurious-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spurious Correlation\n",
    "\n",
    "**Coauthored by Samuel (Siyang)  Li, Thomas Sargent, and Natasha Watkins**\n",
    "\n",
    "This notebook illustrates the phenomenon of **spurious correlation** between two uncorrelated but individually highly serially correlated time series \n",
    "\n",
    "The phenomenon surfaces when two conditions occur\n",
    "\n",
    "* the sample size is  small \n",
    "\n",
    "* both series are highly serially correlated\n",
    "\n",
    "We'll proceed by \n",
    "\n",
    "- constructing many simulations of two uncorrelated but individually serially correlated time series \n",
    "\n",
    "- for each simulation, constructing the correlation coefficient between the two series\n",
    "\n",
    "- forming a histogram of the correlation coefficient\n",
    "\n",
    "- taking that histogram as a good approximation of the population distribution of the correlation coefficient\n",
    "\n",
    "In more detail, we construct two time series governed by\n",
    "\n",
    "$$ \\eqalign{ y_{t+1} & = \\rho y_t + \\sigma \\epsilon_{t+1} \\cr\n",
    "             x_{t+1} & = \\rho x_t + \\sigma \\eta_{t+1}, \\quad t=0, \\ldots , T } $$\n",
    "             \n",
    "where\n",
    "\n",
    "* $y_0 = 0, x_0 = 0$\n",
    "\n",
    "* $\\{\\epsilon_{t+1}\\}$ is an i.i.d. process where $\\epsilon_{t+1}$ follows a normal distribution with mean zero and variance $1$\n",
    "\n",
    "* $\\{\\eta_{t+1}\\}$ is an i.i.d. process where $\\eta_{t+1}$ follows a normal distribution with mean zero and variance $1$\n",
    "\n",
    "We construct the sample correlation coefficient between the time series $y_t$ and $x_t$ of length $T$\n",
    "\n",
    "The population value of correlation coefficient is zero\n",
    "\n",
    "We want to study the distribution of the sample correlation coefficient as a function of $\\rho$ and $T$ when\n",
    "$\\sigma > 0$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll begin by importing some useful modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Empirical distribution of correlation coefficient r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now set up a function to generate a panel of simulations of two identical independent AR(1) time series\n",
    "\n",
    "We set the function up so that all arguments are keyword arguments with associated default values\n",
    "\n",
    "- location is the common mathematical expectation of the innovations in the two independent autoregressions\n",
    "\n",
    "- sigma is the common standard deviation of the indepedent innovations in the two autoregressions \n",
    "\n",
    "- rho is the common autoregression coefficient of the two  AR(1) processes\n",
    "\n",
    "- sample_size_series is the length of each of the two time series\n",
    "\n",
    "- simulation is the number of simulations used to generate an empirical distribution of the correlation of the two uncorrelated time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spurious_reg(rho=0, sigma=10, location=0, sample_size_series=300, simulation=5000):\n",
    "    \"\"\"\n",
    "    Generate two independent AR(1) time series with parameters: rho, sigma, location, \n",
    "    sample_size_series(r.v. in one series), simulation. \n",
    "    Output : displays distribution of empirical correlation\n",
    "    \"\"\"\n",
    "        \n",
    "    def generate_time_series():\n",
    "        # Generates a time series given parameters\n",
    "        \n",
    "        x = []  # Array for time series\n",
    "        x.append(np.random.normal(location/(1 - rho), sigma/np.sqrt(1 - rho**2), 1))  # Initial condition\n",
    "        x_temp = x[0]\n",
    "        epsilon = np.random.normal(location, sigma, sample_size_series)  # Random draw\n",
    "        T = range(sample_size_series - 1)\n",
    "        for t in T:\n",
    "            x_temp = x_temp * rho + epsilon[t]  # Find next step in time series\n",
    "            x.append(x_temp)\n",
    "        return x\n",
    "    \n",
    "    r_list = []  # Create list to store correlation coefficients\n",
    "    \n",
    "    for round in range(simulation):    \n",
    "        y = generate_time_series()\n",
    "        x = generate_time_series()\n",
    "        r = stats.pearsonr(y, x)[0]  # Find correlation coefficient\n",
    "        r_list.append(r)    \n",
    "     \n",
    "    fig, ax = plt.subplots()\n",
    "    sns.distplot(r_list, kde=True, rug=False, hist=True, ax=ax)  # Plot distribution of r\n",
    "    ax.set_xlim(-1, 1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparisons of two value of $\\rho$\n",
    "\n",
    "The next two cells we'll compare outcomes with a low $\\rho$ versus a high $\\rho$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHrRJREFUeJzt3XlwJOd5HvDnnXsGM1gAC2Dv5XJ5SpYiUkRRjliRopUsU7JDMjGTUCkltKJkIyV2yZWkHKmUylWVxE6qYkeJLGejyzFFUiIplmlHskKZZBhR4lJYnntwyT2x2AMYDI6ZwWAuzJs/uhucxeKYGXRPH3h+VagdTPd0v9sYPPjm66+/FlUFEREFS8jtAoiIyH4MdyKiAGK4ExEFEMOdiCiAGO5ERAHEcCciCiCGOxFRADHciYgCiOFORBRAEac2PDg4qPv27XNq80REgXTkyJEpVR3a6HYcC/d9+/ZhdHTUqc0TEQWSiJy3YzvsliEiCiCGOxFRADHciYgCiOFORBRADHciogBiuBMRBRDDnYgogBjuREQB1HK4i8gtIvJq01deRH7LyeKIiKgzLV+hqqonAdwGACISBnARwJMO1UXkmIcPjwEA/s4H9rpcCZFzOu2W+SiA06pqy2WyRERkr07D/QEAj9hZCFE3WK12oqBrO9xFJAbgHgCPrbDsoIiMishoNpu1oz4iIupAJy33TwB4WVUnli9Q1UOqOqKqI0NDG56xkoiIOtRJuH8K7JIhIvK0tsJdRFIAfgnA950ph4iI7NDWzTpUtQRgq0O1EBGRTRy7ExORV12aXUAmwbc+BRvf4bSpHL80h4cOj2EwHUMiGsaDH9zndklEjuDcMrSpvHA6BwCYKlZx4nLe5WqInMNwp01jqljBual5HLh1GMloGCevFNwuicgx7JahTePI+RkogJuG05gqVnA6W3S7JCLHsOVOm8bLYzMIi2BnXxJ7+lPIl+uYzJfdLovIEQx32jReuzCLHX0JRMMh7OpLAgDeuDjnclVEzmC406agqjh5pYDtvQkAwI4+49+jF3lSlYKJ4U6bwmShgplSDdu3GKEej4TRn4riFPvdKaAY7rQpWCNjrJY7AAxnEjg1yXCnYGK406ZghfhwU7gPZeI4nS1isaFulUXkGIY7bQqns0X0JiLoiYWXnhvOxFGtNzA+U3KxMiJnMNxpUziTnccNw2mIyNJzw5k4AODtCXbNUPAw3GlTODNVxP7B9FXPDWWMLhqeVKUgYrhT4BXKNUzkK7hhuOeq55OxMLb1xtlyp0BiuFPgnZ2aB4BrWu7Wc2enGO4UPAx3CjxrDpkbhnquWbZvMIXzOZ5QpeBhuFPgncnOIyTA3q2pa5bt29qD3HwV+XLNhcqInMNwp8A7lyuhLxXDE0cuXrNs36DRmj9ndt0QBQXDnQJvbLqEgVRsxWXXm+F+luFOAdNWuItIn4g8LiJvisgJEfnLThVGZJex3DwGelYO970DRlfNuSn2u1OwtNty/68A/lxVbwXwPgAn7C+JyD75cg0zpdqq4Z6IhrFzSwLPnZzscmVEzmo53EWkF8CHAHwDAFS1qqqzThVGZIcxcyTMauEOGP3uU8VKt0oi6op2Wu77AWQBfEtEXhGRr4vItWPLiDzkO4fHAKwf7rn5ardKIuqKdsI9AuD9AL6mqrcDmAfwxeYVROSgiIyKyGg2m7WxTKLOTJuhvXWNcL9+aw9K1UXMlTgckoKjnXAfBzCuqofN7x+HEfZLVPWQqo6o6sjQ0JBdNRJ1bHq+gp5YGPFoeMXlDx8eW7qI6WyOI2YoOFoOd1W9AuCCiNxiPvVRAMcdqYrIJtPz1TW7ZABga9pYfp7hTgESaXP93wTwHRGJATgD4DP2l0Rkn9lSDTvNm2GvZqAnBgHHulOwtBXuqvoqgBGHaiGylapibqGGd+3oXXO9aDiETCKC8ZmFLlVG5DxeoUqBNVuqod5QbElG1123LxXDRYY7BQjDnQLr8lwZAFoK9/5UFOOzvEqVgoPhToF1ec5oibfacr88W+bNsikwGO4UWO213GOoNxQT+bLTZRF1BcOdAuvy3AJCAqQT648b6EsZfwAuzrLfnYKB4U6BdXmujN5EFCGRddftN6cEHp9hvzsFA8OdAuvybBm9LXTJAE0td46YoYBguFNgXcmXW+pvB4yx7oPpOMe6U2Aw3CmQVBWX5xZaDncA2NWfZJ87BQbDnQIpX66jXGugt4WTqZbd/Um23CkwGO4USNbNN1oZKWPZ3We03Bsc604BwHCnQMoWzHCPt94ts7s/iWq9wbsyUSAw3CmQrHDPtNFy39VvzB45zn53CgCGOwXSUrdMvPVwf/XCHACw350CgeFOgZQtVBASIBlb+Q5MK9mSMLpwJuY4BQH5H8OdAmmqWEE6Hmnp6lRLIhpCNCy4wvllKAAY7hRI2UKlrZEyACAi2JKMMtwpEBjuFEjZYgWZNkbKWHoTUVxhtwwFAMOdAmmqUG275Q4AvUmGOwUDw50C56EXz2OyUG5rpIylNxHFZKHMC5nI99p694vIOQAFAIsA6qrKm2WT5yxUF9HQ9sa4W3qTEdQWFdOlKgbTcQeqI+qO9t/9wEdUdcr2SohsUqzUAbQ3xt3Saw6HvDJXZriTr7FbhgKnUDbDvYOWuzWLJPvdye/aDXcF8H9E5IiIHHSiIKKNKlZqANDZaBkr3Dkcknyu3abNXap6SUSGATwtIm+q6vPWQjPwDwLA3r17bSyTqHVFs+XeSZ97Oh6BgC138r+2Wu6qesn8dxLAkwDuXLb8kKqOqOrI0NCQfVUStaFQqSMSEsQj7fc6hkOCdDyyNPEYkV+1/O4XkR4RyViPAXwcwFGnCiPqVLFcRzoRgbQx9UCzdCKC3DzDnfytnc+t2wA8af7CRAA8rKp/7khVRBtQrNQ7Gilj6YlHMFWs2lgRUfe1/BugqmcAvM/BWohsUSjX0Zdq/2SqJR2P8IYd5HscCkmBUyjXlsardyIdjyDHljv5HMOdAqVab2C+uohMcmPdMgu1RXzrhbM2VkbUXQx3CpSs2Z2y0ZY7AMxXFm2picgNDHcKlAnz4qPeDsa4W9Jx4+5N1jQGRH7EcKdAmcxbN8buvOXeY7bcrYuhiPyI4U6BMlkwWu6dXJ1qeadbhuFO/sVwp0CZyJcRknda351YarlXGe7kXwx3CpSJfAWZRLStG2MvFw2HkIiG2C1DvsZwp0CZyJc31CVj6YlFeEKVfI3hToGSLVQ2dDLVko4z3MnfGO4UKBP58oaGQVp64hGU2OdOPsZwp8Co1BcxU6rZ0nJPxcIo8SIm8jGGOwWGNcbdrpb7fLUOVd3wtojcwHCnwJg0b7Bh3SpvI1KxMBrKq1TJvxjuFBiT+Y1fwGRJxYxtzMzXNrwtIjcw3CkwJpbCfeMt956YMb/MTIlT/5I/MdwpMCYKFUTDgpQZzBthbWOa4U4+xXCnwJjIlzGcSWzo6lRLKm51yzDcyZ8Y7hQY2UIFQ5m4LdvqsfrcS+xzJ39iuFMgPHx4DCevFLCt155wj0dDELDlTv7VdriLSFhEXhGRP3OiIKJOFcp1bOtN2LKtkBh99zyhSn7VScv9CwBO2F0I0UbUFhtYqC3aFu6AMRyS4U5+1Va4i8huAL8C4OvOlEPUmYI5Pe+wTX3uAJCKhzHNbhnyqXZb7r8P4LcBNByohahjhbJx4vPYpbxt2+yJRTDLE6rkUy2Hu4j8KoBJVT2yxjoHRWRUREaz2awtBRK1wmq523F1qiUVY8ud/KudlvtdAO4RkXMAHgVwQEQeal5BVQ+p6oiqjgwNDdlYJtHarDlg0hu4vd5yKbPlzsnDyI9aDndV/ZKq7lbVfQAeAPCMqn7ascqI2mDdzNqaE8YOqVgY1cUG5quc+pf8h+PcKRCKlTpSsTDCoY1fnWrpiZvzy7Brhnyoo3BX1edU9VftLoaoU8VK3dYuGaBpZkgOhyQfYsudAmG+UkeP7eFuTh7Gljv5EMOdAqFYWbS95W7NL8PhkORHDHcKhGKlxpY7UROGO/letd5AudawveWeiIWNycPY504+xHAn37Na1naHe0gESU4eRj7FcCffmyoaN8ZOxzd+B6blUrEI76NKvsRwJ9+zwt3uPnfAuJcqW+7kRwx38r1c0ZluGcC43R5PqJIfMdzJ93LzVreMA+HOljv5FMOdfC9XrCISEsQi9r+djW4ZTh5G/sNwJ9/LFitIxyMQsW9eGUsqFkG13kCJk4eRzzDcyfdyxSrSNs7j3sy6kIldM+Q3DHfyvdx8ZWmqALtZI3A4HJL8huFOvpcrVh05mQqw5U7+xXAnX1NV5IpVR8a4A5z2l/yL4U6+li/XUV1sON7nzrHu5DcMd/K1nINTDwBAMhaGCDDDaX/JZxju5Gs5s0XtVLdMSAR9yShvtUe+w3AnX5sqOHd1qqU/FWOfO/kOw518bcqh6X6b9fcw3Ml/GO7ka1afe8qhce4A0J+KYprj3MlnWg53EUmIyEsi8pqIHBORf+tkYUStyBWrSMXCCIfsn3rA0p+KYZYtd/KZdpo7FQAHVLUoIlEAPxGRH6rqiw7VRrSuqWLFsZOploGeGIdCku+03HJXQ9H8Nmp+cao8cpWTV6da+lIxVOoNLHDyMPKRtvrcRSQsIq8CmATwtKoeXrb8oIiMishoNpu1s06iFU3NVxwP94GeKABgml0z5CNthbuqLqrqbQB2A7hTRN6zbPkhVR1R1ZGhoSE76yRakZNTD1j6UjEA4Fh38pWORsuo6iyA5wDcbWs1RG2o1huYW6g5dnWqZaDHDHe23MlH2hktMyQifebjJICPAXjTqcKI1jO9NMY96uh++s2WO0+qkp+083l2B4A/EpEwjD8K31PVP3OmLKL1TTk8r4ylP2X88Zjl/DLkIy2Hu6q+DuB2B2shaovT88pYtiSjEGHLnfyFV6iSb3VjXhkAiIRD2JKMss+dfIXhTr6Vm+9OuAPW5GHsliH/YLiTb+WKVcQjIcQizr+N+1Oc9pf8xfkmD5FDpopVDKbjEHFuXhkAePjwGErVRVTqDUf3Q2QnttzJt6aKFWxNx7qyr1QswpY7+QrDnXwrN1/BYDrelX31xMKcfoB8heFOvpUrVrG1p1st9zDKNU4eRv7BcCdfajTU7JbpUsvdHJFjXThF5HUMd/KlmVIVtUXF9t7uhHtv0rhKdbJQ7sr+iDaK4U6+NGlewDTcm+jK/jIJo+X+xJGLXdkf0UYx3MmXJvJGC3pbl1rumYTRci+UeSET+QPDnXxpMm+23DPdabmnYmGEBCiU613ZH9FGMdzJl6yW+7NvTnZlfyERpOMRhjv5BsOdfGmiUEYyGkYk3L23cCYRRaHCbhnyB4Y7+dJkvoLeZHdnz8gk2HIn/2C4ky9NFCroTTh7B6blMoko8gx38gmGO/nSZL68NIKlWzKJCEqVOuqLnECMvI/hTr7TaCgmC5WlsefdkklEoDBmoyTyOoY7+U5uvorFhqK3y+FudQNZI3WIvIzhTr5zZc4IV2tKgG6xPilYV8cSeVnL4S4ie0TkWRE5ISLHROQLThZGtJoLMyUAxq3vusnq4+f8MuQH7XyurQP4Z6r6sohkABwRkadV9bhDtRGt6MK0O+GejkcgeOfqWCIva7nlrqqXVfVl83EBwAkAu5wqjGg1F2ZK6E1EkIyFu7rfcEiQikfYLUO+0FGfu4jsA3A7gMPLnj8oIqMiMprNZjdeHdEyDx8ew0tnp7FnIOXK/nsTEWTZLUM+0Ha4i0gawBMAfktV883LVPWQqo6o6sjQ0JBdNRJdZXq+hj397oR7JsGWO/lDW+EuIlEYwf4dVf2+MyURrU5VMVuqYs9A0pX9ZxJRDoUkX2hntIwA+AaAE6r6X5wriWh1hUod9Ya61i2zJRnFZKGCGq9SJY9rp+V+F4C/C+CAiLxqfn3SobqIVjQzb1wd6la3TF8yCtV3xtoTeVXLQyFV9ScAxMFaiNY1UzLDfSCJyy4EbJ85/PLi7IJrnx6IWsErVMlXpueN+dR3u9VyTxkXMl2cWXBl/0StYriTr8yUqsjEI0hEuzvG3bLFnPLg0izDnbyN4U6+MjNfRX9PDA8fHnNl/9FwCOl4BBcZ7uRxDHfylZlSFf2p7k4YtlxfKspwJ89juJNv1BYbmFuoob+nu3PKLLclyXAn72O4k29cml1AQ4GBLk8Ytlx/KoZLswtQVVfrIFoLw51843zOmA1yazruah1bklGUaw1Mz/OOTORdDHfyjfPmVL8DLnfLbDX3f878Y0PkRQx38o2x3DwiIen6vVOXG8wYnxzOZIuu1kG0FoY7+ca5XAkDPTGExN0LpftTMUTDgtPZeVfrIFoLw518Y8wMd7eFQ4K9Aym23MnTGO7kC6qKsenSUn+32/YPpXFmii138i6GO/lCtlDBQm3REy13ANg/1IPzuXnUOfUveRTDnXzhnZEy7g6DtNwwmEZtUfG15067XQrRihju5AtLY9w91HIHjE8URF7EcCdfGMvNIyRAX4+788pYbt6eAQBc4S33yKMY7uQL56dL2NmXRCTkjbdsbyKK67amOMcMeZY3flOI1nE+V8J1W71z56OHD48hk4hyXnfyLIY7+cL53Dz2DvS4XcZVdvUlMVOqYa5Uc7sUomsw3Mnz8uUaZkq1pZtje8XOvgQA4OilOZcrIbpWy+EuIt8UkUkROepkQUTLjeW8MWHYcju3JAEARy8y3Ml72mm5fxvA3Q7VQbSq8x4N9554BH3JKI5eyrtdCtE1Wg53VX0ewLSDtRCt6FzOuMzfK2Pcm+3sS+IYW+7kQexzJ887PVnE9t4E4tGw26VcY2dfEmem5lEo86QqeYut4S4iB0VkVERGs9msnZumTezw2WmkXZ7DfTW7+qx+d3bNkLfYGu6qekhVR1R1ZGhoyM5N0ybVaCiyhQqGMt6YU2a5Pf1GuL96YdblSoiuxm4Z8rTL+TKqiw0MezTcU/EIrh/swStjM26XQnSVdoZCPgLgZwBuEZFxEfmsc2URGU5NGjfE8GrLHQBu39OHVy7MQlXdLoVoSTujZT6lqjtUNaqqu1X1G04WRgS8E+7DmYTLlazu9r19yBYqGJ/hVATkHeyWIU87NVlEKhZGT8x7I2Ust+/tBwC8wn538hCGO3na8ct5bOtNQFy+KfZabtmeQSIaYr87eQrDnTyrWm/gxKU8dpvDDb3qsdFxbO9N4JUxttzJOxju5FlvTRRQXWxgV7+3wx0A9gykcPxSHpX6otulEAFguJOHvT5uXNa/u98787ivZk9/CtXFBo5xnhnyCIY7edYbF2fRl4qiP+WNW+utZe+A8QeIXTPkFQx38qzXLszhvbu2ePpkqqU3GUVfMoqfn+XceuQNDHfypJn5Kk5cyeOO6/rdLqVlNw6n8cLpKdQXG26XQsRwJ2/6yakpqAIfutk/cxTdtC2DQrmO18bZNUPuY7iTJz3/VhbJaBjHfXSC8sahNEIC/N+TnBGV3MdwJ89RVTz/dhY3DKcR8kF/uyUZC+O2PX147i2GO7mP4U6e88bFOUzkK7h5OO12KW27+z3b8fr4HE5NFtwuhTY5hjt5zmOj44iEBO/ZtcXtUtqmCoQEeOzIuNul0CbHcCdPKdcW8SevXsQv7OxFwoO31VtPJhHFLdsy+P7LF1HjqBlyEcOdPOVPX7uEfLmOO64bcLuUjt15/VZkCxU8+fJFt0uhTYzhTp5RrTfwlWfexi/s7MUNQz1ul9Oxm7elsasvif/27NtsvZNrGO7kGY+8NIYL0wsYua7fF1elrkZEcODWYVyYXsAjL425XQ5tUgx38oTT2SJ+54dv4sahNG7elnG7nA27dXsG+4d68B9+cAK5YsXtcmgTYriT6yYLZfyjPz6CRDSE++/Y7etWu0VE8Nf+0k5U6w385x+ddLsc2oQY7uSqN6/k8Ynf/38Yy5Xwa3fsRm/S+zNAtmpbbwIfvGEQ3x29gBdOTbldDm0ybYW7iNwtIidF5JSIfNGpoij4VBXf+/kF3PfVF1CtN/CZu/Zh/6D/Llpaz8fetQ03DqXxhUdfweU53kCbuqflcBeRMICvAvgEgHcD+JSIvNupwii45ko1/I0/+Cl++4nX8f69/fiNAzfiuq3+HR2zllgkhE++dwfKtQbu/9rPcOzSnNsl0SYRaWPdOwGcUtUzACAijwK4F8BxJwqj4Gg0FOX6Is7nSvjx8Ql8+6fnMFOq4sCtwzhw67Cv5o/pxLbeBB7+hx/AZ771c/zKV36CO67rx3VbU+hLxrAlGcWOvgT2DqSwdyCFbb0JhEPBPh7UHe2E+y4AF5q+HwfwAXvLoSA5eaWA+776AhZqV99X9K/cNIj37e7DTo/f+NpORy/m8fkP34CfnsnhTLaIU5NFLNQWUa1fPQ5eBIiGQuiJh/HKv/q4S9VSELQT7is1J/SqFUQOAjhoflsRkaOdFtZFgwD8cLYrMHWeB/BQd2pZi+ePp/xrAD6o0+SHOv1QIwDcYsdG2gn3cQB7mr7fDeBS8wqqegjAIQAQkVFVHdlwhQ5jnfZinfZinfbxQ42AUacd22lntMzPAdwkIteLSAzAAwCesqMIIiKyV8std1Wti8hvAPgRgDCAb6rqMccqIyKijrXTLQNV/QGAH7S4+qH2y3EF67QX67QX67SPH2oEbKpTVHX9tYiIyFc4/QARUQBtKNxF5G+KyDERaYjIqmehV5u2wDw5e1hE3haR75onam0nIgMi8rS5n6dFpH+FdT4iIq82fZVF5D5z2bdF5GzTstvcqtNcb7GplqeanvfS8bxNRH5mvj9eF5G/3bTMseO53hQZIhI3j80p81jta1r2JfP5kyLyy3bV1GGd/1REjpvH7i9E5LqmZSv+/F2q89dFJNtUzz9oWvag+R55W0QedLnO32uq8S0RmW1a1pXjKSLfFJFJWWWIuBi+Yv4fXheR9zcta/9YqmrHXwDeBWNM5nMARlZZJwzgNID9AGIAXgPwbnPZ9wA8YD7+QwCf30g9a9T5nwB80Xz8RQC/u876AwCmAaTM778N4H4nauukTgDFVZ73zPEEcDOAm8zHOwFcBtDn5PFc673WtM4/BvCH5uMHAHzXfPxuc/04gOvN7YQdOn6t1PmRpvff56061/r5u1TnrwP47yu8dgDAGfPffvNxv1t1Llv/N2EMCOn28fwQgPcDOLrK8k8C+CGMa4p+EcDhjRzLDbXcVfWEqq43n+nStAWqWgXwKIB7RUQAHADwuLneHwG4byP1rOFec/ut7ud+AD9U1ZJD9aym3TqXeO14qupbqvq2+fgSgEkAQw7VY1nxvbZsnebaHwfwUfPY3QvgUVWtqOpZAKfM7blSp6o+2/T+exHGdSXd1srxXM0vA3haVadVdQbA0wDu9kidnwLwiEO1rEpVn4fRaFzNvQD+lxpeBNAnIjvQ4bHsRp/7StMW7AKwFcCsqtaXPe+Ebap6GQDMf4fXWf8BXPvD//fmR6XfE5G4E0Wi9ToTIjIqIi9aXUfw8PEUkTthtKhONz3txPFc7b224jrmsZqDcexaea1d2t3XZ2G06Cwr/fyd0Gqdv2b+LB8XEetCR08eT7N763oAzzQ93a3juZ7V/h8dHct1h0KKyI8BbF9h0ZdV9U/Wez1Wn7Zg3ekM2rFWnW1uZweA98IYz2/5EoArMALqEIB/AeDfuVjnXlW9JCL7ATwjIm8AyK+wnleO5x8DeFBVrYlUbDuey3e3wnPLj0FX3o/raHlfIvJpACMAPtz09DU/f1U9vdLru1DnnwJ4RFUrIvI5GJ+KDrT4Wru0s68HADyuqs0THnXreK7H1vfmuuGuqh9roai1rDZtwRSMjx0RswV1zXQG7VirThGZEJEdqnrZDJvJNTb1twA8qaq1pm1fNh9WRORbAP65m3Wa3RxQ1TMi8hyA2wE8AY8dTxHpBfC/AfxL82OmtW3bjucy606R0bTOuIhEAGyB8VG5ldfapaV9icjHYPwx/bCqLt2rb5WfvxNh1MqUI7mmb/8ngN9teu1fXfba52yv8J19tfqzewDAP2l+oovHcz2r/T86Opbd6JZZcdoCNc4UPAujfxsAHgTQyieBTjxlbr+V/VzTH2cGmNWvfR8ApyZEW7dOEem3ujFEZBDAXQCOe+14mj/rJ2H0IT62bJlTx7OVKTKaa78fwDPmsXsKwANijKa5HsBNAF6yqa626xSR2wH8DwD3qOpk0/Mr/vxdrHNH07f3ADhhPv4RgI+b9fYD+Diu/jTc1TrNWm+BcULyZ03PdfN4rucpAH/PHDXziwDmzIZQZ8dyg2d//zqMvyoVABMAfmQ+vxPAD5adBX4Lxl/DLzc9vx/GL9ApAI8BiG+knjXq3ArgLwC8bf47YD4/AuDrTevtA3ARQGjZ658B8AaMEHoIQNqtOgF80KzlNfPfz3rxeAL4NIAagFebvm5z+niu9F6D0eVzj/k4YR6bU+ax2t/02i+brzsJ4BNOHLs26vyx+TtlHbun1vv5u1TnfwRwzKznWQC3Nr3275vH+RSAz7hZp/n9vwHwO8te17XjCaPReNn8vRiHcS7lcwA+Zy4XGDdEOm3WMtL02raPJa9QJSIKIF6hSkQUQAx3IqIAYrgTEQUQw52IKIAY7kREAcRwJyIKIIY7EVEAMdyJiALo/wPq+GmH0tGP9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "spurious_reg(0, 10, 0, 300, 5000)  # rho = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For rho = 0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXZyb7QvYFsgMh7LIEAoi4oaL+xF1xK1qXalevt+21P1sfrb331t57f/W2VavUWncRbalLQasILuyr7JAQsgLZCVnIOt/fHxlsjIFMyCRnls/z8eDBzJkzZ945mXzmO9/zPd8jxhiUUkr5F5vVAZRSSg09Lf5KKeWHtPgrpZQf0uKvlFJ+SIu/Ukr5IS3+Sinlh7T4K6WUH9Lir5RSfkiLv1JK+aEAq144Pj7eZGZmWvXySinllbZu3VptjEkY6HYsK/6ZmZls2bLFqpdXSimvJCLF7tiOdvsopZQf0uKvlFJ+SIu/Ukr5IS3+Sinlh7T4K6WUH9Lir5RSfkiLv1JK+SEt/kop5Ye0+CullB+y7AxfpQbDaxtL+rX+rXnpg5REKc+mLX+llPJDWvyVUsoPafFXSik/pMVfKaX8kBZ/pZTyQ1r8lVLKD2nxV0opP6Tj/JXH6+/YfaVU37Tlr5RSfkiLv1JK+SEt/kop5Ye0+CullB/S4q+UUn5Ii79SSvkhLf5KKeWHdJy/8not7Z2U1DZTUttMaW0zzW2dXz4WFRpIWkwoqbFhpMeGEWjX9o5SoMVfebFj9S2sL6xhR2kd7Z0GAZKjQhgWEgiAwVBxooW9R08AEBxgY3JqNNMzYkiLCUVELEz/T3oBGmUFLf7K6xw70cLKXUfJr2wkwCZMSYtmcmo0aTGhBAfav7Z+c2sHpXXN7CyrZ0dpHZuLakmJDuX8MQk4HAabzTM+BJQaSlr8lddobu3go/0VbCysJSTQzmUTkpmREUNY8JnfxmHBAeQkDyMneRhXtY/gi7LjfJ5fzWubSthUVMu3LxjFVeeM0C4h5Ve0+CuvsLu8nrd3lHOyvZO8kXHMH5vYZ9HvTUignbysOHIzYtl9pJ4vSo/z0LIv+M2HB/nW+aO4cXoqIb18e+gvnY9IeTot/sqj1TW1sXRzCTvL6hkRHcI3p2UxPCp0wNu124RzUqP51bWT+Hh/JU+uLuBnf9vN71blc+95WdyWl0H4WXy49Icxhk5jcDggKEC/daihpcVfeax1h6p5cOkOqhtbmT8ukfPHJGJ3c/+8zSbMH5/ExeMSWX+ohqfWFPCfK/bz9JpDXD8tlWumpDAxZdhZHxx2GENdUxvlx09ytL6F6sZWahrbqGtuo63DgXGuFxJoIzY8iLjwYMYmRzJu+DC3fANR6nS0+Cu36U9Xx5lGrHR0OvjtqnyeXF1AVnw4N+WmMSJ64K39MxER5oyOZ87oeLaX1LHk00JeWl/Enz4/zMj4cOaNSWBqejRT02JIiQn92oeQMYaapjaKa5o5XN3E33ceofx4C0frT9La4QDAJhAbHkxceBBZCeGEBNgIsNsQoP5kO3XNbRTXNLGrvJ4Am5CTHMmFOYmD/rMr/6TFX3mUsrpmfrB0B1uL67gpN5WfL5zA37YfGdIMU9Nj+MPt06lvbmfl7qO8t/Moy7aU8sK6IgBEus4fiA0LosNhaG7roLG1g5Z2x5fbCLQLycNCmJIWTUp0KCOiQ0kcFkyA7czdOw5jKK1tZmd5PTtKjrP3SAG5mbFcMj6JiEHuhlL+Rd9NymOs3HWUf/vLThwGfrtoCldPSbE0T1RYIItmprNoZjodnQ4OVDTwRWk9FSdaqGlqpa65nUCbEBoUQHiQnRHRoWTGh5ERF866gpqz6qKyiZARF05GXDjzxybx8f4K1hfWsKv8ODdMS2P8iGFu+4al/JsWf2W51o5OfvneXl7ZUMI5adH8ftFU0uPCrI71FQF2GxNGRDFhRJRL628srB3wa4YG2bly8ghmZMXy1tYyXtlYzMXjErkwJxGbh5ygpryXFn9lqdLaZr7z2jZ2ltXzrXkj+eFlOTrevofEyBDuPW8kf9tezqp9lRyrb+Gm3DTdT2pAXHr3iMgCETkgIgUi8nAvj6eLyGoR2S4iO0XkCvdHVb5mzYFKrnrycw5XN7Hkjun85IpxWtBOI9Bu44bpqVwxaTh7j5zg1Y3FdHQ6+n6iUqfR51+aiNiBp4DLgfHALSIyvsdqPwWWGWOmAouAp90dVPkOhzH85sOD3PXCZpKHhfDud+dy6YRkq2N5PBFh7uh4rp2awsGKRl7bVEKHQz8A1NlxpdtnJlBgjCkEEJGlwNXA3m7rGGCY83YUMLTDM5TXaGrtYNmWUvIrG7lheiq/vHoioUE6nr0/cjNj6TSGt3ccYemmUm6Zme728x+U73Ol+KcApd3ulwF5Pdb5OfAPEfkeEA7M721DInIfcB9AerqOQvA3lQ0tvLiuiIaWDh6/bhI3z0jzmJk1vU1eVhydDsN7O4/y911HWXjOCKsjKS/jSvHv7a/T9Lh/C/CCMeb/ichs4GURmWiM+cp3UmPMEmAJQG5ubs9tKB92uLqJVzYUY7MJ9543kkUzve/D39Pm65kzKp7jze18XlBNYmQws0bGWR1JeRFXin8ZkNbtfipf79a5G1gAYIxZLyIhQDxQ6Y6QyrvtLDvOm1vLiA0PYvHsTGLDg6yO5DMWTEymurGV93YeIS4iiOzESKsjKS/hSvHfDGSLSBZQTtcB3Vt7rFMCXAy8ICLjgBCgyp1BlXfaVlzHX7aVkREXzh2zMjyuf9/TWvP9ZRPh5tw0nv20kNc3lfCdC0YTFxFsdSzlBfoc7WOM6QC+C3wA7KNrVM8eEXlMRBY6V/tX4F4R+QJ4HbjTGKPdOn5uc1Etf9lWxqiECO6ck+lxhd9XBAfauWNWBoLw+qYS2nUIqHKBSyd5GWNWACt6LHu02+29wLnujaa82ebDtSzfUc6YpAhuy8vQ8fuDLCY8iBump/LyhmJW7Dpq+dQYyvPpX6Ryu93l9fxtRzk5SZHcroV/yIwbPozzRsez8XAtO8uOWx1HeTj9q1RuVVjdyLItpaTFhnHLzHQCtPAPqUsnJJMeG8Zft5dT3dBqdRzlwfQvU7nNsfoWXtlQTEx4EN+YlaFXp7KA3SYsmpGGXYTXNpXQ0t5pdSTloXRiN+UWNY2tvLShiCC7jbvmZPZ5fV1vH2XjyaLDgrgpN5UX1xfzi3f38KvrJlsdSXkgbZqpAWvvdPDtV7fR2NLB7bMyiA7TcfxWy0kexrzsBF7fVMrftpdbHUd5IC3+asB++d5eNh6u5bppKaTGeNY8/P7skvFJzMiM4f8u30VBZYPVcZSH0eKvBmTZ5lJeWl/MffNGMiUtxuo4qhu7Tfj9LdMIDbTzwCvbaG7rsDqS8iBa/NVZO3CsgZ+9vZtzR8fxbwvGWh1H9SI5KoTfLppKQVUjjyzfjZ57qU7R4q/OSnNbB995bRuRIYH8781TdUphDzY3O54HLx7D8u3lvL6ptO8nKL+gxV+dlUff3sOhqkZ+u2gKCZE6l4yn+95FozkvO56fv7uH3eX1VsdRHkCHeqoz6m1I5raSOt7aWsaFOYkU1zRTXKPDNj2dzSb8781T+D+//5wHXt3Ke987j6jQQKtjKQtpy1/1S2VDC+/sOEJWfDgXj0u0Oo7qh7iIYJ68dRpHj7fwoze/0P5/P6fFX7msvdPB0k2lBNiFm3LTsOlVuLzO9IwYHr58LP/YW8Fznx22Oo6ykBZ/5bK/7zzKsRMt3JSbpl0GXuzuuVksmJDM4+/v57N8veyGv9Lir1yys+w4m4pqmZedwJgkvVqUNxMR/vvGyYxOiODbr26joLLR6kjKAnrAV/WpprGV5dvLSY8N45LxSVbHUf1wpjmUFp4zgqfXFHDzs+t54PxR3DNv5BAmU1bTlr86o45OB0s3l3ZdLnBGmo7n9yEx4UHcPiuD+pPtvLqphNYOnQHUn2jxV2e0cs8xyo+f5PppqcTohG0+JyMunOunpXK4uokHl+6g06EjgPyFFn91Wu/vPsb6QzWcOyqO8SOGWR1HDZJz0qK5ctJwVu4+xk//tkuHgPoJ7fNXvSqobOSHb35Bakwol01MtjqOGmTnjo4nMz6Mp1YfIjosiB9floPoUF6fpsVffU1DSzvfenkLwQE2bp2ZToBNvyD6gx9emkNtUzt/WHOIAJvw0CVj9APAh2nxV19hjOGHb35BUU0zr9ydx+HqJqsjqSEiIvzHNRMxxvD7jwvodBh+pN8AfJYWf/UVv1tVwAd7KvjpleOYPSpOi7+fsdmE/7x2Ejab8PSaQ3Q6DA9fPlY/AHyQFn/1peXby3jio4NcNy2Fu+dmWR1HWcRmE/796onYRXj200Jqmtr41XWTCLRr958v0eLvA/pzMfRb89J7Xb7+UA0/fmsns0fG8fh1k7Wl5+dsNuGxqycQHxHMEx8dpKqhladvm0Z4sJYMX6Ef5Yr8iga+9fIWMuPCeeaO6QQF6NtCdR0D+MH8bB6/bhKf5Vdxyx83UN3YanUs5Sb6V+7niqqbuO25jQQH2nn+zhk6YZv6mkUz0/njN3I5WNHA9X9YR5EeB/IJ+h3Oj5XVNXPbcxvpcBjeuG8WabFhVkdSFuqr+/CuOVm8uL6IK3/3Ga/dO4tz0qKHJpgaFNry91PH6lu49Y8baWhp56VvziRbZ+pUfUiLDeP+eaMICrCxaMkGPt5fYXUkNQBa/P1QSU0zNz67jtqmNl785kwmpkRZHUl5ifjIYO4/fxSjEyO496WtLN2kl/D0Vlr8/Ux+RQM3PruOhpYOXr0nj6npMVZHUl4mMiSQpffNYu7oeB7+6y6e+PCgzgfkhbT4+5Hy4ye5eckGHAbeuG+29tmqsxYeHMBzi3O5cXoqv12Vz2Pv7dUPAC+jB3z9RFF1Ey+uLyI+IphX78kjMz7c6kjKi506ODwlLZry4yf589oidpXVc83UlK9d2/l055Yoa2nx9wP5FQ28srGYqNBA3rx/NiOiQ62OpHyEiHDlpOEEB9hYfaCKDofhhumpX/sAUJ7HpW4fEVkgIgdEpEBEHj7NOjeJyF4R2SMir7k3pjpb+46e4KUNxcRHBHPfvFFa+JXbiQiXjE/m0vFJ7Cg9zts7jmgXkBfos+UvInbgKeASoAzYLCLvGGP2dlsnG/gJcK4xpk5EEgcrsHLdniP1vL6phBHRodw1J4vQILvVkZQPuyAnkdYOB58crCI00M4CvQ6ER3Ol22cmUGCMKQQQkaXA1cDebuvcCzxljKkDMMZUujuo6p/d5fUs3VxCSnQod52bRUigFn41+C4dn0RLeyef5lcRFmRn3pgEqyOp03Cl+KcApd3ulwF5PdYZAyAiawE78HNjzPtuSaj6bd/REyzdXEJqTBh3zsnUwq+GjIhw1TkjONneyft7jhEfEWx1JHUarvT593bkpmeHXgCQDVwA3AI8JyJfG0coIveJyBYR2VJVVdXfrMoFJTVNLN18qqtHC78aejYRrp+WSkp0KMu2lnKwosHqSKoXrhT/MiCt2/1U4Egv67xtjGk3xhwGDtD1YfAVxpglxphcY0xuQoJ+HXS3yoYWXlxfzLCQQL4xO5NgLfzKIoF2G7fPyiDIbuO+l7ZQ39xudSTVgyvdPpuBbBHJAsqBRcCtPdb5G10t/hdEJJ6ubqBCdwZVZ9bQ0s4L64qw2YQ752QSofOuK4tFhQZyW146z312mBufXcfi2ZkuXSdCzwsYGn22/I0xHcB3gQ+AfcAyY8weEXlMRBY6V/sAqBGRvcBq4EfGmJrBCq2+qtNhWLq5lKbWDhbPziBO+1mVh8iIC+fKycM5WNHIhkItCZ7EpeahMWYFsKLHske73TbAQ85/aoh9uLeCw9VN3Dg9ldQYnZZZeZa8rFgOHGtg5e5jjEyIIGlYiNWRFDq3j9f7cG8Fn+ZXMTMzVidpUx5JRLhuWgrBATaWbSmlw+GwOpJCp3fwaqW1zTy0bAcp0aFcOXm4S8/pz/V+lXKXyJBArpuWyssbilm1r5LLJugJYFbTlr+XcjgMP35rJ8bALTPTCbTrr1J5tnHDhzE9I4bP8qs4Wn/S6jh+TyuGl3ptUwnrC2t45MpxxIYHWR1HKZdcPjGZkEA7b+84gkPn/7GUFn8vVFrbzK9W7GPu6HgWzUjr+wlKeYiwoACumDicktpmthbXWR3Hr2nx9zLGGB7+604AHr9+kkvjppXyJFPTo8mKD+f93cdobO2wOo7f0uLvZf6yrZy1BTX85IpxOqxTeSURYeE5I2jt6OSD3cesjuO3tPh7kYaWdh5fuZ+p6dHcOlPPglTeK2lYCOeOimdbSZ0e/LWIFn8v8rtV+dQ0tfLzqyZgs2l3j/JuF+QkEhJo531t/VtCi7+XKKhs5M9ri7hpeppeeF35hNAgOxeOTSS/spH8Sp35c6hp8fcCxhh+8e4eQoPs/GhBjtVxlHKbWVmxxIQF8v7uYzr0c4hp8fcCqw9U8ll+NQ/OH6MXx1A+JcBu49IJyRytb2FH6XGr4/gVLf4ertNh+PXKA2TGhfGN2RlWx1HK7SalRJESHcqqfRV0OrT1P1R0bh8PdWoOnm3FdRyoaGDRjDTe3FJmcSql3M8mwvxxiby4vpjtJXXcoY2cIaEtfw/W3ungo30VpESHMjElyuo4Sg2aMUmRpMaEsvpAJe2dOuvnUNDi78E2FtZw/GQ7CyYmY9MzeZUPExEuGptIXXM7y7eVWx3HL2jx91At7Z2sPlBFdmIEoxIirI6j1KDLSYokJTqU36/O19b/ENDi76HWHqrmZHsnl+q858pPiAgXj0uktPaktv6HgBZ/D1R/sp21BdWMHz6MlOhQq+MoNWRykiKZlBLF02sKdOTPINPi74FeWFtES7uDi8YmWh1FqSElIjxwwSiKapr5xx6d9mEwafH3MCda2vnT54WMGz6MEdrqV37osgnJZMaF8cwnhzB61u+g0eLvYV5YW8SJlg4u1la/8lN2m3DvvJF8UVbPhsJaq+P4LC3+HuRESzvPfVbIJeOTtNWv/Nr101KJjwjimU8OWR3FZ2nx9yCnWv0/uDjb6ihKWSok0M5d52bxycEq9h09YXUcn6TF30N09fUfZv64JD2bVyng9rwMwoLsLPm00OooPkmLv4d4cW0R9SfbeXC+tvqVAogKC+Sm3DTe23mEyoYWq+P4HC3+HuBESzvPaatfqa9ZPCeT9k7z5USHyn10Vk8PoK1+pf6pZ6EfkxTBnz47TGx4EAG2r7ZXb83Ta1mfLS3+Q6i31ktLeydPrznE2ORIdpbVs7Os3oJkSnmuOaPieWFdEbvLTzBFL2HqNtrtY7H1hTWcbO/k4rFJVkdRyiONTowgPiKI9YeqrY7iU7T4W6ilvZPP86sZmxxJSoyO61eqNzYRZo+Mo7TuJKW1zVbH8Rla/C2krX6lXDMtPYbgABvrC2usjuIztPhbRFv9SrkuONDO9IwYdpXV09DSbnUcn6DF3yIbtNWvVL/MHhmHwxg2Hdb5ftxBi78FWto7+Uxb/Ur1S1xEMGOSItl4uJYOh17pa6BcKv4iskBEDohIgYg8fIb1bhARIyK57ovoe9Y5r9Kl8/Ur1T+zR8XR2NrB7nIdEj1QfRZ/EbEDTwGXA+OBW0RkfC/rRQLfBza6O6QvOdnWyefOq3SlxoRZHUcpr9I17DOYdYf0wO9AudLynwkUGGMKjTFtwFLg6l7W+yXwX4BOwnEGnxVU0druYP447etXqr9sIsweFUeZDvscMFeKfwpQ2u1+mXPZl0RkKpBmjHnPjdl8TmNrB+sKapiUGkVyVIjVcZTyStPSonXYpxu4Uvyll2VfXltNRGzAE8C/9rkhkftEZIuIbKmqqnI9pY/49GAV7Z0OHeGj1AAEB9qZlhHDrvJ6qhtbrY7jtVwp/mVAWrf7qcCRbvcjgYnAGhEpAmYB7/R20NcYs8QYk2uMyU1ISDj71F6o8kQLGwprmJoeTUJksNVxlPJqeVmxdDoMb2wu7Xtl1StXiv9mIFtEskQkCFgEvHPqQWNMvTEm3hiTaYzJBDYAC40xWwYlsZd6anUBDmO4SFv9Sg1YYmQIIxPCeW1jCZ0Ovcj72eiz+BtjOoDvAh8A+4Blxpg9IvKYiCwc7IC+oKyumdc2lTA9I5bY8CCr4yjlE2ZlxVF+/CSr91daHcUruTSlszFmBbCix7JHT7PuBQOP5Vue/LgAQbgwx7+6upQaTOOGDyNpWDAvbShm/nj9Rt1feobvICuqbuLNrWXcmpdOdJi2+pVyF7tNuHVmBp8erKKousnqOF5Hi/8g+92qfALtwrcvHGV1FKV8zi0z0wiwCa9uLLY6itfR4j+ICiobWL6jnMWzM0mM1HH9Srlb4rAQLpuYzLItZZxs67Q6jlfR4j+Invgon7BAO986X1v9Sg2WO2ZlUH+ynXd3Hul7ZfUlLf6DZO+RE/x951HunpulI3yUGkR5WbGMSYrglQ3a9dMfWvwHyW8+PMiwkADuPm+k1VGU8mkiwh2zMthZVs+O0uNWx/EaWvwHwY7S43y0r4L75o0kKjTQ6jhK+bxrpqYQHmTn5fXa+neVFv9B8JsPDxIbHsSd52ZZHUUpvxAZEsi101J4d+cR6prarI7jFbT4u9mmw7V8erCK+88fSUSwS+fQKaXc4I5ZmbR1OFi2Ref7cYUWfzcyxvA//zhAQmQwd8zKtDqOUn4lJzmSmVmxvLKxGIfO99MnbZoOwGsbS75yv6CykU2Ha7lq8nCWby+3KJVS/uuOWRl87/XtfJJfxYU5epnUM9Hi7ybGGD7ce4yo0EBmZMZaHUcpv9CzAdbhcBAZHMDjK/Zz9PhXLyp4a176UEbzeNrt4yYHKhoorTvJRTmJBNh1typlhQCbjdzMWA5WNFCrB37PSKuUGxhj+GhvBbHhQUzLiLE6jlJ+bWZWLCJdgy/U6Wnxd4M9R05wpL6Fi8cmYrf1dtVLpdRQiQoNZNzwYWwprqW902F1HI+lxX+AHMbw0b4KEiKCOSct2uo4SikgLyuO5rZOdpfXWx3FY2nxH6CdZfVUNrRy8bhEbKKtfqU8waiEcOIjgtlQWGN1FI+lxX8AOh2GVfsqSB4WwsSUKKvjKKWcRIRZI2MprTtJ+fGTVsfxSFr8B2BHaR01TW1cMj5JW/1KeZipaTEE2oWN2vrvlRb/s9TW4WDV/kpSY0IZmxxpdRylVA+hQXampEXzRdlxvdBLL7T4n6U3tpRyvLmd+eOSEG31K+WR8rLiaO80bC2pszqKx9HifxZa2jt58uN8MuLCyE6MsDqOUuo0RkSHkh4bxsbCGp3vpwct/mfh1Y0lVJxo5RJt9Svl8WaNjKWmqY3PC6qtjuJRtPj3U1NrB39YU8C5o+MYmaCtfqU83cQRUYQHB/DiuiKro3gULf799OL6Iqob23jokhyroyilXBBgt5GXFcuq/ZUcqmq0Oo7H0OLfDyda2nn2k0IuzElgus7ho5TXyMuKJSjAxp/XHrY6isfQ4t8Pz39+mPqT7drqV8rLRIYEcs2UEby1tUwv8+ikxd9FdU1t/OmzwyyYkMykVD2bVylv8825WbS0O3htU0nfK/sBLf4uWvJZIY1tHfzLJWOsjqKUOgtjk4cxd3Q8L64roq1DZ/vU4u+C6sZWXlhbxFWTR5CjZ/Mq5bXuPi+LyoZW3tt5xOooltPi74I/rDlEa0cnD87PtjqKUmoAzs9OIDsxgiWfFmKMf5/0pdfw7aHnNUHrT7bz4roipqbFsKGwlg2FenUgpbyVzSbcf/4o/vXNL1hzoIoLx/rvRd615d+HNQcqcRjj128SpXzJwikjGBEVwh/WHLI6iqW0+J9BXVMbW4rqyM2MJTY8yOo4Sik3CLTbuHfeSDYV1bK12H+/yWvxP4OPD1QiAhfmaKtfKV9y84w0YsIC+cOaQqujWMal4i8iC0TkgIgUiMjDvTz+kIjsFZGdIrJKRDLcH3VoVTe0sr2kjrysWKJCA62Oo5Ryo7CgABbPyeSjfRUcrGiwOo4l+iz+ImIHngIuB8YDt4jI+B6rbQdyjTGTgbeA/3J30KG2an8Fdpswb0yC1VGUUoNg8exMwoLsPLW6wOoolnCl5T8TKDDGFBpj2oClwNXdVzDGrDbGNDvvbgBS3RtzaFWcaGFnWT2zR8YTGaKtfqV8UUx4EN+Ynck7XxyhoNL/JnxzpfinAKXd7pc5l53O3cDKgYSy2kf7KggKsDEvO97qKEqpQXTveVmEBtr5/cf5VkcZcq4U/96uVtLr2REicjuQC/z3aR6/T0S2iMiWqqoq11MOoSPHT7LnyAnOHR1PWLCeBqGUL4uLCPbb1r8rxb8MSOt2PxX42rnRIjIfeARYaIxp7W1DxpglxphcY0xuQoJn9qV/tK+C0EA7547SVr9S/sBfW/+uFP/NQLaIZIlIELAIeKf7CiIyFXiWrsJf6f6YQ2NbSR37jzVwXnY8oUF2q+MopYbAV1v//jPyp89+DWNMh4h8F/gAsAPPG2P2iMhjwBZjzDt0dfNEAG86r2lbYoxZOIi5B8Vv/nGQ8CA7s0fFWR1FKeVmPadu6S42PIhAu40Hl+7g1rwMbs1LH8Jk1nCpU9sYswJY0WPZo91uz3dzriG3tqCazwuquWLScIIDtNWvlD+JCA7gvNHxrNpfSWltc99P8AF6hi/gcBgeX7mflOhQ8rJirY6jlLLA3NHxhAcH8P6eY34x46cWf2DF7qPsKq/nXy4ZQ6Bdd4lS/ig40M5FOQkcrm5izUHPHI3oTn5f6do7HfzPBwfISYrk2qlnOn1BKeXrZmR1TeL465X76XT4duvf74v/G5tLKapp5keX5WC39XZKg1LKXwTYbFwyPon9xxpYvr3c6jiDyq+Lf3NbB79dlU9uRgwXj9OZO5VSMCklinPSovn1+/tpbO2wOs6g8evi/+e1RVQ1tPLw5WNxDlFVSvk5mwi/WDiBqoZWnz7xy2+Lf11TG8+sOcT8cUnkZuoIH6XUP01Ji+bG6ak8//lhCqt8c9oHvy3+T68poKmtgx8vyLE6ilLKA/14wVhCAuz8+9/3WR1lUPhl8S8/fpIX1xVz3bS0HR7yAAAM4ElEQVRUxiRFWh1HKeWBEiKD+cH8bD7eX8lHeyusjuN2fjFtZc/Tut/aWobDGEbGh5/xlG+llH9bPCeTZVtKefTt3cwaFUeED83063ct//K6k2wvqWPWyDiiw/Si7Eqp0wu02/jVdZM5eqKF//nggNVx3Mqvir8xhvd2HiEsOICLxurQTqVU36ZnxHDHrAxeXF/E9pI6q+O4jV8V/y/K6imubeay8UmEBOrkbUop1/zoshySh4Xw8F920dbhsDqOW/hN8W/rcPD+7qOkRIcyLSPG6jhKKS8SGRLIL6+eyIGKBp70kQu++87Riz58crCSEy0d3DIzHZue0KWUOoPTDQSZmhbNkx/n09bhID02DMBr5/73i5Z/5YkWPs2vZkpaNBlx4VbHUUp5qavOGcGw0ECWbSmltaPT6jgD4vPF3+EwLN9RTpDdxuUTk62Oo5TyYiGBdm6cnkZdUxsrdh21Os6A+Hzxf21TCcU1zVwxKZnIkECr4yilvFxWfDjnZcezuaiOPUfqrY5z1ny6+B+rb+HXK/czMiGcael6kFcp5R7zxyWRGhPKW1vLOFzdZHWcs+Kzxd8Yw8/e3k1bp4Nrp6TorJ1KKbcJsNu+HDzywCtbOdnmff3/Plv8l24u5cO9Ffzw0hziIoKtjqOU8jExYUHcPCONAxUNPLJ8l9dd99cni39BZQO/eHcPc0fHc/fcLKvjKKV81JikSH5wcTZ/3V7O82uLrI7TLz43zr+1o5Pvvb6DsKAAfnPTOdj00oxKqUH0/Yuy2Xf0BP/+972kxYRy6QTvGFXocy3/X63Yz76jJ/jvGyaTOCzE6jhKKR9nswn/e/NUJqdG8/2l2/mi9LjVkVziU8X/jc0lvLCuiLvOzeTicUlWx1FK+YnQIDvPfSOX+Ihg7n5xM8U1nj8CyGeK//pDNTyyfDfnZcfzyBXjrI6jlPIzCZHBvHDXDDochluWbKC0ttnqSGfkE8X/cHUTD7y6lYy4MJ68dRoBdp/4sZRSXmZ0YiSv3pNHU1snizz8A8Drq2RZXTOLn9+EAM/fOYOoUD2LVyllnQkjonj1njwaWtq55Y8bKKnxzA8Ary7+pbXN3PzsBuqa23j+zhk6aZtSyiNMTInilXvyaGzt4Nqn17LNAy8C47XFv7imiZufXU9jawev3pPHVJ2+QSnlQSanRvPXB+YQERLALUs2sNLDJoLzynH+6w5V873XttNpDK/ek8fElCirIyml/NTp5v4/5ba8DF7ZUMwDr25jXnYCzy3OJSjA+na39Qn6wRjDkk8PcftzG4kOC+St++do4VdKebSI4ADunpvFjMwYPs2v4oZn1nnEZHBeU/zL6pq596Ut/OeK/Vw2IZm3vzuX0YkRVsdSSqk+BdptXDs1lVtnplNc08yVv/uMZz45ZOkFYTy++Ld1OHhqdQHzf/MJawtq+OmV43j6tmlEBHtlj5VSyo9NTIni/QfPY86oOB5fuZ9Ln/iUD/Ycs2RSOI+toCda2nljUyl/XnuYI/UtLJiQzM+uGk9KdCjQdz+bUkp5ouFRoTy3eAafHqzil+/t5Vsvb2VsciR3z81i4ZQRBAfYhySHS8VfRBYAvwXswHPGmMd7PB4MvARMB2qAm40xRf0N09HpYOPhWlbsOsrbO47Q2NrBrJGxPH79ZOaNSejv5pRSymPNG5PAyh+cx/Lt5fzp88P86K2d/Pr9/VwxaThXTBrOjMxY7IM4MWWfxV9E7MBTwCVAGbBZRN4xxuztttrdQJ0xZrSILAJ+Ddzc17Zb2js5cKyBbSV1bC2uY21BNXXN7YQG2rlsQhJ3zx3JpFQ9oKuU8k0Bdhs35qZxw/RU1hbU8OrGYt7YXMpL64uJjwhiZlYsMzJjyc2IZXRiBKFB7vtW4ErLfyZQYIwpBBCRpcDVQPfifzXwc+ftt4AnRUTMGTqyDhxrYNyj73NqjRFRIVyQk8hlE5I5f0yCW39IpZTyZCLC3Ox45mbH09Tawcf7K1m1r4LNRXWs2HXMuQ6kxoS67TVdKf4pQGm3+2VA3unWMcZ0iEg9EAdUn26joUF2vn9RNmOSIpmWEc3wKPf9UEop5a3CgwO46pwRXHXOCADKj59kR8lxCiobKahq5HM3vY4rxb+3TqeeLXpX1kFE7gPuc95tfejSnN0uvL7V4jnDh5gH0Zzu4w0ZQXO625DkvG3gm8gZ+CZcK/5lQFq3+6nAkdOsUyYiAUAUUNtzQ8aYJcASABHZYozJPZvQQ0lzupc35PSGjKA53c2bcrpjO66M898MZItIlogEAYuAd3qs8w6w2Hn7BuDjM/X3K6WUslafLX9nH/53gQ/oGur5vDFmj4g8BmwxxrwD/Al4WUQK6GrxLxrM0EoppQbGpXH+xpgVwIoeyx7tdrsFuLGfr72kn+tbRXO6lzfk9IaMoDndza9yivbOKKWU//H4uX2UUkq536AWfxG5UUT2iIhDRE57FF1EFojIAREpEJGHuy3PEpGNIpIvIm84DzgPRs5YEfnQ+TofisjXrgwjIheKyI5u/1pE5BrnYy+IyOFuj02xKqdzvc5uWd7ptnzQ96eL+3KKiKx3vjd2isjN3R4b1H15uvdat8eDnfumwLmvMrs99hPn8gMicpk7c51FzodEZK9z/60SkYxuj/X6+7co550iUtUtzz3dHlvsfJ/ki8jins8d4pxPdMt4UESOd3tsSPaniDwvIpUi0usQeOnyO+fPsFNEpnV7rP/70hgzaP+AcXSNSV0D5J5mHTtwCBgJBAFfAOOdjy0DFjlvPwM8MEg5/wt42Hn7YeDXfawfS9eB7TDn/ReAGwZzX/YnJ9B4muWDvj9dyQiMAbKdt0cAR4Howd6XZ3qvdVvn28AzztuLgDect8c71w8GspzbsVuY88Ju778HTuU80+/fopx3Ak/28txYoND5f4zzdoxVOXus/z26BrYM9f6cB0wDdp/m8SuAlXSdVzUL2DiQfTmoLX9jzD5jzIE+Vvty+ghjTBuwFLhaRAS4iK7pIgBeBK4ZpKhXO7fv6uvcAKw0xgz1lZn7m/NLQ7g/+8xojDlojMl33j4CVAJDMXNfr++1Hut0z/8WcLFz310NLDXGtBpjDgMFzu1ZktMYs7rb+28DXeffDDVX9ufpXAZ8aIypNcbUAR8CCzwk5y3A64OU5bSMMZ/Sy/lR3VwNvGS6bACiRWQ4Z7kvPaHPv7fpI1Lomh7iuDGmo8fywZBkjDkK4Pw/sY/1F/H1N8d/OL+KPSFds5wOBldzhojIFhHZcKpriqHbn/3alyIyk67W2KFuiwdrX57uvdbrOs59dWqqEleeO5Q5u7ubrhbhKb39/geDqzmvd/4+3xKRUyeMeuT+dHafZQEfd1s8VPuzL6f7Oc5qXw54Pn8R+QhI7uWhR4wxb7uyiV6WmTMsPytnytnP7QwHJtF13sMpPwGO0VXElgD/BjxmYc50Y8wRERkJfCwiu4ATvax3VvvTzfvyZWCxMcbhXOy2fdnbS/ayzNWpStz6fuyDy68lIrcDucD53RZ/7fdvjDnU2/OHIOe7wOvGmFYRuZ+ub1UXufhcd+nPay0C3jLGdL/E1lDtz7649b054OJvjJk/wE2cbvqIarq+1gQ4W2C9TSvhsjPlFJEKERlujDnqLEiVZ9jUTcByY0x7t20fdd5sFZE/Az+0MqezKwVjTKGIrAGmAn/BTfvTHRlFZBjwd+Cnzq+wp7bttn3Zi4FMVeLKc4cyJyIyn64P3PONMa2nlp/m9z8YxarPnMaYmm53/0jXdO+nnntBj+eucXvCf76Wq7+7RcB3ui8Ywv3Zl9P9HGe1Lz2h26fX6SNM15GM1XT1r0PX9BGufJM4G92np+jrdb7WH+gscqf61a8BBmvCuj5zikjMqa4SEYkHzgX2DuH+dCVjELCcrv7LN3s8Npj7ciBTlbwDLJKu0UBZQDawyY3Z+pVTRKYCzwILjTGV3Zb3+vu3MOfwbncXAvuctz8ALnXmjQEu5avfpoc0pzNrDl0HTNd3WzaU+7Mv7wDfcI76mQXUOxtLZ7cvB/no9bV0fSq1AhXAB87lI4AVPY5iH6Tr0/SRbstH0vUHVgC8CQQPUs44YBWQ7/w/1rk8l64rl51aLxMoB2w9nv8xsIuuQvUKEGFVTmCOM8sXzv/vHsr96WLG24F2YEe3f1OGYl/29l6jq1tpofN2iHPfFDj31chuz33E+bwDwOWD/LfTV86PnH9Tp/bfO339/i3K+StgjzPPamBst+d+07mfC4C7rMzpvP9z4PEezxuy/UlXo/Ko82+jjK5jOfcD9zsfF7ourHXImSW323P7vS/1DF+llPJDntDto5RSaohp8VdKKT+kxV8ppfyQFn+llPJDWvyVUsoPafFXSik/pMVfKaX8kBZ/pZTyQ/8fNNRScAsd16YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "spurious_reg(0.99, 10, 0, 300, 5000)  # rho = .99"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we change the series to length 2000 when $\\rho $ is high?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0XPWZ5//3U6V9X63NlmUb7wu2UdssAQwBYgIBMtATm6QbuslhyDI9J9OzJJP+Jd0k3ZP0MpnuTjqJkxCS9IEEnAQcYhazGAhesIx3GduyvGixdkvWLlXV8/tDZU9hJKskVdWtUj2vc+qo6i5VH11V6an7/d77vaKqGGOMMQAupwMYY4yJHlYUjDHGXGJFwRhjzCVWFIwxxlxiRcEYY8wlVhSMMcZcYkXBGGPMJVYUjDHGXGJFwRhjzCUJTgcYTUFBgVZUVDgdwxhjYsbevXvbVLVwqs8TlUWhoqKCqqoqp2MYY0zMEJEzoXgeaz4yxhhziRUFY4wxl1hRMMYYc4kVBWOMMZdYUTDGGHOJFQVjjDGXWFEwxhhziRUFY4wxl1hRMMYYc0lUntFsTDR5avfZUac/uLY8wkmMCT/bUzDGGHOJ7SkYE2Kj7VnYXoWJFVYUjJkgj8/H6bY+/vb31XT0DpOUICwpzeaTq8rISLaPlIlt9g42JkjDXh/vnupg+/FWegc9JLiEjOQEhr0+nn63jm++UM2N8wu5ZWEhIuJ0XGMmZdyiICJPAHcDLaq6bJT5/x34dMDzLQYKVbVDRE4D3YAX8KhqZaiCGxNJNS09/GZfPZ19w8wrTOe6uWVcNSODpAQXqkr9+X7ePN7Kq0ebaeke4P7VM0l0W5ediT3B7Ck8CXwX+PloM1X1H4B/ABCRTwBfUtWOgEVuUdW2KeY0JqSCbfcf9HjZcqCBXbUdFGQk88hH5jCvMOMDy4gIs/LS+PTact483sor1c30D3l56PoKXLbHYGLMuEVBVd8SkYogn28j8PRUAhkTLVq7B3ns3/ey98x5bpiXzx1Li6/47V9EWLdwBqlJbp7f38hrR1u4fUlRBBMbM3Uh61MQkTRgPfDFgMkKvCIiCvxQVTeF6vWMCaez7X1s/NEu2nsH2bimnOVl2UGvu3ZOPvUd/bxxrIXyvDQWFmeGMakxoRXKRs9PAO9c1nR0g6quBu4EviAiN421sog8KiJVIlLV2toawljGTEz9+ZGC0Dvk4dn/dP2ECsJFn7i6lOKsFH79Xj0Dw94wpDQmPEJZFDZwWdORqjb6f7YAvwXWjLWyqm5S1UpVrSwsnPK1p42ZlM6+IT79491cGBjm3x9Zy/KZEy8IAEkJLv7D6jJ6Bj28cawlxCmNCZ+QNB+JSDZwM/CZgGnpgEtVu/337wAeD8XrGRMOPp/yl88coLGzn18+ei3LJrGHEGhmbhrXlOeyo6adf3n1BAWZyR9axk5qM9Fm3D0FEXka2AksFJF6EXlERB4TkccCFvsk8Iqq9gZMKwL+ICIHgHeB36vqS6EMb0wo/ejtWl57v4W/umsJ18zOC8lz3rG0iAS3sPXwuZA8nzHhFszRRxuDWOZJRg5dDZxWC1w92WDGRNK5rn7+bftJ7lpewp9eNztkz5uZksjNCwp5pbqZhvP9lOWmhuy5jQkHO7vGxD1VZcuBRrJTE/m7Ty4P+dnI187NJyXRxfbj1rdgop8Nc2Hi3oH6Ts609/EfVpXx+0Ohb+ZJSXRz3dx83jjWSvOFAYqyUkL+GsaEiu0pmLg2OOzlxUNNzMxNZfXs3LC9zvXzCkh0C28et8OtTXSzomDi2s7adroHPdy9ojSsQ1KkJyewpiKPg/WdXBgYDtvrGDNVVhRM3BoY9vL2iTYWFmVSnpcW9te7dm4+PoU9pzvGX9gYh1hRMHFrV207/cNePrp4RkReLz8jmfkzMthzqgOvTyPymsZMlBUFE5cC9xJm5oZ/L+Gia+fmc2HAw9FzFyL2msZMhBUFE5f2nO6gf9jLrYsis5dw0cLiTHLSEtlV2x7R1zUmWFYUTNwZ8vh4p6aNuQXpzIpAX0IglwhrKvKobeulvWcwoq9tTDCsKJi4s+VAIxcGPNw435mBF1eV5yLAvrpOR17fmCuxomDiiqqy6a2TFGUls6AoY/wVwiA7NZGrZmTw3tnz+KzD2UQZKwomrmw/1srx5h5unF8Y8uEsJmJ1eS6dfcPsOmV9Cya6WFEwceWHb52kJDuFFZO8TkKoLCnNIjnBxeaqekdzGHM5Kwombhyo62RXbQd/fsMcElzOvvUT3S5WzMzhxcNN9A15HM1iTCArCiZubHqrlsyUBDasmeV0FACunplN/7CX147a6KkmelhRMHHhTHsvLx4+x2eunU1mSqLTcQCoKEinMDOZFw42Oh3FmEusKJi48P3tJ0lwu/iz6yucjnKJS4S7lpfwxrFWegatCclEB7uegpkWntp9dtTpD64tp/58H5v31vPpteXMiLJrGdy9ooQnd5zm1epm7ltV5nQcY2xPwUx/399+EpcIj62b53SUD1ldnktxVoo1IZmoYUXBTGvnuvp5tqqeBypnUpIdfddHdrmEu1aU8ObxVrr67ToLxnnjFgUReUJEWkTk8Bjz14lIl4js99++FjBvvYgcE5EaEflyKIMbE4z/88pxAD53c/TtJVx094oShr3Ktupmp6MYE9SewpPA+nGWeVtVV/pvjwOIiBv4HnAnsATYKCJLphLWmIlo7Oxn83v1PHxDRcQHvpuIlbNyKMtJtSYkExXGLQqq+hYwmUtFrQFqVLVWVYeAXwL3TuJ5jJkwVWXroXPkpiXxhVuucjrOFYkId68o4Q8n2jjfO+R0HBPnQtWncJ2IHBCRF0VkqX9aGVAXsEy9f9qoRORREakSkarWVru4uZmao+e6qW3r5Uu3zSc7NTrOS7iSu1eU4vEpLx9pcjqKiXOhKArvAbNV9WrgX4Hn/NNHG21szCEhVXWTqlaqamVhoTNDGpvpwePz8eLhcxRmJrNxTbnTcYKyrCyL2flp/P7QOaejmDg35aKgqhdUtcd/fyuQKCIFjOwZBI4nMBOwRlMTdrtrO2jvHeLjy4pJcMfGAXYXm5B2nGy3i+8YR035EyMixeIfg1hE1vifsx3YA8wXkTkikgRsALZM9fWMuZK+IQ+vv9/C/BkZLCjKdDrOhNy5rASvT3n1qB2FZJwz7hnNIvI0sA4oEJF64OtAIoCq/gB4APiciHiAfmCDqirgEZEvAi8DbuAJVT0Slt/CGL/X329hYNjLnctKHL1ewmQsLc1iZm4qLx1u4lN/FBvNXmb6GbcoqOrGceZ/F/juGPO2AlsnF82YiWnrHmRXbTuVFXkUZ0fXcBZjuXx4jor8dN460caFgWGyomTgPhNfYqPB1ZggvHj4HIluF7ctnuF0lElbWpqF16e88b4Np22cYUXBTAsnW3s42tTNugWFUTM09mTMyksjMyWBlw7boanGGVYUTMxTHTm+Pzs1keuvKnA6zpS4RFhSksX2Y630D3mdjmPikBUFE/NePtJE/fl+bls8g8QYOQT1SpaWjlyR7c3jdhKnibzY/wSZuObx+vjHV45TmJHMylm5TscJiTkF6eSkJdrZzcYRdpEdE1MuP1pn75nz1LT08OCactyu2DoEdSxul3D74iJeOtLEkMdHUoJ9dzORY0XBxCyfKm8eb6EkO4WlpVmjLjPWFdmi3fplxTy7t54dJ9tYtzB2j6Yysce+gpiYdaTxAm09Q9y8oDDmTlQbzw1XFZCe5LYmJBNxVhRMTFJV3jreSn56EsvKsp2OE3IpiW5uXVzEK0ea8frGHEfSmJCzomBiUk1rDw2d/dy0oBDXNNtLuGj90mLae4fYc3oylzMxZnKsT8HEpHdq2shMSWDVrByno4TFU7vPMujxkuAS/vm1E9S29gLw4FobE8mEl+0pmJjT1j3I8eYe1s7Ji5mhsScjOcHN/KJMqhsv4FNrQjKRMX0/UWba2nWqHbcIf1SR53SUsFtamkVX/zAN5/udjmLihBUFE1MGPV72njnPsrKsmB7jKFiLi7NwyciRVsZEghUFE1P2ne1k0OPjunmxPcZRsFKT3MwrzOBIYxdqTUgmAqwomJiy53QHpdkpzMpNdTpKxCwpzaK9d4jmC3aZThN+VhRMzDjc0MW5rgGuqcibdierXcmSkiwEONLY5XQUEwesKJiYsXlvPQku4eqZ0+9ktSvJTElkdn6a9SuYiLCiYGLCoMfLc/sbWFySRVpS/J1es7Q0m6YLA5xq63U6ipnmxi0KIvKEiLSIyOEx5n9aRA76bztE5OqAeadF5JCI7BeRqlAGN/Hl1eoWOvuGqZw9PYbHnqiLA/7ZFdlMuAWzp/AksP4K808BN6vqCuAbwKbL5t+iqitVtXJyEY2BZ6rqKMlOYd6MDKejOCInLYmynFResgHyTJiNWxRU9S1gzMFXVHWHqp73P9wFzAxRNmMAONfVz9snWnngmpnTdpyjYCwtzeJAXSeNnXYimwmfUPcpPAK8GPBYgVdEZK+IPBri1zJx4jfvNeBTeOCa+P6+sax0pIPdhtM24RSyoiAitzBSFP5nwOQbVHU1cCfwBRG56QrrPyoiVSJS1dpq16Y1I1SVZ6vqWDsnj9n56U7HcVRBZjILijKsX8GEVUiKgoisAH4M3Kuq7Renq2qj/2cL8FtgzVjPoaqbVLVSVSsLCwtDEctMA3tOn+d0ex//sXKW01Giwvqlxew53UFbj53IZsJjykVBRMqB3wB/oqrHA6ani0jmxfvAHcCoRzAZM5Znq+rISE7gzuXFTkeJCuuXleBT2Fbd7HQUM02Ne8C3iDwNrAMKRKQe+DqQCKCqPwC+BuQD/+Y/y9TjP9KoCPitf1oC8JSqvhSG38FMQ0/tPsuw18eWA40sK83muX2NTkeKCotLMinPS+Olw01sXGPXVjChN25RUNWN48z/LPDZUabXAld/eA1jgnP03AUGPT5Wlk/PC+lMhoiwflkxP33nFF39w2SnTv+RYk1k2RnNJmodqOskKyWBOQXx3cEc6KndZ3GJMOxVvvlCNU/tPut0JDPNWFEwUalv0MOx5m5WzMyJ63MTRjMzN5WslAQbC8mEhRUFE5UONXbhU1g5Ta/BPBUuEZaUZnOipZshj8/pOGaasaJgotL+uk5mZCZTkp3idJSotLQ0i2Gvcry52+koZpqxomCiTl1HH2fa+1g5KyeurpswERX56aQlue0aCybkrCiYqLPlwMjhp1fPtKajsbhdwpKSLN5v6mbQ43U6jplGrCiYqKKqPLevgdn5aeSmJzkdJ6otLc1i0ONjR037+AsbEyQrCiaqVJ+7wImWHutgDsK8wgySE1w2FpIJKSsKJqo8t6+BBJewvDS+Lrk5GQluF4uKM3mlugmP145CMqFhRcFEDZ9PeeHgOdYtLCQtOf4uuTkZS0uzOd83zLunx7zkiTETYkXBRI19dec51zXA3StKnY4SMxYUZZKS6OJla0IyIWJFwTjuqd1neWr3Wf7pleMkuISO3iGnI8WMpAQXNy8o5KUjTfh86nQcMw1YUTBRwafK4YYu/zdft9NxYsr6ZcU0XxhkX12n01HMNGBFwUSFs+19XBjwsLzMOpgn6tZFRSQluHjhoA0vbqbOioKJCocaukhwCYuKM52OEnOyUxO5deEMfneg0Y5CMlNmRcE4zqfK4cYuFhZnkmxNR5Ny36oy2nqGeOeknchmpsaKgnHcmfY+ugc8LLOmo0m7ZVEhWSkJPL+vwekoJsZZUTCOs6ajqUtOcPPx5SW8dKSJviGP03FMDLOiYBzl9SlHGvxNRwnWdDQV960qo2/Iy7bqZqejmBgWVFEQkSdEpEVEDo8xX0TkX0SkRkQOisjqgHkPicgJ/+2hUAU308Oe0x10D9pRR6GwpiKP0uwUnt9vRyGZyQt2T+FJYP0V5t8JzPffHgW+DyAiecDXgbXAGuDrIpI72bBm+vn9wXMkuoVFxVlOR4l5Lpdwz8oy3jzeSnvPoNNxTIwKqiio6lvAlQZXuRf4uY7YBeSISAnwMWCbqnao6nlgG1cuLiaO+HzKS0eaWFiUSVKCtWSGwn2rSvH6lN8fOud0FBOjQvVJLAPqAh7X+6eNNd0Y9tV10to9yBIbETVkFhVnsag4k+fsKCQzSaEqCqNdM1GvMP3DTyDyqIhUiUhVa2triGKZaLatupkEl7CwyI46CqX7VpXx3tlOzrT3Oh3FxKBQjU9cD8wKeDwTaPRPX3fZ9O2jPYGqbgI2AVRWVtrIXnFgW3UT187NJzXJjjqaiqd2n/3AY59PEeC5fY38l9vmOxPKxKxQ7SlsAf7UfxTStUCXqp4DXgbuEJFcfwfzHf5pJs6dbO3hZGsvty8pcjrKtJOTlsScwnR+s68eVft+ZSYm2ENSnwZ2AgtFpF5EHhGRx0TkMf8iW4FaoAb4EfB5AFXtAL4B7PHfHvdPM3Hu4rH0t1lRCItrynM5095H1ZnzTkcxMSao5iNV3TjOfAW+MMa8J4AnJh7NTGfbqptZWppFWU6q01GmpaWl2Ww9dI7NVfX8UUWe03FMDLFrHpqI6+gd4r2z5/mLW629O1ySElwsLM7iuf0NLC7J+sAhvw+uLXcwmYl2dnC4ibi3T7SiCrcsmuF0lGlt9ewcBj0+qs9dcDqKiSG2p2Ai5uJRMs9U1ZGe5OZwQxfVjfYPK1wq8tPJTUvkvbPnWTkrx+k4JkbYnoKJKJ8qx5u7mV+UiUtGO43FhIpLhFXluZxs6aGrf9jpOCZGWFEwEdVwvp++Ia+dsBYhq2bloMD+s3YUkgmOFQUTUceauxFg/owMp6PEhfyMZCry09h7ttPOWTBBsaJgIup4czez8tJIS7burEhZXZ5LW88g9ef7nY5iYoAVBRMxfUMeGs73M7/I9hIiaVlZNolu4T1rQjJBsKJgIuZUWy8KXFVoRSGSUhLdLC3N5kB9J8Nen9NxTJSzomAipra1l0S3UJZrZzFH2qryHAaGfbzf1O10FBPlrCiYiDnZ2kNFfjoJLnvbRdq8wgyyUhJ4z8ZCMuOwT6eJiNbuQVq6B5lrTUeOuHjOwomWblq6B5yOY6KYFQUTEbtq2wGYV5jucJL4tbo8F5/C8/sanY5iopgVBRMRO062k5LooiTb+hOcUpiZzKzcVDbvtessmLFZUTARsfNkGxX56bhdNrSFk1bPzuVYczdHbMwpMwYrCibsGjv7Od3exzzrT3DcirIckhJcbN5b73QUE6WsKJiw23lypD9hrvUnOC41yc3ti4vYcqCRIY+ds2A+zIqCCbudte3kpiVSlJXidBQD3H9NGR29Q2w/1uJ0FBOFrCiYsFJVdp5s57p5+TZUdpS4aX4hBRnJ1oRkRhVUURCR9SJyTERqROTLo8z/jojs99+Oi0hnwDxvwLwtoQxvot/Zjj4aOvu5bm6+01GMX4LbxSdXlfLGsRY6eoecjmOizLhFQUTcwPeAO4ElwEYRWRK4jKp+SVVXqupK4F+B3wTM7r84T1XvCWF2EwMu9idcN6/A4SQm0P3XzGTYq2zZ3+B0FBNlgtlTWAPUqGqtqg4BvwTuvcLyG4GnQxHOxL4dJ9spzEy2k9aizKLiLJaVZbH5PWtCMh8UzKD2ZUBdwON6YO1oC4rIbGAO8HrA5BQRqQI8wLdU9blJZjUx4uK1mFWVN95vYW5hOk+/WzfOWibS7l89k7/5XTXHmrpZWGxXwjMjgtlTGK13cKzTITcAm1XVGzCtXFUrgQeB/ysi80Z9EZFHRaRKRKpaW1uDiGWiXWv3IN2DHjs/IUrdc3UpCS7h17a3YAIEUxTqgVkBj2cCYw2esoHLmo5UtdH/sxbYDqwabUVV3aSqlapaWVhYGEQsE+1OtvUC2CB4USo/I5l1CwvZsr8Rn8+GvTAjgikKe4D5IjJHRJIY+cf/oaOIRGQhkAvsDJiWKyLJ/vsFwA1AdSiCm+hX29pDTloiuWmJTkcxY/jE1aU0XRhgz+kOp6OYKDFun4KqekTki8DLgBt4QlWPiMjjQJWqXiwQG4Ff6gdH2loM/FBEfIwUoG+pqhWFOOBTpba1lyUlWYidnxBVLvb5AAx5fCS6hX/adpxn/tN1DqYy0SKoq6er6lZg62XTvnbZ478eZb0dwPIp5DMxqqlrgP5hrw1tEeWSElwsLsnicEMXw14fiW47nzXe2TvAhEWt9SfEjBVlOfQNeXmnps3pKCYKWFEwYVHb2kN+ehLZqdafEO0WFGWQkuhiywG7+I6xomDCwOtTTrX12qGoMSLB7WJpaTavHGlmYNg7/gpmWrOiYEKusbOfQY/P+hNiyIqZ2fQMemzkVGNFwYRebWsPYP0JsWRuQQYFGUnWhGSsKJjQq23rpSgrmYzkoA5uM1HA7RLuWl7Ca0db6B4YdjqOcZB9ak1IDXl8nG7vpXJ2ntNRzASlJLoZ9Ph4/HfVrCrPvTT9wbXlDqYykWZ7CiakDtR3MuxVGxU1Bs3KSyMnNZGD9V1ORzEOsqJgQmpHTTsCzCmw/oRY4xJhWVk2NS099A/ZUUjxyoqCCakdJ9soyUkhNcntdBQzCcvKsvGq8n7TBaejGIdYUTAhMzDsZd/ZTubZXkLMmpmbSlZKAocbrSjEKysKJmT2nO5gyOuzQ1FjmEuEpWXZnGjuZtBOZItLVhRMyLx9oo0kt4s5BdbJHMuWlWbj8SnHmrudjmIcYEXBhMxbx1uprMglKcHeVrFsdn4aGcnWhBSv7NNrQqKle4D3m7r5yPwCp6OYKXKJsKQ0i2NNFxjy+JyOYyLMioIJiYvDLt803y6lOh0sK81m2KucaLEmpHhjRcGExNsn2shLT2JJSZbTUUwIzClIJy3JzeEGO5Et3lhRMFOmqrx9oo0brirA5bJLb04HbpewpCSL95u6GfTYUUjxxIqCmbL3m7pp7R7kxqusP2E6WVaWzaDHxx9O2BXZ4okVBTNlrx1tBmDdQutPmE7mFqaTkujixcNNTkcxERRUURCR9SJyTERqROTLo8x/WERaRWS///bZgHkPicgJ/+2hUIY30eHVoy1cPTObGVkpTkcxIZTgcrG4OItt1c0Me+0opHgxblEQETfwPeBOYAmwUUSWjLLor1R1pf/2Y/+6ecDXgbXAGuDrIpI7yromRrV2D3KgvpOPLi5yOooJg6Wl2XT1D7PzZLvTUUyEBLOnsAaoUdVaVR0CfgncG+TzfwzYpqodqnoe2Aasn1xUE43eONaCKnx08Qyno5gwmF+UQVqS25qQ4kgwRaEMqAt4XO+fdrn7ReSgiGwWkVkTXBcReVREqkSkqrW1NYhYJhq8drSZ0uwUOxR1mkp0u7hl0Qy2VTfh9anTcUwEBFMURjvG8PJ3x++AClVdAbwK/GwC645MVN2kqpWqWllYaB2WsWBg2MvbJ9q4dfEMROxQ1OnqzmXFtPUMsed0h9NRTAQEUxTqgVkBj2cCH7i6t6q2q+qg/+GPgGuCXdfErndq2ugb8nKb9SdMa7csnEFygouXrAkpLgRTFPYA80VkjogkARuALYELiEhJwMN7gKP++y8Dd4hIrr+D+Q7/NDMNvHDwHNmpiVw/z85PmM7SkxO4aUEhLx1uwmdNSNNewngLqKpHRL7IyD9zN/CEqh4RkceBKlXdAvyFiNwDeIAO4GH/uh0i8g1GCgvA46pq+6DTwMCwl23Vzdy1vMRGRZ3mntp9lpzURJouDPD3Lx+jPC+NB9eWOx3LhMm4RQFAVbcCWy+b9rWA+18BvjLGuk8AT0who4lC24+10jPo4e6rS8Zf2MS8RcVZuEU40tBFeV6a03FMGNlXPDMpLxxsJC89ievm5jsdxURAapKbeTPSOdzYhao1IU1nVhTMhPUNeXjtaAvrlxWT4La3ULxYWprN+b5hznUNOB3FhFFQzUfGBNpW3Uz/sJfURDdP7T7rdBwTIYtLsnhuXwOHG2047enMvuaZCXumqo5Zeal2LeY4k5GcwJyCdI402GU6pzMrCmZC6jr6eKemnT++ZhYuO2Et7iwty6a1Z5ATzXZFtunKioKZkGf31iMC918z0+koxgFL/cOZ2FhI05cVBRM0r0/ZXFXHjfMLKctJdTqOcUBWaiLleWlWFKYxKwomaG+daKWxa4BPVc4af2EzbS0rzeLouQucae91OooJAysKJmg/23Gawsxkbl9iYx3Fs6Wl2YA1IU1XVhRMUP751RNsP9bKipnZbN5bb4eixrHc9CSWl2VbUZimrCiYoOysbcPtEtZU5DkdxUSBO5cXc6Cuk7qOPqejmBCzomDG1dU/zHtnOrl6ZjaZKYlOxzFR4BMrSgHYcsBGwp9urCiYcT397lmGvD6usyGyjd+svDQqZ+fy3L4GGwtpmrGiYK5oYNjLj98+xVUzMuwwVPMB964q40RLD+832Yls04kVBXNFz+6tp61nkJsX2CVSzQfdtbyEBJfw3P4Gp6OYELKiYMbk8fr44ZsnWVWew1wb58hcJi89iZsWFPK7/Y12RbZpxIqCGdPz+xupP9/PF9Zdhdg4R2YUn1xVRmPXAO+cbHM6igkRKwpmVMNeH//82gmWlmZx66IZTscxUer2JUXkpCXyqz11TkcxIWLXUzCj2ry3nrMdffzkoUpcLttLMB8UePLi4pIsXjzcxI/fquWzN811MJUJhaD2FERkvYgcE5EaEfnyKPP/q4hUi8hBEXlNRGYHzPOKyH7/bUsow5vwGPR4+dfXTrByVo7tJZhxVc7OxetT9tV1Oh3FhMC4ewoi4ga+B9wO1AN7RGSLqlYHLLYPqFTVPhH5HPD3wKf88/pVdWWIc5swemr3WRq7Bvj2AyusL8GMqyQ7lZm5qVSd6UBV7T0T44LZU1gD1KhqraoOAb8E7g1cQFXfUNWL57vvAmyw/RjV1T/Mv7x2guvn5fORq+xkNROcytl5NF8YpOrMeaejmCkKpiiUAYG9SPX+aWN5BHgx4HGKiFSJyC4RuW8SGU0E/dv2Gjr7h/nqXYvtG58J2spZOaQmuvnJ26ecjmKmKJiO5tH+M4x6ULKIfAaoBG4OmFyuqo0iMhd4XUQOqerJUdZ9FHgUoLy8PIhYJtTqOvr4ydunWDUrhwPO0kjNAAANxklEQVR1XRyoswu0m+AkJbhYMyePV6qbqOvoY1ZemtORzCQFs6dQDwReVWUm8KFRsETkNuCrwD2qOnhxuqo2+n/WAtuBVaO9iKpuUtVKVa0sLLSzZ53wd1uPIgK3Lyl2OoqJQdfOzcclwpM7TjsdxUxBMEVhDzBfROaISBKwAfjAUUQisgr4ISMFoSVgeq6IJPvvFwA3AIEd1CZKvHm8lRcPN7Fu4QyyU20kVDNx2amJfHx5Cb/aU8eFgWGn45hJGrcoqKoH+CLwMnAUeEZVj4jI4yJyj3+xfwAygGcvO/R0MVAlIgeAN4BvXXbUkokCA8Nevv78YeYWpHOjdS6bKXj0prn0DHqsbyGGBXXymqpuBbZeNu1rAfdvG2O9HcDyqQQ04ff97Sc53d7HLx5ZQ11Hv9NxTAxbVpbN+qXF/OQPp3j4+gpy05OcjmQmyIa5iHOHG7r43hs13LuylBvnW1+Ombov3b6A3iEPm96udTqKmQQrCnFs0OPlL585QF56En9zz1Kn45hpYmFxJp9YUcqT75ym+cKA03HMBFlRiGP/+PIxjjV38+37V5CTZrv5JnT+8o4FeFX5xgvWhRhrrCjEqf/1m0P86O1TrJ2Tx7muAZ7affYDg5wZMxWz89P5wrqreOHgOd463up0HDMBVhTi0MnWHn79Xj2zclO5a3mJ03HMNPXYurnMLUjn/3v+MAPDXqfjmCBZUYgzLRcG+LOf7sHtEjauKSfBbW8BEx7JCW6+ed8yzrT38Te/O+J0HBMku55CHOnqH+ZPn3iXtp5BHrquwvoRTMiN1gR584JCnn63jrVz8rlv1ZWGTTPRwIpCnOjoHeLPfvouJ1t7+OnDazjb0Tf+SsaEwG2Li+gd9PC/fnuIhcWZLC7JcjqSuQJrO4gDDZ39PPCDHbzf1M33P30NH5lvZy2byHG7hO8+uJqslEQeeuJd6uwLSVSzojDNvXuqg3u/+w6t3YP84pG13LakyOlIJg4VZ6fw80fWMOjx8Sc/2W3nL0QxKwrTlKrykz+cYsOmnYDy5zfMoaalxw49NY5ZUJTJEw//Ea3dg9z//R3UtvY4HcmMworCNNTSPcDDP93DN16oZmFRJp9fdxVFWSlOxzKGa2bn8vSj19I/5OWBH+xkz+kOpyOZy4jqqNfLcVRlZaVWVVU5HSOmXPz2f/TcBX79Xj1DHh8fX17C2jl5dgU1E3Xaugf52c7TdPmv8vfw9RX2Pp0iEdmrqpVTfR47+miaGPL42HroHO+e7qAkO4VPVc5ihu0dmChVkJnM59ddxc7aNv7md9W8U9PG331yub1no4A1H00Dh+q7+O4bNbx7uoMb5xfwuZvn2YfLRL3UJDeb/qSSv7prMW+faOP277zFb/fVE42tF/HEmo9i2MCwl399/QQ/eLOW9CQ3D1wzi6tmZDgdy5gJa+0e5Nfv1XO2o4/bFhfxzfuWUZxtX2wmIlTNR1YUYtTeMx38j80HOdnaywPXzGRRcSZpSdYaaGKXT5V3atp4/f0WXCI8dvM8Hr1pLqlJbqejxYRQFQVrPoox53uH+Przh3ngBzsZGPbxsz9fwz/+8dVWEEzMc4lw4/xCtn3pZm5ZVMh3Xj3Orf+0nef2NeDzRd+X1+nK9hRiRFffML/YdZofvllLz6CHtXPz+NiSYpIT7VuUmZ5OtfXy+0ONNHYOUJyVwlfvWszHl5fgdtlRSqOJaPORiKwH/hlwAz9W1W9dNj8Z+DlwDdAOfEpVT/vnfQV4BPACf6GqL4/3elYURvh8yntnz/Pc/gZ+vbeB/mEvH100g2Vl2XbegYkLPlUO1HWy/VgrrT2DlOWk8uDacu5fPdP6HC4TsaIgIm7gOHA7UA/sATaqanXAMp8HVqjqYyKyAfikqn5KRJYATwNrgFLgVWCBql5xcPV4LAqqSvOFQWrbeqhp6aHq9Hl21bbT0j1IUoKLu1eU8NmPzGVJaZadkWzijk+Vgowkfr7zDDtOtgOwujyHjy4uYu2cPJbPzCY5Ib73miNZFK4D/lpVP+Z//BUAVf3fAcu87F9mp4gkAE1AIfDlwGUDl7vSa8ZaUfD5FJ8qXlV8Phj2+egf8tI35KVvyHPpfu+gh46+ITp6hmjvHeJ83xAdvUOcaO6hvXeQYe//+1tkJidQUZDO4pIsFhdnWjORMX5tPYMcaujiSEMXjV0jYyi5XcKcgnRm56VRkJFMQWYSBRnJZKcmkpTgIsnt+sDPxICfCS7B7RIS3P6fLpf/p3zgZ7SfXBfJk9fKgLqAx/XA2rGWUVWPiHQB+f7puy5bN6YHVP/vzx5gy4FGfKr4FLyT7ADLTE4gNz2JvPQkctISmVeYTn5G8sgbOiOJ7NTEqH8TGuOEgoxkblk4g1sWzqB30MOsvFSONF7g/aZu6s/3c7Chi47eoUl/Nsfivlgc/I9FQBC+/cAK7rm6NKSv5aRgisJo/5ku39pjLRPMuiNPIPIo8Kj/4aCIHA4im5MKgDanQwTBcoaW5QytmM957zcjnGRsC0PxJMEUhXpgVsDjmUDjGMvU+5uPsoGOINcFQFU3AZsARKQqFLtB4RQLGcFyhprlDC3LGToiEpI292DOU9gDzBeROSKSBGwAtly2zBbgIf/9B4DXdaSzYguwQUSSRWQOMB94NxTBjTHGhN64ewr+PoIvAi8zckjqE6p6REQeB6pUdQvwE+AXIlLDyB7CBv+6R0TkGaAa8ABfGO/II2OMMc4J6jRYVd0KbL1s2tcC7g8AfzzGun8L/O0Ec22a4PJOiIWMYDlDzXKGluUMnZBkjMozmo0xxjjDxj4yxhhziSNFQUT+WESOiIhPRMbs0ReR9SJyTERqROTLAdPniMhuETkhIr/yd4CHI2eeiGzzv842EckdZZlbRGR/wG1ARO7zz3tSRE4FzFvpVE7/ct6ALFsCpkfT9lwpIjv974+DIvKpgHlh3Z5jvd8C5if7t0+Nf3tVBMz7in/6MRH5WChzTTDjfxWRav+2e01EZgfMG/Xv71DOh0WkNSDPZwPmPeR/j5wQkYcuXzfCOb8TkPG4iHQGzIvI9hSRJ0SkRcY4TF9G/Iv/dzgoIqsD5k18W6pqxG/AYkaOqd0OVI6xjBs4CcwFkoADwBL/vGeADf77PwA+F6acfw982X//y8C3x1k+j5GO9jT/4yeBByKwPYPKCfSMMT1qtiewAJjvv18KnANywr09r/R+C1jm88AP/Pc3AL/y31/iXz4ZmON/HrdDGW8JeP997mLGK/39Hcr5MPDdUdbNA2r9P3P993OdynnZ8v+ZkQNtIr09bwJWA4fHmP9x4EVGzgu7Ftg9lW3pyJ6Cqh5V1WPjLLYGqFHVWlUdAn4J3CsiAtwKbPYv9zPgvjBFvdf//MG+zgPAi6raF6Y8Y5lozkuibXuq6nFVPeG/3wi0MDJkSriN+n67bJnA/JuBj/q3373AL1V1UFVPATX+54t4RlV9I+D9t4uRc4MiLZhtOZaPAdtUtUNVzwPbgPVRknMjI2O5RZSqvsXIl82x3Av8XEfsAnJEpIRJbsto7lMYbXiNMkaGz+hUVc9l08OhSFXPAfh/zhhn+Q18+E3zt/5duu/IyGiy4RBszhQRqRKRXRebuIji7Skiaxj5BncyYHK4tudY77dRl/Fvr8DhXMZbN1IZAz3CyDfIi0b7+4dDsDnv9/8tN4vIxZNcI7UtJ/Ra/ma4OcDrAZMjtT3HM9bvMaltGbYrs4jIq0DxKLO+qqrPB/MUo0yb0NAZwbhSzgk+TwmwnJHzOS76CiODAyYxcrjY/wQedzBnuao2ishc4HUROQRcGGW5aNmevwAeUlWff3LItudoLznKtJAP5zJFExk25jNAJXBzwOQP/f1V9eRo60cg5++Ap1V1UEQeY2QP7NYg1w2VibzWBmCzfvA8q0htz/GE9H0ZtqKgqrdN8SnGGiKjjZHdowT/t7Uxh84IxpVyikiziJSo6jn/P6mWKzzVfwR+q6rDAc99zn93UER+Cvw3J3P6m2NQ1VoR2Q6sAn5NlG1PEckCfg/8lX93+OJzh2x7jiIiw7lEICMichsjRfhmVR28OH2Mv384/omNm1NV2wMe/gj4dsC66y5bd3vIE/6/1wr277YB+ELghAhuz/GM9XtMaltGc/PRqMNr6EgPyhuMtN/DyPAawex5TEbg8B3jvc6H2hv9//guttvfB4RrkL9xc4pI7sXmFhEpAG4AqqNte/r/1r9lpI302cvmhXN7xsJwLuNmFJFVwA+Be1S1JWD6qH//MGQMNmdJwMN7gKP++y8Dd/jz5gJ38MG974jm9GddyEhH7c6AaZHcnuPZAvyp/yika4Eu/xeoyW3LSPSej9Jb/klGqtgg0Ay87J9eCmy9rFf9OCPV96sB0+cy8qGrAZ4FksOUMx94DTjh/5nnn17JyBXoLi5XATQArsvWfx04xMg/r38HMpzKCVzvz3LA//ORaNyewGeAYWB/wG1lJLbnaO83Rpqn7vHfT/Fvnxr/9pobsO5X/esdA+4M42dnvIyv+j9TF7fdlvH+/g7l/N/AEX+eN4BFAev+uX8b1wB/5mRO/+O/Br512XoR256MfNk85/9c1DPSV/QY8Jh/vgDf8/8Ohwg4onMy29LOaDbGGHNJNDcfGWOMiTArCsYYYy6xomCMMeYSKwrGGGMusaJgjDHmEisKxhhjLrGiYIwx5hIrCsYYYy75/wHjPosY6RXpkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "spurious_reg(0.99, 10, 0, 2000, 5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Try other values that you want\n",
    "\n",
    "Now let's use the sliders provided by widgets to experiment\n",
    "\n",
    "(Please feel free to edit the following cell in order to change the range of admissible values of $T$ and $\\rho$)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f4bc249a7bd04a13b4fabd454ccc23fc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='rho', max=0.999, step=0.01), IntSlider(value=20, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import interactive, fixed, IntSlider\n",
    "\n",
    "interactive(spurious_reg, \n",
    "            rho=(0, 0.999, 0.01),\n",
    "            sigma=fixed(10), \n",
    "            location=fixed(0), \n",
    "            sample_size_series=IntSlider(min=20, max=300, step=1, description='T'), \n",
    "            simulation=fixed(1000))"
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}