presidential.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:227e67016d3c4d454ef1c0a33aa3827d37e643b7c3269ef4c8a2fcf8c5967a27"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Forecasting Presidential Elections (BDA 15.2) #\n",
      "\n",
      "In this notebook, we learn about the advantage of using hierarchical linear model as opposed to non-hierarchical model, by building a forecasting model for presidential elections.  \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "% pylab inline\n",
      "from pandas import *\n",
      "\n",
      "np.random.seed(13571)\n",
      "\n",
      "# define sampler for inverse chi-squared distribution\n",
      "def r_inv_chisquare(size, df, scale_sq):\n",
      "    return df * scale_sq / np.random.chisquare(df, size)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again, the data is from BDA book's homepage.  Somehow the below command does not work if engine=\"c\", which is a default value.  Very frustrating!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# I could not figure out how to do this in \"c\" engine.  how frustrating!\n",
      "presidental_df = read_csv(\"http://www.stat.columbia.edu/~gelman/book/data/presidential.asc\",\n",
      "                          skiprows=28, delimiter=r\"\\s+\", engine=\"python\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "presidental_df.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Dvote</th>\n",
        "      <th>year</th>\n",
        "      <th>state</th>\n",
        "      <th>evotes</th>\n",
        "      <th>constant</th>\n",
        "      <th>n1</th>\n",
        "      <th>n2</th>\n",
        "      <th>n3</th>\n",
        "      <th>n4</th>\n",
        "      <th>s1</th>\n",
        "      <th>...</th>\n",
        "      <th>s6</th>\n",
        "      <th>s7</th>\n",
        "      <th>s8</th>\n",
        "      <th>s9</th>\n",
        "      <th>r1</th>\n",
        "      <th>r2</th>\n",
        "      <th>r3</th>\n",
        "      <th>r4</th>\n",
        "      <th>r5</th>\n",
        "      <th>r6</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0</th>\n",
        "      <td> 0.462334</td>\n",
        "      <td> 1992</td>\n",
        "      <td> 1</td>\n",
        "      <td>  9</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0.575</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.345</td>\n",
        "      <td>-0.058</td>\n",
        "      <td>...</td>\n",
        "      <td>-0.608585</td>\n",
        "      <td>-0.5169</td>\n",
        "      <td>-0.110199</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td> 0.436794</td>\n",
        "      <td> 1992</td>\n",
        "      <td> 2</td>\n",
        "      <td>  3</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0.575</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.345</td>\n",
        "      <td>-0.084</td>\n",
        "      <td>...</td>\n",
        "      <td>-0.608585</td>\n",
        "      <td> 0.2181</td>\n",
        "      <td>-0.170000</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td> 0.489230</td>\n",
        "      <td> 1992</td>\n",
        "      <td> 3</td>\n",
        "      <td>  8</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0.575</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.345</td>\n",
        "      <td>-0.068</td>\n",
        "      <td>...</td>\n",
        "      <td> 0.691415</td>\n",
        "      <td>-0.3627</td>\n",
        "      <td>-0.170000</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td> 0.598887</td>\n",
        "      <td> 1992</td>\n",
        "      <td> 4</td>\n",
        "      <td>  6</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0.575</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.345</td>\n",
        "      <td>-0.033</td>\n",
        "      <td>...</td>\n",
        "      <td>-1.408585</td>\n",
        "      <td>-0.2121</td>\n",
        "      <td>-0.050335</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td> 0.590312</td>\n",
        "      <td> 1992</td>\n",
        "      <td> 5</td>\n",
        "      <td> 54</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0.575</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.37</td>\n",
        "      <td>-0.345</td>\n",
        "      <td> 0.021</td>\n",
        "      <td>...</td>\n",
        "      <td> 0.991415</td>\n",
        "      <td> 0.1496</td>\n",
        "      <td> 0.109646</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "      <td> 0</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "<p>5 rows \u00d7 24 columns</p>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "      Dvote  year  state  evotes  constant     n1    n2    n3     n4     s1  \\\n",
        "0  0.462334  1992      1       9         1  0.575 -0.37 -0.37 -0.345 -0.058   \n",
        "1  0.436794  1992      2       3         1  0.575 -0.37 -0.37 -0.345 -0.084   \n",
        "2  0.489230  1992      3       8         1  0.575 -0.37 -0.37 -0.345 -0.068   \n",
        "3  0.598887  1992      4       6         1  0.575 -0.37 -0.37 -0.345 -0.033   \n",
        "4  0.590312  1992      5      54         1  0.575 -0.37 -0.37 -0.345  0.021   \n",
        "\n",
        "  ...        s6      s7        s8  s9  r1  r2  r3  r4  r5  r6  \n",
        "0 ... -0.608585 -0.5169 -0.110199   0   1   0   0   0   0   0  \n",
        "1 ... -0.608585  0.2181 -0.170000   0   0   0   0   0   0   0  \n",
        "2 ...  0.691415 -0.3627 -0.170000   0   0   0   0   0   0   0  \n",
        "3 ... -1.408585 -0.2121 -0.050335   0   1   0   0   0   0   0  \n",
        "4 ...  0.991415  0.1496  0.109646   0   0   0   0   0   0   0  \n",
        "\n",
        "[5 rows x 24 columns]"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The data looks like above.  For detailed description of each column, please look at skipped rows at the top of the file and the corresponding book chapter in BDA book; here we just consider ``constant``, ``n1`` to ``r6`` columns as given covariates.\n",
      "\n",
      "Here, we want to reproduce Figure 15.1, but unfortunately we cannot infer mapping from state indices to actual state names, so I am just showing points here.  The below plot shows how the share of democratic vote in 1984 is related to that in 1988; almost linear relationship is observed."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df_1984 = presidental_df.loc[(presidental_df['year'] == 1984),('Dvote','state')]\n",
      "df_1984.rename(columns={'Dvote':'Dvote1984'}, inplace=True)\n",
      "df_1988 = presidental_df.loc[(presidental_df['year'] == 1988),('Dvote','state')]\n",
      "df_1988.rename(columns={'Dvote':'Dvote1988'}, inplace=True)\n",
      "df_1984_1988 = merge(df_1984, df_1988, on=('state'))\n",
      "fig, ax = plt.subplots()\n",
      "df_1984_1988.plot(x='Dvote1984',y='Dvote1988', kind='scatter', ax=ax)\n",
      "ax.set_xlabel('Dem vote by state (1984)')\n",
      "ax.set_ylabel('Dem vote by state (1988)')\n",
      "ax.set_aspect(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAAEPCAYAAADVt9ZXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX2cHFWV97+/JCQMisQgIpsnOrxqECUBCQi6jGsmE55n\nJRJYXUSWqAs8uqDI8BBeZIkiy4tE47sgPARdPwsukBhXmUmCJLzIO0xASVCE0cCigjGKEiWQs39U\n9UxN091TM11Vfaf7fD+f+qSrum7Vr6t7Tu4995xzZWY4juM4EeMaLcBxHCck3Cg6juMkcKPoOI6T\nwI2i4zhOAjeKjuM4CdwoOo7jJMjVKEqaK2mDpJ9LWljlnA5JD0r6iaQ1I2nrOI6TNcorTlHSeOBR\nYDbwFHAvcKyZrU+cMxm4A+gysyclvcbMnk3T1nEcJw/y7CnOAh4zs34z2wpcC8wrO+cDwA1m9iSA\nmT07graO4ziZk6dRnApsTOw/GR9LsjcwRdItku6TdPwI2jqO42TOhByvnWZcvh1wAPBuYAfgTkl3\npWzrOI6TOXkaxaeAaYn9aUQ9viQbgWfNbAuwRdKtwP7xecO1RZIbT8dxhmBmqqd9nsPn+4C9JbVL\nmgi8H1hRds73gHdIGi9pB+Bg4JGUbQEws6C2888/v+Eaxoou1+Sast6yILeeopm9KOkUoBcYD1xl\nZuslnRy/f7mZbZDUAzwEbAO+aWaPAFRqm5fWLOnv72+0hIqEqMs1pcM1FUuew2fM7CbgprJjl5ft\nXwZclqat4zhO3nhGS8YsWLCg0RIqEqIu15QO11QsuQVvF4EkG8v6HcfJFklYwBMtLcmaNWsaLaEi\nIepyTelwTcXiRtFxHCeBD58dx2kafPjsOI6TMW4UMyZUX0uIulxTOlxTsbhRdBzHSeA+Rcdxmgb3\nKTqO42SMG8WMCdXXEqIu15QO11QsbhQdx3ESuE/RcZymwX2KjuM4GeNGMWNC9bWEqMs1pcM1FYsb\nRcdxnATuU3Qcp2lwn6LjOE7GuFHMmFB9LSHqck3pcE3F4kbRcRwngfsUHcdpGtyn6DiOkzFuFDMm\nVF9LiLpcUzpcU7G4UXScFqW3t5c5c45mzpyj6e3tbbScYHCfouO0IL29vRx11Als2XIJAG1tC1m2\n7Bq6uroarKw+svApulF0nBZkzpyjWbXqSOCE+Mg1dHauYOXKGxopq258oiVAQvW1hKjLNaXDNRXL\nhEYLcByneLq7T+L2209gy5Zov61tId3d1zRWVCDkOnyWNBdYAowHrjSzS8re7wC+BzweH7rRzC6I\n3+sH/gi8BGw1s1kVru/DZ8cZJb29vSxefAUQGcmx7k+EwH2KksYDjwKzgaeAe4FjzWx94pwO4HQz\nO7JC+yeAA81sU417uFF0Wp5mNG6jJXSf4izgMTPrN7OtwLXAvArn1foAdX24RhCqryVEXa4pHbU0\nlWaRV606klWrjuSoo04oJLwmxOeUFXkaxanAxsT+k/GxJAYcKmmdpB9K2rfsvdWS7pN0Yo46Hafh\njDZmcPHiK+KwmhOAKMSm1Gt0RkeeEy1pxrUPANPM7HlJRwDLgX3i9w4zs6cl7QKskrTBzG7LS2xW\ndHR0NFpCRULU5ZoiymMGb7/9hCExg/6ciiVPo/gUMC2xP42otziAmT2XeH2TpK9JmmJmm8zs6fj4\nM5KWEQ3HX2YUFyxYQHt7OwCTJ09mxowZA19YqYvv+74f8v5gb+8NAAO9vUmTJg3bfvbsg7n99oXx\nLPJ6Jk68gu7u/wjq8+W539fXx+bNmwHo7+8nE8wsl43I4P4CaAcmAn3A9LJzdmVwsmcW0B+/3gHY\nMX79CuAOYE6Fe1ho3HLLLY2WUJEQdbmmiM7O+QZLDSzellpn5/zUmnp6eqyzc751ds63np6enNWm\n09QoYptQl+3KradoZi9KOgXoJQrJucrM1ks6OX7/cuAY4KOSXgSeB/4xbv464EZJJeP6HTNbmZdW\nx2kk9cYMdnV1tfSMc9Z4mp/jBICH1WRD0HGKReBG0XGyoVmMcuhxii1JqPFbIepyTenIW9NoYh1D\nfE5Z4bnPjtPiDI11hC1bomNjtbdYLz58dpwWp5nKiPnw2XGcuunuPom2toXANcA18ez3SVXPb/aK\n3W4UMyZUX0uIulxTOvLW1NXVxbJlUe+ws3NFzQrcg/7HvQvNtS4S9yk6jpM61nFo9k1HU/ofvaeY\nMaHmhIaoyzWlI0RN0NFoAblRs6coqQ34e+CdwN8AW4CHgR+Y2U/zl+c4Tki0QsXuqj1FSZ8myjl+\nO3AXcDnwXaJK2BdLWi3prYWoHEOE6JOCMHW5pnSEpKnkfzzwwKuG9T+OVWr1FO8BFlWJeVksaVeG\nVsFxHKcF6OrqYtKkSYEO6+vH4xQdx2kaco1TlLSnpKslfVbSjpK+Kemnkv5TUns9N3Ucp3q8X7PH\nAQZPtZpiRAVdPwqcDfwUOAN4PfAR4Ef11izLYsPrKaYmRF2trKmnp8fa2naN6ygutba2Xa2np6fi\n8UsuuaQQTSMhxO/OLP96iq80s68DSPqYmV0WH79K0qn5mGjHaQ2q5RtHr4ce/+53r+LMM8982TWa\npbJNaNQyitskvRHYCWiTdJCZ3Stpbzy+sSqhOp9D1OWa0jFlyi4vOzbcui55E+JzyopaRvFM4PvA\nNuC9wNlxCM5OgK+u5zh1UCveL00coFe2yY+qPT4zu9nM9jGzN5nZ7WY2HzgE2M3MlhcncWwRUkxZ\nkhB1tbKmavnGlY6vW7cuuImXEL+7rBguo+Vvgd+Y2aOS3kEUyP0I8IMixDlOM1Mt3zh5vLe3l/PO\nu5gXXvgCMDhMboXMkkZRNU5R0heBg4DtgB7g3cBNwOFAn5mdUZTIanicotPs1Kp16BMtLyeLOMVa\nPcVOYD+gjWgN56lm9mdJFxMtV9pwo+g4rYyv4pcPtWaRtwFGlOts8VY6XpclbmZC9bWEqMs1DU93\n90lMnPhJ0haALYrQnlOW1Oop/pAogHt74Crgu5LuIho+31qANsdpebq6urjggrNYvXoFAN3dzVeA\nITRq5j5LOpQoQvxOSXsBRwG/BK43s20FaayK+xQdx0lS2LrPkl4NvGRmf6znZlnjRtEZS/jESP7k\nXRBiqqRvSfoD8Dvgp5I2Slokabt6btrMhOprCVFXK2kazdrKeWuqhxA1ZUWtiZZ/B64GJgPHADcA\n04lCdL6avzTHaR6GZqBE6XmlXqMTFrXiFNeZ2f6J/QfM7ID49aNm9saCNFbFh8/OWKGZ1lYOmbzj\nFJ+VdDzwI+Bo4In4puPwkBzHGRGegTJ2qDV8/jBwJNALHAycEh9/NXBOmotLmitpg6SfS1pY4f0O\nSX+Q9GC8fSpt21AJ1dcSoq5W0jSStZWL0lQPIWrKiqo9RTP7JfAPFY7/Drh+uAtLGg98BZhNlBFz\nr6QVZra+7NS1ZnbkKNs6zpjBM1DGBqNao0XSv5rZZ4Y55+3A+WY2N94/C8DMLk6c0wF0m9l7Rto2\nPu4+RcdxBsg1JGcY0tRTnApsTOw/GR9LYsChktZJ+qGkfUfQ1nEcJ3OqDp8lPVejXVuKa6fpwj0A\nTDOz5yUdASwH9knRboAFCxbQ3t4OwOTJk5kxY8ZAVeCS36PI/b6+Pk477bSG3b/aftIHFIIegCVL\nljT8+yrfD/H7Kx0LRU9Iv6e+vj42b94MQH9/P5lQbfEW4FfA66q8t3G4xV+ICtL2JPbPBhYO0+YJ\nYEratvjCVakJUZdrSodrSg8ZLFxVK07xQuB7ZnZPhfcuNbOXr6Qz9JwJwKNEdRj/G7gHONYSkyWS\ndgV+a2YmaRbwXTNrT9M2bm/V9DuO03rkGqdoZufWeK+mQYzPeVHSKUQhPeOBq8xsvaST4/cvJ8qU\n+aikF4HngX+s1Tb9x3IcxxkdtXKf22s1lDRO0rRa55jZTWb2RjPby8wuio9dHhtEzOyrZrafmc0w\ns0PN7K5abccCocZvhajLNaXDNRVLrdnnyyTdKOmfJL1Z0mslvV7SuyVdANxBlAvtOE5Kent7g1uE\nyhnKcPUU3wwcBxwK7AZsAdYTLVx1vZn9pQiR1XCfojOWKF+rua1tYSZrNXtJskEKq6cYKm4UnbFE\nHkUh8jK0Y5VGBm87VQjV1xKiLteUjlqahitJltdwPcTnlBU11312HCc7iq6UU96LLK0Z3aq9yLT4\n8NlxRkDSf3f44Qewdu0DQHpfXtb+v1rD51as4Zh3PcXSTcYRTbbsbmafkfR6okyXlwV1O04zU26A\nVq36OFEZgLek7oVlXSmnVJJs0NB6T7Buhkt5Ab4BfA3YEO9PAe6rN5Umiw1P80tNiLpC0tTT02Od\nnfPtwAPfaT09PRXP6eycb7DUwOJtqcH8gdednfNz0Tba59TT02NtbbvGOpdaW9uuVT9bUZryhgzS\n/NJMtBxsZh8jCsfBzDYRrdPiOE1BclGp++8/dESLSoVMPYVtW5rhrCZwN1Gq3YPx/i6l143eCLCn\n6Iw9KvUAK/X6ynte8CqD7sx7Yc7oIYOeYprZ5y8Dy4DXSvo3onzlT9Vu4jjNR7n/7vDDz4wnWp5w\nX14zkcZyEqXznRJv0+u1xFltBNhTDNXXEqKuojWV/IadnfMHenU9PT02c+bhNm7cznGvb2FwvT7/\n7tJDET1FSd82s+OJ0vvKjznOmKBSzN65557KhRd+eeDYuHGfZI89pvKVr3ivr5UZNk5R0oNmNjOx\nPwF4yMz2rdGsEDxO0UlLpZi9KVMuYNOm82ilOL5mJ9c0P0nnxEsSvEXSc6UN+C2wop6bOk5RlNLc\n7r9/HfBwo+U4Y4CqRtHM/s3MdgQuM7MdE9sUMzurQI1jilBzQkPUlbemZKhN1CP8JnAGcA1tbQs5\n/fQP0da2ELhm4Njs2Qfnqmk0tOJ310iG9Sma2VmSXg3sDWyfOH5rnsIcp8RoU+OGFkuImDLlAg48\ncHC2+G1ve9uQbJBJkyZlrt8ZW6TxKZ4IfByYBjxItKjUnWb2d/nLq437FJufekpjtWLub6tTVOmw\nTwCzgH4zexcwE/hDPTd1nLQMVxqrFt3dJ71seNzdfVKOap1mII1R/IuZbQGQtL2ZbQDemK+ssUuo\nvpYQdeWtaTRpbq34nEZDiJqyIk1Gy5OxT3E5sErS74H+XFU5Tky9NQizrkrjND8jqqcoqQN4FdFC\n9S/kJSot7lNsDRqxBomvezI2KWSNlkrZK6FktLhRbA5CM0C+7snYpaiJlv3KbjoBOLCemzYzofpa\nQtS1Zs2aIbGEq1Yd2fCyXWvWrKlrcicvTaERoqas8IwWp6GEZoAcJ83w+eJQM1h8+Dz2CTGW0IfP\nY5eifIrvAPrM7E+SjieKU/yimf2ynhtngRvFsU+oBig0P6eTjiyMYpqahQ8TDbP3J8poOQVYW2/N\nsiw2vJ5iakLUVdJUqc5hozWFhGtKDwVV3n7RzLZJei/wVTO7UtKH0xhcSXOBJUTLGVxpZpdUOe8g\n4E7g/WZ2Q3ysH/gj8BKw1cxmpbmnM/bwWEInJNIMn28FeoAPAe8EniEaTr9lmHbjgUeB2cBTwL3A\nsWa2vsJ5q4DngasTRvEJ4ECLFsqqdg8bTr/jOK1DUSE57wf+AnzYzH4NTAU+l6LdLOAxM+s3s63A\ntcC8CuedClxPZGzLqc834DiOM0KGNYpm9rSZfd7Mbov3f2Vm30px7anAxsT+k/GxASRNJTKUXy/d\nLnlrYLWk++JKPWOCUOO3QtRVhKZSkdk5c45OFf/Yqs9ppISoKSvS+BRHS5px7RLgLDMzSWJoz/Aw\nM3ta0i5EOdcbSoY5yYIFC2hvbwdg8uTJzJgxg46ODmDwiytyv6+vr6H3H0v7fX19uV7/0ksv5bzz\nLuaFF74AwNq1x3LBBWdx5plnVm0f4vdXIhQ9Ie339fWxefNmAPr7+8mCEeU+j+jC0iHAIjObG++f\nDWxLTrZIepxBQ/gaIr/iiWa2ouxa5wN/MrPFZcfdp+hUpZ4YSA/JGZsU4lOUdKSkNL7Hcu4D9pbU\nLmkikW9yiLEzsz3MbHcz253Ir/hRM1shaQdJO8b3fwUwB19gw6mLXuAb3H//umGH0aGlHjrFknai\n5TFJl0p6U9oLm9mLRDGNvcAjwHVmtl7SyZJOHqb564DbJPUBdwP/ZWYr0967kYTqawlRV96aBovM\nngF8EPi/bNp03oCRq+RvXOO5z6kIUVNWpFmj5ThJOwHHAkslGXA18B9m9twwbW8Cbio7dnmVcz+U\neP04MGN4+Y5TnVKR2Q984F/YtOkySsPoLVvg7LMvYsOGDUPWgV62zNdocRg+o8UGs0deA3wS+CWR\noXsM+Hi90eP1bASY0eKER2fnfIOlBhZvS23KlD1fdqyzc76ZRRk2bW27xu8vtba2XRueaeOkgwwy\nWtL4FOdJWgasAbYDDjKzI4C3AqfnYagdJ0sqrdXyhje87mXn3X//OubMORpgxMsYOE3EcFaT6Jf0\nt1Xem12vVa5nI8CeYqg5oSHqKlJTeX51eW8QXmXQbbAwuJ5hq393I4Eicp/N7ARJu0maB2wD7rUo\nswUzW52XsXZam6xDYirlVy9bdg2LF1/B/fevY9OmE4HLgDVs2TKdxYuv8N5hqzKc1QT+GfgVUY/x\nW0Q+xY/Ua42z2Aiwp+gMZTQVcIr26VXyOZb8i87Yggx6imkMz8+AnRP7OwM/q/fGWWxuFMNmtMat\naCPlEyvNQxZGMU2c4rPAnxL7f4qPORUINX4ra11pcoqHi/dr5LNK6ofBiZUDD7wquImVEH9TIWrK\niqo+RUnd8cvHgLslLY/35wEP5S3MCZfyatmlGL+sDEm9az0PRzX9K1fewJo1awZya53WpGrus6RF\nDBZ1UPlrM/t07uqGwXOfG0PanOJ6lhrIM/c4xHVhnGzIIve5ak/RzBbVc2HHKWWUDBq39L3JLKtx\nlxtYx6lJvU7JRm4EONESavxWlrqympgo4llV0vrZz362qv4Qvz/XlB4KWqPFcYZQTw+waIZO9kR5\nz2vXrhgz+p3iya2eYhG4T9EZDvcfthZF1VN8o6SbJf003n+rpE/Vc1PHKYpKec/uV3RqkSZO8ZvA\nOcAL8f7DRGXEnAqEGr8Voq4iNJWG+mmLO7TqcxopIWrKijQ+xR3M7O5oCZXIiylpa76yHCc7fF1p\nZySkWff5JqJlSP/TzGZKOoYo9/mIIgTWwn2KTpKi11XxdVzCIwufYpqwlz2Bm4kWlfpv4A6gvd5p\n7yw2AgzJGSuMplBDyBSdv+z50mFCQQUhdo//fSXwquSxRm8hGsVQ47eSukL5g87yWWVVRCKtpiKL\nVoT4mwpRk1k2RjHNRMuNsfX5k5n9MT52fV3dU6ehhLYwk+OERK2CENOBfYGdJM1nMP/5VcD2xcgb\ne4RaTCBEXVlqyqqIRFpNeRetGI2mIglRU1bUKggxDzgKeA9D12t+DrjWzH6cv7za+ETL6KinUEPI\n+ESLU9REy6H1jtHz2nCfYmrKdYUw0RLis3JN6QhRk1k2PsU0cYob49X83hHv3wp8wsyerMsaOw3F\nY/ccpzJp4hRXA98B/j0+dBxwnJl15qxtWHz47DhOkiyGz2mM4joz23+4Y43AjaLjOEkKKQgB/E7S\n8ZLGS5og6YP4Gi1VCTUnNERdrikdrqlY0hjFDwPvA34NPA38A/ChPEU5TlakWWDLcZKkGT7vYmbP\njOri0lxgCTAeuNLMLqly3kHAncD7zeyGtG19+OzUollDj5zqFOVT/DnwBHAdcKOZ/T6luPHAo8Bs\n4CngXuBYM1tf4bxVRLnVV5vZDSNo60bRqYoXmG09CvEpmtnewHnAfsD9kv5L0vEprj0LeMzM+s1s\nK3At0fKo5ZxKlDb4zCjaBkeovpYQdWWpKathcrM/p6wIUVNWpFqjxczuJlr7+ULgC0RljL89TLOp\nwMbE/pPAwckTJE0lMnZ/BxzE4DKqw7Z1nBLV1nFOm4qXzEyZPfvgpk5hc4ZnWKMoaSeidL/3A3sB\ny4gM2HCkGdcuAc4yM1NUxbbU7U09Jl6wYAHt7e0ATJ48mRkzZgz8qEv/mxW9X6JR96+039HREZSe\nEsnF50fSvre3l3PPvRCAbdu2xQbxDQADBS7OOedUFi06ndWroyzV2bNPZ9KkSUPufc8997Bo0efj\n9utZu/Zi9t9/f7q6uhr+fELeD+X31NfXx+bNmwHo7+8nE4ZLeSHyJy4B3k7sg0yzAYcAPYn9s4GF\nZec8Hl//CaKc6t8AR6Zpa4Gm+RVNCOl6RVNe+kyabPAmg/kGPSMq41VkCTAnfyionuK4UV046oX+\nAmgHJgJ9wPQa518NzB9J2xCNYpE5oSOpixhirupoNVUyZHBI/O9rbOLEyan/g3j5tRYGZxSb6bvL\nmyyMYpqJlm0j73+Cmb0InAL0Ao8A15nZekknSzp5NG1Ho6OZaURdxN7eXg44oIOdd96LAw54R0Cx\nf39D9Bwu481v3j912E35an8TJ17hq/21OvVa1UZuBNhTLJKih349PT02ceIuAz3TkfbKstSR7CHD\na+Jh8+ieQSu6IJoVihg+h7y1ulEselmBasPWRgw3S4Zs5szDbeLEyQ1fWsEJgyyM4rDDZ0l7SPqC\npGWSvh9vK4Zr16oUGb81kjWNQ4wrq0dTV1cXK1fewAMPrGHFimtTr+ucp6a8cE3FkiZOcTlwJfB9\noORf9DSSQCiyLmJ390msXXs8L7xQOnIGEye+SHf3olTtyytVJ8Nj6sFrQzqZMlxXEri73u5oXhst\nOnxupA+sp6fHZs483KZM2dNmzjws9f1DWUHQaW7IYPicJvf5OKKg7ZXAXxPG9IHcLHVKWjH3eawW\nOfA8ZKcIiqqnuB9wInAxsDixORXI29cy2jCcEH1AmzaNqvhSroT4nFxTsaTxKb4P2MPMXhj2TMep\nQqU85Pe97/TGinKcCqQZPi8HTjaz3xQjKT0+fB47w2fwJUGd/CmqnuJa4K1ENQ1LPkUzsyPruXEW\ntKJRhOY1Ls36uZziKGrd5454Ozz5ut4Zniw2Apx9DjUnNERdSU2hzE6H/pxCIURNZgUFb5vZGqAf\n2C5+fQ/wYF2W2Gk4oa1d0og8bsepRJrh80lEs89TzGxPSfsAXzezdxchsBatOnyulxD9kh6y42RB\nFsPnNLPP/0K0PMBdAGb2M0mvreemTmMZ2iuDLVuiY400immrZDtO3qSJU/yrmQ0EbUuagKf5VSXU\n+K0QdSU1jSSPuyhNoeCaiiVNT3GtpHOBHSR1Ah8jyoN2xiih9so8h9kJgTQ+xfHAR4A58aFeonWY\nG95bdJ/i6PHwF6cZKSROMb7RawHM7Lf13Cxr3Cg6jpMk19xnRSyS9CzRwvSPSnpW0vnxyntOBUbi\naykyLKaarkaG5oTol3JN6QhRU1bU8il+EjgMOMjMnoCo4Czwjfi9z+cvr3mptlZxkcPYEDQ4TmhU\nHT5L6gM6zeyZsuO7AKvMbEYB+moylofPIcTlhaDBcbIk79JhE8oNIkB8LM2steM4zpijllHcOsr3\nWpq0vpbypTWjsJj8ltaspKtoDWk0NRrXlI4QNWVFrR7fWyU9V+W9tjzEtBKlYOXBsJjifXkhaHCc\n0EgVkhMqY9mn6DhO9hS1HIHjOE7L4EYxY0L1tYSoyzWlwzUVixtFx3GcBO5TdBynaQjepyhprqQN\nkn4uaWGF9+dJWifpQUn3Sjos8V6/pIfi9+7JU6fjOE6J3IxiXF3nK8BcYF/gWEnTy05bbWb7m9lM\n4MPAlYn3DOgws5lmNisvnVkTqq8lRF2uKR2uqVjy7CnOAh4zs34z2wpcC8xLnmBmf07svhLYVnYN\nLzzhOE6h5OZTlHQM0GVmJ8b7HwQONrNTy857L3AR8Frgf5vZ3fHxx4E/AC8Bl5vZNyvcw32KjuMM\nELpPMZW1MrPlZjYdeC/w2cRbh8XD6iOAf5H0zhw0Oo7jDCHPwg5PAdMS+9OAJ6udbGa3SdpD0hQz\n22RmT8fHn5G0jGg4flt5uwULFtDe3g7A5MmTmTFjBh0dHcCg36PI/b6+Pk477bSG3b/aftIHFIIe\ngCVLljT8+yrfD/H7Kx0LRU9Iv6e+vj42b94MQH9/P5lQ78LR1TYig/sLoB2YCPQB08vO2ZPBIfwB\nwMb49Q7AjvHrVwB3AHMq3KPigtiNJNRFwkPU5ZrS4ZrSE9uEumxXrnGKko4AlgDjgavM7CJJJ8fW\n7HJJZwL/RFR1Zwtwhpn9OC5me2N8mQnAd8zsogrXtzz1O44ztihsjZZQcaPoOE6S0CdaWpJQ47dC\n1OWa0uGaisWNouM4TgIfPjuO0zT48NlxHCdj3ChmTKi+lhB1uaZ0uKZicaPoOI6TwH2KjuM0De5T\ndBzHyRg3ihkTqq8lRF2uKR2uqVjcKDqO4yRwn2JG9Pb2JhaVP8kXlXecBuC5z4EYxd7eXo466gS2\nbLkEgLa2hSxbdo0bRscpGJ9oCYTFi6+IDeIJwBvYsuWSgV5jKIToA3JN6XBNxeJG0XEcJ4EPnzPA\nh8+OEwbuUwzEKIJPtDhOCLhPMSC6urpYufIGzjnn1CANYog+INeUDtdULG4UHcdxEvjw2XGcpsGH\nz47jOBnjRjFjQvW1hKjLNaXDNRWLG0XHcZwE7lN0HKdpcJ+i4zhOxrhRzJhQfS0h6nJN6XBNxeJG\n0XEcJ4H7FB3HaRrcp+g4jpMxuRpFSXMlbZD0c0kLK7w/T9I6SQ9KulfSYWnbhkqovpYQdbmmdLim\nYsnNKEoaD3wFmAvsCxwraXrZaavNbH8zmwl8GLhyBG2DpK+vr9ESKhKiLteUDtdULHn2FGcBj5lZ\nv5ltBa4F5iVPMLM/J3ZfCWxL2zZUNm/e3GgJFQlRl2tKh2sqljyN4lRgY2L/yfjYECS9V9J64L+I\neoup2zqO42RNnkYx1bSwmS03s+nAe4HP5qinEPr7+xstoSIh6nJN6XBNxZJbSI6kQ4BFZjY33j8b\n2GZml9Ro8wvgIGCfNG0leTyO4zhDqDckZ0JWQipwH7C3pHbgv4H3A8cmT5C0J/C4mZmkA4CJZrZJ\n0rBtof7IELeYAAAIeklEQVQP7ziOU05uRtHMXpR0CtALjAeuMrP1kk6O378cOBr4J0lbgS1Exq9q\n27y0Oo7jlBjTGS2O4zhZE2xGS4rA7+PiwO+HJN0h6a1p2zZIU398/EFJ9xSoqfAA+To15fKc0uhK\nnHeQpBclHT3StgVratRvqkPSH+L7PijpUyP9PAVoOi/x3siek5kFtxENmR8D2oHtgD5getk5bwd2\nil/PBe5K27ZoTfH+E8CUBjynVyRevwVYH8Bzqqgpr+c0ks8bn/cjohCxoxv9rKppavBvqgNYMdrP\nU6Sm0TynUHuKaQK/7zSzP8S7dwP/K23bBmgqkfXEUIgB8vVoKpHHBFraz3sqcD3wzCjaFqmpROG/\nqRr3bfRzqvUsUj+nUI3iSIO3PwL8cJRti9AEUdzmakn3SToxAz2pNRUcIF+PJsjnOaXSJWkq0R/b\n1xNaUrVtgKbS60b8pgw4NHaB/FDSviNoW7Sm0nupn1OeITn1kHr2R9K7iP6oSn6pvGaO6tEEcJiZ\nPS1pF2CVpA1mdlsRmsxsObBc0juJAuQ767xvnpryeE5pdS0BzjIzkyQGexeN/E1V0wSN+009AEwz\ns+clHQEsJ4otzot6NY3oOYXaU3wKmJbYn0b0v8MQ4omMbwJHmtnvR9K2YE2Y2dPxv88Ay4iGBIVo\nSmi4DdhD0pT4vIY9pyqa8npOaXUdCFwr6QmicLGvSToyZduiNTXsN2Vmz5nZ8/Hrm4DtGv2bqqFp\n5M8pSydtVhtRD/YXRI7ViVR2rL6eyPl6yEjbNkDTDsCO8etXAHcAcwrStCeDoVcHABsDeE7VNOXy\nnEbzeYGrgfmNflY1NDXyN7Vr4vubBfQ3+jnV0DTi5xTk8NnSBX7/K/Bq4OvRqIKtZjarWttGagJe\nB9wYH5sAfMfMVhakqdAA+Xo0kdNzGoGuEbVtpCYa+5s6BviopBeB54F/rNW2kZoYxXPy4G3HcZwE\nofoUHcdxGoIbRcdxnARuFB3HcRK4UXQcx0ngRtFxHCeBG0XHcZwEbhSbCEkvxeWRfiKpT9LpcWpY\no3XN0wiXqJW0SFJ3xjrOyfK8Cu2uV1QtHkkXSvqVpOfKznmDpJvjHN1b4tzm0nuXxt/dI5K+WOH6\nX0peT9L/kfTp0Wh1quNGsbl43sxmmtl+RLnERwDnN1gTwFFE63ePhDwCaM/O+LwBJL0ZGGdm/fGh\n71E5newyYKmZ7Q98Brgobn8ocChRKbX9gIMkHZ64/tuAySSei5n9AHiPpLaR6nWq40axSbEoz/Mk\n4BQASeMlfU7SPXEv5aT4eIektZKWS/qFpIsUFcu9W1Fhzj2S15U0TtITknZKHPu5pF0ktUv6UXz9\n1ZKmxX/s7wE+F/did5e0p6Sb4qolt0p6Y5WPsb+kH0v6maR/ju91jaR5iXt/p5QLnDi2W3zdByU9\nLOkdki4G2uJj347PWx5r+Ini6ilVzvtg/DwelPQNSZX+bo4jKkJQev73mNmvK5w3nag2IsAaBktg\nGbA9MAloI6ob+OvSdwdcCpzJy0tgrQH+vsrzc0ZDvXmJvoWzAc9VOPZ74LVEBvLc+Ngk4F6iXNKO\n+JxdifJKnwTOj8/7OPCFCtdcAiyIXx8MrIxffx84Pn79IWBZ/HogZzfevxnYK9H+5gr3WESU4zoJ\n2Bn4FbAb8LeJ6+4EPE7UQ0u2PR04J349DnhlpecDvDr+tw14OLH/XOKc6cAKYHy8/7XSZyy71hrg\nzcN9J8B3gI/Hr+cT1ZIs3fey+LvYDFyQaPMJ4BNVrvcB4EuN/u010xZk7rOTC3OAt0g6Jt5/FbAX\nsBW418x+AwPLzJZyQ38CvKvCta4jyvNeSpRjel18/BCi9bsB/p2od1NC8fVfSVSh/D8T7s6JFe5h\nwHIz+yvwV0m3ALPM7HuSvibpNUT5rtebWXmR2nuB/y9pu/ga6yo/Ej4hqaR3GrA3UF6u/t1ElWru\ni/W2EffgytiNykVgyzkD+IqkBcCtRBVgXpK0F/AmojqBIipx1Utk9I8BOqr4h58B/ibFfZ2UuFFs\nYuKh70tm9tv47+kUM1tVdk4H8NfEoW2J/W1U/o3cBewVG6Z5RL6xgUtWkVPyhY0DNpvZzBF8lNJ1\nS8bvW8DxRIUkFrzsRma3KarT+PfAUkmfN7NvD7lY9LnfTVTR6C+x0d2+yr2vMbPhJl+21Gif1PY0\nUUGM0n8QR5vZH+PiBndZXP5K0k1E/3mU/vN6LL7EDpJ+ZmalWoHbx/d2MsJ9ik2KooKa3wC+HB/q\nBT4maUL8/j6SdhjNtS0aty0DvgA8YoN1I3/MYHWS44h6QgDPEf1xY2Z/BJ4o9VgVMbDAV/IjAPMk\nTZK0M3A4UQ8Qoh7qabGUDRU+++uBZ8zsSuAqoGSAt5Y+f6zn97FBfBNRL5cK590MHBM/TyRNia9f\nznqinmZNJO2c8EmeHesD+CVweOz73S7+vI+Y2Q/NbDcz293MdieaTEsWdN2HaOjvZIQbxeaiNEHw\nE2AV0MNgL+5K4BHgAUkPE5W3n0DUg6s201vrveuIDN91iWOnAh+StC5+7xPx8WuB/yfpfkm7x+99\nRFIf0RB9yERJ4t4PAbcAdwKfsXjiwsx+G3+Wq6to6wD6JD0A/ANQCm+5AngonkDpASZIeoRoBvjO\nRPuB8ywqffUpYGX8uVYSlaMq5wfxfYGB8JqNRN/JRkn/Gr/1LmCDpEeBXYAL4+PXE9UMfJjIl9pn\n0exypedS/lkrneeMEi8d5ow54h7uQ8BMM3tuuPOLQNL2RAb8sAo+zrzuuStRfcDZRdyvVfCeojOm\nkDSbqJf4pVAMIoCZ/YUoJjSLhZrSMo1opt3JEO8pOo7jJPCeouM4TgI3io7jOAncKDqO4yRwo+g4\njpPAjaLjOE4CN4qO4zgJ/gdWqdv8J/CK2wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1034d0ad0>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Below, the same analysis done for 1972 & 1976 year pair.  Here it is shown that the pattern is not very linear, but it is also pointed out in the text that outlier states are in Jimmy Carter's home region, hinting us that we may have to take home regions as features of the model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df_1972 = presidental_df.loc[(presidental_df['year'] == 1972),('Dvote','state')]\n",
      "df_1972.rename(columns={'Dvote':'Dvote1972'}, inplace=True)\n",
      "df_1976 = presidental_df.loc[(presidental_df['year'] == 1976),('Dvote','state')]\n",
      "df_1976.rename(columns={'Dvote':'Dvote1976'}, inplace=True)\n",
      "df_1972_1976 = merge(df_1972, df_1976, on=('state'))\n",
      "fig, ax = plt.subplots()\n",
      "df_1972_1976.plot(x='Dvote1972',y='Dvote1976', kind='scatter', ax=ax)\n",
      "ax.set_xlabel('Dem vote by state (1972)')\n",
      "ax.set_ylabel('Dem vote by state (1976)')\n",
      "ax.set_aspect(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEPCAYAAADWANh/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXucXFWV77+/JAQ76tCEIDKZaCsPQUASuAYUHeKQToV7\nNRHCiICYRoXMHdEZp5kJyEgiM9whaITxxUMcOiN8hLmEhKgknQ6mgSgQImleSRCEzCSAIxhUvLSS\nkHX/OKc61ZXqqlOP03XO6fX9fM4ntffZe5+16nSt7L322nvLzHAcx8kio5otgOM4Tly4gXMcJ7O4\ngXMcJ7O4gXMcJ7O4gXMcJ7O4gXMcJ7PEauAkzZS0RdJTkuaXuH+RpI3h9ZikXZJao9R1HMephOKK\ng5M0GngSmA48BzwEnGVmm4co/2Hgb81serV1HcdxShFnD24q8LSZbTWzncCtwOwy5c8Gvl9jXcdx\nnL2I08BNBLYVpLeHeXshaRyQA5ZWW9dxHGco4jRw1Yx9PwKsM7Pf1FDXcRynJGNibPs5YFJBehJB\nT6wUH2fP8DRyXUluCB1nhGBmqrZOnD24DcBhktokjQXOBFYUF5K0H/DnwJ3V1gUws2G7FixYMKzP\nG+7L9UvvlWXdzGrvx8TWgzOzXZIuBLqB0cB3zWyzpHnh/evDoh8Fus2sv1LduGSNytatW5stQqy4\nfukly7rVQ5xDVMxsJbCyKO/6ovQSYEmUuo7jONXgKxmqoKOjo9kixIrrl16yrFs9xBboOxxIsjTL\n7zhONCRhCZtkyBy9vb3NFiFWXL/0kmXd6sENnOM4mcWHqI7jJB4fojqO4xThBq4Ksu7ncP3SS5Z1\nqwc3cI7jZBb3wTmOk3jcB+c4jlOEG7gqiNvP0d3dzYwZc5gxYw7d3d2xPqsUWffjZFm/LOtWD7Gu\nRXWi093dzWmnzaW/fxEA69bNZdmyJeRyuSZL5jjpxX1wCWHGjDn09MwC5oY5S2hvX8Hq1UvLVXOc\nEYH74BzHcYpwA1cFcfo5OjsvoKVlPsHOUUtoaZlPZ+cFsT2vFFn342RZvyzrVg/ug0sIuVyOZcuW\nsHjxDQB0drr/zXHqxX1wjuMkHvfBOY7jFOEGrgqy7udw/dJLlnWrh1gNnKSZkrZIekrS/CHKTJO0\nUdLjknoL8rdKejS8tz5OOR3HySax+eAkjQaeBKYTnHP6EHBW4elYklqBnwA5M9suaYKZvRTeexY4\n3sx2lHmG++AcZwSQRB/cVOBpM9tqZjuBW4HZRWXOBpaa2XaAvHEroGqFnPpp9pIxx2kUcRq4icC2\ngvT2MK+Qw4DxktZK2iDp3IJ7BqwJ88+PUc7IZN3P0dvbO7BkrKdnFj09szjttLmZMXJZfn9Z1q0e\n4oyDizJ23Ac4DjgFGAfcL+kBM3sK+ICZPS/pQKBH0hYzu6+4gY6ODtra2gBobW1l8uTJTJs2Ddjz\n0huV7uvra2h7SUv39fVx8813hOth5wK99Pd3sHjxDeRyuabL5+9v5KR7e3vp6uoCGPh910KcPrgT\ngYVmNjNMXwLsNrNFBWXmAy1mtjBM3wisMrPbi9paAPzezBYX5bsPrsH4mlgniSTRB7cBOExSm6Sx\nwJnAiqIydwIfkDRa0jjgBGCTpHGS3gwg6Y3ADOCxGGV1QpKwZMxxGkVsBs7MdgEXAt3AJuA2M9ss\naZ6keWGZLcAq4FHgQeA7ZrYJeCtwn6S+MP+HZrY6LlmjknU/R29v78CSsfb2FbS3r8jUlk1Zfn9Z\n1q0eYl2LamYrgZVFedcXpb8KfLUo7xlgcpyyOUOTy+UyY9SckY2vRXUcJ/Ek0QfnOI7TVNzAVUHW\n/RyuX3rJsm714AbOcZzM4j44x3ESj/vgHMdxinADVwVZ93O4fukly7rVgxs4x3Eyi/vgHMdJPO6D\ncxzHKcINXBVk3c/h+qWXLOtWD27gHMfJLO6Dcxwn8bgPznEcpwg3cFWQdT+H65desqxbPbiBcxwn\ns7gPznGcxFOrD67sjr6SWoAPAx8E/hToJzgb4Udm9kQtgjqO4wwXQw5RJX2Z4NT59wEPANcD/wG8\nDlwpaY2k9wyLlAkh634O1y+9ZFm3eijXg1tvZguGuLdY0kHApHKNS5oJXAOMBm4sPDKwoMw04GqC\nM1JfMrNpUes6juOUI85zUUcDTwLTgeeAh4CzzGxzQZlWgl5izsy2S5pgZi9FqRvWdx+c44wAGh4H\nJ+nC8FR5JB0q6V5Jv5H0oKRjIrQ9FXjazLaa2U7gVmB2UZmzgaVmth3AzF6qoq7jOE5ZyoWJ/LWZ\nvRh+/jrBMHJ/YD5wXYS2JwLbCtLbw7xCDgPGS1oraYOkc6uoO+xk3c/h+qWXLOtWD+V8cKMLPh9o\nZsvCz735U+crEGXsuA9wHHAKMA64X9IDEesC0NHRQVtbGwCtra1MnjyZadOmBYKGL71R6b6+voa2\nl7S06+fppKR7e3vp6uoCGPh918KQPjhJVxD0mi4HPk4QInIH8BfAHDP7cNmGpROBhWY2M0xfAuwu\nnCyQNB9oMbOFYfpGgpPut1eqG+a7D85xRgC1+uDKTjJIOg/4K+AQYF8Cw7McuNLMfltBoDEEEwWn\nAM8D69l7kuEI4JtALmz/QeBM4OeV6ob13cA5zggglsX2ZnaTmZ1gZhPM7M1mdqSZXVLJuIV1dwEX\nAt3AJuA2M9ssaZ6keWGZLQQ9tkcJjNt3zGzTUHWrVa7RZN3P4fqllyzrVg9D+uAkzQJWm9kfam3c\nzFYCK4vyri9KfxX4apS6juM41VDOB9cPvArcBXwf6Daz14dRtor4ENVxRgZxDFG3AIcC9wEXAS9I\nuk7SyTXK6DiOM6xU8sG9bGY3mNlfAMcCm4FFkraVq5dVsu7ncP3SS5Z1q4fI+8GZ2Qtm9q9mdiLB\n7iKOEwvd3d3MmDGHGTPm0N3d3WxxHNL7Tsr54D5kZmuHWZ6qcB9c9uju7ua00+bS3x+EPLa0zGfZ\nsiXkcrkmSzZyScI7iSUOrqDx/YHXzex3tQgXF27gsseMGXPo6ZkFzA1zltDevoLVq5c2U6wRTRLe\nSRyL7SdK+ndJvwV+DTwhaZukhZL2qUfYtJJ1P4frl16yrFs9lPPB3QzcBLQCZwBLgSMJ1o9+K37R\nnJFIZ+cFtLTMB5YAS2hpmU9n5wXNFmtEk+Z3Us4H94iZHVuQftjMjgs/P2lm7xomGYfEh6jZpLu7\nm8WLbwCCH5f735pPs99Jw31wku4GuoAfA3OAk81sjqRRwBYzO7wOeRuCGzjHGRnEEej7KWAWwXrQ\nEwjWhkKwJ9wXq5YwA2Tdz+H6pZcs61YPQxo4M/tPM/tLMzvazM4xsxfC/F+b2e3DJ6KTp5pYpFrj\nltIa75Rm/DuPETOr+gIuq6Veo69A/JHBqlWrrKXlIIMugy5raTnIVq1aVXfZRtRzase/82iEv/Xq\nbURNlWBbLfUafY0kA9fefnr4I7Dw6rL29tPrLtuIek7t+HcejVoNXLk4uFeGuoCDY+9aJpCs+zl2\n7HixcqEUk+X3l2Xd6qHcmQwvA1PN7JfFN0bqYvtm0tl5AevWzaW/P0gHsUhL6i5byMc+9mE2bZpf\ndT2ndmp9V040Kp3JcKeZrS9x7yoz+4e4havESAsTqSYWqda4pWbHO40k8t/1Sy/9GtjFhAkH+Xc+\nBLGuRU0qI83ANRI3ZM0lCQvY00StBq6cA7+tnPOOIMRkUi2Ov0ZdDPMkw9q1a4f1eXEx1MxdVvQb\niiTp1+jJhSTpFgc0epIB+KqkOyR9UtJRkt4i6W2STpH0T8BPCNamlrO6MyVtkfRUeERg8f1pkn4r\naWN4fang3lZJj4b5ew2TndpZvPiGsOcwFwh6EfnenONkiSEnGczsDElHAecQrGo4mOBs1M3Aj4Ar\nrMyBNJJGExwJOB14DnhI0grb+3Sse8xsVikRgGlmtqMaheIkf0BtVnH9ho9GTy4kSbckUW4WFTN7\ngtqXZU0FnjazrQCSbgVmExjIQsqNq6sfczsV8Zm75pPL5Vi2bEmBH9T9b3EQecvyGpgIFIaTbA/z\nCjHg/ZIekXSXpHcX3VsjaYOk82OUMzJZiTXK/7ja21fQ3r5iwLmdFf2GImn65XI5Vq9eyurVS+s2\nbknTLSmU7cHVSZTpzYcJJipelXQqsBzI71Jykpm9IOlAoEfSFjO7r7iBjo4O2traAGhtbWXy5MkD\n3fX8S29Uuq+vr6HtNTOdy+XYd999M6tfqXTW9ctSure3l66uLoCB33ctxBYmIulEYKGZzQzTlwC7\nzWxRmTrPAscX+90kLQB+b2aLi/ItLvkdx0kOcWyXlG94lKRzJV0Wpt8maWqEtjcAh0lqkzQWOBNY\nUdT2QZIUfp5KYHB3SBon6c1h/huBGcBjVWnmOM6IJ4oP7tvA+4Czw/Tvw7yymNkugj3kuoFNwG1m\ntlnSPEnzwmJnAI9J6gOuAT4e5r8VuC/MfxD4oZmtjqhTbGTdz+H6pZcs61YPUXxwJ5jZFEkbAcIe\nVqRDZ8xsJbCyKO/6gs/fosT5Dmb2DDA5yjMcx3GGoqIPTtKDwPuBDaGhOxBYbWZThkPAcrgPLt34\ncjEnKrX64KL04L4BLAPeIun/EAwr/7HaBzlOIcVrMdetm+trMZ2GU9EHZ2Y3A/OBfwGeB2ab2X/E\nLVgSybqfYzj1a8ZysSy/vyzrVg9RZlG/Z2abzeyb4bVZ0veGQzhn+Oju7uaiiy5LxbkAfoaBE5lK\nq/GBjUXpMcCmWlb2N/piBG1ZHifNOBfAz41wqoFGn8lAsAb1FWBX+G/+2gFcWcvDGn25gSvNqlWr\nrL39dGtvPz3Sj79Z5wJUK6dZdbLW0r6TTBpu4GyPEUmEMRtCtsZ8exFJw55btfRw9hiNtYk/+CSq\ngSv1PSxatKgJEg8PafjbrIfYDFzQNvsT7A7y5/mrloc1+nIDtze19Mb2GIP5iR/2RTXgpb6H44//\nYBMkHh7S8LdZD7UauCiTDOcD9wKrgS8TrExYWI/fL61kdc+tPbuLPDVod5EkMtROKFEYP/7Agc9Z\nm6jI6t9m3VSygMDjQAvQF6aPAJbVYk0bfeE+uL1wJ3xAue/Bv6P0QVw9OOAPZtYPIOkNZrYFeFcc\nxjbppCHWqJ4eThr0i9rzKvU95LeHyuKW7Wl4d80gykqG7ZL2J9irrUfSy8DWWKVy6iKXyzVtiBnn\n8qtqVz8Ufw9uBEYg1XT3gGnALGBsLd3FRl/4EDVRxD30q3UCpThUxIeo6YMah6gVe3DhSoZzQ2vS\nm88Dzo3H5DppZfDQD/r7g7xG9eKCA5KjU67H5+chjAyiDFGPLkxIGgMcH484yaa3tzfTs1VJ1q+7\nu5snnngEuGggb+zYv6ezc+hVg3sb3M0DBreZw/g4SPK7ayZDTjJI+qKkV4BjJL2Sv4BfUbQzb9bJ\nO7YvuuiyTIQUxEVn5wW0tMwHlgBLwtO6LmhI24sX38Brr10D3Ezw53cdRx11eCKNVNZCUFJNpTEs\nI3wlg/trqiOu5VH1BTAnf42tUx5iXKr1AeBN4edzga8Bb6/lYY2+hsPANWudpjOYehbnD+d6VP97\niYdaDVyUOLhrgVclHQv8HfAM8O+N7Uemhd5mCxArSQ6jqDW+r/Ds0XwcXBZJ8rtrJlEmGXaZ2W5J\nHwW+ZWY3SvpUlMYlzSQ4TGY0cKMVHRkoaRpwJ4HRBFhqZv8cpe5wMfgU+M20tHT5KfBNIg0TA4P/\nXgj9kP730jQqdfEI1qF+EXiK4LSr0cBjEeqNBp4G2oB9gD7gyKIy04AVtdS1YRqimvm2O2khKe8p\nKXJkCWL0wR1MMDT9YJh+G/DJCPXeB6wqSF8MXFxUZhrwg1rq2jAaOCf5NNu570YtXmo1cFHOZHjB\nzL5mZveF6f8ysyg+uInAtoL09jBvUPPA+yU9IukuSe+uou6wk3U/Rxr1y4dknH32ZyuuL41Lv3xA\ncU/PLHp6ZnHaaXOHPTwkje9uOIjig6uVKOf5PQxMMrNXJZ1KsN718Goe0tHRQVtbGwCtra1Mnjx5\nIOAx/9Ible7r62toe0lLp02/q666ii996Upee+1qgvOQNhNMBAX3d+x4cVAAbCn91q9fz5o1DwIw\nffoJTJ06tWp59gQUvx1gwLjmJzWS8n2lKd3b20tXVxfAwO+7Jmrp9kW5gBMZPMy8BJhfoc6zwPio\ndcnYEDUpw5ykyFGJwSEZqwwmVDVEbdSwtlxoSFq+y6RDjD64WcCoqhsOeoe/IJgoGEvpSYaD2HP4\n9FRga9S6ljED12wfUtLkiMLehqXTxo8/pOFnUVQyUkN9Z2n6LpNOnAbuFoIwjquAI6pqHE4FniSY\nEb0kzJsHzAs/f5ZgQ80+4KfAieXqlmg/ru+zJHFuC52EANG1a9cmQo5SNGJXkPz7y7c1fvwhFXWN\n+oxS8g3nd+lblpe+KvrgzOwcSfsBZwFdkgy4Cfi+mb1Soe5KYGVR3vUFn78FfCtqXWdksmdXkE8A\nP+Huu8/h8su/wKWXXlr1riCDdxh5B/D5gXulYtai7pCShhi9EUlUSwhMAL4A/CeB4Xka+HwtVrVR\nFz5EzawchQQ9oU6DPXKNGrV/g/xl5Ye19fTCkvhdphViHKLOBpYRDCX/AXhLmD+O0GfWrCtLBs4s\nOQ7ppMiRJzAyJzZkuFetwarXSCXtu0wrcRq4JQxxTCAwvZaHNuoabgOXdT9HUvVbtWqVjRp1QN0G\nbu3atTUZrDQYqaS+u0ZRq4GL4oObK+lgSbOB3cBDZvbL8N6a+gfJjlOeXC7H5Zd/gcsu+wK7dwd5\ntazxzMe8HXHEEcB3mDDhoEh+O/evpZd8iMbQBaTPAJcBawEBJwOXm9l34xevPJKskvxOdqjnQJvi\n7ctbWuYn+vxXZzCSMDNVXS+Cgfs58D4z+3WYPgC438yqWnEQB27gskdcp3LNmDGHnp5Z5GdDIdh6\nafXqpQ1p34mXWg1clP3gXgJ+X5D+fZg34sj6er9m6xf/ms7NDWwrWTT73SWVIX1wkjrDj08DD0pa\nHqZnA4/GLZgz8ojzVK7Ozgu4556zeO21IwEYNeoLnHxyZ4Va8Z7z6sRPuR7cm4E3ESyZWk6weN4Y\nvEHliCLrpxZlWb9cLsdll3UyalQncB27d3+KK674RtkeYhJ2CYlKlt9dXdQy9ZqUi4zFwY10Bodw\ndNqoUQfYlCknNyw0o9oYuKQuWxuJEOOZDE5I1v0czdYvf+7ClCnfYdSof2P37sVs3Hhew3pOO3a8\n2AApk0mz311ScQPnJIpcLseECQexe/fVlNu8shY+9rEPV3Vua5znvDrDQ8UwkSTjYSLZJM6Qjmon\nDXySIRnEGQf3LuDbwFvN7ChJ7wFmWXj6VTNxA5dNGhGU64YpW9Rq4KI48u8FTgA2hmkBT9Ti8Gv0\nha9FbSiLFi2qe81lo9Zt1tPOUOtNs/z+sqybWe2TDFGMyIbw340FeX21PKzRlxu4xrFq1SobO3b/\nurb2Scr2QEPNfmb5/WVZN7PaDVyUQ2delHRoQVfxDOCFqruKGSDLsUaLF98QHt5Se5BtnIG6jSDL\n7y/LutVDFAN3IXAD8C5JzxMcDHNOrFI5Th346fLOAJW6eMA7wn/fBPxJYV6zL3yI2jCyNETNy1Ls\nw8vy+8uybmbxDlHvAKaYWeGC+9uB4ytVlDQTuAYYDdxoZouGKPde4H7gTDNbGuZtBX4HvA7sNLOp\nEWR1aiSXy/FP/3Qxa9asAKKdb1CqjWrPSIgL38PNgTJhIpKOBN4NfAW4iGD21IA/Af7ezI4q27A0\nmuBUrOnAc8BDwFlmtrlEuR7gVeCmAgP3LHC8me0o8wwbSn7HcbJDrWEi5XpwhwMfAfYL/83zCnB+\nhLanAk+b2dZQwFsJdiIp3rPmcwQ9wveWaKP6uBfHcZyQIZdqmdmdZtYBfMTMziu4Pm9mP43Q9kRg\nW0F6e5g3gKSJBEbv2vxjC0UA1kjaICmKQY2drK/3c/3SS5Z1q4coPrhtkpYBHwjT9wJ/Y2bbK9SL\nMna8BrjYzEySGNxjO8nMXpB0INAjaYuZ3VfcQEdHB21tbQC0trYyefLkgSnz/EtvVLqvr6+h7SUt\n7fp5Oinp3t5eurq6AAZ+37UQZanWGoLT7W8Os84BzjGz9gr1TgQWmtnMMH0JsLtwokHSM+wxahMI\n/HDnm9mKorYWAL83s8VF+e6DSym+lMqphjiXaj0SJa9EmTEEm2W2AWOBPuDIMuVvAk4PP48D3hx+\nfiPwE2BGiTpVTzc7zac4nGTUqP1typSTEnskn9N8iHE/uF9LOlfSaEljJH2CCGcymNkugiDhbmAT\ncJuZbZY0T9K8CtXfCtwnqQ94EPihma2OIGusZN3PMVz6DV7xMJfdu69m48bX99r3rbu7mxkz5jBj\nxpyG7AeX5feXZd3qIYoP7lPAN4CvhemfAudFadzMVgIri/KuH6LseQWfnwEmR3mGkxX+lP7+vxpY\n2lW8o8i6dXP9mD+naqL44A40s0Ruheo+uHRSbLwgv6nkLwf2ffNj/pxC4jw28KeSVkv6tKT9a5DN\ncQaxZ2vym8JDYD4B/NJ3zHUaTkUDZ2aHAV8CjgZ+JumHks6NXbIEknU/x3Dql8vlePjhXu666xba\n25+lvX3FoCFoHNuFZ/n9ZVm3eojig8PMHiQ4G/UK4GqCv7rvxSmYMzIYas1okta1Ouklig9uP+A0\n4EzgUGAZwYzoz+IXrzzug3OckUGcZzI8S3DY823AA0myKG7gHGdkEOckwyFm9rdmdv9ItyZZ93O4\nfukly7rVQ5RJht3DIYjjOE6j8XNRHcdJPHEOUR3HcVJJRQMn6Z2Srpa0TNIPwmtFpXpZJOt+Dtcv\nvWRZt3qIEge3HLgR+AGQ98f5uNBxnMQTJUzkQTM7YZjkqQr3wSUH39/NiZM44+DOIQjwXQ38MZ9v\nZg9X+7BG4wYuGRQvnm9pme87fzgNJc5JhqMJDpm5ElhccI04su7nqFW/4v3d+vsXDfTmkkSW31+W\ndauHKD64jwHvNLPX4hbGcRynkUQZoi4H5pnZfw+PSNHxIWoy8CGqEzdx+uDuAd5DcHBz3gdnZjar\naikbjBu45OCTDE6cxGngpoUfjT0nYJmZ3VPtwxrNcBu43t7egSPOsojrl16yrBvEOMlgZr3AVmCf\n8PN6YGNEoWZK2iLpKUnzy5R7r6RdkuZUW9dpPI0+7MVxmkWUHtwFBLOo483sEEmHA9ea2SkV6o0G\nngSmA88RDHHPMrPNJcr1EJyJepOZLa2irg9RG4z705wkEmeYyGcJTrX/HYCZ/Rx4S4R6U4GnzWyr\nme0EbgVmlyj3OeB24MUa6joNJi0hH44ThSgG7o9mNhDgK2kM0ZZqTQS2FaS3h3kDSJpIYLiuDbPy\n7Vas2wyyHmvk+qWXLOtWD1Hi4O6RdCkwTlI78NcE61IrEcUIXgNcbGYmSRRMYkSoC0BHRwdtbW0A\ntLa2Mnny5AFna/6lNyrd19fX0Paame7u7ubSS68A4IorLiWXy9HX18f06Sewbt18+vsBNjN27A10\ndn6/6fI2Ip2l95f1dG9vL11dXQADv+9aiOKDGw18GpgRZnUDN1Zyfkk6EVhoZjPD9CXAbjNbVFDm\nGfYYtQkEfrjzgV9Vqhvmuw+uBir52Tzkw0kasYWJhI2/BcDMflWFQGMIJgpOAZ4nmH3da6KgoPxN\nwA/M7I6odd3A1YYfquykjYZPMihgoaSXCIzNk5JekrQgHE6Wxcx2ARcS9Pg2EZzEtVnSPEnzaqkb\nXa14yLqfw/VLL1nWrR7K+eC+AJwEvNfMnoVg80vguvDe1yo1bmYrgZVFedcPUfa8SnWdxtDZeQHr\n1s0N/WyEhyovaa5QjhMDQw5RJfUB7Wb2YlH+gUCPmU0eBvnK4kPU2nE/m5MmGu6Dk/S4mR1d7b3h\nxA1cdNygOWkmjkDfnTXeyyxp9XPkZ017embR0zOL006bW3IJVlr1i0qW9cuybvVQzsC9R9IrpS7g\nmOES0KmfLKxO8PWxTi0MOclgZqOHU5A0kOXdGiC5+hXH7a1bN7em9bFJ1a8RZFm3eoiyksFJOWmf\nNR3cA4X+/iDP/YhOJfzg5ypIq58jl8uxbFkQzNvevmLI3k9a9YtKlvXLsm714D24EUIul0ttjyft\nPVCneURaqpVUPExk5OBhLiObWNeiJhU3cI4zMohzw0snJOt+DtcvvWRZt3pwA+c4TmbxIarjOInH\nh6iO4zhFuIGrgqz7OVy/9JJl3erBDZzjOJnFfXCO4yQe98E5juMUEauBkzRT0hZJT0maX+L+bEmP\nSNoo6SFJJxXc2yrp0fDe+jjljErW/RyuX3rJsm71ENta1PC4wW8C04HngIckrSg6PGaNmd0Zlj8G\n+A/gyPCeAdPMbEdcMjqOk21i88FJeh+woOBs04sBzOzKMuVvNLOjwvSzwP8ws1+XeYb74IbA1246\nWSKJPriJwLaC9PYwbxCSPippM/BD4FMFtwxYI2mDpPNjlDNzRN2i3HGyTpwGLlLXysyWm9mRwEeB\nfy64dZKZTQFOBT4r6YMxyFgVafFz1LpFeVr0q5Us65dl3eohzv3gngMmFaQnEfTiSmJm90l6p6Tx\nZrbDzF4I81+UtAyYCtxXXK+jo4O2tjYAWltbmTx58sD2zfmX3qh0X19fQ9uLK72Hwems6Jf19+fp\nafT29tLV1QUw8PuuhTh9cGOAJ4FTgOeB9cBZhZMMkg4BnjEzk3QccKeZTZI0DhhtZq9IeiOwGviy\nma0ueob74EpQfIZBS8v8ms4wcJykUKsPLrYenJntknQh0A2MBr5rZpslzQvvXw/MAT4paSfQD5wZ\nVn8rcIekvIy3FBs3Z2jyW5TvmWRw4+aMTHwlQxX09vZm+vQi1y+9ZFk3SOYsquM4TlPxHpzjOInH\ne3CO4zhFuIGrgqzHGrl+6SXLutWDGzjHcTKL++Acx0k87oNzHMcpwg1cFWTdz+H6pZcs61YPbuAc\nx8ks7oNzHCfxuA/OcRynCDdwVZB1P4frl16yrFs9uIFzHCezuA/OcZzE4z44x3GcItzAVUHW/Ryu\nX3rJsm7adCo1AAAJLklEQVT14AbOcZzM4j44x3ESj/vgHMdxiojVwEmaKWmLpKckzS9xf7akRyRt\nlPSQpJOi1m0GWfdzuH7pJcu61UNsBk7SaOCbwEzg3cBZko4sKrbGzI4ND3j+FHBjFXWHnfy5mlnF\n9UsvWdatHuLswU0FnjazrWa2E7gVmF1YwMz+X0HyTcDuqHWbwW9+85tmixArrl96ybJu9RCngZsI\nbCtIbw/zBiHpo5I2Az8k6MVFrus4jlOOOA1cpOlNM1tuZkcCHwX+OUZ56mbr1q3NFiFWXL/0kmXd\n6iG2MBFJJwILzWxmmL4E2G1mi8rU+QXwXuDwKHUleYyI44wQagkTGROHICEbgMMktQHPA2cCZxUW\nkHQI8IyZmaTjgLFmtkNSxbpQm8KO44wcYjNwZrZL0oVANzAa+K6ZbZY0L7x/PTAH+KSknUA/gSEb\nsm5csjqOk01SvZLBcRynHL6SoQQRApTPCQOUH5X0E0nvaYactVJPAHbSiRogLum9knZJOn045auX\nCO9umqTfhu9uo6R/bIactRLl/YU6bpT0uKTesg2amV8FF8GQ+GmgDdgH6AOOLCrzPmC/8PNM4IFm\ny91g/d5Y8PkYYHOz5W6UbgXlfkwQmjSn2XI3+N1NA1Y0W9YY9WsFngD+LExPKNem9+D2JkqA8v1m\n9tsw+SDwZ8MsYz3UE4CddKIGiH8OuB14cTiFawBR9Uvr5FsU/c4GlprZdgAze6lcg27g9qbaIONP\nA3fFKlFjqScAO+lU1E3SRIIfzbVhVpqc0FHenQHvD10Md0l697BJVz9R9DsMGC9praQNks4t12Cc\nYSJpJfIfvKQPEfz4U+OjoooAbGC5pA8SBGC3xypVY4ii2zXAxWZmkkS6ejtR9HsYmGRmr0o6FVhO\nEFeaBqLotw9wHHAKMA64X9IDZvZUqcJu4PbmOWBSQXoSwf8kgwgnFr4DzDSzl4dJtkYQSb88Znaf\npHdKGm9mO2KXrj6i6HY8cGtg25gAnCppp5mtGB4R66Kifmb2SsHnlZK+nZJ3B9He3zbgJTPrB/ol\n3QscC5Q0cE13LCbtIjD6vyBwdI6ltKPzbQTO0BObLW9M+h3CnhCi44BtzZa7UboVlb8JOL3Zcjf4\n3R1U8O6mAlubLXeD9TsCWEMwITEOeAx491Bteg+uCIsWoHwZsD9wbdgT2GlmU5slczVE1K9kAHbS\niahbaomo3xnA/5a0C3gV+HjTBK6SKPqZ2RZJq4BHCSa/vmNmm4Zq0wN9HcfJLD6L6jhOZnED5zhO\nZnED5zhOZnED5zhOZnED5zhOZnED5zhOZnEDl1IkvV6wZUyfpL8Llx41W67Z1R7xKGmhpM4Gy/HF\nRpYrUe/2cMdpJF0h6b8kvVJU5u2S7g7Xha4N18Ei6UMF2xltlNQvaVZ475Zwu6DHJH1X0pgw/39J\n+nItso5k3MCll1fNbIqZHU2wTvRUYEGTZQI4jeAs22qIIxjzkgaXG0DSUcAoM9saZt1JsGqgmK8C\nXWZ2LHA58C8AZrY2fHdTgL8gCMhdHda52cyOMLNjgBbgM2GdHwEfkdRSrbwjGTdwGcDMXgQuAC6E\n4OBsSV+RtD7sPVwQ5k+TdI+k5ZJ+Ielfws07Hww373xnYbuSRkl6VtJ+BXlPSTpQUpukH4ftr5E0\nSdL7gY8AXwl7Ju+QdIikleHOD/dKetcQahwr6aeSfi7pM+GzlkiaXfDsW/I9nYK8g8N2N4a9ng9I\nuhJoCfO+F5ZbHsrwuKTzw7xS5T4Rfh8bJV0nqdRv5ByCRez573+9mf2yRLkjCfadA+il9NZGfwnc\nZWZ/CNtaWXDvIQZvxdULfLhEG85QNHv9mV81r9t7pUTey8BbCIzdpWHevgQ/lDaCzRBfJlivOJZg\nIfOCsNzngatLtHkN0BF+PgFYHX7+AXBu+Pk8YFn4edD6TuBu4NCC+neXeMZCgnWH+wIHAP8FHAz8\neUG7+wHPEPScCuv+HfDF8PMo4E2lvh9g//DfFoL1i/sXlyMwSCuA0WH623kdi9rqBY6q9E6AW4DP\nh59PJ1hatH9RmR8D/7NEW/sAPwNOKsg7G/h6s//20nT5WtRsMgM4RtIZYfpPgEOBncBDZvbfMHBM\nY35o9DjwoRJt3Uaw9raLYF3jbWH+iQRn2QLcDFxVUEdh+28i2P34/xa4B8eWeIYBy83sj8AfJa0F\npprZneFuGBMI1ljebmbFm28+BPybpH3CNh4p/ZXwN5Ly8k4i2FdsfVGZUwh2G9kQytsClOqZHUy0\nzTIvAr4pqQO4l2C3jNfzNyUdDBxNsPaymG8D95jZTwryXgT+NMJznRA3cBkhHF6+bma/Cn+cF5pZ\nT1GZacAfC7J2F6R3U/rv4QHg0NDIzCbwJQ00OYQ4eZ/aKOA3FviaqkHs2UX434FzCRb8d+z1oGA7\npw8SDN26JH3NzL43qLFA71MIdn/5Q2hA3zDEs5eYWaWJh/4y9Qtle4Fg44K8sZ9jZr8rKPIx4A4z\ne72wnqQFwAFmdn5Rk28In+1ExH1wGUDSgcB1wDfCrG7grwtm4A6XNK6Wti0YGy0DrgY22Z69737K\nnp0qziHooQC8QtBjJPwxP5vvSSqg1AE9AmZL2lfSAcDJBD0zCHqOfxuKsqWE7m8DXjSzG4HvAnlj\nujOvfyjPy6FxO4Kg90mJcncDZ4TfJ5LGh+0Xs5mgB1gWSQcU+PAuCeUr5Czg+0V1PkPQAz+7RJOH\nEwyvnYi4gUsveef440APsIo9vasbgU3Aw5IeI9ieewxBz2qoGcty924jMGK3FeR9DjhP0iPhvb8J\n828F/l7SzyS9I7z3aUl9BMPgQZMEBc9+FFgL3A9cbqHT3sx+Fepy0xCyTQP6JD1M4LD/1zD/BuDR\ncPJgFTBG0iaCmcz7C+oPlLPg7N1/BFaHeq0G3lrimT8KnwuApKskbSN4J9skXRbe+hCwRdKTwIHA\nFQV12oCJZnZPUdvXEvhR79fep2JNC5/tRMS3S3ISTdjzfBSYYgW71TYTSW8gMMYnlfAJxvXMg4Bb\nzGz6cDwvK3gPzkkskqYT9N6+nhTjBmBBSMcCyh9G1GgmEcwYO1XgPTjHcTKL9+Acx8ksbuAcx8ks\nbuAcx8ksbuAcx8ksbuAcx8ksbuAcx8ks/x+FkpDP6189QgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1022e5590>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let us do the actual regression analysis to form a model.  We confine ourselves to years before 1992, and drop every row that contains missing values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df_before1988 = presidental_df[presidental_df['year'] <= 1988].dropna(axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I generally do not like one-letter variables, but when doing regression analysis it is very tempting...!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# let us predict percentages instead of fractions\n",
      "X = matrix(df_before1988.iloc[:,4:].as_matrix())\n",
      "y = matrix(df_before1988.loc[:,('Dvote')].as_matrix()).T * 100"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For linear regression, posterior distribution can be simply described by MLE."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def linear_regression(X,y):\n",
      "    XtX = X.T * X\n",
      "    Xty = X.T * y\n",
      "    Q,R = numpy.linalg.qr(X)\n",
      "    R_inv = np.linalg.inv(R)\n",
      "    beta_hat = np.linalg.solve(R, Q.T * y)\n",
      "    n = X.shape[0]\n",
      "    k = X.shape[1]\n",
      "    s_sq = (np.asarray(y - X * beta_hat) ** 2).sum() / (n-k)\n",
      "    return beta_hat, s_sq\n",
      "\n",
      "beta_hat, s_sq = linear_regression(X,y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we will be using national RMSE (root mean squared error) as test statistic for model checking.  Please refer to the book for the precise definition."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rmse_nationwide(X, y, beta):\n",
      "    return sqrt(np.mean(Series(np.squeeze(np.asarray(y - X * beta)), \n",
      "                               np.asarray(df_before1988['year'])).groupby(level=0).mean() ** 2))\n",
      "rmse_nationwide(X, y, beta_hat)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "1.0091015470474629"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we sample 200 $\\beta, \\sigma^2$ parameters from the posterior distribution, generate new dependent variable $y^{\\text{rep}}$, and compare how the test statistic changes as $y$ is switched to $y^{\\text{rep}}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sample_num = 200\n",
      "n = X.shape[0]; k = X.shape[1]\n",
      "Q,R = numpy.linalg.qr(X)\n",
      "R_inv = np.linalg.inv(R)\n",
      "\n",
      "original_rmses = []\n",
      "replicated_rmses = []\n",
      "\n",
      "for sample_index in range(sample_num):\n",
      "    sigma_sq = r_inv_chisquare(1, n-k, s_sq)[0]\n",
      "    sampled_beta = beta_hat + R_inv * np.asmatrix(np.random.normal(size=k)).T * sqrt(sigma_sq)\n",
      "    original_rmses.append(rmse_nationwide(X, y, sampled_beta))\n",
      "    sampled_y = X * sampled_beta + sqrt(sigma_sq) * np.asmatrix(np.random.normal(size=n)).T\n",
      "    replicated_rmses.append(rmse_nationwide(X, sampled_y, sampled_beta))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.scatter(original_rmses, replicated_rmses, marker='.', s=1)\n",
      "plt.xlim(0,1.5)\n",
      "plt.ylim(0,1.5)\n",
      "plt.plot([0,3], [0,3], color='k', linestyle='-', linewidth=1)\n",
      "plt.axes().set_aspect(1)\n",
      "plt.xlabel(r'$T(y,\\theta)$')\n",
      "plt.ylabel(r'$T(y^{rep},\\theta)$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "<matplotlib.text.Text at 0x1061f7590>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAEPCAYAAACDYTshAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGlNJREFUeJzt3X+wnVV97/H3pxilVAGZjLQTxAiktvaK1Q4xiCQHcaZH\nxh9XamrxR+VHwekMkT86Ldppm+P13nGwc3vR0EtpJMh1vGRGsN7QQbzU3p1AGihYIbQGm0OgQLjD\nVEu53mq9ofneP579nDzZ2T+effazn/38+LxmMux99jp7raP7fM5a61nPWooIzMyK8BOzboCZNYcD\nxcwK40Axs8I4UMysMA4UMyvMS2bdgOWS5MtTZjMSEer39Vr3UCKi9H+bN292vQ2s0/WO/nfddddx\n1llnDf2drHWgmFk5PvvZz7J161Y6nc7Qcg4UMxsqGyarVq0aWtaBMqa5uTnX28A6XW9/44QJgCLq\nObcpKeradrM6GBQmkogmTsqa2XSM2zNJOVDM7CjLDRNwoJhZxiRhAg4UM+uaNEzAgWJmFBMmMINA\nkbRN0nOSHh1R7hxJL0q6uKy2mbVRUWECs+mh3ALMDysg6TjgOuBuoO/lKTObXJFhAjMIlIi4F3h+\nRLFNwO3AP06/RWbtVHSYQAXnUCStAt4L3Nj9klevmRVsGmEC1dy+4HrgExERksSQIc/CwsLS47m5\nuZktYTark3HDpNPpjLwpMDWTpfeSVgN3RsQb+rx2gCMhshL4IXBlROzoKeel92ZjKqJnMmzpfeV6\nKBFxRvpY0i0kwbNjyLeYWQ7TGuZklR4okm4DNgArJT0NbAZWAETETWW3x6wNyggT8N3GZo1XdJj4\nbmOzliqrZ5JyoJg1VNlhAg4Us0aaRZiAA8WscWYVJuBAMWuUWYYJOFDMGmPWYQIOFLNGqEKYgAPF\nrPaqEibgQDGrtSqFCThQzGqramECDhSzWqpimIADxax2qhom4EAxq5Uqhwk4UMxqo+phAg4Us1qo\nQ5iAA8Ws8uoSJuBAMau0OoUJOFDMKqtuYQIOFLNKqmOYgAPFrHLqGibgQDGrlDqHCcwgUCRtk/Sc\npEcHvP4hSY9I2itpt6Szy26j2SzUPUxgNj2UW4D5Ia8fANZHxNnAp4E/LaVVZjPUhDCBGQRKRNwL\nPD/k9T0R8UL36QPAaaU0zGxGmhImUP05lCuAu2bdCLNpaVKYQAXPNk5JugC4HDhvUJmFhYWlx3Nz\nc8zNzU29XWZFqUuYdDodOp1OrrIzOYpU0mqSQ9DfMOD1s4GvAvMRsTigjI8itdqqS5j0U6ujSCWd\nThImHx4UJmZ1VucwGaX0Hoqk24ANwErgOWAzsAIgIm6S9AXgfcBT3W85FBFr+7yPeyhWO00Ik2E9\nlJkMeYrgQLG6aUKYQM2GPGZN1JQwGcWBYjZlbQkTcKCYTVWbwgQcKGZT07YwAQeK2VS0MUzAgWJW\nuLaGCThQzArV5jABB4pZYdoeJuBAMSuEwyThQDGbkMPkCAeK2QQcJkdzoJgtk8PkWA4Us2VwmPTn\nQDEbk8NkMAeK2RgcJsM5UMxycpiM5kAxy8Fhko8DxWwEh0l+DhSzIRwm43GgmA3gMBmfA8WsD4fJ\n8pQeKJK2SXpO0qNDynxe0n5Jj0h6U5ntM3OYLN8seii3APODXpR0EXBWRKwBrgJuLKthZg6TyZQe\nKBFxL/D8kCLvAW7tln0AOFnSqWW0zdrNYTK5Ks6hrAKezjx/BjhtRm2xlnCYFOMls27AAL2nkvU9\nInBhYWHp8dzcHHNzc9NrkTWWw2S4TqdDp9PJVXYmR5FKWg3cGRFv6PPanwCdiNjeff4YsCEinusp\n56NIbWIOk/HV7SjSHcCvA0haB/xzb5iYFcFhUrzShzySbgM2ACslPQ1sBlYARMRNEXGXpIskLQL/\nAlxWdhut+Rwm0zGTIU8RPOSx5XKYTKZuQx6zqXGYTJcDxVrDYTJ9DhRrBYdJORwo1ngOk/I4UKzR\nHCblcqBYYzlMyudAsUZymMyGA8Uax2EyOw4UaxSHyWw5UKwxHCaz50CxRnCYVIMDxWrPYVIdDhSr\nNYdJtThQrLYcJtXjQLFacphUU64NliStAN4PnEuy3+sJwGHgh8Be4MsR8a/TaqRZlsOkukZusCTp\nHGA9cE9E7O3z+lnARcDeiOhMo5ED2uUNllrIYTJ7wzZYyhMob4iIgaf8ZcqdARyMiB8vr5njcaC0\nj8OkGiYNlDNJhjsnkZyXsyciHi68lWNyoLSLw6Q6Jt0C8u3AHcA3gQuB/yjpW5I+VGAbzQZymNRH\nnkD5CeDlEfFNkrN03gW8FTgs6TeXU6mkeUmPdQ9Ev7bP6ydJulPSw5L+VtKly6nH6s9hUi95hjwC\nrgHeBfwAuBk4APwYuDgi/nCsCqXjgO8C7wAOAg8Cl0TEvkyZ3wVeERGflLSyW/7UiHgxU8ZDnoZz\nmFTTsCHPyMvG3d/a6yX9Mcnw51zg3cD3gduW0Z61wGJEPNlt3HbgvcC+TJnDwIndxycC38+GiTWf\nw6Sech/0FRGHgG90/02i32Hob+kpcwNwp6RngVcAvzphnVYjDpP6WvbJgd3zie8BrgJ+CviLnIvb\n8oxT5oG/iYgLuleZ7pH0xoj4QbaQD0tvHodJ9ZR2WLqkVRFxcMzvWQcsRMR89/kngcMRcV2mzJ8D\nn4mI3d3n3wSujYiHMmU8h9IwDpN6mNrJgeOGSddDwBpJqyW9FPgAyQHpWU+RTNoi6VTgdSQTwdZQ\nDpNmGDnkkTQP/DTw3yPi/01aYUS8KOlqkrmY44CbI2KfpI91X78J+DTwRUl7Se4d+p2I+KdJ67Zq\nqmuYbNmyFYBNm66ccUuqI88cygvAk8BHga1FVBoRXwe+3vO1mzKP/zfwy0XUZdVW1zBJedh9tDzr\nUP4D8Dzwf4A/q0pPwXMo9Vf3MMmraT2ZSedQdgBfJtmq4PNFNszaqy1hkmrLH7+xrvJIOr4q+564\nh1JfbQuTpplopayki4DHSe7n+VbRjbN2qUuYjDNMadqQZhJ5hjyvAl5OshzebNnqEiapcXrA7i0n\n8kzK/hxwLfCXEfGlUlqVg4c89VK3MJlE03ssEw15gHdGxGUFt8lapA1h0hsiO3fed9TztsgTKKdJ\n+ijJvTrLWRlrLdakMNmyZSu7du1m/frz+gZF2mNOX2tjDzrPkGc9ycK28yJiOdsVTIWHPNXXpDCB\nJFB27ryP9evPI9kmqH09EJhwHUpE7IqIpyLiNknbJd0s6XJJryu+qdYUTQmTLVu2Lg1nADZseBsf\n//hVQDt7IKOMtX1BRPyapDXAOuAaSRcCXwN+v4j7fKwZ6hIm2aCAwb2NbHBkH6e9FDti3IVt67rf\ns6f7fCPwCPDuiPjP02niwLZ4yFNBdQkTSAIlItIu/FLPI8/3bdv2Jc4887XcfvutU25l9Ux0jEbP\nG/0ecAh4M8lS/KeADsmitzsnb2p+DpTqqVOYjKP3Ck7b51ImvWyc9TWS8MhuhvQbwBMTtM8aoKlh\nkspeBk7/wZFejiWG9lAkvYxk9/nvjXwj6fSIeKrIxo2ozz2UimhCmIxajJYGx7BhUdMXtKWWfZWn\ne6zoOkkflPSTA978lZKuAl4zeVOtbpoQJqlBf6DSoMgzx9L2P3J51qF8C/gcsBo4BTgeWAH8G8k8\nyjPA1oh4YaotPbZd7qHMWJPCBAb3MIb1Tob1SpraY5l0DuVPI+K/Fdwmq7mmhUmq35L5UYEw7A9b\n2/7o5bnbeO2gRWySXlVwe6wGmhgmaW9iw4a3jRUCmzZdiaRj1rSk2rZWJU+gvAb4nKQDku6R9J8k\nvVvSTwO/MuX2WcU0MUxSae/k4x+/6pgVssPs2rV76Xt7ta2HkmfI85V0A+nuKtm1wIXAJ4CzgRun\n1zyrkiaHSb8b+nrDYNCcyPr15/UNjuyl5X7f10R5AuXfS9oWEYciYj+wn2SPWST99rgVdo/luJ7k\nCI0vZNe0ZMrMAf+FZPL3exExN249Vqwmh0kqGwDpXcW9YTAsOFL9AqQtPZU8gXIN8D5J+yPi2z2v\n/c9xKpN0HMm5xe8ADgIPStoREfsyZU4G/hj45Yh4RtLKceqw4rUhTHrt33+AiDhqTiUbEKN6Hen3\njHPJuQlGBkpE/D3w9wNee2TM+tYCixHxJICk7cB7gX2ZMh8E7oiIZ7p1jFxUZ9PTxjDJ9kaGhcig\nXkdvyLSldwITnm08dmXS+0l6Hld2n38YeEtEbMqUSYc6vwC8Avhcv60nvQ5l+toUJnnmObI3E44q\n22RF3sszqTwJsILk5sMLgROAPZLu787fHGVhYWHp8dzcHHNzc8W00loVJqmI6Bss/W4ObNMfs06n\nQ6fTyVW27B7KOmAhIua7zz8JHO652fBa4CcjYqH7/AvA3RFxe897uYcyJW0Mk1R2VWz2snGe7Q3a\ncjWnSj2Uh4A1klYDzwIfAC7pKfM/gBu6E7gvA94C/FGJbWy1NoZJNgj6zX+kQbJx46UAfOUrXxz4\nXnkvNTdVqYESES9Kuhr4Bsll45sjYp+kj3VfvykiHpN0N7CX5CygrRHxnTLb2VZtDJNU71WZ3nDZ\nsmUri4sHOOOM1w58jzw7vjVdqUOeInnIU6w2h0nWoBsBN268lIhY2qGt37zKsOdNUqUhj1WQw6R/\nAGS/1m81bLZXs2vXbs4//61Hzbu08WweB0rLOUyODoRe/Ra1ZZ+n35u9nJzt4Qy6ctRUHvK0mMMk\nMWy/k34TsdleSLpZdbrRdb9Nq/Ps9lYnHvLYMRwmRwzb9jGdiO23SlYSZ5752qXNqvOunG0y91Ba\nyGEyes4kfZ5dk5Id1qSGhdGw1+tsopMDrVkcJkek8xu9C9iyssvss/MivXug9Ns/pd8fvHH2Wakj\nD3laxGFyRNpzSC8H79q1G0jmStKJVjg2FDZtunIpELLbG/SWHdYzaXLP2oHSEg6T/tLLwffe+1cc\nPnz05krp3MjOnffx+ONPkM2B889/6zFDoexWkP3Wo6SP03JNHA45UFrAYTJY+kvdewVmw4a39f2F\nz86pZEOh91jT3u/JPk7LNHGexYHScA6T5cmGwK5du3n88Se4/PKPHPXL33uaYL/NlIbtnt/Eu5Yd\nKA3mMMlv0D08Gzdeyp49DxIR7Nx531FzLXDsfMi4q2Obtiu+A6WhHCbj69dbWFw8wMqVp3DFFb9+\nzFzLoNWz2ffJu1VkU3gdSgM5TPLLc5Nf7yrXfhOtvY+zZZu0Sha8DqVVHCbj6/3DtG3bl5aW3APc\ne+9f9V1jkr2PZ9B7ZdevtIF7KA3iMMlv2FBk48ZLOXw4uOOOW9myZevSvEjvlZ8m9j7yGNZDcaA0\nhMNkPOOGQRoskli//rxW7HsyiG8ObDiHyfgGHc7V7/Xsucc33/ylY47X6N0XJRs4beNAqTmHSXEG\nLUxLX8vey5PqDY70gLBBrzedhzw15jCZnUFXdNKNmtJgauL8ioc8DeQwKU6erQx6vw7915Ck9/iM\nep+mcqDUkMOkeIOGOf3s3HkfGza8bSk4hoVP9rJyG0Kl9ECRNA9cT3KMxheyh3z1lDsH2AP8akR8\ntcQmVprDpHj9ftEHXcVJ/9u7xL7fMRzAUcOfXk3svZQaKN3Du24A3gEcBB6UtCMi9vUpdx1wN9Cs\nmx0m4DCZrkG/4P0Wq8Hg+3Z6eyRbtmwdeM9O0+YBy+6hrAUWI+JJAEnbgfcC+3rKbQJuB84ptXUV\n5jBZvlHzIemy+X473w/qPQwKlXHu2WlSzyRVdqCsAp7OPH+G5KjRJZJWkYTM20kCpVkRvgwOk8kN\n6glkv947oTpKdriznIldB8rk8oTD9cAnIiKU9BMHDnkWFhaWHs/NzTE3Nzdp+yrHYTK5Ub2MfmXG\nORlw48ZLWVw8wGWXffiYOvIEWa+qBU6n06HT6eQqW3agHARenXn+apJeStYvAdu7Y86VwDslHYqI\nHb1vlg2UJnKYzFZ2ONN7EyAcCYXFxQN915yMCrJhwVGluZXeP9af+tSnBpYtO1AeAtZIWg08C3wA\nuCRbICLOSB9LugW4s1+YNJ3DZHayS+3TX+zsUaTZPWQBLr/8I8sOgKbNrZQaKBHxoqSrgW+QXDa+\nOSL2SfpY9/WbymxPVTlMZm/Y3rHZI0dTo3Ze69cbqXNwDFL6OpSI+Drw9Z6v9Q2SiLislEZViMNk\n9kZdNu43vMnTQ6nSMGZafC9PhThMrA58L08NOEyaY9yrNMO2TqgbbwFZAQ6T6up3dGie40RH9Z4H\nHYFa9163eygz5jCpvnFuHIT8PYwm3jToOZQZcpjM1ji715fZlqrzrvcV5DCphvSP0q5du5cWsmW/\nPq48w6FRbakzD3lmwGFSDdneQHbh2qS9hOUEQx16Jnl4yFMyh4nVnYc8FeEwsaZzoJTEYWJt4EAp\ngcPE2sKBMmUOk/qZ5EpN2zlQpshhUl91nPCvAl/lmRKHiTWVr/KUzGFibeVAKZjDxNrMgVIgh0nz\necJ2OAdKQRwm7VHlubtZc6AUwGHSLqP2j20zB8qEHCbt4x7KYL5sPAGHibVR5S4bS5qX9Jik/ZKu\n7fP6hyQ9ImmvpN2Szp5FO4dxmJgdq/QeiqTjgO8C7yA5SfBB4JKI2Jcpcy7wnYh4QdI8sBAR63re\nZ2Y9FIeJtVnVeihrgcWIeDIiDgHbSQ5HXxIReyLihe7TB4DTSm7jQA4Ts8FmESirgKczz5/pfm2Q\nK4C7ptqinBwmZsPNYgvI3OMUSRcAlwPn9Xs9e1h674HORXOYWFt1Oh06nU6usrOYQ1lHMicy333+\nSeBwRFzXU+5s4KvAfEQs9nmf0uZQHCZmR1RtDuUhYI2k1ZJeCnwA2JEtIOl0kjD5cL8wKZPDxCy/\nWRyW/qKkq4FvAMcBN0fEPkkf675+E/AHwCuBG7urEg9FxNqy2+owMRuPF7YN4DAx669qQ57Kc5iY\nLY8DpYfDxGz5HCgZDhOzyThQuhwmZpNzoOAwMStK6wPFYWJWnFYHisPErFitDRSHiVnxWhkoDhOz\n6WhdoDhMzKanVYHiMDGbrtYEisPEbPpaESgOE7NyND5QHCZm5Wl0oDhMzMrV2EBxmJiVr5GB4jAx\nm43GBYrDxGx2GhUoDhOz2WpMoDhMzGavkoeld8t8vvv6I5LeNOz9HCZm1VB6oHQPS78BmAdeD1wi\n6ed7ylwEnBURa4CrgBsHvV/ZYZL3BDXXW686XW8xKnlYOvAe4FaAiHgAOFnSqb1vNIueSZP+z69q\nvW36WZtWb1UPS+9X5rTeN/Iwx6xaZhEoeU/n6j1I6Jjvc5iYVUslD0uX9CdAJyK2d58/BmyIiOcy\nZep55KFZAww6ObD0s43JHJYOPEtyWPolPWV2AFcD27sB9M/ZMIHBP5CZzU4lD0uPiLskXSRpEfgX\n4LKy22lm46vtYelmVj2VXylb9CK4ouqV9KFufXsl7ZZ09rTrzJQ7R9KLki6etM689Uqak/RtSX8r\nqVNGvZJOknSnpIe79V5aQJ3bJD0n6dEhZabxeRpa75Q+TyN/1m654j5PEVHZfyRDokVgNbACeBj4\n+Z4yFwF3dR+/Bbi/pHrPBU7qPp6ftN48dWbK/SXw58CvlPSzngz8HXBa9/nKkur9XeAzaZ3A94GX\nTFjv+cCbgEcHvF745ylnvYV+nvLUOY3PU9V7KIUtgiu63ojYExEvdJ8+QJ91MkXX2bUJuB34xwnr\nG6feDwJ3RMQzABHxvZLqPQyc2H18IvD9iHhxkkoj4l7g+SFFpvF5GlnvFD5PeX5WKPjzVPVAKWwR\n3BTqzboCuGvadUpaRfJLl96KUMQEWJ6fdQ1wiqT/JekhSR8pqd4bgNdLehZ4BLimgHqX066Jf7nH\nVMTnaaRpfJ5mcdl4HIUtgptSvUi6ALgcOK+EOq8HPhERIUkc+3NPq94VwJuBC4ETgD2S7o+I/VOu\ndx74m4i4QNKZwD2S3hgRP5ig3jyK/jzlr7i4z1MehX+eqh4oB4FXZ56/muQvxrAyp3W/Nu166U6c\nbQXmI2JU17KIOn+JZG0OJHMK75R0KCJ2TLnep4HvRcSPgB9J2gW8EZgkUPLUeynwGYCIeFzSE8Dr\nSNYyTcs0Pk+5FPx5yqP4z1MRE07T+kcSeI+TTNy9lNGTsusoZjIrT72nk0wqrivrZ+0pfwtwcUk/\n688Bf0EygXcC8Cjw+hLq/a/A5u7jU0kC55QCfubV5JuULeTzlLPeQj9Peeqcxuep0j2UmNEiuDz1\nAn8AvBK4sZvwhyJi7ZTrLFzO/40fk3Q3sJdkonRrRHxn2vUCnwa+KGkvSXf8dyLinyapV9JtwAZg\npaSngc0kQ7qpfZ7y1EvBn6ecdRbOC9vMrDBVv8pjZjXiQDGzwjhQzKwwDhQzK4wDxcwK40Axs8I4\nUMysMA4UMyuMA8WmTtLLZt0GK0ell95bPUj6WeBzwNeAjcA/AH9NcurjZuB+4Mc53kfd7zlMsunP\nCcDxEfHgdFpuRXOgWBF+EXhPRByS9D7gsxHxXUk/Ak6M/BsyXQN0gP8LvD8irpP0W5IejmQTJqs4\nD3msCPszv/A/GxHf7T5eDfxZnjeQtAJ4V0Q8DLwGOKn70j3A+wtsq02RA8UmFhHfBpC0hmRLgtQp\nkeyhgqR/J+n3lJyzhKTtPW/zduAHkj4K/CbdXdMiYi/JfqtWAw4UK9JakrmT1PGZxy8HDpFMlawh\n2Rog61ySLQxu7X7fnsxrPtStJhwoVqRzSCZgUyvSBxFxP/DmiNhDsnHR7p7v/RngQPeK0M90hz6p\nE6bUXiuYA8WKdA6QvSLzbz2v/7D733PpBoqSI2khOSLjx8DFwB/1fN/hIhtp0+NAsYlJeqOk3wbO\nBt4n6VXdl37YU/QpSRuBC7tXgVaRTLoC3EYSJidHxJd7vq/3fayifNnYJhYRj5Acc/GHPS89I+mV\nEfG8pN8guSR8kGS9ChFxUNJV3cePkuxVexRJZ5FsPWk14B6KTdNWkoVuAE+QTMy+G/j9TJmfGvEe\nFwG9PRarKO8pa1Ml6XzgHyLiqQGvHx8R/zrgtTOA0yOiM8UmWoEcKFZZkl4WESOX7Ft1OFDMrDCe\nQzGzwjhQzKwwDhQzK4wDxcwK40Axs8I4UMysMP8fw/CLsXreYxIAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106127090>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Above plot corresponds to Figure 15.2 in the book.  It can be seen that realized test variable is much smaller than that from replicated data; in other words, we might be underestimating the magnitude of prediction error.  Authors argue that this might be due to correlation of dependent variable within observations in the same year.  This motivates us to move on to hierarchical model.\n",
      "\n",
      "Note that the scale of this plot is not the same to that of Figure 15.2.  Please let me know if you find errors in my code!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    }
   ],
   "metadata": {}
  }
 ]
}