Statsmodels example: Weighted Least Squares

example_wls.ipynb

{
 "metadata": {
  "name": "example_wls"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Weighted Least Squares"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy import stats\n",
      "import statsmodels.api as sm\n",
      "import matplotlib.pyplot as plt\n",
      "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
      "from statsmodels.iolib.table import (SimpleTable, default_txt_fmt)\n",
      "np.random.seed(1024)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## WLS Estimation\n",
      "\n",
      "### Artificial data: Heteroscedasticity 2 groups \n",
      "\n",
      "Model assumptions:\n",
      "\n",
      " * Misspecificaion: true model is quadratic, estimate only linear\n",
      " * Independent noise/error term\n",
      " * Two groups for error variance, low and high variance groups"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nsample = 50\n",
      "x = np.linspace(0, 20, nsample)\n",
      "X = np.c_[x, (x - 5)**2, np.ones(nsample)]\n",
      "beta = [0.5, -0.01, 5.]\n",
      "sig = 0.5\n",
      "w = np.ones(nsample)\n",
      "w[nsample * 6/10:] = 3\n",
      "y_true = np.dot(X, beta)\n",
      "e = np.random.normal(size=nsample)\n",
      "y = y_true + sig * w * e \n",
      "X = X[:,[0,2]]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### WLS knowing the true variance ratio of heteroscedasticity"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mod_wls = sm.WLS(y, X, weights=1./w)\n",
      "res_wls = mod_wls.fit()\n",
      "print res_wls.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            WLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.847\n",
        "Model:                            WLS   Adj. R-squared:                  0.844\n",
        "Method:                 Least Squares   F-statistic:                     265.9\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           3.32e-21\n",
        "Time:                        20:53:31   Log-Likelihood:                -46.062\n",
        "No. Observations:                  50   AIC:                             96.12\n",
        "Df Residuals:                      48   BIC:                             99.95\n",
        "Df Model:                           1                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1             0.4379      0.020     22.088      0.000         0.398     0.478\n",
        "const          5.2726      0.185     28.488      0.000         4.900     5.645\n",
        "==============================================================================\n",
        "Omnibus:                        5.040   Durbin-Watson:                   2.242\n",
        "Prob(Omnibus):                  0.080   Jarque-Bera (JB):                6.431\n",
        "Skew:                           0.024   Prob(JB):                       0.0401\n",
        "Kurtosis:                       4.756   Cond. No.                         17.0\n",
        "==============================================================================\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## OLS vs. WLS\n",
      "\n",
      "Estimate an OLS model for comparison: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "res_ols = sm.OLS(y, X).fit()\n",
      "print res_ols.params\n",
      "print res_wls.params"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.43486879  5.24256099]\n",
        "[ 0.43794441  5.27260714]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare the WLS standard errors to  heteroscedasticity corrected OLS standard errors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "se = np.vstack([[res_wls.bse], [res_ols.bse], [res_ols.HC0_se], [res_ols.HC1_se], [res_ols.HC2_se], [res_ols.HC3_se]])\n",
      "se = np.round(se,4)\n",
      "colnames = ['x1', 'const']\n",
      "rownames = ['WLS', 'OLS', 'OLS_HC0', 'OLS_HC1', 'OLS_HC3', 'OLS_HC3']\n",
      "tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)\n",
      "print tabl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "=====================\n",
        "          x1   const \n",
        "---------------------\n",
        "WLS     0.0198 0.1851\n",
        "OLS     0.0233 0.2707\n",
        "OLS_HC0 0.0281 0.194 \n",
        "OLS_HC1 0.0287 0.198 \n",
        "OLS_HC3 0.029  0.2003\n",
        "OLS_HC3  0.03  0.207 \n",
        "---------------------\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate OLS prediction interval:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "covb = res_ols.cov_params()\n",
      "prediction_var = res_ols.mse_resid + (X * np.dot(covb,X.T).T).sum(1)\n",
      "prediction_std = np.sqrt(prediction_var)\n",
      "tppf = stats.t.ppf(0.975, res_ols.df_resid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Draw a plot to compare predicted values in WLS and OLS:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prstd, iv_l, iv_u = wls_prediction_std(res_wls)\n",
      "plt.figure()\n",
      "plt.plot(x, y, 'o', x, y_true, 'b-')\n",
      "plt.plot(x, res_ols.fittedvalues, 'r--')\n",
      "plt.plot(x, res_ols.fittedvalues + tppf * prediction_std, 'r--')\n",
      "plt.plot(x, res_ols.fittedvalues - tppf * prediction_std, 'r--')\n",
      "plt.plot(x, res_wls.fittedvalues, 'g--.')\n",
      "plt.plot(x, iv_u, 'g--')\n",
      "plt.plot(x, iv_l, 'g--')\n",
      "plt.title('blue: true, red: OLS, green: WLS')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<matplotlib.text.Text at 0x4b23b90>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEICAYAAACpqsStAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdYFFfXwH80KSrFAhZEFLuiJvYKsRF7i+0FbIktiZpq\n8hqNphhNYuIbNQ0TjbBgN7FgiwWwlxgFsaAIWLCBIChS935/zOcK0pfdBfT+nmee3Z25c+/ZYTlz\n59xTjIQQAolEIpGUS4xLWwCJRCKRaI9U4hKJRFKOkUpcIpFIyjFSiUskEkk5RipxiUQiKcdIJS6R\nSCTlGKnESwlnZ2f27duX57GgoCDq1KljYInKJ/JaSV50pBIvJYyMjDAyMiptMTA2Nubq1aulLYbB\nEELw7bff0qhRI6ysrKhbty6zZ88mPT1d02b8+PHMnTs3z/O3bNlC69atsbGxoXr16vTs2ZPo6GgD\nSW8YGjduzPr16zWfDx8+jLGxca591tbWqNVq/vjjD7p165ZnX+Hh4fTp04eqVatiZ2dH27Zt2blz\np96/w4uEVOISCor3yszMNKAkoFar9dr/jBkzWLFiBX5+fjx8+JCdO3eyb98+Ro4cqWmT3w32ypUr\njBs3jiVLlvDgwQOioqJ46623MDEx0bmchr7u2XFzcyMkJETzOSQkhCZNmuTa17lzZ4yNC1YhAwcO\nxMPDgzt37nD37l2WLl2KtbW13mR/EZFKvBQ5ceIEzZs3p0qVKkycOJG0tLQ82z07W352prh9+3Za\nt26NnZ0dXbp0ISwsrEjjd+/eHYBWrVpRuXJlNmzYQFBQEI6OjnzzzTfUrFmTiRMnsnr16lwzrewy\npaWl8cEHH1C3bl1q1KjBtGnTSE1NLZIM48ePZ9q0afTr149KlSoRFBREbGwsw4cPx97envr167Ns\n2TJN+8ePHzN+/HiqVKlC8+bNOXnyZJHGAbh8+TI///wzAQEBdOjQAWNjY5o1a8amTZvYtWsXQUFB\nmrZ53djOnDlDvXr1eOWVVwCoVKkSw4YNK7I55/Tp07z00ktYW1szcuRIRo0apfk7PnvdX3/9dYQQ\nLFq0iAYNGlCtWjVGjRpFQkKCpr9jx47RuXNn7OzsaN26NcHBwZpj7u7ufPrpp3Tt2hVra2s8PDyI\nj48vkpzdu3fPobAPHTrERx99lGPfwYMHNb+f/IiLiyM6OppJkyZhamqKmZkZnTt3pkuXLkWSQ1I0\npBIvJYQQBAQEsGfPHiIjI4mIiODLL78s0rnZZ4r//vsvr7/+OitWrOD+/ftMmTKFQYMGkZGRAcBb\nb73FW2+9lWc/T/4pQ0NDSU5OZsSIEQDcuXOHhIQErl27ho+PT4EzdYCPP/6YK1eucPbsWa5cucLN\nmzf5/PPPi/RdANasWcPcuXN5+PAhnTp1YuDAgbz00kvExsayb98+/ve//7Fnzx4APvvsM6Kiorh6\n9Sq7d+9m9erVOWbNBX3fffv2UadOHdq2bZtjv6OjIx07dtSMkR9t2rTh4sWLvPfeewQFBfHw4cMi\nf8f09HSGDh3KxIkTSUhIYMyYMfz11185ZM9+3X/99VeWLl3K1q1bCQkJ4datW9jZ2Wm+282bNxkw\nYACffvopCQkJLF68mOHDh+dQ1GvWrOGPP/7g7t27pKens3jxYs2xVq1asXbt2jxl7datG+Hh4SQm\nJqJWqzl16hSjRo0iMTFRs+/IkSOFKvGqVavSoEEDPD092bJlC3fu3Cny9ZIUAyEpFZydncWvv/6q\n+bxjxw7h4uIihBDiwIEDwtHRUXPMyMhIREZGaj6PHz9ezJ07VwghxNSpUzXvn9C4cWMRHBxcJDme\n7fvAgQOiQoUKIi0tTbNv1apVomvXrnmep1arRcWKFXP0ceTIEVGvXr0ijT9+/Hgxbtw4zedjx44J\nJyenHG2++uorMWHCBCGEEPXr1xe7d+/WHPPx8clxrQriiy++EB07dszz2OjRo8XkyZM1Ms2ZMyfP\ndseOHRMjR44U1atXFxYWFmL8+PHi4cOHhY4dHBwsateunWNf165dNX+7vK5706ZNxb59+zSfY2Nj\nhZmZmcjMzBSLFi0S3t7eOfrz8PAQq1evFkII4e7uLhYsWKA59tNPP4lXX321UDmf4OzsLLZs2SJO\nnz4tunTpIoRQrtGTfZaWliI9PV0Ikffv4wk3btwQb7/9tnBxcRHGxsaie/fu4vLly0WWQ1I4pqV9\nE3mRyf4Y7uTkRGxsbLH7iImJwdfXN4fJISMjg1u3bmktV/Xq1alQoUKR2t67d4+UlBTatGmj2SeE\nKJZt29HRUfM+JiaG2NhY7OzsNPuysrI0s77Y2Nhc162oVKtWLd/rEhsbS/369Qvto0OHDqxbtw5A\nM0NdsGABX331VYHnxcbGUrt27Rz7njXDPHvdo6OjGTp0aA67s6mpKXfu3CEmJoYNGzawbds2zbHM\nzEx69Oih+VyjRg3Ne0tLy2I9OTwxqTg5OWmufdeuXTX7OnTogJmZWaH91K5dW/PbvHHjBpMnT2bs\n2LEcOXKkyLJICkaaU0qRa9eu5Xhfq1atPNtZWVmRkpKi+ZxdETk5OfHJJ5+QkJCg2R4+fMioUaO0\nluvZRb2KFSvmGP/27dua99WqVcPS0pLz589rxk9MTCQpKUmr8ZycnKhXr16O75OUlMT27dsBqFmz\nZq7rVlR69OjB9evXc9nRr1+/zvHjx+nZs2eeMuVH27ZtGTp0KOfOnSu0bc2aNbl582aOfc/K/uyY\nTk5O7Nq1K8e1SElJoVatWjg5OeHt7Z3jWHJyMrNmzSpUlqLwRIkfPHhQsx7SrVs3zb7CTCl54ejo\nyJtvvlmk6yUpOlKJlxJCCH788Udu3rzJ/fv3WbBgAaNHj86zbevWrfH39ycrK4tdu3blWGCaNGkS\nv/zyCydOnEAIwaNHjwgMDCzyrMvBwYHIyMgC27Rq1Yrw8HDOnj1Lamoq8+fP1xwzNjZm0qRJvPPO\nO9y7dw9Q7LXZ7cvGxsY5ZH72OmSnffv2VK5cmW+++YbHjx+TlZXFuXPnOHXqFAAjR45k4cKFJCYm\ncuPGjRxPIIXRqFEjpk6diqenJ8ePHycrK4vw8HCGDx9O7969NbNYIQSZmZmkpqZqtvT0dA4fPsxv\nv/2m+Z4XL15k27ZtdOrUCVAWJ/Pz1ujcuTMmJiYsX76czMxMtmzZUuii7NSpU5k9e7ZG2d+7d4+t\nW7cC4OXlxbZt29izZw9ZWVmkpqYSFBSU40bx7LUtDt27d+f06dOEhIRoFiJdXV25evUqBw4cyKXE\nhRCkpaXluGaJiYnMmzePyMhI1Go1cXFxrFy5UnO9JLpBKvFSwsjICE9PT/r06YOLiwsNGzZkzpw5\nOY4/4YcffmDbtm3Y2dkREBDA0KFDNcfatGnDihUrePvtt6lSpQoNGzbE19dXc/60adOYNm1avnLM\nnz+fcePGYWdnx8aNG/N0r2vUqBGffvopvXr1onHjxnTr1i1Hm6+//poGDRrQsWNHbGxs6N27NxER\nEYAyy61cuTKurq75XofsfRkbG7N9+3bOnDlD/fr1qV69OpMnT9bM7OfNm0fdunWpV68er776KmPH\njs1xfmHfd/ny5bzxxht4eXlRuXJl+vbtS48ePdi0aVMOmRYtWoSVlZVm69WrF7a2tmzduhVXV1fN\nucOGDdPMfq9fv56v54WZmRmbN2/m999/x87ODn9/fwYMGJDDfPLsdZ85cyaDBg2iT58+WFtb06lT\nJ06cOAEos9otW7bw1VdfYW9vj5OTE999910OxZ29v2evc4sWLVizZk2+16lhw4bY29tTs2ZNjUug\nkZERHTp0IDk5mc6dO+fo+8iRI1haWmquV8WKFTE2NiYmJoZevXphY2ODq6srlpaW/PHHH/mOKyk+\nRqIkt2uJpBD8/f05f/48CxYsKG1R9M6kSZMYOXIkvXv3LlL7Dh068OabbzJu3Dg9SyZ5nilQiU+c\nOJHAwEDs7e01vscnTpzg7bffJiMjA1NTU3766SfatWtnMIElkvJKSEgIjRo1olq1avj7+/Pmm29y\n9epVHBwcSls0STmmQHPKhAkT2LVrV459s2bN4osvvuDff//l888/19lCikTyvHPp0iVNUNaSJUvY\nuHGjVOCSElOgi2G3bt1y5YWoWbMmDx48ACAxMTGX25REIsmbSZMmMWnSpNIWQ/KcUahNPDo6moED\nB2rMKTExMXTt2hUjIyPUajVHjx7N5e9aFhI7SSQSSXmkuMuUxfZOef3111m6dCnXrl1jyZIlTJw4\nMV9B5Kabbd68eaUuw/O0yespr2VZ3bSh2Er8xIkTGhe31157TePyJJFIJBLDU2wl3qBBA022tP37\n99OoUSOdCyWRSCSSolHgwuaYMWMIDg4mLi6OOnXq8Pnnn+Pj48Nbb71FWloalpaW+Pj4GErWFxZ3\nd/fSFuG5Ql5P3SGvZemjl2AfIyMjre07EolE8jwSGBjC0qV7SEszxdw8kxkz+tC/f870BdroTpnF\nUCKRSPRMYGAIM2fuJjLyaeRyZOQnALkUeXGRuVMkEolEzyxduieHAgeIjFzAsmV/l7hvqcQlEolE\nz6Sl5W30SE0teX1WqcQlEolEz5ib51342sIiq8R9SyUukUgkembGjD64uHySY5+Ly2ymTy9axsuC\nkN4pEolEYgACA0NYtuxvUlNNsLDIYvr03jrxTpFKXCKRSMoI2uhOaU6RSCSScoxU4hKJRFKOkUpc\nIpFIyjFSiUskEkk5RipxiUQiKcdIJS6RSCTlGKnEJRKJpBwjlbhEIpGUY6QSl0gkknKMVOISiURS\njpFKXCKRSMoxUolLJBJJGSA+JV6r8wpU4hMnTsTBwQFXV9cc+5ctW0bTpk1p0aIFH330kVYDSyQS\nyYtORlYGWy5uYdi6YdRfWl+rPgpU4hMmTGDXrl059h04cICtW7cSGhrKuXPn+OCDD7QaWCKRSF5k\nNoRvoPb3tVl8dDH9Gvbj2jvXtOqn0FS00dHRDBw4kLCwMABGjhzJ1KlT6dGjR/6dylS0EolEUiDR\nidFkqjNpUKWBZp9Bqt1fvnyZkJAQZs+ejYWFBYsXL6Zt27a52s2fP1/z3t3dHXd39+IOJZFIJOWa\n1MxUjlw/Qo96uSe9zrbOBAUFoQpSlWiMYivxzMxMEhISOHbsGCdPnmTkyJFcvXo1V7vsSlwikUhe\nFIQQnLh5gtVnV7M+fD2ta7Smq1NXKphUyNXWvVs33FNTYf16+PlnPvvss2KPV2wl7ujoyLBhwwBo\n164dxsbGxMfHU7Vq1WIPLpFIJM8Tv576lSXHlqAWasa2Gss/k/+hrm3d3A3Dw2H1alCpwNERxo0D\ntVqrMYutxIcMGcL+/ftxc3MjIiKC9PR0qcAlEokEqF6xOisHr6STYyeMjIzybvTOO7BxI3h5wb59\n0LRpicYscGFzzJgxBAcHEx8fj729PZ9//jleXl5MnDiRM2fOUKFCBb777rtc9m65sCmRSJ5X1ELN\n3Ud3qVGphnYd3L0LVauCiUmuQ7JQskQikeiJyPuR+Ib64nfWj/a127P2tbV5NxQCTp+G8+fB27tY\nYxjEO0UikUheFDLVmaw+s5rVZ1dzMe4iY1zHsHHkRl6q8VLuxrGxio179Wp4/BimTjWIjFKJSyQS\nST6YGJlw6tYp3u/0Pn0b9s3TwwQhYNgwCA5WXn/5Bbp2hfxs4jpGmlMkEokEyFJnYWKc205dJA4f\nhpdeAiurEsmgje6UCbAkEskLS3xKPD+e+JEOv3Xgu6PfFdw4Kkqxc+dFly4lVuDaIpW4RCJ5ocjI\nymDrpa0MXz8cl6UuHLp+iM/cP+O9Tu/lbpyUBCtXgpsbtG8PR44YXuBCkOYUiUTyQhF+N5wp26cw\nrtU4RjYfiY2FTe5Gt2/D++9DYCC4uyvBOP37Q4U8bOI6RLoYSiQSiS5ITYXff4dRo6BaNYMNK23i\nEonkhSc1M5V159bRz78fYXfCCm587x6kpOTeb2EBb71lUAWuLVKJSySSco8QgiPXjzBl+xRqf1+b\n3/79jf+4/geXKi65G6elwebNMHgwNGwIJ08aXmAdIv3EJRJJuWfp8aX8fOpnxrUax5kpZ6hjUyd3\no0uXYNkyWLcOmjdX7NwqFVSubHiBdYi0iUskknJPelY6ZsZm+SedAggJUQJyvL3B2dlgshUHubAp\nkUieS9RCzYGoA+y4soPFvRcXrKzT0/XuRaIv5MKmRCJ5rrgUd4nZ+2bj/D9nPvz7Q+pY1yFDnZG7\noVoN+/fD+PFKfu7kZIPLWlrImbhEIimTvLH1DQIvB+Lp6sm4VuNwdXDN3ejSJfD1BT8/Jb3r2LEw\nZgzU0DJNbCkjzSkSieS54dqDa9SqXAtT4wL8Lz75RPE2GTsWWrY0nHB6QipxiURSbhBCcOb2Ga4n\nXWdQ40GlIkNgYAhLl+4hLc0Uc/NMZszoQ//+3UtFFpD5xCUSSTngVvIt/MP88T3rS1JaEu90fCfv\nhkLAP/8o5pJbt2DDBp3KERgYwsyZu4mMXKDZFxn5CUCpKvLiImfiEonEIGRkZTBk3RCOXD/C0CZD\nGddqHN3qdsPY6Bn/ihs3wN9fUd6pqYqpxNsb6tfXqTweHnPYs+fLPPbPZdeuL3Q6VlHR+Ux84sSJ\nBAYGYm9vT1hYzvDV7777jg8//JC4uDiqVKlSfGklEskLhZmJGTPaz2D9a+upWKFi3o3Uanj1VSW1\n66+/Kq96Kq6Qlpa3+ktN1TKnOKVjnilQiU+YMIHp06czduzYHPuvX7/O33//Td26dfUqnEQiKX9E\nJUQBUM+uXq5jHg08Cj7Z2BjCwgxSFcfcPDPP/RYWWVr1V1rmmQL9xLt164adnV2u/e+99x7ffPON\n3oSSSCTli6S0JH4//TvdV3Wn/W/tOX7zeP6NL16E2bMVt8C8MFBZsxkz+uDi8kmOfS4us5k+vbdW\n/S1duieHAgeIjFzAsmV/F+n8ydsmazVusRc2t2zZgqOjIy0LceeZP3++5r27uzvu7u7FHUoikZRx\nohOjmb1vNjsu7+CVeq/wXqf36NewX+5alHFxsHatYue+fh08PaFt29IR+v95MjtetmwuqakmWFhk\nMX36q1rPmrUxzwQFBREUFATAnjN7tBq3WEo8JSWFr776ir//fnpnyc8In12JSySS55OKZhXpXKcz\nS/supZpVPmlbL16EDh2Uogqffw69eoFp8R3j9GFv7t+/u85MHUUxz0zeNpmI+AiszKwIGB5Al25d\neFjrIdGJ0TRr2IyYLTHFHrdYVzIyMpLo6GhatWoFwI0bN2jTpg0nTpzA3t6+2INLJJLywd1Hd6li\nWSVX4E31itV5u/3bBZ/cuLEy+7a21nr8suYOmNcNZcaMPkRGfpJDRsU886rmc0R8BMExwQC0W9GO\nB6kPaGhTn0lZrQgY748dxXcSKZYSd3V15c6dO5rP9erV459//pHeKRLJc0haZhrbIrbhe9aXkJgQ\nQiaE0NIhHzNqTIyS1tXLC551eDAyKpECh4LszXMNrsTzu6H88IMHP/zgUaB5xtLUEoAKJhUYbtqS\nN86k0GD7EehVGzy184opcGFzzJgxdO7cmYiICOrUqcOqVatyHC8wk5hEIimXhN0JY1rgNGp9X4uf\nTv7E8KbDuf7u9dwKPCkJVq2CV16BNm3g5k29yaQPd0BtKWgBs3//7ji9dQfGB2Hs/Q9deua8Zmte\nW8OrJo25/XNFFgXcpYHbUIiOhk2btL7RFTgTX7NmTYEnX716VatBJRJJ2SUqMQrHyo6cnnyaurb5\nuBEHBMCbbyoKfMYM6NcPzM31JpOu3QFLQmE3lEvxlwiJCQEUG/j6Ees1bWwtbNnZezV42kO93C6Y\n2iDD7iWSF5RMdWaeyaUGNR5UeC6Tnj3hyhWD1aAsir3ZUGhuKAMnQ9UIyLCCTf6kV7nJe7vf49j1\nYwA0tKqDz0Cf3B106KBTeaQSl0heILLUWeyP2o9vqC8HYw5yZcaV/LME3r4Nu3YpObqfxcFBr3I+\ni67dAUuC5oZSNQKclUVKo/dqEFnRjh6RXTl0oSPfcgQf+37YWtjqXR6ZO0UieQE4f+88vmd9UYWq\nqFGpBt4tvRnjOgb7is94lT1+DFu2KP7cR4/CkCHwyy96NZWURwIDQxj/90Ti7CKxTq7Bkpt1Gb8n\nAuOWrZTF3ddeAxubYvcrU9FKJJI8mbp9KjYWNni39KaFfYu8G82bB0uXQvv2ShHhIUPAysqwgpYx\nnvh1W5paMrnNZHrV70Vlc6WwcmJqIpO3TcZnoA+267Yo6wNOTiUaTypxiUSiPUePKu6BtWqVtiRl\nhnYr2nEq9hQAdhVsODzwL5q2cNfbeDKfuETyAiKE4MTNE/iG+iKE4Kf+P+XfOC5OcWnLK+S9Uye9\nyVgcDJ0J8NkoSlsLW3Ze3sns/bO5cO8CAC1SKnPwJ4Gtw1XQoxLXBqnEJZJyyrUH11CFqvA964ta\nqBnXahxeLb1yN0xLg8BAxc4dFARTp5Z63pL8KI3IzOxRlE9cAmulGPNdWC1abYpi2jAbfNrMw/aK\nN1SurBcZSoJU4hJJOSQ1M5Vuq7rRr2E/Vg1eRUfHjrmD79LTYeZMpSKOq6tSXMHXt8TRk/rEkJGZ\naqHm3N1zWJkpdv92tdppXAJb1WgNTm5w8lfWOzrqdFxdI5W4RFIOsTC1IGpmVO6qONmpUAFat4aP\nP84dCl9G0VehhinbZ/LQPAFTYcrszh8QV015inEwr8pOr528uXuGskD5xCXQwQFmzdJ6TEMilbhE\nUkYJvxuOb6gv3Z26079R/1zHNQo8IQEyM6F69dydTJmiZyl1i74KNdzsZgM1zgDwQdh0vEU7Ai/V\nxXVvGHS7lyOqsrxRYO4UiURiWO4+ussPx36gjU8bPFQeGGFEk2pNcjfMyIBt22DECHB2hp07DS6r\nPihJoYbJ2ybj/oc7/fz7kZiaCGQzz2QoJhOHOGsuf2fF9F/v4DpkilLP09VV91/EgMiZuERSRjhy\n/Qj9/PsxsPFAFvVcRI96PTAxfsaMcP06fPutUmChUSOlgLCPD+RRgas8UpLIzLwWKDXmmU0B1Bww\ngKHbX6Vr6us07vArQV55LAKXQ6QSl0jKCO1qteP6u9c1wSR5kpUFVasqPt0uLoYTTscU5EaoTaGG\niPgIbiTdAKCaRVXNAqXGPJNqy62Nh3jifNmqFBJn5UdcHPz221lWr76m1flSiUskZQQzEzPMTMyU\nD0lJijvbsx4nzs5KZGU5piRuhNl9ulcNXsWmC5vwC/Uj6v5Vhpm24EpCAut2VcT2Q8UDpywlznr4\nEMLDITQUzp1TtvBwePgwk6wsW1JTW2nVr4zYlEjKCpmZ8Pffihvgjh1w4oRSFec5w8NjDnv2fJnH\n/rns2vVFgee6/+GuMZm81vQ1qiQ8ZvDpFPqs/wfTzl0V89KgQTnSBQQGhrBs2d/ZzDO99Ro8lJWl\nJHgMDVW2sDBlu3ULmjZVTPAtWiivzZvDxIlz+PvvJ9dDRmxKJOWP8HD47TdYs0bJMe3tDcuXK2aT\n55CiuBFmn3H7D/OnUoVKmJmY5fDpXjFoBbZzF0AzR7i4Jt/Mirqso/ksycmKoj57Fs6cUV7Dw8He\nHlq2VDZPT0VhN2iQd2nR9PSSqWGpxCWSUuRS3CXq/XOCCpUrQ0iIslj5nFMUN8Lsi5R1ltThB4//\n8XqbNwgYHvA06ZSFrbLIayBu3VIU9b//Pt1u31Zm061aKdvYsYrCLk4Cw/yuR1GR5hSJxFBkZYHJ\n09lmpjqTRssaETA8gI6OHUtRMMOSl03cxWU2P/zwKl17tmJd+Do+2vsRiamJ2Btbs+aMC6/YvYTR\nb78bRD4hlJKhp0/DP/8or//+q3h1vvRSzq1Roxx/Uq3IeT2kOUUiKVtkZsLeveDnp2iC8HDNYuWf\nF/6kVuVaBlXghk4ulRdPxpsa+DLJFe5jhikLXvmB/v27czzmCHuPqPgpshmbH51gxeMu2I6ZAAMH\n6kUWIZR8YKdOKQr7idI2N1fKhr78spJq5uWXwdEx9zqzLsjuVrl7d/HPL3AmPnHiRAIDA7G3tycs\nLAyADz/8kO3bt1OhQgVcXFxYtWoVNs88O8iZuOSF58wZZYFyzRolx7S3N4wapYmqFELQ4bcOzO42\nmyFNhhhEpLxnwEqVdkMrciEEr6x+RWMyGdFshBI1mZSkKOyRI5XrpcPyb0IotZxPnVK2kyeVV0tL\nRWG3baso65dfhpo1dTZssdB5PvGDBw9SqVIlxo4dq1Hif//9Nz179sTY2JiPP/4YgEWLFpVYEInk\nuWLs2KfKOw8Pk4MxB3l96+tceOtC7oAePVESrxBteTbN64PUB6hCVajCVNSoWIOgmCDa1WrHHu89\nOi9llpCgKOkTJ55umZnQrp2ytW2rKO/80qeXxlOLzvOJd+vWjejo6Bz7evd+Gv7aoUMHNm3aVKwB\nJZIXAl/fAg8vPrqY9zq9ZzAFDvpJLlUY2RcoGy1rhFqdxUiTlqwMsaCJ53imNKueM/GUlqSnK14i\nx47B8eOKwo6NVWbV7dsrFdOWLlXuq0UxiZRGSlxtKZFNfOXKlYwZMybPY/Pnz9e8d3d3x93dvSRD\nSSRlh+x2bgcH+P77YncxofUE+rj00YNw+aPr5FJF4YlLoK1xRb6/VI8Rmy9g7m4HE2ZAv36sNx9X\n7D6FULIPHD+uKO1jxxTrVf360LEjuLkpCQibNs3bpa8oGColblBQEEFBQSXrRBRCVFSUaNGiRa79\nX375pRg2bFie5xShW4mkfKFWC3H6tBDvvitEjRpCtGsnxNKlQty9W9qSFZnt24OFi8tsoahBZXNx\n+a/Yvj24RP1O2jpJuK1yE31VfUXC44QcxxIeJ4gRSzqLhO7thfj5ZyHi44vdf2qqEEePCvHdd0K8\n9poQtWoJYW8vxKBBQnz1lRD79gnx4EGJvkIu3Nzm5bhOTzY3t3m6HegZtNGdWt2n/vjjD3bs2MG+\nfftKdgeRSMoL9+8rC22jRinVccpwJGV+ttySJJcqiLC7YRy7cQx4mnjqCbYWtqx/5zC8U/T+7t6F\nw4eV7ej9Zq3PAAAgAElEQVRRZZbdqBF07gyDB8PXXysxUfrwFHmCwZ9aYmO1rm1abCW+a9cuvv32\nW4KDg7GwsNBqUImk3FG1Kly6pF/NoQMKs+VqG7347AKliZGJJm/JqZtKIeEmj6z41W1xsfoVAi5e\nfKq0Dx9WlHinTtClC3zxhWLTrlSp2CKXCIPkXLl7F9avB39/iIri0fkzWnVToHfKmDFjCA4OJi4u\nDgcHBz777DMWLlxIeno6VapUAaBTp0789FPOwqzSO0VS7niSt8TPDyZOhF69dD5EamYqFqb6nfjo\nywMle86SEc1GEJ0QRY1HRniHGdFtx3lm/McWn65fYztwBJiZ5dtPRobih33oEBw8qLxWrqwo7K5d\nldfmzcG4DFQ60FvOlS1b4JdflMeMAQOUuPxevbj68DouVVx0652yZs2aXPsmTpxYPIElkrKKEIpG\n8fNT8nM7OysugS+/rLMh7j66S0BYAL5nfWlfuz2/DPhFZ33nhb7Km5379wZUAevkmgyvMI5hp/7E\n7NJlQlt14M1m3bh/uhKjws8zw/xoDkWXkqIsPB48qGQVOHFCMYV06wajR8OPP0Lt2lqLplf0lXNF\nxMeDtzdGGzdCxYqa/fXt6mvVn4zYlLy4bNig1J/08tJp3pL0rHS2XtrK6rOrORhzkEGNB/Ft729x\nd3bXSf8FURJbbnaTyfd9vifwciDXIu4S+KUp8TdPwYDJJG334ZO/vqXy9/9BDDbNZbq5fPkz/vnH\njtRUV4KDlYRQrq7QvTu8+64y035O6lcUjFqtmEtq1NDsunL/Cv6h/qgeqVj76lraZFPgJUHmTpG8\nuGRmKokvdGznfpj+kBEbRjC6+WiGNxtOpQqGM+gWlJeksFll91XdOXjtIABmxqZ4V2hH/Y3JzDkV\nlquth8dchBB5mm7s7KJ4++16dO+u2LZ1pKvKB2Fhio17zRro0oV7K35gXfg6VKEqohKjGN1iNF6u\nXrSt1RajPH53Oo/Y1BapxCVlgowM2L1bmXH/+iu8IAvxhdlyn12ktLWw5fbD29T5vg6ZIpP6aVYc\n/MOEWm79mRVuwrdhqlxj1KlziMREZ5KTHXMdc3ObT1DQfH1+xbJFeroSKxAQAImJ8J//aPLP/nzy\nZw5fP4xXSy961e+FqXHBxg+dR2xKJOUOIZSkGH5+sG7d0zqUOpxUJKUlsen8JnxDfZnWdhojm4/U\nWd+6oDBbbl61KGtUdODin7X4b+fH+Lw8D9uL3mBtzele8/Pso0KFdBo33sWpU2/kOqbP4KEyiZmZ\noryXL1dWZ7Otyk5rN41p7abpdXipxCXPF2++Cfv2KYr72DEljE8HZKmz2Ht1L76hvgRGBOLm7Mb0\n9tPp37C/Tvo3BHce3mHNuTWohRpQCis8qUWJkREuIefwt6jMqVOwbyns3w9Hj87F3PwaaWlOmn6e\nmGegETNnlo3SZwYhJUUxwVkrpd+EEJy+dRpVmIpDLoc43u0rjI0M71YjzSmS54sHD5R/Mh3bubdc\n3MKXB79kbMuxjG4xmuoVq+u0/+JSWHKmJyYTc1NzRjYbyeaLmzkSc5hBNOKteqNYbHIcn4E+2Jjb\ncuGCkkVg714IDla8R3r2hB49FC+SgwfzN88YuvSZwXmSYiEgALZtg+XLierXGf8wf/zD/EnPSsfT\n1RNPV08aVyt5AJi0iUuef9LSIDBQydr/7rsGG1YIkedCVGlQlJSy2f26HUxs+C6iHkN2XKXiqwO5\nN3o6O+I7sHev8tBSoQL07q24xvfoocmW+2ITFQVLligmubp1FRv3qFFQowbD1g2jZuWaeLl60dGx\no05/F9ImLnk+UavhyBFQqWDjRqVwoY7jFVIzU9l2aRsB5wJYNXhVrqx6ZUWBQ7bkTAMnQ9UIyLAi\nclMAy5Z9p1HiTxJPtbprzN7LL3Gz5US+GD2U7UGViN2hKOvevWH+fMXiVIa+XtkgNVWJ0j10CBo2\nzHFo86jNpSRU3kglLinbqNWK0hZCsXOfPq3kE9UBQgiO3jiK71lfNpzfwEs1XmJsq7GYm5jrpH99\noQnoqRoBzspsmwGTSb3XDFAu1X8bBXAjchLO4Z/jfKopbdIUpb1ypZJDu6QlxZ4b7t3L8eiRpc4i\nOCYYVYSK6l2q8/UzCrwsIpW4pGxjbAw7dkCdOjqfLr6/5312XtnJ2JZjOTPlDHVs6gC6Lwag6/5M\nLVKh2UaodgEA8xQrKu75mATnK4wbp3hVWlvb4uGxgT7vg5+7EtpuSBnLNElJsHmz4s998iRcuECo\n8T1UoSoCwgKwr2iPp6sno1uMLm1Ji0YJsibmi566lTyvPHokxJo1QgSXLCVqcXmY9lCo1eoc+/JO\n1zpb63StJe3v2TSv3xz6RlT+orKo83o18VbbpqL/SFPxk8VoUc/ooujY8Z746SchIiMNK2O5Yf9+\nIUaOFMLaWslju26dEI8eiXuP7on6P9QXH+/9WJy7c65URdRGd0olLikdMjOF2LtXiHHjhLC1FcLD\nQ/kn0yHpmeli+6XtYuHBhUU+p0+fT/LMI+3hMUcrGUran9sqN8F8BPMRg/xGiBVvzhcXqliL3aav\nigkWP4mGtf8W7dr9Jv78M0Qr+XQhY7nB11eIX34RIi4u16Fnb+alhTa6U5pTJIbn3Dnw8FDySnh7\nw6JFOXJMlAQhBGdun8E31JeAsADq29VnQusJRT5f1wmkitLfsxGUaqHmYtxFOtbuTPojZYHSKqEd\nwR/44PDyQyw/mYr7KAdW6ihxVGmUbdMr9+/D/2dZBSWXzc7LO1FZbGVa12n0qFo11yllaeG6uEgl\nLjE8DRsqaV+bNdN51z18exCVEMXYVmM5OOEgjaoWL6mVrosBFKW/7BGUzZY350FKMu1jOhK+bQ9V\nawfQeMBkvhvoQ59vbTEz020x4aLKWOa5fl3JV+LvD1WqIPbv58j1I6jCVGwI30Bz++Z4uXrxck3d\nZagsM+j+gUCaUyRCiMREIVauVF4NyJX4KyJLnaX1+bouYVZYf5lZWaLmN06C+YgKs83F2x0aieuV\nbEV45zdEdGSm1t9DlzKWWbKyhPDxEcLNTYgqVYSYNEmIoCAhsrKE31k/0XR5U/FVyFciOiG6tCUt\nMtroThnsI9Ed6emKa4RKBbt2Kc7I33+vhADqCLVQExQdhLGRsd5Su+o6CjEwMISpge+QXOE+Zpjy\nv25LqWHbj23bYPt26ODQi/h6wSyN6onztIlUGD4QLC11+I2KJmO5jLz84AMlv22/fmD+1DU0S52F\nsZFxuTOTyIhNSenh6wvvvw9Nmij5uUeMyGGXLCkX7l3AL9QPVaiKqlZVmd11NiOaj9BZ//qmo09X\njt86DIBZxAheurKegQOVwi6tHh3BqFFDGSqZH2q1krckW422R+mP+OviX/x58U/8h/mzd/fx58JF\nUipxSelx5Yri062jhFNPiE2OZfDawcQmx+Lp6ol3S29cHVyLdG5p+D5nX6T87KXV/Bp4jK3Rftyr\nuhlMsnA2a8cuzz00rqt72/ZzhRBKRYknubnffZfMd2ey7+o+VGEqtl3aRhenLni5emEZbc8H7+4v\nMA1BeUEr3akjU04O9NStpLS5e1eIv/4y6JBZ6iyxN3KvyMwqnn24NHyft24NFpWnN9S4BDKngmjw\nTgOxeFAjcbWujRjxuatIeJygt/GfC+7dE+Krr4Ro3lyIunWFmD1biPBwIYQQnps8RfsV7cXSY0vF\nnYd3NKc8Ty6S2ujOAs+YMGGCsLe3Fy1atNDsi4+PF7169RINGzYUvXv3FgkJuX+UUok/Rzx6JMTa\ntUIMGCCEjY0QY8cqC0o6RK1Wi4MxB8Xdh3d11qeh/rHT0oTYs0eI/v1vChOTRIFnX8F8hN17FiKk\nhpm41bGbEJs2CfH4sU7HfW6JihJiyhQhDh7M9Tt7nJH3NXRzm5fn39rNbZ7+5dUx2ujOApPfTpgw\ngV27duXYt2jRInr37k1ERAQ9e/Zk0aJFxZv6S8oP77+vVLFduVKxcV+/DqtX66wU+ZX7V5gXNI8G\nyxowZfsUohKjdNIv6M/3OTAwhNpT21LxreaYv9ERu9q3mDsXIiLCyMqygU0BcG4EfX76niG3bzHe\nxg2GDXthqgplJzAwBA+PObi7z8fDYw6BgSFPDz5+nKtQR8LjBHzi97DgP3VyFVcAsDDN+xo+Fy6S\nJaEwLR8VFZVjJt64cWNx+/ZtIYQQt27dEo0bN9bJ3URSBtm7V4jYWJ13ezDmoOj0Wydh/629mLlz\npjh185TOI+Z0PROPjxfi3XcvCsuKZwXTXDUmk8pvOIvt24Ofq9mgLsjLnNWw/sfi6JffCzFhghKl\nGxYmUjNSxabzm8TQtUOF9UJrMWL9CLEjYkeJxyoXLpJ5oI3uLHawz507d3BwcADAwcGBO3fu5Nlu\n/vz5mvfu7u64u7trcYuR6J1r15QIt9atcx/r2VMvQ1a3qs7sbrPxcPHAzMRML2PMmNGHyEjtq85M\n3jaZc7cieBBnRfWQAE5euI14eSGp03ZhZpJMBtDilgk1VYtYdv1vzM3zXox6YWaDz6BJlws04QKv\n8ztjrq7h8ddqmP8hfPklmTXsafBDPRpUaYCXqxerBq/CxsKm2GM9WbxctmxuNhfJwgtDlwWCgoII\nCgoqWSeFaflnZ+K2trY5jtvZ2enkbiIxIAkJQqxY8TRI4vvvdT6EWq0WYXfCdN5vcdi+PVh4eMwR\nbm7zhIfHnCLNzGJilMthPeNpzpJO348Qa7b9T3j3ryZOVa0o/rAYKGq91lmYWtzRzLafp9mgLsj+\nZOKJn/iS2aIp4bmeTO6n3C8dAcso2ujOYs/EHRwcuH37NjVq1ODWrVvY29uX7C4iMRzx8TBlihLy\n3rs3vPMO9O2bI0iipEQnRqMKVeEX6gfAqUmnqGxeSB5UPVFYweAnREbCpk1KvYnIqEwGDzSlUUsr\nTj1Q6lDu8PbB9kYct6N20yven0TsYOPT8y0sssr1bFBnZGQoRYPJZqeudBv/FnFwpwdENcPpmScT\nO0s7Q0v53FFsJT5o0CBWr17NRx99xOrVqxkyZIg+5JLoA1tbJbpkxQqw0+0/z4bwDSw7sYwLcRcY\n1XwUvkN8aV+7vc4i5nTt833lCvxnzWQu3I0g7ZEl7pUnYz1+M5VTgvB5L4qH6QFM3jYZn4E+SpWf\nBrY0/OZjqs5cTGI+Jpqi3jSeK1JTlbBTlQouXYLz53mUkUJLTyMO1nXhcdV4uDgYYro/30WUS5EC\ng33GjBlDcHAwcXFxODg48PnnnzN48GBGjhzJtWvXcHZ2Zv369dja5i5lVUC3En1z/jw4OCjlpQzE\nmrA1WJlZ0bdhXyqYVNBp30WpKVkUrlyBDRuU7eZNEJPacc/sFAB2RpbMP+/A6KPJ2B8/l29WRUOG\np5fpQg1BQeDnB3/+CS+/rNSgHDaMkw8j6O3Xmy5OXWie2YbTa9PITLEsX6H8pYiM2HyRuXVLiWxT\nqeDOHeV9d93+wwghiH8cTzWrajrt9wn5KS0Pjzns2fNlrvYeHnPZteuLAvscEzCZE5ERxN2ywnx7\nAK8NsGXECOXS1PvGgevpd2lxz5iD9wdj6/k69OmjMQmUJrq6cemNadOgQQMYPVpxQ/1/0jLTSExN\nxKGSQykKV36REZsvIseOCdG7t+KyNWGC4haYqdvsd9EJ0WJByALRZHkTMSBggE77fkJBEZbFdd+7\ndk2IxYuFaNdOCLNJTxcoX1s3Ike7BL8VYsTXbUXCnRi9fKeSUGaiEPMI7IpJjBELDy4UD1IfGFaW\nFwBtdKfMJ17eqVwZ3ngDtmzRaea7LHUWf5z5A79QP87dPcdrzV7j90G/08mxk87GyE52l7QnREYu\nYNmyuUVy37t9WzGTrFsHFy7A4CFqRswKJjbmCjcfKguUKwb55Djf1usN1vOG7r+MDijVQg2JicpK\nr0oFrq6wdCmJqYlsPL8RVaiKsLthjGg2gpSMFKzNrfUvj6RApBIvDwgBYWFK1fdnadZML8UVjI2M\n+ff2v8zsMJN+Dfthblp0DxZtbLkFKa0PP+yRp8/3hAkDmDnzEr/cnEeGdSxmRmq8e71Jl0/PsSbM\nl9MXYep5OGFrh+9He5QFynKCwaMQ09OVgtQqleK91KsXzJgB/fqx7Pgy5hyYQ6/6vbT6PUj0i1Ti\nZZmICCWLm0oFFSoolbmzpePUBUIIMtQZuRYjjYyMWN5vebH7y8uWGxn5CUCBirwgpZXdfe/RIwuS\nkhpiYfE+EyfaYGR0kfQRt6HuQdIB3/TjTPK3Z/vhJFp2GARjvRWFZJrzp16mFw0pebBSsXnwAJYt\ngzFjcnkvDW4yGK+WXtIdsKyie6uOtImXGF9fxaBbo4YQ77wjxKlTQug4LD27nfvrQ1/rrF9tbbkF\nBcukpwuxfbsQLu9MEiavu4lqM/uKn1YmiFde+Vxp+/9Jp5jUTiy0eEN802KQEElJxRyr7FV31yZY\nqUjk8VuKiIsQvmd8ddO/RGu00Z1yJl4WsbCAL79UKuOY6u5PlJyWzPrw9Ro794jmI3Ru59bWlvts\nsIy5uZqePUcSGNia8eOhUSMw6hdBlm0wcYBK3Q8zdR/l5E0BMGAybPfhv6m2uFWdz4eV8w8wKsj+\nXpZm4zr1O8/uvbRwIXh4cO/RPdaFr0MVqiI6MRrvVt4IIcpdNZwXHanES4uMDMVZ2dk597ER+qlY\nE5scS+DlQJ3YNfMzR5TEltu/f3fq1euOSgUBAXDjhuJ+vG5vBAfu+7H48HEAHFPM+D2qJTOfjJVq\nCxvXF3ksQy8aqoWa4OhgnG2dqWenu1J1hZKcrPhxq1SKKW7IEPj2W3B3Z9xf49hycQv9G/Vnnts8\nerv0xtRYqoNyie4fCKQ5JV/UaiGOHhXi7beFqF5diPHj9TSMWudZAbNTkDlCmxwisbFKvpKXXhLC\natQkUWeum+jyc19xPyVB7L0QKBy+sBHvTnEWwU0rihEfOImEnX8KkZmpdb4SQ7nvhd0JEx/9/ZFw\n/N5RtPq5lQiKCtJp/4WiUgkxcKAQ69YJkZKS49CRa0dEUmr+JidJ6aCN7pRK3BCkpQnx6adCuLgI\n0aiREJ9/LsSVKzofJiohSnwR/IVovKyxOHPrjM77f0JhSrAottzkZCH8/ITo00dxcR8/XnFxd1v1\n1K97xPoRIiP6qsjw6K0opIcPc/Wjjd1Y38mqjl4/Klr/0lrU/q62mPX3LBF6O1Qn/RaXc3fOieM3\njpfK2BLt0EZ3yucnQ2Bmpti2166FNm1AhzbHhMcJbDi/AVWoigtxFxjZfCSrBq+ipUMe7og6ojBz\nRH623KwsJVp72o7JRCVFYFfJioVj/flo+QU2XfanTecvsdpsBSh+3T4DfTC1sIVde/KVRRu7sb6T\nVTnZOPFdn+9wq+uGibEe/bqvXlW8l7ZuhZAQsLTkVvIt1pxbg1+oH/ce3WOe2zza126vPxkkpY5U\n4rrk4UPF1v1scikjI5g7Vy9Drj67mkPXDvF+p/f1krckL4pr9z5/XkmzoVKBuXkysX2OkekYxj3g\n7chaOEda433OGCwHEDD8mcRTeqKki4aZ6kyCo4PpUa9HroXAWpVrUatyrQLPf5zxmJ1XdjKs6bDi\nDRwXB+vXKxfzyhUYNQqWL+dO5gPGqYZy/OZxhjQZYpibiKRsoIcnghfLnPLE/23MGCGsrRUbwXNO\nUcwR9+4JsXSpEG3bClGzphAffijEjz+eFC4us4XR2y6C+QjbD8zE1vpmIqbXq0Ls36/z2p26Rq1W\ni39i/xHv7npX1FhcQ7Rf0V7ruqD3U+4Lx+8dxd7IvcU70dNT+a1t36789v6fzKxMsf7cevEo/ZFW\n8kjKBtroTpkAS1uiomDxYiXWu2FDxY1i5EioptvkUJH3I1GFqjh47SB7vPdgbKSb+pZFoaCAmLyy\n+fXp052dO5UynNuYjF2DCOrWsmLnGwFUrWirSWQ1ym4B13stx2Hbl+xKHUN3j4WFJrIqbfzO+rHo\n8CIeZzzGq6UXXi29aFS1UYn63HVlF1O3TyV0WmiRwteFEJy9fYY6Nk5UtTJchkqJ4dBGd0pzirZk\nZkKtWnD8ONTTrdtYfEo868PXowpTcTn+MqNbjGZhz4UYYTj/3cIiL7ObI86eVRT3xInK/WzsWEFk\n1j+cvXuau4mKDXz9iPUaW/q6hE9gwyeafg2SD6SE1Laujc8AHzrX6awzP+pXG7xKb5fevL/nfVYM\nXKHsPHdOMZWkpcGSJQDcSLqBf6g/fqF+PEx/SMDwADpbddaJDJLyj5yJF8adO2Bvr9PFyMLo49eH\nKpZV8G7pTR+XPnqrQ1kQhaV/vX8fBvpMJvRmBJmPrXirRgD9hydyOFmF3+lV3HhwnRQyaFezLXvG\n/o2thW2JUsoagoysDK7cv0LT6k0NNmZSWhItlzfn58xX6bvmpFJ9ydMTvLw4YpPEnP1zOHvnLMOb\nDse7pTddnLoY9GlMYljkTFxXJCXB5s3KjOiff5TkU46OBht+t9fuUo+ay88DJTa2MaNGwe7dYD4l\ngofVggFYW6k5f2xPYtQNG1YfSqZx19FMaXUdH69NmgVKg+cDKQJCCE7cPIEqTMW6c+voXrc7G0du\nLPxEHWGNOSvXpbK12xH6LlkObm5grChpm7vhvNnuTQY0GoCFqYXBZJKUL+RMPDsHDsBPP8GePfDK\nK8qMaMAAnaZ4BbgYdxFVqIqKZhX5b7f/6rTv/Chuwqf8Zs3W1jdYtMiR0aPBc0c/dl7ZSbta7fj+\ngDntUmww9xqvXDOLvJWOISvjFIQQggUHF+AX6odaqPFu6Y1XSy/q29U3uCxkZYFJ2TcpSfSPrOxT\nUtasUWbhI0ZAlSo67frOwzusPbcWVZiKm0k3+Y/rfxjXahyuDq46HScvtKkSs3nzIaZNi+Vuh71Q\nNQIyLKl2yYWPPmjGB6PeBCAxNfGpO2AFa80Msryw5OgSOtfprNNaoLkQAk6dUnwse/WCQYNyHC7r\n2RQlhkUq8aISH2/Q+pP3H9+n4bKGDGg0AC9XL3rU62FQ/93CbNHZFUlaWnWqVh3K0aO1cHa+z8Wu\nHUmpchkAS7UpPxkPYPy8Pw0me0lJz0onOS3Z8N4cUVFP0whnZYGXF4wfD3XrapqU+RJsEoNjUJv4\nwoULUalUGBsb4+rqyqpVqzA3L8OJ4q9dU2ba/v5K1feQEIMNXcWyCrHvxZZaIv2CIiwDA0N4++0g\nol3vQq0IyLDCZm9tPvvanG3GX5NxPRKAJveNOZLhhd34aYYUXSuEEBy9cRS/UD82hG9gVpdZzOoy\ny3ACBAUpT3OjRsEff0CHDnkujJeXbIqSso1WSjw6OpoVK1Zw4cIFzM3NGTVqFGvXrmXcuHG6lq9k\nCKEkuFepIDwcXnsNli+Hrl11PIwg9E4oqjAVw5oMo1Od3KlddaXAtXn8zi/C8vHj2kyeXInY2E/B\n3R2clUXKB2mV2LvOmbeS/kVV91VmdE7E590/sbW218l30BfxKfEsO7EMVagKU2NTvFt6c2ryKZxt\nnQ0rSNeuEBtbaMHlUi3BJnlu0EqJW1tbY2ZmRkpKCiYmJqSkpFA7W8XrMoORkfJY+9570Lcv6PhJ\n4Yn/ripMRVJaEp6unoWGW5cEbavm5OUVYmZ2h5s3PTGq6QuJTpCh5CzhZjvY7kNyh/8x9MAtsLJi\nfT79ljXUQs39x/dZM3wNbWu11a+d+/Bh5alu4ULlyS47RcwBb/ASbJLnEq2UeJUqVXj//fdxcnLC\n0tISDw8PevXqlaPN/PnzNe/d3d1xd3cviZwFk5UFjx6BdR5RbwsX6mXIvy7+xetbX2d40+H82O9H\nujp11bv/rraP3/36defcORs+Oz2cx1b3MDNNxc2pLXetD3DzRhSt72VwJlthBVJtFUViZaXX76Mt\naZlpmBqb5lpXqF6xOkv7LtXfwJcuKU91/v6Kx5K3d4m6K4sulxLDEhQURFBQUMk60Sa+/8qVK6Jp\n06YiLi5OZGRkiCFDhgiVSqU5rmW3xUOtFuLECSFmzlTKmH3zjf7HzEZKeop4nPHYoGO6uc3LMwWs\nm9u8PNvHxyt5ups0UTbHz1/WpHl1nmUuDrSyEZEDhoredSbrLS2rrlCr1eLwtcNi6vapourXVcXR\n60cNK8Cnnyq/s3ffFeL0aZ2Vy9NbCTZJuUQb3anVTPzUqVN07tyZqv/v4TFs2DCOHDmCp6dnye4o\nReHuXfj5Z2U2pFYrvtxBQdC4sU6HEUJw7MYxNpzfwNe9vs4VNWlpVjLfcV3atp88fgcGhjBl+0we\nmD4mPaUapn/9xdC+1fDxUcy0bv/L4EYStHlsx96OP2L75WtgZsbMwBCM9ZSWtaTEJMaw6syqHHbu\nfyb/Q13buoWfrEtmzlQyUeqwXB7ouASb5MVEm7vFmTNnRPPmzUVKSopQq9Vi7NixYvny5SW6mxSZ\n69eFmD5diGPHdF48WAilYOynBz4V9X+oLxovayy+DP5SJKcl63QMbQv1FpQ9cP36Q6JatS2C8U+L\nKliNa5Gjz4SEW2KE/xCR8DhBp99Hn2wI3yCm75guTtw4oddqRSIrS4h9+4RYvFh/Y0gkhaCN7tTa\nT/ybb75h9erVGBsb8/LLL/Pbb79h9v+r8TrxE09KgkqVDBpA8s6ud1h7bi2jW4zGq6UXbWq20cvi\nWElyiDwb8div3xDOnWvDyvXRZDX9CyP3OQiLRxjFtkb4HsDD7bsykZekMERpFeg9d04JxAkIUDJQ\njhsH77xjeDkkEgzsJz5r1ixmzdKx7216OuzcqZhKdu+GQ4fAVf8RjU/4uOvHLO6zWO8FY0viWta/\nf3c2pas4GhHBrWtWnNrkir3bx1R4Mxj3CzaM3GzK/1o24Nr2jSSk2pZpdzXx/yYrv1A/9kTu4fxb\n5w1S1EJD795w4YJiktu5E1q0MNzYEomOKBsJsE6fhl9/hU2boFkz5Z/q5591HlWZpc7iQPQBYhJj\nePtMIEMAABg/SURBVP3l13Mdr1Gphk7Hyw9tXcsiI+GXX0D1IIKM2sFgD+16hDFzZzy2Nxuw/MZ3\nvEFvsiKe/lnLorva1YSrqEJV+IX6YWxkzNiWY9k7dq9hFTgoeXLq15d5SyTlmrKhxK9eBRcXRZk7\nOem0ayEEZ++cRRWqIiAsgFqVazG5zWSdjlFciuJaNnnbZCLiI7A0tWJcpQBW/mbCv8crM2ECdG5n\nRXCsUodyj/vv2L7nTGDIv1yeuZusyL759llWWHBwAZamlqiGqvSbtyQrC/bvVxYjX3kl9/GGDfUz\nrkRiQAybO+XBA7Cx0fVw+ZKlzqLDbx2IS4nDq6UXnq6eBs0VXRCFZfPrssKdI7FKBKVZmjUtqcPB\nueewtHwm8VS2OpRlJUNgqXP2rOLPHRAANWvCf/8Lw4eXtlQSSaGUzQRY9+4pJcz8/SEhQQl/N+AC\n1oV7F2hcrXG5SaR/8pSa1xev5FydjxGV4rF5bMzP+60Y0Wc6pl9+VdriFYoQglOxp/AN9UUIwfJ+\nyw03+I0b0K+fMlnw9FSCcZqWjZu2RFIUypYSX7NGWfU/fFj5x/L0VBaSKujW7pmWmcaOyzuoVbkW\nHRw76LRvffPEZGJhYsUwEcCqn225EvmYll3r0ffGA9bWq07lrf/jZq2TLFnat0zPqmMSY/AP88f3\nrC+Z6kxNfm6XKi6GE0KtVhbDu3Ytd2lxJRIoa0p8zhxo0gQGD1ZcBXWIWqg5dO0Q/mH+bDy/EVd7\nV+Z2n0vP+j11Oo4uKCiop5OPO8duKSYT+3sjWNFnPcuXz+XR369yllY84ul1KyslzPIiNTOVxssb\n069hP8a2HEtHx476s3NnZsK+ffDSS0rZPInkOaJslWf7Qj8KJ/ROKAPXDMTa3Brvlt6cmXKGOjZ1\n9DJWSdEkrGp2T5Pm9cJ/79M2dBfH43eQkJEAVaCFXTsOfuSDrQV8/70JR+iSq6+y7CpoYWpB1Mwo\n/ZmshFDs3H5+SjphR0fFTUcqcYmkjHinFIMGVRqwbcw2Wjq0LG1RCkWTsKqbGzgr+ctv1jem2enK\n/HDFlPavf8oHDodyLFCWxcx2Tzx8fM/60qNeDwY0GpCrjd4U+L59SvBNcrJSWGH/fuUJTyKRAGVU\niSelJbH5wmZGNBtBxQoVcxyzMrMqFwocngb1mJFGBuCUYMTbq4Zx07kir53+HUxMWM+MHOeUpcx2\nscmx+If64xvqS3JaMl4tvWhevblhhXBygh9/lHZuiSQfyowST89KZ/eV3fiH+bPzyk7cnd3pUa9H\nLiVeVnmySGllZkXA8AAiw22JiBgMgFi7lToDBnN3+yZmpdbCo+bcfANMntjLl5VyQqqDMQcZtHYQ\nw5oOY3nf5XSr201/s+3MTDh5EjrlLqZBw4bSn1siKYAyUWNz5b8r+WjvRzSu2hhPV09GNh9p+JqI\nJSAxNZEOv3UgIj4CgKq3R2AVuJ5eva5y4EAA0dFzNG1dXGbzww9lJ0tgfmRkZZChzsDKTE85xZ+1\nc9epAwcOlNkc5hKJIShb3inF6PZS3CXMTc0NUkarJNXFs8+2Vw9ZzbEbx/A9/it7ovZh/iiDexWz\nqHa/Cd92OorXa7aYmpbtAJwL9y7gF+rHrC6zcgQN6Z2VK2HJkqd2bi8vaeeWSChr3inPcCv5Fidu\nnmBwk8G5jjWupttc4PlRWHmzwhR8RHwEwTGKS2DvXzpjdf0eo49l8Er4MHZU6sjhgX/geuMVqvcJ\nxdS0u6bfsqK0Ae49usfac2vxDfXlZtJNPFt6kpGVYVgh7OyUWqfdukk7t0RSQvSqxJPTktl8YTP+\nYf6cjD3JyOYjGdR4UOmkHKXg8mZAofUrn5gWmlRuR6edi0k7dodjvVtx5PafxMRMhxXTOQBciyi8\n7mVpsOToEj4L/owBjQawoMcCetbrmavEmc7IylKKBdfJw/1z6FD9jCmRvIhokbe8UAAxedtkYbPQ\nRgxaM0isO7dOpKSn6GOoYlFQebM+fT5RPg+cpBRW8OotaLdM1JvaQGSps0R6uhAr/BKE3aQRon6z\nBPHjj0I8fCienvfM5uExp7S/bi5uJt0USalJ+h0kNFSIDz4QolYtITw99TuWRPKcoY1K1ttMfFTz\nUSzosYBqVtX0NUSxKcgHOzXVFIwzwfEYOIQBYG9/gI/2WvC/RY/5388VqVfPlj/eX8+AX55aAUqS\nG1wfXE24ypHrR/Bq6ZXrWK3KtfQzaGYm/PCDskh5/75i4967V+YtkUgMgN6UeI96PfTVtdZofLCb\n3YOqSgSl8z9NmD59CEuX7qHm6A4kWZ/nEeAYW5VaviqmZ3ZjeGhFNm+Gtm1z91kWgnMepD5gw/kN\n+J715ULcBTxdPfF09TSc2crEREl0tmQJuLlJO7dEYkDKjJ+4IXhio/beN5YEmxgAHNs/pn//7wHY\nN+kUEQmr2DFgCze2/0ayeRi//nSOCRPyT6xV2sE5bwa+iX+YP73q9+L9Tu/Tt2Ff/RVXUKshNTW3\nG6CRESxapJ8xJRJJgbxQSvxh+kPia0dRq54lCfeVogrbvP9EiP9r7+6joqz2BY5/x0AtsCt5AMmx\nQ06arwwEaqdjihrhxa4vQWoqoMnxZm+aZt08ecTq+O7xElnLOtpiPOrVfFeUTGsM7Zh6RA20sAKP\nEhBLUV4URdj3jyfJFwQZ5plh5PdZy7Vi5pm9N896+s3DfvZv/+Cuu3rzhU83vi9swoOHfPB/7H0m\nT+7NwIE174zo7OScEV1H8E7fd/RdV3/smDZV8o9/wOuvw8sv69eXEKJObF4nfu7cOeLi4sjIyMBg\nMLBs2TIeffRRrVF7FEqup98q49zN+ODxrD9oYctPKTx++i5GqW6sD2/L++EfsX1DSxYs0GYAXnsN\nhg+3+2659ZZbnEvBhQLHbjdQWKgFbosFcnNh5EiIiXFozVMhGhuHrhOfOHEiERERrF27litXrlBa\nWmprUzYn4NT0uWvXdH9zfCfx+5qxsNMIfMb9N0WdenLqIwOPdNaevS1cqG117qSVj9W6UH6BTd9t\nwnLUwr7T+3iz15uODeK//KKlws+eDf36SR1KIRoom4L4+fPnSU1NJSkpSWvEzY3/sLHsWm0JOLV+\n7pptXjNfu1j1uatrugMu38fu4ERazogi90xT/icBPh4I4eGwZYu2LXVDcr7sPFN2TGH98fX0aNOD\n6IBo1g1bp2/6O9z8Dfbww9qduBCiQbMpiGdlZeHt7c3YsWM5cuQIwcHBJCQkcM81D7weeqg3PXua\naN/+94SGhhIaGlptWzUl4NSURVn1ud694fepAGRf8iQx8XMGDuzNysiVVXUo80+2ZOqLsG6dVmDo\n4EF48EFbfnP9eTb1xOxr5u2+b+u3JBDgxAktSC9fDps3yzSJEE5gtVqxWq31a8SWBekHDhxQbm5u\nav/+/UoppSZOnKimT59+3YJ1UMpkmqa2bt1dY1s1JeBs3bpbmUzTrnvdZJqmNm/5UgUMjlEtBj+t\nmkxvoohHtXyxtaJ5oerTZ0ZV2/v2KTV0qFLe3krNmKFUQYEtv60+Ci8WqvNl5x3b6ZkzSn3wgVKP\nPqqUj49SkyYpdeiQUpWVjh2HEKJatoRkmxb0Go1GjEYj3bt3ByAqKopDhw7ddJx2R/15jW3VtM66\n6m77v8bDmFAYFUFuZXdiv/pPLrdZxZRfdvD6ByPxSu/PuaXHoKwlzZpVsGOHNo07fDj07QtZWRAf\nD7+rR95RcvJXhIe/RWhoPOHhb5Gc/FWd2yivKCc5M5lhnw7D/3/9sWZbbR+QLZYvB6sV3npLKyq8\naJE2n9SQHgYIIerEpumU1q1b07ZtWzIzM+nQoQM7d+6kS5fqiwXUlrlY0zrr+fO/0F5olQn+2kNK\ntysXmbvFiHl4HDH5hXx/Zg6s1Q7z9V3Fjz++zuTJ8MYbMGIEuLvb8htez9Z5+6uyCrNI3J/Iym9X\nYrrPRExADEueWoLX3V71H1xdTJyo/RNC3DFsXp2SmJjIqFGjuHz5MiaTiU8++aTa42rLXLwaBJ9P\nfoTipmdxx40ZfeZz0T8f9+aXtIPKf51rz+lO0aYNrOuzkD/Fv8HC7l+RkBBPdnYwOTmP4eUVwbx5\nLYCveP/9HSxdWvftZqtT27x9bc6VncOjqQepY1Np30rHAgcnT2pruXfvhpQUyZwUohGwOYibzWYO\nHDhQ4zG3m7k4cGBvTAX3svtkGgDjjkTRJ6Upfxq6gKwf/syP61bCU+Nh60eY2szj5ZcHUFIC33/f\nm4yM3gQEwN//ru1sum1b/e6aq3O7+6NcrrhcbbZkkF8QQX46LYMpKtKe2Fos8O23MGwYvP22TJEI\n0UjolrEZHj79tjMXl+9ZzDcnvwbAWNyEz9RoOj8/CQID8TCmkpi4kLKCzjTvs5AxY57i4MHHGDNG\nm+/euvX6ZYL1vWuuTk3z9kop9vx7D5ajFtYfX0/6hHT8WvjZ1I9NoqLAwwNeeQUiIqBZM8f1LYRw\nOt2CeErKOze9dmMdyqvVZDrtPsaOf/ch0XyJj6asp6Xnb08grxZVyM2Fv/0NXngBhgyBPXu0pcw3\n0mNXwerm7R8IeIFWkRcxvWeiuVtzYs2xHH3+qGMDOMD27ZKII0Qj5rC9Uy5ducSh3EP8K/dfgBbQ\n1zyzBoCQPy8G4PFqPnfyJMybp5VhjI6Gw4e1Aui3oseugtXtj+IXVYlnG0/WmNcQ7Bes346BeXmw\ncqU2PfLqqze/LwFciEZNtyAesSKCFU+v4FjBMZbvfo/VmRsxXLoCHnBvsR+RTWNr/PyJE1rG96ZN\nMH48fPcd+PjU3q9euwo6tMzaxYvaL26xwD//qf3pERfnmL6FEC5Ft0LJxIMHTWlb3ISRR5vg/t3D\nWAoSOf5Uwq8PKOeTkBB+U2BMT4dZs+Dzz+Gll7QN8+67r27926M4ccYvGSQdSeJo/lFSRqfUbQD1\nUVwM7dpBcLD2p8eQIdqctxDijtegqt37v9GMZYV9CH1mKhHzviDl81k3HRcePr1q7vzQIXj3Xfj6\na5g8GSZMgBYt7D2ymhWUFrAqfRWWIxbySvIYHTCaWHMsnbwdXKHm7Nm6f3MJIVxeg6p2nzYjl5a/\nJrNcfHdPtceUld3FN9/AO+9oc91Tp2rLnG+sOeAog/5vECYvE7P7z6bfg/30KyJ89iysWQPdu2t3\n3DeSAC6EuE26BfGW12Qj3uphY3r6KIYN07Ir166F5s31Gs3t2fvcXpoYdEqQKS/XEnCSkrS5ogED\noGfNBSeEEKI2Dknpe+WVJzGZ/nzda25uZxk5UnHihLZs0BEBPKcoh7l75rLk4JJq39ctgO/dC23a\naMtswsO1JTerVze8fXCFEC5Htznxa5tVCmbPPsKCBfdSVnYvJpOVd9/1YfDg6hYV2teF8gtsOL4B\ny1ELB3IOENk5kvGPjKd7m+66912lqEgrJGwyOa5PIYTLaVAPNpVSKAW7dsHMmVqhmLfegmefBTcH\nrU4vKC3g4fcfpkebHsSaYxnScQh3u9+tT2elpdqywGeesc+uW0KIRqdBBfHPPlPMnAlnzsD06dqO\ngs7ISykoLcDbw1ufxisrtc2mLBbYuBH+8AdYuhT8HJy1KYS4IzSoIN6xo+Ivf9H2Y9IzeBddKuLT\njE95rO1jjl0KaLFo305eXhAbq/2J0bq14/oXQtxxGtQSw/R0/YJ3RWUFu7J2kXQkieTMZPo+2JeQ\n+0P06exWunTRypqZzY7tVwghruGQB5v2ZM22Mnr9aFp7tibWHMuz3Z7ld/fUo2RPTcrLISMDAgP1\naV8IIa7RoKZT9AriBaUF5Jfm09Wnqy7to5SWPmqxaLtuhYRAcrLszy2E0J0tsbNBln4prygn5YcU\nKlXlTe95e3jrF8ATEqBrV22FiZeXtvnUtm0SwIUQDZbDtqK9HUfyjpB0JImV366knVc7HvF7BB+P\n29i60F6aNYMlS+CPf5TALYRwCQ1iOmVNxhpmpc7i7MWzxJhjiDHH0KFVB3sPS1NZCefPa3faQgjR\ngDh0TryiooKQkBCMRiNbtmyp10B2/bQLg8FAqH+ofqnvP/ygzXMvX65t77pokT79CCGEjRy6xDAh\nIYHOnTtTXFx8W8crpcgtyeX+Fvff9F7/dv1tHUbNLl7UgnZSkhbER46E9etltYkQ4o5h023v6dOn\n2bZtG3FxcbV+a+QW5zJ/73y6fdiNqDVRNg3SZkppGZVvvgmnT2t330FBMt8thLhj2HQn/uqrrzJ/\n/nyKiopuecywF4eRlpfG6aLT9Avtx4ejP6TXA71sHmitKiuhyQ3fSffcAytW6NenEELUg9VqxWq1\n1quNOgfxrVu34uPjQ1BQUI2dm542Mdh7MEM7DeUed52qPOTna0WEk5K0u+3hw/XpRwghdBAaGkpo\naGjVzzNnzqxzG3V+sDlt2jSWL1+Om5sbZWVlFBUVERkZicVi+a1RHZN9uHQJtmzRAndqKgweDDEx\n0LfvzXfiQgjhQhyesbl7924WLFhQ79UpdZKcDAsXaptORUaCp6c+/QghhIM5ZQMsg6MfEg4cqP0T\nQgjRMJJ9rlNSoi0DXLFC27tEigYLIRoJ1907pbISvvwSxoyBtm21SvBxceDh4eyRCSFEg9Yw7sSn\nTIGdO7V57pEjpbiCEKJRct2taMvKHFPuXgghGrCGO51y5Qps3w5//Wv170sAF0IIm+gbxNPTYepU\neOABiI+HVq107U4IIRob/fYTj4iAo0chOhq++AI6dtStKyGEaKz0mxPPzIR27fQtdS+EEHcQ132w\nKYQQogE/2BRCCKELCeJCCOHCJIgLIYQLkyAuhBAuTIK4EEK4MAniQgjhwiSICyGEC5MgLoQQLkyC\nuBBCuDAJ4kII4cIkiLsAq9Xq7CHcUeR82o+cS+ezKYifOnWKvn370qVLF7p27cp7771n73GJa8j/\nKPYl59N+5Fw6n01b0bq7u7No0SICAwMpKSkhODiYsLAwOnXqZO/xCSGEqIFNd+KtW7cmMDAQAE9P\nTzp16sTPP/9s14EJIYSoXb23os3OzqZPnz5kZGTg6empNWow2GVwQgjR2NQ1JNersk9JSQlRUVEk\nJCRUBXBbBiGEEMI2Nq9OKS8vJzIyktGjRzNkyBB7jkkIIcRtsmk6RSlFbGwsrVq1YtGiRXqMSwgh\nxG2wKYjv2bOH3r17ExAQUDX/PXv2bAYMGGD3AQohhLg1m6ZTevXqRWVlJYcPHyYtLY20tLSqAJ6S\nkkLHjh1p3749c+fOtetgGyN/f38CAgIICgqiR48ezh6Oy3nuuefw9fWlW7duVa+dPXuWsLAwOnTo\nwJNPPsm5c+ecOELXUd25jI+Px2g0EhQURFBQECkpKU4coWu5Vb5NXa9Pu2ZsVlRU8NJLL5GSksKx\nY8dYtWoVx48ft2cXjY7BYMBqtZKWlsb+/fudPRyXM3bs2JsCy5w5cwgLCyMzM5P+/fszZ84cJ43O\ntVR3Lg0GA5MnT77pZk7U7mq+TUZGBvv27WPx4sUcP368ztenXYP4/v37eeihh/D398fd3Z0RI0aw\nadMme3bRKMlqH9s9/vjjeHl5Xffa5s2biY2NBSA2NpaNGzc6Y2gup7pzCXJ92qq6fJucnJw6X592\nDeI5OTm0bdu26mej0UhOTo49u2h0DAYDTzzxBCEhIXz88cfOHs4dIT8/H19fXwB8fX3Jz8938ohc\nW2JiImazmXHjxsnUlI2ys7NJS0ujZ8+edb4+7RrEJcnH/vbu3UtaWhrbt29n8eLFpKamOntIdxSD\nwSDXbT1MmDCBrKwsDh8+jJ+fH1OmTHH2kFxOSUkJkZGRJCQk0KJFi+veu53r065BvE2bNpw6darq\n51OnTmE0Gu3ZRaPj5+cHgLe3N0OHDpV5cTvw9fUlLy8PgNzcXHx8fJw8Itfl4+NTFWji4uLk+qyj\nq/k20dHRVfk2db0+7RrEQ0JCOHHiBNnZ2Vy+fJnVq1czaNAge3bRqFy4cIHi4mIASktL2bFjx3Ur\nA4RtBg0aRFJSEgBJSUmSrFYPubm5Vf+9YcMGuT7rQCnFuHHj6Ny5M5MmTap6vc7Xp7Kzbdu2qQ4d\nOiiTyaRmzZpl7+YblZ9++kmZzWZlNptVly5d5HzaYMSIEcrPz0+5u7sro9Goli1bps6cOaP69++v\n2rdvr8LCwlRhYaGzh+kSbjyXS5cuVdHR0apbt24qICBADR48WOXl5Tl7mC4jNTVVGQwGZTabVWBg\noAoMDFTbt2+v8/VZ7w2whBBCOI9U9hFCCBcmQVwIIVyYBHEhhHBhEsSFEMKFSRAXQggXJkFcCCFc\n2P8D5WiaMzPf6nEAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Feasible Weighted Least Squares (2-stage FWLS)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "resid1 = res_ols.resid[w==1.]\n",
      "var1 = resid1.var(ddof=int(res_ols.df_model)+1)\n",
      "resid2 = res_ols.resid[w!=1.]\n",
      "var2 = resid2.var(ddof=int(res_ols.df_model)+1)\n",
      "w_est = w.copy()\n",
      "w_est[w!=1.] = np.sqrt(var2) / np.sqrt(var1)\n",
      "res_fwls = sm.WLS(y, X, 1./w_est).fit()\n",
      "print res_fwls.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            WLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.876\n",
        "Model:                            WLS   Adj. R-squared:                  0.873\n",
        "Method:                 Least Squares   F-statistic:                     338.3\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           2.25e-23\n",
        "Time:                        20:53:32   Log-Likelihood:                -43.178\n",
        "No. Observations:                  50   AIC:                             90.36\n",
        "Df Residuals:                      48   BIC:                             94.18\n",
        "Df Model:                           1                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1             0.4390      0.019     22.520      0.000         0.400     0.478\n",
        "const          5.2710      0.177     29.828      0.000         4.916     5.626\n",
        "==============================================================================\n",
        "Omnibus:                        4.076   Durbin-Watson:                   2.251\n",
        "Prob(Omnibus):                  0.130   Jarque-Bera (JB):                4.336\n",
        "Skew:                           0.003   Prob(JB):                        0.114\n",
        "Kurtosis:                       4.443   Cond. No.                         16.5\n",
        "==============================================================================\n"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}

Hi, I was curious how you determined the weight for the WLS. I've seen somewhere w=1/(x^2), but what if the x data contain 0. Thanks.

I have not looked at this since I posted 8 years ago. I'm unlikely to find an answer for you. Sorry.

No worries, thanks anyway.