{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weighted Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WLS Estimation\n",
    "\n",
    "### Artificial data: Heteroscedasticity 2 groups \n",
    "\n",
    "Model assumptions:\n",
    "\n",
    " * Misspecification: 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",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, (x - 5)**2))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5., 0.5, -0.01]\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,1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WLS knowing the true variance ratio of heteroscedasticity\n",
    "\n",
    "In this example, `w` is the standard deviation of the error.  `WLS` requires that the weights are proportional to the inverse of the error variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.927\n",
      "Model:                            WLS   Adj. R-squared:                  0.926\n",
      "Method:                 Least Squares   F-statistic:                     613.2\n",
      "Date:                Tue, 17 Dec 2019   Prob (F-statistic):           5.44e-29\n",
      "Time:                        23:40:04   Log-Likelihood:                -51.136\n",
      "No. Observations:                  50   AIC:                             106.3\n",
      "Df Residuals:                      48   BIC:                             110.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2469      0.143     36.790      0.000       4.960       5.534\n",
      "x1             0.4466      0.018     24.764      0.000       0.410       0.483\n",
      "==============================================================================\n",
      "Omnibus:                        0.407   Durbin-Watson:                   2.317\n",
      "Prob(Omnibus):                  0.816   Jarque-Bera (JB):                0.103\n",
      "Skew:                          -0.104   Prob(JB):                        0.950\n",
      "Kurtosis:                       3.075   Cond. No.                         14.6\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "mod_wls = sm.WLS(y, X, weights=1./(w ** 2))\n",
    "res_wls = mod_wls.fit()\n",
    "print(res_wls.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS vs. WLS\n",
    "\n",
    "Estimate an OLS model for comparison: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.24256099 0.43486879]\n",
      "[5.24685499 0.44658241]\n"
     ]
    }
   ],
   "source": [
    "res_ols = sm.OLS(y, X).fit()\n",
    "print(res_ols.params)\n",
    "print(res_wls.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the WLS standard errors to  heteroscedasticity corrected OLS standard errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================\n",
      "          x1   const \n",
      "---------------------\n",
      "WLS     0.1426  0.018\n",
      "OLS     0.2707 0.0233\n",
      "OLS_HC0  0.194 0.0281\n",
      "OLS_HC1  0.198 0.0287\n",
      "OLS_HC3 0.2003  0.029\n",
      "OLS_HC3  0.207   0.03\n",
      "---------------------\n"
     ]
    }
   ],
   "source": [
    "se = np.vstack([[res_wls.bse], [res_ols.bse], [res_ols.HC0_se], \n",
    "                [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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate OLS prediction interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "prstd_ols, iv_l_ols, iv_u_ols = wls_prediction_std(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare predicted values in WLS and OLS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFnCAYAAAB+YZr1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1xV9f/A8de5l+lAxC2498ZEE81EcWXu1FJzlZnr6x717furbLhwZZlpZrkqt+VKza3gwE3ugQg4EEVF5r33/P74CGY5EC7ci76fjweP8Nxzz/lwNd7ns95vTdd1hBBCCJG1DLZugBBCCPEikgAshBBC2IAEYCGEEMIGJAALIYQQNiABWAghhLABCcBCCCGEDTw1AGuaNk/TtOuapoX87Zi3pml7NU07omlasKZpdTK3mUIIIcTzJS094J+AFv84NgkYq+u6N/Dx/T8LIYQQIo0cnnaCrus7NU0r+c/DgNv97/MAkWm5Wf78+fWSJf95KSGEEOL5dPDgwRu6rhd41GtPDcCPMRTYqGnaZFQvut7jTtQ0rS/QF6B48eIEBwen85ZCCCFE9qJp2qXHvZbeRVj9gWG6rhcDhgE/PO5EXdfn6Lruo+u6T4ECj3wIEEIIIV446Q3APYGV979fBsgiLCGEEOIZpDcARwIN73/fGDhrneYIIYQQL4anzgFrmvYL4Afk1zQtHPgEeA/4StM0ByCB+3O86ZGcnEx4eDgJCQnpvUS24eLigpeXF46OjrZuihBCCBtLyyroLo95qZY1GhAeHk7u3LkpWbIkmqZZ45J2Sdd1oqOjCQ8Pp1SpUrZujhBCCBuzeSashIQE8uXL91wHXwBN08iXL98L0dMXQgjxdDYPwMBzH3xTvCg/pxBCiKeziwAshBBCvGiyXQBefTiC+hO2UuqDddSfsJXVhyMyfE2j0Yi3tzdVqlShRo0aTJ06FYvF8sT3hIaG8vPPP2f43kIIIV5M6c2EZROrD0fw4crjxCebAYiIiefDlccBaFfTM93XdXV15ciRIwBcv36drl27cvv2bcaOHfvY96QE4K5du6b7vkIIIbLG6sMRBGw8TWRMPEXdXRnVvEKG4oY1ZKsecMDG06nBN0V8spmAjaetdo+CBQsyZ84cvvnmG3RdJzQ0lAYNGvDSSy/x0ksvERgYCMAHH3zArl278Pb2Ztq0aY89TwghhG2ldN4iYuLRedB5s8YIakZkqx5wZEz8Mx1Pr9KlS2OxWLh+/ToFCxZk8+bNuLi4cPbsWbp06UJwcDATJkxg8uTJrF27FoC4uLhHnieEEMK2ntR5s2UvOFsF4KLurkQ8ItgWdXe1+r10XQdUopBBgwZx5MgRjEYjZ86ceeT5aT1PCCFE1sqqztuzylZD0KOaV8DV0fjQMVdHI6OaV7DqfS5cuIDRaKRgwYJMmzaNQoUKcfToUYKDg0lKSnrke9J6nhBCiKz1uE5aZnTenkW2CsDtanoyvkM1PN1d0QBPd1fGd6hm1SGEqKgo+vXrx6BBg9A0jdu3b1OkSBEMBgMLFy7EbFbDGLlz5+bu3bup73vceUIIIWwrqzpvzypbDUGDCsLWHrOPj4/H29ub5ORkHBwc6N69O8OHDwdgwIABvPHGGyxbtoxGjRqRM2dOAKpXr46DgwM1atSgV69ejz1PCCGEbaXEDHtbBa2lzHVmBR8fH/2fC5NOnjxJpUqVsqwNtvai/bxCCPEi0zTtoK7rPo96LVsNQQshhBDPCwnAQgghhA1IABZCCCFsQAKwEEIIYQMSgIUQQgjgbPTZLL2fBGAhhBAvtCNXj9Dwp4ZUnFmRczfPZdl9s90+YGuLjo7G398fgKtXr2I0GilQoAAA+/fvx8nJyZbNE0IIkQlMFhO34m9RIGcBXB1cuRRzianNplIkV5Esa8MLH4Dz5cuXWorw008/JVeuXIwcOfKhc3RdR9d1DAYZMBBCiOwsPjmeH4/8SEBgADUL12TlmyupkL8C5wefx2gwPv0CVmRXAXjoULgfC63G2xumT3/29507d4527drxyiuvsG/fPlavXk2NGjWIiYkB4Ndff+XPP/9k7ty5XLt2jf79+xMWFobBYGDGjBnUrVvXuj+IEEKIdLsVf4tZwbOYvnc6UXFR1PWqSy/vXqmvZ3XwBTsLwPbmxIkT/Pjjj3z33XeYTKbHnjd48GBGjx5N3bp1CQ0NpVWrVoSEhGRhS4UQQjzJ9L3T+WznZ7Qo24IPX/mQBsUboGmaTdtkVwE4PT3VzFSmTBlq16791PP+/PNPTp8+nfrnW7duER8fj6urbSttCCHEi+r0jdMEBAbQrmI7WpVvxX9e/g8dKnWgRuEatm5aKrsKwPbm7wUVDAYDf8+bnZCQkPq9ruuyYEsIIezA/oj9TNwzkVUnV+Hs4Ey1gtUAyJ8jP/lz5Ldx6x4mq4rSyGAwkDdvXs6ePYvFYmHVqlWprzVp0oSZM2em/vmItSeyhRBCPFXv33rz8tyX2XpxK/9t8F8uDb3EkLpDbN2sx5IA/AwmTpxIixYt8Pf3x8vLK/X4zJkz2bNnD9WrV6dy5cp8//33NmylEEK8GMwWM8tPLCfBpEYk/Uv5M7npZMKGhvFF4y8omLOgjVv4ZFKOMIu9aD+vEEJYW4IpgflH5hMQGMD5W+dZ0G4B3Wt0t3WzHulJ5QhlDlgIIUS2YLKYmBw4mel7p3Pt3jXqeNYhoGkAbSu2tXXT0kUCsBBCCLsWlxxHDsccGDUjK06uwLuwN2Pqj8GvpJ/NtxJlhARgIYQQduls9FkmB05m+cnlnBl0hnw58rG953ZyOuV8+puzAQnAQggh7EpwZDAT90xkxYkVOBmd6O3dG5NFJUN6XoIvSAAWQghhRy7eukjt72vj5uzGB698wOCXB1M4V2FbNytTSAAWQghhM2aLmZUnVxJyPYSxjcZSKm8plnZcSrMyzcjjksfWzctUT90HrGnaPE3TrmuaFvKP4//RNO20pml/aZo2KfOamPnCw8Np27Yt5cqVo0yZMgwZMoSkpCS2b99Oq1at/nX+2rVrqVmzJjVq1KBy5crMnj3bBq0WQojsK8GUwJyDc6g4syKdl3dm6Ymlqft5O1Xp9NwHX0hbIo6fgBZ/P6BpWiOgLVBd1/UqwGTrNy1r6LpOhw4daNeuHWfPnuXMmTPExsby0UcfPfL85ORk+vbty5o1azh69CiHDx/Gz88vaxsthBDZ2PbQ7ZT6qhTvr30fdxd3lnVaRkj/EFwcXGzdtCz11CFoXdd3appW8h+H+wMTdF1PvH/Odau16FHBrHNnGDAA4uKgZct/v96rl/q6cQM6dnz4te3bn3i7rVu34uLiQu/evQEwGo1MmzaNUqVK0ahRo3+df/fuXUwmE/ny5QPA2dmZChUqPP3nEkKIF9jV2Kvcir9FpQKVKJ+vPDUL12SE7wgal2qcrbcSZUR6U1GWBxpomrZP07QdmqY9tmSQpml9NU0L1jQtOCoqKp23yzx//fUXtWrVeuiYm5sbxYsX59y5c/8638PDgzZt2lCiRAm6dOnC4sWLsVgsWdVcIYTIVs7dPEe/tf0oOb0kA9cPBKBo7qKs77Ye/9L+9hF8szAj5N+ldxGWA5AXqAvUBpZqmlZaf0ReS13X5wBzQKWifOqVn9RjzZHjya/nz//UHu8j2vfIfwCPOw4wd+5cjh8/zp9//snkyZPZvHkzP/300zPdVwghnmdHrx5l3O5xLD+xHEeDIz1r9GRU/VG2btbDEhNh0SJVC/f336FUqSy9fXp7wOHASl3ZD1gA+6rzlEZVqlThn/mp79y5w+XLlylTpsxj31etWjWGDRvG5s2bWbFiRWY3Uwgh7J6u61h0NSK45eIW/jj3B6PrjSZ0aCizW8+mrEdZG7fwvpgYmDhRBdw+fcDBQU1hZrH0BuDVQGMATdPKA05A1rfeCvz9/YmLi2PBggUAmM1mRowYQa9evciRI8e/zo+NjWX733rZR44coUSJElnVXCGEsDspVYnqzK3DgqPqd2k/n36EDQ1jfJPx9rWPNy4OypaFDz6AqlVh0yY4dAhqP3YmNdOkZRvSL0AQUEHTtHBN094F5gGl729N+hXo+ajh5+xA0zRWrVrFsmXLKFeuHOXLl8fFxYVx48YBsGXLFry8vFK/Dh8+zKRJk6hQoQLe3t588sknMvwshHgh/X0rUadlnYhJiMHN2Q2AHI457Gcr0V9/qR4vqKnMceNU0N20CZo2BRvNQ0s5wiz2ov28QojnV/NFzdl0fhM+RX0YU38M7Su2x2gw2rpZiq7Dzp0QEADr1qnAe+YMeHpmaTOeVI4wvUPQQgghXjBXY6/y0ZaPiEmIAWBM/TH82f1P9vfZT8fKHe0n+J45A3Xrqm2t+/fDZ59BWFiWB9+nkVSUQgghnujczXNMDpzMT0d+ItmSTK2itehQqQONSzW2ddMeiI+Hy5ehfHkoUgQsFpg1C3r2BFdXW7fukSQACyGEeKQkcxLdV3VP3UrUy7sXI+uNtJ/VzADR0fDtt/D111CwIBw7Brlzw4EDtm7ZU0kAFkIIkUrXdU5EnaBKwSo4GZ3QdZ3R9UYzpO4Q+1rNHBoKU6bAvHlqZfPrr8OoUTZbUJUeEoCFEEJgtphZcXIFE/dM5Ni1Y5wffJ7ieYqztNNSWzftYRYLGAywezfMng3dusHIkVCliq1b9sxkEZYQQrzA/r6V6M3lbxKbFMus12dRKGchWzftAV2HjRuhSROYfL/2z5tvwsWL8OOP2TL4ggRghg0bxvTp01P/3Lx5c/r06ZP65xEjRjB16lSqVq36r/fu3buXl19+GW9vbypVqsSnn36aFU0WQgiruRZ7jYHrB+Lu4s7yTss5MeAEfV7qg7ODs62bBsnJsHAheHtDixZw8iTkzatec3S0u1XNz+qFD8D16tUjMDAQAIvFwo0bN/jrr79SXw8MDKR+/fqPfG/Pnj2ZM2cOR44cISQkhM6dO2dJm4UQIr2u3L3CmM1j6LSsEwAl3EtwrN8x9vfZzxuV37CfrUSgqtz16AEmk+rpXrwI771n61ZZjd3NAfv95PevY52rdGZA7QHEJcfRcvG/yxH28u5FL+9e3Ii7QcelD5cj3N5r+xPvV79+fYYNGwaoykhVq1blypUr3Lp1ixw5cnDy5Enypjxx/cP169cpUqQIoMoYVq5cOQ0/oRBCZL0z0WcI2BPAgmMLMFlMdK7SmSRzEk5GJyoVsJPkQFevwowZ0L8/FCsGQ4ZA167w2mtq3vc5Y3cBOKsVLVoUBwcHwsLCCAwMxNfXl4iICIKCgsiTJw/Vq1fHycnpke8dNmwYFSpUwM/PjxYtWtCzZ09cXF6sgtJCCPu34sQKOi3rhJPRiXdrvssI3xGU8Xh8sZksd+qUWtG8YIHq7VaooPbv1qlj65ZlKklFCXTr1o3WrVuzYcMGhg8fTkREBIGBgeTJk4fo6Gj69etHq1atCAkJ+dd7z58/z6ZNm/j111/RNO2hQg2PYg8/rxDi+abrOpsvbMaoGfEv7U9MQgxTAqcwqM4gCuWyo8VVFgt07gwrVoCLixpyHjFCFUt4TkgqyqdImQc+fvw4VatWpW7dugQFBT1x/jdFmTJl6N+/P1u2bOHo0aNER0dnUauFEOJhJouJJSFLqDWnFs0XNScgMAAAdxd3Pm/8uX0EX7MZdu1S3xsMKmvVxx/DpUsqc9VzFHyfRgIwah547dq1eHh4YDQa8fDwICYmhqCgIHx9fR/7vnXr1pEygnD27FmMRiPu7u5Z1WwhhEi1/MRyKnxTgbdWvEVcchw/tPmB3976zdbNeiA+Xu3brVQJXn1VZawClcFq7FiVxeoF88LPAQNUq1aNGzdu0LVr14eOxcbGkj9/fmJjYzl9+jReXl6pr0+bNo0VK1YwbNgwcuTIgYODA4sXL8ZotKMVhEKI59qt+Fs4GZ3I6ZST2KRYCuQowOSmk2lbsS0GzU76V3fuqIVVX38N16+Djw8sWQKyaFXmgLPai/bzCiGsL/xOONOCpjHn0BzG+o1luO9wLLoFDQ3NXlIxJiervbq3bkGJEvDKKzB6NDRsmK3SRWbUk+aApQcshBDZxMmokwQEBrDo2CIsuoU3q75J09JNAeynx3vokKrBe+6cKgWYNy+cPw8FCti6ZXZHArAQQmQTA9cPZG/4Xt6v9T7DfYdTKm8pWzdJ0XXYtEkF3i1bVDWi99+HxES1ulmC7yPZRQDWdd1+hk0yUVYO9wshsjdd19lwbgPT9k5jfrv5FM1dlFmvz8LD1YMCOe0soC1frrYTFS0KkyZB376QJ4+tW2X3bB6AXVxciI6OJl++fM91ENZ1nejoaEnUIYR4omRzMkv+WsKkPZM4fv04Xm5enL95nqK5i1IhfwVbN0+5cwe+/x7y51cJM9q0UUk03nwTHpO4SPybzQOwl5cX4eHhREVF2bopmc7FxeWhldRCCPF38cnxVJtVjfO3zlO5QGV+avsTXap1wcloJ0EtMhK++kptJ7p9G7p3VwHY2Vl9L56JzQOwo6MjpUrZyTyGEEJksei4aP449wfdqnfD1dGVXt69qFGoBq+Xf91+FlaBShX54YcqkUbHjjBqlNpSJNLN5gFYCCFeRJdiLjE1aCpzD88lPjmeV4q/Qgn3Evzv1f/ZummKrsOOHVCxIhQuDFWrqrnd4cOhdGlbt+65YEePV0II8fyLuBNBj1U9KPt1Wb4N/paOlTtyrP8xSriXsHXTFJMJli5VhRAaNVLpIQGaN4dvvpHga0XSAxZCiEym6zq3E2/j7uKOi4MLm85vYlDtQQzzHUbxPMVt3bwHZs1SW4kuXoRy5eC771Q9XpEpJAALIUQmsegW1p5Zy4TdE0i2JLO/z37y5chH2LAw+1lYdecOuLmp77dtg0KF1HxvmzYgqXUzlQxBCyGElSWZk/jx8I9U/bYqbX9ty5XYK/Ss0ROLbgGwj+B79iz066fmd0+cUMd++gkCA6F9ewm+WUB6wEIIYWWLji3i3d/fpUahGizusJjOVTrjYLCTX7dBQWqYefVqtWe3Rw/IlUu9liOHbdv2grGTfxFCCJF9XYu9xox9MyiXrxy9vHvRtVpXPHN70qxMM/tKMBQTA/7+Kj3kRx/BoEFqyFnYhARgIYRIpwu3LjA5cDI/HvmRRFMig18eDICLgwvNyza3ceuAhASVoWrHDli0CNzdYf16qF0bcubMsmasPhxBwMbTRMbEU9TdlVHNK9CupmeW3d9eSQAWQoh0+GzHZ4zdMRYHgwM9qvdgZL2R9pMqMjoavv1WbRu6fh1q1VJlAT08wM8vS5uy+nAEH648TnyyGYCImHg+XHkc4IUPwrIISwgh0kDXdbZc2ELUPZU2t3bR2oz0HcnFIRf5vs339hN8d++G4sXh449Vpqpt2+DAARV8bSBg4+nU4JsiPtlMwMbTNmnPkwRdDmL8rvEEXQ7KkvtJD1gIIZ7AbDGz8uRKJu6ZyMErB/m80ef879X/8Vq513it3Gu2bp5y4IDq4TZrpnq7vXrBgAFQpYqtW0ZkTPwzHbeVwLBAGi9ojMliwsnoxJYeW/At5pup95QALIQQj/H9we+ZFDiJczfPUdajLLNbzaZHDTtJTGGxwIYNakXzjh3w0ksqALu6wsyZtm5dqqLurkQ8ItgWdXe12j2sMcc8YtMIEs2JgNpGtj10e6YH4KcOQWuaNk/TtOuapoU84rWRmqbpmqblz5zmCSFE1opPfhAs1p9bj7uLO8s6LePUwFP0rdUXFwc7KCm6bp3KzdyqFVy4oBJnbNtm61Y90qjmFXB1fHhPsaujkVHNrTNknzLHHBETj86DOebVhyOe+L67iXeZFjSNC7cuAPBerfdwMjhh1Iw4GZ3wK+lnlfY9SVp6wD8B3wAL/n5Q07RiQFMgzPrNEkKIrBVxJ4Jpe6cx99Bc9r+3n/L5yrOg3QJyOeWyj61Et26BwaAK3Scmqj28ixZB587g6Gjr1j1WSk80s1ZBP2mO+Z/3CLocxNozawm/E87vZ34nJiEGHZ3hvsN5p+Y7VMpfie2h2/Er6ZfpvV9IQwDWdX2npmklH/HSNGA08JuV2ySEEFnmZNRJAgIDWHRsEWbdzJtV3sSoqR5bbufcNm4dEBoK06fD3LkwYgSMHQvt2qlsVfbwYJAG7Wp6ZtqK57TOMQeGBdJwfkNMFhMAfiX8mNBkAi97vZx6jm8x3ywJvCnSNQesaVobIELX9aNPezLUNK0v0BegeHE7SjouhHjh3Um8g8/3Plh0C31r9WWE7whK5bWT+uSHDqn53WXLVKDt0kXV4QXVExbA0+eYQ66HULVgVXZc2oHZonrKRs1IszLNHgq+tvDMf4uapuUAPgI+Tsv5uq7P0XXdR9d1nwIFCjzr7YQQwmp0XWf92fUM2TAEADdnN35941fChobxTctvbB98df3B9599puZ6hw1T87wLFkC1arZrm5161Byzi6NGI+9wGv7UkGqzqrEvfB9+Jf1wcXD59xxvQoLKgZ2UlOVtT08PuAxQCkjp/XoBhzRNq6Pr+lVrNk4IIawh2ZzMryG/MilwEiHXQyjmVoyPXv2IgjkL0rpCa1s3T83p/vwzTJsGK1aoUoAzZqj53jx5MvXW2T1LVUpbP96wkstx+8npohOXI5gv952mmFsxpjWfRuUClcntnJstPbY8mON1LqOG82fOhKgoyJsX2rbN0rY/cwDWdf04UDDlz5qmhQI+uq7fsGK7hBDCKo5ePUqbX9sQdjuMKgWqsKDdAt6q+haORjtYuBQTo2ruzpgBV66oHm50tArAWTBll12yVD3tIaFgvkuc08eQ5JREjNlMaYfSLGy/kDervPnQ37NvMV98PaqrUYUFC9SDT6tWam69YcMs/7meGoA1TfsF8APya5oWDnyi6/oPmd0wIYRIr6h7UVy6fQmfoj6U9ShL9ULVmdlyJi3LtcSg2cn8aUKCCrQ3bkCTJmoYtGnTLF1Y9SwriG3lSQ8JL5c1MGPfDH4+/jNJ5iTMuhmjZuTdmu/ydvW3H1xE1+HSJShZUlV8OnwYevZUgbhiRRv8VEpaVkF3ecrrJa3WGiGEyICLty4yJWgK8w7Po3ie4pwceJKcTjlZ02WNrZumHDoEa9bAJ5+oikSTJkHNmuDtbZPmZIcsVY96SLhjusT7a78mRvsTk8WEXwk/ouKiSDIn4WR0olHJRurE5GS1iG3KFDh/Hi5fhty5Ye9eu6h3LJmwhBDZ3omoE3yx8wuW/rUUg2age/XujKw30j727+o6/PEHTJ4MW7eqAPDuu+DlBb1727RpWZGlKqNSHgYSDSdJMBxH03Nwy3E2mtmRfrXfZbjvcMp6lCXoctCD+V23yurz/uorCA9XvdyAALV3Guwi+IIEYCFENqXrOiaLCUejIyHXQ1hzZg1D6w5lWN1heLrZx/ApJ0+qRBkhIeDpqXq8fftm+sKqtBrVvMJDw7tg3SxV1lAkjzMnY38mxnEhYEHDgVzmFlTM8S7fvv5G6nm+xXzx9aqrhvAPHoRRo6BRIzXH/tprdrl1SwKwECJbMVvMrD61mol7JtK6fGv+r+H/8UalN2hauil5XfPaunlqYVVoqBpWLl4c8uWD+fPhrbce9MDsRGZnqcqIJHMSPx//mXDn8cQknQEd0EDXTbhqhfioRd0HJ+/fr4aZc+dWCUtq1YJTp6CC/TxIPIoEYCFEtpBgSmDB0QVMDpzM2ZtnKetRltJ5SwNgNBhtH3wvXXqQscrTU/V+c+aE7dtt266nyMwsVem18uRKBm8YTMTdCGoUqkG9IsNZdmYmFj0Zg+bI8IbtaVe9MKxerQLv7t1qVOE//3lwETsPviABWAiRTfT5vQ+Ljy+mVpFaLO24lA6VOmA02MFcXkgIjBsHS5c+yFg1YkS2SRNpL9acXsPeiL20KtcKdxd3KuavyLy282hauimapjHkcseH8zR/+qnax1uihHrweecd1QPORjT975lXMpmPj48eHBycZfcTQmRfkXcjmb53Ov19+lMqbymOXD3CzfibNCrZyPaLqywWtcLW2RmWLIH33oP334fBg6FYMdu2LZs5deMUozaPYu2ZtWhouDi4PLoW79Wr8M030Lw5NGigRhz27YMOHcDBfvuSmqYd1HXd51Gv2W+rhRAvrNikWGp8V4Ob8TepXKAypfKWwruwbbbqPCQxERYvVsOeXbvCRx/BG29AixZ2s7Aquwi6HMTEPRP57fRvOBgc0NDQ0f9dizckRH3eP/+sHnpy5lQBuEQJ9ZWN2d+yMCHEC+/wlcPciLvBko5L6OXdy9bNUaUAx49XiRzefVeV/6tcWb3m4CDBNw2CLgcxbtc4gi4HATAreBa7wnbxScNP+O2t3x6dp7lnT5UdbOlSNcpw5gx8+KHtfggrkx6wEMLuBEeqqaoGxRvYuCX3vfsurFqlhj8XLgR/f5njfQY7QnfQdGFTki3JOBud2dZzGwFNA5j1+ixyOuUEUHmaz/+J33kLvkVqqze+/LJaTNWvH3h42PAnyBwSgIUQdif4SjDF3IpRKFch2zTgwAE17Dl+PJQqpRb7fPopVK9um/ZkU7cTbjP74Gy+2PkFyZZkAJItyQ8PMQPcvInvwm34fjNL5cQuVEvlaB4wwEYtzxoSgIUQdmdg7YG0r9g+a29qsajyf5Mnw86d4Oam5nlLlZIygOlgspioOqsq4XfC8Snqw/FrxzFZTA8PMcfGwgcfwI8/QlwcNGv2ICf2C0ACsBDC7tQrVi9rb2gyqeQNx46pVcxTp6phZze3rG1HJsiKcoMpaSBLuJfgZNRJPmv0GQ4GByY2mUiFfBWoVbTWg1SRJRriq91fKZ4jh3rY6dwZhg9/4R50ZBuSEMKuXLx1kdPRp1MLqGea6GhYu1Yt9AH48ksoXRo6dlSLrJ4D/6wkBCrV5PgO1awWhIMuB9FofiMSzYkAOBmdONrvKBXz/6PKkMmk5tEnT4azZ1VhhJw51XE73kaUUU/ahiSroIUQdmX5ieW8tvg17iXdy5wbnDsHAweqnm6vXioYgNpS1KXLcxN84cnlBq0h7HYYb614KzX4amiMqjfq4eB7964qilCunOrpRkfD558/KIjwHAffp5EALISwK/OqPP8AACAASURBVAciD1DKvRT5cuSz7oUvX1Z7dsuXV+kiu3RRe0zLlbPufexIZpQbTDQlEnI9BIDCuQrj7uKOo8ERo2bExcGF18u9rk5MGV396y8YOlSl51y5Ek6fVg9ALpk4upFNvLiPHkIIuxQcGYxP0UeO2D07sxkiI1VvN08eVY/3ww9h0CAoUsQ697Bj1io3GHQ5iA3nNnAj7garTq3CqBm5MORC6nDzQ6UAb7jAh2+rz3vmTKhbF44elRXkjyABWAhhN6LjorkYc5F+Pv0ydqG4OFWBaOpUlS7y+HG1oOrcObupBZsVrFFu8PfTv/PG0jcwWUwA1PGsw7jG43A0PBiq9/V8Gd9jN6HnR7BtG+TKpR5yUkjwfSQJwEIIu3HwykEAahetnb4LXLumel3ffqvmGuvUUXVhdV0lzniBgi9krNygRbdg0AysOrkqNfgaNSPtKrTDv7R/6nmrD0dwc+SHvLN1Idfd8nN9yEdU/XQkuLtnzg/1HJFV0EIIu2GymDgRdYKyHmXJ4Zgj7W9MCbALF6pVzW3awMiRUL++ZKx6BrqusztsN5MCJ1GlQBUmNJlAYFgg/gv9STYn42R0UoUSXMrCrFnsKOFNv3NOFLgeTs3IU6yr2ABHF2errrLO7qQYgxAiW3AwOFC9UBqHK3Vd7SGdPBkaNlQB9803VfrC8uUzt6HPmd1hu/ku+DuOXD3CX1F/kT9HfvxK+AFQr3g9tvbYquZ4DaXw/eInWLAAEhI41ewd4mt2ICxvEcLyqjl10/1V1tkpAMfHw+HDsH+/+vPQoVlzXwnAQgi78d8t/6VF2Ra8WuLVx59kMsGKFSpV5IEDkD+/yqAE4OQkwfcZpezjNVlMaGiM8B3BZ40+e2gEwreYL74fz4Gf/qvm1Hv0gGHDmDD/wiOvmZFV1pnNbIYTJ1SwTfk6flwdB/X8JgFYCPFCuRp7lfG7x5M/R/4nB+B33lFDzeXLw3ffqWDg+myrep8X6c1yFZMQw3fB39G6fGu2h27HolsAMGgG8rnmU8E3ORl+/x3atVNz59WrwyefqPzMBQsCUNT9ilVWWWemK1dg794HXwcPwr37W8zz5FHLBD74AGrXVl9Fi2Zd2yQACyHswsHIxyzACg+Hr79Wq2qLFVN7SDt2VMn6DS9uKoN/ZrmKiInnw5XHAR4bhMPvhDN973RmH5xNbFIsRs2IX0k/nI3OJJmTVJ7m/D4QEAAzZqjPft06aNkShg371/WsscramhIT1VBySrANCoKwMPWaoyN4e6vntzp11FfZsg/+Ca0+HEGnBZmbsvOfJAALIexCcGQwGho1i9RUB44eVcPMv/yiCiVUrQrdu6sxQvHELFf/DBxBl4MYvnE4ByIPAPBm1TcZVW8U3oW9gfulAM9sxG/DSXzrdFBFEho1glmzoEWLx7YhI6usreHqVQgMfPB18CAkJanXihdXW5CHDlX/rVnz8bk/0vMwYw0SgIUQdiH4SjCVClQil9FV9bg2bFC5ggcOhCFDVFUikeppWa50XedA5AFMZhNNFjYhwZSA0WBkSccldKjU4cEbrl1Tc7xF68DAKtC2rSqM8NJLaWpHu5qeWRJwzWaVVGvPngcB98L9KWhnZ/DxgcGDoV499Yz2LEPJz/IwY00SgIUQtpeYyNWI0/iU9VXzjRUrqpXNfftC3ry2bp1delyWqyJ5nFh5ciWT9kxiX8Q++tTsQ5I5CR0dXdc5feO0imZr1qgRhtOn4dIlNY9+7JhayGYHEhLUAqndu2HXLhVw79xRrxUqpHaYDRigAu5LL6kgnF6ZkbIzLSQACyFs5+ZNmD0bvv6aA1eukHR8qTo+dapt25UNpMy/xphDSDAcx9lSCc0hklDHtbyx9CKl85bm25bfUil/JRYfX/xgjvdQNLxTUWUFK1EC/vvfBxe1YfC9dUv1blMCbnDwg+HkKlVU6u769dVXqVLW3d5trZSdz0oCsBAi60VFqYo4P/zwUCF2pyo1bN2ybKNdTU9O3TzIR7v/h0VPxqA54urgQkW3snz12kQ6VOqA0aAyf23psUXt473tjm/rAWoF0pIl0KGDzaoRXb+uAu2OHWo797Fjamu3o6MaTh46FF55RfVw81m5Lsc/2WoxmQRgIUTWuX1b7f1wcIBFi6BTJxg+nK/itrEv4icW0xTJW5U2YbfDWBM6EQuJoIGmmRj0cj/G+49HS+kehoTA1Kn45s2L75QpKsIF11FjtlmcISwiQgXalIB78qQ6niOHCrJjx8Krr6png6zeVWarxWQSgIUQmctsVvtJp0xRvd2DB9W8bkpBdmDDotFcib3yIHCIxzp+7TgBgQH8EvILFosFo6Z6uU5GJ9pWaKseYDZtUp/3pk0qmg0cqN6saVCrVpa08+pVVZdh+3b135Syy25uqmfbs6ea5n/pJfuYds6qxWR/JwFYCJE5/l6R6Nw5KFlSjSuazaoHfD/46rpOcGQw7Sq2s2177dTfS/3p6NSfV5+cjjkZVHsQQ+sOJfJu5INSgMV84eOP1fB+4cLwxRfQr1/mj+GihpS3b38QcE+dUsfd3FTPtl8/FXC9vV+4mhiPJQFYCJE5li5Vy1Tr1FHft2//yPnGS7cvER0fbb0awFaS3ixT1rQ7bDf+C1QhBBcHFzZ338yMFjPoVr0bHq4eAJSw5MZ393bI4QTFgG7d1Cqlrl0ztjT4Ke7cUcPJW7bA1q0qnSOoSoQNGqiEF40aqYBro2lmuycfixDCOk6cUL3dl15SgbdLF5Vq6CkViYIjVYU0ewrAtkrMkCI+OZ75R+fz0ZaPSDKrpcBJ5iR2XtrJhw0+VCedPw/Tp8O8eWq0QdfV8HKFCurLyhISVGapLVvU14EDajDDxUUNKXfpogJurVpqIZV4uqcGYE3T5gGtgOu6rle9fywAaA0kAeeB3rqux2RmQ4UQdkjX1Xjj5MkqcYaLi0pBBKr39corT72EUTPysufLVCtYLZMbm3a2SswAsOLECvqv609UXBSV8lciNjkWs8WsthCV9FMn9emjAq+Dg+rxDh8O1az7+VksamXypk2webPaHpSQoIaPa9dW+ZP9/cHX9/EZpsSTpaUH/BPwDbDgb8c2Ax/qum7SNG0i8CEwxvrNE0LYtX79YM4clZz/s8+gf39VnegZtK/UnvaV2mdSA9MnKxMzBF0OYvWp1dQrVo+2Fdvi5eaFT1EfxtQfw6slXmVv+F62X9iKX6QTvp7303CWKwcffqjyYxcpYrW2hIerYLt5M/z5p9otBlC5Mrz/vgq4r76qFrJbkz0M99vCUwOwrus7NU0r+Y9jm/72x71AR+s2Swhhl2Ji4PvvVa+raFE17li7Nrz9drq6QbquY9EtqftV7UVWJWZYeHQhvX/rjVk3Y9SM7Oq9C99ivqzvtl6dcPcuviv34zt9LoSGgkd1aN4cxlinv3PvnprH3bhR9XRTFk4VKqS2ZjdtCk2agGcmxkJbD/fbkjVKibwDbHjci5qm9dU0LVjTtOColMcpIUT2EhqqquEUKwajR8Pateq4n58aDk3nGOS5m+fIMyEPa06vsVpTrWFU8wq4Oj78UGDNxAw7Qnfw2uLX6LG6B2b9wVD39tDt6pvYWBVkixVTK8c9PWHlShUNM0DX1bByQIC6lIcHvP66eqYqUULNJBw9qkr4LVqktgplZvCFJw/3P+8ytAhL07SPABOw+HHn6Lo+B5gD4OPjo2fkfkKILGaxqN7tkiWqbtubb8KIEaq0jBUERwZzL/kexfMUt8r1rCUzEjOYLebUnv6Sv5Zw6Moh3q/1PguOLniQJjLv/c/VxQWWL1fd0BEjMlQB6uZN1btN6eVGRqrjVavCf/6jOtQNGthuHtdWeZgfYjLZZKl2uu+oaVpP1OIsf13XJbAK8bywWGDfPrW6xmBQ+3VHjFClZry8rHqr4MhgXBxcqFygslWvaw3WSMwQdDmIzRc2czfxLqtOrWJ+u/nUL16fLxp/wdTmU3FxcKFn9e5s3zQHv7Uh+M58By5eVAvYQkLSlRLKYlE1cTdsgPXr1V+lxaJynzRtqqoLNmuW+T3btLJVHmZAdfVnzVLrGHbtUnPrWShdAVjTtBaoRVcNdV2Ps26ThBA2kZI4Y9o0lbboxAmoVEmNT2aS4CvBeBf2xtH4/O1b2XhuI61+aYXJYgKgUv5Kqa95uHqoJcU/zcV36lR8T55UETElUQk8U/C9dUv1bjdsgD/+gGvX1M4vHx/43//gtdfUVL09JsCwSR7mQ4fUFq5ff1W931atIDk58+73GGnZhvQL4Afk1zQtHPgEterZGdh8P3XcXl3X+2ViO4UQmeXWLRV0v/0WoqPVb+pff8303oDZYubQlUP0qtErU+/zOJm58taiW3hr+VupwdegGXi7+tvUL17/wUmBgfDeeypTRUpe7DTmZNR1VRt33To1HR8YqHq5Hh5qSPm119R/Cxa0yo+TqbI8D3NMjEo+7eCgVvH/5z9Z3vNNoWXl6LGPj48eHBycZfcTQjxBQoKa+LtxQ6WJbNJEDTW/8kqWJOqPS45j0p5JNCjeAP/S/pl+v7/758pbUL2u8R2qpfsX/5GrR5h/ZD6Tm03GaDAyYfcExu4YS7I5GSejE1v85uG7YJtKCzlunIqiQUFqqD8Nn3dCgtpyvXatCryXLqnjL70ELVuqrzp17LOXa1N37qg90/v3w88/q2MbN6p5dXf3TL+9pmkHdV1/ZJYZCcBCvEh0XeUNnDJFrc4JClK//KOjsyRfsL2oP2HrI+cdPd1d2fNB4zRfJzAskB8O/0BIVAj7I/aTyykXge8EUq2QSooRFBbI9u0/4rf+JL5L9qi53YED1eefBpGRKuCuXauyT8XFqepBTZuq1cstW9rPXK7dOX8evv5aBd+7d9WD5bp1Kjl1FnpSAJZUlEK8CJKS1LDy1Klqn0mhQiqJQ0phBBsE33M3z1EgRwHyuFg5q0MaWGPl7W+nfqP9kvboqE5Mv1r9GOc/jryueVPP8Z37B76fz1XJST75RKXofMK4sK6rv541a1QBqZT+SsmSKrdyq1aqoIFknnqK9evVh2U0wltvwZAhakLczkgAFuJFsGgRvPuuSmn0ww8qUb+Nf4v3WNUDB4MDO3vvzPJ7p3flbVxyHCeiTuBT1IeQ6yGpwdeoGSmepzh5kwwwc7JaZly9utq25eUF3bs/clHV6sMRTFx3hovHcqBd9sQUWpgbVx3QNDVCOm4ctG4NVapkefne7CUxUT1g5s4NHTqop5SPP4a+fVXCGDslAViI59GFC2qVp7e36jp16aJ+ETVvbhe/yU0WE4evHqa/T3+b3P9ZV95Gx0Uz88BMvt7/NQCXh12mcanGuO5yVXt4DY74/X4UXi+mhjsTE1ltzkfAxmtExnhS9KughxYWxcTA57Nu8v0iA7HnGqAnOaA5mshZOpqBfZ34vwF5KVQo8z+HbO/6dfjuO7WA8No1aNdOBeCcOeHTT23duqeSACzE8yQoSM0vrlql9vCOHq2Ou7qqDaB2IjgymARTArWL1rbJ/dO68nb1ydVMCpzEoSuHSDQn0qp8K8bUH4Oz0RnfYr5s6bGF7dOH4Lf8IL4Ry1MTlazWCv0rveLIn86wYWlOLh50Z9s2MJk8MORMIGelSFzLXcOl+A0MjhYOu7pSqFDa56FfWAEB8H//p3q/LVuqLVwZyBSWbE5mW+g2mpVpZsVGPpkEYCGeFwMHqp6AuzuMGqW2V9jhCp3pe6czevNocjjmoEGJBjZrx5MSbZgtZvZH7OetFW+RaE7EqBlZ1H4R3ap3U/PmmzZBs2b4FvPF1+116NxQJSopVgyAgAlbiU82k3wjF3FnCxF3tjBJV9w5D5Qvrxabz4/Yg1PRmH8NSGRpBqjsxGJRc7t166o59bJloXdvNb9bsWKGLr3x3EbeW/Mel+9c5mi/o1QvVN1KjX4yCcBCZFexsfDjj6rXVbAgtG2rfhH17q2qotuRAxEH8HLzokjuIlQrWI33XnqPEfVG4OVm3cxaGaHrOlsubmHSnkmUz1cez9yeqft4AcJunFMPONOmwblzqlyQv79aXJV6DZXj4a/fi3PvdGFMN9Xfg1ORW7g3PEWOclc5/b0fALsmJBLxiCKuWZIBKjuJjVUJYr76SiWImTxZPcG0b6++0iniTgTJlmRKupfE082TMh5lmPX6LKoWrGrFxj+ZBGAhspuICLW9YvZsNZno5KRqxTVrpr7sREpAm7B7AlsubmFM/TFMaDIB/9L+Wb7v93GCLgex9eJWzLqZ1adWc/jqYQrnKszr5V6njmcdnIxOao7XouH3nylw4u6DRCUNGwKqQxwYqGolrFwJYWGAoTQuxW7iVisU13JXccidCKhtTilskgEqO9F1VXR49my4fVutSvvlF3jjjQxd9vi140wOmszPx3+mY+WO/PLGL1QtWJVtPbdZqeFpJwFYiOzCbFYrmRcvVsNx7durnoCvr61b9i8rT65k3K5xHLxykCK5ihDQNIC+tfraulkPCbochP8CfxJMCejoFHcrztzWc3m7+ts4OzjD7dtqjvf8FvxGzMC3XD34TiUqMZk1tm9X9RJWr1brf5yd1fPP2LFgKHGN8VuPPjG4ZnkGqOxA1+HMGahQQS0WPH1aLRwcNkwNPWfAjtAdjN89no3nN5LTMScDfAYwtO5QKzU8fSQAC2HPLBY1punjo/Y0mkxqL+mQIVC6tK1b9xCTxYSDQf1KWX5iOXcS7/B96+/pXr27Cmh24kbcDWbun0lUXBRJ5iR0dAyagfd93ufdmu+ooeUpU+DECXzPn8e3mC/sGEqycy62boXlfdUat+hotdi2ZUvVKWvZUu2CUYrg5m55anC1RsGH50JyMqxYoYb3DxxQQ/ylS6tjGUjtlWxOxsHggKZprD+7niNXj/Bl4y/p59NP5eO2NV3Xs+yrVq1auhAiDeLidH3OHF2vVEnXNU3Xz5yxdYse63bCbX3S7kl60SlF9WNXj+m6rus3427qJrPJxi172IWbF/RB6wbprl+46nyK3n9tf931C1fdONaou37hqgd+9z9dr15d10HXCxfW9S++0BNu3tPXrtX1Xr10PW9e9VLu3Lretauur1ql/ppEBty+resTJui6l5f6cMuV0/VvvtH12NiMXTbhtj4lcIpebGoxfcPZDanH4pPjrdHqZwIE64+JidIDFsKexMSo/bvffgtRUWof74IFqlq6nbkWe42v9n3Ftwe+5XbibZqUboJFtwA8lA3KVoIuB7E9dDt+Jf1YcHQBcw7NwagZ6V69OyPqjaBygcp0r95dnXPNBd9Ow6FqVUxz5vFnwa78usqZ1aXU9GOePGqNW8eOKg2kZKLKoKQktXbh3j21lahBA1UWsGVLtX0unSLuRDBj3wxmH5zN7cTbNCzREDdnlXoy5b/2RHJBC2EPEhPVJOK1ayrvoL+/mt/187OLxBn/lGhKxGuaF9Fx0XSs3JEx9cdQq2gtWzcr1YQ/f+O/ezqh62YMmiNNS3SmhmcRBr88GE83T5UnePr01BSRyUk6h6fv4LuTDVm1WiMmRgXd9u1VkaImTdJcqEg8jq6rhNbTp6uk1lu3quMREVbZLqfrOuW+LsfFmIt0rNyRkb4jqe1pm33mfye5oIWwRymFEaZOVdmTdu5UOZovXbLLOnJHrh5hScgSxvmPw9nBmW9bfkuNwjUon6+8rZuWymQxMWbdbL46+H/oWjJoYNGTOXTBSL8ag/EMuQRTBsOqVegODlx+vT+fvwcrV2rcvOmHm5vq6b75purpStC1goQEVYVo+nQ4flz92x4wQC0qNBrTHXx1XWdb6DbmHZ7HvLbzcDI6Maf1HEq6l6R0XvtaH/E4EoCFyGr/LIxQsKBKopHyC8mOgq+u6+y4tIMJuyew8fxGcjvlps9LfSjjUYZOVTrZunmpEk2JzDk4h6l7pxIaE4pRL4D69WZBwwGjqQq3Rn4AWxeRnDsvW1/6kJGhgwhZXYRcuaBNGxV0mzWT4WWr++EHVfijenW1b71LFzXak07J5mSWnVjG5MDJHL56mEI5C3H6xmmqFapG41LZK4OYBGAhstqPP6pC4JUrw9y50K2bXf7WD40J5a3lb7EvYh8FcxZkXONx9K/dH3eXzK+hmhYpe3gbl2qMT1EfpgRNwcvNi7hr3XGx1EYzH6NM1G/ccX2V8PiX+eFaBc6712VaTC8sf+WkdWv49E017fiIOgkivY4cUb1df39VhKJHD/Vv3QrTKeF3wqk/rz5ht8OomL8ic1vPpVv1brg42N//P2khAViIzHbunMriU7u2+mXUrZtaVGUnhRH+LsmcxPmb56lUoBJFchXBweDArNdn0bNGT1wd7SdKLftrGV1XdMWkm3Dd5cqWHls48N4BCuQsQLsPl9Bi+3y6HP6DPEn3+D/X1/givgHXDBbcXzMwp4vq8T7YMiQyzGJRRYunTYPt21XR4qr3M0rlzg2NGqX70pF3Izly9Qgty7XEM7cnTUs3pW2Ftrxe/nUMWvoXbNkDCcBCZAZdhz171DDz6tUP19zNlcuuCiMA3E28y5yDc5i2dxoOBgfO/ucszg7O7H5nt9XvtfpwRLqTTxyMPMikwEks+2tZainAJHMS20O341vMl3s9B7B80fdoFgsreIMpDOdYwbIUrhLCpFEedPez39J02VrnzmrPbrFiMGkS9OkDeTO2Ej7keghTgqaw+NhicjnlInJEJC4OLsxtM9dKjbY9CcBCZIa+fdXwsocH/Pe/ao63SBFbt+pfou5FMWPfDL458A0xCTE0KtmIBoXfpeGknVy5nWD17EyrD0f8q0rQhyuPAzz1HsGRwdT+vjZuzm50q9aN5SeXk2xOxgkH7h5viP+n0GSrO64M5Ley/QmtoGEqGcZLXsfu/wzWDb4ZeZDI9sLC1Fa5MWNUoO3bVy0X79ABHB0zdOlj147xwZ8fsOHcBnI45uD9Wu8zzHdYth1mfhIJwEJYw507arFJjx6qp9uhA9SsCT17qnRJdkbXdTRNIyg8iC93fUn7Su0ZU38Mkdc97wfIBODZAmRaBGw8/VB6RoD4ZDMBG0//6/q7Lu3im/3fUCR3Eaa3mE6tIrX4oc0PdKzcEecEJ/x3eHHyrx9odzqKD8KTiSwLiZ+Mo0MXGJqa8TFzVmhn5EEiW9u7Vw0zr1ih/ly3rqrBm8Ec5CaLiTuJd/Bw9cBsMXPwykE+b/Q5/X36ky9HPis03D5JABYiIy5dghkz4Pvv1VYiDw8VdF97zdYte6SjV48ycc9EynmUY2yjsbQq34rTg05TLl85AOov3JrmAJkejyu19/fjsUmx/G/r/5ixb0ZqmshOlTtRr1h9KkZ1ZvsrU6l7cCa99Oscd/DmiP9UpizxpZZv1k2pP8uDxHMhPl4tqgoKUhukhw1T5S6LF8/QZWOTYvnh0A9M2zuNV0u8yoL2C6hZpCaXh13Gyfj87wGTACxEephM8PbbKhs/qD0sw4apnM12JmUr0cQ9E/nj3B/kcsrF6HqjATBohtTgC2kLkBlR1N2ViEdcK6UE34oTK3hvzXvcSriV+pqGxqdz/uT8z/WJvOhAKLOILOrDpcEjqDm8EdUcs34hW2Z/TnYhJgZ27YLWrdUy8WrVoGtX6NUrw+Uur9y9wtf7v2ZW8CxiEmJoULwBnat0Tn39RQi+ANl7CZkQWclshn371PcODipl3vDhcPGiqlBkh8EX4KOtH9FofiMOXTnEl42/JGxoGP/X8P8eee7jatFaq0btqOYVcHV8OLm+0fEa3V9R83vl8pWjYcmGTPGbjSOuaBYDxiSdvvPnUL60ie8XuJAr/BTeEeuoPaYxDjYIvpD5n5NNnTunerdeXmoq5fp1dXz2bLWf1wq1pgMCA5iwewL+pfzZ++5edvbeSavyrTJ83exGesBCPE1K4fvp0yE0VKUxLFlSZfexQ4mmRBYfX8wrxV+hfL7ydKrcieJ5iqdpK1Fm16hNGZ79eMNKQuM2ojlFEKsfY09UV95JqseFoOoYf17K4XXL+a5Ica6VPE3da3nw6dSHTl8kQQ4HII9V2pIRz2Ut37NnYeRIWLNGPWB27QpDh2Y4MYyu6+y8tJOAwACG1R2Gf2l/RtcfzcDaAynjUcZKjc+eJAAL8Tg3bsDkyQ8K3/v6woQJqmdgh+4k3kndShR5N5KPX/2YsY3GUrNITWoWqZmma2RFjdrLSSv5yzQUi6MFdGjp2R2n7RPx7Ks+8m7um1ho7kqCQwVcus9WyRzsLFPGc1PLNzFRFf3w8lKf8YED8L//qVSRhQtn6NImi4mVJ1cyOXAyByIPkD9Hfq7fU73pwrkydu3nhRRjEOKf4uJUIoHISChTRs2BDRtml4XvU3y+43OmBE3hduJtGpdqzJj6Y2hauimanST6+Htd1gY/NmB32P39xRYjebYOZ2ygifxVCuI+/gOaN7XgsG2zSsacgco44gmiouC779RWoqpVYfNmddxkUr1fK6g/rz6BlwMp51GOEb4j6FGjh10lc8kqUoxBiKexWGDdOlWI3WBQRRKKFlWVWjzsoHD3I1y+fRkvNy80TSMqLoqmZZoyut5ou6gAk1IKsI5nHY5fP87UoKnMem0uiSeakbz5EyjdBgxJOFtgXdgU6mkaWsP+8DqAQWUJE9Z36pRKDrNwoSqS0KKFerhMkYHgey32GvMOz2NkvZE4Gh0ZWHsgI31H0qZCG4wG49Mv8AKSACxebPfuqXq706apObBixWDIEBWQDQa7DL6Hrhxi4p6JLD+xnO09t9OgRAO+avGV3fR2gy4H0XhBYxJNianZqjyTX6Vbp9zcDoEiRZowI+YtYu/+iF9UDnw7D4DBg9VnL6xP19W/Z6NRze8uXKj2qw8ZonI0Z9CpG6eYGjSVBUcXkGROol6xejQs2ZCu1bpaofHPNwnA4sU2d65aaFK7tqpQZIVMPplB13W2XNzCjVLGLgAAIABJREFUxD0T+fPCn7g5uzGq3qjULUT2EnwBtoVuI8GkEnmgA8H9uLtxCpOrzaf8nHzU710ehyMDYGdVlbLQzf4KpT8X4uNh0SK1ePCDD9Rcer9+0Lu3qoOcRo/L+HUr/ha9fuvF76d/x8XBhV7evRjuO9yuylPaOwnA4sVy9KgagmvaVO3j7d0batWC+vXtrjDC3yWYEuiyogsOBgcmNpnI+7XeJ4+L7VcDp9gfsZ/vD87lDZeZbPulEXg5g8GEweLIHJcEersVx3A4Gq5/AQ4fqS1bdrptK9u7elXN7c6apVa1eXs/yEOeO/czVaH4Z8av8JhYhq9cC7SijXcRbsXf4uNXP2ZgnYEUzGk/ZTSzCwnA4vlnscAff6jAu2WLSg1Zvbp6zc0NXnnFtu17hPjkeOYfnc/qU6tZ13Udro6ubHp7E5UKVHpiTtysyk8cdDmIbaHbcDY6s/z4GvZe3YEh0Z25Pw4gb6Ivb/TcRsWk4by+IRjf0Pmqyv2IEepBR2Su1q3h4MEHiwcbNkz3w2VKxi8LCdwzbuGOw2os2h0m/OFJu5qe7Oi1w65GX7IbCcDi+ffWW7BsGXh6qm1EfftmuFJLZrkVf4tZwbP4at9XXL93nTqedbh27xpFcxd96lairMpPnDLHmzrMHFsAdk+lQe53+W+/v4hpFs5XO+JpsbIwkYVasnnaSJq2bWC1+4u/sVhg/Xq1VW7xYvVAOWOGGmIuV+7p73+KyzFXueOwlrsO67Bod3CyVMA9qRdXEyyAfU19ZEcSgMXz5+pVtcVi6FBwd1fDzO3aqWotdji/m+L4tePUm1eP2KRYmpdpzpj6Y/Ar6ZfmX3KZnZ/4buJd1hw4wncbdpPgkKzy6FkMNHIZwJKO+Siw8FX49ChdL04monBFvmjcBwDX4HuMLx6R/fbI2rN792D+fFVn+swZtY/3zBk1rG+F7XIW3YJBM+DudovLyb/gan4ZN1N7nC1V0NDwfB4yftmBpwZgTdPmAa2A67quV71/zANYApQEQoHOuq7fetw1hMgSx4+r1cyLF0Nyshpm7tDBbgsjAJyMOsnZm2dpU6ENlQtUpk/NPvT07ol3Ye9nvlZm5ScOvXGVIT/PYP31WZjMFgzLVmHs5oROEs4GA1/+PpMCR25AlSp8+cYogvM/nN3ouS5SYAvXr0PFinDrFtSpA7/8Am+8YZWHy6DLQQQEBlAgRwFmt57NZ6+1Y8TKXJiTHyzayvYZv+xIWna5/wT8s3r4B8AWXdfLAVvu/1kI20hIUOXQqleHJUvgvffg9GkVfO1U4OVA2v7alsrfVmbQ+kGYLWaMBiPTWkxLV/AF6+YnDrocxDu/jKbMR20pNaMEv9+cgHOkP/1zb+LHr0pR2Xki+RLeZPUiR8q4VFZz7MePM7dsQ5Ic/h0InqsiBVlo9eEI6k/YSuteXzGp/TBWH45QqSEHD4Y9e1R5wLfeylDwtegWfjv1G6/Me4V68+qxPXQ7Xm4q21u7mp5M6dAET3dXNMDT3ZXxHarJw5SVPLUHrOv6Tk3TSv7jcFvA7/7384HtwBgrtkuIJ4uPV6XRGjcGFxcoVAjGj1fzu3a4dzfFvvB9jNw8kt1hu/Fw9eDjVz9mUJ1BVklUYI38xPHx8PmCnUyMaIFFSwRHC5532/JlowC6N4ni6idjSZh6nP9n77zDqizfOP55OXAYoqKi4gBxizNXCo5wp+VeubXUyr3IykzTSgXEPXJr5kpRs7LMgQvcuLeCorhQEGSe8f7+eBz5y8E4cDj6fK7LKwTO+9yc8L3f537u+/sd98lcHJQSfNHiQ5Jz5WFSgYq0UZTXuh1JUs+mI9fZ7ruIgJBAat04wwP7XDQq+x4AbcaPN9k63wV9x4Q9E3B3cmfG+zP4uOrHOGqfGS60qVpEJtxMIr1nwAVVVb0FoKrqLUVRXtp/rihKf6A/gFsGvSMlkudGLGJi4Pp1KFRIiAtkU1IMKSToEnCyc0Jn1HH94XWmN5vOJ9U+ee5Gl1Eyok989pyRr5ds5Y8YX/RWj8AlBayMaBQNA0vloNcPPeHAARztc7KxyvvY6lNIsrEjxj4X/KvE/EaaFJiDvXup2aYzbR7c4kauAkxs8AlrqzTjETYZLuffT7jPvCPzaFS8EZ6unvR+pzce+T3oUL4D1layLSgryfR3W1XVBcACEFrQmb2e5A0lIgLGjXt2vtuypbACzKBgfGbyKOURC48uJOBAAC3LtGTuB3Op61aXK0OuZNqNLi27leRk+DUwhe83r+ZCPj8ocAYHW1faunfk9zvnSDGkoEWD94+rQFsCZs3C86oL8dr/7maflJjfGJMCc3D9utAhL1cOXF2JdMjD9/V6sa2MJ4Z/VUjSW84Piw4jICSAJceXkKBLwKga8XT1pHie4hTPU9xUP4UkDaT3LnBHUZRCj3e/hYC7pgxKIgHEiEV0tBARUFXYsEEoJw0dCmWyr9rO3fi7zDw4k7mH5xKdFM17xd6jddnWT79u7l3G+oMhzPotiJO/eRNTcAvUm4SLUolvGvxMf/e62MyZT0i+fgTVdsHbrR6eVaOhRQvQaHCavJP415SYZckyjRw6JGbU168XDYNbtoC7O8MGzjJZOX/Qn4OYd2QeGkVDt8rdGOk5kooFKpoiekkGSO+d4DegFzD58X83mywiieSJhN60aeDmJhp83Nzg1i3hUpTNGbtzLAuPLaStR1u+8PqCWkVrmTskDAYxLjp+1e8cK90WrI1YtbNlbLlfqF19K80TXFACAmB1HzAa8fz0Uzx9vhIvLvbsOrLEbEK2boXvv4fgYDG/O3w4DB789MsZea+NqpHtV7fTqHgjNFYaSuQpwSjPUQypNYQiueTDUXYhNWNIqxENV86KotwAxiES7zpFUT4BrgMdMzNIyVvCnTvifHfuXCGhV7WqkItUVaHkk02T79HIo/gG+zKi9ghqFa3FmPpjGOE5grLO5k9Kd+7A4sUwe/UFbhX3h6pLQTGAAoqSgn3R87T4JQkmTBAKYQMHigpD8ReXJGWJOYPExorfY2trCA0VD5XTp8PHH/9HIjI973WyPplVp1bhH+LP2Xtn2dh5I23KtWGE54hM/bEk6UP6AUvMz5MEO2UKfPXVs/Pd+vWzrT6zqqr8c/UffPf7siNsB7lsczG3xVy6Ve6W7muaSkZSVWHfPtGntn496Jp9CtUXorWypXnp5vx9ZSs6QwpajZYdvXbieQPYu1d0kDs5pTt+ySsIDxcKVYsWwcKF0LmzGJ+zsREuRRkkWZ/MtAPTmHlwJrce3aJKwSqM8hpF5wqdsdFkX/GZtwHpByzJfhiN8Pff4uyrTx/o2hU+/RTats3W57sgkm+Tn5uwI2wHhRwLmcQcwRQykvHxokdtyqr9XHVaiMPFT/j883o4NqqIJuc3DCrZhQLL1xPyaxBBuZPwbtYTT1dPcMUk6kmSFxASIo5SNmwQ9pYdOz6zALR7uaZ3aknQJeBg44C1lTVLQpdQsUBFlrdZTuMSjaVMpAUgE7Aka0lMFCND06fDuXPC9N7w+IzLySnb7sASdAmsO7OOnlV6YqVY0c6jHV0qdqF75e7YWttm+PoZkZG8fFlU7RcvSyH2nQng/SMoKvrqq/modxCeroPF+eJP1SExEc8WLfAcORIaNMhw3OkhQZfA8uPL2R+xn/kfzk/XKFZWmU5kCKNRPFzeuQM+PjBokJCMNAHHbx/HP9ifHWE7uDLkCg42Dhzpf4RcttLa0ZKQCViStTRvDrt3i/Pdn3+GTp1AqzV3VC/lfsJ9Zh+azaxDs7ifeJ9iuYvRoHgDBtQcYNJ10ioj+cTgafZs2LotBSvPmdgOmAY2kU+/x2A0EBQeJHa5KSnQrZtIxCYwYc8ILVe3ZGfYTgDae7SnrUfbNL0+q0wn0szDh6LEvHKlKOk7OoozAHd38XEGeXLs4Rfsx/ar23HUOtK/Wn+S9ck42DjI5GuBpEaKUiJJPydOiNJyXJz4+5gxEBQk7NK6d8+2yTc2OZYhW4fgNt2N8bvH4+Xqxd4+e/F29073NZ/IChb/8g/qTN4pZAUfk1oZyZgYUbUvUwY+aJ1EaCh8O9aa4u0X4VmmHAGN/bBXtGiMoNUZ8E58rJEze7Y4ezRD8r14/yKD/hxETFIMAN/U+4btPbZjb21PUHhQmq/3qmqBWQgLE8YfRYvCqFGio/nu48nMihVNknwBjkQeodnKZpy5e4bJjSYTMTyCqc2mksc+ezp7SV6P3AFLTI/RKEYsAgJg507R9dmlC3h7Q5Mm5o7ulUQnRpPHPg8ONg78dfkvOlXoxCjPUVQoUCFD133dru11IyfnzsGsWcIAJ6HEanK0CiBn3jDODAsnr6MjI6J3kXvpahg7g9rGFIJqOOPd8GM863YRF8vi80BVVdkfsR//YH9+u/AbWo2WD8t8yPul3qdBcVH6ruNWh6BrQWm+dmaZTqSLy5ehbFlxvtu5s6gwVK9ukkvHJcex8NhC4pLjGOc9jhqFa7Cx80aal2pukmMPifmRCVhiWmJioHZtYYZgAf67IJLFrvBd+O735cSdE4QNDcPO2o7TA06j1Zhmh/66M94XjZyMbFIW65tFaDoa/vkHbEoEk3fIlyTY7SUesDZaczBiJ809WpHbJidMnAhVquA5cjaeH3wgkoIZSNAl0GhFIw7cOEBe+7x8U/8bBtYcSEHHgs99n18TP3LY5Ejz9c2qN63XQ2Cg6Gr+4gsoVUp0N7dpI37fTUBkXCQzD85k/pH5PEx+yPul3kdVVRRFoU25NiZZQ5I9kAlYknEiI4WYQIcOoomqUSMhG9mhQ7b23zUYDWw4twHf/b4cvXWUgjkKMqTWEAxGkShNlXwhdbu2J4n44UNYuhS+6AxXroj7+uDvTzJLX4dYa3sUvYKKimo0cPzb/jRf+6Eoc547ZzZpzviUeIIjgmlSsgkONg5UzF+R7pW60/ud3uTQvjjJptf1ySxiIE/Od2fOFJKRFSuKUTlrazE7bSKWH19Ovy39MKgG2nu0Z5TXKN4t8q7Jri/JXsgELEk/oaFixGLNGjHL2KQJ5M4Nc+aYO7JUsT9iP53Xd6Z03tIs+HABPar0wM765aMhGem8Tc2u7dIl+HpeCL+dCCLlshelal2h84hb/NxvDDY2lalzahX5j1/mw0vjSUFFa1TxLtFIzJM6OJgl+d55dIc5h+cw5/AcYpNjiRgegYujCwtbLUzV61edWoVRNdK9cvdUr5nlYiAbNkDv3vDoEbz3njgL+OADk8zvqqrKnmt7yOeQj4oFKlKraC36V+/P8NrDKZm35OsvILFopBCHJO2cPCnUkoKCxM7r44+FP2nJ7H3DeJD4gLmH52JUjXz73rdPu0qfyPW9iv8/wwWx60qtN+rLXv9j20rkjinC9Omw5XgI9GwE1kmgiH+XtYvWZl+ffSK+wEBo357fPfIwo2ZJYop3YkzrrmlKPKYa37kZe5MJuyew/MRyUgwptCrbilFeo6jjWidN86dNf27K7Ue3Ofn5yTTHkGmoqpjfdXQUHtOXLsF334kdb7VqJllCb9QTeC4Qv2A/jkQeoVeVXixrs8wk15ZkL6QQhyTjxMfDgwfg6iokC8PDwc9PmCNk09ndJ1yLuca0A9NYdGwR8bp4OlXo9PRMrWnJpqm6RkbmdOG/uzaXHDmorq/CuN55OHUK8ueHSp/P4ZTVs11yn7IfsfhsKZR582HgQDYXrcaO9mP4o8S7wh0niTSN32R0fEdVVR6lPCKnbU4MqoFVp1fR+53eDK89PN2ym97u3ozZOYaohCicHZzTdQ2T8eR8NyAADh4UjYOrVkHp0mK0yEQsO76MCbsnEBYTRqm8pZj3wTx6VellsutLLAeZgCWv5sYNMcKyYAHUqSOcWkqWFIeTZmrySQs/HfmJgX8ORFEUulbqyijPUVQqWCnN1zFF522bqkWoXagI8+bBvFlw4B6UqXOO73+CkT092HKlNV02rAFVRWtU6DdmA0qYXjSxAb47r3Kz1POKVWl5CEjvQ4TBaGDT+U34BfuR2y43f3f/G7fcbtweeful57up5clY1+7w3bQv3z5D18oQixaJJrbr10Vj1Zw50Mt0SfFe/D3yOeTDSrHi8oPLuDi6MLXpVFqVbfXa6ovkzSX730El5uHYMSEPWby42Ok2bCh0mp+QTZOvqqrsCtvFuXvnADHqMrTWUK4OucryNsvTlXwh9XO6L+PUKVGpd/UKYcKuSbh8+BNes1pzsUl5Qp3GYmcHHSt0ZK+xFxO3G9nxiwbPZn3h/HmYPx/I+ENAWl+foEtg7uG5lJ1dlg6/diAqIYrWZVvz5Ngqo8kXoGbhmjjYOKRrHjjDhIc/U2G7cUP8rm/eLN7zAQNEpSeDXLp/ic9+/wy36W78fvF3AMZ7jyf4k2DaerSVyfctR+6AJc94cjPSaOCPP+D334U92uDBL3XHyS4YjAYCzwXiG+zLkcgj9KvWjwUtF1CxQEWmNpv62te/7mw0PZ23qgrbtsHUqWKMyLZUCMaeDUBJ5hSQKy4X4+uOZeD1gmLn5eaGZ8NeeFoXh88+A+fnS7IZHb9J6+vnHp6Lzz8+1CpSiymNp9CmXBuTJwwbjQ113eoSERth0uu+kpAQUWYODBRKVW3bwtixMH686ZaICMEv2I9N5zdho7GhZ+WelM8vRFDM7QctyT7IJiyJUKlasgRmzBCORB07Cts0EKo+2Zzlx5czcc9ErkRfoXTe0vh4+by2o/nfpLbBKrUNTElJ4ugwIADOnAGXIskMGmBNUg1ffggZg4qKFQrfahoybt45Mcb1ww/w9dcmiTO9r794/yIBIQE0Kt6IjhU6Ep0YzZl7Z9LcWJVWkvXJmS8sYTA8O989cED0LXz2mXi4LFzYtEsZDZSaVYqHSQ8ZUHMAg94dhIujecbDJOZHNmFJXsy1a88s0mJjxRlv/vzia9k88T5IfEAeuzwoisKZe2dwdnDGr4lfus7UUns2+m/BjBcRFSUsAGfPFkqEFWvE0GnGT+zVzaB8izm4OHoz9bAdKboktHqVpot2QPnGwrC3WbPXxpnR8ZuXvT5/vnDarh3E5vOb0Wq0FHcS1Y489nmo61Y3VdfOCJmafA2GZ+NCX30lFMFmzxbnuyaSiEzSJ7Hy5EpWnlzJX93/ws7ajk2dN1Eyb8l0GU1I3h7kDvhtRVWFmMCFC2LHO3w4vJv9B/7DY8IJCAlgcehiNnTawPul3ifFkIKNlU26d2nFv/yDF/0rUICwyR+89vVXrohx6CVLIDFfCCVa/IZ75escjt1CXEocTUo0YWKx3tSq35WQiBCCAgbjHV8Az4GToEqVdMVsKnpv6s3yE8vJa5+XATXEbu3/Fauygu6B3XF3cuf7ht+b5oJPHi43bYLTp8HeXpz5urqaZH4XhGzpvCPzmHlwJnfi71DVpSq/dvxVzu9KnkPugCVixGLDBpElAgNFg8miRUJA3tXV3NG9ltBbofgF+7HuzDqsFCu6Ve5GiTwlgIwrVqX3bPXgQfD3F2+ntTU0/TiE7UUacdWQyNUoaFK8MVOsmlJ17kYI6Qah5fF8xxPPgMNZrs38hERdIj+f/JkuFbuQ0zYnrcq2okbhGvR5p49JmqrSS1RCFMdvH894Aj54UJSZN2wQf+/USVR37O2FK5GJuPLgClXmVyFeF0+zks3w8fKhYfGG0oNXkiZkAn7TiYkRLjizZkFEhBixCAsTu18LMWE3GA20WtOKh0kPGV57OENrD6VoLtP4qkLaGqyMRtGb5u8vHOdy5VbpMnofCRXm807RsmzdnQKABoUG649Sdct2KFFCvP+lSomLmOEmHZUQxZxDc5h9eDZRCVHYWdvRs0pP2nm0y/JYXoS3uzdf7fiKu/F3KZCjQPouEhoqdMhz5xaiGYMHm/ThMvRWKKfvnqZHlR6UyFOC4bWH07FCRyoXrGyyNSRvFzIBv8lcvy7s5+LjhRPRnDlCQi+bjhA9QW/U8+uZX1l5aiUbO29Eq9GyodMGyuQrg5Od6UU/UnO2mpwstBj8/cWUiqubkT6+mzmd25dfbh3AOdKZluUaoNVoSTGkoE0x4K0vChsWQevWJit7phWdQcfQv4ay7PgyEvWJfFjmQ3y8fKjnVs8s8byMJ/PAe67toUP5Dql7UVycEM2OjYVvvoF33hEe061bQ86cJolLVVW2XdmGX7AfO8J2UMixEB9V/AgbjQ0TG040yRqStxeZgN8kVFWY3Z8/Lzo83dzAx0fckN5Jn/B9VhKfEs+S0CUEHAggPCaccs7luP7wOqXylsp0QfqXNVg9fAg//QTTp8OtW1C6QQhtPv+Do/oVLI2LoLi2OHMqf0XvzddIXjabaV2nEJFyhErGUtz54WMwk0F8eEw47k7u2GhsuPzgMl0qdmGk18inozDZjeqFqpPDJgdB4UGvT8AREaKisGCB+B/UpIn43VcU4TFtIvZf38+APwdw8s5JCucsjG9jX/pX74+NJvsajEgsC5mA3wRSUoQhwrRpcPy4KLt9/LEwux83ztzRpYqr0VepubAmDxIfUMe1DjPen8GHZT7ESjHPbv3WLTGVNW+e2GDVbxZN60nLWX7ja67GiDLzhCI9+GrtTay3T0Jv78D6ik1IflCE3DYluE7aZCJNgVE18sfFP/AL9uPgzYOEDw2nUM5C/NX9L7O9j6nFRmND/+r9cXdyf/U3Ll4sHi5VVbhtDR8OtWqZLI645DjiUuIonLMwue1yY1SNLG29lK6VuprUHUsiAdkFbfn8/Tf06SMyRvnyMGyY2AXYZ4E3aga5/OAyp+6coq1HW1RVZfjfw+lYviN13OqYLaYLF0SZecUK0bfW4qMIcr8/nc03FpCsT8aoGjGoBjRYMXG7ka+uFoYhQ3g/oRznk//7PFvEyZ79XzbM1JiT9En8cvIX/EP8OR91HrfcbgyvPZx+1fqZtbHKJBgM4tC9eHFhjHDunGgeHDIEihUz2TL/9uBtXro5q9uvBniqGS6RpBfZBf2mcfasOMctV+7ZjWnpUmja1GzdtWnh0M1D+O73JfBcIM4OznxQ5gO0Gi3T35+eaWu+TkTjyBGYPFl0NGu10KFfGInvjmfLtVWoYSoflW5L04tGPlO2kKKIzmvv3qOhx1eg1XLhyz9euG5atKLTS2RcJP1/70/lgpX5pd0vdCzf0WLLpEn6JB6lPMJZtYdly0Tt//Jl+PRTIcnp4SGkxUzE2Xtn8Q/2Z+XJlU89eEfUHvH06zL5SjITmYAthSe6htOmiV1vp06wdi2UKQN//WXu6FLF0cijjNw2kt3XdpPbNjdf1v2Swe8OzvTS3stcgFQVcscUYdIk2L4dHMoFU3fMP4zu0JQShZ14d1EgA0t1Zfg+A8UmB0JiIqUHtCToI0+83b3xdH3WRZ5Rmci0cC3mGtMPTOd2/G1Wt19NiTwlCP00lEoFKll0wjCqRtymudExuSRzAi5AdLSYTV+7FtqZrlv7SdVPURSWHV/GmtNrpAevxCzIBGwJrF4N338vdr4uLsK15dNPzR1VqkgxpBCXHEc+h3woisLV6KsENA2gb7W+5LQ1Tafq6/h/pStVhftnnOm5zJFHN6BgIQPNvp/CNv1Y9mLkyO9T2NFzB7dSBuPYbbLYEvfoAcOH41m+PC8a3kqPVnRa+fcstKIodKvUDYPRgMZKY/mjMKdPY1W+PNULVyfo0hFo0ECMEnl5mayq80Qv3C/YjwkNJvB+qfcZXWc0X9T5wvxWiJK3EpmAsyu3bkGBAmJ85exZsLUVB5OdOomPszkPkx6y4OgCph+cTuMSjVneZjnVClUjbGhYljvAPCkDqwaF+LOFiT1YEt39nGjy3adbwAIOavz5O/rS0+9PMaQQFB6EZ91GMNZGOOMUfLU6VEZlIl/H8uPL6b25Nzm1ORlWexhDaw3FNXf2F1B5JUajqN4EBMCOHbBlC97FvPny8l/cXTEv/fPA/0eCLoFlx5cREBLAlegrlMxTEr1RD0A+h3wmWUMiSQ+yCSu7ceyYKDOvXQu//ipGiFJSwMbGIs53b8beZMbBGfx09Cdik2NpVLwRo+uMpknJJmaLqfbEIC7tdSbm+j2M+Q+gialJnuL5sKs5gxvG5VQv8A5tY1z4IeFvUhQVrcaGHR/vfq7EnNXoDDrWnlmLi6MLjUs05n7CfRaHLqZ/9f6ZMgudpeh04nw3IECMzBUWjWz078+hhEvUWlSLtR3W0qlCpwwvpaoq1RdUJ/R2KLWK1MLHyydTXJ0kkpfxqiYsVFXNsj/Vq1dXJS/AYFDVjRtVtX59VQVVzZFDVQcPVtWrV80dWZoZtnWYavWdlfrR+o/Uo5FHzRpLbKyq+vqqqlM+vUrRYJVv7FTGoTLOWi32TYC6NPiYuvOrj1RjTkdVBTW45Tvqjwt6qMHX9pkv5qRYdWrwVNU1wFVlPGrXDV3NFovJSU4W/9XpVLVYMVWtWlVVV6589nlVVXUGner4o6P6+e+fp3uZy/cvq6P/Ga0m68V1N57bqO4J36MajcaMRC+RpAvgiPqSnCh3wObEaBTdzAaDaKbS64V8Xt++wi4tm6OqKnuu7cEv2I8RniNoWLwhtx/dJlGXSPE85vMPjo4WOg3Tp4uPa7c6Rfi7XbmtP/04cOhe/kt+7jRJnO0aDOK8scaLH1KzipkHZ/Ltrm95mPyQ94q9h4+XD81LN8/2M7yv5exZsdvdvl3MednaiiMWF5cXVnVWn1pNOedyVC1UNU3LHL55GL9gPzac24BG0bCr1y6zjrRJJCDHkLIf166JDPHbb3DyJNjZic5md3eh6p/NMRgNbDy/Ed/9vhyOPEx+h/zcjb8LkCm+p6n14b1zR1Tv584VKoWtWoFVq0/ZdGMBdqod1lihGo1oDTDApaJ40fLlZpXmPHfvHMWciuFg44CDjQNNSzZllNeoTFf+ynRUFXbuFCNDW7eK3/E3tZBzAAAgAElEQVTevYUsqq0tFCr00pd2qdQlTUvdT7hP+3Xtn3bX+3j5MKTWEArnNK3Pr0RiarL/3f5NQVUhJERkiMBA8eTfoYMwS3BxeSbUbwF4L/dm3/V9lMpbinkfzKNXlV7Y22SO8MfLRojgWePTzZvg5wfztoSQUnQnlXtoWNx3ODWq2rL0UDVqRLfi8+VnuBBzhaB3cuNdvyee77QUC5gh+aqqyr7r+/AN9uX3i78zu/lsBr47kL7V+tK3Wt8sjydTCAmBxo1F89rEiUK9yjl1ncY6g47x237h99BkYmOLvvChK8WQwsk7J6lRuAZ57fPiqHVkatOp9KvWL8u66yWSjCJL0FnFkSNQs6YoLffvD4MGWYQNIAgnnWXHlzGs9jCsraz55eQv2FnbZUkzS53JO184X1vEyZ5VXRoyZYpQJ9QXDYIeTTEqOgDG1R/L+AYT4MED8T6XLw8jR4qHHjNVGVRVfToGc/DmQZwdnBlUcxAD3x1o+WMwDx4IoQyDAcaOFQ+cGzZAy5Zp7trfcPQanbZ44GDwJp9uICBGuia1q0QDD0d+OvoTMw7OIC45jhsjbpDLNldm/EQSiUmQJWhz8OCBsAFMToZvv4Xq1eGXX0RXcw7LkAe88uAKASEBLD2+lER9ItULVadB8QZ0q9wty2J4kZKULtqBk1tLUmosoNFRdbAvJ3L/SIpRJF9FBd28xeD9HeTNC6dOCcUwM3WRG1UjVooViqIwNWQq9xLuMafFHHq/0xsHGwezxGQyLl8WVZ1lyyAhQQhmPDFG6JBKV6P/I+CfK9gaK5Bkderp5x7pohj4xwji/tpKXEocjYo34os6X5BTK3e7EsslQwlYUZThQF9ABU4BfVRVTTJFYBbLhQtCxX/5cnFDatXq2Q2pa1dzR5cqYpJi6LelH4HnArG2sqZ7pe5mc9L5t8KU7n4OHh4oRfyZwmAXx6DPwMfHmgYb1pDrTj4eahMxoqI1gjband8OXKGVZynhx2sG7ifcZ+7huSwKXcThfocpkKMA6zutp2COgm/GGMysWTB0qBiR69ZNGCNUqpThy0bGJGJrXYlEmyPoicIaZwzKA27p1/ORRyd8vHyoVqiaCX4AicS8pPsATFGUIsAQoIaqqhUBDfCRqQKzSGbOFPrMS5bARx+JBqvNmy1ifteoGrl0X4hR5LLNxc3Ym3zh9QXhQ8NZ3HpxhpLvptCb1Jm8k+Jf/kGdyTvZFHoz1a/1aVYWq4e5uB1sIPJ8MPE229H2bY/D18WY4BuNm5tCxyMfcm9KBJvW5KBp2DsUTRjP0upfM2X39XTHnBHCosMY/Odg3Ka78W3Qt1QsUJHY5FgACucsbLnJV68X8+mnH3eTN2gAY8aIpsIlS0ySfEE8dNkZxLXu2H4NgFYtSQ3btaxuv1omX8kbQ0ZL0NaAvaIoOsABiMx4SBZEYqIoK1evDlWrQqNGMH48fP65ULGyAJL1yaw6tQr/EH9uxd3i+vDrOGod2f/xfpPoCqemieplnD0La32LELY/HHp1Ak0yKICi5VPexbhmNfQZwMZC76JrMZzfPeqRbP1MVzorjBD+n9uPblNmdhkUFLpV7sYoz1FUKFAhy+MwKbGxwoFoxgy4fl04bk2bBhUrij8mxqdZWb4MTCRO3xgbtSggzoDHvG8+YRSJJDNIdwJWVfWmoij+wHUgEdimquq2//8+RVH6A/0B3Nzc0rtc9uLWLTHrMn8+REUJ0/uqVaFCBfHHAniY9PBpM0tkXCSVC1ZmZvOZ2GpEw4ypRP3/X4cZIFFnwO/vCy9NwGfOiMbZdevEcXnNYRs4bJ0s4lJhdLDKhG374JOy0AecnXOzoVKj/1wnM4wQ/h9VVdl2ZRuHbh5i7HtjcXF0YWHLhTQp0YQiubLGBzhTmTBBjBLFxsJ778Hs2fDBB5m65DNZz68yRdZTIskupLsLWlGUPMAGoDMQA/wKrFdVdeXLXvNGdEEPGyaSr14PH34ozr28vS2izAzP/E0P3jhI7cW1aVyiMT5ePjQp0SRTnHSKf/kHL/oNU4Cwyc/fyE+fFvf7XzclYvvuMqo1vM6WoZO4GB+M99L6GAwGtEbYcbc5ngMnQZUqwH932fCsazazbto6g441p9fgH+LPyTsncc3lytmBZ3HUOmbKelnK8ePivVUU+PprCA/PFkIlEoklklld0I2BMFVV7z1eJBDwAl6agC0Sg0HYADZrJmZG8+UTTkRDhkDp0uaOLtWcvHMSv2A/cmpzMveDudQqWotzA89Rzrlcpq6bGpu+M2dgiG8IOyP+xDrPbRy+3kyCcg8rh3LkThmEp6snQeX9CDqxGe/2I/Cs2uq5a2W2EcL/s//6frps6EJEbAQV8ldgWetldKnUJdNtFTMVoxG2bBG73b17hUlCs2bwww8W83ApkVgaGUnA14HaiqI4IErQjQAL397+i9hYYXI/cyZcvfrshjR2rLkjSzWqqrIjbAd+wX5su7KNHDY5GPTuoKdfz+zkC6+26Tt3Dr77DtYGB0OvBlA8Bb0CFfSuzNzuQr0D51HUX+CLL/DsOBzPjsNfuk6bqkUytUR5+9Ft7ifcp0KBCpTKW4pyzuWY98E8y5eKTEkRg9TTpsGlS+DmJmQjPR+ft8rkK5FkGhk5Az6oKMp64BigB0KBBaYKzGzExYm53cWLxcdeXjBlimiwsjB+3Psj3+z6BhdHF35s+COf1fiMPPZ50nSN1MpAvowX7U67lq3Ar/4FWbXzBFpbAw71A0nQ6EEBKyN0DoqgQkJ1lPWzoU2bNMVrai5EXcA/2J8VJ1dQq0gt9vTZQ0HHgmzr8Z92B8tCpxPjQyB2uYULw5o10L69RcihSiRvAlIJC8ScbmQkFCkiznbLloXatcV5b82a5o4u1cQlx7Ho2CK8XL2oVbQWV6OvsitsF90rd8fWOu0ewqY+W710CSZMVPkleCdKXV+MxbdRMLk6Ch9xRzsGRU3B2mhF7dufoy/Zjv1fNkzzGqbiSOQRftj7A5vPb8bW2pY+7/RhhOcISuW1HMnQF3L6tNjt7t4t2sy12lcaI0gkkowhlbBeRkqKmGucPl0k4GvXxA3p7FmLML1/wq24W8w8OJP5R+cTkxTDV3W/olbRWpTIU4ISedIvQpGeDuYXERYGg6aEsPXufNSiB6DHRfJr8jD8Sin6rzlKh+79IP8PJCsnsVUrcy2fB4oZRoiMqhGjasTaypqQiBD2XNvD2PpjGfjuQJOZw5sFVRVORFOnCtMPe3thjJCQIH7fX2GMIJFIMo+3MwFHRYkRojlz4PZtIZ4xfry4UYFFJd8v/vmCGQdnoDfqaefRjlGeo6hVtJZJrv2yOdrUztfeuAHjvk9g2bajGLs2A5ckFEVl9Nl8jAu8j10Be+Z69yUqhxO2xqLY4vH0tVkxQvSEZH0yK0+uxD/EnxG1R9Cvej/6VuvLx1U/JofWMmRDX8mePdC0qdjlfv+9MEbIl8/cUUkkbz1vVwI2GECjgRMnRDNVs2ai0appU7Na0qUFVVXZH7Efz6KeaKw0ODs407dqX0Z4jqBk3pImXSs1Hcwv4vZt+HbyfZacmoOh+ixKdvYmXJuCQVWxMkAu1Qa7JSugc2cKn7mHLvAUvKBJK7OJSYrhpyNiFvrWo1tUdan6dHY3s9ydsoToaPGAaWUFo0dD/fpiqLpVK4t6uJRI3nTe/ARsNIoO5unThWpPQAA0bCg0m8uUMXd0qUZv1BN4LhD/YH8ORx5mQ6cNtPNoxxd1vsi0NV/Vwfwi7t+HMX7hLDoTgKHyYqifQPNHxehw7TSDKmhJMaSg1drg7b8B3LyArB8h+jdt17YlKDyIJiWasKLtChoVb5Qps9BZxpUr4vd8yRJRXu7YUXxeUZ59LJFIsg1vbgJ+9AhWrBDyeRcvigarlo89YBXFYpJviiGFn478xLQD0wiLCaN03tLM/2A+zUs1z/S1U5scHz6E4QEh/LI/iBSPpWiqhdH9biG+2pBA+Ye3oEcPPD6aR1BkCN7u3ni6ev5nnaxIuGfunmH6gelMbjyZfA75+LHhj9hZ21G1UNVMXzvTmTFDiMJYWwvTjxEjoHJlc0clkUhewZvbBf3pp7BggehiHj5cWKM9GbuwAFIMKWg1WgxGA2Vnl6VAjgL4ePnQqmyrbCPmHx+vMnTGDpaFTcRQ8DBYp6BVrFi3WkfrqHwwYAAMHChM2c2Eqqrsvb4X3/2+/HHpDxxsHAjsFEizUs3MFpNJMBhg40ZR1SlXTqhXrVsnfKYLFzZ3dBKJ5DGv6oJ+MxKwqsKBA6L8Nno0VKsmdr337ok5XgsqK56POk9ASABbL2/lwqALONg4EJUQla0M2xOS9Aycu56VV33R5w/FVqdFZ6PDiIpG0TDRsSVfffYLOJjX6zZRl0jDFQ05cOMAzg7ODHl3CANqDiCfgwU3ID16JPoWpk0T7eUjRojuZolEki15c8eQdDpYv14k3kOHwMlJCAlUqyZKzBZSZlZVlX3X9+Ef4s9vF37DztqOXlV6kahLxMHGwezJNyQihKDwIOq5eXNyRzmGnq2OPmcYhW3yMG6rLWVvJdO8t4YUDWg1Wrw7fmG25JukT2Lf9X00LtEYext7qrpUpWflnvR6pxcONuZ9IMgwP/wA/v4QEyOUqvz9oXVrc0clkUjSieUmYKNROBCdOSMS7Zw50KuXsM+xMI7dOkb9ZfXJZ5+Pb+t/m63mTkMiQmi4oiHJ+hTQ26Iu20HDWi4MOnONVpceounYiaDBPSh17QoRCUdwdajBnSg3cM3aOKMTo5l3ZB4zDs4gKiGK8KHhuOZ2Ze4Hc7M2EFNz7pwoMSuK6HJr3BhGjhRCMRKJxKKx7BL0okXivOv99y1mjAggQZfAsuPLiE6MZkz9Maiqyroz62hZtmW22qVdeXCVlku7cS7ugLAvMmroVmgiPxeqgLJnNwwdyqb7mix3Ivo3d+PvMmXfFBYcW8CjlEc0K9mM0XVG4+3ubbkdzaoqDECmToV//oEdO0Tnvqpa1HGKRCJ5G86ALYS78XeZfWg2cw/P5X7ifRq4N2B7z+3ZTsz/avRVPv31a7ZH/opiBCvFCIDWyoYdH+9+rou5zuSdL5wVLuJkn6lSksn6ZGytbbkZe5Mys8s8FSGp4lIl09bMdHQ6WLlSjMqdPi0UqoYMEQ2FedKm4S2RSLIHb+4ZsAWx8uRK+v7WlxRDCq3KtsLHywcvV69ssUsLiQhhV/guahSugWtKU4Z/B/vcNzP4mJbRB5IIr1WOPW2q4P3BwP+MEGVULSstqKrKnmt78A32JVmfzPae2ymSqwg3R9zEyc7J5OtlGU8EYoxG+OorKFAAli2DLl2EVKREInkjkQk4k3iiWJXHLg8VClSgRuEa9KrSixGeIyjrnPkqT6nlx+2BfLO/M6qqB9UKZdleHKO9OOjaGA93K2x+H0mRevWo85IHhfSqZaUFg9HA5gub8d3vy8GbB8nvkJ8htYZgVI1YKVaWm3yvXBHdzLt2CXU2W1s4eFBYAmaDBzOJRJK5ZK/a5xuAwWhg/dn1eC72pN7SevgG+wLCe/enlj9lm+QbnxJP3w0TGLu3FyqPrQAxUrDWUmZvukXlExux+WOzkDF8RTLwaVYWe5vn55JNLSW58NhC2q9rT1RCFPM+mMe1Ydf4pv432a50n2qCg0W3funSYlb93XeF9SVAsWIy+UokbwlyB2xClh1fxsQ9E7kafZWSeUoyp8UcelXpZe6wXkjAnhksPj2O0ndtueYMBgWs0ODg6sZPh87Rs2HqHHIyQ0oyOjGa+UfmUzpfaTqU70C3St3IZ5+Pdh7tso0ISbrZvRu8vcWZ7pdfSuEMieQtRibgDHI3/i7ODs5YKVZcun+JAjkK4NfEj9ZlW2ebZBESEcKGcxsIiw6jS4VuRIe0Y+3EduxkHAXvFGNCk5rsqAzWvIPBxiPN57emkpK8EXuDaSHTnnY0D6w5kA7lO5DTNicdK1iolnF8vBDOMBpFQ1W9erB4MXTqBI6O5o5OIpGYEdkFnU6eKFatOLGCXzv+SsuyLdEZdNhospfc5dLQpfTb0g+DagAVWh1347fN1/D0hIJF1hJaIsd/Sp6Z3cH8IibunsjEPRMxqkY+qvgRPl4+lt3RfOsWzJolXImio6FFC/jjD3NHJZFIshjZBW0inugK+wf7s+XiFuys7ej9Tm/K5y8PkO2Sb/8t/Vl4bCGogAIaI+R6lJfN65Jp2cGWzcfrcj4Nbkem5EmTWqUClchtl5tyzuX4rMZnjPAcgbuTe6avn6nMny92u3o9tGkDo0YJSVSJRCL5FxbaxWIeDKqBnht7EnIjhHHvjeP6sOvM/3C+yX1404vOoGPVqVXEp8QDUO2MHQMPgp0erAwKGo0d/RfNpVVHWxRFlI4ntatEESd7FMTON7MFNIyqkU3nN1FnSR3qLa3HomOLAOhYoSMzm8+0zOSrqrB9u+hqBiGF2q+f0CMPDJTJVyKRvBBZgn4Fj1IesTR0KatPr2Znr53YWdtx6s4pSuYtmW0Uq0IiQvj7yt/EJMWw8Wwg1+MimOU0gNPn5rBmYRx9bZaTNNqD/O8domnp/1oBZhWqqrIkdAl+wX5cuH8Bdyd3RnmOok/VPtnmvUwzOh2sXSs0mU+cgKFDhS65RCKRPEaWoNPIrbhbzDo0i3lH5hGTFIOXqxe3H93G3cmdSgUrmWydTaE3M9Q9vPf6Xhoub4jeqAegyj0Nc/6BB1cfsRj4fGBORo8dRP78AI1MFndaeGKrqCgKa86swd7GntXtV9OhfAesrSz412/WLJgyBW7ehPLlRWNV167mjkoikVgQFnwHzBwuRF2g8vzK6Aw62nq0ZaTnSLxcTV9C3BR68zkN5ZsxiXwVeArgtUn4QeID8trnZd+1fU+Tr5URqpwpxw8XF1K0gyfnJkGpUo+T/GLTjQillltxt5hxcAaLQxdzrP8xXHO7sq7DOpzsnLKF+le6iIwU8pCKIkwSypaFhQuhWTOL0iKXSCTZg7c+Aauqyq7wXVyNvkrfan0pk68M494bR6cKnSiVt1Smrev394Xnmp8AEnUG/P6+8NIEefjmYfyC/fjj/G9c/fgE7xXzxlbVolN1qAZbQh0WsiDE86lRTkaSfHq5EHUB/2B/Vpxcgd6op0P5Dk8fEvLYW6ie8dGjwhhh3ToICoK6dWHmTLB+6//5SCSSDPDW3kF0Bh2/nv0V/2B/Qm+HUiJPCfq80weNlYav632d6eunRkN5U+hNxv4ZSFjSZhTrazziMrlTrBhy0Mjdh7/yzbFvSL4UhHONIEZ28Gb0RM/nJorSk+QzQlRCFJXmVUJjpeGTqp8w0nNktmlQSzNGI/z5p0i8QUGQMycMHw4lSoivy+QrkUgyyFt5F/n78t/029KPiNgIyjmXY2HLhXSv3D1LhTNep6G8KfQmwwLXcU0zGjQ6UGHQQRgZVoa/nL+m+tTO5HaG2eM96d/fE5sXTEBltlGCqqpsvbyVvdf2MqnxJJwdnPml3S+85/5etvEzTjfJydCnD9jbiyTcty/kymXuqCQSyRvEW5OAb8TeQG/U4+7kTpFcRSiZtyRzP5hLi9ItzKIp7NOs7At9dAc1Ksr0A9OZ/udfPFQLAgah02yE7Q4dmHt5NTbh1owYLYxzcud++RqZZZSgM+hYe2Ytvvt9OXX3FK65XBlddzROdk6Wq1gVFQXz5sHffwu5SHt72LkTypXjhU83EolEkkHe+AR8/PZxpoZMZc3pNXQs35FV7VdRsUBFdvXaZda4npSAv90aSETCEVzsy1Cx+B0GbltGtD6O+hFwP+doHua1RjUaMBq0nD80ghxl73D6jyK4u79+jZcl+YwIbRyJPEL7de25/vA6FfJXYHmb5XxU8SO0Ggu1zbt0STgSLVsGiYnQvDk8eAD580Ml03W8SyQSyf/zxibgHVd3MGnfJHaE7cBR68igmoMYWnuoucN6joLO17msjibZJpkYvZHzF6HdOfA5nYvzBVsxOWcjrHY2xpDnMDYJVcnX0EAJjwu4u6fu/NZURglRCVFExkVSuWBlSuUthYezB3NazDFb9cBkBAeLhiobG+jeHUaMgAoVzB2VRCJ5S3ijEnCyPvnpzOlfl//iXNQ5pjSeQv/q/bOdZ+yhm4dYcXIFKYYUjBhRVBh+zomp3j9y1qcXM4dacWGXHdZ54slfxgn70ndw0D7Cp1nadmUZMUoIjwlnavBUFocuppxzOY72P4qTnRN/df8rXdczOwYDbN4MsbHQuzfUqgWTJkGvXuDiYu7oJBLJW8YboYT1IPEB84/MZ9ahWSxtvZT3S71PbHIsdtZ22ao0+qRpyXfnRHbfPkCt+DyczJ0kxCoUa9a32sGW+XVYsED0+7T9OIZzeUK5/SghS2d4z9w9w4/7fmTt6bVYKVZ0r9ydUV6jnmpeWxzx8aLEPG2akIusVQtCQqTvrkQiyXTeWCWsq9FXmX5gOotDF5OgS6Bpyabktc8LQC7b7NGxGhIRQlB4EFaKFSsPLeR03BWKPoSAgwp93Rtx+rPBbL++nzsHvOlS35P4eBg4EL79FpydnYAGWRKnqqoYVSMaKw0Hbx7ktwu/Maz2MIbVHkbRXEWzJIZMYfVq4bn74IFIvJMnQ9u2MvlKJBKzY7EJ2Kgaabi8IZFxkXSt1JURniOoXLCyucN6ju1Xt9NqdStSDClYqeD6wMCKA3Z8VPdTbFaPQHV1IzIQlvrUJywMPvwQ/PxE421WYTAa2HR+E77BvnSp2IVhtYfRvXJ32pZra7nCGefOibndokXB1RXq13/mSCQTr0QiySZYbAK2UqxY0XYFpfKWonDOwuYO5zluP7rNrP0BBByYTjJ6VFRQNHzi0pwe21ZD7tyEhsKwHrBnj2i2/ecfaNw462JM1iez4sQK/EP8uXj/IiXylMDFUZyDajVatPbZp3SfKlRVvJn+/vD772LXO2uWaLKqW9fc0UkkEsl/yFACVhTFCVgEVES4zn6sqmqIKQJLDfWL1U/3azNqhPAiLt6/yNSd37P8zCpSMOAdBsHuVug1Vmg1Whp0H8udpNyMGQlLlkC+fMI6tm9f0GSdBggAXQO7EngukGqFqrG2w1rae7TPUiESk7JhgygtHzkixoe++w4+/9zcUUkkEskryegOeAbwl6qqHRRF0QIW4StnCo3kTaE3n87wujrUYELzdkzZ0phQXQR9jsMI+4aUHvgtIcVtCLq2G68i3uxZ5UmzHyApSUy8fPMNOGVRc/atuFvMPDiTobWH4uLogo+XD5/X+JxGxRtZpjlCQgI4PP51++sv0dn800/Qo4cQ0ZBIJJJsTrq7oBVFyQWcAEqoqbxIdvEDrjN55wsVooo42bP/y4avff2m0JsMDVzDdc2XgB4FG9wMU5hkd4+GkTcpOHQMlCkDiMroxo3g4wNXr0LLlkLZsHRpU/9UL+bS/Uv4Bfux/MRy9EY9K9qsoFvlblmzeGYQGSmMEObPF6pVtWqJ5OvoKB2JJBJJtiOzuqBLAPeApYqiVAGOAkNVVY3/v8X7A/0B3NzcMrCc6ciIRnKKIYXRW6YRZZwH1sLlB1XHQ/U4szW96DLnWQI/eRKGDYNdu4S+w7Zt0KSJSX6E12IwGugW2I11Z9ah1Wjp804fRnmNylSHp0zl9Gnx5PLLL2Ket337Z9rMUqNZIpFYIBlJwNZANWCwqqoHFUWZAXwJjP33N6mqugBYAGIHnIH1npLR89v0aiTHJT6kgm8xInhIiUcQkRv0Vgpgg52x0tMEfv++GCOaP1+UmGfPhk8/zXwDHVVVOX77OFULVUVjpSGXbS6+rPslQ2sNpaBjwcxdPDNJShKdzMnJ8Nln4qnmiSuRRCKRWCgZSQk3gBuqqh58/Pf1iAScqZji/DY1GslPznivJezFWaMwtU0AbaoWofe1POS77Myxgt3Z5pCDJM1p7IyVsDV6UCinA7Nni+QbGwsDBoh+oLx5TfgGvAC9Uc+GsxuYsn8KobdDOTPgDOXzl2dBywWZu3BmodPBr7/Cpk2wZg3Y2UFgIFSunPlvpkQikWQR6U7AqqreVhQlQlGUsqqqXgAaAWdNF9qLMYXH7es0kjeF3mRA4Gxua3xRrY3EAt+sLAqMYMKcs2w6/4D5gaew1Rmw1Qt1KOON/Nw8WJXBl6FhQ5gxAypWNNmP/UKS9EksDV2Kf4g/V6OvUjZfWRa1XETJPBbqwRsbC4sWwfTpEBEhBqIjI8U8r7e3uaOTSCQSk5LRouhg4JfHHdBXgT4ZD+nVmMrj9mUayVejLjF6VStu5TgvPqGAYoQHmhuPk3zD5xL4tXBI3FeRB2cK4O4uJmIyW2hJVVUURWHdkUsM+ns41sZilNV+x49efWhXzTXzFs5MTpwQZebYWJFs580TzkSysUoikbyhZCgBq6p6HHhhd1dmkRket0bVyP2E++TPkR/7ezHEqefpeMaOjR469FYqKNao1jWfS/JNyxbh2MYi+C4TM7zffw8jR4pqaWYRGRfJtJBpnLhzgs8rLGHSHzdx0c/GWnUhKVlhzMYzWClWWaIXbRJOnoRr10RreIUKwpGoTx+okaW/UhKJRGIWLG574dOsLPY2zwtGpNfjds+1PXRY8QElx+elwwRRLy7kUZMWCdM5UnwtzvrJOOm7UzDlB2yNHhR2skdVYf168PCAiROhXTu4cAHGjMm85Hsh6gJ9f+tL8RnFCTgQgLODM1P+PkmizoCNWggFsd1+UorP1qiqkP1q1gyqVBED0aoqOtTmzJHJVyKRvDVYnBSlKTxuY5Nj+Xrdp8y9sgZVEeXiXnEeqCkpKFotH/brwP7AUxh1HtgaPQCR5DuXqkDjxrBzp+gH+vlnUTV9GaZQ2/rr8l+0+KUFtta29K3al5FeIymRpwTFv/zjhd+f1lJ8lrJzp0i4J04I+78ffxTt4ZYoBCKRSCQZxOISMGTM47yPm24AABIpSURBVBZg2aJBzIla8/TvVlYabLv0QNFqn14fniX5AnaOFLxUnRGTHMmVK3VjRent1lZVlX+u/kOyPpmWZVvyXrH3+M77Oz6t8SkFchR4+n2ZUYrPFGJjRVdzvnzi7zqd0OHs2hVsbc0bm0QikZgRiytBp4fzkSfp61+flSu/AODj9t+zWNMeext7NIoGrUaLt7v3c69pU7UIe79oyLdlPyB83nv8vtqRTz6BixeFXeDrZnpf1a39IvRGPWtPr6X6guo0W9kM32BfAOxt7Bn73tjnki+YthSfKdy4IeS/XF3FATlAgwZw6pQ455XJVyKRvOVY5A44tSzYHcD0oMmc4x52Oih97RF0B0cXNz7+Zj0ej716vd298XT1fO61R48KQ50DB8DTE/78E6pXT/3aaenW3nx+MyO3jeRK9BXK5ivL4laL6Vbp1XKRpijFZwonTwpHotWrxdlup05CnxlEqVmWmyUSiQR4gxNwq4nl2WI8BypYqwobKk+kRcevn/seT1fP/yTeBw9EQ9VPP0GBArBsmcgfaZ2GeV2JOCYpBgWF3Ha5MapG8jnkw6+JH63LtcZKSd1iGS3FmwxVfZZYp04VohmDBsHQoeDubtbQJBKJJLvyxpSgk3VJLF37FdFREQDksc+DogIKqBorThSyeuXuy2gUGhBlysDChTBkiOhu7tUrfaOoLysRf/KeE6P/GY3bNDcCQgIAaFOuDQc+OUBbj7apTr7ZAp0OVq6EqlXh2DHxuUmThIjGtGky+UokEskrsNgdcMjj8nGNglUJ3f4zMyJ+JdJOR8qNcD4duZpPO/nx68+NSTGkvPCM998cOSLOdQ8dEt7tc+aILueM8P8l4jy57uNc+G8+37EevVFPpwqdaOvRFsDy7ABjY8VTyvTp4qzXw0N8DqBwYfPGJpFIJBZCuu0I04Op7AhDIkJotKIRSfpE1Me73Ca3HPDx6Evjvj+i5Mjx9PtedsYLwjTh669FLilYEPz8oFu3zDmmbLu2LVsvbX3qSlQyr4XKRRoMwgjh+nWhWOXjA++/LxWrJBKJ5AVklh2h2QgKDyLFkMLj3MtA5w+YNfa355KAmMFNJDKmMr87JeLT7ObTXanRKCZhRo+Ghw+Fuc748aZztVNVlR1hO/Dd78ucFnMona80U5tOZd4H83BxdDHNIlnJyZPCFOGHH4Ts16RJolYvRTMkEokk3VhkAvZ290ar0YrysrWWrq3G/Cf5vmwGtxhFGDBAdDfXqyfKzZUqmSYug9HAxvMbmbxvMkdvHcXF0YWr0Vcpna80JfJYmH2eqsKOHaKj+e+/wcFBHIiXLStmeCUSiUSSISwyAXu6erKj546XlpdfNIMb/0jh84FG7h4EZ2dYsUJID5uq3GwwGqj6U1VO3T1F6bylWfDhAnpU6YGddSaKQ2cWYWFCY/P4cVGb/+EH4cMrrQAlEonEZFhkAoYXjxA94d+ztqoK8WeLEL2rHMYEWwYOELoQTk4ZjyE2OZbfLvxG98rd0Vhp6FmlJ+5O7rQt1xaNleb1F8hOxMUJlZHq1aFIEaFctWiROBTPTIcJiUQieUux2AT8Kp7M4Kbcc+TBPxVJjsiHtlA05fqcYvasmhm+/t34u8w4MIM5h+fwMPkhVQpWoVLBSozyGmWC6LOYyEiYORPmzwdHRwgPB60Wtm83d2QSiUTyRvNGJuDB9csx0CeJBwfcsdLqydvsJPlrRPJd+4wd9t5PuM+4oHEsDl1Msj6Zdh7tGF1nNJUKmugQOSu5dEk0U61cKTqb27UTHc2v09iUSCQSiUl44+62W7bAl4ML8+AaFKgeia3XGVwLa/BpVindqlGPUh7hqHXEztqOwHOBdK3YlS/qfEFZ52yiu5xaVFWIZ2i1cOWK6Gzu3x+GD4eSFjoWJZFIJBbKG5OAr18XyoebNglv9z17oF69wkD6hSH2Xd/H5H2TuRJ9hdOfnyaHNgdXhlzB3iabOQ69Dr0eNmwQHc0NGoCvr/DjjYh45lIkkUgkkizF4tUTdDqRVzw8xLTM5MlCFbFevfRdz6ga+f3i79RdUpd6S+tx8OZBulXqhs6oA7Cs5BsfD7NmiZndjz4SalVPZq4URSZfiUQiMSMWvQPev19Mx5w+DS1bil6ijMoP/37xd1qvaY1bbjdmNZ/Fx1U/xsHGwSTxZjlDh8LixeDlBQEB0KqVVKySSCSSbIJFSlEC/PijcC1ydRWbvNat03edBF0CS0OXYqOxoX/1/uiNegLPBdK2XFtsNDYmiTXLuHhRuBENGiR2uufPC3snLy9zRyaRSCRvJa+SorTY7dATGeKzZ9OXfKMTo/lhzw+4T3dn0NZB/HnpTwCsrazpVKGTZSXf/fuhTRsoVw6WL3/mTFSunEy+EolEkk2x2BK0l1f6c8uy48sYvHUwj1Ie0bxUc76q+xV13eqaNsCsQFVFM9U//wiVqm++EbvfAgXMHZlEIpFIXoPFJuC0cun+JRy1jhTKWYgSeUrQskxLRtcZTRWXKuYOLW0kJsLmzdC5s2ikatRInO326QOPXaAkEolEkv2x2DPg1HLs1jEm75vM+rPrGfTuIGY2n5ml65uMqCiYOxdmz4Z79yA4GDxfLMUpkUgkkuzBG2dHmBp2h+/mx30/su3KNnLZ5mJ0ndEMrT3U3GGlnehoGDtW+CcmJkKLFuLwu3Ztc0cmkUgkkgzwRiVgVVVRHtsbLTuxjOO3jzOp0SQ+r/E5ue1ymzm6NPLggTjXtbcX8l6dO8OoUUJlRCKRSCQWzxtRgtYZdKw6tYop+6ewou0KahSuQVRCFDlscliWcIbRCFu3gp8fXLsm9JqtrSE5GWxtzR2dRCKRSNLIGzmGBGKGd+bBmZScWZLem3tjo7EhUSesCJ0dnC0n+SYnw9KlYnb3ww//1969x1hRnnEc/z6sUAOliIVSiligGk0bETZ441aj9QI10molGLCkkBAFFRNNAU28Ro00rdoGRWtJaYst1hZr8FIJVhsToSByM1hZFMJWChQVi41R2ad/vLN6cpjZPewezjtz/H2SzZmdec/yPLxn5tl5Z3ZeaGqCWbPCIyRBxVdEpA4Vdgj6YMtBhj44lG3vbmP08aNZeNFCxp0w7tMh6EJZsQKmTQsFePHi8NjIbt1iRyUiIkdQYQtwQ5cGbvn2LQzuPbh4f8Pb3Az33Qd9+sDcueHGqpUrw0QJRfwFQkREDlthCzDAFadeETuEw7NpU5g54tFHw0M0ZswI67t0gXPOiRubiIjUVKGvARfKXXfB0KHw+OMwc2a4zvvAA7GjEhGRSAp9BpxrrXPwNjbCiSfCeeeFu5yvukrTAIqISOfPgM2swcxeNbPl1Qio8Mrn4F20KKw/7bTwrGYVXxERoTpD0LOBLVX4OcV3991w/PFw7bXQvz8sWwZ33hk7KhERyaFOFWAzOw74LvBIdcIpoB07wg1VEO5uHjMGXnrpsykCu+gyu4iIHKqz1eE+4MdAS1YDM5thZmvNbO3evXs7+c/lyOrVcOmlMHhwKLYQhp6feAJGjYobm4iI5F6HC7CZXQTscfdX2mrn7g+7+wh3H9G3b9+O/nP50NICy5fD2LFhMoTnn4d588JNVqCzXRERqVhn7oIeBVxsZuOBo4Evmdnv3H1KdULLoQ8/DPPudu8O994L06dDz56xoxIRkQLqcAF293nAPAAzOxu4oe6K7/798NBD8NRT4Wy3e/fwevLJ0LVr7OhERKTANGaaprk5zLk7cCDMmROK7b59Ydspp6j4iohIp1XlQRzu/gLwQjV+VnRr1sDIkeHO5okTwxy8jY2xoxIRkTqjJ2G5w4svwu7dYdL7xka48cZwrXfQoNjRiYhInfr8DkEfPBiey3zGGWEWojvuCMW4oQFuu03FV0REjqjPZwF+9lk46SS47DJ47z1YuDAMPWsqQBERqZHPzxD0vn1hgoR+/aBHjzAX7/z5MGFCOOsVERGpofo/A37rLbjmmvCM5ttvD+tGj4aXX4ZLLlHxFRGRKOr3DPjVV8MZ7mOPhSI7eXKYhxc01CwiItHVVwF2/6y4LlgQHqBx/fUwezYMGBA3NhERkRL1MQT98cewZAkMHx4mSYAwDeDOneEsWMVXRERyptgF+MABuP9+OOEEmDIFPvoIPvggbOvXD3r1ihufiIhIhuIOQbe0wLBhsG1bmIN3wQIYP14zEomISCEUtwB36RIenjFoEJx1VuxoREREDktxCzDA5ZfHjkBERKRDNF4rIiISgQqwiIhIBCrAIiIiEagAi4iIRKACLCIiEoEKsIiISAQqwCIiIhGoAIuIiESgAiwiIhKBCrCIiEgEKsAiIiIRqACLiIhEoAIsIiISgbl77f4xs73Ajir+yD7Af6r482JSLvlTL3mAcsmresmlXvKA6ufydXfvm7ahpgW42sxsrbuPiB1HNSiX/KmXPEC55FW95FIveUBtc9EQtIiISAQqwCIiIhEUvQA/HDuAKlIu+VMveYByyat6yaVe8oAa5lLoa8AiIiJFVfQzYBERkUJSARYREYmgEAXYzC40s3+aWZOZzU3Z/gUzW5psX21mg2ofZfvMbKCZ/c3MtpjZa2Y2O6XN2Wa238zWJ183x4i1Ema23cw2JXGuTdluZvbzpF82mlljjDjbYmYnlfxfrzez983surI2ue0TM1tkZnvMbHPJumPNbIWZbU1ee2e8d2rSZquZTa1d1OkycvmJmb2efH6WmdkxGe9t87NYaxm53Gpm/yr5HI3PeG+bx7tayshjaUkO281sfcZ789YnqcffqPuLu+f6C2gAtgFDgG7ABuCbZW1mAguT5UnA0thxZ+TSH2hMlnsCb6TkcjawPHasFeazHejTxvbxwDOAAWcCq2PH3E4+DcC/CX84X4g+AcYCjcDmknXzgbnJ8lzgnpT3HQu8mbz2TpZ75zCX84GjkuV70nJJtrX5WcxJLrcCN7TzvnaPd7HzKNv+U+DmgvRJ6vE35v5ShDPg04Emd3/T3T8C/gBMKGszAVicLD8OnGtmVsMYK+Luu9x9XbL8X2ALMCBuVEfUBOA3HqwCjjGz/rGDasO5wDZ3r+bT2o4od/878E7Z6tL9YTHwvZS3XgCscPd33P1dYAVw4RELtAJpubj7c+7+SfLtKuC4mgfWARn9UolKjnc101YeyTF2IvD7mgbVQW0cf6PtL0UowAOAnSXfN3No0fq0TbKz7ge+XJPoOigZJh8OrE7ZfJaZbTCzZ8zsWzUN7PA48JyZvWJmM1K2V9J3eTKJ7INJUfoEoJ+774Jw0AG+ktKmaH0DMI0wopKmvc9iXlydDKcvyhjqLFK/jAF2u/vWjO257ZOy42+0/aUIBTjtTLb8b6cqaZMbZvZF4E/Ade7+ftnmdYQh0FOBXwBP1Dq+wzDK3RuBccAsMxtbtr0w/WJm3YCLgT+mbC5Sn1SqMH0DYGY3AZ8ASzKatPdZzIMHgW8Aw4BdhOHbckXql8tp++w3l33SzvE3820p6zrdL0UowM3AwJLvjwPezmpjZkcBvejY8M8RZ2ZdCZ2/xN3/XL7d3d939wPJ8tNAVzPrU+MwK+Lubyeve4BlhOGzUpX0XV6MA9a5++7yDUXqk8Tu1qH+5HVPSpvC9E1yw8tFwGRPLsiVq+CzGJ2773b3g+7eAvyS9BgL0S/JcfYSYGlWmzz2ScbxN9r+UoQCvAY40cwGJ2cpk4Any9o8CbTelfYD4PmsHTWm5JrJr4At7v6zjDZfbb1+bWanE/poX+2irIyZ9TCznq3LhJtlNpc1exL4oQVnAvtbh3pyKPO3+aL0SYnS/WEq8JeUNn8Fzjez3slQ6PnJulwxswuBOcDF7v6/jDaVfBajK7v/4fukx1jJ8S4PvgO87u7NaRvz2CdtHH/j7S+x70yr5ItwN+0bhLsDb0rW3U7YKQGOJgwdNgH/AIbEjjkjj9GEYYuNwPrkazxwJXBl0uZq4DXC3Y+rgJGx487IZUgS44Yk3tZ+Kc3FgAVJv20CRsSOOyOX7oSC2qtkXSH6hPBLwy7gY8Jv6dMJ9z+sBLYmr8cmbUcAj5S8d1qyzzQBP8ppLk2Ea2+t+0vrXzt8DXi6rc9iDnP5bbIfbCQc9PuX55J8f8jxLk95JOt/3bp/lLTNe59kHX+j7S96FKWIiEgERRiCFhERqTsqwCIiIhGoAIuIiESgAiwiIhKBCrCIiEgEKsAiIiIRqACLiIhE8H+A6bAV6WYWCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res_wls)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "# OLS\n",
    "ax.plot(x, res_ols.fittedvalues, 'r--')\n",
    "ax.plot(x, iv_u_ols, 'r--', label=\"OLS\")\n",
    "ax.plot(x, iv_l_ols, 'r--')\n",
    "# WLS\n",
    "ax.plot(x, res_wls.fittedvalues, 'g--.')\n",
    "ax.plot(x, iv_u, 'g--', label=\"WLS\")\n",
    "ax.plot(x, iv_l, 'g--')\n",
    "ax.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feasible Weighted Least Squares (2-stage FWLS)\n",
    "\n",
    "Like `w`, `w_est` is proportional to the standard deviation, and so must be squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.931\n",
      "Model:                            WLS   Adj. R-squared:                  0.929\n",
      "Method:                 Least Squares   F-statistic:                     646.7\n",
      "Date:                Tue, 17 Dec 2019   Prob (F-statistic):           1.66e-29\n",
      "Time:                        23:40:04   Log-Likelihood:                -50.716\n",
      "No. Observations:                  50   AIC:                             105.4\n",
      "Df Residuals:                      48   BIC:                             109.3\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2363      0.135     38.720      0.000       4.964       5.508\n",
      "x1             0.4492      0.018     25.431      0.000       0.414       0.485\n",
      "==============================================================================\n",
      "Omnibus:                        0.247   Durbin-Watson:                   2.343\n",
      "Prob(Omnibus):                  0.884   Jarque-Bera (JB):                0.179\n",
      "Skew:                          -0.136   Prob(JB):                        0.915\n",
      "Kurtosis:                       2.893   Cond. No.                         14.3\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "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 ** 2))).fit()\n",
    "print(res_fwls.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}