sandwich....ipynb

{
 "metadata": {
  "name": "sandwich..."
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import pandas\nfrom IPython.display import HTML\nimport statsmodels\nimport statsmodels.api as sm",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 282
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "OK. Laten we artificiele data genereren, met \n\n\n$\\beta = [1, -1]$\n\n$ X = \n\\begin{bmatrix}\n1 & 0.5 \\\\\\\n0.5 & 1 \\\\\\\n-0.5 & -0.5 \\\\\\\n1 & 0.5 \\\\\\\n0.5 & 1 \\\\\\\n-0.5 & -0.5 \\\\\\\n1 & 0.5 \\\\\\\n0.5 & 1 \\\\\\\n-0.5 & -0.5 \\\\\\\n\\end{bmatrix}$\n \n$Y = \\beta X + \\epsilon$\n\n$\\epsilon_t \\sim 0.5 \\times N(0, 1) + 0.5 \\times \\epsilon_{t-1}$"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "X = np.array([[1, 0.5],\n              [0.5 ,1],\n              [-0.5, 0.5]])\n\nX = np.tile(X, (3, 1))\n\nbeta_real = np.array([[1, -1]])\nY_real = np.dot(X, beta_real.T)\n\nar = 0.5\nnoise = np.random.randn(X.shape[0] + 1)\nnoise = [ar * noise[i-1] + (1-ar) * noise[i] for i in np.arange(1, len(noise))]\nnoise = np.array(noise)[:, np.newaxis]\n\nY = Y_real + noise\n\nplt.figure(figsize=(5, 5))\nplt.subplot(121)\nplt.title('Y', fontsize=24)\nplt.imshow(Y, interpolation='nearest')\nplt.axis('off')\nplt.subplot(122)\nplt.title('X', fontsize=24)\nplt.imshow(X, interpolation='nearest')\n_ = plt.axis('off')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAOUAAAFICAYAAACvPDmSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACNRJREFUeJzt3U+IlfUex/Hv78ykNaVcI0qSu8iFmxb9oyKoSKLaFQR3\nGy2ShiBQkuRCQ2NBF/rDLZIWktJGhGzVokUQl6BFBUIYRCWERpYQ1ZVw4Yj97uJyvURWM85wno+H\n12s38zzM94uc9zx6OPJrvfdeQIzR0AsAvyZKCCNKCCNKCCNKCCNKCDOWKO+5554ajUZ1yy231C+/\n/PKH9x48eLCmp6drNBrV/v37x7EeF7hDhw7VqlWrajQa1d69e//w3u3bt9doNKr169fXTz/9NKYN\nl6iPwZEjR/qaNWt6a60///zzv3vfwsJCv+6663prrT/wwAPjWI0JMTc311trfd26df2777475z0f\nffRRn5qa6qPRqL/11ltj3nDxxhJl773v2rWrt9b6JZdc0r/88stz3vPss8/+6R8snMvCwkK/9tpr\ne2utP/jgg7+5furUqT+8nmRsUfbe+5133tlba/2OO+74zbXPPvusr169uo9Go753795xrsWE+PDD\nD/vU1FRvrfUDBw786tr/nqSXX355P378+EAbLs5Yozx8+HCfmZnprbX+6quvnv3+mTNn+m233dZb\na/2+++4b50pMmK1bt/bWWr/qqqv6jz/+2Hvv/ZNPPukXXXRRH41G/Y033hh4wz831ih77/3FF1/s\nrbV+2WWX9SNHjvTee3/55Zd7a62vXbu2f/311+NeiQly8uTJvnHjxt5a6w899FA/ffp0v/HGGy+o\nX/hjj/LMmTP91ltv7a21fu+99/avvvqqX3rppb211l977bVxr8MEeu+993prrbfW+v33399ba33N\nmjX96NGjQ6+2KGOPsvf///uxtdY3bNjQW2v9rrvuGmIVJtQjjzxyNszWWt+1a9fQKy1a632Y/7r1\n3HPP1VNPPVVVVTMzM3Xo0KHauHHjEKswgQ4ePFg333xzVVVt2rSpPv/884E3WrzBPtGzY8eOuuKK\nK6qqanZ2VpCsmN57bd++/ezXhw8frg8++GDAjZZmsCinpqZqZmamqqrWrl071BpMoN27d9f7779f\nF198cW3evLl677Vly5ZaWFgYerVF8dlXJso333xTTz75ZFVVzc3N1f79+2vdunX1xRdf1DPPPDPw\ndosjSibK7Oxs/fzzz3X99dfXjh076sorr6yXXnqpqqpeeOGF+vTTTwfe8M+Jkomxb9++euedd2p6\nerr27NlTo9F/X94PP/xw3X333XX69OnasmVLDfTe5qKJkonw/fff19atW6uq6oknnqgbbrjhV9d3\n795dMzMz9fHHH9crr7wyxIqLJkomwuOPP14//PBDbdq0qXbu3Pmb69dcc83Z78/NzdXRo0fHveKi\niZIL3ttvv11vvvlmjUajev3112vVqlXnvG/btm1100031cmTJ2t2dnbMWy7eoFG21qq1NuQKXOBO\nnDhRjz32WLXW6tFHH63bb7/9d+8djUa1Z8+emp6ernfffbf27ds3xk0Xb7BP9ADn5q+vEEaUEEaU\nEGZ6uT/gH7VtJfY4L3+vfw42m8WbH/DNvPljg42uuvr83q7xpIQwooQwooQwooQwooQwooQwooQw\nooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQwooQw\nooQwooQwooQwooQwyz5evW1eqVWWrv9ruNkswbcDHoW3YbDRNX+eaXlSQhhRQhhRQhhRQhhRQhhR\nQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhRQhhR\nQhhRQhhRQhhRQhhRQhhRQpjlH4XX5ldolaXrfbjZLN6gr5FjOwebXVc7Cg8mgighjCghjCghjCgh\njCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCghjCgh\njCghjCghjCghjCghjCghjCghzLKPwqsDbYVWOQ9/W97qjMd8G+41Mn9ssNGOwoNJIUoII0oII0oI\nI0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oII0oI\nI0oII0oII0oII0oII0oII0oII0oII0oIs+zzKf9y6vhK7bJk/169frDZLMG3A55PuWGw0TV/nml5\nUkIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIY\nUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUIYUUKY6eX+gBMXD3fMWS3rED/GpW14erDZ/djO\nwWafL09KCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNK\nCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCLPso/D6X9evxB7nO33A2SzW03Xh\nHUc3JE9KCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNK\nCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCCNKCNN678s6T26+tZXaZemzl7c64/Lt\ngK+RDYONPu/XpyclhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBEl\nhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhBElhFn2UXht20qtsnT9n8PN\nZvFamx9sdj+2c7DZdbWj8GAiiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLC\niBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCiBLCLPsoPGBleVJC\nGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCGFFCmP8Ah5AQ\nZqLlzRsAAAAASUVORK5CYII=\n",
       "text": "<matplotlib.figure.Figure at 0xd5ac8d0>"
      }
     ],
     "prompt_number": 283
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Now let's calculate OLS using\n\n$\\beta = (X^T X)^{-1} X^T Y$\n\nThe explained timecourse\n\n$\\hat{Y} = X \\beta$\n\nThe residual\n\n$\\hat{r} = Y - \\hat{Y}$\n\nThe OLS variance\n\n$\\hat{V} = r^T r (X^T X)^{-1}$"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#OLS\ncalc_beta = np.linalg.pinv(X.T.dot(X)).dot(X.T)\nbeta = calc_beta.dot(Y)\nprint 'beta:\\n%s' % beta # more or less right...\nY_ = X.dot(beta)\nresid = Y - Y_\n\n\n\nols_variance = np.linalg.pinv(X.T.dot(X)) * resid.T.dot(resid) / (X.shape[0] - X.shape[1])\nprint '\\nols_variance:\\n%s' % ols_variance\n\nplt.figure(figsize=(12,5))\nax = plt.subplot(111)\nplt.plot(Y, label='$Y$', lw=2, c='g')\nplt.plot(Y_, label='$\\hat{Y}$', lw=2, ls='--', c='b')\nplt.plot(resid, label='$\\hat{r}$', lw=2, ls='--', c='r')\n\n#plt.fill_between(np.arange(9), Y.ravel(),  Y_.ravel(), color='purple', alpha=0.2) #np.zeros(9), \n\nplt.grid()\nax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.05),\n          ncol=3, fancybox=True, shadow=True, prop={'size':20})",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "beta:\n[[ 0.77988065]\n [-0.83997098]]\n\nols_variance:\n[[ 0.09358095 -0.04679048]\n [-0.04679048  0.09358095]]\n"
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 284,
       "text": "<matplotlib.legend.Legend at 0x11ac7b50>"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAE1CAYAAAAVjX0SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VMXXwPHvpgGhd5AWOqGDFOldijTpCAI/VBAF9RVF\nQBQQFAQVpRdFROlFmoAkklBDL9ICCAQIgQDSIaRs5v1jJLQEstm7uZvN+TxPHt3s3b0nhy1nZ8/M\nWJRSCiGEEEIIIUQcN7MDEEIIIYQQwtlIkSyEEEIIIcQTpEgWQgghhBDiCVIkCyGEEEII8QQpkoUQ\nQgghhHiCh9kBCOFqgoOD+eabb1izZg3h4eFmhyOEEEmSO3duWrZsyUcffUSpUqXMDkeIZGfXSHLv\n3r3JnTs35cqVi/f6wMBAMmfOTKVKlahUqRKjR4+253RCOL3g4GAaNGiAj48PQUFBxMTEoJSSH/mR\nH/lJUT8xMTEEBQXh4+NDgwYNCA4ONvvlVYhkZ1FKJXmd5C1btpAhQwZ69OjBoUOHnro+MDCQ7777\njlWrVtkVpBApxZtvvomPjw/Dhg0zOxQhhDDEqFGjOHToEHPmzMHb29vscIRINnaNJNepU4esWbM+\n8xg7anAhUpw1a9bQrVs3s8MQQgjDdO/eHX9/f+bMmUNYWJjZ4QiRbOwaSQYICQmhVatW8Y4kb9q0\niXbt2pE/f37y5cvHN998Q+nSpZ8OwmKxJwQhnEpMTAzu7u5mhyGEEIawWq14enri6+uLm5sbR44c\nkQEwkaIl9vHr0NUtKleuzPnz5zl48CADBgygbdu2CR5rdv+Vq/wMHz7c9Bhc6cfWfAJSIAshXIq7\nuztKKYYNG8abb75JSEiIvBc5+XuR/CT8YwuHFskZM2aM619q3rw50dHRXLt2zZGnTPVCQkLMDsGl\nSD6FEOIhpRQ3btx47nHy2mksyac5HFokh4eHx1Xtu3btQilFtmzZHHlKIYQQQjiIxWIhNjbW7DCE\nSBZ2rZPctWtXNm3axNWrVylQoAAjR44kOjoagL59+7J06VKmTZuGh4cH3t7eLFy40JCgRcJ69epl\ndgguRfIphBC2k9dOY0k+zWH3xD1DgrBYbO4TEcIZyWNZCOGKLBYL8+fP5/Lly9SrV4+KFSuaHZIQ\nSWLL+7RsS+1iAgMDzQ7BpUg+hXho8uTJTJ061ewwRAogr53GknyaQ4pkIYQQzzV16lT279/Pnj17\nmDFjhtnhCCGEw0mR7GLq169vdgguRfIpBPz2229cvXqVn376idmzZxMWFiZzTMQzyWunsSSf5rBr\n4p4QQgjX1759e9KlSxd3eeTIkURERJgYkRBCOJ6MJLsY6VsyluRTCB4rkJ/1OyEekNdOY0k+zSFF\nshBCCJFKyORLIRJPimQXI31LxpJ8itRs165d1K5dm/Tp0+Pm5kbmzJlp3LjxUzunhoWFUaJECdzc\n3EibNi116tTh0qVLJkUtEpKcky/ltdNYkk9zSJEshDCcFFeuoVq1amzdupVZs2YBujfZ39//qZ1T\nX3jhBUaNGsXLL7/M+fPn2bJlC3ny5DEjZJEAmXwphO1kMxEXExgYKJ84DWRrPuWx/Lj58+fTvXt3\nevXqxezZs+M9ZtGiRfz888/8+uuv5MyZM5kjFIlx7949cuXKhaenJ+Hh4Xh5eT12/ZkzZxg8eDDz\n5s3Dw0PmgzujiIiIp/rI4/tdQmzdTETei4wl+TSObCYihHAKbdu2xdvbm99//52oqKinrj9z5gzL\nly9nzZo1UiA7MW9vb9q2bcvNmzdZvXr1Y9fdunWLgQMHMmPGDCmQnZhMvhTCdlIkuxj5pGksyad9\npLhyHd26dQNg3rx5cb+zWq3069ePcePGkSVLFrNCE05IXjuNJfk0hxTJQgiHkuLKNbz88stkz56d\ntWvXcvPmTQAGDRpEnz59KFasmMnRiedZtmwZAwYMoG3btkRERPDLL78wePBgunfvzu7du80OTwin\nJMM3Lkb6loyVXPm0jLQ4/BzxUcMd3z/9ZHGVOXNmly6uLAn8UybUAmfv8cnVAu/u7k7nzp2ZOnUq\nS5YsITIykjJlylCvXr3kCSC5udA/ZGRkJLt27WLSpEkUL16cjh07MmTIEFq1akWlSpXIli0bVatW\nNfSc8l5kLMmnOWQkWQjhUA+Kq6ioKJYsWcKUKVNcu7hyYQ++Ffjqq68ICwujd+/eJkckEmPTpk3U\nrFmTqKgoQkNDqVy5MrVq1eLu3btky5aNDh06mB2iEE5JVrcQwkDyWI5fUFAQtWrVwsfHh65du/Ll\nl1+aHZJIgtjYWLy9vSlfvjy7du0yOxyRSDt37sTX15cDBw5Qv359/v77b8qWLWvTfdi6uoUQzsqW\n92lptxBCOFz16tXx8vIiR44cUiCnYEePHiUqKopmzZqZHYqwQfXq1QHYuHEjOXLksLlAFiK1knYL\nFyP7uxtL8mkMKa5cw+bNmwGkVSaFCggISLZ/O3ntNJbk0xxSJAshHE6KK9ewefNmPD09qVmzptmh\nCBtFRESwc+dOmfwlhA2kSHYx8gJoLMmnMaS4SvmsViubN2+mYsWKsglFCrR9+3aioqKS7TVNXjuN\nJfk0hxTJQgiHkuIqZTt16hRNmjShTJkyhIeHExwcTMOGDZk7d67ZoQkbXLx4kWrVqlGmTBmzQxEi\nxZAi2cVI35KxJJ9JJ8WVayhatCh+fn4EBwdjtVq5efMmGzdupEePHmaHJmzQvXt3duzYkWznk9dO\nY0k+zSGrWwghHOJBcSWEEEKkRLJOshAGkseyEMIVyTrJwlXY8j4t7RZCCCGEEEI8QYpkFyN9S8aS\nfAohhO3ktdNYkk9zSJEshBBCCCHEE6QnWQgDWSwWYmJicHd3NzsUIYQwhNVqxdPTk3nz5klPskjx\npCdZCJPkzp2bc+fOmR2GEEIY5ty5c2TNmtXsMFKf8+fh5s3Hf2e1QliYOfGkQlIkuxjpWzKWrfls\n2bIl8+bNc0wwQghhgt9++41KlSrZdBt5L7LTzZvQrBnUrAnnz+t8Xr8OzZtDgwZPF8/CIaRIFsJA\nH330EVOmTGHUqFGcOXMGq9VqdkhCCGEzq9XKmTNnGDVqFD/88APNmzcHQCmFp6enydG5uOho6NgR\njh7VlzNm1P/18oLwcDhxAnr2hNhY82JMJaQnWQiDBQcH8+WXX/LHH39w48YNeWwLIVIci8VC1qxZ\nqVSpEs2bN+eFF14A4PLly3Tt2pVcuXKZHKGLUgr69oVZsyBnTti5EwoXfnj9qVNQpQrcuAGjRsGw\nYebFmkLZUnNKkSyEAyilWLp0KRcvXiRHjhxYLBazQxJCiCRTSnH16lVy5cpFp06dcHOTL6IdYvx4\nGDQI0qaFgAB46aWnj1m3Dl55Rf//H3/oFgyRaLbUnLIttYsJDAykfv36ZofhMpKaT4vFQsuWLdmw\nYQMhISFSJP/n4sWL5M2b1+wwXIbk0ziSy2eLjY2lUKFCNG3aNFEFsrwXJYFScPas/v9ff32sQH4s\nn82bwxdfwGefwY8/SpHsQFIkC+Eg6dKlo02bNty5c4fr169LfzKwe/duqlatanYYLkPyaRzJZcLc\n3d3JmjUrGTJkMDsU12axwKRJut/4eY/FoUMhXz54/fXkiS2VknYLIYQQQgiRKsg6yUIIIYQQQthB\nimQXI2tTGkvyaSzJp7Ekn8aRXBpL8pkI0dHw77+JOlTyaQ67iuTevXuTO3duypUrl+Ax7733HsWL\nF6dChQrs37/fntMJIYQQQqR8SkH//lCtGhw/bsx9Xr4MnTpBaKgx9yfs60nesmULGTJkoEePHhw6\ndOip69euXcvkyZNZu3YtO3fu5P3332fHjh1PByE9yUIIIYRILb75Bj7++NlLvdmqRw+9Kka1arB5\nM6RJY/99uqBk60muU6fOM/dzX7VqFT179gSgevXq3Lhxg/DwcHtOKYQQQgiRci1frtdCBpg715gC\nGeC776BQIdi1CwYMMOY+UzmHLgF34cIFChQoEHc5f/78hIaGkjt37qeO7dWrFz4+PgBkyZKFihUr\nxq0J+KAXRy4///KjfUvOEE9Kvyz5lHw682XJp3GXH/zOWeJJ6Zcf/M5Z4nGay9Onw/vvU18pGDOG\nwJw5ITDQmHzmyEHg0KEwYAD1Z82CqlUJLF7cuf5+Ey4fOHCAGzduABASEoJNlJ3OnDmjypYtG+91\nLVu2VFu3bo273KhRI7V3796njjMgDPGfgIAAs0NwKZJPY0k+jSX5NI7k0liSzwR89ZVSoNQbbygV\nG5vom9mUzzlz9Dm8vJQ6fNj2GF2cLTWn3eskh4SE0KpVq3h7kt9++23q169Ply5dAChVqhSbNm16\naiRZepKFEEIIkSqsWgXNmoGXl+PO0b8/ZMkCI0eCu7vjzpMCOc06ya1bt2bu3LkA7NixgyxZssTb\naiGEEEIIkSq0bu3YAhn0zn2jR0uBbCe7iuSuXbtSs2ZNjh8/ToECBZg9ezYzZsxgxowZALRo0YIi\nRYpQrFgx+vbty9SpUw0JWiTs0f4lYT/Jp7Ekn8aSfBpHcmksyaexbM6nxeKQOFIbuybuLViw4LnH\nTJ482Z5TCCGEEEKkTBERkC6d2VGIJLK7J9mQIKQnWQghhBCuZMUKeO893YNcsaLZ0UBYGNy8Cb6+\nZkdiKqfpSRZCCCGESHX27IHXXoPz52HDBrOjgWPH4MUXoWVLuH7d7GhSDCmSXYz0gRlL8mksyaex\nJJ/GkVwaK1Xn89w5aNVKt1r87396Zz072Z1PHx/ImxdOn4Zu3SA21u6YUgMpkoUQQgghjHDrFrzy\nCly6BA0bwvTpzjGJLl06vdNf9uywbp1eGk48l/QkCyGEEEIYYdEi6NIFSpWC7dsha1azI3qcn59e\nozk2Flau1MvRpTLSkyyEEEIIkdw6d4bFi+GPP5yvQAZo0gS++kqPbp88aXY0Tk+KZBeTqvvAHEDy\naSzJp7Ekn8aRXBorVeezY0coUsTQuzQ0n4MGwc6dMHCgcffpoqRIFkIIIYRILSwWqFrV7ChSBOlJ\nFkIIIYRIithYcJPxxpREepKFEEIIIRzp3DmoVAmCgsyOxBiyLNxTpEh2Mam6D8wBJJ/GknwaS/Jp\nHMmlsVw+n7du6Y05/v4bRo1y+Okcns+AAChbFs6edex5UhgpkoUQQgghEismRq9icegQlCwJ8+aZ\nHZH9vvtO78rXrp3eBEUA0pMshBBCCJE4SsG778K0aZAjh14lwuCVLExx7ZqezHf6NPTqBbNnO8cm\nKA5gS80pRbIQQgghRGIcOAAvvgienrBxI9SsaXZExjl4EGrU0CPJU6dCv35mR+QQMnEvFXP5PrBk\nJvk0luTTWJJP40gujeWy+axYUe9UN3dushbIyZLPChXgxx/1/3/4od5aO5XzMDsAIYQQQogUo2VL\nsyNwnNdegxMnoHZtyJPH7GhMJ+0WQgghhBAiVZB2CyGEEEIIIewgRbKLcdk+MJNIPo0l+TSW5NM4\nkktjuUQ+Y2J0+4G/v9mRuEY+UyApkoUQQgghHqUUvP8+LFgAPXrAvXtmR2SuP/+EffvMjiLZSU+y\nEEIIIcSjvv8e/u//wMtLL/VWq5ZpoVhjrZy9eRafLD64WUwY21y3Dl55BQoWhL17IXv25I/BQNKT\nLIQQQgiRFCtX6iXQAObMMaVAPn39NDP3zqTjko7k+iYXRScW5b117yV7HAA0aABVqugtq7t2BavV\nnDhMIEWyi5G+JWNJPo0l+TSW5NM4kktjpdh83rypd5xTCkaN0kVhMvj33r8sPbqUvmv6UnRiUYpO\nLErfNX1ZenQp1yKuQQhM3T2VXRd2JUs8j0mbFpYtg5w5wc8Phg1L/hhMIuskCyGEEEIAZM6sC8IV\nK+DTTx12mvsx99l2bhv+Z/zxO+XHvov7UDxsAciaNisNCzekSZEmNC7SmGE/D2PhnYW8u/Zddryx\nA3c3d4fFFq8CBWDRImjSBMaO1SPL7dsnbwwmkJ5kIYQQQggHilWxHLx0EP/T/vid9mPLuS3cj7kf\nd72Xuxe1C9amceHGNC7SmMp5Kz9WCN+JukOpyaW4cPsCM1rOoM+Lfcz4M+C77+Drr2HJEqhb15wY\n7GRLzSlFshBCCCGEwc7eOBtXFP915i+u3rv62PUV81SkcZHGNCnShNoFa+Pt6f3M+1t8ZDGdl3Ym\nW7psHO9/nBzeORwZfvyUgn//hRwmnNsgUiSnYoGBgdSvX9/sMFyG5NNYkk9jST6NI7k0VmrM5/WI\n6wSEBOB/2h//0/6cvHbysesLZCpAk6JNaFKkCQ0LNyRX+lyJvu/AwEDq1atHk1+b8NeZv3ir8lvM\nbDXT6D8hVbCl5pSeZCGEEEKkTlOnQv780Lq1zTeNjIkkKDQobrR4T9geYlVs3PWZ02SmYeGGNC6i\nWyiKZyuOxWJJcqgWi4XJLSZTflp5ftz3I29WfpNq+aol+f7E88lIshBCCCFSn1WroG1bsFggOBiK\nF3/m4UopDl0+FFcUbz67mXvRDzcZ8XTzpGaBmnEtFC++8CIebsaPRQ72H8zX276mygtVzJnEF59r\n1yBbNrOjSBRptxBCCCGESMi+fVCnjt5J74sv4LPP4j0s9FYofqf88D+jWygu37382PXlcpWLK4rr\nFKpDBq8MDg/9TtQdfKf4EnorlOmvTKdvlb4OP2eCYmN1/iZNgt27oUgR82JJJCmSU7HU2AfmSJJP\nY0k+jSX5NI7k0lhOnc/z56F6dbh4UW85PWeOHk0Gbt6/SWBIYNzSbMf/Pf7YTfNlzBdXFDcq0og8\nGfIkS8hP5nPJkSV0WtqJrGmzcmLACXMm8YEuktu2hdWroUIF2L4dvJ89AdFs0pMshBBCCPEkpaBT\nJ10g16tH9LQp7Di3Na4o3nVhF1b1cEe5jF4ZaVC4AY0LN6ZJ0SaUzF7Srr5io3Qo3YHGRRrjf9qf\nIX8NYVarWeYE4uYGv/4KVavCwYPw1lvw229xHzpSOhlJFkIIIUSqoJTizPoFpBn6GQP7FeWPf4O4\nE3Un7noPNw9eyv9SXFFc9YWqeLp7mhhxwoKvBlN+WnmiY6PZ8cYOquevbl4wR47o0fm7d+H77+H9\n982L5Tmk3UIIIYQQAgi7HRa3LJv/aX8u3rkICvhvsLN0ztJxLRT1CtUjY5qMpsZriweT+F7M+yI7\n39xp7iS+JUv0KH3BgnoiZLp05sXyDFIkp2JO3QeWAkk+jSX5NJbk0ziSS2OZmc/bkbfZdHZT3CoU\nR68cfez6PBnyxG333KhwI/JlymdKnLZIKJ+PTuKb9so03q7ydvIH96hZs6BVK8iTPL3aSSE9yUII\nIYRIFaKt0ewO2x1XFO8I3UFMbEzc9ek901Pfp37caHHpnKWdoq/YCBm8MjCh6QQ6LunI0L+G0t63\nPTnT5zQvoLfeMu/cDiAjyUIIIYRIMZRSHP/3eNzSbAFnArgddTvueneLO9XyVaNxkca0jShE+ete\neHR73cSIHUspRdPfmuJ32o83Kr3Bj61/NDskpybtFkIIIYRwGeF3wnVP8X/rFYfeCn3s+pLZS9Kk\naBMaF25MfZ/6ZE6bGUJD9WSysDBYuTJJu+qlFMevHqfctHJEx0YT9EYQL+V/yeyQnJYtNaebvSdb\nv349pUqVonjx4nz99ddPXR8YGEjmzJmpVKkSlSpVYvTo0faeUjxDYGCg2SG4FMmnsSSfxpJ8Gkdy\naSx783k36i7rTq5j4IaBlJ9Wnjzf5qH7792Zc2AOobdCyZU+F6+Ve43ZrWdz7oNzBPcPZlLzSbQp\n1UYXyLdvQ8uWukCuWxeaNjXmDzPJ8/JZMkdJBtYcCMC7a9/FGmt95vHJxmqFTz/VG42kQHb1JFut\nVvr374+/vz/58uWjatWqtG7dGl9f38eOq1evHqtWrbIrUCGEEEK4ppjYGPaG7cXvtB/+p/3Zfn47\n0bHRcden80hHPZ96cUuzlc1VFjdLAuN8MTHQpYtet7dECfj9d0iTJpn+EvMMqzOMeX/PY9/Ffczc\nO5N+VfuZHRJMnQpffaXXUt6zB3LlMjsim9jVbhEUFMTIkSNZv349AGPHjgVg8ODBcccEBgby7bff\nsnr16oSDkHYLIYQQItVQSvHPtX/iiuKNZzZyM/Jm3PVuFjeqvFAlbrJdjfw1SOORyEL344/hm28g\ne3bYsQOKFXPQX+F8lh1dRoclHciSNgsn+p8wdxIfQFQU1K8PQUHQoAFs2AAe5q4ZkWyrW1y4cIEC\nBQrEXc6fPz87d+58Kpjt27dToUIF8uXLxzfffEPp0qWfuq9evXrh4+MDQJYsWahYsWLccicPvmaQ\ny3JZLstluSyX5XLKvHzl7hUmLZ7E3rC9HE5/mHM3z0EImg8Uy1aM0ndL82LeFxnQaQBZ02XVtw+B\nND5pEn++UqWoX7QozJlDYGgohIY6xd+fHJezhWejSlQV9rCHIX8NoXum7ubGt307fPgh9fv3h4AA\nArt1g379kjWeAwcOcOPGDQBCQkKwhV0jycuWLWP9+vXMmqW3Q/ztt9/YuXMnkyZNijvm9u3buLu7\n4+3tzbp163j//fc5ceLE40HISLJhAgMD4x4cwn6ST2NJPo0l+TSO5NJYgYGBVKtVja3ntsaNFh+4\ndOCxY7Kny07jIo3jfnyy+BgXQHQ0eDrnTnlJYcvj88S/Jyg7taxzTeLbulWPJMfEwNKl0L69aaEk\n20hyvnz5OH/+fNzl8+fPkz9//seOyZjx4c41zZs355133uHatWtky5bNnlMLIYQQwglFW6MZs3UM\nm7ZuItIaGff7tB5pqVOwTlwLRYU8FRLuK7aXCxXItiqRvQQf1fyIMVvH8M4f77D7rd3m7sQHULs2\nTJgAK1boiZQphF0jyTExMZQsWZK//vqLF154gWrVqrFgwYLHJu6Fh4eTK1cuLBYLu3btolOnTk8N\nd8tIshBCCOEaJu6cyPvr38eChcp5K8ctzVarYC3SeqQ1O7xU4W7UXXyn+HL+1nmmtJjCO1XfMTsk\nUApiY8Hd3II92UaSPTw8mDx5Mk2bNsVqtfLGG2/g6+vLjBkzAOjbty9Lly5l2rRpeHh44O3tzcKF\nC+05pRBCCCGc1JW7VxgeOByA3zv/TptSbRx7wjt3YO1a6NTJsedJYdJ7pef7Zt/TfnF7Pt34KR1L\ndzR/Ep/FYnqBbCvZTMTFSF+dsSSfxpJ8GkvyaRzJpTH6runLzL0zqRpVlZ2jdzp2+2erFdq2hTVr\n9GoWAwc67lwmS8rjUylFs3nN2HBqA70r9ean1j85JrgUJlk3ExFCCCGE2HdxH7P2zsLDzYP+1fo7\ntkAG+PBDXSBny+bSu+kllcViYVLzSXi6eTJ7/2yCzgeZHdLToqPh0CGzo0iQjCQLIYQQwi5KKer8\nXIdt57fxYY0P+fblbx17wkmT4L33wMsL/P2hTh3Hni8F+3Tjp3y15Ssq5qnInrf2mD+J74Hbt6FV\nKzhwQO/IV7x4spxWRpKFEEIIkWwWHF7AtvPbyJU+F5/X/dyxJ1u7Fj74QP//Tz9JgfwcQ2sPpWDm\nghy4dIDpe6abHc5DGTLoDV9u3oRXX9X95U5GimQX82AhbWEMyaexJJ/GknwaR3KZdHei7jDIbxAA\nYxqNIXPazI7NZ9GiULgwDB8O3bs77jxOxJ58pvdKz/dNvwdgWMAwLt+9bFBUdrJYYM4cKFUKjhyB\nN97QK2A4ESmShRBCCJFkY7aO4cLtC1R5oQq9KvZy/AlLloS9e3WRLBKlbam2NC3alBv3bzDYf7DZ\n4TyUMaNeOzljRli8GL77zuyIHiM9yUIIIYRIklPXTlF6ammirFFs772dGgVqmB2SSMDJf09SdlpZ\noqxRbOu9jZoFapod0kMrVuiWi5YtYeVKcHPcGK70JAshnNKNG3rVJiGEaxi4YSBR1ih6VOghBbKT\nK569OB/X/BiAd9e+S0xsjMkRPaJtW/Dz08WyAwtkWzlNJJ06wb//mh1Fyid9dcaSfBpnxw7ImzeQ\nChXgkd3shR3k8WkcyaXtNpzawMrjK8nglYGxjcY+dp2h+QwIcLpe1eRmVD6H1nHSSXwAjRs73WYj\nTlMkL1miJ6ieO2d2JEIIo125Ah07wv37en5G375mRySEsEe0NZr3178PwGd1PyNvxrwAREbqnYcN\nM2kSNGz4cDULYRdvT++Hk/g2OtEkPiflND3JZcsqDh+GF16A9euhXDmzoxJCGMFqhebN9TdpL72k\nJ6WPHw/58pkdmRAiqSYETeDDDR9SLFsxDvc7TBqPNCgFPXvCtWswd67e4+OB69cha1YbT/LHH3qT\nkNhY+PXXVLOShaMppWgxvwXr/1lPr4q9+LnNz2aHlDCl9CoYBkqRPclbtkDduhAWpkeUN282OyIh\nhBGOHoXt2yFHDv2N0fz5UiALkZJdvnuZEZtGAPB90+9J45EGgJkzdS0bEAAXLz48fv9+KFIEZs2y\n4SQHDkDnzrpATkVLvSUHi8XCxGYT8XL3Ys6BOWw7t83skOIXFqa/RQgyb6dApymSs2SBP/+E9u31\nutJbt5odUcokfXXGknzar1w5vZnS8uXwzz+BZofjUuTxaRzJZeIN/WsotyJv0bxYc14p8QoAe/bo\nDfBAF8tXrgTGHb9xo56026cPfPFFItqLL1zQqxzcvQvduslSbxj/+CyevTiDaum1rZ1uEt8D06ZB\nYCB06ACXLpkSgtMUyQBp08KiRfpnyBCzoxFCGMXX99mbYsXGwrJlqX5ujhBOb0/YHmbvn42nmycT\nmk4AdHtFhw4QFQX9+um69lEDB+p6x81N17vvvPOcVW6sVj1yVru23lHP4K/bhTak9hAKZS7EwfCD\nzjeJD+Dzz/UbR1iYXt0hOjrZQ3CanmQnCEMIYZJPPoFx4/RI05Qp4OFhdkRCiCcppag1uxZBoUF8\nXPNjxjUZB+hBrbFjoWpV3TqZJk38t//9d+jaVU/u69IFFix4xslu3oSYGL1tsXCYFcEreHXRq2RO\nk5nj/Y+TO0Nus0N63KVL8OKLulB+7z344Qe77zJF9iQLIVKvmjX1N0kzZ+oRqYgIsyMSQjxp3qF5\nBIUGkTvGwwHAAAAgAElEQVR9bobVHRb3+y++gKFD9ZyDhApk0HtF+PvrCX1t2jznZJkzS4GcDNqU\nbEPzYs25GXmTT/w/MTucp+XJo79m9PSEiRPhr7+S9fQppkg+dQr+9z+4d8/sSJyb9NUZS/Jpu4UL\ndbEb3wf1hPLZpo1+88yaVW+21KSJ/gpXPJs8Po0juXy225G3GeSne1i/bvw1mdJkirvO0xO+/BIK\nFXp4fEL5rF1bv5936eLIaF2Pox6fFouFic31JL5fDv7inJP4XnoJJk/Wn8YaNEjWU6eIIlkpeO01\nmDNHrzUtm44I4ZyOHIE33tDrIG/caNtta9XSE3bz54dt22RZVCGcyZdbvuTinYtUy1eN1yu8btd9\nZcnyxC+OH5cJCSYqlq0Yn9TSo8hOO4mvTx/47LNk340vxfQkHzsGzZrpzUZKldIrYRQsmEwBCiGe\n6/Zt3ZN4/LherWnu3KTNtwkNhQED9Gh0zpzGxymEsM3Jf09SdlpZoqxR7HxzJ9XyVTPuzg8e1MPL\nrVtz64efyZTDy7j7Fol2L/oepaeU5uzNs0xsNpEB1QeYHZLDuGRPsq+vXmu1bFkIDoYaNeDQIbOj\nEkKAHgR66y1dIJcpA9OnJ31Cev78eoKPFMhCOIcPN3xIlDWKXhV7US1fNX78Ea5eNeCOw8L0Um93\n7nDlCvgU92TVKgPuV9jM29ObH5rpSXHDAoYRfifc5IicQ4opkkFvQPDopiNbtpgdkfORvjpjST4T\nZ+pUvXRjxox6jkX69PEfJ/k0luTTOJLL+K07uY41J9aQ0SsjYxqNYckS/YG4Rg29SkVCnpvPu3eh\nVSv91VGtWnxX5ieu37Dw6qs2bjqSSiTH47N1yda0KN6CW5G3GOQ/yOHns9v58w6fyJeiimR4uOnI\nvHl6rUUhhPkaN9bf8syeDSVLOuYcMTGwbp1j7lsI8bQoaxQf/KknBwyvN5ybF/LQu7e+bsCAZ69k\n8UxWq55otG8fFC0KK1bw1XdpGT5cr5me6E1HhKEe7MSXxj0Ncw/OZes5J97V7exZqFxZL5kSHOyw\n06SYnmQhhHOLjtaz3B2lb1/dp/z55zBihOwvIISjfbv9Wz7y+4iS2UsS1ONv6tT04sgRva/DwoV2\nPAevX4eXX9bLXOzYASVKxF01fTq8+64ulj/4ACZMMOZvEYn3ecDnjNo8ivK5y7O3z1483Jxw4Xql\n9KLbixbpkZlduyBTpuffDhftSRZCODdHFsgAVaroic1ffAFvv61HloUQjnHpziVGbhoJwISm3/Pe\nu7pALlUKfvzRzg+pWbPCpk0QEPBYgQz6ub10KWTIAI0a2XEOkWSDaw/GJ4sPf4f/zdTdU80OJ34W\ni96NsWxZPRmmZ0/9ycpgLlUkHz6sP3mm5jdP6aszluTTWDOWzuDCrQtJuu1bb8Hy5bLpyKPk8Wkc\nyeXjhvw1hNtRt2lZoiVNizYjb14912DpUj334Hmem09vb6hQId6rXn0VzpzRc/qElpyPz0cn8X0W\n8BmX7lxKtnPbJH16/aaQOTOsWKG3fTSYyxTJMTHQrp3esbBdO9l0RAhHOnfO9n7BdSfX8fYfb1Nq\nSimWHV2WpPM+uenIsGHPv40Qwja7LuxizoE5eLl7MaHpBNzc9LbxJ07o1WuSQ44cyXMeEb9WJVrx\nSvFX9CQ+PyeexFe8OPz2G3h4gJfxywe6VE9yUJD+5Hntmp55u3q17GophNEezJdo1kx/25U27fNv\nczfqLmWmluHszbNxv/uo5keMaTQmSf1uR4/CoEF6Am/mzDbfXAiRgFgVS42farDrwi4+qfUJYxsb\nMDp37Zp+orq7231Xd+8mvHqOMNapa6coM7UMkdZINvfaTJ1CdcwOKWEhIeDjk6hDU21Pco0aeseu\nggV1wVy7th7xEkIYIzJStzlcuwY3biT+g/vwwOGcvXmWinkq8t3L3+Hh5sE327+h0dxGSfoqr3Rp\nWLNGCmQhjPbrwV/ZdWEXeTPk5dM6n9p/h3fv6n3m27fX/2+HDRugSBG9I6dwvKLZivJJbSffie+B\nRBbItnKpIhkebjpSrpxeFWSrE69g4gjSV2csyefjPvwQ9uyBQoXg118Tt0Povov7mLBjAm4WN/pm\n78v/1fg/AnoGkDdDXjaf3UzlGZWde6khJyaPT+NILuFW5C0+8ddF0ei648mYJhHNxwkIDAx8fKm3\nw4ftnkQwdy5cvqyXnExtm46Y9fgcXEtP4jt0+RBTdk0xJQYzuVyRDHrTkc2bdZvKa6+ZHY0QrmH+\nfL1piJeXnryTLdvzbxMTG0Of1X2IVbEMqDaAUjlKAVC7YG329d1HvUL1uHjnIg1+acD3O763u+0q\nKkpPmhdC2G705tGE3w3H99YARnZ4jaAgO+/w4491NZs1K/zxh92NxnPm6DWU799HNh1JJuk80zGx\n2UQAPg/8nIu3L5ockQ2MmJymnICThCGESEBsrFI1ayoFSk2blvjbTQiaoBiBKvBdAXXr/q2nro+2\nRquPNnykGIFiBKrzks7qduTtJMVotSrVtatSbm62xSiEUCr4SrDy/MJT8WF+lSVblAKlhg+34w4n\nTtQvGJ6eSgUEGBSlfi0aPlzfNSg1bpxhdy2eoeX8looRqO7Lu5sdSuKsW6dUrlxKbd781FW21Jwu\nOZIshDCWxQJ+fnrkpm/fxN3m3M1zDNuol5+Y0mJKvF/derh5ML7JeJZ2XEpGr4wsOrKIarOqEXzV\n9h2ULBbdbhUbC/36wfDhsmOXEIn1f3/+H9HRilxrA7hxzZOXX4bPPrPjDh9MCJo1C+rXNyJEQD/P\nR4zQm46kS6fXTxeO90OzH0jjnobf/v6NzWc3mx3O8/31l+7N6dgRLiRt2VFw0XaLZ9m927W3u5S+\nOmNJPh/y9oY330zcJgJKKd5d+y53o+/S3rc9rUq2AhLOZ/vS7dn91m5K5yzNsavHqDqrKkuPLrUp\nPotFv6nPnJl6Nh2Rx6dxUnMu/zjxB+v+WYfXxu+5HFyM/Pl1u6I9i1EEVqqkK9mePY0L9BF9+8Lp\n09CggUPu3umY/fgskrUIg2sPBvQkvmhrtKnxPNeYMfrBER6uZ5tHRibpblJVkXzvnl5ndfhw/Wbv\nym+eQphp2bFlrDmxhkxpMjGx+cRE3aZkjpLsfHMnXcp24U7UHTou6cjADQNtfjF+6y34/feHm46M\nH5+Uv0CI1CEyJpIP/vwALpchatu7eHjA4sWQM+dzbhgRAb/8Am+8Ef/1L7yQ+K+dkihPHofevXjC\nJ7U+oXCWwhy+fJgpu518Ep+Hh96yukABvfX5Bx8k6W5cap3kxFi9Gjp31s/vli11Dr29k+XUQqQK\nN+7fwHeKL5fuXGJqi6n0q9rPptsrpZi0axIDNwwkJjaGOgXrsKjDIvJmzGvT/WzfDl9+qd/wZV1V\nIeI3bts4PvH/BN8cvnyZ/2+uhHvQp88zbvDPP3qEePZsuH5d/27XLqhaNVniTYz79xO3fruw3ZoT\na2i1oBUZvTJyvP9xm1+Xk92ePXo94JgYOHQIfH1tqjlTXZEMj2868tJLer1V2XREiIeuXIGRI/U3\nVonZgvZR/f7ox/Q906lZoCZb/rcFN0vSvrDadm4bnZZ2Iux2GHky5GFxh8XOvZi9ECnMxdsXKTG5\nBHei7vBn9z95uejLz77BW2/Bjz8+vFylCrzzjh55cpLRpuXL9VKVa9fq9dSF8VotaMWaE2voXr47\nv776q9nhPN/8+ZA/P9StC6TizUQS69FNR3bu1CPxrsLsviVXkxrz+WBp0ylToH9/22677dw2pu+Z\njoebBzNazniqQLYln7UK1mJfn33U96nPpTuXaPBLAyYETUjWD9TOLjU+Ph0lNeZy8F+DuRN1hzYl\n2zy/QAa9BXDatNC7t57gs3s3/O9/8RbIZuRTKb1M5dmzevDQlTYdcabHZ4qbxPfaa3EFsq1SZZEM\nDzcdmTcPXnnF7GiEcB4jR4K/v+5J/OqrxN8uyhpFnzX6e9pBtQZRNldZu2PJnSE3fq/78XHNj7Eq\nKx9u+JAuy7pwO/J2ku/z3j39DZwQqVnQ+SDmHpxLGvc0fNf0u4dXKJXwagBvvw1hYXo/eidcVsJi\n0S2VrVvrTpDGjWHlSrOjcj1FshZhSO0hQAqZxGcP+5eiW6dKliypihUrpsaOHRvvMQMGDFDFihVT\n5cuXV/v27XvqegPCEEIYYO1avfaom5tS/v623Xb0ptGKEahiE4upe1H3DI9t2dFlKuNXGRUjUKUm\nl1JHLx+1+T6io5Vq1UqpNGmU+v13w0MUIkWwxlqV76fdFb1rqaF/DdW/vH5dqe+/V6pUKaUKF1Yq\nJsbcIO0QHa1Unz4PX8vmzDE7ItdzL+qeKvJDEcUI1HfbvzM7HJvYUnPaNZJstVrp378/69ev5+jR\noyxYsIBjx449dszatWv5559/OHnyJDNnzqRfP9sm8Qghkse5c9C9u/7/L76ARo0Sf9uT/55k1OZR\nAEx/ZTrpPNMZHl8733bs6bOHMjnLEHw1mGo/VmPJkSU23YfFoifdR0ZC+/Z6/pEQqc3kwAUcmzIK\n5gTS7EgLvdzTCy/oFQCCg/UTJCTE7DCTzMNDP7eHD9fdIaVKmR2R63l0J77hgcNT1k58NrCrSN61\naxfFihXDx8cHT09PunTpwsonvttYtWoVPf9bJ7F69ercuHGD8PBwe07rcJs2wbRpZkeRNM7Ut+QK\nUlM+c+aEdu2gRQsYMiTxt1NK0XdNXyKtkfSs0JNGRRKuru3NZ4nsJdj55k66lu3Knag7dFraiQ//\n/DDRX/e5u+vn9hdfuMamI6np8eloqSWX1+/d5ON3csMNH4qWuUGtXz7Q7RMREfqT8dKlukAuWtSu\n85idzwebjhw/DtWrmxqKIczOZ3xeKfEKrUq04nbUbT72+9jscBzCw54bX7hwgQIFCsRdzp8/Pzt3\n7nzuMaGhoeTOnfux43r16oWPjw8AWbJkoWLFitT/b5eeBw+O5Lh89Sq0bBnInTtw4UJ9Ro2CTZuS\n7/xyWS6beXnWLPDzC2Tz5sTffuhPQwnYFkB23+x88/I3Do939/bdvJXtLWo2r8n//fl/TFg4Af+N\n/vw57E/yZsz73Ntv2hRInTowc2Z93n4bvvgikIgIGDfOMfHK5ZRx+QFnicdRl+u0m05UcHU80t/E\nf1V2Nv/cAgr7UH/0aChZUh+/bZvL5POffwL55x/nyX9Kz+eTl7tm6Mqf5/9kHvN4q/JbqBDlVPEF\nBgZy4MABbty4AUCIjd+Q2LUE3LJly1i/fj2zZs0C4LfffmPnzp1MmjQp7phWrVoxePBgatWqBUDj\nxo0ZN24clStXfhhEMi8B9zyzZ0OfPnqWf+/eMGOG/vpGCPG4K3evUGpKKa5FXOOXtr/Qo0KPZD1/\n0PkgOizpELdM3KIOi6hbKPGzmFet0rvmLl0KadI4MFAhzBITo2ezXbvG7MyNeaOTHrSa/Ntp3n2t\nmMnBmScqCry8zI7CNXyx6QuGBw6nTM4y7O+7H093T7NDeqZkWwIuX758nD9/Pu7y+fPnyZ8//zOP\nCQ0NJV++fPac1uF694YVK/S+8LNnw6uv6hnxQojHDdwwkGsR12hUuBGvl3892c9fo0AN9vfdTwOf\nBly6c4mGvzTk2+3fJvoFsHVrXShLgSxczsWLuq/IxwfatUMNGsSP2yZD2htU7rw2VRfIv/wCL74I\noaFmR+IaBtUaRJGsRThy5QiTdk16/g1SELuK5CpVqnDy5ElCQkKIiopi0aJFtG7d+rFjWrduzdy5\ncwHYsWMHWbJkearVwhm1bAkbN0K2bLB+Pezda3ZEifPkVzPCPq6cz9u34datpN/e/7Q/v/79K2k9\n0jK95XQsFstzb+OIfOZKn4sNr2/gk1qfYFVWPvL7iE5LOyV6mbhEhO20XPnxmdxcJpdWq97co2BB\n3XB/4QIUL86RPm05kP4bMn1Qmz9mOL5J11nzGRMDEybA4cNQsyYcPWp2RInjrPkESOuRNm4S34jA\nEYTdDjM5IuPYVSR7eHgwefJkmjZtSunSpencuTO+vr7MmDGDGTNmANCiRQuKFClCsWLF6Nu3L1On\nTjUk8OTw0kt6MfL586GObPQlXIhS8MYbeqnTI0dsv31EdARvr3kbgM/rfk6xbOaOSnm4eTC28Vh+\n7/w7mdJkYunRpVT7sRpHryTtHfDmTXhioR4hUgZ3d/0JWCk9E9fPj/uHD9AmbyARnjCqzdvkyZTT\n7ChN4+GhB8Bq1YLz511v0xGzvFLiFVqXbO1yk/hS5bbUQqR2P/ygV3vKmFFvrFGihG23H/rXUMZs\nHUPZXGXZ12efU/Wgnfj3BO0Xt+fw5cOk90zPT61/onPZzom+fWSkXuFj3z7diiEfkIXTioyMv1co\nOFg/uf9rbRyzZQxDNw5NMT2jySEiArp21ZuNpE0Ly5bp571IujPXz1B6amnux9wnoGcA9X3qmx1S\nvGRbaiFEgrZvh48+0v8/e7btBfKh8EOM3z4eCxZmtpzpdG+4JbKXYMcbO+hWrht3o+/SZVkX/u/P\n/0v0MnFKQebMcOMGNGkCv//u4ICT6tYtOHUKol14tyvxtIgI3VRbvbqeYR6fUqUgXz6io+HCrQt8\nueVLQG8n7GzPV7OkS6cn7Pbpoz9nPLIIl0iiwlkLM7T2UMB1duKTIjmJ1q7VTzBn48x9SymRq+Xz\n8mXo1En35X3wAXToYNvtY1Usfdb0ISY2hn5V+1GjQA2bbp9c+UzvlZ5fX/2Vyc0n4+nmyfc7vqfB\nLw0S1SuXNi0sWaJ34I2M1Dlymk1HevWCypX1ZInMmQksVkyPGKaUSRNOzOmf6//8oz/d5s+vHwe7\ndsGGDXqZhnjcvQtVq0Kr/pu5G3WXdr7tnrmGudGcPp883HRk/34oV87saJ4tJeQT4ONaH1M0a1GO\nXjnqEpP4pEhOgjNnoGNHXWykoBZrIViyRM/jqVkTxo2z/fbT90xnR+gO8mbIy1cNvzI+QANZLBbe\nrfYum3ptIl/GfGw7v43KMyqzKWTTc2/r7q6f249uOrJqlQOCPHBA/6OMHw/9++sZw+XK6RHihI7f\nvx+uX9fVfI4c+hNPQl8HjBun7z8kJOXumCJ0j3H58vDtt3Dtml6a4aef9OMknnXMlIK+feHgQdi/\nriJesVn5psk3JgTu/CwWKFzY7ChcR1qPtExs/nAnvpQ+iU96kpNAKRgzBj79VF/+9FMYNSplz5IX\nqcf8+VCvXly7YqJduHUB3ym+3I66zdKOS2lfur1jAnSAy3cv03VZVzae2Yi7xZ2xjccysMbARK3I\nMWsWrFsHixfbsF661aqX4Dp7Vv80bAh58jx9XI0asGPH07/384PGjZ/+/bZt4OkJhQpBrlz6Ref2\nbT2a/KRbtyBLlofFcY4ceqZm1ap61QN390T+McIpvP22HjXu10//Gz7DtGnwzjvg5hVB7Fsv8ln7\nDnzR4ItkCtR1xMTIHglJ1XZhW1YeX0nXsl2Z336+2eE8xpaaU4pkO8imIyI1ab+4PcuPLad1ydas\n6LwiUQWmM4mJjeGzgM8Yu3UsAO192zO7zWwypcn03NsqlcgPwQMH6kXWz59/vFd41Spo1erp44cO\n1cto+PjowrdQIf3/JUuCt3ei/q4EXbumZ2ju3q1/rl7Vvy9YUBfuT4qJ0QV31qz2nVckjVIQGAgZ\nMjy3CH6W3bv1ig1RUUD71yhQayvB/YPx9rTz8ZTKTJ4Mixbpp648JWwXciME3ym+3I+5z8YeG2lQ\nuIHZIcWxqeZUTsBJwkiS1auVSpdOqTRplDpwwOxolAoICDA7BJci+dRWHFuhGIHK8FUGde7GuSTf\njzPkc8WxFSrTmEyKEaiSk0qqw+GHEz543z6lZs1Satgwpbp3V6pOHaUKFlRq6dL4j+/VSyld7iiV\nK5dS1aop1bGjUlu3OuRvSXQ+Y2OVOnNGqcWLlfrll/iP2b1bx12smFJduyr13XdKbdmi1J07RoXr\n1Ex7bF6/rtQPPyjl66vz36JFku8qNlapypX13aSt8aNiBGrhoYUGBpt4zvBcT6o7d/TTHJQqU0ap\nc0l/yTNMSsznqE2jFCNQpaeUVlExUWaHE8eWmlPGPe30YNOR8HCoUMHsaIQw3u3I2/Rf1x+ALxt+\nSYHMKXsaeJtSbdjbbQsDZ3UgcudxZvpVpkW3ETTtNOTpg2fOjH/W3unTXL2qd+IsWPCR3w8bBp98\non9p70iwkSwWPULt45PwMWfP6mn+//yjfxYs0L9v2lTvqCSMdeWK/iZh/vyHW7q+8IJetSLRX108\nzmLRq7G88vZ2Dr/4DnUL1aVTmU4GB+760qeHrVv1Q//IET2H488/oXRpsyNLWT6q+RFzDszh6JWj\nTNw5kYE1B5odks2k3UIIFzZ6tP6W354PcO+vf5+JOydS9YWqBL0RhLtbCuhltVrj77n97jsYO1YX\nKI/4rAHcGvQe45uMx8v9kYlQCxboAvFBG8R/bRF3sxWgUXMvQkP11WXLOvSvST7R0XorsgctGrt3\nw6uv6h7mJ61fD6tX69aAqlX1smPS55x4ERF6pYpr16BRI91r3Lq17jm3w5HLR6gwvQIKxb4++6iQ\nR0ZvkuraNf1Psm2bbrn44w89jUAk3rqT62gxvwUZvDIQ/G4w+TLZOBnGAaQnWQjBvHnQvbueuxUS\notf+tdWuC7t46ceXcLO4safPHirmqWh4nHY5fFhPcgsJeThJLiREFxxfxbP6xoNdVP6b/KYKFSI4\n/T2GZtjFihJWahaoyeIOi5/7Qn7jhn7z3LJF53flSqhb1yF/ofN67z2Y9MgSTxky6OXpBg7UyRHP\nt2qV7j8vWdKQu1NK0eTXJvx15i/6VenH1Fdk+SV7Pdh0ZPt2XSwXL252RCnPq4teZUXwCrqU7cKC\n9gvMDkd6kp3F0qVKBQYm7zlTYt+SM0up+Tx8WClvb91TN3160u4jKiZKVZhWQTEC9fGGjw2JK1H5\njI1V6uJFpXbsUGrhQqW+/jrhHuDp0x/2AD/606NH/Mf/+69SFy4oZbU+9uug80Eq/3f5FSNQucbn\nUgFnnh/nvXtKvfqqPl2aNEotX/78P81opj4+9+zR/zYdOihVqNDD3M+bF//xx47p3Dspw3MZHa3U\nihVKvfyyUj//bOx9J2D50eWKEaisY7Oqq3evJss5E5JSXzvjEx2t1D//mBtDSs7nmetnVLrR6RQj\nUBtPbzQ7HJtqTlkn2UEOHoTXXtM9TcuWmR2NSE1u34b27XWb4+uvJ7wp1/N8v+N7DoYfxCeLD8Pr\nxfN1e1JZrXDzZvzXLV+ue3nz5oWXXoIuXXSP77x58R9ftapeX3j8eL0e8K5deseUOXPiPz5bNt33\n6fb4S99L+V9iX599NCrciMt3L9N4bmPGbxv/zNGGdOme3nRk27ZE/P2u4sUXYdCgh+swh4fr76Pj\nW7oOYMAAve5gvnzQpo3uBfrzT7hzJ1nDdriLF/WaoIULQ9u2esOP2bMdcqqoKPj5Z72Wd0R0BB9u\n+BCAUQ1Gkd07u0POmRp5eEDRomZHkXL5ZPFhaJ2UuROftFs4iNUK778PU6boyRSTJ+t1K4VwJKV0\nXbl4sd6TYseOpM0fO3P9DGWmliEiJoJ13dbRrFizpAd1/z5s3qy3qdywAU6ehJdf1gXVkwIDoUED\nyJ798SXRqlXTf5iDWWOtfBbwGWO2jgGgnW87fm7z8zOXiVMKvvxSr+T2669P1d/igQ4dwN//6Q9I\ne/fqNg1XsG+fnngXE6MvFy+uW3969tQf0Az2wQe6g6hPHyjQbTSfBXxGuVzl2Nd3Hx5uMi/f0ZTS\nH1CkFf/57sfcp9y0cvxz7R/GNxnPRzU/Mi0W6Ul2ErLpiDDDggW6LTQwMOGN2J5FKUXzec3589Sf\n9i8Ev3evbtZ9MHv/gTp1dOH8pKgo/ZMhQ9LPaYCVwSvpuaInNyNvUiJ7CZZ3Wk6ZXGWeeZvYWCmQ\nnys2Vq+c8WBS4IEDuqc8vslqjRvriW0PJgZWqKBX33BmViv4+upPqP366U1kHPSgWLwYOnfWqVuy\nNpyuOwsTERPhdGvSurIxY/SXV/Pn62+WxLOt/2c9zec1J71neoL7B5M/U35T4pCeZCcze7ZS7u5K\npU+v1KlTjj1XSu5bckYpNZ8REUm/7by/5ylGoLKMzaIu3b6UuBtFR8f/+3v3dHN0xYpKDRmiAn74\nQf8uBTj570lVflp5xQiU95feav7f880O6Skp9fH5XGFhT/eZe3rqdacTeqzZyaZc7tun+9vjc/++\nIfE8y7FjSmXIoNMycaJSXZZ2UYxAdVzc0eHnTiyXfWz+58oVpbJm1f8GtWsrde2aY8/nKvl8deGr\nihGozks6mxaDLTWnjHskg//9T89+X74cihQxOxqRGqRNm7TbXYu4xgfrPwDgmybfkDtD7oQPDgnR\n+9+2aqXbI65ff/qYdOngwgXYv1+vNlG+fIoZcimWrRhBbwTxevnXuRd9j9eWv8Z7694jyhqV6Pu4\neFEvIyVslCOH7hWaNAl69NCjszEx+huJ+LY1vX1bD+edPPlwG26jRUTAL7/oXvnKleGnn+I/zsGj\n3Xfv6s6VO3d0B1L5VptZeHghaT3SMr7JeIeeWzyUI4de3SZ/fr2mcp06eqNN8WwTmk4gnUc6Fh1Z\nxMYzG80O57mk3UIIEefNVW/y0/6fqFuoLoE9A+PfenrMGN18e+zY47//4w9o0SJ5Ak1GSimm75nO\n++vfJzo2mhr5a7Ck45JELRNXt66u7davf2LTEWG7W7f0p474lksLCNCtDaDX5KtSRbdoNGgATZrY\nd97QUPj+ez1D7sEnnsyZYcgQPak0mYWH66Wrr1+HHTut1Jv/IgfDDzKi3giG1zdwgq1IlPPnoVkz\nOHpUF8wbNujPdCJhX27+kmEBw/DN4cuBtw88vjZ9MrCl5nSakeR70feef5AQ4inh4cbcz6aQTfy0\n//sx5aMAACAASURBVCe83L2Y0XJG/AUy6NG6Y8cgUya9jMZPP+nRYhcskEG/oPar2o8t/9tC/kz5\nCQoNovLMygScCXjm7e7d04Oax47pHbsOH06mgF1VpkwJryecJo3+RiNPHv3pxN9ff5ibOTP+420Z\nlLl4Eb79VhfIL76oH+9hYaYUyAC5c+v5Bn5+sODELA6GH6Rg5oJ8XOtjU+JJ7QoU0CPKtWvrh5XJ\n0ylShI9qfkSxbMU4dvUYP+z4wexwnslpiuQqM6uw/+J+s8NIdvPnw99/G3d/gYGBxt2ZcPp8rl2r\nV5r68Uf77icyJpK+a/riaYWZ3l0pNfZHWLgw/oM/+EC/S1+9CkuXQu/eelm1RHD2fD5L9fzV2ddn\nH42LNNbLxP3amHHbxiU4IvHCC3puYp06+jNEQnMV7ZGS82momjX1xhxhYXpob/lyveVzpwS2ZJ42\nTT9xOnWCceMgIIDA+FZbAT0iPXy4nqG1Z49+vJu85biXF3hnv8awjcMA+Pblb/H2dKJt0Eldj81s\n2fQI8qZNumh2BFfKZxqPNExqrjciGrlpJKG3Qk2OKGFOUyQfu3qM6j9W5+utX2ONtZodTrLYvFm3\n29Wtq59cQtgiJETvqBcRYedo8pUrbPikA2OnHOf6ODd6DvxFj5z98kv8x5cvD/Xq2b19bkqUM31O\n1ndbz6d1PiVWxfKJ/ye0X9yem/fjX/c5a1b95tmunR7gfPll/bWscBCLRX/n/eqrel2+jh3jP27/\nfv0EWrJEjwg3bAgtWyb8jzNihC6WncjwwOH8G/Ev9X3q0963vdnhpHrp0slayrZoVqwZ7XzbcTf6\nLgM3DDQ7nAQ5TU/yu3+8y5TdUwCoV6gec1+dS8HMrt3Ed/++3uxh6VI9MjBvnp6QIcTzREbqr/f2\n7IFXXtGDaEldaSpk1a/4tOnx8Bdly+rWiVat9ElEvFYfX83rv7/OzcibFM9WnOWdl1M2V9l4j7Va\n9Z4nsbEwfbosA2m6mBjdB7Nnz8Pl6A4e1N+SjBtndnSPiY3V/330+X0o/BCVZlRCoTjQ9wDlcpcz\nJzjxXA8qLHnOP+3czXOUmlyKiJgI/F/3p1GRRsly3hS7TvLak2vpvbI34XfDyZwmM9NemUbXcl3N\nDs+hntx0ZNIkePdds6MSzq5fP11s+fjo/QuyZn3ODUJDISjoqZG1WBVLox/r0mfiNqLr1abHJ/Md\n932hCzp17RTtF7fnYPhBvD29mdlyJt3Kd4v3WNl4wMlFRuo1ujNmNDuSx4weDdu367my2bPriaSN\n5jYiICSA/tX6x31tLZzTkCF6kuWUKfLcj89XW77i042fUipHKQ6+fTBZJvGl6HWSw++Eq1bzWylG\noBiB6rasm7oecd3E6BwvNlapL7/U6y1mzarU5ctJvy9XWUvRWThjPs+d02ukenkptWdPAgdFRSkV\nEKDUoEFKlSv3cK3ZsLDHDpu5Z6ZiBCrX+Fzq2j0HL/SpnDOf9robdVf1/L1n3GtW/7X9VWRMZLKc\n2xXzaRZnzKWfn1IWi/7x89O/W3JkiWIEKvvX2dW/9xJYq9kJOGM+k9upU0qlTatfetu2tW+JeFfN\n5/3o+6r4xOKKEaivt36dLOe0pfR1mp7kB3Klz8XKLiuZ0XIG3p7ezDs0jwrTK7D5rMEzXpyIxaLn\nmMyeDatXQ86cZkcknFmBAvrb4d9+05Pt4/Vg+atx4+DQIUifHtq00cto/efSnUsM8h8EwA/NfiBr\nuucNR4v4eHt683Obn5n+ynS83L2YvGsy9ebUS/RklNDQpzckFCI0FLp21Z9uP/9cb0B4L/peXP/m\n6IajyZbO+K2uhXGKFNGrkGTJAitW6NUIZd30xz06ie+LTV843SQ+p2q3eNKJf0/QbXk39oTtwYKF\nT2p/wsj6I5N9TT0hnFJ0tP56OH36p6/r31+vHdu8uf6pXfupTQ66LuvKwsMLaVasGWtfW5vwkm8i\n0XZd2EWHxR04f+s8Ob1zsrDDQhoWbpjg8VeuQK1a+oPx6tV6lrwQUVFQv77ukHr5Zb2Kjbs7jAwc\nyYhNI6iQuwJ7++zF3U2+v08JjhzRaymHhkLp0vDnn3p+qXiow+IOLDu2jE5lOrGowyKHnivF9iTH\nJ9oazchNIxmzdQyxKpbKeSszr908SuUolcxRCuEEwsL0zhRr1+ohiiFDYPDgp4+Ljn7m6hPrTq6j\nxfwWeHt6c7jfYQpnLezAoFOXq/eu8tqy1/A77YebxY2vGn7FoFqD4v0QcuKEHiE8f15vQCCbjgiA\nH37QcwgLFNBzDnLkgLM3zlJqSinux9xnU69N1C1U1+wwhQ0ebDoSGal7zHPlMjsi53Lu5jl8p/hy\nL/oefq/70bhIY4edK0VuJpIQT3dPRjcczaZem/DJ4sO+i/uoPKMy03ZPSzW79P3yC5w9m7hjXWkt\nRWfgNPncuBEqVoR8+eCNN2DZMt06ceRI/Mc/o0C+G3WXfn/0A2Bk/ZHJWiA7TT4dKId3DtZ1W8ew\nusOIVbEM/mswry56Nd5l4kqU0KOFZcsmbdOR1JDP5OJMuXznHfj4Y1i8WBfIAB/7fcz9mPt0LtM5\nRRTIzpRPZ/Bg0xE/v6QVyK6ez4KZCzKsjl73u//a/kRZo0yOSHP6IvmB2gVrc6DvAXpU6EFETATv\nrH2HVgtaEX7HoO3GnNTq1dCrl37zPHTI7GiEGbZvh6UbMuklqry99XquU6bA6dN6yruNRmwawdmb\nZ6mYpyIfvPSBAyIW7m7ujGowitVdV5MlbRZWHl9JlVlVOBT+9JM4X77HNx2pXVv/V6Renp56OsFL\nL+nLgSGBLDm6hHQe6RjfZLy5wYkky5ZN72Ej4vdhjQ8pkb0Ex/89zoSgCWaHA6SAdov4LD6ymL5r\n+nLj/g1yeudkdpvZtCzR0oERmufGDWjbVm82kjkzrFyp93EQLiQmBnbsgHXrIDhYjxL/5/JlqFwZ\nwi7E4jd4I42G14a0aZN8qv0X91N1VlUUih1v7KBqPufaIMEVnb5+mnaL2nEw/CDpPNIxs9VMupfv\n/tRx9+9Dt256QwInW6pXmCgmNobKMypz6PIhRjUYxbC6w8wOSRhMKVlH+YENpzbQ9LemeHt6E/xu\nMAUyG78kqUv1JCck9FYoPVf0ZOOZjQC8XeVtp9ya0wiy6YgLUgrmztW9xRs26E9DD5w6BUWKYLVC\n06bw1196cldAgH2b3FljrVT/sTp7L+7l/erv8/3/t3ffUU1kXxzAv6GJgIC9AIqdIgqIXbGAvfde\nERRE17Kurq4/sWHvgn1dO/aCiooFRUBYRBQVG4o0UVFAmkCS9/tj1gooZcIkcD/ncDyByczlOiE3\nM/e9121D0X8Pki8Z2dzdr39C/wEAODVzwvqu63MMQpZIuDfLwi4MQ0oetyA3OHs5w1DXEI+cHqGs\nalmhQyI8YoxbG6FqVW4WEyqWvw7iG2wyGEcHH+V9/yWqJzkv+tr68B7tjTWd10BNWQ3bgrfBcrsl\n7sTdETo03qmrAx4e3AspK4tbfCQtLfdtS3rfUnGTWT5FImDNGq7pMCmJa0797Tdu5JaeHgBuJdyr\nV7n+tSNHir4K9Oagzbjz+g4MtA2wpOOSov8OhVBaz8+yqmXxd5+/sb3Xdqgpq8H9X3dY77FGdHL0\nd9spKxesQC6t+ZQFIXPp6QkkJOT8/vv091hwfQEAYF2XdQpVINO5mT937wLbt3N/7x0duQ/KuSlN\n+VzXdR00VDVw7NExeEd4CxqLwhbJAKAkUsKs1rMQNDEIJpVN8OT9E7Tc3RLLfZdDIs3jTFNQysrc\nanyrVgHnzuU+6xeRM2/ecFeLHz/O/eezZnH/qc+fA0+eABs2cJeOy5TBhQvcSltKSsDhw1/q5kKL\nSo7CX9e427RuPdxQrox8rSpWGohEIjg0dcCt8bdQU6cmAmMDYbnDEldfXP3lc6Ojua4cUvIEBQED\nB3JtVT/Oobvg+gIkfkqETW0b9DPqJ0yARKYsLYGTJ7mLYdu3c3eJMzKEjkpYNXVqYoE19+HQ2csZ\nmeJMwWJR2HaLH2VkZ2Du1bnYFLgJANCuZjvs678PhrqGPERISD5IJNw7npcX10Zx57+7GgsWAIsX\nF2hXT59yfyyHDeMWmikKxhj6ePTBuafnMNB4II4POV60HZIiS0hPwMiTI3E54jKUREpY2nEp5rSd\nAyVRzusWMTHcwF0LC+4Dk0bJ6ygrtd6/54qkqCjuTuGWLV9/di/+Hix3WEIEEe5NvgfTKqbCBUpk\nzs8P6N2bW8K6TRvu7kL5Ury+U5YkC423NsaT90+w3GY55rbNZarTQioVPcl5ufT8EsadGYf41Hho\nl9GGWw83jDQbSQslENn7PLnpZ+rq3IoAEyYAgwcXeHcZGdz6H0XtTz3+6DgGHxsM7TLaCJ8Sjhrl\nahRth4QXEqkEi24swpKbXOtLn4Z9sLffXuiq63633b//cvOrfvjAFcu06EjJIJVyE9V4eQHNm3Mz\nnHxe74cxhg57O+Dmq5uY1mIaNnbbKGywpFh8XnREV5ebLk5X99fPKcm8I7zR5UAX3gfxlYqe5Lx0\nrdcVYY5h6GfUDx8zP2L0qdEYcXIEEjMShQ5Npvbs4XrafHx8uPuyEgk3IoAUSY4+MIkEiIzMfeMu\nXbipCZyduSvJ799z74CFKJABoGzZohfISZ+SMM1rGgBghc0KwQvk0tRX9yvKSspY3HExzg0/B111\nXZx9chZWO6xw/83977Zr1gy4dYubZ9Xfn5siLiqK+xnlkz/Fnctly7g/DxUrAseOfb8g5rFHx3Dz\n1U1U0qgEl/YuxRoXX+jcLDhTU+417uWVs0AujfnsXLczBpkMQnp2OmZenilIDCqCHFXGKmlUwskh\nJ/H33b/x28Xf4PHAA35RftjXfx86GHYQOjze7d/PXaxs2IBhc9o4IPablUeUlbmvBw+A+vVzPrlD\nB+4dV1kZUFH5+u+ZM0CtWjm3t7MD4uO/7vfzc9avB6pXz7m9qyt3CezH7Z2duXeHH3l4AKmpObfv\n0QMol0sf7b//cqMZf9y+YcMcyzAD4ArXz3n59vdVUcl7WPG7d9w6ol5e3L9ly3I5+3F7IyPg2TO5\nGp7859U/8Tr1NVobtMYkq0lCh0Ny0bNBT9xxuIOBRwciND4ULXe1xPZe2zG6yegv2xgbc4uOdOvG\nvZQLuugIkT8ZGdyH4IMHv19lMS0rDb9f/h0A4NrJFeXLluJ77qWQAf8znim0dV3W4cKzCzj+6Dgu\nR1xGl7pdijcAJgdkGcbThKes+c7mDC5gIhcRm315NvuU/Ulmx5MZqZSxqCjGTp9mbOFCxuLjv/wo\nNpYxMzPGAMaCVVswqUjEmEjEfePz1/Pnue+3du3vt/vV9nXq5L79s2e5b1+3bsG2l/X+89r+6dPc\nt2/aNGcua9f+Lv98ePeO+y/mk1+UH4MLmMpiFRb2JozfnRPepWels/GnxzO4gMEFzPGcY46/VR8+\nMNauHWNLlggUJOHVkyc5v7fg2gIGFzCLbRZMLBEXf1BELvH9/qBIVviuYHABa7C5AS/1W0FqzhLX\nk5ybbEk2lvouxdKbSyFlUphXM8fBAQdhUtlEZsfkjZsbcPYsEBLy/RxBZ89yXf7/+bzoyJMbr5FV\nrhKmzVLFX/OkUIaEaxFQU8v93n1MDLeYvETytU1DLAZMTHK/EnvjBnel9/O2n7fv3RvQ0sq5/e7d\n3EiEH7efPj33xsr587lZIX7cfvNmbiLJHw0dyv0OP27v5QXo6+fc3sqKW6nux9/36VOgTp2c29eu\nDcTFcVfcu3fnvho04PVq8ceP3C11S0tg587c01hQWZIsWG63xMN3DzGv3Tws67Ss6DslMscYw66Q\nXXD24pZlba7XHMcHH/+uFy8ri5sOUI5uWBCevEx8CWM3Y2RKMnFr/C20qdlG6JCIHJBKubvFGhrc\nNHFmZkJHVLz4HsRXoJqzyCU5D4orDP8of1ZnYx0GFzD1pepsc+BmJhX641l2NmMPHnCXg3MzfvzX\nK5gVKjBma8vYH38wFpbzymBGBmPt219nAGMtWsg47lLi+tmzXGJlRCplbPBg7r/XzIyxtDR+9rv0\nxlIGF7B6m+qx9Kx0fnbKg+vXrwsdgkL4N/ZfVmt9LQYXsEqrKjHvCO9ct/s2n6mpjB09ylhmZjEF\nWcLIw7k54MgABhewESdGCB1KkclDPkuKy5cZA65/KQVatWJszx7+3i8UgXeEN4MLmMYyDfYq6VWR\n9lWQmrPQw4I+fPiAzp07o0GDBujSpQuSvl0x7BuGhoZo3LgxLCws0Lx588IejhetDFohdFIoxpuP\nxyfxJ0z1mooeh3ogPjW++IJ48YK7ujplCtCqFaCtDTRqxM2nmxsHB+DUKW6wWEIC4O0NrFzJPecH\n6urAwoXcJq6uue9O+PsGCqZcuSItA/0rmzZxg3bKleNWVORjeq9n7599mTFhW89tCrUAAeFY1bDC\nHYc76Fq3KxLSE9D1QFe4+rpCyqR5PufoUWDIEO4Gypw53PTbRH7ktUjEZ1dfXMXJ8JPQVNXEKlta\nl5x81bkzsGsX4OTElQwBAcD48YC1tdCRFR/bOrYYbDKYG8R3qRgH8RW2Ep89ezZbuXIlY4yxFStW\nsDlz5uS6naGhIXv//v1P91WEMArt2MNjrPyK8l+u1JwOP83vAfK6Qr16dc5+WENDxtas4ff4efj9\nd8b69GHs3DnGxNTuJig/P8ZUVLhT4MQJfvYplUpZp72dGFzAxp4ay89OiWDEEjFbeH0hE7mIGFzA\neh/qzRIzEnPd9sQJxho1+v5Pi60tY7duFXPQJIfwcMaMjPL+v8iWZDNTN1MGF7BlN5cVb3BEoaSm\nMrZ7N2PNmzO2YoXQ0RSv6ORoprlMk8EF7OKzi4XeT0FqzkJXpw0bNmTx/w1eev36NWvYsGGu2xka\nGrKEhISfByFQ10dMcgyz3Wf7ZaCM/Vl7lpqZWvAdJSUx5uPD2Lp1jI0axZipKdcmkZuAAMaGDWNs\n1SrGrlxh7BcfIPgkFjNWrdrXN1ADA8YWLWIsOrrYQiDf6NmT+3+YOZO/ff5z9x8GF7CKKyuyd2nv\n+NsxEdT5p+e/fKivu7EuC30dmut2Uilj/v6MjRvHWNmy3Pl14UIxB0u+k5LCmIkJ938xalTu22y6\nvYnBBazOxjosI1t27V2kZMnrQtfly7l2ZJYIK2+tZHABq7+pfqEH8RWk5iz0wL3y5csjMTHx89Vo\nVKhQ4cvjb9WpUwc6OjpQVlbGpEmTYG9vn2MbkUiEsWPHwtDQEACgq6sLc3NzdOjQAcDX+QFl8VjK\npJjqPhU7QnZAXFOM+hXqY0b1GTCuZJy//V27Bh8bG+7xf7+PD/eLo0NEhMzj//Hxt3Mp5vbzt2+B\nBQt8cO4cEBfH/VxZ2QcnTgB9+8o+PkV7/Kt8FuVxixYdsGkT0LSpD1RUir4/02amMHYzxvvw95jb\ndi6W2y2XeX4K+liW+Szpj2s1qYWBRwfi7u27UFNWw86pO1Ez8evcYT9ub27eAceOAXXr+kBJKefP\nra07QElJfn4/oR9//h6f+2cM6NzZB1evAsbGHRAUBAQHf7/96YunMfrUaKTWSMXpoaehE68jF/mQ\nx3yW5sefv/er7a9d88HIkUB8fAe0agVYW/ugQwegWzf5+n0K+9j7mjcmnp2IqPJRcO3kilaSVr98\nfmho6JeW4MjISOzdu5efgXu2trasUaNGOb7OnDnDdHV1v9u2fPnyue4jLi6OMcbY27dvWZMmTdjN\nmzeLVNXLyv34+6yReyMGFzDlRcpsic9ilh0VydjZs9zl1rwu971+zZiaGmNWVow5ODC2bRtjQUEy\nHez1M/kdLCGRcBeyhwxhrF8/2cakyBRp8Mnok6MZXMBs9toIPyA1D4qUT3mUkZ3B7M7Yfbn71Wd5\nn0JdeXz7ljF9fW4McF6zJZY2sjg33dy4K8iamow9epT7Ng6eDgwuYF32d5Hb121h0GudX/nN58eP\njDk5Maat/fWusY4OY87OJWdQ7+dBfGWXli3UIL6C1JyFvpJsZGQEHx8fVKtWDa9fv0bHjh3x+PHj\nnz5n0aJF0NLSwqxZs777vqyngMuvT+JPWHR2FqznuMPyNVA17ZsfqqlxU5+pquZ8YnZ27t9XEFJp\n7rPDhYUBr15xs54pKxd/XCT/rry4gs77O0NdRR1hjmGoV6Ge0CERGdodshtTLkxBpiQT9SrUw+bu\nm9GtXrd8P3/nTm5M8Gc2Ntzjfv24P3Wk6OLiuFklMzOBw4eBYcNybnP39V003dEUykrKuD/5Powr\nGxd/oKRESkvjBvNu3w4EBnLTjAYFCR0Vf4YeH4qjD49igPEAnBhyokDPLZZlqfv06YO9e/cCAPbu\n3Yt+/frl2CY9PR0pKSkAgLS0NFy+fBlmQk/wJ5EAjx8Dhw7lGG6srqKO5f23wPZdOVRNAxLVAZ+6\nyngwpjvYnj1cNZkbBS6QgbyXPl69mpv+uHZtYNEibjpiIn8ysjMw+dxkAMAC6wVUIJcCdpZ28Jvg\nB5PKJnj+4Tm6H+yOQUcHITo5Ol/PnziRW/523DhuAcmrV7kpx3+4fkGKoEYNbsaaefNyL5AZY5jq\nNRUMDFObT6UCmfBKU5ObAeP2bSA0lFsUNzdicfHGxZe1XdZCU1UTJ8NP4tLzS7I7UIGvU//n/fv3\nzMbGhtWvX5917tyZJSZyI65jY2NZjx49GGOMRUREsCZNmrAmTZowU1NT5urqWuRL34Vy6BB3r6F1\na8Y0NL7egwgPz337GzfYh4d32ECPAV9uaw4+Opi9Ty++QXaFxectrg0bvl+gTkmJmxkjIoK3Q8g9\nvvKZkcHYvHmMJSfzsrvv/HnlTwYXsEbujVimWL7vp9EtWH55X/Vmq26t+jLiW3OZJlt5a2WBzoPE\nRMY2b+ZmxggMlGGwcq64z81D9w8xuIBVWV2FJWUkFeuxiwO91vklq3zOnMlYy5aKOe/yqlurCjWI\nryA1p/DNwIynIjk9nbFPeSSpTZvv50WqWZNrxP3F8E+pVMr23N3DtFy1GFzA9NbqsasvrhY9Vhni\n+4X0be+yqipj6urc0rilBV/5nDSJO/W6deNld1/cj7/PVBarMJGLiPlH+fO7cxmgN05+fc5ndHI0\nG3R00JcP9cZbjNn1l9cLtC+pNO+ZK3fsKPm9y8V5bqZmpjK9tXoMLmC77uwqtuMWJ3qt80sW+ZRK\nGatTJ2fv8v37vB9KJjLFmcxoi1GBp04sSM2pmMtSp6Rw9w9CQr5+hYdzqzL0759z+3/+4ZY6trQE\nLCyASpUKFF/EhwiMPjUaATEBAIBZrWZhWadlKKOSy7LNJdjbt1xPU69eOX/2uXOFepdz2r8fGDOG\nW+Xb3587DfkgZVK0+bsNbsfchlMzJ7j1cONnx0RhXXp+Cc5eznj+gVtJZKTZSKzpsgbVtKoVep8R\nEUC9/zp4qHeZH39d+wvLfJfBqoYVAicGQklU6M5HQorkx95lgFvyPiaGaxmSd1dfXIXtfluUVSmL\n8CnhqKVb65fPKUjNqZhF8qRJwI4d339PSQlYtw747Td+g/uPWCqGq68rFt9YDAmToHHVxjg04BBM\nq5jK5HiK5swZwNmZ63W0s+NW/SLAgwdA8+ZARgb3R+jbwVJF5f6vO6ZcmILqWtURPiUcOuo6/O2c\nKKxP4k9Y7bcarrdc8Un8CdpltLG041I4NnOEipJKgff34gWwZAlw5Ah3HgNA5cpc//KcOTwHr8Ci\no7mFUdu1+/l2ER8iYOJugixJFvwn+KOVQatiiY+QX7l3jxvUm5AAeHgIHU3+DTs+DEceHkF/o/44\nOfTkL7cvUM1ZkEvbsvIljLg4bim4xYsZ69+fMXf33J+wezdjlpaMTZzIbXP7drE10wREB7C6G+sy\nuICVWVKGbby9kUmkkmI5dn4IdYvLzu773uXevRnz9FT8Vf2Kks/kZMYaNOByMmZM3reyCyMmOYZp\nL9dmcAE7/vA4fzuWMboFy6+f5TPiQwTrebDnlxYM823mLCA6oNDHSkxkbMsWxszMuHP6jz8KvSu5\nVJRzMzOT6+tUVmbs+C9ejn0P92VwARt9cnShj6cI6LXOL3nI59278tm7HJMc82Vchtczr19uX5DS\nV37u8VSvzl3b79UL+N//gFOngCtXct92wgTgzh3uI4+jI9CiBaChUSxhttRvibuT7sLOwg6Zkkz8\ndvE39DjYA69TXhfL8eXVjh3cCPkhQ7iWC09PbmaM48eFjkw4ampAx46AmRmwdSt3C4sv0y5Ow8fM\nj+jTsA8GGA/gb8ekxKhTvg48h3vi9NDTqKVTC6HxoWi1uxUmnp2IhPSEAu9PVxeYMoW72hQQAEyd\nmvt2n682lya//87NIlCjBtC+fd7bXY64jDNPzkBLTQsrbFcUX4CE8GDDBm7GjBo1uDvHYWFCR8TR\n09bDwvYLAQBTvaYiU5zJ277lp90CALS1uYbNz1/NmwP16wsdXp5Ohp+Evac9PmR8QMWyFbGz9070\nN86lJ7qUefuWawM/cQK4cQNQVxc6ImGlpADlyvG3v7NPzqKvR19oqWnhkdMjGOgY8LdzUiKlZ6dj\nme8yrPZbjWxpNiqUrYAVNitgZ2nHaz8sY4CVFVdQT5pUOnqXPTyA4cO5mUB9fblrNrnJlmSjybYm\nCE8IxwqbFZjTlnpViGI5eBDYvPlr7zIAtGwJ7NoFmArceZolyYL5NnOEJ4RjacelmG89P89tFbMn\n+flzbkLevCbtlVNxKXEYf2Y8LkdcBgDYWdhhQ7cN0FLTEjgy+ZWWxrWPjxsHGFB9VyApmSkwcTdB\nzMcYbOy2EdNaTBM6JKJAHic8hvMFZ1x9eRUA0EKvBdx7usOyOj+jSV+9AoyNv+9dHj8esLf/Oviv\nJAkP5xZpSEvjigdn57y33XB7A2ZcmoF6FerhgeODUjfwm5Qc9+5xd48PHOAWy4mLAypUEDoq9jlL\nWwAAIABJREFU4NrLa7DZZ4OyKmXxaMojGOoa5rpdsSwmwru6dRWuQAaAGuVqwGukFzZ224gyymWw\n++5umG8zR2BM4K+fLAPfrvMurzw8uI4aQ0OuJePcuRzrusgNecvnX9f/QszHGDSr0QxTmk0ROpwC\nk7d8KrqC5tOokhG8R3vDY6AHqmtVR2BsIJrtbIapXlOR9CmpyPHUqsW9YW7ZwrUZvXsHrFoF9OzJ\nXWWWZ4U5N9+84WatGTaMa0XJy9u0t1jow90OXt91fakokOm1zi95ymeTJoCbG/dav3Qp9wJZIgHS\n04s3rk61O2FYo2HIEGdgxqUZvOxT8apSOaQkUsK0FtMQ7BCMxlUbIyIxAm3+boPFNxZDLFXQ5Wxk\nyMzsa+/yuXNcoWxoyM3gp8gkEuDjR9ntPyg2CJsDN0NZpIwdvXdAWYnm2yMFJxKJMLTRUDx2fowZ\nLWdABBG2BG1Bwy0Nsf/e/oLNNJSLH3uXx4/nrrDy2ZMvLzp04GYg3bnz57/f/Gvz8THzI7rX646e\n9XsWW3yEyJKmZt49+BcuCNO7vKbzGmipaeH049PweuZV5P3JT7uF8GHwIlOcifnX5mNtwFoAQCv9\nVtjffz/qVqgrcGTy5+1bYO9e7rbN8+fAxYtA165CR1V4f/3FTZN14gTQuDG/+86WZKPZzma49+Ye\nZreejVWdV/F7AFJqhb0Jg9MFJ9yKugUAsK5lDbcebmhUpZHMj+3hwRWX/fuX3N7l4LhgNN/ZHMpK\nynjg+AANKzUUOiRCZG72bGDNmq+PW7bkpkAdOlT28yys8V+D2d6zUbd8XTxwegB1le8HRilmT7Lw\nYfDq6ourGHt6LGJTYqGlpoVN3TZhnPk4iEri5ZQikkqBmzcBa+vcO26SkrirU/Ls/HluYhYlJW5S\nlo4d+d3/ar/V+OPKHzDUNcQDxwfQVNPk9wCkVGOMYd+9fZjtPRvv0t9BRUkF01tOx8L2C2U2vkIq\nBerU4fqYK1fmxijY28v1WO0CY4yhzd9tEBATgN9b/47VnVcLHRIhxeb+fe4i2P79X++y7tsHjB4t\n2+NmS7Jhvt0cj949wpKOS/CX9V/f/Vwxe5JLGJs6NrjveB+DTQYjNSsVE85OwOBjg/E+/b1MjytP\nfUv5paTE3bbMrUB+84abHbB3b25aOXExd6/kJ5+RkV9f9EuX8l8gv0x8+aWfcWvPrQpdICvi+SnP\n+MqnSCTCWPOxeOL8BI5WjpBIJVjjvwZGW4xw7OExmVzEEIu5q02fe5dXrwYaNOBW9UtL4/1wv5Sf\nXEqlBdvnwbCDCIgJQFXNqlhgvaBwgSkoeq3zSxHz2bgxNz4hLg7Yswfo0gUYNEj2x1VVVsWW7lsA\nAMt8lyEyKbLQ+6IiWYYqlK2AI4OOYG+/vSinVg4nwk+g8bbG8I7wFjo0hREQwPX6njsH9OnDTYDi\n4sKtbiUPPn3iXvSJidyVZL5XIGOMwfG8IzLEGRjeaDi61evG7wEI+Ub5suXh3tMdQfZBaFajGWJT\nYjHk+BB0O9gNz94/4/VYamo5e5fLlgVSU7leR3kjlXJ/g1asyF+xnJKZgj+8/wAArLBdAe0y2jKO\nkBD5pKnJ3Sm6dIl7jf8oPR2YMYO78syXjrU7Ynij4fgk/lSkQXzUblFMXia+xKhTo+Af7Q8AmN5y\nOpbbLM/RK0Ny+rF3GeCWv965U9i4AK4Hedgwrni/cwcoX57f/R8OO4wRJ0dAV10Xj6c8RlWtqvwe\ngJA8SKQS7AzZiXlX5yHxUyLUlNXwR5s/MK/tPJRVzeWdjgdJSdxVJxOT3H9Wtiw3m4QQFi8GFi4E\nKlYEHj4Eqv7ipfjn1T+x4tYKNNdrjgC7AF7noyakJPnnH+5DMsDNMz5pEje4v6gfluNS4tBwS0Ok\nZqXi/Ijz6FG/BwDqSZZbYqkYK26tgIuPCyRMArMqZjg44CDMqpoJHZpCkEoBHx+uWP79d27RAnlw\n8CA3N6wlP1PNfvEh4wOMthjhXfo77Oq9C3aWdvwegJB8eJf2DnOuzMGe0D0AgNq6tbGp+yb0atCr\nWOOYOZPrbfw873Jx9i5fvgx0++8mzsWL3G3jn3n+4TlM3U2RJcnCbbvbaKGfxwojhBA8fsy1ZXzb\nu6ytDaxdy10QK4q1/mvxu/fv3w3io55kOaWipIK/rP+Cv50/6leoj7C3YbDaaYX1AeshZQVsdsuD\nIvYt5ZeSEtCpEzciPq8CecECfnuX85PPkSP5L5AB4A/vP/Au/R2sa1ljgsUE/g8ggJJ8fgqhOPJZ\nWbMy/u77N26Nv4XGVRvjZdJL9D7cG309+hap16+gQkOBhITve5ePHOEWM+BDXrmMjgZGjODmeV64\n8NcFMgDMvDQTWZIsjDMfV2oLZHqt86sk59PI6Gvv8t9/czNhfPwI6OkVfd/TWkyDSWUTRCRGYLVf\nwQfOUpEsgOZ6zREyKQT2lvbIkmRh5uWZ6HqgK2I/xgodmkJ7/pwbOCePvcsFdSPyBnbf3Q01ZTVs\n77WdZkUhgmtTsw3uONzB+q7rUU6tHM4+OQsTNxMsu7kMmWKeKtWfuHoVuH37a+/ytWvcctBv3sj2\nuE5OwPv33PSUC/Ix9s7rmRc8n3qinFo5LLdZLtvgCClBNDW513dAANefnNcH0piY/O9TVVkVbj3c\nAACut1zxMvFlgWKidguBnXl8BhM9JyIhPQEVylbAjl47MNBkoNBhKaTERGD3bq4d49l/Y4yUlLgB\nA7t3CxpagWSKM9FkWxM8ef8ELu1dsLDDQqFDIuQ7cSlxmHV5FjweeAAAGlRsALcebrCtY1ssx09O\n5tqcnj0D1q+X7bGio4FZswB3d6BSpZ9vmyXJgtlWMzx9/xSrO6/G761/l21whJQy794B+vqAhUXB\nepdHnBiBww8Oo2/Dvjgz/Az1JCuS+NR4jD8zHhefXwQAjDMfh03dNqFcmXICR6aYGPvau3ziBPcG\nt5ynCzp793LT2lhY8LO/3Lj4uGDRjUUwqmSE0EmhpWIJW6KYrr64CmcvZzxOeAwAGGI6BOu6rIOe\nNg/3SYsgMBA4fpxbvKA4e5c/9z82qNgAYY5hUFMuoSukECKQq1eBgQO5D8oA17s8ahRXMP9sEa9v\nB/HBBdSTrEiqaVXDhREXsLn7ZqirqOOf0H9gvt38y0wYBVGS+5bySyTi5io+fBiIjeUG/OQmKurX\nvcvf5tPPjxtE0Lo191xZCH8XjuW3uIp+e6/tJa5ApvOTX0Ln06aODe5NvoflNsuhoaqBow+PwsjN\nCOsC1iFbki1YXG5u3GpfDRp8Hcfwq97louYyPjUei24sAgBs6Lqh1BfIQp+bJQ3lk2Nj83Xe5c+9\ny+7uwMaNP39ejXI1sKjDogIfj4pkOSESieDc3Bl3HO7AvJo5XiS+QLs97bDQZ6GgbzaKrnJl7utH\njHELlOS3d/ntW+62jljM9SjWrMl/rFImxaRzk5AlycJEy4mwrmXN/0EI4Zmashrmtp2LR06P0M+o\nH1KzUjHr8ixY7rCE7ytfQWJydgYmTOCWv71+netd1tcHbt2S3THnXZ2HlKwU9GrQC93rd5fdgQgp\n5TQ0uDbKgABunnVnZ8DR8dfPm9p8KkwrmxboWNRuIYcyxZn4n8//sNpvNRgYWui1wIEBB1CvQj2h\nQysx3r3jrgh/nndZSQno0YO7Pduz5/er/0kk3ACCa9eAtm25f1VV+Y9pV8gu2Hvao4pmFTye8hjl\ny/I86TIhxeDCswuY6jUVLxJfAADGNhmLVZ1XoYpmlWKP5XPv8vbtwNOn3BWon81lzhhw+jTQt2/u\nK4DmJSg2CC12tYCqkioeOj1E/YolaG1tQhTYzJlAo0bA0KFc7/KNyBvoULsD9SSXBNdfXseY02MQ\n8zEGmqqa2NBtA+ws7GimA5587l3evh04eRLIzuZuzz5+zLVsfDZ/PuDqClSpAty9C9SowX8s8anx\nMHYzRtKnJBweeBjDGg3j/yCEFJOM7Ays9FuJFbdWIFOSCV11XSzrtAyTmk6CspJyscfDGPDyJVCn\nTs6ficXcz+rX59o0nJ25N1QPj/ztW8qkaL27NQJjAzGnzRyssF3Bb/CEkEKJiADq/Xdt8XPvsoMD\nYG5O8ySXCB1rd8T9yfcxrNEwpGWnwd7THgOODkBCekKez6G+pfz73Lvs4cH1Lq9ezRXE3xbIPj4+\naNwY0NHhtpNFgQwAMy7NQNKnJHSr1w1DTYfK5iBygM5PfslrPsuqloVLBxc8cHqAbvW6IelTEqZc\nmIIWu1ogKDao2OMRiXIvkAHgwgXuw3GTJj6Y8d/qtf365X/f++/tR2BsIKprVcf8dvOLHmwJIa/n\npqKifBZc9epc73KrVl97l83NC7YPKpLlXPmy5XFowCEc6H8A2mW0cfrxaZhtNcOl55eEDq1EqVyZ\nW8VvzJicPxs6FIiM5ApqWfB65gWPBx7QUNWAew93ulNASox6FerhwogLODHkBPS19XHn9R203NUS\nk85NwoeMD0KHBwB49Yrrcbx/n7ub5OzMLTWfHx8zP2Lu1bkAgJW2K2lGIkLkyOfeZX//r73LOjoF\n2we1WyiQyKRIjDk1Br5R3GCYaS2mYYXNCpRVLStwZKSw0rLSYOpuilfJr2heVVKipWalYsnNJVgX\nsA5iqRiVNCphpe1KjDMfByWRsNdrkpOBQ4e4sQpz5wJq+ZyY4g/vP7DafzVa6reE3wQ/wX8PQsjP\npacDmpr5rzmpSFYwEqkEq/xW4X8+/4NYKoZpZVMcHHAQTao1ETo0UgizvWdjjf8amFczx7/2/0JF\nSUXokAiRqUfvHsHpvBNuvLoBAGht0BruPdwV7m/Y0/dP0ci9EcRSMQInBqKZXjOhQyKE5ENBak76\n2KtglJWU8We7PxFgF4AGFRvg4buHaL6rOdb4r4GUSalviWeyzOfd13exPmA9lERK2NFrR6kokOn8\n5Jci5tOksgmuj72OA/0PoKpmVfhH+8NyhyWmX5yOj5kfBYuroLmccWkGsqXZmGAxgQrkXCjiuSnP\nKJ/CoCJZQVnVsEKIQwgmW01GliQLs71nw3afLd6mvRU6NJIPEqkEDuccIGESTG0+ld5kSakiEokw\nsvFIPHF+gmktpgEANgZuRMMtDXE47LDc31k8//Q8Ljy7AO0y2nC1cRU6HEKIjFC7RQlw7uk5TDgz\nAe/S30FXXRfbem7DENMhNABMjm28vRHTL02HgbYBHjo9pAE/pFQLjQ+F03knBMQEAAA6GnaEWw83\nGFc2FjiynDLFmTDbaoZnH55hXZd1mNFqhtAhEUIKoCA1JxXJJcSb1DewO2uH88/OAwAsqlnA0coR\nI8xGQFNNU+DoyLeikqNg4maCtOw0nB12Fr0b9hY6JEIEJ2VS7Lm7B3OuzMH7jPdQVVLFzFYzscB6\ngVz9DVvltwpzrsyBUSUj3J98H6rKMlhZiBAiM9STXApV1aoKz+GemFFtBippVMLd+LtwOOcAvXV6\n+O3ib3ic8FjoEBUS331gjDE4X3BGWnYaBhoPLHUFMvXV8ask5VNJpAQ7Szs8cX4Ch6YOEEvFWOm3\nEibuJjgVfkrmF1Lyk8vXKa+x5OYSAMDGbhupQP6JknRuygPKpzCoSC5BRCIR+jTsg+gZ0djffz9a\nG7RGcmYyNgVugrGbMWz22eD4o+PIlmQLHWqpdTL8JDyfekK7jDY2dd8kdDiEyJ2KGhWxvdd2BNgF\nwKKaBaKSozDg6AD0PNQTER8iBI1t7tW5SM1KRZ+GfdClbhdBYyGEyB61W5RwofGh2Bq8FQfuH0B6\ndjoAoLpWdTg0dYC9pT30tPUEjrD0SP6UDGM3Y7xOfQ33Hu5wbOYodEiEyDWJVIJtwdsw/9p8JGcm\no4xyGcxtOxdz286Fuop6scZyO+Y2Wu1uBTVlNTxyeoS6FeoW6/EJIfygnmSSQ/KnZOy7tw9bg7ci\nPCEcAKAsUkZfo75wsnJCp9qdaKCfjDmdd8LW4K1opd8KtybcooUHCMmnN6lvMNt7Nvbf3w8AqFu+\nLjZ334zu9bsXy/GlTIoWu1ogOC4Y89rNw7JOy4rluIQQ/lFPcimWV9+SjroOpraYiodOD3F97HUM\nNhkMkUiEk+EnYbvfFsZuxth4eyOSPiUVb8Byjq8+MP9of2wN3goVJRXs6L2j1BbI1FfHr9KSz6pa\nVbGv/z7cGHcDppVNEZEYgR6HemDAkQGISo7i5Rg/y+U/of8gOC4YNcrVwJ9t/+TleCVdaTk3iwvl\nUxil8526FBOJROhg2AFHBx9F1PQoLO6wGHrl9PDk/RNMvzQdNdbWwMSzExHyOkToUEuMLEkWHDwd\nAAB/tPkDjao0EjgiQhSTdS1r3J10F6s7r4amqiZOPT4FYzdjrLy1ElmSLJkcM/lTMv68yhXGqzuv\nhpaalkyOQwiRP9RuQSCWiuH5xBNbg7fC+4X3l++30GsBRytHDDEdgrKqZQWMULG5+rpi/rX5qFeh\nHu5Pvk+5JIQHMR9jMPPSTBx7dAwAYFzJGG493NCxdkdejzPr8iysC1iHNgZt4Dvel9rSCFFwxdJu\ncezYMZiamkJZWRkhIXlfdbx48SKMjIxQv359rFy5srCHIzKkoqSC/sb9cXn0ZTxxfoIZLWdAV10X\ngbGBGHdmHPTX62O292zBR5Yromfvn2HxjcUAgG09t1GBTAhP9LX1cXTwUVwadQn1K9RHeEI4Ou3r\nhJEnR+J1ymtejvE44TE2BW6CCCJs6r6JCmRCSplCF8lmZmY4deoUrK2t89xGIpHA2dkZFy9exKNH\nj3D48GGEh4cX9pAkH4rat9SgYgOs67oOsTNjsbvPbjSt3hQfMj5gjf8a1NtcD90OdMPZJ2chkUr4\nCVjOFSWfjDFMPj8ZmZJMjGkyBjZ1bPgLTEFRXx2/KJ9Al7pdEOYYhiUdl0BdRR2Hwg7ByM0IG29v\nhFgqzvd+fswlYwy/XfwNYqkYEy0nwrK6Jc+Rl2x0bvKL8imMQhfJRkZGaNCgwU+3CQoKQr169WBo\naAhVVVUMGzYMZ86cKewhSTHSUNXABIsJCHYIRtDEIIwzHwd1FXVciriEvh59UWdTHbj6uuJN6huh\nQ5Vb++/vx7WX11CxbEWs7bJW6HAIKbHKqJTBX9Z/4ZHTI/Ru0BsfMz9i+qXpsNphBf9o/0Lt0/Op\nJy5HXIZOGR2azYKQUkpFljuPjY2FgYHBl8f6+voIDAzMddtx48bB0NAQAKCrqwtzc3N06NABwNdP\nUPT41487dOjA+/7TnqVhrM5YrJ25Fv+E/oN1h9chKjIK85Pnw8XHBe1YO/Rt2BdTh0yFSCSSq3wU\n9XFh85n8KRkz788EANhXtMeDoAdy8fsI/VgW52dpfkz5/P5x7fK1MbP6TLQQt8DO9ztx7809tPlf\nG3Sv1x17p+9FZc3K+dpfliQLMx7OAACM1hmNh/8+lIvfjx7TY3pc8MehoaFISuJm7oqMjERB/HTg\nXufOnREfH5/j+66urujdm1tOt2PHjli7di0sLXPeijpx4gQuXryInTt3AgAOHDiAwMBAbN68+fsg\naOCeQpEyKa68uAL3f93h+dQTUiYFADSq0giOVo4Y1XgUtMtoCxylsMaeHot99/bBprYNvEd7Uy8j\nIcUsPTsdrr6uWOW3CtnSbJRXL4/lNsth39T+l1MwLvddjnnX5sGksglCJ4XS8tOElCC8Ddzz9vZG\nWFhYjq/PBfKv6OnpITo6+svj6Oho6Ovr5+u5pHA+f4qSJSWRErrU7YLTw04j8rdI/GX9F6pqVsWD\ntw8w5cIU6K3Tg9N5J4S9CZN5LLJWmHxeeXEF++7tg7qKOrb12kYF8jeK4/wsTSifedNQ1cDSTksR\n5hgG2zq2SPyUiMnnJ6PV7la4E3cnx/afcxn7MRbLfLn2io3dNlKBXEh0bvKL8imMn3+czqe8KnIr\nKys8e/YMkZGRyMrKwpEjR9CnTx8+DknkhIGOAZZ0XIKoGVE4MugI2tdqj9SsVGwN3orG2xqj3Z52\nOBx2GJniTKFDLRYZ2RmYfG4yAGCB9QLUq1BP4IgIKd0aVmqIy6Mu48igI6hRrgaCYoPQbGczTLkw\nBYkZiTm2n3NlDtKy09DfqD9s69gKEDEhRF4Uep7kU6dOYdq0aUhISICOjg4sLCzg5eWFuLg42Nvb\n4/z58wAALy8vTJ8+HRKJBHZ2dvjzz5yrFVG7Rcny8O1DbA3ein339iElKwUAUEWzCiZaToSDpQNq\n6dYSOELZmX9tPlx9XdGoSiPccbgDNWU1oUMihPwnJTMFi24swobbGyBhElTWqIw1XdZgdOPREIlE\n8IvyQ9s9bVFGuQzCp4SjdvnaQodMCOFZQWpOWkyEyExKZgoOhR2Ce7A77r+5D4Br1ehZvyecmjmh\nS90uJWp55rA3YbDcYQmJVAK/CX5oZdBK6JAIIbkIexOGKRemwDfKFwDQrmY7bO6+GRPOTkDI6xD8\nZf0XlnRcInCUhBBZKJbFRIh8kqe+pXJlymGS1SSETgrFrfG3MNJsJFSUVOD51BPdD3ZH/c31sdpv\nNRLSE4QONU/5zaeUSeFwzgFiqRiOzRypQM6DPJ2fJQHls3DMqprhxrgb2NtvL6poVoFvlC/M/zRH\nyOsQ6GvrY26buUKHqPDo3OQX5VMYVCQTmROJRGhTsw0ODDiA6BnRWG6zHLV0auFF4gv8ceUP6K/T\nx9jTY3E75rbC3lHYFrwNt2Nuo7pWdbh2chU6HELIL4hEIoxpMgZPnJ9gSrMpX76/uvNqaKppChgZ\nIUReULsFEYREKsHF5xfhHuwOr2deYOD+/y2qWcCpmROGNxquMG9UcSlxMHYzxsfMjzg++DgGmgwU\nOiRCSAHdf3Mf0cnR6FG/B81IQ0gJRj3JRKG8SHyB7Xe2Y3fIbrzPeA8A0Cmjg3Hm4zDZajKMKhkJ\nHOHPDTo6CCfCT6BPwz44PfQ0vcESQgghcop6kksxRexbqlO+DlbarkTMzBjs778frfRbITkzGRsD\nN8LYzRg2+2xw4tEJZEuyiz22X+Xz7JOzOBF+AlpqWtjSfQsVyL+giOenPKN88odyyS/KJ78on8Kg\nIpnIDXUVdYxqPAr+dv64O+kuHJo6QENVA9deXsOgY4NguNEQLj4uiP0YK3SoALjZO6Zc4HoZl3Va\nBgMdg188gxBCCCGKgtotiFxL/pSMfff2wT3YHY8THgMAlEXK6GfUD07NnNDRsKNgV2+nX5yOjYEb\n0axGMwTYBUBZSVmQOAghhBCSP9STTEocxhhuvLoB93/dcerxKYilYgBAw4oN4WjliLHmY6Grrlts\n8fwb+y9a7GoBJZESgh2CYV7NvNiOTQghhJDCoZ7kUqyk9i2JRCJ0MOyAo4OP4tX0V1jUYRH0yunh\nyfsnmH5pOmqsrQF7T3uEvA7h9bi55TNbkg17T3swMMxsNZMK5AIoqeenUCif/KFc8ovyyS/KpzCo\nSCYKp0a5Gvhf+/8hcnokTg45Cds6tsgQZ2BXyC403dEULXe1xL57+/BJ/Ekmx99wewPuvbkHQ11D\nLGy/UCbHIIQQQoiwqN2ClAhPEp5g251t+Cf0HyR9SgIAVCxbEeMtxmNy08moW6EuL8d5mfgSpu6m\nyBBnwGukF7rV68bLfgkhhBAie9STTEqt9Ox0eDzwgPu/7rjz+s6X73er1w1OVk7oUb9HoQfYMcbQ\n41APXHx+EcMbDcehgYf4CpsQQgghxYB6kkux0t63pKGqgQkWExDsEIygiUEYZz4O6irquPj8Ivp4\n9EGdTXXg6uuKN6lv8rW/b/Pp8cADF59fhK66LtZ3XS+j36BkK+3nJ98on/yhXPKL8skvyqcwqEgm\nJVYzvWbY03cPYmbEYE3nNahXoR6ikqMw/9p8GKw3wIgTI+D7yjdfnyg/ZHzA9EvTAQBrOq9BVa2q\nsg6fEEIIIQKidgtSakiZFFdeXIH7v+7wfOoJKZMCABpVaQQnKyeMajwK5cqUy/W5E89OxO67u2Fd\nyxo+Y31oZT1CCCFEAVFPMiG/EJUchZ0hO7Hzzk68SeNaL7TUtDC68Wg4WjnCrKrZl21vvrqJ9v+0\nh5qyGu5NvgejSkZChU0IIYSQIqCe5FKM+pbyp6ZOTSzpuARRM6LgMdAD7Wu1R2pWKrYGb0XjbY1h\nvccaHg88cMH7Ahw8HQAA89rOowK5iOj85Bflkz+US35RPvlF+RSGitABECIkNWU1DG00FEMbDcXD\ntw+xNXgr9t3bB98oX/hG+UItWg1ZBlkwqmSEuW3nCh0uIYQQQooJtVsQ8oOUzBQcDDsI93/dEfY2\nDABwY9wNWNeyFjgyQgghhBQF9SQTwgPGGAJjAyGRStCmZhuhwyGEEEJIEVFPcilGfUv8EYlE+PT8\nExXIPKLzk1+UT/5QLvlF+eQX5VMYVCQTQgghhBDyA2q3IIQQQgghpQK1WxBCCCGEEFIEVCSXMNS3\nxC/KJ78on/yifPKHcskvyie/KJ/CoCKZEEIIIYSQH1BPMiGEEEIIKRWoJ5kQQgghhJAioCK5hKG+\nJX5RPvlF+eQX5ZM/lEt+UT75RfkUBhXJhBBCCCGE/IB6kgkhhBBCSKlAPcmEEEIIIYQUARXJJQz1\nLfGL8skvyie/KJ/8oVzyi/LJL8qnMKhIJoQQQggh5AfUk0wIIYQQQkoF6kkmhBBCCCGkCKhILmGo\nb4lflE9+UT75RfnkD+WSX5RPflE+hUFFcgkTGhoqdAglCuWTX5RPflE++UO55Bflk1+UT2EUukg+\nduwYTE1NoaysjJCQkDy3MzQ0ROPGjWFhYYHmzZsX9nAkn5KSkoQOoUShfPKL8skvyid/KJf8onzy\ni/IpDJXCPtHMzAynTp3CpEmTfrqdSCSCj48PKlSoUNhDEUIIIYQQUqwKXSQbGRnle1uauaL4REZG\nCh1CiUL55Bflk1+UT/5QLvlF+eQX5VMYRZ4CrmPHjli7di0sLS1z/XmdOnWgo6MDZWWCoL0LAAAG\nRUlEQVRlTJo0Cfb29jmDEImKEgIhhBBCCCH5kt/S96dXkjt37oz4+Pgc33d1dUXv3r3zdQA/Pz9U\nr14d7969Q+fOnWFkZIR27doVKlhCCCGEEEKKw0+LZG9v7yIfoHr16gCAypUro3///ggKCspRJBNC\nCCGEECJPeJkCLq8rwenp6UhJSQEApKWl4fLlyzAzM+PjkIQQQgghhMhMoYvkU6dOwcDAALdv30bP\nnj3RvXt3AEBcXBx69uwJAIiPj0e7du1gbm6OFi1aoFevXujSpQs/kRNCCCGEECIjRR64V1QXL17E\n9OnTIZFIMHHiRMyZM0fIcBTahAkTcP78eVSpUgVhYWFCh6PQoqOjMWbMGLx9+xYikQgODg6YNm2a\n0GEprE+fPqF9+/bIzMxEVlYW+vbti+XLlwsdlkKTSCSwsrKCvr4+PD09hQ5HoRkaGkJbWxvKyspQ\nVVVFUFCQ0CEptKSkJEycOBEPHz6ESCTC33//jZYtWwodlkJ68uQJhg0b9uXxixcvsGTJEno/KqTl\ny5fjwIEDUFJSgpmZGfbs2YMyZcrkub2gRbJEIkHDhg1x5coV6OnpoVmzZjh8+DCMjY2FCkmh+fr6\nQktLC2PGjKEiuYji4+MRHx8Pc3NzpKamomnTpjh9+jSdm0WQnp4ODQ0NiMVitG3bFmvWrEHbtm2F\nDkthrVu3Dnfu3EFKSgrOnj0rdDgKrXbt2rhz5w7N58+TsWPHon379pgwYQLEYjHS0tKgo6MjdFgK\nTyqVQk9PD0FBQTAwMBA6HIUTGRmJTp06ITw8HGXKlMHQoUPRo0cPjB07Ns/nCLosdVBQEOrVqwdD\nQ0Ooqqpi2LBhOHPmjJAhKbR27dqhfPnyQodRIlSrVg3m5uYAAC0tLRgbGyMuLk7gqBSbhoYGACAr\nKwsSiYQKkiKIiYnBhQsXMHHiRJodiCeUR34kJyfD19cXEyZMAACoqKhQgcyTK1euoG7dulQgF5K2\ntjZUVVWRnp4OsViM9PR06Onp/fQ5ghbJsbGx3/1n6+vrIzY2VsCICMkpMjISd+/eRYsWLYQORaFJ\npVKYm5ujatWq6NixI0xMTIQOSWHNmDEDq1evhpKSoH/CSwyRSARbW1tYWVlh586dQoej0F6+fInK\nlStj/PjxsLS0hL29PdLT04UOq0Tw8PDAiBEjhA5DYVWoUAGzZs1CzZo1UaNGDejq6sLW1vanzxH0\nLywtIkLkXWpqKgYNGoSNGzdCS0tL6HAUmpKSEkJDQxETE4ObN2/Cx8dH6JAU0rlz51ClShVYWFjQ\n1U+e+Pn54e7du/Dy8oKbmxt8fX2FDklhicVihISEwMnJCSEhIdDU1MSKFSuEDkvhZWVlwdPTE4MH\nDxY6FIUVERGBDRs2IDIyEnFxcUhNTcXBgwd/+hxBi2Q9PT1ER0d/eRwdHQ19fX0BIyLkq+zsbAwc\nOBCjRo1Cv379hA6nxNDR0UHPnj0RHBwsdCgKyd/fH2fPnkXt2rUxfPhwXLt2DWPGjBE6LIWW23z+\npHD09fWhr6+PZs2aAQAGDRqEkJAQgaNSfF5eXmjatCkqV64sdCgKKzg4GK1bt0bFihWhoqKCAQMG\nwN/f/6fPEbRItrKywrNnzxAZGYmsrCwcOXIEffr0ETIkQgBw/Yl2dnYwMTHB9OnThQ5H4SUkJCAp\nKQkAkJGRAW9vb1hYWAgclWJydXVFdHQ0Xr58CQ8PD3Tq1An79u0TOiyFRfP586tatWowMDDA06dP\nAXB9tKampgJHpfgOHz6M4cOHCx2GQjMyMsLt27eRkZEBxhiuXLnyy7a/n664J2sqKirYsmULunbt\nColEAjs7O5o9oAiGDx+OGzdu4P379zAwMMDixYsxfvx4ocNSSH5+fjhw4AAaN278pZhbvnw5unXr\nJnBkiun169cYO3YspFIppFIpRo8eDRsbG6HDKhGoba1o3rx5g/79+wPgWgVGjhxJ8/kX0ebNmzFy\n5EhkZWWhbt262LNnj9AhKbS0tDRcuXKF+uWLqEmTJhgzZgysrKygpKQES0tLODg4/PQ5gs+TTAgh\nhBBCiLyhodGEEEIIIYT8gIpkQgghhBBCfkBFMiGEEEIIIT+gIpkQQgghhJAfUJFMCCGEEELID6hI\nJoQQQggh5Af/B8lCcNd3RwTuAAAAAElFTkSuQmCC\n",
       "text": "<matplotlib.figure.Figure at 0xe888e90>"
      }
     ],
     "prompt_number": 284
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "resid",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 285,
       "text": "array([[-0.60184156],\n       [-0.42663825],\n       [-0.43852397],\n       [-0.18071406],\n       [-0.07025626],\n       [-0.25620732],\n       [ 0.11442828],\n       [ 1.16502185],\n       [ 0.02660395]])"
      }
     ],
     "prompt_number": 285
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Laten we eens vergelijken met standaard OLS-fitting library"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from statsmodels.formula.api import ols\ndata = pandas.DataFrame({'Y': Y.squeeze(), 'X1': X[:,0], 'X2': X[:, 1]})\nmodel = ols('Y ~ X1 + X2- 1', data)\nprint'variance:\\n%s' % str(1/(model.fit().t() / model.fit().params) ** 2)\nmodel.fit().summary()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "variance:\nX1    0.093581\nX2    0.093581\ndtype: float64\n"
      },
      {
       "html": "<table class=\"simpletable\">\n<caption>OLS Regression Results</caption>\n<tr>\n  <th>Dep. Variable:</th>            <td>Y</td>        <th>  R-squared:         </th> <td>   0.435</td>\n</tr>\n<tr>\n  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.355</td>\n</tr>\n<tr>\n  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   5.394</td>\n</tr>\n<tr>\n  <th>Date:</th>             <td>Wed, 14 Aug 2013</td> <th>  Prob (F-statistic):</th>  <td>0.0532</td> \n</tr>\n<tr>\n  <th>Time:</th>                 <td>14:43:59</td>     <th>  Log-Likelihood:    </th> <td> -6.4531</td>\n</tr>\n<tr>\n  <th>No. Observations:</th>      <td>     9</td>      <th>  AIC:               </th> <td>   16.91</td>\n</tr>\n<tr>\n  <th>Df Residuals:</th>          <td>     7</td>      <th>  BIC:               </th> <td>   17.30</td>\n</tr>\n<tr>\n  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n</tr>\n</table>\n<table class=\"simpletable\">\n<tr>\n   <td></td>     <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n</tr>\n<tr>\n  <th>X1</th> <td>    0.7799</td> <td>    0.306</td> <td>    2.549</td> <td> 0.038</td> <td>    0.057     1.503</td>\n</tr>\n<tr>\n  <th>X2</th> <td>   -0.8400</td> <td>    0.306</td> <td>   -2.746</td> <td> 0.029</td> <td>   -1.563    -0.117</td>\n</tr>\n</table>\n<table class=\"simpletable\">\n<tr>\n  <th>Omnibus:</th>       <td>12.044</td> <th>  Durbin-Watson:     </th> <td>   1.213</td>\n</tr>\n<tr>\n  <th>Prob(Omnibus):</th> <td> 0.002</td> <th>  Jarque-Bera (JB):  </th> <td>   4.871</td>\n</tr>\n<tr>\n  <th>Skew:</th>          <td> 1.572</td> <th>  Prob(JB):          </th> <td>  0.0876</td>\n</tr>\n<tr>\n  <th>Kurtosis:</th>      <td> 4.762</td> <th>  Cond. No.          </th> <td>    1.73</td>\n</tr>\n</table>",
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 286,
       "text": "<class 'statsmodels.iolib.summary.Summary'>\n\"\"\"\n                            OLS Regression Results                            \n==============================================================================\nDep. Variable:                      Y   R-squared:                       0.435\nModel:                            OLS   Adj. R-squared:                  0.355\nMethod:                 Least Squares   F-statistic:                     5.394\nDate:                Wed, 14 Aug 2013   Prob (F-statistic):             0.0532\nTime:                        14:43:59   Log-Likelihood:                -6.4531\nNo. Observations:                   9   AIC:                             16.91\nDf Residuals:                       7   BIC:                             17.30\nDf Model:                           1                                         \n==============================================================================\n                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n------------------------------------------------------------------------------\nX1             0.7799      0.306      2.549      0.038         0.057     1.503\nX2            -0.8400      0.306     -2.746      0.029        -1.563    -0.117\n==============================================================================\nOmnibus:                       12.044   Durbin-Watson:                   1.213\nProb(Omnibus):                  0.002   Jarque-Bera (JB):                4.871\nSkew:                           1.572   Prob(JB):                       0.0876\nKurtosis:                       4.762   Cond. No.                         1.73\n==============================================================================\n\"\"\""
      }
     ],
     "prompt_number": 286
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "OK. Dit lijkt allemaal te kloppen."
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Nu de sandwich.\n\n$\\hat{V}_W  = (X^TX)^{-1}X^T W X(X^TX)^{-1}$\n\nen\n\n$ W = \\hat{r}\\hat{r}^T$"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Eerst maar eens W, dit is een covariance matrix\nW = resid.dot(resid.T)\nplt.imshow(W, interpolation='nearest')\nplt.colorbar()\nplt.title('$W$', fontsize=24)\n\nprint 'diag(W):\\n%s' % np.diag(W)\n\nprint '\\nAlles op de diagonaal van W > 0? %s' % (np.all(np.diag(W) > 0))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "diag(W):\n[  3.62213265e-01   1.82020194e-01   1.92303276e-01   3.26575722e-02\n   4.93594248e-03   6.56421884e-02   1.30938317e-02   1.35727591e+00\n   7.07770087e-04]\n\nAlles op de diagonaal van W > 0? True\n"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAERCAYAAAAXPEBjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9QVNf9N/D3ld38MEQQo4C79IsKzgIKoiCPVSqmIokd\nqT9IBoOPjkFC03ES86PRzjxJ0SdJZUxqjDQdY9Tqg1VaOwnU0I0yEbX4xY2RxFjNV2JhsoAwWsWA\nxADrPn9QVxFY9nL33oVz36+ZO+PuPWfPOTP68fy690hOp9MJIqIhbpivK0BE5A0MZkQkBAYzIhIC\ngxkRCYHBjIiEwGBGREJgMCMiITCYEZEQGMwEs3XrVlgsFgwbNsx1TZgwAZs2bXKlWbt2LSZMmOC6\nf//992PatGnYuXNnt9+aOXMm7rvvPgwbNgxGoxFJSUm4cuWK1k0i8ojEJwDENG3aNFRVVSEuLg5V\nVVW9pjGZTGhsbMSBAwewaNGiXtPs3bsXeXl5sFqtmDBhgppVJlKEPTNBjRkzBgDw8MMP95nGaDTC\n6XTCYDD0maaqqgr79+9nIKNBj8FMUEFBQW7v19XVoampCQBw9erVXtPU1NSgtbUV06ZN83r9iLyN\nwUxQ/QWzdevWISQkBABw7dq1XtO88cYbeO2117xeNyI1MJgJauTIkX3e+8c//gGDwYAZM2YA6L1n\ndurUKYSGhmLs2LGq1ZHImxjMBNVXz+zWrVtYu3YtXn/9dVfA6y2Ybdq0CWvXrlW1jkTexGAmqL6C\n2fbt25GcnAyz2exKc+8w8+DBg/jJT34Cf39/1etJ5C19L2PRkNZbMLt27Ro2b94Mm80GABg1ahSA\n7j0zh8OBnTt34i9/+Ys2FSXyEvbMBNVbMPvNb36DX/ziFxgxYkS3NHcHs127dmHZsmXw8/PTpqJE\nXsKemaDuDWZnz56F1WrFuXPneqS5Pcxsa2vDJ598wl4ZDUnsmQnq3tXM559/Hhs2bOi2QfbeYebv\nfvc7vPjii9pVksiLGMwEdXfP7K9//StaWlqQmZnZa5rr16+jqakJFy9edG3XIBpqOMwUlNFohL+/\nP77//nv86le/6vEQOXAnmDkcDrzyyit49dVXta4mkdewZyawkSNH4vPPP0d0dDRSUlJ63L+79xYQ\nEICIiAgNa0fkXXxrhsDi4+Nx9uxZfPnll4iOju41TWBgICRJwjfffOOaQyMaijjMFNioUaOwYsWK\nPgPZ7TQ5OTkMZDTksWcmsA8//BDJycl45JFH+kxTWFiIJ554Avfff7+GNSPyPgYzIhICFwCISAgM\nZkQkBEULAJIkeaseRDQASmaJgoKC+nwxZ29GjhzZ51uJBwNFc2aSJCHfuVpWnsN5NqTmTZeV5xz6\nXo3zFgMcHqc9nVeKqXnzZZdxBHNk55Gr4Xqox2k7N+bDsE7+O8tiA76SnUcu2yezPU9cmAcsy5Nf\nyGMDyCPXWQ/LeC8P+KWHaW+bJCkKZpIkQU6JeVAWPNXGrRlEOiZSABCpLUQkk9HXFfAizYPZ+BST\n1kV6XWhKpK+r4BXDZs30dRW8IzbF1zVQLjHFJ8WK1JvRvC0TGMwGjWGzZvm6Ct7BYDZgD/qkVHWI\nFJiJSCYOM4lICCIFAG6aJdIxo4zrXk8//TSCg4MxefLkXn977969iIuLQ2xsLGbOnIkzZ8647oWH\nhyM2Nhbx8fGYPl3eVq2+9BvMrFYrLBYLIiMjkZ+f75VCiWhwMMi47rVy5UpYrdY+f3v8+PE4duwY\nzpw5g1dffRXPPPOM654kSSgvL0dVVZXrtDCl3AYzh8OB1atXuw7C2LdvH86fP++VgonI95T0zJKT\nk3ucNXG3GTNmICAgAACQlJSEurq6bve9vQHXbTCz2WyIiIhAeHg4jEYjMjMzUVxc7NUKEJHvuAte\nXwM4cNelxI4dOzB//p2nZiRJwty5c5GQkIDt27cr/PUubuf/6uvrERYW5vpsNptx8uTJbmkO593p\nIo5PMQmx9YJoUPqsvOvyIndbM6b/57pt3wDLOHLkCHbu3ImKigrXdxUVFQgNDcXly5eRmpoKi8WC\n5OTkAZbQxW0w8+RBcrnPWRLRACWmdN+P9t56xT+p9mrmmTNnkJOTA6vV2m1IGhra9Qzx6NGjsWjR\nIthsNsXBzO0w02QywW63uz7b7XaYzWZFBRLR4KFkzqw/3377LRYvXozCwsJuh+W0tbWhpaUFAHDj\nxg0cOnSozxVROdwG5oSEBFRXV6O2thZjx45FUVER9u0baGeTiAYbJT2zpUuX4ujRo7hy5QrCwsKw\nfv16dHR0AAByc3OxYcMGXLt2Dc8++yyAruMPbTYbGhsbsXjxYgBAZ2cnsrKyMG/ePKVNcd8Wg8GA\ngoICpKWlweFwIDs7G1FRUYoLJaLBQckTAP11bD744AN88MEHPb4fP348vvjiCwUl967fwPz444/j\n8ccf93rBROR7Ij0BIFJbiEgmPptJRELgWzOISAjsmRGREEQKACK1hYhkMsqJAJ2qVcMrGMyIdMzA\nYEZEIjD6+boG3qM4mGlxpmUtxqlehp8G/+204GHVy2i/eb/qZVwJeET1MuCvfhHAi1oUMqjJ6pkN\ncgI1hYjkMqr/f59mGMyI9EygCCBQU4hINoEigEBNISLZBIoAAjWFiGTjaiYRCUGgCCBQU4hINq5m\nEpEQBIoAAjWFiGQTKAL0e6I5EQnMT8Z1j6effhrBwcFuDyN57rnnEBkZibi4OFRVVbm+t1qtsFgs\niIyMRH5+vleawmBGpGcGGdc9Vq5cCavV2udPl5aW4ptvvkF1dTXef/9918EmDocDq1evhtVqxblz\n57Bv3z6cP39ecVMYzIj0TEEwS05O7nYW5r1KSkqwYsUKAEBSUhKam5vR2NgIm82GiIgIhIeHw2g0\nIjMzE8XFxV5pChHplZsIUH616xqo+vp6hIWFuT6bzWbU19ejoaGhx/cnT54ceEH/wWBGpGdutmak\nhHZdt62/KP/nnU6n/EwDxGBGpGcqRgCTyQS73e76XFdXB7PZjI6Ojm7f2+12mM1mxeVxzoxIzxSs\nZvYnPT0de/bsAQBUVlYiMDAQwcHBSEhIQHV1NWpra9He3o6ioiKkp6crbgp7ZkR6piACLF26FEeP\nHsWVK1cQFhaG9evXo6OjAwCQm5uL+fPno7S0FBEREXjooYewa9euriINBhQUFCAtLQ0OhwPZ2dmI\niopS3BTJqWBQK0kSVjjfU1yJ/ojyptmvEKt6Gf9uGqV6GeHBtaqX8a+KGNXLwKzv1C/j7Aj1fnuS\npGhOSpIkOP+3jPT/T9s5MLnYMyPSM741g4iEIFAEEKgpRCTbA76ugPcwmBHpGYeZRCQEgSKAQE0h\nItkEigCKm2KAwxv1cEuLbROitMPPIEY7tCgC+F6DMlTcmuENHGYSkRAEigACNYWIZBMoAgjUFCKS\njQeaEJEQBIoAAjWFiGQTKAII1BQiko2rmUQkBIEigEBNISLZBIoAbt80a7fbMWfOHMTExGDSpEl4\n9913taoXEWlBxTfNas1tXDYajdi8eTOmTJmC1tZWTJs2DampqV55KyQRDQICvTXDbc8sJCQEU6ZM\nAQD4+/sjKioKDQ0NmlSMiDSg4NzMwcbjKtbW1qKqqgpJSUndvj+dV+r6c2hKJEJTIr1XOyK647Py\nrsubhsDw0VMeBbPW1lZkZGRgy5Yt8Pf373Zvat58VSpGRPdITOm6bntvvfLfVNjjslqtWLNmDRwO\nB1atWoW1a9d2u//WW29h7969AIDOzk6cP38eV65cQWBgIMLDwzFixAj4+fnBaDTCZrMpqku/Teno\n6MCSJUuwbNkyLFy4UFFhRDTIKAhmDocDq1evRllZGUwmExITE5Gent5tTv3ll1/Gyy+/DAA4ePAg\n3nnnHQQGBgLoOlClvLwcQUFBippwm9umOJ1OZGdnIzo6GmvWrPFKgUQ0iLgZZpZ/CZSf6fu+zWZD\nREQEwsPDAQCZmZkoLi7uc4HwT3/6E5YuXdrtO2+e9uR2AaCiogKFhYU4cuQI4uPjER8fD6vV6rXC\nicjHHuj7SkkC8nLuXPeqr69HWFiY67PZbEZ9fX2vxbS1teGTTz7BkiVLXN9JkoS5c+ciISEB27dv\nV9wUtz2zWbNm4datW4oLIaJBSsECgCRJHqf929/+hlmzZrmGmEBXZyk0NBSXL19GamoqLBYLkpOT\nB1wftz0zIhKcgq0ZJpMJdrvd9dlut8NsNvdazP79+3sMMUNDQwEAo0ePxqJFixQvADCYEemZgmCW\nkJCA6upq1NbWor29HUVFRUhPT++R7vr16zh27Bh+/vOfu75ra2tDS0sLAODGjRs4dOgQJk+erLgp\nRKRXCiKAwWBAQUEB0tLS4HA4kJ2djaioKGzbtg0AkJubCwD46KOPkJaWhgcffNCVt6mpCYsWLQLQ\ntWUjKysL8+bNG3hlAEhOBcsJkiQh27lVUQU8UYNw1cvQ4kCTM1D2P48n/v3vR1QvI3xUjeplVB+N\nU70MpDSpX8bZYPV+e5KkaDVQkiQ4v5KRfrJ3Vx+9jT0zIj0TKAII1BQiko1nABCREASKAIqbcgRz\nvFEPt1rwsOplaHGwrRbzWR2tD/afSCG7Iaz/REoF9p9EsZdVnM+6bbAHi8FePxkEagoRySZQBBCo\nKUQkl1NvrwAiIjE5BIoAAjWFiORiMCMiIfxw/30yUrerVg9vYDAj0jGHnziTZgxmRDrmEOgQAAYz\nIh3rZDAjIhE4BAoB4rSEiGTjMJOIhMBgRkRC+AFytmYMbgxmRDrGOTMiEoJIw0weaEKkYw74eXz1\nxmq1wmKxIDIyEvn5+T3ul5eXIyAgwHXu7uuvv+5xXrnYMyPSMSX7zBwOB1avXo2ysjKYTCYkJiYi\nPT29x4nms2fPRklJyYDyysGeGZGOOWDw+LqXzWZDREQEwsPDYTQakZmZieLi4h7pejsExdO8crBn\nRqRj7ubMqspbUFXe0uf9+vp6hIXdeeuw2WzGyZMnu6WRJAknTpxAXFwcTCYT3nrrLURHR3uUVy4G\nMyIda3ezNSMmZRRiUka5Pu9cf6nbfUmS+v39qVOnwm63Y/jw4fj73/+OhQsX4sKFCwOvsBscZhLp\nWCf8PL7uZTKZYLfbXZ/tdjvMZnO3NA8//DCGDx8OAHj88cfR0dGBq1evwmw295tXLgYzIh1TMmeW\nkJCA6upq1NbWor29HUVFRUhPT++WpqmpyTVnZrPZ4HQ6ERQU5FFeuTjMJNIxJfvMDAYDCgoKkJaW\nBofDgezsbERFRWHbtm0AgNzcXBw4cAB/+MMfYDAYMHz4cOzfv99tXiUkp4Lz1iVJwnjnWUUV8ASP\nmvOcFkfNPRDY96Swt9ysDVK9DBSqXwRWqfjbFqnXlUJPSZKEYuc8j9P/XDqkqDy1sWdGpGN8n9ld\nGq6HeqMebrXfVP8MeT+DQ/UytOg1oVP9/59utg5XvQw8oH4RSNGgjEGuHer/29IKe2ZEOibSs5kM\nZkQ6xmEmEQmBrwAiIiFwmElEQmAwIyIhMJgRkRB+EGhrhkfPZjocDsTHx2PBggVq14eINKT0TbOD\niUc9sy1btiA6OhotLeo/xkJE2hkKQcpT/fbM6urqUFpailWrVg3q57KISD4lrwAabPrtmb3wwgvY\ntGkTvvvuu17vd268cxDBsFkzMWzWLO/VjojusJV3XV6km31mBw8exJgxYxAfH4/y8vLef2DdWjXq\nRUT3mp7Sdd1WsF7xT4o0zHQbzE6cOIGSkhKUlpbi5s2b+O6777B8+XLs2bNHq/oRkYpECmZu58ze\nfPNN2O121NTUYP/+/Xj00UcZyIgE8gPu8/ga7GQNmD05wICIhg7dzJndbfbs2Zg9e7aadSEijelm\nmElEYlO6adZqtcJisSAyMhL5+fk97u/duxdxcXGIjY3FzJkzcebMGde98PBwxMbGIj4+HtOnT1fc\nFnH6mEQkm5L9Yw6HA6tXr0ZZWRlMJhMSExORnp7e7WCS8ePH49ixYwgICIDVasUzzzyDyspKAF3T\nVuXl5QgK8s55DwxmRDqmZM7MZrMhIiIC4eHhAIDMzEwUFxd3C2YzZsxw/TkpKQl1dXXdfsObG/EZ\nzIh0zN2cWWP5BTSW9336eH19PcLCwlyfzWYzTp482Wf6HTt2YP78+a7PkiRh7ty58PPzQ25uLnJy\ncmTWvjsGMyIda3ez5SIoZRKCUia5Pn+x/uNu9+Xsbjhy5Ah27tyJiooK13cVFRUIDQ3F5cuXkZqa\nCovFguTkZBm1744LAEQ6puTZTJPJBLvd7vpst9thNpt7pDtz5gxycnJQUlKCkSNHur4PDe062W30\n6NFYtGgRbDaborYwmBHpmAMGj697JSQkoLq6GrW1tWhvb0dRURHS09O7pfn222+xePFiFBYWIiIi\nwvV9W1ub6y08N27cwKFDhzB58mRFbVE8zIwN+ErpT/TrSoD6J4FrcaK53RDWfyKFtDjTMijk36qX\ncfULk+ploED9IvCOBmUooGSfmcFgQEFBAdLS0uBwOJCdnY2oqChs27YNAJCbm4sNGzbg2rVrePbZ\nZwEARqMRNpsNjY2NWLx4MQCgs7MTWVlZmDfP89PVeyM5FSwnSJKE6c5yRRXwxBUIEsyuM5h5SpNg\n9n/UL0LVYGaRFK0GSpKEBc4/e5z+b9KTg/o1YFwAINKxofCeMk8xmBHpmC6fzSQi8bjbmjHUMJgR\n6RiHmUQkBA4ziUgIIr0CiMGMSMcYzIhICAxmRCSEH3C/r6vgNQxmRDrGnhkRCYHBjIiEwH1mRCQE\n7jMjIiFwmElEQmAwIyIh/NDOB82JSACOTnFCgDgtISLZHJ3iDDN5oAmRjjk6/Ty+emO1WmGxWBAZ\nGYn8/Pxe0zz33HOIjIxEXFwcqqqqZOWVg8GMSMc6O/w8vu7lcDiwevVqWK1WnDt3Dvv27cP58+e7\npSktLcU333yD6upqvP/++66DTTzJKxeDGZGO3XIYPL7uZbPZEBERgfDwcBiNRmRmZqK4uLhbmpKS\nEqxYsQIAkJSUhObmZjQ2NnqUVy7OmRHpmbs5s/8+ClQe7fN2fX09wsLunDhmNptx8uTJftPU19ej\noaGh37xyMZgR6dlNNyEg/qdd123v/N9utyVJ8qgIrY6nUxzMbJ/M9kY93PNXvwgNjs0EAjUo4wH1\ni9DkTMtHOlQvYtLfv1C9jLMXE1QvQxEFf+9NJhPsdrvrs91uh9lsdpumrq4OZrMZHR0d/eaVi3Nm\nRHrWKeO6R0JCAqqrq1FbW4v29nYUFRUhPT29W5r09HTs2bMHAFBZWYnAwEAEBwd7lFcuDjOJ9ExB\nz8xgMKCgoABpaWlwOBzIzs5GVFQUtm3bBgDIzc3F/PnzUVpaioiICDz00EPYtWuX27xKSE4FA1pJ\nkgCrBuNhDjM9p8EwE60alKHFMPO/hvgwM2KYovkoSZKAShn5/5ek2fzXQLBnRqRnDl9XwHsYzIj0\nTIsRiUYYzIj07KavK+A9DGZEesaeGREJQaBg1u8+s+bmZmRkZCAqKgrR0dGorKzUol5EpAUF+8wG\nm357Zs8//zzmz5+PAwcOoLOzEzdu3NCiXkSkBfV3wGjGbTC7fv06jh8/jt27d3clNhgQEBCgScWI\nSAN62ZpRU1OD0aNHY+XKlfjyyy8xbdo0bNmyBcOHD7+TqDDvzp9jU7ouIvK+k+XAyb7fYjEgQ2D4\n6Cm3TwCcOnUKM2bMwIkTJ5CYmIg1a9ZgxIgR2LBhQ1dmPgEgD58A8ByfAOifN54A2CYjf+7gfgLA\n7QKA2WyG2WxGYmIiACAjIwOnT5/WpGJEpAGBFgDcBrOQkBCEhYXhwoULAICysjLExMRoUjEi0oBA\nwazf1cytW7ciKysL7e3tmDBhguupdyISwBAIUp7qN5jFxcXhs88+06IuRKQ1vWzNICLB6WVrBhEJ\njg+aE5EQ9DRnRkQCE2jOjAeaEOmZQ8Ylw9WrV5GamoqJEydi3rx5aG5u7pHGbrdjzpw5iImJwaRJ\nk/Duu++67uXl5cFsNiM+Ph7x8fGwWq39lslgRqRnKu0z27hxI1JTU3HhwgX89Kc/xcaNG3ukMRqN\n2Lx5M/75z3+isrISv//97/H1118D6Ho64cUXX0RVVRWqqqrw2GOP9VsmgxmRnqkUzEpKSrBixQoA\nwIoVK/DRRx/1SBMSEoIpU6YAAPz9/REVFYX6+nrXfbmPTimfM3ssT/FP9O9FDcr4Xv0iXg5Wv4wU\n9YtAgfpFaHFA7xJpuupl/M+V66r9tlemu9z9yKXyrmsAmpqaEBzc9fc9ODgYTU1NbtPX1taiqqoK\nSUlJru+2bt2KPXv2ICEhAW+//TYCA90/3MwFACI9+8HNvaCUruu20+u73U5NTUVjY2OPbG+88Ua3\nz5IkdT3U3ofW1lZkZGRgy5Yt8PfveqvEs88+i9deew0A8Oqrr+Kll17Cjh073DaFwYxIzxRszTh8\n+HCf94KDg9HY2IiQkBBcunQJY8aM6TVdR0cHlixZgmXLlmHhwoWu7+9Ov2rVKixYsKDf+nDOjEjP\nOmRcMqSnp7te6rp79+5ugeo2p9OJ7OxsREdHY82aNd3uXbp0yfXnDz/8EJMnT+63TAYzIj1TaWvG\nunXrcPjwYUycOBGffvop1q1bBwBoaGjAz372MwBARUUFCgsLceTIkR5bMNauXYvY2FjExcXh6NGj\n2Lx5c79lcphJpGcqPQEQFBSEsrKyHt+PHTsWH3/8MQBg1qxZuHXrVq/59+zZI7tMBjMiPePjTEQk\nBIEeZ2IwI9Izd1szhhgGMyI94zCTiITAYSYRCYFvmiUiIXCYSURCYDAjIiFwzoyIhMCtGUQkBA4z\niUgIHGYSkRC4NYOIhMBhJhEJgcGMiITAOTMiEoJAPTO+NpuIhMBgRkRed/XqVaSmpmLixImYN28e\nmpube00XHh6O2NhYxMfHY/r06bLz301yyj02+O7MkgScHXB2/eGg3nMG9f9eGQNbVC+jo/lh9X48\nYpjsU7/v1nWWpZz8ksflvfLKK3jkkUfwyiuvID8/H9euXcPGjRt7pBs3bhw+//xzBAUFDSh/t9ox\nmGmIwcxzDGb980owa5eR4z6Py7NYLDh69Kjr/MyUlBR8/fXXPdKNGzcOp06dwqhRowaU/27850Wk\na+5WAI4BOD6gX21qakJwcDCArgOBm5qaek0nSRLmzp0LPz8/5ObmIicnR1b+uzGYEemau70ZM/5z\n3fZGt7upqalobGzskeuNN7qnkyTpP73AnioqKhAaGorLly8jNTUVFosFycnJHue/G4MZka59P+Cc\nhw8f7vPe7eFhSEgILl26hDFjxvSaLjQ0FAAwevRoLFq0CJ999hmSk5M9zn83rmYS6VqHjMtz6enp\n2L17NwBg9+7dWLhwYY80bW1taGnpmre8ceMGDh06hEmTJnmc/15cANAS+8Ge4wJA/7yyAFAjI8c4\nj8u7evUqnnzySXz77bcIDw/Hn//8ZwQGBqKhoQE5OTn4+OOP8a9//QuLFy8GAHR2diIrKwu//vWv\n3eZ32x4GMw0xmHmOwax/XglmF2TkmKioPLX1O8z87W9/i5iYGEyePBlPPfUUfvhBoFdTEulep4xr\ncHMbzGpra7F9+3acPn0aX331FRwOB/bv369V3YhIderMmfmC24HPiBEjYDQa0dbWBj8/P7S1tcFk\nMmlVNyJS3cBXMwcbt8EsKCgIL730En70ox/hwQcfRFpaGubOnds90Xt5d/6cmNJ1EZH3nSwHTh71\n8o8O/uGjp9wuAFy8eBELFizA8ePHERAQgCeeeAIZGRnIysrqyswFAHm4AOA5LgD0zysLAEdk5Jgz\ndBcATp06hR//+McYNWoUDAYDFi9ejBMnTmhVNyJSnU4WACwWCyorK/H999/D6XSirKwM0dHRWtWN\niFSnkwWAuLg4LF++HAkJCRg2bBimTp2KZ555Rqu6EZHqBn+Py1PcNKslzpl5jnNm/fPKnNlfZeRY\nMqjnzPjPi0jXdLI1g4hEN/jnwjzFYEaka+LMmWn/CqDPyjUv0utEaAMA2Mp9XQPvOFnu6xoo57M2\niLOayWA2ECK0ARAomHl7V7wP+KwN4uwz4zCTSNcGf4/LUwxmRLo2+HtcnlK+z4yIfEb5PjPPjRw5\nElevXh1weWpTFMyIiAYLHmhCREJgMCMiITCYEZEQNA1mVqsVFosFkZGRyM/P17Jor7Hb7ZgzZw5i\nYmIwadIkvPvuu76u0oA5HA7Ex8djwYIFvq7KgDU3NyMjIwNRUVGIjo5GZWWlr6skGw8N8hKnRjo7\nO50TJkxw1tTUONvb251xcXHOc+fOaVW811y6dMlZVVXldDqdzpaWFufEiROHZDucTqfz7bffdj71\n1FPOBQsW+LoqA7Z8+XLnjh07nE6n09nR0eFsbm72cY3kqampcY4bN8558+ZNp9PpdD755JPOP/7x\njz6u1dCkWc/MZrMhIiIC4eHhMBqNyMzMRHFxsVbFe01ISAimTJkCAPD390dUVBQaGhp8XCv56urq\nUFpailWrVg3q17q4c/36dRw/fhxPP/00AMBgMCAgIMDHtZLn7kODOjs7eWiQApoFs/r6eoSFhbk+\nm81m1NfXa1W8Kmpra1FVVYWkpCRfV0W2F154AZs2bcKwYUN32rSmpgajR4/GypUrMXXqVOTk5KCt\nrc3X1ZLl7kODxo4di8DAwJ6HBpFHNPubLNoG29bWVmRkZGDLli3w9/f3dXVkOXjwIMaMGYP4+Pgh\n2ysDgM7OTpw+fRq//OUvcfr0aTz00EPYuHGjr6sly8WLF/HOO++gtrYWDQ0NaG1txd69e31drSFJ\ns2BmMplgt9tdn+12O8xms1bFe1VHRweWLFmCZcuWYeHChb6ujmwnTpxASUkJxo0bh6VLl+LTTz/F\n8uXLfV0t2cxmM8xmMxITEwEAGRkZOH36tI9rJQ8PDfIezYJZQkICqqurUVtbi/b2dhQVFSE9PV2r\n4r3G6XQiOzsb0dHRWLNmja+rMyBvvvkm7HY7ampqsH//fjz66KPYs2ePr6slW0hICMLCwnDhwgUA\nQFlZGWIi9M9oAAAAhklEQVRiYnxcK3l4aJD3aPagucFgQEFBAdLS0uBwOJCdnY2oqCitiveaiooK\nFBYWIjY2FvHx8QC6ltYfe+wxH9ds4IbyFMDWrVuRlZWF9vZ2TJgwAbt27fJ1lWThoUHew2cziUgI\nQ3cpi4joLgxmRCQEBjMiEgKDGREJgcGMiITAYEZEQvj/sVzU5xpBx40AAAAASUVORK5CYII=\n",
       "text": "<matplotlib.figure.Figure at 0x11da6410>"
      }
     ],
     "prompt_number": 287
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Nu de HAC-sandwich\ntop_sandwich = np.linalg.pinv(X.T.dot(X)).dot(X.T)\nbottom_sandwich = X.dot(np.linalg.pinv(X.T.dot(X)))\n\nsandwich = top_sandwich.dot(W).dot(bottom_sandwich)\nprint 'Sandwich:\\n%s' % sandwich\n\nprint '\\nTop sandwich:\\n%s' % pandas.DataFrame(top_sandwich)\nprint '\\nBottom sandwich:\\n%s' % pandas.DataFrame(bottom_sandwich)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Sandwich:\n[[  3.55377379e-19   8.97478465e-19]\n [ -3.61400724e-20  -1.56004646e-17]]\n\nTop sandwich:\n          0             1         2         3             4         5  \\\n0  0.222222  5.551115e-17 -0.222222  0.222222  5.551115e-17 -0.222222   \n1  0.000000  2.222222e-01  0.222222  0.000000  2.222222e-01  0.222222   \n\n          6             7         8  \n0  0.222222  5.551115e-17 -0.222222  \n1  0.000000  2.222222e-01  0.222222  \n\nBottom sandwich:\n          0             1\n0  0.222222  5.551115e-17\n1  0.000000  2.222222e-01\n2 -0.222222  2.222222e-01\n3  0.222222  5.551115e-17\n4  0.000000  2.222222e-01\n5 -0.222222  2.222222e-01\n6  0.222222  5.551115e-17\n7  0.000000  2.222222e-01\n8 -0.222222  2.222222e-01\n"
      }
     ],
     "prompt_number": 288
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Dit lijkt niet te kloppen... Hele kleine getallen en nog negatieven op de diagonaal ook! Wat als we de HC sandwich doen met alleen de diagonaal, zoals in de meeste teksten (HC_3)"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Nu de HC-sandwich\nW_ = np.diag(np.diag(W))\nplt.imshow(W_, interpolation='nearest')\nplt.title('W_', fontsize=24)\nplt.colorbar()\ntop_sandwich = np.linalg.pinv(X.T.dot(X)).dot(X.T)\nbottom_sandwich = X.dot(np.linalg.pinv(X.T.dot(X)))\n\nsandwich = top_sandwich.dot(W_).dot(bottom_sandwich)\nprint 'Sandwich:\\n%s' % sandwich\n\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Sandwich:\n[[ 0.0329194  -0.012773  ]\n [-0.012773    0.08903137]]\n"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAEQCAYAAADoGVLaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHQFJREFUeJzt3X9QVNfZB/DvAttIoY1BEZC1s0ZUQAFXUWLjRsyAGjsh\ngkxHQyQxioypUTP9I80ffVv7zhgzaaaxZpJpMmrLmEjeqa2/hhLGmlWD1TVC1QTaqGEb8EfERBTY\nKOy67x+4W5Yfd+/dXe7u2fv9zNwZ2Hvu3od0fHrOueeeR+dyuVwgIgpzUaEOgIhIDiYrIhICkxUR\nCYHJioiEwGRFREJgsiIiITBZCaq9vR1RUVGIiorCgQMHhm23bt06T7u//vWvw7Z78cUXERUVhays\nrJEIlyhgTFaCSkxMRHp6OgDg2LFjw7brf05Ou/z8/OAESBRkTFYCmz9/PoDhk9A333yD5uZmJCUl\nSbbr6OjAZ599BgB47LHHRiBSosAxWQnMnVj++c9/oru7e9D548ePAwB+8pOfYOrUqTh79iw6OzuH\nbOdyuaDT6TwJkCjcMFkJzJ1YnE4n6uvrB513Jyuz2Yx58+bh3r17ku2mTJmCcePGjWDERP5jshLY\n+PHj8fDDD8Plcg05xDt27Bh0Oh3MZjPMZrPns6HaARwCUniLCXUAFJj58+fjyy+/HJSEurq60NjY\niOTkZE9CAwYnK7vdjoaGBs93EYUr9qwE5+4Nffrpp7h7967n8xMnTuDevXueHtWkSZOQnJyMM2fO\n4M6dO17tHA4H56so7DFZCc6dYO7evYtTp055PnfPQ/Uf2pnN5mHbGY1GpKamqhEykV+YrARnNBph\nMBgGzVu5f3b3rABg3rx5Xuf6/8xeFYU7JqsI4O49uRNPT08PrFYrRo8e7bUifeAke09Pj6eXxWRF\n4Y7JKgK4E83JkyfhdDphtVpx9+5dPProo17tsrOzER8fj1OnTsHhcOD06dO4c+cOdDodnwRS2GOy\nigDuRNPd3Y0zZ854ra/qLyoqCnPnzkVXVxcaGho87VJTUzFx4kR1gyZSiMkqAkydOhXjxo3zzFtJ\nrZvqP2TkfBWJhMkqQriTkMViwT/+8Q/ExsYiNzd3UDv3JLvFYsGJEye8riUKZ0xWEcLdO6qtrcXt\n27eRl5eHmJjBa37z8vKg1+s97bi+ikTBZBUh3L2je/fuARg8X+U2atQozJo1y9MuKSkJU6ZMUSdI\nogAwWUWI6dOnIyEhATqdzvM+4HDc53y1IwonOhY5JSIRsGdFREJgsiIivzz//PNISkoadt/+/fv3\nIycnByaTCbNmzcKRI0c854xGI7Kzs2EymTBnzhxZ9+MwkIj8cvz4ccTHx6O8vBznz58fdL67uxtx\ncXEAgPPnz6O4uBgXL14EAEycOBFnzpxBQkKC7PsFtJ+VTqcL5HIiClAgfY2EhATcvHlTdvuHHnoI\n3377red3s9kMm802bHt3ogL69lcbO3as13mlsQe++d6bCv9j1f4aWPxrZddsUth+xFkA5Ic4hmCw\nYPDfYQFwFICS/yP6n+CE4zcLxP/fwwLlf8PmgO548+ZN/FrivO3+4WZRkNjc9u3bh1deeQVXr15F\nXV2d53OdToeCggJER0ejsrISFRUVPr+LO4XSAPkQ/x8+ySWVANLuH24WP75/6dKlWLp0KY4fP46V\nK1fi3//+NwCgvr4eKSkpaG9vR2FhIdLT030uo+EEO5GG6RUcgTCbzXA4HPjmm28AACkpKQD66l8W\nFxfDarX6/A71k1Vavuq3DD5jqAMIEmOoAwgSY6gDCAJjSO4ao+BQ6tKlS555Kfc+/2PGjIHdbveU\nhOvu7kZdXZ2sSuDqDwOZrMKIMdQBBIkx1AEEgTEkd40N4NoVK1bg6NGjuHHjBiZMmIDNmzejt7cX\nAFBZWYm9e/eiqqoKer0e8fHxqK6uBgBcu3YNJSUlAACHw4GysjIsXLjQ5/0CWrqg0+mUT7D7I+wm\n2InCweaAngbqdDq8p6B9BQJ7+hgoTrATaZhICUCkWIkoyAKdOFeTzwn22tpapKenY/LkyXjttdfU\niImIVDKSE+zBJpmsnE4n1q9fj9raWjQ1NWHPnj1obm5WKzYiGmFqLV0IBslkZbVakZaWBqPRCL1e\nj+XLl2P//v1qxUZEI0ykZCXZu7t8+TImTJjg+d1gMHhV8wXQ9/qMW1p+hCxNIApHNni/ABO4QJYu\nqE0yWcl6UVnpe35E5CcjvNdjHQ34G8NhLkouyVhTU1PR2trq+b21tRUGg2HEgyIidYTD8E4uyTmr\n3NxcXLhwATabDT09Pfjwww9RVFSkVmxENMJEehooGUNMTAzeeustLFq0CE6nE6tXr0ZGRoZasRHR\nCBOpZ+UzYT7xxBN44okn1IiFiFQWDj0muUSKlYiCLKJ6VkQUuSJm6QIRRTb2rIhICCIlAJFiJaIg\n0yvJAI4RC0MW7sFOpGExMfKPgXwVOX3//feRk5OD7OxsPProozh37pznnD+7uTBZEWmYPlr+MdCq\nVatQW1s77Hc//PDDOHbsGM6dO4df/vKXWLt2LQD/d3MJfBioxpbDo1S4xx0V7kEUZobqMbkddQDH\nJIZ+voqczp071/NzXl4e2traAHjv5gLAs5uLrwXnnLMi0jD9A8OfK3gAKOj3+/9e9/8+O3bswJIl\nSwDI3M1lCExWRFqmQgb4+OOPsXPnTtTX1wOQuZvLEJisiLRshDPAuXPnUFFRgdraWjz00EMA/N/N\nhRPsRFo2gtsufPXVVygpKcHu3buRlvbfQvT+7ubCnhWRlg3xlE8uX0VOf/Ob3+DmzZtYt24dAECv\n18Nqtfq9m0vgRU7xK38vl49PA4mGEHiRU5eCHZ90zSxySkShIvE0MNwwWRFpmUAZQKBQiSjoBMoA\nAoVKREEXwAS72pisiLRMoAwgUKhEFHQCZQCBQiWioBMoAwgUKhEFHZcuEJEQBMoAAoVKREHHp4FE\nJASBMoBAoRJR0AmUAQQKlYiCjsNAIhKCQBlAoFCJKOhGhToA+bhTKJGWRSs4BvBVNxAANmzYgMmT\nJyMnJweNjY2ez41GI7Kzs2EymTBnzhxZoTJZEWlZANsa+6obWFNTg4sXL+LChQt49913PTuGAn0b\n/1ksFjQ2NsJqtcoKlcmKSMsCSFZms9lTBGIoBw4cwLPPPgugr25gR0cHvv76a895pbuOijFnpcaW\nw6NVuEeHCvcgUkLiaaDlP4DlK/+/eqj6gJcvX0ZSUhJ0Oh0KCgoQHR2NyspKVFRU+Pw+MZIVEY0M\niQyQP6nvcNt8XPnXD9d7+uSTTzB+/Hi0t7ejsLAQ6enpMJvNkt/FYSCRlo1gKa6B9QHb2tqQmpoK\nABg/fjwAIDExEcXFxbLmrZisiLTsAQWHQkVFRaiqqgIAnDx5EqNHj0ZSUhLsdjs6OzsBAN3d3air\nq5N8oujGYSCRlgWQAXzVDVyyZAlqamqQlpaGuLg47Nq1CwBw7do1lJSUAAAcDgfKysqwcOFCn/cT\no26gGjjBTsIJQt3AbQrab2TdQCIKFb4bSERCECgDCBQqEQWdQBlA8mlga2srFixYgGnTpmH69On4\n/e9/r1ZcRKSGAN4NVJtkXtXr9fjd736HGTNmoKurC7NmzUJhYSEyMjLUio+IRlKk7LqQnJyMGTNm\nAADi4+ORkZGBK1euqBIYEalgBBeFBpvsEGw2GxobG5GXlzfgjKXfz8b7BxEFn+3+EURhMLyTS1ay\n6urqQmlpKbZt24b4+PgBZ/ODHxURDcEI787A0cC/Mgx6THL5DLW3txfLli3DM888g6VLl6oRExGp\nJVKSlcvlwurVq5GZmYlNmzapFRMRqUWgYaDkBHt9fT12796Njz/+GCaTCSaTSXJnQCISzCgFR4hJ\n9qzmzZuHe/fuqRULEalNoJ6VQCNWIgo6gTKAQKESUdAJlAEECpWIgk6gDCBQqEQUdALNWXFbYyIt\nC/B1m9raWqSnp2Py5Ml47bXXBp2/efMmiouLkZOTg7y8PHz++eeyrx2IyYpIywLYg93pdGL9+vWo\nra1FU1MT9uzZg+bmZq82W7ZswcyZM3H27FlUVVVh48aNsq8diMmKSMsC6FlZrVakpaXBaDRCr9dj\n+fLl2L9/v1eb5uZmLFiwAAAwdepU2Gw2XL9+Xda1Q4VKgEr7o6eqcI/LKtyDIoZEBrCcACz/GP78\nUEVMT5065dUmJycHf/nLXzBv3jxYrVb85z//QVtbm6xrFYRKRBFPqsjpY32H2+Y3vM/3FYyR9otf\n/AIbN26EyWRCVlYWTCYToqOjZV2rIFQiinSuAJ4GDixi2traCoPB4NXmBz/4AXbu3On5feLEiZg0\naRK+++47n9cOxDkrIg1zxsg/BsrNzcWFCxdgs9nQ09ODDz/8EEVFRV5tbt26hZ6eHgDAe++9h/nz\n5yM+Pl7WtQOxZ0WkYUMlIbliYmLw1ltvYdGiRXA6nVi9ejUyMjLwhz/8AUBfodOmpiY899xz0Ol0\nmD59Onbs2CF5rRQWOVUVJ9gpmAIvcnrL8T3Z7R+M6WGRUyIKDWe0OEvYmayINMwp0Ps2TFZEGuZg\nsiIiETgFSgHiREpEQcdhIBEJgcmKiIRwF/KXLoQakxWRhnHOioiEwGEgEQmByYqIhMB1VkQkBM5Z\nEZEQOAwkIiH0cOkCEYmAc1ZEJASR5qy4rTGRhjkRLfsYiq9Cpb/97W9hMpk8BSNiYmLQ0dEBADAa\njcjOzobJZMKcOXN8xsqdQlXFnUIpmALfKXS/a6Hs9k/p6rzu53Q6MXXqVBw+fBipqamYPXs29uzZ\nM+z2xIcOHcKbb76Jw4cPA+grHnHmzBkkJCTIur84fUAiCjqpOavPLN/gc8u3w57vX6gUgKdQ6XDJ\n6oMPPsCKFSu8PlOSbJmsVMVeD4WXnqHqwt83JX88puSP9/z+f5svep1XUqjUbrfjo48+wttvv+35\nTKfToaCgANHR0aisrERFRYVkrExWRBoWyDorJYVKDx48iHnz5mH06NGez+rr65GSkoL29nYUFhYi\nPT0dZrN52O/gBDuRhjkQLfsYSE6RU7fq6upBQ8CUlBQAQGJiIoqLi2G1WiVjZbIi0jAnYmQfA8kt\nVHrr1i0cO3YMTz31lOczu92Ozs5OAEB3dzfq6uqQlZUlGSuHgUQaFsgwUE6RUwDYt28fFi1ahNjY\nWM+1X3/9NYqLiwEADocDZWVlWLhQ+skkly4QCSvwpQvvuJ6T3X6d7o8sckpEocEXmYlICHclli6E\nG1kT7E6nEyaTCU8++eRIx0NEKgr0dRs1yepZbdu2DZmZmZ7ZeyKKDOGQhOTy2bNqa2tDTU0N1qxZ\nE9LJNSIKvkDWWanNZ8/qpZdewuuvv47bt28P08LS72fj/YOIgs92/wgekbaIkYz00KFDGDduHEwm\nEywWyzCt8oMeFBENxQjvzsDRgL9RpGGgZLI6ceIEDhw4gJqaGty5cwe3b99GeXk5qqqq1IqPiEaQ\nSMlKcs5qy5YtaG1tRUtLC6qrq/H4448zURFFkLv4nuwj1BQNWJW8ZU1E4S9i5qz6mz9/PubPnz+S\nsRCRykQaBoqTVoko6JisiEgI4bB+Si4mKyINi8g5KyKKPCINA7lTKJGG9eB7so+h+KobCAAWiwUm\nkwnTp09Hfn6+omv7Y8+KSMMCmbNyOp1Yv369V93AoqIir1JcHR0d+NnPfoaPPvoIBoMBN27ckH3t\nQOxZEWlYIHuw968bqNfrPXUD+/vggw+wbNkyTyGJsWPHyr52IPasyA95Ktxj6PpzFFxSc1Y3LE24\nYWka9rycuoEXLlxAb28vFixYgM7OTmzcuBErV65UVHPQjcmKSMOkktVD+Vl4KP+/FWf+tXmv13k5\nb7T09vaioaEBf//732G32zF37lw88sgjfr0Nw2RFpGGBzFnJqRs4YcIEjB07FrGxsYiNjcVjjz2G\ns2fPwmAwyK456MY5KyING+m6gU899RQ++eQTOJ1O2O12nDp1CpmZmbJrDvbHnhWRhg23JEEOOXUD\n09PTsXjxYmRnZyMqKgoVFRXIzMwEgCGvlcK6geQHTrCHh8DrBs5z1clu/4luIesGElFo8HUbIhKC\nSK/bMFkRaRiTFREJgcmKiIQgUvl4JisiDWPPioiEwGRFRELgtsZEJASusyIiIXAYSERCYLIiIiHc\n7Ql9WXi5mKyINMzpECcFiBMpEQWd08FhIBEJgMmKiITg6BUnWXFbYyINu+eMkX0MRW6h0tOnTyMm\nJgZ79/636ITRaER2djZMJhPmzJnjM1b2rIi0LIBhoNxCpU6nEy+//DIWL17s9blOp4PFYkFCQoKs\n+7FnRaRld2LkHwPILVS6fft2lJaWIjExcdA5Jdsks2dFflBjf/RYFe7xnQr3CHMOiXNWC3DaMuxp\nOYVKL1++jP379+PIkSM4ffq0V71AnU6HgoICREdHo7KyEhUVFZKhMlkRaZlUspqZ33e4vb3Z67Sc\nQqWbNm3C1q1bodPp4HK5vHpS9fX1SElJQXt7OwoLC5Geng6z2TzsdzFZEWmZVLLyQU6R0zNnzmD5\n8uUAgBs3buBvf/sb9Ho9ioqKkJKSAgBITExEcXExrFarZLLinBWRlvUqOAaQU6j0yy+/REtLC1pa\nWlBaWop33nkHRUVFsNvt6OzsBAB0d3ejrq4OWVlZg2/SD3tWRFrm9P9SOUVOh3Pt2jWUlJQAABwO\nB8rKyrBw4ULJ+7HIKYUpTrD7FniRUxxVcP18HYucElGI3Al1APIxWRFpWQAT7GpjsiLSMoGSlc+n\ngR0dHSgtLUVGRgYyMzNx8uRJNeIiIjU4FBwh5rNntXHjRixZsgR//vOf4XA40N3drUZcRKSGIZYk\nhCvJZHXr1i0cP34cf/rTn/oax8TgwQcfVCUwIlJBAEsX1CaZrFpaWpCYmIhVq1bh7NmzmDVrFrZt\n24bvf//7/VpZ+v1svH8QUfDZ7h9BFAbDO7kk56wcDgcaGhrwwgsvoKGhAXFxcdi6deuAVvn9DuNI\nxEhEAPr+feX3O4LgjoIjxCSTlcFggMFgwOzZswEApaWlaGhoUCUwIlKBQBPskskqOTkZEyZMwBdf\nfAEAOHz4MKZNm6ZKYESkAoGSlc+ngdu3b0dZWRl6enowadIk7Nq1S424iEgNYZCE5PKZrHJycnD6\n9Gk1YiEitUXK0gUiinCRsnSBiCJcGDzlk4vJikjLImnOiogimEBzVtzWmEjLnAqOIfgqcrp//37k\n5OTAZDJh1qxZOHLkiOxrB+JOoRSmuFOob0HYKfRFBddv994p1Ol0YurUqV5FTvfs2eNV5LS7uxtx\ncXEAgPPnz6O4uBgXL16Ude1AHAYSaZnUnNVlC3DFMuzp/kVOAXiKnPZPOO5EBQBdXV0YO3as7GsH\nYrKiMDXyvZ5fYbPvRgHaHO4jD6k5q3H5fYfbp97/veQUOQWAffv24ZVXXsHVq1dRV1en6Nr+OGdF\npGV3FRwDyClyCgBLly5Fc3MzDh48iJUrV/o9dGXPikjLRrjIaX9msxkOhwPffvstDAaDomsB9qyI\ntG2Ei5xeunTJ05Ny79gyZswYWdcOxJ4VkZaNcJHTvXv3oqqqCnq9HvHx8aiurpa8VgqXLpBmiT/B\nHoSlC08quP4gi5wSUajwdRsiEoJAr9swWRFp2RBLEsIVkxWRlnEYSERC4DCQiITAnUKJSAgcBhKR\nEJisiEgInLMiIiFw6QIRCYHDQCISAoeBRCQELl0gIiFwGEhEQhAoWXGnUCItC2CnUMB37b9//etf\nmDt3LkaNGoU33njD65zRaER2djZMJhPmzJnjM1T2rIi0LICeldPpxPr1671q/xUVFXnt+DlmzBhs\n374d+/btG3S9TqeDxWJBQkKCrPuxZ0VEfulf+0+v13tq//WXmJiI3Nxc6PX6Ib9Dyc6j7FkR0TAs\n94+h+VP7rz+dToeCggJER0ejsrISFRUVku2ZrEizwr4Aacjl3z/cvPesl1s3cDj19fVISUlBe3s7\nCgsLkZ6eDrPZPGx7DgOJNM3/GXaldQMHSklJAdA3VCwuLobVapVsz2RFpGkOBYc3JbX/Bs5N2e12\ndHZ2AgC6u7tRV1eHrKwsyUg5DCTSNP/ft5FTN/DatWuYPXs2bt++jaioKGzbtg1NTU24fv06SkpK\nAAAOhwNlZWVYuHCh5P1YN5BIWEGoG4hrCq5IZt1AIgoVcd5kZrIi0jRx3rdhsiLSNHF6Vj6fBr76\n6quYNm0asrKy8PTTT+PuXYG2FiQiH/x/Gqg2yWRls9nw3nvvoaGhAefPn4fT6UR1dbVasRHRiAvw\nTWYVSQ4Df/jDH0Kv18NutyM6Ohp2ux2pqalqxUZEI+67UAcgm2SySkhIwM9//nP86Ec/QmxsLBYt\nWoSCgoIBrSz9fjbeP4go+Gz3j2AK/fBOLslh4KVLl/Dmm2/CZrPhypUr6Orqwvvvvz+gVX6/wzgS\nMRIRgL5/X/kY/M5eIMQZBkomq08//RQ//vGPMWbMGMTExKCkpAQnTpxQKzYiGnERMsGenp6OkydP\n4rvvvoPL5cLhw4eRmZmpVmxENOLE6VlJzlnl5OSgvLwcubm5iIqKwsyZM7F27Vq1YiOiERf6HpNc\nfDeQSFjBeDdwr4IrlvHdQCIKlQhZukBEkS70c1FyMVkRaZo4c1Yh2CnUpv4tg84W6gCCxBbqAILE\nFuoAgsAWovuK8zSQycovtlAHECS2UAcQJLZQBxAEthDdN7B1Vr6KnALAhg0bMHnyZOTk5KCxsVHR\ntf1xD3YiTfO/Z+UuclpbW4umpibs2bMHzc3NXm1qampw8eJFXLhwAe+++y7WrVsn+9qBmKyINM3/\nnpWcIqcHDhzAs88+CwDIy8tDR0cHrl27JuvagYIwwb7Zd5NBjgZ+25CLhL8B4N8RTkLxN7wsu2V8\nfLzX73KKnA7V5vLly7hy5YriAqkBJatQLhAjosAE+u9XbpHTYOUJLl0gIr/IKXI6sE1bWxsMBgN6\ne3sVF0jlnBUR+UVOkdOioiJUVVUBAE6ePInRo0cjKSlJUYFUN/asiMgvcoqcLlmyBDU1NUhLS0Nc\nXBx27dolea2UgF5kVqq2thabNm2C0+nEmjVr8PLL8if3wkVrayvKy8tx/fp16HQ6rF27Fhs2bAh1\nWH5xOp3Izc2FwWDAwYMHQx2OXzo6OrBmzRp8/vnn0Ol02LlzJx555JFQh6XIq6++it27dyMqKgpZ\nWVnYtWsXHnjggVCHFX5cKnE4HK5Jkya5WlpaXD09Pa6cnBxXU1OTWrcPmqtXr7oaGxtdLpfL1dnZ\n6ZoyZYqQf4fL5XK98cYbrqefftr15JNPhjoUv5WXl7t27Njhcrlcrt7eXldHR0eII1KmpaXFNXHi\nRNedO3dcLpfL9dOf/tT1xz/+McRRhSfV5qz8WVcRjpKTkzFjxgwAfY9yMzIycOXKlRBHpVxbWxtq\namqwZs0aYZ/q3rp1C8ePH8fzzz8PoG9o8eCDD4Y4KmX6F2VxOBwsyiJBtWQ13HoLkdlsNjQ2NiIv\nLy/UoSj20ksv4fXXX0dUlLjPWFpaWpCYmIhVq1Zh5syZqKiogN1uD3VYivQvyjJ+/HiMHj16iKIs\nBKiYrOSuyRBFV1cXSktLsW3btkGL5cLdoUOHMG7cOJhMJmF7VQDgcDjQ0NCAF154AQ0NDYiLi8PW\nrVtDHZYi8oqyEKBispKzJkMUvb29WLZsGZ555hksXbo01OEoduLECRw4cAATJ07EihUrcOTIEZSX\nl4c6LMUMBgMMBgNmz54NACgtLUVDQ0OIo1KGRVnkUy1Z+bOuIhy5XC6sXr0amZmZ2LRpU6jD8cuW\nLVvQ2tqKlpYWVFdX4/HHH/eshRFJcnIyJkyYgC+++AIAcPjwYUybNi3EUSnDoizyqbbOyp91FeGo\nvr4eu3fvRnZ2NkwmE4C+R8+LFy8OcWT+E3mIvn37dpSVlaGnpweTJk3yrOMRBYuyyKfqOisiIn+J\n+yiIiDSFyYqIhMBkRURCYLIiIiEwWRGREJisiEgI/w+dHHblTqYDrQAAAABJRU5ErkJggg==\n",
       "text": "<matplotlib.figure.Figure at 0xf4a5fd0>"
      }
     ],
     "prompt_number": 289
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "R..."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%load_ext rmagic",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "The rmagic extension is already loaded. To reload it, use:\n  %reload_ext rmagic\n"
      }
     ],
     "prompt_number": 290
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "X1 = X[:, 0]\nX2 = X[:, 1]",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 291
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%R -i X1 -i X2 -i Y\nmodel = lm('Y ~ X1 + X2 - 1')\nprint(summary(model))\n\nlibrary(sandwich)\n\nhc_models = c(\"HC3\", \"const\", \"HC\", \"HC0\", \"HC1\", \"HC2\", \"HC4\", \"HC4m\", \"HC5\")\n\n\nfor (i in 1:length(hc_models)) {\n    print(vcovHC(model, type=hc_models[i]))\n}",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": "\nCall:\nlm(formula = \"Y ~ X1 + X2 - 1\")\n\nResiduals:\n    Min      1Q  Median      3Q     Max \n-0.6018 -0.4266 -0.1807  0.0266  1.1650 \n\nCoefficients:\n   Estimate Std. Error t value Pr(>|t|)  \nX1   0.7799     0.3059   2.549   0.0381 *\nX2  -0.8400     0.3059  -2.746   0.0287 *\n---\nSignif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1\n\nResidual standard error: 0.562 on 7 degrees of freedom\nMultiple R-squared:  0.5728,\tAdjusted R-squared:  0.4507 \nF-statistic: 4.692 on 2 and 7 DF,  p-value: 0.05097\n\n            X1          X2\nX1  0.05441779 -0.02111455\nX2 -0.02111455  0.14717431\n            X1          X2\nX1  0.09358095 -0.04679048\nX2 -0.04679048  0.09358095\n           X1          X2\nX1  0.0329194 -0.01277300\nX2 -0.0127730  0.08903137\n           X1          X2\nX1  0.0329194 -0.01277300\nX2 -0.0127730  0.08903137\n            X1          X2\nX1  0.04232495 -0.01642243\nX2 -0.01642243  0.11446891\n            X1          X2\nX1  0.04232495 -0.01642243\nX2 -0.01642243  0.11446891\n            X1          X2\nX1  0.04232495 -0.01642243\nX2 -0.01642243  0.11446891\n            X1          X2\nX1  0.05441779 -0.02111455\nX2 -0.02111455  0.14717431\n            X1          X2\nX1  0.03732709 -0.01448322\nX2 -0.01448322  0.10095209\n"
      }
     ],
     "prompt_number": 292
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Groter, maar nog steeds veel _kleiner_ dan de OLS_vars..."
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Misschien werkt het alleen bij een groot aantal replicaties?"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "X = np.array([[1, 0.5],\n              [0.5 ,1],\n              [-0.5, 0.5]])\n\nX = np.tile(X, (3, 1))\n\nbeta_real = np.array([[1, -1]])\nY_real = np.dot(X, beta_real.T)\n\nar = 0.5\n\nresids = []\nbetas = []\nWs = []\n\nfor i in np.arange(500):\n    noise = np.random.randn(X.shape[0] + 1)\n    noise = [ar * noise[i-1] + (1-ar) * noise[i] for i in np.arange(1, len(noise))]\n    noise = np.array(noise)[:, np.newaxis]\n    Y = Y_real + noise\n    calc_beta = np.linalg.pinv(X.T.dot(X)).dot(X.T)\n    beta = calc_beta.dot(Y)\n    Y_ = X.dot(beta)\n    resid = Y - Y_\n    \n    \n    betas.append(beta)\n    resids.append(resid)\n    Ws.append(resid.dot(resid.T))\n\nW = np.sum(Ws, 0) / (len(Ws) - 1)\nplt.imshow(W, interpolation='nearest')\nprint 'beta:\\n%s' % np.mean(betas, 0)\nprint '\\nols_var:\\n%s' % str(np.mean([np.sum(r**2) / (X.shape[0] - X.shape[1]) for r in resids]) * np.linalg.pinv(X.T.dot(X)))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "beta:\n[[ 0.98694206]\n [-1.0181403 ]]\n\nols_var:\n[[ 0.13605541 -0.0680277 ]\n [-0.0680277   0.13605541]]\n"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAAD5CAYAAAAURMgdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADU5JREFUeJzt3X9slAWex/HPYEuiNNbggV079egVvM4U2g6U4LlGfpym\nYrYskZGDlvTEon8QN8CZ8NdtUrMJQoiRArlcYgxIQMvGP4Rj2ea2QdYu3YJmusYELxCdSYZW/rmk\naK067fjcHyLC7nWmM/M8z9Av71cySZs+ne+X6DsznT59JuA4jiMApswo9gIA3EfYgEGEDRhE2IBB\nhA0YVFLINwcCAbf2AJCHyX6pVVDYkuQ05nZ851WpsyK373ln8Je5fUMeWn/33tQPfrtTau3Mfcgv\nzub+PTn7eQ7H/kbSr3Mfcbg09+/J1XO/z+HgY5Lacp8RXJ379+SoIvn5lI4b7exSWee2nO77aqBm\n0q/xVBwwiLABg3wPe0WZ3xM9sGhFsTdwyePFXsAli4q9QMFmrljm6v0Rdj7MhL282Au4pL7YCxRs\n5opHXL0/nooDBhE2YBBhAwZlDbunp0e1tbVasGCB9uzZ48dOAAqUMex0Oq2XXnpJPT09unjxot55\n5x19+umnfu0GIE8Zw75w4YLmz5+vefPmqbS0VBs2bNCJEyf82g1AnjKeUjo0NKSqqqobnweDQZ0/\nf/6WYzqv/vTxijIjv84CbkOpswNKnT2f/UBlCXsqf+SR63nfAPIzc8Ujt/y+e/SV/ZMem/GpeGVl\npZLJ5I3Pk8mkgsGgCysC8FLGsJuamnT58mUlEgmlUikdP35ca9as8Ws3AHnK+FS8pKREBw8eVHNz\ns9LptDo6OhQKhfzaDUCesv499urVq7V6tfd/twrAPZx5BhhE2IBBhA0YRNiAQYQNGETYgEGEDRgU\nKORtdAOBgN52vL/md/Q+7/+i7O2R9Z7PeC7Q7fkMqceHGbU+zPjGhxkP+DBjtof3PWPSNwzgERsw\niLABgwgbMIiwAYMIGzCIsAGDCBswiLABgwgbMIiwAYMIGzCIsAGDCBswiLABgwgbMIiwAYMIGzCI\nsAGDCBswiLABgwgbMIiwAYMIGzCIsAGDsr7xfTatv3vPjT0ySo38i+czrgZ+6/mMvziDns9oDPR5\nPkP6kw8z1vkwI+/3ysjBKz7M+Fs8YgMGETZgEGEDBhE2YBBhAwYRNmAQYQMGETZgUMawk8mkVq5c\nqbq6Oi1cuFD79+/3ay8ABch45llpaalef/11NTY2anR0VEuWLNGTTz6pUCjk134A8pDxEbuiokKN\njY2SpLKyMoVCIQ0PD/uyGID8Tflc8UQiocHBQS1btuzWL7zd+dPHi1b8cAPggcT1W3ZTCnt0dFTR\naFRdXV0qKyu79YutnTmtBiBf867ffvTHSY/M+qr4+Pi41q1bp02bNmnt2rUFrwbAexnDdhxHHR0d\nCofD2r59u187AShQxrDPnTuno0eP6v3331ckElEkElFPT49fuwHIU8afsR977DF9//33fu0CwCWc\neQYYRNiAQYQNGETYgEGEDRhE2IBBhA0YFHAcJ++rpgcCAUnvu7jOZJZ7PuEvzj96PmMicNnzGZ87\nv/B8xvrQf3k+Q//zB+9n6AHvR+yr9+6+twc0Wb48YgMGETZgEGEDBhE2YBBhAwYRNmAQYQMGETZg\nEGEDBhE2YBBhAwYRNmAQYQMGETZgEGEDBhE2YBBhAwYRNmAQYQMGETZgEGEDBhE2YBBhAwYRNmCQ\nC28YkHJxncn0+jBjsecTfuts8XzG54FTns+41/lXz2dsDfzK8xn+/H/l5b9jFm8YANxJCBswiLAB\ngwgbMIiwAYMIGzCIsAGDphR2Op1WJBJRS0uL1/sAcMGUwu7q6lI4HL5+QgqA213WsK9cuaLTp09r\ny5Ytk57lAuD2UpLtgB07dmjv3r368ssvJzniNzd9/Lik5a4sBuCvfSCpb0pHZgz71KlTmjt3riKR\niM6ePTvJUb/ObTcAeXr8+u1HuyY9MuNT8f7+fp08eVLV1dXauHGjzpw5o/b2dpeWBOCVjGHv2rVL\nyWRS8Xhc3d3dWrVqlY4cOeLXbgDylNPvsXlVHJgesr549qPly5dr+XJeGAOmA848AwwibMAgwgYM\nImzAIMIGDCJswCDCBgwq/Lrih334i6/n4t7PUMz7EbXrPB/xH58+5/mMuYG3PJ/xc6fc8xk/C3zo\n+QxvPcx1xYE7CWEDBhE2YBBhAwYRNmAQYQMGETZgEGEDBhE2YBBhAwYRNmAQYQMGETZgEGEDBhE2\nYBBhAwYRNmAQYQMGETZgEGEDBhE2YBBhAwYRNmAQYQMGFf6GATrt4jqT+XsfZoR9mPEHH2bM9nzC\nF84/ez7jPwPXPJ/xd06H5zN+NXTQuzsP3s0bBgB3EsIGDCJswCDCBgwibMAgwgYMImzAIMIGDMoa\n9sjIiKLRqEKhkMLhsAYGBvzYC0ABSrIdsG3bNj399NN69913NTExoa+//tqPvQAUIGPY165dU19f\nn956660fDi4pUXl5uS+LAchfxrDj8bjmzJmjzZs36+OPP9aSJUvU1dWle+6556ajjt308SJJ9Z4s\nCtzx/vzBD7cpyPgz9sTEhGKxmLZu3apYLKZZs2Zp9+7df3VU2003ogY880+PS//27z/dMsgYdjAY\nVDAY1NKlSyVJ0WhUsVjMvUUBeCJj2BUVFaqqqtKlS5ckSb29vaqrq/NlMQD5y/qq+IEDB9TW1qZU\nKqWamhodOnTIj70AFCBr2A0NDfrwww/92AWASzjzDDCIsAGDCBswiLABgwgbMIiwAYMIGzAo6++x\nswqudmGNLK78r/czlPf7JuTgAR9m/N7zCT8LeH9ewwFnj+czmgNvej7jPSfp2X2vzfA1HrEBgwgb\nMIiwAYMIGzCIsAGDCBswiLABgwgbMIiwAYMIGzCIsAGDCBswiLABgwgbMIiwAYMIGzCIsAGDCBsw\niLABgwgbMIiwAYMIGzCIsAGDCBswKOA4Tt5Xyg8EAqpwPnNzn//X1UC15zOkV7wfsa/T+xnbx7yf\noSHvR1yp8nzEe5W/9HxGfeC/Pbvvf5A0Wb48YgMGETZgEGEDBhE2YBBhAwYRNmAQYQMGZQ371Vdf\nVV1dnRYtWqTW1lZ99913fuwFoAAZw04kEnrjjTcUi8X0ySefKJ1Oq7u726/dAOSpJNMX7733XpWW\nlmpsbEx33XWXxsbGVFlZ6dduAPKUMezZs2fr5Zdf1kMPPaS7775bzc3NeuKJJ245ZrSz68bHM1cs\n08wVj3izKXCHG5B0forHZgz7s88+0759+5RIJFReXq5nn31Wx44dU1tb241jyjq3FbAqgKl65Prt\nR12THagsP2N/9NFHevTRR3X//ferpKREzzzzjPr7+11ZEoB3MoZdW1urgYEBffPNN3IcR729vQqH\nw37tBiBPGcNuaGhQe3u7mpqaVF9fL0l68cUXfVkMQP4y/owtSTt37tTOnTv92AWASzjzDDCIsAGD\nCBswyPewU2cH/B7pgUSxF3DJB8VewB1/nv7/DrerKELYUz135naWKPYCLukr9gLuMBC221XwVBww\niLABgwq+rjiA4pks36wnqORzpwCKi6figEGEDRhE2IBBvobd09Oj2tpaLViwQHv27PFztGuSyaRW\nrlypuro6LVy4UPv37y/2SnlLp9OKRCJqaWkp9ip5GxkZUTQaVSgUUjgc1sDA9DsBypMLhjo+mZiY\ncGpqapx4PO6kUimnoaHBuXjxol/jXfPFF184g4ODjuM4zldffeU8/PDD0/Lf4TiO89prrzmtra1O\nS0tLsVfJW3t7u/Pmm286juM44+PjzsjISJE3yk08Hneqq6udb7/91nEcx1m/fr1z+PDhgu/Xt0fs\nCxcuaP78+Zo3b55KS0u1YcMGnThxwq/xrqmoqFBjY6MkqaysTKFQSMPDw0XeKndXrlzR6dOntWXL\nlmn7241r166pr69Pzz//vCSppKRE5eXlRd4qNzdfMHRiYsK1C4b6FvbQ0JCqqn56z+NgMKihIR/e\nZ9lDiURCg4ODWrZsWbFXydmOHTu0d+9ezZgxfV9micfjmjNnjjZv3qzFixfrhRde0NiYH+8P7p6b\nLxj64IMP6r777vubC4bmw7f/qtZOZhkdHVU0GlVXV5fKysqKvU5OTp06pblz5yoSiUzbR2tJmpiY\nUCwW09atWxWLxTRr1izt3r272Gvl5OYLhg4PD2t0dFTHjh0r+H59C7uyslLJZPLG58lkUsFg0K/x\nrhofH9e6deu0adMmrV27ttjr5Ky/v18nT55UdXW1Nm7cqDNnzqi9vb3Ya+UsGAwqGAxq6dKlkqRo\nNKpYLFbkrXLj1QVDfQu7qalJly9fViKRUCqV0vHjx7VmzRq/xrvGcRx1dHQoHA5r+/btxV4nL7t2\n7VIymVQ8Hld3d7dWrVqlI0eOFHutnFVUVKiqqkqXLl2SJPX29qqurq7IW+XGqwuGFnRKaU6DSkp0\n8OBBNTc3K51Oq6OjQ6FQyK/xrjl37pyOHj2q+vp6RSIRST/8uuKpp54q8mb5m84/Jh04cEBtbW1K\npVKqqanRoUOHir1STm6+YOiMGTO0ePFiVy4YWtAfgQC4PU3fl0QBTIqwAYMIGzCIsAGDCBswiLAB\ng/4PajCfiIKeDVcAAAAASUVORK5CYII=\n",
       "text": "<matplotlib.figure.Figure at 0x11da8210>"
      }
     ],
     "prompt_number": 293
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Beta wordt uitstekende geschat.\nWe  zien de autocorrelatie duidelijk terug in de one-off-diagonal van W...\n\n\nSandwich?"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "top_sandwich = np.linalg.pinv(X.T.dot(X)).dot(X.T)\nbottom_sandwich = X.dot(np.linalg.pinv(X.T.dot(X)))\n\nsandwich = top_sandwich.dot(W).dot(bottom_sandwich)\nsandwich",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 294,
       "text": "array([[  2.58281051e-17,   7.70988212e-18],\n       [  1.28176790e-17,   4.85722573e-17]])"
      }
     ],
     "prompt_number": 294
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Het lijkt ook niet in de volgorde van dots te zitten, links-naar-rechts gaat goed\nnp.linalg.pinv(X.T.dot(X)).dot(X.T).dot(W).dot(X).dot(np.linalg.pinv(X.T.dot(X)))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 295,
       "text": "array([[  2.58281051e-17,   7.70988212e-18],\n       [  1.28176790e-17,   4.85722573e-17]])"
      }
     ],
     "prompt_number": 295
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Nee..."
    }
   ],
   "metadata": {}
  }
 ]
}