{ "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": [ "