{
 "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:                Sun, 16 Aug 2020   Prob (F-statistic):           5.44e-29\n",
      "Time:                        18:00:44   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+YZr1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAACRnklEQVR4nOzdd1hURxfA4d9l6aKiYkUQe+8VjYq99941UWOL3Wg+001sIPbejcbYUKPG2LGBBcWOXRTBiiAibcv9/hjRmFgQFnbReZ/HJ7rs3jug2bMzc+YcRVVVJEmSJElKWxamHoAkSZIkfYpkAJYkSZIkE5ABWJIkSZJMQAZgSZIkSTIBGYAlSZIkyQRkAJYkSZIkE3hvAFYUZZmiKA8VRbnwj8fKKYpyTFGUM4qiBCiKUiV1hylJkiRJH5ekzIBXAI3/9dhU4CdVVcsB37/4syRJkiRJSWT5vieoqnpIURS3fz8MZHrx+8xAWFJu5uTkpLq5/ftSkiRJkvRxOnXq1GNVVbO/6WvvDcBvMRzYpSiKF2IWXf1tT1QUpT/QH8DV1ZWAgIBk3lKSJEmS0hdFUW6/7WvJTcIaCIxQVdUFGAEsfdsTVVVdpKpqJVVVK2XP/sYPAZIkSZL0yUluAO4F+Lz4/QZAJmFJkiRJ0gdIbgAOA2q/+H1d4JpxhiNJkiRJn4b37gErirIW8ACcFEW5C/wA9ANmKopiCcTxYo83ObRaLXfv3iUuLi65l0g3bG1tyZs3L1ZWVqYeiiRJkmRiScmC7vKWL1U0xgDu3r1LxowZcXNzQ1EUY1zSLKmqSnh4OHfv3iV//vymHo4kSZJkYiavhBUXF0e2bNk+6uALoCgK2bJl+yRm+pIkSdL7mTwAAx998E30qXyfkiRJ0vuZRQCWJEmSpE9NugvAWwJDqTF5P/nH7aDG5P1sCQxN8TU1Gg3lypWjZMmSlC1bFm9vbwwGwztfExwczO+//57ie0uSJEmfpuRWwjKJLYGhfONznlitHoDQyFi+8TkPQOvyzsm+rp2dHWfOnAHg4cOHdO3aladPn/LTTz+99TWJAbhr167Jvq8kSZKUNrYEhuK56wphkbHkcbRjTKOiKYobxpCuZsCeu668DL6JYrV6PHddMdo9cuTIwaJFi5gzZw6qqhIcHEzNmjWpUKECFSpUwM/PD4Bx48Zx+PBhypUrx/Tp09/6PEmSJMm0EidvoZGxqLyavBljBTUl0tUMOCwy9oMeT64CBQpgMBh4+PAhOXLkYM+ePdja2nLt2jW6dOlCQEAAkydPxsvLi+3btwMQExPzxudJkiRJpvWuyZspZ8HpKgDncbQj9A3BNo+jndHvpaoqIAqFDBkyhDNnzqDRaLh69eobn5/U50mSJElpK60mbx8qXS1Bj2lUFDsrzWuP2VlpGNOoqFHvc/PmTTQaDTly5GD69OnkzJmTs2fPEhAQQEJCwhtfk9TnSZIkSWnrbZO01Ji8fYh0FYBbl3dmUtvSODvaoQDOjnZMalvaqEsIjx49YsCAAQwZMgRFUXj69Cm5c+fGwsKC3377Db1eLGNkzJiRZ8+evXzd254nSZIkmVZaTd4+VLpaggYRhI29Zh8bG0u5cuXQarVYWlrSo0cPRo4cCcCgQYNo164dGzZsoE6dOmTIkAGAMmXKYGlpSdmyZendu/dbnydJkiSZVmLMMLcsaCVxrzMtVKpUSf13YlJQUBDFixdPszGY2qf2/UqSJH3KFEU5papqpTd9LV0tQUuSJEnSx0IGYEmSJEkyARmAJUmSJMkEZACWJEmSJBOQAViSJEmSgGvh19L0fjIAS5IkSZ+0M/fPUHtFbYrNLcb1J9fT7L7p7hywsYWHh1OvXj0A7t+/j0ajIXv27ACcOHECa2trUw5PkiRJSgU6g46I2AiyZ8iOnaUdtyNv493Qm9wOudNsDJ98AM6WLdvLVoQ//vgjDg4OjB49+uXXdTodlpaf/I9JkiTpoxCrjWX5meV4+nlSPld5fDr5UNSpKDeG3kBjoXn/BYzIrCLL8OHwIhYaTblyMGPGh72md+/eZM2alcDAQCpUqEDGjBlfC8ylSpVi+/btuLm5sXr1ambNmkVCQgJVq1Zl3rx5aDRp+5coSZIkvVtEbATzA+Yz49gMHsU8olreavQu1/vl19M6+ILcA36rq1evsnfvXqZNm/bW5wQFBbFu3TqOHj36sgvSmjVr0nCUkiRJUlLMODaD8fvHUzFPRQ72Pojf5360LNrSpGMyqxnwh85UU1OHDh3eO5Pdt28fp06donLlyoCoKZ0jR460GJ4kSZL0DlceX8HTz5PWxVrTvEhzvqr6FW2Lt6VsrrKmHtpLZhWAzck/mylYWlpiMBhe/jkuLg4QPYN79erFpEmT0nx8kiRJ0n+dCD3BlKNT2By0GRtLG0rnKA2Ak70TTvZOJh7d6+QSdBK4ublx+vRpAE6fPs2tW7cAqFevHhs3buThw4cAPHnyhNu3b5tsnJIkSZ+yPlv7UHVJVfbf2s//av6P28NvM6zaMFMP663kDDgJ2rVrx6pVqyhXrhyVK1emSJEiAJQoUYJffvmFhg0bYjAYsLKyYu7cueTLl8/EI5YkSfr46Q16Nl/eTPMizbG1tKVe/nqUyl6K/hX7k9Emo6mH916yHWEa+9S+X0mSJGOL08Wx8sxKPP08uRFxg1WtV9GjbA9TD+uN3tWOUM6AJUmSpHRBZ9Dh5efFjGMzePD8AVWcq+DZwJNWxVqZemjJIgOwJEmSZNZitDHYW9mjUTRsCtpEuVzlGFtjLB5uHiiKYurhJZsMwJIkSZJZuhZ+DS8/LzYGbeTqkKtks8+Gby9fMlhneP+L0wEZgCVJkiSzEhAWwJSjU9h0aRPWGmv6lOuDzqAD+GiCL8gALEmSJJmRWxG3qLy4MplsMjHus3EMrTqUXA65TD2sVCEDsCRJkmQyeoMenyAfLjy8wE91fiJ/lvysb7+ehgUbktk2s6mHl6reW4hDUZRliqI8VBTlwr8e/0pRlCuKolxUFGVq6g0x9d29e5dWrVpRuHBhChYsyLBhw0hISMDX15fmzZv/5/nbt2+nfPnylC1blhIlSrBw4UITjFqSJCn9itPFsejUIorNLUbHjR1Zf2k9cTpRZbBDyQ4fffCFpFXCWgE0/ucDiqLUAVoBZVRVLQl4GX9oaUNVVdq2bUvr1q25du0aV69eJTo6mvHjx7/x+Vqtlv79+7Nt2zbOnj1LYGAgHh4eaTtoSZKkdMw32Jf8M/Pz5fYvcbR1ZEOHDVwYeAFbS1tTDy1NvXcJWlXVQ4qiuP3r4YHAZFVV418856HRRvSmYNaxIwwaBDEx0LTpf7/eu7f49fgxtG//+td8fd95u/3792Nra0ufPn0A0Gg0TJ8+nfz581OnTp3/PP/Zs2fodDqyZcsGgI2NDUWLFn3/9yVJkvQJux99n4jYCIpnL06RbEUon6s8o9xHUTd/3XR9lCglklsLughQU1GU44qiHFQUpfLbnqgoSn9FUQIURQl49OhRMm+Xei5evEjFihVfeyxTpky4urpy/fr1/zw/a9astGzZknz58tGlSxfWrFnzWqMGSZIk6ZXrT64zYPsA3Ga4MfivwQDkyZiHv7r9Rb0C9cwj+KZhRch/Sm4SliWQBagGVAbWK4pSQH1DXUtVVRcBi0CUonzvld81Y7W3f/fXnZzeO+N9w/je+A/gbY8DLFmyhPPnz7N37168vLzYs2cPK1as+KD7SpIkfczO3j/LxCMT2XhpI1YWVvQq24sxNcaYelivi4+H1atFL9w//4T8+dP09smdAd8FfFThBGAAzKvPUxKVLFmSf9enjoqKIiQkhIIFC771daVLl2bEiBHs2bOHTZs2pfYwJUmSzJ6qqhhUsSK479Y+/r7+N19X/5rg4cEsbLGQQlkLmXiEL0RGwpQpIuD27QuWlmILM40lNwBvAeoCKIpSBLAG0n70RlCvXj1iYmJYtWoVAHq9nlGjRtG7d2/s7e3/8/zo6Gh8/zHLPnPmjOx+JEnSJ01v0LPx0kaqLKnCqrPivXRApQHcGX6HSfUnmdc53pgYKFQIxo2DUqVg9244fRoqv3UnNdUk5RjSWsAfKKooyl1FUb4AlgEFXhxN+gPo9abl5/RAURQ2b97Mhg0bKFy4MEWKFMHW1paJEycCsG/fPvLmzfvyV2BgIFOnTqVo0aKUK1eOH374QS4/S5L0SfrnUaIOGzoQGRdJJptMANhb2ZvPUaKLF8WMF8RW5sSJIuju3g0NGoCJ9qFlO8I09ql9v5IkfbwarW7E7hu7qZSnEmNrjKVNsTZoLDSmHpagqnDoEHh6wo4dIvBevQrOzmk6jHe1I0zuErQkSZL0ibkffZ/x+8YTGRcJwNgaY9nbYy8n+p6gfYn25hN8r16FatXEsdYTJ+Dnn+HOnTQPvu8jS1FKkiRJ73T9yXW8/LxYcWYFWoOWinkq0rZ4W+rmr2vqob0SGwshIVCkCOTODQYDzJ8PvXqBnZ2pR/dGMgBLkiRJb5SgT6DH5h4vjxL1Lteb0dVHm082M0B4OMybB7NnQ44ccO4cZMwIJ0+aemTvJQOwJEmS9JKqqlx6dImSOUpirbFGVVW+rv41w6oNM69s5uBgmDYNli0Tmc3NmsGYMSZLqEoOGYAlSZIk9AY9m4I2MeXoFM49OMeNoTdwzezK+g7rTT201xkMYGEBR47AwoXQrRuMHg0lS5p6ZB9MJmFJkiR9wv55lKjTxk5EJ0Qzv9l8cmbIaeqhvaKqsGsX1K8PXi96/3TqBLduwfLl6TL4ggzAjBgxghkzZrz8c6NGjejbt+/LP48aNQpvb29KlSr1n9ceO3aMqlWrUq5cOYoXL86PP/6YBiOWJEkyngfRDxj812AcbR3Z2GEjlwZdom+FvthY2ph6aKDVwm+/Qbly0LgxBAVBlizia1ZWZpfV/KE++QBcvXp1/Pz8ADAYDDx+/JiLFy++/Lqfnx81atR442t79erFokWLOHPmDBcuXKBjx45pMmZJkqTkuvfsHmP3jKXDhg4A5HPMx7kB5zjR9wTtSrQzn6NEILrc9ewJOp2Y6d66Bf36mXpURmN2e8AeKzz+81jHkh0ZVHkQMdoYmq75bzvC3uV607tcbx7HPKb9+tfbEfr29n3n/WrUqMGIESMA0RmpVKlS3Lt3j4iICOzt7QkKCiJL4ieuf3n48CG5c+cGRBvDEiVKJOE7lCRJSntXw6/iedSTVedWoTPo6FiyIwn6BKw11hTPbibFge7fh1mzYOBAcHGBYcOga1do0kTs+35kzC4Ap7U8efJgaWnJnTt38PPzw93dndDQUPz9/cmcOTNlypTB2tr6ja8dMWIERYsWxcPDg8aNG9OrVy9sbT+thtKSJJm/TZc20WFDB6w11nxR/gtGuY+iYNa3N5tJc5cvi4zmVavEbLdoUXF+t0oVU48sVclSlEC3bt1o0aIFO3fuZOTIkYSGhuLn50fmzJkJDw9nwIABNG/enAsXLvzntTdu3GD37t388ccfKIryWqOGNzGH71eSpI+bqqrsubkHjaKhXoF6RMZFMs1vGkOqDCGngxklVxkM0LEjbNoEtrZiyXnUKNEs4SMhS1G+R+I+8Pnz5ylVqhTVqlXD39//nfu/iQoWLMjAgQPZt28fZ8+eJTw8PI1GLUmS9DqdQce6C+uouKgijVY3wtPPEwBHW0cm1J1gHsFXr4fDh8XvLSxE1arvv4fbt0Xlqo8o+L6PDMCIfeDt27eTNWtWNBoNWbNmJTIyEn9/f9zd3d/6uh07dpC4gnDt2jU0Gg2Ojo5pNGpJkqRXNl7aSNE5Rem8qTMx2hiWtlzK1s5bTT2sV2Jjxbnd4sWhVi1RsQpEBauffhJVrD4xn/weMEDp0qV5/PgxXbt2fe2x6OhonJyciI6O5sqVK+TNm/fl16dPn86mTZsYMWIE9vb2WFpasmbNGjQaM8oglCTpoxYRG4G1xpoM1hmITogmu312vBp40apYKywUM5lfRUWJxKrZs+HhQ6hUCdatA5m0KveA09qn9v1KkmR8d6PuMt1/OotOL+Inj58Y6T4Sg2pAQUExl1KMWq04qxsRAfnywWefwddfQ+3a6apcZEq9aw9YzoAlSZLSiaBHQXj6ebL63GoMqoFOpTrRoEADAPOZ8Z4+LXrwXr8uWgFmyQI3bkD27KYemdmRAViSJCmdGPzXYI7dPcaXFb9kpPtI8mfJb+ohCaoKu3eLwLtvn+hG9OWXEB8vsptl8H0jswjAqqqaz7JJKkrL5X5JktI3VVXZeX0n049NZ2XrleTJmIf5zeaT1S4r2TOYWUDbuFEcJ8qTB6ZOhf79IXNmU4/K7Jk8ANva2hIeHk62bNk+6iCsqirh4eGyUIckSe+k1WtZd3EdU49O5fzD8+TNlJcbT26QJ2MeijoVNfXwhKgoWLwYnJxEwYyWLUURjU6d4C2Fi6T/MnkAzps3L3fv3uXRo0emHkqqs7W1fS2TWpIk6Z9itbGUnl+aGxE3KJG9BCtaraBL6S5Ya8wkqIWFwcyZ4jjR06fQo4cIwDY24vfSBzF5ALaysiJ/fjPZx5AkSUpj4THh/H39b7qV6YadlR29y/WmbM6yNCvSzHwSq0CUivzmG1FIo317GDNGHCmSks3kAViSJOlTdDvyNt7+3iwJXEKsNpbPXD8jn2M+vq31ramHJqgqHDwIxYpBrlxQqpTY2x05EgoUMPXoPgpm9PFKkiTp4xcaFUrPzT0pNLsQ8wLm0b5Ee84NPEc+x3ymHpqg08H69aIRQp06ojwkQKNGMGeODL5GJGfAkiRJqUxVVZ7GP8XR1hFbS1t239jNkMpDGOE+AtfMrqYe3ivz54ujRLduQeHCsGCB6McrpQoZgCVJklKJQTWw/ep2Jh+ZjNag5UTfE2Szz8adEXfMJ7EqKgoyZRK/P3AAcuYU+70tW4IsrZuq5BK0JEmSkSXoE1geuJxS80rR6o9W3Iu+R6+yvTCoBgDzCL7XrsGAAWJ/99Il8diKFeDnB23ayOCbBuQMWJIkychWn1vNF39+QdmcZVnTdg0dS3bE0sJM3m79/cUy85Yt4sxuz57g4CC+Zm9v0qF9aszkX4QkSVL69SD6AbOOz6JwtsL0LtebrqW74pzRmYYFG5pXgaHISKhXT5SHHD8ehgwRS86SScgALEmSlEw3I27i5efF8jPLidfFM7TqUABsLW1pVKiRiUcHxMWJClUHD8Lq1eDoCH/9BZUrQ4YMaTaMLYGheO66QlhkLHkc7RjTqCityzun2f3NlQzAkiRJyfDzwZ/56eBPWFpY0rNMT0ZXH20+pSLDw2HePHFs6OFDqFhRtAXMmhU8PNJ0KFsCQ/nG5zyxWj0AoZGxfONzHuCTD8IyCUuSJCkJVFVl3819PHouyuZWzlOZ0e6juTXsFotbLjaf4HvkCLi6wvffi0pVBw7AyZMi+JqA564rL4NvolitHs9dV0wynnfxD/Fn0uFJ+If4p8n95AxYkiTpHfQGPT5BPkw5OoVT904xoc4Evq31LU0KN6FJ4SamHp5w8qSY4TZsKGa7vXvDoEFQsqSpR0ZYZOwHPW4qfnf8qLuqLjqDDmuNNft67sPdxT1V7ykDsCRJ0lssPrWYqX5Tuf7kOoWyFmJh84X0LGsmhSkMBti5U2Q0HzwIFSqIAGxnB3Pnmnp0L+VxtCP0DcE2j6Od0e5hjD3mUbtHEa+PB8QxMt9g31QPwO9dglYUZZmiKA8VRbnwhq+NVhRFVRTFKXWGJ0mSlLZita+CxV/X/8LR1pENHTZwefBl+lfsj62lGbQU3bFD1GZu3hxu3hSFMw4cMPWo3mhMo6LYWb1+ptjOSsOYRsZZsk/cYw6NjEXl1R7zlsDQd77uWfwzpvtP52bETQD6VeyHtYU1GkWDtcYaDzcPo4zvXZIyA14BzAFW/fNBRVFcgAbAHeMPS5IkKW2FRoUy/dh0lpxewol+JyiSrQirWq/CwdrBPI4SRUSAhYVodB8fL87wrl4NHTuClZWpR/dWiTPR1MqCftce87/v4R/iz/ar27kbdZc/r/5JZFwkKioj3UfyefnPKe5UHN9gXzzcPFJ99gtJCMCqqh5SFMXtDV+aDnwNbDX2oCRJktJK0KMgPP08WX1uNXpVT6eSndAoYsaW0SajiUcHBAfDjBmwZAmMGgU//QStW4tqVebwwSAJWpd3TrWM56TuMfvd8aP2ytroDDoAPPJ5MLn+ZKrmrfryOe4u7mkSeBMlaw9YUZSWQKiqqmff98lQUZT+QH8AV1czKjouSdInLyo+ikqLK2FQDfSv2J9R7qPIn8VM+pOfPi32dzdsEIG2SxfRhxfETFgC3r/HfOHhBUrlKMXB2wfRG8RMWaNoaFiw4WvB1xQ++G9RURR7YDzwfVKer6rqIlVVK6mqWil79uwfejtJkiSjUVWVv679xbCdwwDIZJOJP9r9wZ3hd5jTdI7pg6+qvvr9zz+Lvd4RI8Q+76pVULq06cZmpt60x2xrpVCn3F1qr6hN6fmlOX73OB5uHtha2v53jzcuTtTATkhI87EnZwZcEMgPJM5+8wKnFUWpoqrqfWMOTpIkyRi0ei1/XPiDqX5TufDwAi6ZXBhfazw5MuSgRdEWph6e2NP9/XeYPh02bRKtAGfNEvu9mTOn6q3Te5WqxLF+v9OHkJgTZLBVibEP4NfjV3DJ5ML0RtMpkb0EGW0ysq/nvld7vDYFxXL+3Lnw6BFkyQKtWqXp2D84AKuqeh7IkfhnRVGCgUqqqj424rgkSZKM4uz9s7T8oyV3nt6hZPaSrGq9is6lOmOlMYPEpchI0XN31iy4d0/McMPDRQBOgy279FKl6n0fEnJku811dSwJ1glE6vUUsCzAb21+o1PJTq/9Pbu7uOOetYxYVVi1Snzwad5c7K3Xrp3m39d7A7CiKGsBD8BJUZS7wA+qqi5N7YFJkiQl16Pnj7j99DaV8lSiUNZClMlZhrlN59K0cFMsFDPZP42LE4H28WOoX18sgzZokKaJVR+SQWwq7/qQULWQBbOOz+L387+ToE9Ar+rRKBq+KP8F3ct0f3URVYXbt8HNTXR8CgyEXr1EIC5WzATflZCULOgu7/m6m9FGI0mSlAK3Im4xzX8aywKX4ZrZlaDBQWSwzsC2LttMPTTh9GnYtg1++EF0JJo6FcqXh3LlTDKc9FCl6k0fEqJ0t/ly+2wilb3oDDo88nnwKOYRCfoErDXW1HGrI56o1YoktmnT4MYNCAmBjBnh2DGz6HcsK2FJkpTuXXp0iV8O/cL6i+uxUCzoUaYHo6uPNo/zu6oKf/8NXl6wf78IAF98AXnzQp8+Jh1aWlSpSqnEDwPxFkHEWZxHUe2JsFqIordiQOUvGOk+kkJZC+Ef4v9qfzdTCfHznjkT7t4Vs1xPT3F2Gswi+IIMwJIkpVOqqqIz6LDSWHHh4QW2Xd3G8GrDGVFtBM6ZzGP5lKAgUSjjwgVwdhYz3v79Uz2xKqnGNCr62vIuGLdKlTHkzmxDUPTvRFr9BhhQsMRB35hi9l8wr1m7l89zd3HHPW81sYR/6hSMGQN16og99iZNzPLolgzAkiSlK3qDni2XtzDl6BRaFGnBd7W/o13xdjQo0IAsdllMPTyRWBUcLJaVXV0hWzZYuRI6d341AzMTqV2lKiUS9An8fv537tpMIjLhKqiAAqqqw07JyfjG1V49+cQJscycMaMoWFKxIly+DEXN54PEm8gALElSuhCni2PV2VV4+Xlx7ck1CmUtRIEsBQDQWGhMH3xv335VscrZWcx+M2QAX1/Tjus9UrNKVXL5BPkwdOdQQp+FUjZnWarnHsmGq3MxqFosFCtG1m5D6zK5YMsWEXiPHBGrCl999eoiZh58QQZgSZLSib5/9mXN+TVUzF2R9e3X07Z4WzQWZrCXd+ECTJwI69e/qlg1alS6KRNpLrZd2cax0GM0L9wcR1tHijkVY1mrZTQo0ABFURgW0v71Os0//ijO8ebLJz74fP65mAGnI4r6z8orqaxSpUpqQEBAmt1PkqT0K+xZGDOOzWBgpYHkz5KfM/fP8CT2CXXc6pg+ucpgEBm2Njawbh306wdffglDh4KLi2nHls5cfnyZMXvGsP3qdhQUbC1t39yL9/59mDMHGjWCmjXFisPx49C2LVia71xSUZRTqqpWetPXzHfUkiR9sqIToim7oCxPYp9QInsJ8mfJT7lc5Uw9LFG4Yc0asezZtSuMHw/t2kHjxmaTWJVe+If4M+XoFLZe2YqlhSUKCirqf3vxXrggft6//y4+9GTIIAJwvnziVzpmfmlhkiR98gLvBfI45jHr2q+jd7neph6OaAU4aZIo5PDFF6L9X4kS4muWljL4JoF/iD8TD0/EP8QfgPkB8zl85zA/1P6BrZ23vrlOc69eojrY+vVileHqVfjmG9N9E0YmZ8CSJJmdgDCxVVXTtaaJR/LCF1/A5s1i+fO336BePbnH+wEOBh+kwW8N0Bq02GhsONDrAJ4NPJnfbD4ZrDMAiDrNN/biccOAe+7K4oVVq4pkqgEDIGtWE34HqUMGYEmSzE7AvQBcMrmQ0yGnaQZw8qRY9pw0CfLnF8k+P/4IZcqYZjzp1NO4pyw8tZBfDv2C1qAFQGvQvr7EDPDkCe6/HcB9znxREztnRVGjedAgE408bcgALEmS2RlceTBtirVJ25saDKL9n5cXHDoEmTKJfd78+WUbwGTQGXSUml+Ku1F3qZSnEucfnEdn0L2+xBwdDePGwfLlEBMDDRu+qon9CZABWJIks1PdpXra3lCnE8Ubzp0TWcze3mLZOVOmtB1HKkiLdoOJZSDzOeYj6FEQP9f5GUsLS6bUn0LRbEWpmKfiq1KR+WrjrrzIFLe3Fx92OnaEkSM/uQ868hiSJElm5VbELa6EX3nZQD3VhIfD9u0i0Qfg11+hQAFo314kWX0E/t1JCESpyUltSxstCPuH+FNnZR3i9fEAWGusOTvgLMWc/tVlSKcT++heXnDtmmiMkCGDeNyMjxGl1LuOIcksaEmSzMrGSxtpsqYJzxOep84Nrl+HwYPFTLd3bxEMQBwp6tLlowm+8O52g8Zw5+kdOm/q/DL4KiiMqT7m9eD77JloilC4sJjphofDhAmvGiJ8xMH3fWQAliTJrJwMO0l+x/xks89m3AuHhIgzu0WKiHKRXbqIM6aFCxv3PmYkNdoNxuviufDwAgC5HHLhaOuIlYUVGkWDraUtzQo3E09MXF29eBGGDxflOX184MoV8QHINhVXN9KJT/ejhyRJZikgLIBKed64Yvfh9HoICxOz3cyZRT/eb76BIUMgd27j3MOMGavdoH+IPzuv7+RxzGM2X96MRtFwc9jNl8vNr7UCfGwL33QXP++5c6FaNTh7VmaQv4EMwJIkmY3wmHBuRd5iQKUBKbtQTIzoQOTtLcpFnj8vEqquXzebXrBpwRjtBv+88ift1rdDZ9ABUMW5ChPrTsTK4tVSvbtzVdzPPYFe4+HAAXBwEB9yEsng+0YyAEuSZDZO3TsFQOU8lZN3gQcPxKxr3jyx11iliugLq6qicMYnFHwhZe0GDaoBC8WCzUGbXwZfjaKhddHW1CtQ7+XztgSG8mT0N3y+/zceZnLi4bDxlPpxNDg6psr39DGRWdCSJJkNnUHHpUeXKJS1EPZW9kl/YWKA/e03kdXcsiWMHg01asiKVR9AVVWO3DnCVL+plMxeksn1J+N3x496v9VDq9dirbEWjRJsC8H8+RzMV44B163J/vAu5cMus6NYTaxsbYyaZZ3eyWYMkiSlC5YWlpTJmcTlSlUVZ0i9vKB2bRFwO3US5QuLFEndgX5kjtw5woKABZy5f4aLjy7iZO+ERz4PAKq7Vmd/z/1ij9ciP+6/rIBVqyAujssNPye2fFvuZMnNnSxiT133Iss6PQXg2FgIDIQTJ8Sfhw9Pm/vKACxJktn4377/0bhQY2rlq/X2J+l0sGmTKBV58iQ4OYkKSgDW1jL4fqDEc7w6gw4FhVHuo/i5zs+vrUC4u7jj/v0iWPE/safesyeMGMHklTffeM2UZFmnNr0eLl0SwTbx1/nz4nEQn99kAJYk6ZNyP/o+k45Mwsne6d0B+PPPxVJzkSKwYIEIBnYfltX7sUhulavIuEgWBCygRZEW+Ab7YlANAFgoFmSzyyaCr1YLf/4JrVuLvfMyZeCHH0R95hw5AMjjeM8oWdap6d49OHbs1a9Tp+D5iyPmmTOLNIFx46ByZfErT560G5sMwJIkmYVTYW9JwLp7F2bPFlm1Li7iDGn79qJYv8WnW8rg31WuQiNj+cbnPMBbg/DdqLvMODaDhacWEp0QjUbR4OHmgY3GhgR9gqjT7FQJPD1h1izxs9+xA5o2hREj/nM9Y2RZG1N8vFhKTgy2/v5w5474mpUVlCsnPr9VqSJ+FSr06p/QlsBQOqxK3ZKd/yYDsCRJZiEgLAAFhfK5y4sHzp4Vy8xr14pGCaVKQY8eYo1QemeVq38HDv8Qf0buGsnJsJMAdCrViTHVx1AuVzngRSvAq7vw2BmEe5W2oklCnTowfz40bvzWMaQky9oY7t8HP79Xv06dgoQE8TVXV3EEefhw8d/y5d9e+yM5H2aMQQZgSZLMQsC9AIpnL46Dxk7MuHbuFLWCBw+GYcNEVyLppfdVuVJVlZNhJ9HpddT/rT5xujg0FhrWtV9H2+JtX73gwQOxx5unCgwuCa1aicYIFSokaRytyzunScDV60VRraNHXwXcmy+2oG1soFIlGDoUqlcXn9E+ZCn5Qz7MGJMMwJIkmV58PPdDr1CpkLvYbyxWTGQ29+8PWbKYenRm6W1VrnJntsYnyIepR6dyPPQ4fcv3JUGfgIqKqqpceXxFRLNt28QKw5UrcPu22Ec/d04kspmBuDiRIHXkCBw+LAJuVJT4Ws6c4oTZoEEi4FaoIIJwcqVGyc6kkAFYkiTTefIEFi6E2bM5ee8eCefXi8e9vU07rnQgcf81Un+BOIvz2BiKo1iGEWy1nXbrb1EgSwHmNZ1HcafirDm/5tUe7+lw+LyYqAqWLx/873+vLmrC4BsRIWa3iQE3IODVcnLJkqJ0d40a4lf+/MY93m2skp0fSgZgSZLS3qNHoiPO0qWvNWK3LlnW1CNLN1qXd+byk1OMP/ItBlWLhWKFnaUtxTIVYmaTKbQt3haNhaj8ta/nPnGO96kj7i0GiQykdeugbVuTdSN6+FAE2oMHxXHuc+fE0W4rK7GcPHw4fPaZmOFmM3Jfjn8zVTKZDMCSJKWdp0/F2Q9LS1i9Gjp0gJEjmRlzgOOhK1hDA2TdqqS58/QO24KnYCAeFFAUHUOqDmBSvUkoidPDCxfA2xv3LFlwnzZNRLiAKmLNNo0rhIWGikCbGHCDgsTj9vYiyP70E9SqJT4bpPWpMlMlk8kALElS6tLrxXnSadPEbPfUKbGvm9iQHdi5+mvuRd97FTiktzr/4Dyefp6svbAWg8GARhGzXGuNNa2KthIfYPbsET/vXbtENBs8WLxYUaBixTQZ5/37oi+Dr6/4b2Lb5UyZxMy2Vy+xzV+hgnlsO6dVMtk/yQAsSVLq+GdHouvXwc1NrCvq9WIG/CL4qqpKQFgArYu1NuVozdY/W/2pqNRYVoMMVhkYUnkIw6sNJ+xZ2KtWgC7u8P33Ynk/Vy745RcYMCD113ARS8q+vq8C7uXL4vFMmcTMdsAAEXDLlfvkemK8lQzAkiSljvXrRZpqlSri923avHG/8fbT24THhhuvB7CRJLfKlDEduXOEeqtEIwRbS1v29NjDrMaz6FamG1ntsgKQz5AR9yO+YG8NLkC3biJLqWvXlKUGv0dUlFhO3rcP9u8X5RxBdCKsWVMUvKhTRwRcE20zmz35Y5EkyTguXRKz3QoVRODt0kWUGnpPR6KAMNEhzZwCsKkKMySK1cay8uxKxu8bT4JepAIn6BM4dPsQ39T8Rjzpxg2YPh2WLxerDaoqlpeLFhW/jCwuTlSW2rdP/Dp5Uixm2NqKJeUuXUTArVhRJFJJ7/feAKwoyjKgOfBQVdVSLx7zBFoACcANoI+qqpGpOE5JksyRqor1Ri8vUTjD1laUIAIx+/rss/deQqNoqOpcldI5SqfyYJPOVIUZADZd2sTAHQN5FPOI4k7FidZGozfoxREiNw/xpL59YdkyMbXs1k0Uziht3J+fwSAyk3fvFlvKR46IIKzRiJrJ48ZBvXrg7v72ClPSuyVlBrwCmAOs+sdje4BvVFXVKYoyBfgGGGv84UmSZNYGDIBFi0Rx/p9/hoEDRXeiD9CmeBvaFG+TSgNMnrQszOAf4s+Wy1uo7lKdVsVakTdTXirlqcTYGmOpla8Wx+4ew/fmfjzCrHF3flGGs3Bh+OYbUR87d26jjeXuXRFs9+yBvXvFaTGAEiXgyy9FwK1VSySyG5M5LPebwnsDsKqqhxRFcfvXY7v/8cdjQHsjj0uSJHMUGQmLF4tZV548Yt2xcmXo3j1Z0yBVVTGohpfnVc1FWhVm+O3sb/TZ2ge9qkejaDjc5zDuLu781e0v8YRnz3D3OYH7jCUQHAxZy0CjRjDWOPOd58/FPu6uXWKmm5g4lTOnOJrdoAHUrw/OqRgLTb3cb0rGaCXyObDzbV9UFKW/oigBiqIEPEr8OCVJUvoSHCy64bi4wNdfw/bt4nEPD7Ecmsw1yOtPrpN5cma2XdlmtKEaw5hGRbGzev1DgTELMxwMPkiTNU3ouaUnevXVUrdvsK/4TXS0CLIuLiJz3NkZfHxENEwBVRXLyp6e4lJZs0KzZuIzVb58Yifh7FnRwm/1anFUKDWDL7x7uf9jl6IkLEVRxgM6YM3bnqOq6iJgEUClSpXUlNxPkqQ0ZjCI2e26daJvW6dOMGqUaC1jBAFhATzXPsc1s6tRrmcsqVGYQW/Qv5zpr7u4jtP3TvNlxS9ZdXbVqzKRWV78XG1tYeNGMQ0dNSpFHaCePBGz28RZbliYeLxUKfjqKzGhrlnTdPu4pqrD/BqdziSp2sm+o6IovRDJWfVUVZWBVZI+FgYDHD8usmssLMR53VGjRKuZvHmNequAsABsLW0pkb2EUa9rDMYozOAf4s+em3t4Fv+MzZc3s7L1Smq41uCXur/g3cgbW0tbepXpge/uRXhsv4D73M/h1i2RwHbhQrJKQhkMoifuzp3w11/ir9JgELVPGjQQ3QUbNkz9mW1SmaoOMyCm+vPnizyGw4fF3noaSlYAVhSlMSLpqraqqjHGHZIkSSaRWDhj+nRRtujSJSheXKxPppKAewGUy1UOK83Hd25l1/VdNF/bHJ1BB0Bxp+Ivv5bVLqtIKV6xBHdvb9yDgkRETCxUAh8UfCMixOx25074+2948ECc/KpUCb79Fpo0EVv15lgAwyR1mE+fhhkz4I8/xOy3eXPQalPvfm+RlGNIawEPwElRlLvAD4isZxtgz4vSccdUVR2QiuOUJCm1RESIoDtvHoSHi3fqP/5I9dmA3qDn9L3T9C7bO1Xv8zapmXlrUA103tj5ZfC1UCzoXqY7NVxrvHqSnx/06ycqVSTWxU5iTUZVFb1xd+wQ2/F+fmKWmzWrWFJu0kT8N0cOo3w7qSrN6zBHRori05aWIov/q6/SfOabSEnL1eNKlSqpAQEBaXY/SZLeIS5ObPw9fizKRNavL5aaP/ssTQr1x2hjmHp0KjVda1KvQL1Uv98//TvzFsSsa1Lb0sl+4z9z/wwrz6zEq6EXGgsNk49M5qeDP6HVa7HWWLPPYxnuqw6IspATJ4oo6u8vlvqT8POOixNHrrdvF4H39m3xeIUK0LSp+FWlinnOck0qKkqcmT5xAn7/XTy2a5fYV3d0TPXbK4pySlXVN1aZkQFYkj4lqirqBk6bJrJz/P3Fm394eJrUCzYXNSbvf+O+o7OjHUfH1U3ydfzu+LE0cCkXHl3gROgJHKwd8Pvcj9I5RVEM/zt++Poux+OvINzXHRV7u4MHi59/EoSFiYC7fbuoPhUTI7oHNWggspebNjWfvVyzc+MGzJ4tgu+zZ+KD5Y4dojh1GnpXAJalKCXpU5CQIJaVvb3FOZOcOUURh8TGCCYIvtefXCe7fXYy2xq5qkMSGCPzduvlrbRZ1wYVMYkZUHEAE+tNJItdlpfPcV/yN+4TlojiJD/8IEp0vmNdWFXFX8+2baKBVOJ8xc1N1FZu3lw0NJCVp97jr7/ED0ujgc6dYdgwsSFuZmQAlqRPwerV8MUXoqTR0qWiUL+J38V7bu6JpYUlh/ocSvN7JzfzNkYbw6VHl6iUpxIXHl54GXw1igbXzK5kSbCAuV4izbhMGXFsK29e6NHjjUlVWwJDmbLjKrfO2aOEOKMLzsXj+5YoilghnTgRWrSAkiXTvH1v+hIfLz5gZswIbduKTynffw/9+4uCMWZKBmBJ+hjdvCmyPMuVE1OnLl3EG1GjRmbxTq4z6Ai8H8jASgNNcv8PzbwNjwln7sm5zD4xG4CQESHUzV8Xu8N24gyvhRUef56FZi5iuTM+ni36bHjuekBYpDN5Zvq/llgUGQkT5j9h8WoLoq/XRE2wRLHSkaFAOIP7W/PdoCzkzJnqP4b07+FDWLBAJBA+eACtW4sAnCED/PijqUf3XjIAS9LHxN9f7C9u3izO8H79tXjczk4cADUTAWEBxOniqJynsknun9TM2y1BW5jqN5XT904Tr4+neZHmjK0xFhuNDe4u7uzruQ/fGcPw2HgK99CNLwuVbFFy/qe84ugVV9m5PgO3Tjly4ADodFmxyBBHhuJh2BV+gK3rYyysDATa2ZEzZ9L3oT9Znp7w3Xdi9tu0qTjClYJKYVq9lgPBB2hYsKHxxvgeMgBL0sdi8GAxE3B0hDFjxPEKM8zQmXFsBl/v+Rp7K3tq5qtpsnG8q9CG3qDnROgJOm/qTLw+Ho2iYXWb1XQr003sm+/eDQ0b4u7ijnumZtCxtihU4uICgOfk/cRq9WgfOxBzLScx13KRcM+RG0CRIiLZfGXoUazzRP5nQSJNK0ClJwaD2NutVk3sqRcqBH36iP3dYsVSdOld13fRb1s/QqJCODvgLGVyljHSoN9NBmBJSq+io0Uv2E6dRGJPq1bijahPH9EV3YycDD1J3kx5yZ0xN6VzlKZfhX6Mqj6KvJmMW1krJVRVZd+tfUw9OpUi2YrgnNH55TlegDuPr4sPONOnw/Xrol1QvXoiuerlNUSNh4t/uvL8Si50T8Tfg3XuCBxrX8a+8H2uLPYA4PDkeEIj/zuONKkAlZ5ER4sCMTNnigIxXl7iE0ybNuJXMoVGhaI1aHFzdMM5kzMFsxZkfrP5lMpRyoiDfzcZgCUpvQkNFccrFi4Um4nW1qJXXMOG4peZSAxok49MZt+tfYytMZbJ9SdTr0C9ND/3+zb+If7sv7Ufvapny+UtBN4PJJdDLpoVbkYV5ypYa6zFHq9BweOraXDp2atCJbVrA2JC7OcneiX4+MCdO4BFAWxdnpCpYjB2he9jmTEeEMecEpmkAlR6oqqi6fDChfD0qchKW7sW2rVL0WXPPziPl78Xv5//nfYl2rO23VpK5SjFgV4HjDTwpJMBWJLSC71eZDKvWSOW49q0ETMBd3dTj+w/fIJ8mHh4IqfunSK3Q248G3jSv2J/Uw/rNf4h/tRbVY84XRwqKq6ZXFnSYgndy3THxtIGnj4Ve7w39uExahbuhavDAlGoRKdX8PUV/RK2bBH5PzY24vPPTz+BRb4HTNp/9p3BNc0rQKUHqgpXr0LRoiJZ8MoVkTg4YoRYek6Bg8EHmXRkErtu7CKDVQYGVRrE8GrDjTPuZJIBWJLMmcEg1jQrVRJnGnU6cZZ02DAoUMDUo3uNzqDD0kK8pWy8tJGo+CgWt1hMjzI9REAzE49jHjP3xFwexTwiQZ+AioqFYsGXlb7ki/Kfi6XladPg0iXcb9zA3cUdDg5Ha+PA/v2wsb/IcQsPF8m2TZuKSVnTpuIUjJCbTI6G9wZXYzR8+ChotbBpk1jeP3lSLPEXKCAeS0FpL61ei6WFJYqi8Ne1vzhz/wy/1v2VAZUGiHrcpqaqapr9qlixoipJUhLExKjqokWqWry4qiqKql69auoRvdXTuKfq1CNT1TzT8qjn7p9TVVVVn8Q8UXV6nYlH9rqbT26qQ3YMUe1+sVP5EXXg9oGq3S92quYnjWr3i53qt+BbVS1TRlVBVXPlUtVfflHjnjxXt29X1d69VTVLFvGljBlVtWtXVd28Wfw1SSnw9KmqTp6sqnnzih9u4cKqOmeOqkZHp+yycU/VaX7TVBdvF3XntZ0vH4vVxhpj1B8ECFDfEhPlDFiSzElkpDi/O28ePHokzvGuWiW6pZuZB9EPmHl8JvNOzuNp/FPqF6iPQTUAvFYNylT8Q/zxDfbFw82DVWdXsej0IjSKhh5lejCq+ihKZC9BjzI9xHMe2OLeYSSUKoVu0TL25ujKH5tt2JJfbD9mzixy3Nq3F2UgZSWqFEpIELkLz5+Lo0Q1a4q2gE2biuNzyRQaFcqs47NYeGohT+OfUjtfbTLZiNKTif81J7IWtCSZg/h4sYn44IGoO1ivntjf9fAwi8IZ/xaviyfv9LyEx4TTvkR7xtYYS8U8FU09rJcm793K/452QFX1WChWNMjXkbLOuRladSjOmZxFneAZM16WiNQmqATOOMiCoNps3qIQGSmCbps2oklR/fpJblQkvY2qioLWM2aIotb794vHQ0ONclxOVVUKzy7MrchbtC/RntHuo6nsbJpz5v8ka0FLkjlKbIzg7S2qJx06JGo0375tln3kztw/w7oL65hYbyI2ljbMazqPsrnKUiRbEVMP7SWdQcfYHQuZeeo7VEULChhULadvahhQdijOF27DtKGweTOqpSUhzQYyoR/4+Cg8eeJBpkxiptupk5jpyqBrBHFxogvRjBlw/rz4tz1okEgq1GiSHXxVVeVA8AGWBS5jWatlWGusWdRiEW6ObhTIYl75EW8jA7AkpbV/N0bIkUMU0Uh8QzKj4KuqKgdvH2TykcnsurGLjNYZ6VuhLwWzFqRDyQ6mHt5L8bp4Fp1ahPcxb4Ijg9Go2RFvbwYULNHoShIxehzsX402Yxb2V/iGMcGDOb8lDw4O0LKlCLoNG8rlZaNbulQ0/ihTRpxb79JFrPYkk1avZcOlDXj5eRF4P5CcGXJy5fEVSucsTd386auCmAzAkpTWli8XjcBLlIAlS6BbN7N81w+ODKbzxs4cDz1Ojgw5mFh3IgMrD8TR1tHUQwNeneGtm78ulfJUYpr/NPJmykvMgx7YGiqj6M9R8NFWouxqcTe2KksfFOWGYzWmR/bGcDEDLVrAD53EtuMb+iRIyXXmjJjt1qsnmlD07Cn+rRthO+Vu1F1qLKvBnad3KOZUjCUtltCtTDdsLc3v/5+kkAFYklLb9euiik/lyuLNqFs3kVRlJo0R/ilBn8CNJzconr04uR1yY2lhyfxm8+lVthd2VuYTpTZc3EDXTV3RqTrsDtuxr+c+TvY7SfYM2Wn9zToa+66kS+DfZE54zvd2jZkQW5MHFgYcm1iwqIuY8b46MiSlmMEgmhZPnw6+vqJpcakXFaUyZoQ6dZJ96bBnYZy5f4amhZvinNGZBgUa0KpoK5oVaYaFkvyELXMgA7AkpQZVhaNHxTLzli2v99x1cDCrxggAz+KfsejUIqYfm46lhSXXvrqGjaUNRz4/YvR7bQkMTXbxiVNhp5jqN5UNFze8bAWYoE/AN9gXdxd3nvcaxMbVi1EMBjbRDm9GcDZHYXKVvMDUMVnp4WG+renStY4dxZldFxeYOhX69oUsKcuEv/DwAtP8p7Hm3BocrB0IGxWGraUtS1ouMdKgTU8GYElKDf37i+XlrFnhf/8Te7y5c5t6VP/x6PkjZh2fxZyTc4iMi6SOWx1q5vqC2lMPce9pnNGrM20JDP1Pl6BvfM4DvPceAWEBVF5cmUw2mehWuhsbgzai1WuxxpJn52tT70eov98ROwaztdBAgosq6NzuUCHvuRffg3GDb0o+SKR7d+6Io3Jjx4pA27+/SBdv2xasrFJ06XMPzjFu7zh2Xt+JvZU9X1b8khHuI9LtMvO7yAAsScYQFSWSTXr2FDPdtm2hfHno1UuUSzIzqqqiKAr+d/359fCvtCnehrE1xhL20PlFgIwDPixAJoXnriuvlWcEiNXq8dx15T/XP3z7MHNOzCF3xtzMaDyDirkrsrTlUtqXaI9NnDX1DuYl6OJSWl95xLi7WsIKQfwPE2nbBYa/rPiYOhnaKfkgka4dOyaWmTdtEn+uVk304E1hDXKdQUdUfBRZ7bKiN+g5de8UE+pMYGClgWSzz5bycZspGYAlKSVu34ZZs2DxYnGUKGtWEXSbNDH1yN7o7P2zTDk6hcJZC/NTnZ9oXqQ5V4ZcoXC2wgDU+G1/kgNkcryt1d4/H49OiObb/d8y6/isl2UiO5ToQHWXGhR71BHfz7ypdmouvdWHnLcsx5l63kxb505F97TbUv+QDxIfhdhYkVTl7y8OSI8YIdpdurqm6LLRCdEsPb2U6cemUytfLVa1WUX53OUJGRGCtebjPwMmA7AkJYdOB927i2r8IM6wjBghajabmcSjRFOOTuHv63/jYO3A19W/BsBCsXgZfCFpATIl8jjaEfqGayW24Nt0aRP9tvUjIi7i5dcUFH5ctJcbv9cg7JYlwcwnLE8lbg8dRfmRdShtlfaJbKn9czILkZFw+DC0aCHSxEuXhq5doXfvFLe7vPfsHrNPzGZ+wHwi4yKp6VqTjiU7vvz6pxB8AdJ3CpkkpSW9Ho4fF7+3tBQl80aOhFu3RIciMwy+AOP3j6fOyjqcvneaX+v+yp3hd/iu9ndvfO7betEaq0ftmEZFsbN6vbi+xuoBPT4T+3uFsxWmtlttpnksxAo7FIMFmgSV/isXUaSAjsWrbHG4e5lyoTuoPLYuliYIvpD6PyeTun5dzG7z5hVbKQ8fiscXLhTneY3Qa9rTz5PJRyZTL389jn1xjEN9DtG8SPMUXze9kTNgSXqfxMb3M2ZAcLAoY+jmJqr7mKF4XTxrzq/hM9fPKJKtCB1KdMA1s2uSjhKldo/axOXZ73f6EByzC8U6lGj1HEcfdeXzhOrc9C+D5vf1BO7YyILcrjxwu0K1B5mp1KEvHX5JAHtLILNRxpISH2Uv32vXYPRo2LZNfMDs2hWGD09xYRhVVTl0+xCefp6MqDaCegXq8XWNrxlceTAFsxY0ztjTKRmAJeltHj8GL69Xje/d3WHyZDEzMENR8VEvjxKFPQvj+1rf81Odnyifuzzlc5dP0jXSokdtSIIPF3XDMVgZQIWmzj2w9p2Cc3/xI+/muJvf9F2JsyyKbY+FopiDmVXK+Gh6+cbHi6YfefOKn/HJk/Dtt6JUZK5cKbq0zqDDJ8gHLz8vToadxMneiYfPxWw6l0PKrv2xkM0YJOnfYmJEIYGwMChYUOyBjRhhlo3vE004OIFp/tN4Gv+UuvnrMrbGWBoUaIBiJoU+/tmXtebymhy58+J8sUFD5v0j+clPh1PJHDhOGkejBgYsD+wRxZhT0BlHeodHj2DBAnGUqFQp2LNHPK7TidmvEdRYVgO/ED8KZy3MKPdR9Czb06yKuaQV2YxBkt7HYIAdO0QjdgsL0SQhTx7RqSWrGTTufoOQpyHkzZQXRVF4FPOIBgUb8HX1r82iA0xiK8AqzlU4//A83v7ezG+yhPhLDdHu+QEKtASLBGwMsOPONKorCkrtgdAMwEJUCZOM7/JlURzmt99Ek4TGjcWHy0QpCL4Poh+wLHAZo6uPxkpjxeDKgxntPpqWRVuisdC8/wKfIBmApU/b8+ei3+706WIPzMUFhg0TAdnCwiyD7+l7p5lydAobL23Et5cvNfPVZGbjmWYz2/UP8afuqrrE6+JfVqty1taiW4eMPL0AuXPXZ1ZkZ6KfLcfjkT3uHQfB0KHiZy8Zn6qKf88ajdjf/e03cV592DBRozmFLj++jLe/N6vOriJBn0B1l+rUdqtN19JdjTD4j5sMwNKnbckSkWhSubLoUGSESj6pQVVV9t3ax5SjU9h7cy+ZbDIxpvqYl0eIzCX4AhwIPkCcThTyQAUCBvBs1zS8Sq+kyKJs1OhTBMszg+BQKVGyMJP5NUr/KMTGwurVInlw3Dixlz5gAPTpI/ogJ9HbKn5FxEbQe2tv/rzyJ7aWtvQu15uR7iPNqj2luZMBWPq0nD0rluAaNBDnePv0gYoVoUYNs2uM8E9xuji6bOqCpYUlU+pP4cuKX5LZ1vTZwIlOhJ5g8akltLOdy4G1dSCvDVjosDBYscg2jj6ZXLEIDIeHv4DleHFky0yPbaV79++Lvd3580VWW7lyr+qQZ8z4QV0o/l3x625kNCN9tgPNaVkuNxGxEXxf63sGVxlMjgzm00YzvZABWPr4GQzw998i8O7bJ0pDlikjvpYpE3z2mWnH9wax2lhWnl3Jlstb2NF1B3ZWduzuvpvi2Yu/syZuWtUn9g/x50DwAWw0Nmw8v41j9w9iEe/IkuWDyBLvTrteByiWMJJmOwNwD14putyPGiU+6Eipq0ULOHXqVfJg7drJ/nCZWPHLQBzPNfuIstyCQYli8t/OtC7vzMHeB81q9SW9kQFY+vh17gwbNoCzszhG1L9/iju1pJaI2AjmB8xn5vGZPHz+kCrOVXjw/AF5MuZ571GitKpPnLjH+3KZOTo7HPGmZsYv+N+Ai0Q2vMvMg7E09slFWM6m7Jk+mgatahrt/tI/GAzw11/iqNyaNeID5axZYom5cOH3v/49QiLvE2W5nWeWOzAoUVgbiuKY0Jv7cQbAvLY+0iMZgKWPz/374ojF8OHg6CiWmVu3Ft1azHB/N9H5B+epvqw60QnRNCrYiLE1xuLh5pHkN7nUrk/8LP4Z206eYcHOI8RZakUdPYMFdWwHsa59NrL/Vgt+PEvXW16E5irGL3X7AmAX8JxJrqHp74ysOXv+HFauFH2mr14V53ivXhXL+kY4LmdQDVgoFjhmiiBEuxY7fVUy6dpgYyiJgoLzx1Dxywy8NwArirIMaA48VFW11IvHsgLrADcgGOioqmrE264hSWni/HmRzbxmDWi1Ypm5bVuzbYwAEPQoiGtPrtGyaEtKZC9B3/J96VWuF+Vylfvga6VWfeLgx/cZ9vss/no4H53egMWGzWi6WaOSgI2FBb/+OZfsZx5DyZL82m4MAU6vVzf6qJsUmMLDh1CsGEREQJUqsHYttGtnlA+X/iH+ePp5kt0+OwtbLOTnJq0Z5eOAXvsqaSvdV/wyI0k55b4C+Hf38HHAPlVVCwP7XvxZkkwjLk60QytTBtatg3794MoVEXzNlF+IH63+aEWJeSUY8tcQ9AY9GgsN0xtPT1bwBePWJ/YP8efztV9TcHwr8s/Kx59PJmMTVo+BGXezfGZ+SthMIVtcJ7astqKgbQmxx37+PEsK1SbB8r+B4KNqUpCGtgSGUmPyflr0nsnUNiPYEhgqSkMOHQpHj4r2gJ07pyj4GlQDWy9v5bNln1F9WXV8g33Jm0lUe2td3plpbevj7GiHAjg72jGpbWn5YcpI3jsDVlX1kKIobv96uBXg8eL3KwFfYKwxByZJ7xQbK1qj1a0LtraQMydMmiT2d83w7G6i43ePM3rPaI7cOUJWu6x8X+t7hlQZYpRCBcaoTxwbCxNWHWJKaGMMSjxYGXB+1opf63jSo/4j7v/wE3He5/nhi3nYKwX4umlz4jNlYVKOUrRWlPd2O5KSbkvAHfZOXYK3vw9V717kiV0m6hWtDUDrH3802n1+8v2Jnw/9jJujGzMbz+Tz8p/jYP2q4ULr8s4y4KaS5O4B51RV9R6Aqqr3FEV5a/65oij9gf4ArinsHSlJrx2xiIyEO3cgd25RXMBMJegTiNHG4GjriNag5c7TO8xoNIMvKnzx2htdSqWkPvGlIAP/W7aTHZFT0VlEQ64EsDCgUTQMLpSBXr/2hGPHcLDLyOayjbHRJRBnZUukXSb4xxLzR9mkwBQOH6Zy6060fnKPu5lyMKHOF6wr24horFK8nB8eE878gPnUy18Pdxd3epfrTfHsxWlfoj2WFjItKC2l+k9bVdVFwCIQtaBT+37SRyokBH744dX+bosWohVgCgvGp6bohGgWn1qM9zFvWhRpwbxm8/jM9TNuDL2Ram90HzJbiY+HDT4J/LJ1LVeyeUKOi9jbuNDGrQPbHwSRoE/AGg0eE38H6wIwezbuN3Px3Pq/s9nEJeaPpkmBKdy5I+qQFysGLi6E2Wfhl5q92F3EHf0/VkiSu5x/K+IW3v7eLDuzjBhtDAbVgLuLO/mz5Cd/lvzG+i6kD5Dcd4EHiqLkfjH7zQ08NOagJAkQRywiIkQRAVWFTZtE5aRhw6CI+Vbbefj8IbOOz2LeyXlExEVQO19tWhVt9fLrpp5lbDzuz+w/fTn3pweRObdBzUnkUkrzbZ3f6O/2GVZzF+CfrR++1XLh4VoT9/IR0LQpaDQ4Tt7P8/csMcslyw904oQ4o75xo0gY3LYN3NwYPni20Zbzh/w1hPkB89EoGrqV6cYo91GUylHKGKOXUiC57wR/Ar2AyS/+u9VoI5KkxBJ606eDq6tI8HF1hXv3RJciM/fd/u9YfHoxbYq34evqX1M1b1VTDwm9XhwX/fH37Zwu3AYsDVi0teG7YmuoVnEnTWJyoXh7w9o+YDDg/uWXuI/5Rrw436vryCVmI9q5E375Bfz8xPndESPgq69efjklP2uDamDvzb3Uy18PjYWGAlkKMNp9NEOrDsU5k/xwZC6ScgxpLSLhyklRlLvAD4jAu15RlC+AO0CH1Byk9Il48EDs786bJ0rolS8vykWqqqjkY6bB91TYKab6TWVktZFUzVuV8bXGM9J9JEWdTB+UHjyApUthztor3MvvBeWXg6IHBRQlAbu8l2m6Jg5+/llUCBs8WKww5H/zkqRcYk6hqCjx79jSEgIDxYfKGTPg88//UyIyOT/reF08v5//HS9/Ly49usTmTptpXaw1I91HpuZ3JSWT7AcsmV5igJ0yBb755tX+bq1aZlufWVVV9tzcw9SjU9l3ax+ZbDIxr+k8upXpluxrGquMpKrCkSMiT23jRtA2+hIqLsbawoYmhZuw68ZOtPoErDXW7Ou1H/e7wOHDIoPc0THZ45feIThYVKhasgQWL4ZOncTxOSsr0aUoheJ18Uw/Np1Zx2dxL/oeZXOWZXT10XQq2QkrjfkWn/kUyH7AkvkxGGDXLrH31acPdO0KX34JbdqY9f4uiODb4LcG7Lu1j9wOuY3SHMEYZSSfPxc5alN+P8pNx8XYX/2CgQNr4lCvFJqM3zKkYBdyrNyI/wZffDPH4dGoJ+4u7uCCUaonSW/g7y+2UjZtEu0tO3R41QLQ9u01vZMqRhuDvZU9lhaWLAtcRqkcpVjZeiX1C9SXZSLTARmApbQVGyuODM2YAUFBoum9/sUel6Oj2c7AYrQxrL+4np5le2KhWNC2eFu6lOpC9zLdsbG0SfH1U1JG8vp1sWq/dEUCUeV+Bo+JoKjoKq6lc29f3F2+EvuLCytCbCzuTZviPmoU1KmT4nEnR4w2hpVnVnI05CgLmi9I1lGstGo6kSIGg/hw+eABjBkDQ4aIkpFGcOb+Gbz8vNh3ax83ht7A3sqegP4BZLKRrR3TExmApbTVpAkcPCj2d3/7DTp2BGtrU4/qrcJjwplzYg6zT8wmPDacfJnzUSd/HQZVHmTU+3xoGcnEBk9z5sDO3QlYuM/CZtAMsAp9+Ry9QY9vsK+Y5SYkQLduIhAboQl7SrRY24L9t/YD0K54O9oUb/NBr0+rphMf7OlTscS8erVY0ndwEHsAbm7i9ymUuO3h6efJ3pt7cbB2oH+F/sTr4rG3spfBNx1KSilKSUq+s2fF0vKzZ+LP48eDr69ol9a9u9kG36j4KIbuHIrrDFd+PPgj1V2qc7jPYTzcPJJ9zcSygvnH7aDG5P2irOALSS0jGRkpVu2LFIFmreIIDITvv7Mkf7sluBcpind9T+wUazQGsNbq8Yh9USNnzhyx92iC4Hs1/CpD/hpCZFwkAN/W/Ja9PfZiZ2mHb7DvB1/vXasFJnHrlmj8kTcvjB4tMpofvjiZWaqUUYIvQEBYAI1WN+Liw4tMrjeZkBEhTGs0jSx25tnZS3o/OQOWjM9gEEcsvL1h/36R9dmlC3h4QIMGph7dO0XERpDFLgv2Vvb8ff1vOpbsyGj30ZTMUTJF133frO19R06CgmD2bNEAJ6bAWjK09CZj1ltcHB5MVgcHRkYcIPOy3+HbmVRTE/Ct5IRH3c9xr9lVXCyN9wNVVeVoyFG8/Lz488qfWGusaV6kOY0LNaZOfrH0XcO1Br63fT/42qnVdCJZrl+HokXF/m6nTmKFoWJFo1z6WfwzFp9ezLP4Z/zg8QOV8lRic6fNNCnUxCjbHpLpyQAsGVdkJFSrJpohpIP+uyCCxYHgA0w9OpWzD85ya9gtbC1tuTDoAtYa48zQ37fH+6YjJ6MaFMUy1JmGY2HPHrAq4EfWoeOIsT3Mc8DSYMnxkP00Kd6SzFYZxZnSsmVxHzUX92bNRFAwgRhtDPVW1ePY3WNktcvKt7W+ZXDlweR0yPna8zwbeJLBKsMHX9+k9aZ1OvDxEVnNX38NhQqJ7ObWrcW/dyMIexbGrOOzWBCwgKfxT2lcqDGqqqIoCq2LtTbKPSTzIAOwlHJhYaKYQPv2IomqXj1RNrJ9e7Puv6s36NkUtImpR6dy6t4pcmbIydCqQ9EbRKA0VvCFpM3aEgPx06ewfDl83Qlu3BDv61/9co7ZuhpEWdqh6BRUVFSDnjPf96fJuuZimTMoyGSlOZ8nPMcvxI8GBRtgb2VPqeyl6F66O73L9SaD9ZuDbHK7PpmkGEji/u6sWaJkZKlS4qicpaU4O20kK8+spN+2fuhVPe2Kt2N09dFUca5itOtL5kUGYCn5AgPFEYs//hBnGRs0gMyZYe5cU48sSY6GHKXTxk4UzlqYRc0X0aNsD2wt3340JCWZt0mZtV27Bv+b78+fZ31JuF6dQlVv0GnkPX7rNx4rqzLUOP872c9cp/m1H0lAxdqg4lGgnjhPam9vkuD7IPoBc0/OZe7JuUTFRxEyIoRcDrlY3HJxkl7/+/nfMagGupfpnuR7pnkxkE2boHdviI6G2rXFXkCzZkY5v6uqKoduHyKbfTZK5ShF1bxV6V+xPyOqjaBg1oLvv4CUrslCHNKHO3dOVEvy9RUzr88/F/1JC5r3G8aT2CfMOzkPg2rg+9rfv8wqTSzX9y7/3sMFMetKam/Ut71+YpvSZI50ZsYM2HbGH3rWA8s4UMT/l9XyVuNInyNifD4+0K4d24tnYWblgkTm78j4Vl0/KPAY6/hOaFQoPx/8mZVnV5KgT6Bl0ZaMrj6aGi41Puj8acPfGnI/+j7nBp774DGkGlUV53cdHESP6WvX4KefxIy3QgWj3EJn0OET5IOnnycBYQH0KtuLFa1XGOXaknmRhTiklHv+HJ48ARcXUbIwOBg8PUVzBDM9u5voduRtph+bzpLTS3iufU7Hkh1f7qk1LNgwSddIyTld+O+sLVeGDFTUleWH3lk4fx6yZ4fSA+dy3uLVLLlP0c4svVQIZf4CGDyYrXkrsK/deHYUqCK648TxQcdvUnp8R1VVohOiyWiTEb2q5/cLv9O7XG9GVBuR7LKbHm4ejN8/nscxj3Gyd0rWNYwmcX/X2xuOHxeJg7//DoULi6NFRrLizAp+PvgztyJvUShrIeY3m0+vsr2Mdn0p/ZABWHq3u3fFEZZFi6BGDdGppWBBsTlpoiSfD7EwYCGD/xqMoih0Ld2V0e6jKZ2z9AdfxxiZt63LO1MttzPz58P82XDsERSpEcQvC2FUz+Jsu9GKLpv+AFXF2qDQ738bUYL1IokNmLr/JqGFXq9Y9SEfApL7IUJv0LPl8hY8/TzJbJuZXd134ZrZlfuj7r91fzepEo91HQw+SLsS7VJ0rRRZsgQmTBD7u4UKiW2UXsYLio+ePyKbfTYsFAuuP7lOLodcTGs4jZZFW7539UX6eJn/O6hkGqdPi/KQ+fOLmW7duqJOcyIzDb6qqnLg1gGCHgUB4qjLsKrDuDn0Jitbr0xW8IWkn9N9m/PnxUq9S3V/fj4wiVzNF1J9diuuNihBoON32NpCh5IdOGzoxYS9Bvat0eDeuB9cvgwLFgAp/xDwoa+P0cYw7+Q8is4pSvsN7Xkc85hWRVuRuG2V0uALUDlPZeyt7JN1HjjFgoNfVWG7e1f8W9+6VfzMBw0SKz0pdC38GgO2D8B1hivbr24H4EePH/H7wo82xdvI4PuJkzNg6ZXENyONBnbsgO3bRXu0r756a3ccc6E36PEJ8mGq31QCwgLoV6Efi1osolSOUkxrNO29r3/f3mhyMm9VFXbvhmnTxDEi60J+GHrWBSWe80CmZ5n48bPvGHwnp5h5ubriXrcX7pb5YcAAcHp9STalx28+9PXzTs5jzJ4xVHWuypT6U2hdrLXRA4aVxorPXD8jJCrEqNd9J39/sczs4yMqVbVpA999Bz/+aLxbhPjj6efJlstbsNJY0bNMT0pkF0VQTN0PWjIfMglLElWqli2DmTNFR6IOHUTbNBBVfczcyjMrmXBoAjciblA4a2HGVB/z3ozmf0pqglVSE5ji4sTWobc3XLwIuZzjGTLIkriKU/j12LeoqFig8L2mLj/MDxLHuH79Ff73P6OMM7mvvxp+FW9/b+rlr0eHkh2IiI3g4qOLH5xY9aHidfGpX1hCr3+1v3vsmMhbGDBAfLjMk8e4tzLoKTS7EE/jnjKo8iCGVBlCLgfTHA+TTE8mYUlvdvv2qxZpUVFijzd7dvE1Mw+8T2KfkMU2C4qicPHRRZzsnfBs4JmsPbWk7o3+s2DGmzx+LFoAzpkjKhGWqhRJx5kLOaydSYmmc8nlUIdpAbYkaOOw1qk0XLIPStQXDXsbNXrvOFN6/OZtr8+eLZg264aw9fJWrDXW5HcUqx1Z7LLwmetnSbp2SqRq8NXrXx0X+uYbURFszhyxv2ukEpFxujhWn1vN6nOr+bv739ha2rKl0xYKZi2YrEYT0qdDzoA/VaoqiglcuSJmvCNGQBXzP/AfHBmMt783SwOXsqnjJhoXakyCPgErC6tkz9Lyj9vBm/4vUIBbk5u99/U3bojj0MuWQWw2fwo0/RO3Mnc4GbWNZwnPaFCgARPy9aZqra74h/jj6/0VHs9z4D54EpQtm6wxG0vvLb1ZeXYlWe2yMqiSmK39u2JVWuju0x03Rzd+qfuLcS6Y+OFyyxa4cAHs7MSer4uLUc7vgihbOj9gPrOOz+LB8weUz1WeDR02yPO70mvkDFgSRyw2bRJRwsdHJJgsWSIKyLu4mHp07xV4LxBPP0/WX1yPhWJBtzLdKJClAJDyilXJ3Vs9fhy8vMSP09ISGn7uz17netzUx3LzMTTIX58pFg0pP28z+HeDwBK4l3PH3ftkmtdmThSrjeW3c7/RpVQXMtpkpGXRllTKU4k+5foYJakquR7HPObM/TMpD8DHj4tl5k2bxJ87dhSrO3Z2oiuRkdx4coOyC8ryXPucRgUbMab6GOrmryt78EofRAbgj11kpOiCM3s2hISIIxa3bonZbzppwq436Gn5R0uexj1lRLURDKs2jLyZjNNXFT4swcpgELlpXl6i41ymzCpdxh4hpuQCyuUtys6DCQBoUKiz8RTlt+2FAgXEz79QIXERE7xJP455zNwTc5lzcg6PYx5ja2lLz7I9aVu8bZqP5U083Dz4Zt83PHz+kBwZciTvIoGBog555syiaMZXXxn1w2XgvUAuPLxAj7I9KJClACOqjaBDyQ6UyVnGaPeQPi0yAH/M7twR7eeePxediObOFSX0zPQIUSKdQceGixtYfX41mzttxlpjzaaOmyiSrQiOto5Gv19S9lbj40UtBi8vcUrFxdVAn6lbuZB5KmvuHcMpzIkWxepgrbEmQZ+AdYIeD11e2LQEWrUy2rLnh9LqtQz7exgrzqwgVhdL8yLNGVN9DDVda5pkPG+TeB740O1DtC/RPmkvevZMFM2OioJvv4Vy5USP6VatIGNGo4xLVVV239iNp58n+27tI7dDbjqX6oyVxooJdScY5R7Sp0sG4I+Jqopm95cviwxPV1cYM0a8IZUrZ+rRvdfzhOcsC1yG9zFvgiODKeZUjDtP71Aoa6FUL0j/tgSrp09h4UKYMQPu3YPCdfxpPXAHp3SrWP4shPzW+Zlb5ht6b71N/Io5TO86hZCEAEobCvHg18/BRA3igyODcXN0w0pjxfUn1+lSqgujqo96eRTG3FTMXZEMVhnwDfZ9fwAOCRErCosWib+gBg3Ev31FET2mjeTonaMM+msQ5x6cI0/GPEytP5X+FftjpTHfBiNS+iID8McgIUE0RJg+Hc6cEctun38umt3/8IOpR5ckNyNuUnlxZZ7EPqGGSw1mNp5J8yLNsVBMM1u/d0+cypo/X0ywajWKoNWklay8+z9uRopl5p+de/DNulAs905CZ2fPxlINiH/iTGarAtzhw8pEGoNBNbDj6g48/Tw5Hnqc4GHB5M6Ym7+7/22yn2NSWWms6F+xP26Obu9+4tKl4sOlqopuWyNGQNWqRhvHs/hnPEt4Rp6MechsmxmDamB5q+V0Ld3VqN2xJAlkFnT6t2sX9OkjIkaJEjB8uJgF2KVBb9QUuv7kOucfnKdN8TaoqsqIXSPoUKIDNVxrmGxMV66IZeZVq0TeWtPOIWRuPIOtdxcRr4vHoBrQq3o0WDBhr4FvbuaBoUNpHFOMy/H//Tzr7GjH0XF1U3XMcbo41pxbg5e/F5cfX8Y1sysjqo2gX4V+Jk2sMgq9Xmy6588vGiMEBYnkwaFDIV8+o93mnz14mxRuwtp2awFe1gyXpOSSWdAfm0uXxD5usWKv3piWL4eGDU2WXfshToSeYOrRqfgE+eBk70SzIs2w1lgzo/GMVLvn+4poBATA5Mkio9naGtr3u0VslR/Zdvt31FsqnQu3oeFVAwOUbSQoIvPao/dY6PENWFtzZdyON973Q2pFJ1fYszD6b+9PmZxlWNN2DR1KdEi3y6RxujiiE6JxUu1gxQqx9n/9Onz5pSjJWby4KC1mJJceXcLLz4vV51a/7ME7strIl1+XwVdKTTIApxeJdQ2nTxez3o4dYd06KFIE/v7b1KNLklNhpxi1exQHbx8ks01mxn02jq+qfJXqS3tv6wKkqpA50plJk2DvXrAv5sdn4/cwtn1DCuRxpMoSHwYX6sqII3ryTfaB2FgKD2qBb2d3PNw8cHd5lUWe0jKRH+J25G1mHJvB/ef3WdtuLQWyFCDwy0BK5yidrgOGQTXgOt2VDvEFmet9BSIixNn0deugrfGytRNX/RRFYcWZFfxx4Q/Zg1cyCRmA04O1a+GXX8TMN1cu0bXlyy9NPaokSdAn8Cz+Gdnss6EoCjcjbuLd0Ju+FfqS0cY4marv8+9KV6oK4Red6LnCgei7kDO3nka/TGG37jsOYyBg+xT29dzHvYSvcOg2WUyJe/SAESNwL1GCNx3eSk6t6A/1z7PQiqLQrXQ39AY9GgtN+j8Kc+ECFiVKUDFPRXyvBUCdOuIoUfXqRlvVSawX7unnyc91fqZxocaMrTGWr2t8bfpWiNInSQZgc3XvHuTIIY6vXLoENjZiY7JjR/F7M/c07imLTi1ixvEZ1C9Qn5WtV1IhdwVuDbuV5h1gEpeBVb3C80t5iDpeEG14RjTZwunmvYjjGi92RVx7+fwEfQK+wb64f1YPvrMSnXFyvrs6VErLRL7PyjMr6b21NxmtMzK82nCGVR2GS2bzL6DyTgaDWL3x9oZ9+2DbNjzyeTDu+t88XDU/+eeB/yVGG8OKMyvw9vfmRsQNCmYpiM6gAyCbfTaj3EOSkkMmYZmb06fFMvO6dbBhgzhClJAAVlbpYn83NCqUmcdnsvDUQqLio6iXvx5ja4ylQcEGJhtTtQm+XDvsROSdRxiyH0MTWZks+bNhW3kmdw0rqZijHG0ic/FrzC4SFBVrjRX7Pj/42hJzWtPqtay7uI5cDrmoX6A+4THhLA1cSv+K/VPlLHSa0mrF/q63tzgyl0ckstG/PydirlF1SVXWtV9Hx5IdU3wrVVWpuKgigfcDqepclTHVx6RKVydJept3JWGhqmqa/apYsaIqvYFer6qbN6tqrVqqCqqaIYOqfvWVqt68aeqRfbDhO4erFj9ZqJ03dlZPhZ0y6ViiolR16lRVdcymU8nrp/KtrcoPqPxgqeb71ltd7nda3f9NZ9WQ0UFVQfVrUU6duKiH6nf7iOnGHBelTvObprp4u6j8iNp1U1eTjcXo4uPFf7VaVc2XT1XLl1fV1atfPa6qqlavVR0mOqgDtw9M9m2uh19Xx+4Zq8brxHU3B21WDwUfUg0GQ0pGL0nJAgSob4mJcgZsSgaDyGbW60UylU4nyuf17SvapZk5VVU5dPsQnn6ejHQfSd38dbkffZ9YbSz5s5iuf3BEhKjTMGOG+H21lucJrtKV+7oLLwYO3UuM47eOk8Terl4v9hsrvflDalqZdXwW3x/4nqfxT6mdrzZjqo+hSeEmZn+G970uXRKz3b17xTkvGxuxxZIr1xtXddaeX0sxp2KUz13+g25zMvQknn6ebArahEbRcKDXAZMeaZMkkMeQzM/t2yJC/PknnDsHtrYis9nNTVT1N3N6g57Nlzcz9ehUToadJLt9dh4+fwiQKn1Pk9qH98EDsXo/b56oUtiyJVi0/JItdxdhq9piiQWqwYC1HgblKiVetHKlSUtzBj0KIp9jPuyt7LG3sqdhwYaMrj461St/pTpVhf37xZGhnTvFv/HevUVZVBsbyJ37rS/tUrrLB90qPCacduvbvcyuH1N9DEOrDiVPRuP2+ZUkYzP/d/uPhaqCv7+IED4+4pN/+/aiWUKuXK8K9acDHis9OHLnCIWyFmJ+s/n0KtsLO6vUKfzxtiNE8CrxKTQUPD1h/jZ/EvLup3QPC5b1HUml8jYsP1GBShEtGbjyIlcib+BbLjMetXriXq6FuIEJgq+qqhy5c4SpflPZfnU7c5rMYXCVwfSt0Je+Ffqm+XhShb8/1K8vktcmTBDVq5ySlmms1Wv5cfcatgfGExWV940fuhL0CZx7cI5KeSqR1S4rDtYOTGs4jX4V+qVZdr0kpZRcgk4rAQFQubJYWu7fH4YMSRdtAEF00llxZgXDqw3H0sKSNefWYGtpmybJLDUm73/j+VpnRzt+71KXKVNEdUJdXl/o0RCDogXgh1rf8WOdn+HJE/FzLlECRo0SH3pMtMqgqurLYzDHQ4/jZO/EkMpDGFxlcPo/BvPkiSiUodfDd9+JD5ybNkGLFh+ctb/p1G06biuOvd6DbNrBgDjSNaltaeoUd2DhqYXMPD6TZ/HPuDvyLplsMqXGdyRJRiGXoE3hyRPRBjA+Hr7/HipWhDVrRFZzhvRRHvDGkxt4+3uz/MxyYnWxVMxdkTr569CtTLc0G8ObKklpI+w5t7Mghb4DNFrKfzWVs5knkmAQwVdRQTt/KXj8BFmzwvnzomKYibLIDaoBC8UCRVGY5j+NRzGPmNt0Lr3L9cbeyt4kYzKa69fFqs6KFRATIwpmJDZGaJ/Erkb/4r3nBjaGksRZnH/5WLT2MYN3jOTZ3zt5lvCMevnr8XWNr8loLWe7UvqVogCsKMoIoC+gAueBPqqqxhljYOnWlSuiiv/KleINqWXLV29IXbuaenRJEhkXSb9t/fAJ8sHSwpLupbubrJPOPytMacMz8PRYIZ5fzAO2zxgyAMaMsaTOpj/I9CAbT61jMaBibQDrCDf+PHaDlu6FRD9eEwiPCWfeyXksCVzCyX4nyZEhBxs7biRnhpwfxzGY2bNh2DBxRK5bN9EYoXTpFF82LDIWG8vSxFoFoOMxljihV55wT7eRzsU7Mqb6GCrkrmCEb0CSTCvZG2CKojgDQ4FKqqqWAjRAZ2MNLF2aNUvUZ162DDp3FglWW7emi/O7BtXAtXBRjCKTTSZCo0L5uvrXBA8LZmmrpSkKvlsCQ6kxeT/5x+2gxuT9bAkMTfJrxzQqisXTTNz30xN22Y/nVnux7tsO+//l4+epEbi6KnQIaM6jKSFs+SMDDW+VI2/Mjyyv+D+mHLyT7DGnxK2IW3z111e4znDle9/vKZWjFFHxUQDkyZgn/QZfnU6cT7/wIpu8Th0YP14kFS5bZpTgC+JDl61eXOuBzf8AsFYLUslmHWvbrZXBV/popHQJ2hKwUxRFC9gDYSkfUjoSGyuWlStWhPLloV49+PFHGDhQVLFKB+J18fx+/ne8/L249+wed0bcwcHagaOfHzVKXeGkJFG9zaVLsG6qM7eOBkOvjqCJBwWwsOFLKmP4Yy30GcTm3FXQNh3B9uI1ibd8VVc6LRoh/Nv96PsUmVMEBYVuZbox2n00JXOUTPNxGFVUlOhANHMm3LkjOm5Nnw6lSolfRjamUVHG+cTyTFcfKzUvIPaAxzc2XWEUSUoNyQ7AqqqGKoriBdwBYoHdqqru/vfzFEXpD/QHcHV1Te7tzMu9e+Ksy4IF8PixaHpfvjyULCl+pQNP456+TGYJexZGmZxlmNVkFjYakTBjrKL+/67DDBCr1eO568pbA/DFiyJxdv16sV1eefgmTlrGi3GpMPaIgZ93H4HPi0AfcHLKzKbS9f5zndRohPBvqqqy+8ZuToSe4Lva35HLIReLWyymQYEGOGdKmz7Aqernn8VRoqgoqF0b5syBZs1S9Zavynp+kyplPSXJXCQ7C1pRlCzAJqATEAlsADaqqrr6ba/5KLKghw8XwVeng+bNxb6Xh0e6WGaGV/1Nj989TrWl1ahfoD5jqo+hQYEGqdJJJ/+4HbzpX5gC3Jr8+hv5hQvi/X7DllhsqqygQt07bBs2iasx/ngsq4ler8faAPseNsF98CQoWxb47ywbXmXNptabtlav5Y8Lf+Dl78W5B+dwyeTCpcGXcLB2SJX7pakzZ8TPVlHgf/+D4GCzKFQiSelRamVB1wduqar66MVNfIDqwFsDcLqk14s2gI0aiTOj2bKJTkRDh0LhwqYeXZKde3AOTz9PMlpnZF6zeVTNW5WgwUEUcyqWqvdNSpu+ixdh6FR/9of8hWWW+9j/bysxyiMs7IuROWEI7i7u+JbwxPfsVjzajcS9fMvXrpXajRD+7eido3TZ1IWQqBBKZi/JilYr6FK6S6q3VUxVBgNs2yZmu4cPiyYJjRrBr7+mmw+XkpTepCQA3wGqKYpij1iCrgek8+ntP0RFiSb3s2bBzZuv3pC++87UI0syVVXZd2sfnn6e7L6xmwxWGRhSZcjLr6d28IV3t+kLCoKffoJ1fn7Qqw7kT0CnQEmdC7P25qLmscso6hr4+mvcO4zAvcOIt96ndXnnVF2ivB99n/CYcErmKEmhrIUo5lSM+c3mp/9SkQkJ4iD19Olw7Rq4uoqyke4v9ltl8JWkVJOSPeDjiqJsBE4DOiAQWGSsgZnMs2fi3O7SpeL31avDlCkiwSqdmXh4It8e+JZcDrmYWHciAyoNIItdlg+6RlLLQL7Nm2anXYuWZINXTn7ffxZrGz32tXyI0ehAAQsDdPINoWRMRZSNc6B16w8ar7FdeXwFLz8vVp1bRVXnqhzqc4icDjnZ3eM/6Q7pi1Yrjg+BmOXmyQN//AHt2qWLcqiS9DGQlbBAnNMNCwNnZ7G3W7QoVKsm9nsrVzb16JLsWfwzlpxeQnWX6lTNW5WbETc5cOsA3ct0x8byw3sIG3tv9do1+HmCyhq//SifTcWQfzc54yui0JkH1uNR1AQsDRZUuz8QXcG2HB1X94PvYSwBYQH8evhXtl7eio2lDX3K9WGk+0gKZU0/JUPf6MIFMds9eFCkmVtbv7MxgiRJKSMrYb1NQoI41zhjhgjAt2+LN6RLl9JF0/tE957dY9bxWSw4tYDIuEi++ewbquatSoEsBSiQJflFKJKTwfwmt27BkCn+7Hy4ADXvMehxleyaLIy4UYj+f5yiffd+kP1X4pVz2KhluJ2tOIoJjhAZVAMG1YClhSX+If4cun2I72p9x+Aqg43WHN4kVFV0Ipo2TTT9sLMTjRFiYsS/93c0RpAkKfV8mgH48WNxhGjuXLh/XxTP+PFH8UYF6Sr4fr3na2Yen4nOoKNt8baMdh9N1bxVjXLtt52jTer52rt34YdfYlix+xSGro0gVxyKojL2UjZ+8AnHNocd8zz68jiDIzaGvNhQ/OVr0+IIUaJ4XTyrz63Gy9+LkdVG0q9iP/pW6Mvn5T8ng3X6KBv6TocOQcOGYpb7yy+iMUK2bKYelSR98j6tAKzXg0YDZ8+KZKpGjUSiVcOGJm1J9yFUVeVoyFHc87qjsdDgZO9E3/J9Gek+koJZCxr1XknJYH6T+/fh+8nhLDs/F33F2RTs5EGwdQJ6VcVCD5lUK2yXrYJOnchz8RFan/PwhiSt1BYZF8nCAHEW+l70PcrnKv/y7G5qdXdKExER4gOmhQWMHQu1aolD1S1bpqsPl5L0sfv4A7DBIDKYZ8wQVXu8vaFuXVGzuUgRU48uyXQGHT5BPnj5eXEy7CSbOm6ibfG2fF3j61S757symN8kPBzGewaz5KI3+jJLoVYMTaLz0f72BYaUtCZBn4C1tRUeXpvAtTqQ9keI/qnNujb4BvvSoEADVrVZRb389VLlLHSauXFD/DtftkwsL3foIB5XlFe/lyTJbHy8ATg6GlatEuXzrl4VCVYtXvSAVZR0E3wT9AksDFjI9GPTuRV5i8JZC7Og2QKaFGqS6vdOanB8+hRGePuz5qgvCcWXo6lwi+4Pc/PNphhKPL0HPXpQvPN8fMP88XDzwN3F/T/3SYuAe/HhRWYcm8Hk+pPJZp+NiXUnYmtpS/nc5VP93qlu5kxRFMbSUjT9GDkSypQx9agkSXqHjzcL+ssvYdEikcU8YoRojZZ47CIdSNAnYK2xRm/QU3ROUXJkyMGY6mNoWbSl2RTzf/5cZdjMfay4NQF9zpNgmYC1YsH6tVpaPc4GgwbB4MGiKbuJqKrK4TuHmXp0Kjuu7cDeyh6fjj40KtTIZGMyCr0eNm8WqzrFionqVevXiz7TefKYenSSJL3wrizojyMAqyocOyaW38aOhQoVxKz30SNxjjcdLStefnwZb39vdl7fyZUhV7C3sudxzGOzatgeE6dj8LyNrL45FV32QGy01mittBhQ0SgaJji04JsBa8DetL1uY7Wx1F1Vl2N3j+Fk78TQKkMZVHkQ2ezTcQJSdLTIW5g+XaSXjxwpspslSTJLH+8xJK0WNm4UgffECXB0FIUEKlQQS8zpZJlZVVWO3DmCl78Xf175E1tLW3qV7UWsNhZ7K3uTB1//EH98g32p6erBuX3FGHapIrqMt8hjlYUfdtpQ9F48TXprSNCAtcYajw5fmyz4xuniOHLnCPUL1MfOyo7yucrTs0xPepXrhb2VaT8QpNivv4KXF0RGikpVXl7QqpWpRyVJUjKl3wBsMIgORBcvikA7dy706iXa56Qzp++dptaKWmSzy8b3tb43q3On/iH+1F1Vl3hdAuhsUFfso27VXAy5eJuW156i6dAR3696UOj2DUJiAnCxr8SDx67gkrbjjIiNYH7AfGYen8njmMcEDwvGJbML85rNS9uBGFtQkFhiVhSR5Va/PowaJQrFSJKUrqXvJeglS8R+V+PG6eYYEUCMNoYVZ1YQERvB+FrjUVWV9RfX06JoC7Oapd18covmy7sS9OyYaF9k0NAt9wR+y10S5dBBGDaMLeGaNO9E9E8Pnz9kypEpLDq9iOiEaBoVbMTYGmPxcPNIvxnNqioagEybBnv2wL59InNfVdPVdookSZ/CHnA68fD5Q+acmMO8k/MIjw2njlsd9vbca3bF/G9G3OTLDf9jb9gGFANYKAYArC2s2Pf5wdeymGtM3v/Gs8LOjnapWkoyXhePjaUNoVGhFJlT5GURkrK5yqbaPVOdVgurV4ujchcuiApVQ4eKhMIsH1bDW5Ik8/Dx7gGnI6vPrabvn31J0CfQsmhLxlQfQ3WX6mYxS/MP8edA8AEq5amES0JDRvwER9y28tVpa8Ydi+NW1WIcal0Wj2aD/3OEKKXVsj6Eqqocun2IqX5TidfFs7fnXpwzORM6MhRHW0ej3y/NJBaIMRjgm28gRw5YsQK6dBGlIiVJ+ijJAJxKEitWZbHNQskcJamUpxK9yvZipPtIijqlfpWnpJq414dvj3ZCVXWgWqCsOIxDRHWOu9SnuJsFVttHkadmTWq85YNCcqtlfQi9Qc/WK1uZenQqx0OPk90+O0OrDsWgGrBQLNJv8L1xQ2QzHzggqrPZ2MDx46IloBl8MJMkKXWZ19rnR0Bv0LPx0kbcl7pTc3lNpvpNBUTv3YUtFppN8H2e8Jy+m37mu8O9UHnRChADOasuZ86We5Q5uxmrHVtFGcN3BIMxjYpiZ/X6uWRjl5JcfHox7da343HMY+Y3m8/t4bf5tta3Zrd0n2R+fiJbv3BhcVa9ShXR+hIgXz4ZfCXpEyFnwEa04swKJhyawM2ImxTMUpC5TefSq2wvUw/rjbwPzWTphR8o/NCG206gV8ACDfYuriw8EUTPuknrkJMapSQjYiNYELCAwtkK075Ee7qV7kY2u2y0Ld7WbIqQJNvBg+DhIfZ0x42ThTMk6RMmA3AKPXz+ECd7JywUC66FXyNHhhx4NvCkVdFWZhMs/EP82RS0iVsRt+hSshsR/m1ZN6Et+/mBnA/yMaFBJfaWUbCkHHqr4h+8f2usUpJ3o+4y3X/6y4zmwZUH075EezLaZKRDyXRay/j5c1E4w2AQCVU1a8LSpdCxIzg4mHp0kiSZkMyCTqbEilWrzq5iQ4cNtCjaAq1ei5XGvMpdLg9cTr9t/dCrelCh5RlX/tx6G3d3yOW8jtMFMvxnyTO1M5jfZMLBCUw4NAGDaqBzqc6MqT4mfWc037sHs2eLrkQREdC0KezYYepRSZKUxmQWtJEk1hX28vNi29Vt2Fra0rtcb0pkLwFgdsG3/7b+LD69GFRAAY0BMkVnZev6eFq0t2Hrmc8I+oBuR8aUmKRWOkdpMttmpphTMQZUGsBI95G4Obql+v1T1YIFYrar00Hr1jB6tCiJKkmS9A/pNIvFNPSqnp6be+J/158fav/AneF3WNB8gdH78CaXVq9l7fm1PE94DkCFi7YMPg62OrDQK2g0tvRfMo+WHWxQFLF0PKltaZwd7VAQM9/ULqBhUA1subyFGstqUHN5TZacXgJAh5IdmNVkVvoMvqoKe/eKrGYQpVD79RP1yH18ZPCVJOmN5BL0O0QnRLM8cDlrL6xlf6/92Fracv7BeQpmLWg2Fav8Q/zZfWM3kXGRbA7y4XbUHWY7DuJC0Fz+WPyMvlYriRtbnOy1T9Cw8H9bAaYVVVVZFrgMTz9ProRfwc3RjdHuo+lTvo/Z/Cw/mFYL69aJmsxnz8KwYaIuuSRJ0gtyCfoD3Xt2j9knZjM/YD6RcZFUd6nO/ej7uDm6UTpnaaPdZ0tgaIqyhw/fOUzdlXXRGXQAlH2kYc4eeHIzmqXAwMEZGfvdELJnB6hntHF/iMS2ioqi8MfFP7CzsmNtu7W0L9EeS4t0/M9v9myYMgVCQ6FECZFY1bWrqUclSVI6ko7fAVPHlcdXKLOgDFq9ljbF2zDKfRTVXYy/hLglMPS1GsqhkbF843Me4L1B+EnsE7LaZeXI7SMvg6+FAcpeLMavVxeTt707QZOgUKEXQX6p8Y4IJdW9Z/eYeXwmSwOXcrr/aVwyu7C+/XocbR3NovpXsoSFifKQiiKaJBQtCosXQ6NG6aoWuSRJ5uGTD8CqqnIg+AA3I27St0JfimQrwg+1f6BjyY4Uyloo1e7ruevKa8lPALFaPZ67rrw1QJ4MPYmnnyc7Lv/Jzc/PUjufBzaqNVpVi6q3IdB+MYv83V82yklJkE+uK4+v4OXnxapzq9AZdLQv0f7lh4Qsdum0nvGpU6Ixwvr14OsLn30Gs2aB5Sf/v48kSSnwyb6DaPVaNlzagJefF4H3AymQpQB9yvVBY6HhfzX/l+r3T0oN5S2BoXz3lw+34raiWN4mmutkTrBg6HEDD59u4NvT3xJ/zRenSr6Mau/B2Anur50oSk6QT4nHMY8pPb80GgsNX5T/glHuo8wmQe2DGQzw118i8Pr6QsaMMGIEFCggvi6DryRJKfRJvovsur6Lftv6ERIVQjGnYixusZjuZbqnaeGM99VQ3hIYynCf9dzWjAWNFlQYchxG3SrC307/o+K0TmR2gjk/utO/vztWbzgBldqNElRVZef1nRy+fZhJ9SfhZO/EmrZrqO1W22z6GSdbfDz06QN2diII9+0LmTKZelSSJH1EPpkAfDfqLjqDDjdHN5wzOVMwa0HmNZtH08JNTVJTeEyjom/sozukXl5mHJvBjL/+5qmaE9CLOs0G2GvfnnnX12IVbMnIsaJxTubMb79HajVK0Oq1rLu4jqlHp3L+4XlcMrkw9rOxONo6pt+KVY8fw/z5sGuXKBdpZwf790OxYrzx040kSVIKffQB+Mz9M0zzn8YfF/6gQ4kO/N7ud0rlKMWBXgdMOq7EJeDvd/oQEhNALrsilMr/gMG7VxChe0atEAjPOJanWS1RDXoMemsunxhJhqIPuLDDGTe399/jbUE+JYU2AsICaLe+HXee3qFk9pKsbL2SzqU6Y61Jp23zrl0THYlWrIDYWGjSBJ48gezZobTxMt4lSZL+7aMNwPtu7mPSkUnsu7UPB2sHhlQewrBqw0w9rNfkdLrDdXUs8VbxROoMXL4KbYNgzIVMXM7ZkskZ62Gxvz76LCexiilPtrp6ChS/gptb0vZvjdUo4XHMY8KehVEmZxkKZS1EcafizG0612SrB0bj5ycSqqysoHt3GDkSSpY09agkSfpEfFQBOF4X//LM6d/X/ybocRBT6k+hf8X+Ztcz9kToCVadW0WCPgEDBhQVRgQ5Ms1jIpfG9GLWMAuuHLDFMstzshdxxK7wA+ytoxnT6MNmZSlplBAcGcw0v2ksDVxKMadinOp/CkdbR/7u/neyrmdyej1s3QpRUdC7N1StCpMmQa9ekCuXqUcnSdIn5qOohPUk9gkLAhYw+8RslrdaTuNCjYmKj8LW0taslkYTk5am7p/AwfvHqPo8C+cyx4liFYolG1vuY9uCGixaJPJ92nweSVCWQO5Hx6TpGd6LDy8y8chE1l1Yh4ViQfcy3RldffTLmtfpzvPnYol5+nRRLrJqVfD3l313JUlKdR9tJaybETeZcWwGSwOXEqONoWHBhmS1ywpAJhvzyFj1D/HHN9gXC8WC1ScWc+HZDfI+Be/jCn3d6nFhwFfsvXOUB8c86FLLnefPYfBg+P57cHJyBOqkyThVVcWgGtBYaDgeepw/r/zJ8GrDGV5tOHkz5U2TMaSKtWtFz90nT0TgnTwZ2rSRwVeSJJNLtwHYoBqou7IuYc/C6Fq6KyPdR1ImZxlTD+s1e2/upeXaliToE7BQweWJnlXHbOn82ZdYrR2J6uJKmA8sH1OLW7egeXPw9BSJt2lFb9Cz5fIWpvpNpUupLgyvNpzuZbrTplib9Fs4IyhInNvNmxdcXKBWrVcdiWTglSTJTKTbAGyhWLCqzSoKZS1Enox5TD2c19yPvs/so954H5tBPDpUVFA0fJGrCT12r4XMmQkMhOE94NAhkWy7Zw/Ur592Y4zXxbPq7Cq8/L24Gn6VAlkKkMtB7INaa6yxtjOfpfskUVXxw/Tygu3bxax39myRZPXZZ6YenSRJ0n+kKAAriuIILAFKIbrOfq6qqr8RxpUktfLVSvZrU9oI4U2uhl9l2v5fWHnxdxLQ43EL/Nws0GkssNZYU6f7dzyIy8z4UbBsGWTLJlrH9u0LmrSrAQJAV5+u+AT5UCF3Bda1X0e74u3StBCJUW3aJJaWAwLE8aGffoKBA009KkmSpHdK6Qx4JvC3qqrtFUWxBtJFXzlj1EjeEhj68gyvi30lfm7Slinb6hOoDaHPGRhpV5fCg7/HP78VvrcPUt3Zg0O/u9PoV4iLEydevv0WHB1T6Zv8l3vP7jHr+CyGVRtGLodcjKk+hoGVBlIvf7302RwhJgbsX/xz+/tvkdm8cCH06CGKaEiSJJm5ZGdBK4qSCTgLFFCTeBFz6QdcY/L+N1aIcna04+i4uu99/ZbAUIb5/MEdzThAh4IVrvopTLJ9RN2wUHIOGw9FigBiZXTzZhgzBm7ehBYtRGXDwoWN/V292bXwa3j6ebLy7Ep0Bh2rWq+iW5luaXPz1BAWJhohLFggqlZVrSqCr4OD7EgkSZLZSa0s6ALAI2C5oihlgVPAMFVVn//r5v2B/gCurq4puJ3xpKRGcoI+gbHbpvPYMB8sRZcfVC1P1TPM0fSiy9xXAfzcORg+HA4cEPUddu+GBg2M8R28n96gp5tPN9ZfXI+1xpo+5fowuvroVO3wlKouXBCfXNasEed527V7VZtZ1miWJCkdSkkAtgQqAF+pqnpcUZSZwDjgu38+SVXVRcAiEDPgFNzvpZTu3ya3RvKz2KeUnJqPEJ5SIBpCMoPOQgGssDWUfhnAw8PFMaIFC8QS85w58OWXqd9AR1VVztw/Q/nc5dFYaMhkk4lxn41jWNVh5HTImbo3T01xcSKTOT4eBgwQn2oSuxJJkiSlUykJCXeBu6qqHn/x542IAJyqjLF/m5QayYl7vLdjDuOkUZjW2pvW5Z3pfTsL2a47cTpnd3bbZyBOcwFbQ2lsDMXJndGeOXNE8I2KgkGDRD5Q1qxG/AG8gc6gY9OlTUw5OoXA+4FcHHSREtlLsKjFotS9cWrRamHDBtiyBf74A2xtwccHypRJ/R+mJElSGkl2AFZV9b6iKCGKohRVVfUKUA+4ZLyhvZkxety+r0bylsBQBvnM4b5mKqqlgSjg29V5gZH8PPcSWy4/YYHPeWy0emx0ojqU4W52Qo+X56vrULcuzJwJpUoZ7dt+ozhdHMsDl+Pl78XNiJsUzVaUJS2WUDBLOu3BGxUFS5bAjBkQEiIORIeFifO8Hh6mHp0kSZJRpXRR9CtgzYsM6JtAn5QP6d2M1eP2bTWSbz6+xtjfW3Ivw2XxgAKKAZ5o7r4I8nVfC+C3gyH2SCmeXMyBm5s4EZPahZZUVUVRFNYHXGPIrhFYGvJR1PonJlbvQ9sKLql349R09qxYZo6KEsF2/nzRmUgmVkmS9JFKUQBWVfUM8MbsrtSSGj1uDaqB8JhwsmfIjt2jSJ6pl+lw0ZbNxbXoLFRQLFEtK78W5BsWdeb0ZmemrhBneH/5BUaNEqulqSXsWRjT/adz9sFZBpZcxqQdoeTSzcFSzUVcvML4zRexUCzSpF60UZw7B7dvi9TwkiVFR6I+faBSmv6TkiRJMol0N70Y06godlavF4xIbo/bg8EHab+yGQV+zEL7n8V6ce7ilWkaM4OA/Otw0k3GUdednAm/YmMoTh5HO1QVNm6E4sVhwgRo2xauXIHx41Mv+F55fIW+f/Yl/8z8eB/zxsneiSm7zhGr1WOl5kZBTLcTl+LNmqqKsl+NGkHZsuJAtKqKDLW5c2XwlSTpk5HuSlEao8dtVHwU/1v/JXNv/CGWmBXo/cwZNSEBxdqa5v3ac9TnPAZtcWwMxQER5DsVKkn9+rB/v8gH+u03sWr6NsaotvX39b9puqYpNpY29C3fl1HVR1EgSwHyj9vxxud/6FJ8mtq/XwTcs2dF+7+JE0V6eHosBCJJkpRC6S4AQ8p63AKsWDKEuY//ePlnCwsNNl16oFhbv7w+vAryOWwdyHmtIiMnOZApU9KOFSU3W1tVVfbc3EO8Lp4WRVtQO19tfvL4iS8rfUmODDlePi81luJTRVSUyGrOlk38WasVdTi7dgUbG9OOTZIkyYTS3RJ0clwOO0dfr1qsXv01AJ+3+4WlmnbYWdmhUTRYa6zxcPN47TWtyztz+Ou6fF+0GcHza7N9rQNffAFXr4p2ge870/uubO030Rl0rLuwjoqLKtJodSOm+k0FwM7Kju9qf/da8AXjLsWnirt3RfkvFxexQQ5Qpw6cPy/2eWXwlSTpE5cuZ8BJteigNzN8JxPEI2y1UPh2NHQHh1yufP7tRoq/6NXr4eaBu4v7a689dUo01Dl2DNzd4a+/oGLFpN/7Q7K1t17eyqjdo7gRcYOi2YqytOVSupV+d7lIYyzFp4pz50RHorVrxd5ux46iPjOIpWa53CxJkgR8xAG45YQSbDMEgQqWqsKmMhNo2uF/rz3H3cX9P4H3yRORULVwIeTIAStWiPjxoadh3rdEHBkXiYJCZtvMGFQD2eyz4dnAk1bFWmGhJO1mKV2KNxpVfRVYp00TRTOGDIFhw8DNzaRDkyRJMlcfzRJ0vDaO5eu+IeJxCABZ7LKgqIACqsaCs7kt3jn7MhhEDYgiRWDxYhg6VGQ39+qVvKOob1si/qK2I2P3jMV1uive/t4AtC7WmmNfHKNN8TZJDr5mQauF1auhfHk4fVo8NmmSKKIxfboMvpIkSe+QbmfA/i+WjyvlLE/g3t+YGbKBMFst2ru36T/qdwZ08mLDqnok6BPeuMf7TwEBYl/3xAnRu33uXJHlnBL/XiLOkikcpzy7GLhvIzqDjo4lO9KmeBuA9NcOMCpKfEqZMUPs9RYvLh4DyJPHpEOTJElKL5LdjjA5jNWO0D/En3qr6hGni0V9McttcM+er4v3pV7fiSgZMrx83tv2eEE0Tfjf/0QsyZkTPD2hW7fU2aZss64NO6/tfNmVqGDWdFouUq8XjRDu3BEVq8aMgcaNZcUqSZKkN0itdoQm4xvsS4I+gRexl8FOzZj93Z+vBQFxBjeWsMgybHeMZUyj0JezUoNBnIQZOxaePhXNdX780Xhd7VRVZd+tfUw9OpW5TedSOFthpjWcxvxm88nlkMs4N0lL586Jpgi//irKfk2aJNbqZdEMSZKkZEuXAdjDzQNrjbVYXra0pmvL8f8Jvm87g5sPZwYNEtnNNWuK5ebSpY0zLr1Bz+bLm5l8ZDKn7p0il0MubkbcpHC2whTIks7a56kq7NsnMpp37QJ7e7EhXrSoOMMrSZIkpUi6DMDuLu7s67nvrcvLbzqD+zxaYeBgAw+Pg5MTrFolSg8ba7lZb9BTfmF5zj88T+GshVnUfBE9yvbA1jIVi0Onllu3RI3NM2fE2vyvv4o+vLIVoCRJktGkywAMbz5ClOifZ21VFZ5fcibiQDEMMTYMHiTqQjg6pnwMUfFR/HnlT7qX6Y7GQkPPsj1xc3SjTbE2aCw077+AOXn2TFQZqVgRnJ1F5aolS8SmeGp2mJAkSfpEpdsA/C6JZ3ATHjnwZE8p4kOyYZ07gmJ9zjNnduUUX//h84fMPDaTuSfn8jT+KWVzlqV0ztKMrj7aCKNPY2FhMGsWLFgADg4QHAzW1rB3r6lHJkmS9FH7KAPwV7WKMXhMHE+OuWFhrSNro3NkrxTGT+1SttkbHhPOD74/sDRwKfG6eNoWb8vYGmMpndNIm8hp6do1kUy1erXIbG7bVmQ0v6/GpiRJkmQUH9277bZtMO6rPDy5DTkqhmFT/SIueTSMaVQ62VWjohOicbB2wNbSFp8gH7qW6srXNb6mqJOZ1F1OKlUVxTOsreHGDZHZ3L8/jBgBBdPpsShJkqR06qMJwHfuiMqHW7aI3u6HDkHNmnmA5BeGOHLnCJOPTOZGxA0uDLxABusM3Bh6AzsrM+s49D46HWzaJDKa69SBqVNFP96QkFddiiRJkqQ0le6rJ2i1Iq4ULy5Oy0yeLKoi1qyZvOsZVAPbr27ns2WfUXN5TY6HHqdb6W5oDVqA9BV8nz+H2bPFmd3OnUW1qsQzV4oig68kSZIJpesZ8NGj4nTMhQvQooXIJUpp+eHtV7fT6o9WuGZ2ZXaT2Xxe/nPsreyNMt40N2wYLF0K1auDtze0bCkrVkmSJJmJdFmKEmDiRNG1yMVFTPJatUredWK0MSwPXI6Vxor+FfujM+jwCfKhTbE2WGmsjDLWNHP1quhGNGSImOlevizaO1WvbuqRSZIkfZLeVYoy3U6HEssQX7qUvOAbERvBr4d+xW2GG0N2DuGva38BYGlhSceSHdNX8D16FFq3hmLFYOXKV52JihWTwVeSJMlMpdsl6OrVkx9bVpxZwVc7vyI6IZomhZrwzWff8JnrZ8YdYFpQVZFMtWePqFL17bdi9psjh6lHJkmSJL1Hug3AH+pa+DUcrB3InTE3BbIUoEWRFoytMZayucqaemgfJjYWtm6FTp1EIlW9emJvt08feNEFSpIkSTJ/6XYPOKlO3zvN5COT2XhpI0OqDGFWk1lpen+jefwY5s2DOXPg0SPw8wP3N5filCRJkszDR9eOMCkOBh9k4pGJ7L6xm0w2mRhbYyzDqg0z9bA+XEQEfPed6J8YGwtNm4rN72rVTD0ySZIkKQU+qgCsqirKi/ZGK86u4Mz9M0yqN4mBlQaS2TaziUf3gZ48Efu6dnaivFenTjB6tKgyIkmSJKV7H8UStFav5ffzvzPl6BRWtVlFpTyVeBzzmAxWGdJX4QyDAXbuBE9PuH1b1Gu2tIT4eLCxMfXoJEmSpA/0UR5DAnGGd9bxWRScVZDe/2/v7mOkqs44jn8fVkiDpahdSiliwWo0bUTZ4AvyUmJrI5RIi5FI0NJAQnzHRFNAE1+jRptWbYOibbG0xVZrixJ8iUSrjYkQkBfBYBWqxq0UKK1QbFVkn/5xLnUy3LsMu8Oce8bfJ9nMcM+Z3efhzL3P3jN373n8+/Rs6cl/94SlCFt7t6ZTfD/8EB58MPzt7oQJsGkTXHZZuIUkqPiKiDShZKeg93bsZeh9Q9n8r82MOmYU8yfMZ9xx4/4/BZ2UZctg+vRQgBcuDLeN7NUrdlQiInIIJVuAW3q0cMPXb2DIkUPS+xve9na4+25obYU5c8KFVc8+GxZKSPEXCBEROWjJFmCAi06+KHYIB2f9+rByxEMPhZtozJwZtvfoAWedFTc2ERFpqKQ/A07KbbfB0KHw6KNw6aXhc957740dlYiIRJL0GXCp7VuDt60Njj8ezj47XOV8ySVaBlBERLp/BmxmLWa2xsyW1iOg5FWvwbtgQdh+6qnhXs0qviIiQn2moGcBG+vwfdJ3++1wzDFw5ZUwYAAsXgy33ho7KhERKaFuFWAzOxr4NvDz+oSToLffDhdUQbi6efRoePHFT5YI7KGP2UVEZH/drQ53Az8AOoo6mNlMM1tlZqu2b9/ezR9XIitWwHnnwZAhodhCmHp+7DEYOTJqaCIiUn5dLsBmNgHY5u4vd9bP3R9w9+HuPrxfv35d/XHl0NEBS5fCmDFhMYTnnoO5c8NFVqCzXRERqVl3roIeCZxrZuOBzwCfM7PfuPuF9QmthD74IKy727s33HUXzJgBffrEjkpERBLU5QLs7nOBuQBmNha4pumK786dcP/98MQT4Wy3d+/weOKJ0LNn7OhERCRhmjPN094e1twdNAhmzw7FdseO0HbSSSq+IiLSbXW5EYe7Pw88X4/vFd3KlXDmmeHK5smTwxq8bW2xoxIRkSajO2G5wwsvwNatYdH7tja49trwWe/gwbGjExGRJvXpnYLeuzfcl/n008MqRLfcEopxSwvcdJOKr4iIHFKfzgL89NNwwglw/vnw3nswf36YetZSgCIi0iCfninoHTvCAgn9+8Phh4e1eO+8EyZODGe9IiIiDdT8Z8BvvglXXBHu0XzzzWHbqFHw0kswaZKKr4iIRNG8Z8Br1oQz3EceCUV26tSwDi9oqllERKJrrgLs/klxnTcv3EDj6qth1iwYODBubCIiIhWaYwp6zx5YtAiGDQuLJEBYBvCdd8JZsIqviIiUTNoFePduuOceOO44uPBC+OgjeP/90Na/P/TtGzc+ERGRAulOQXd0wCmnwObNYQ3eefNg/HitSCQiIklItwD36BFunjF4MIwYETsaERGRg5JuAQaYMiV2BCIiIl2i+VoREZEIVIBFREQiUAEWERGJQAVYREQkAhVgERGRCFSARUREIlABFhERiUAFWEREJAIVYBERkQhUgEVERCJQARYREYlABVhERCQCFWAREZEIzN0b98PMtgNv1/FbtgL/qOP3i0m5lE+z5AHKpayaJZdmyQPqn8uX3b1fXkNDC3C9mdkqdx8eO456UC7l0yx5gHIpq2bJpVnygMbmoiloERGRCFSARUREIki9AD8QO4A6Ui7l0yx5gHIpq2bJpVnygAbmkvRnwCIiIqlK/QxYREQkSSrAIiIiESRRgM3sHDP7i5ltMrM5Oe1mZj/J2l8xs7YYcR6ImQ0ysz+Z2UYze9XMZuX0GWtmO81sbfZ1fYxYa2Fmb5nZ+izOVTntpR8XMzuh4v96rZntMrOrqvqUdkzMbIGZbTOzDRXbjjKzZWb2RvZ4ZMFrO92vGq0glx+a2WvZ+2exmR1R8NpO34uNVpDLjWb2t4r30fiC15ZmXAryeLgih7fMbG3Ba8s2JrnH36j7i7uX+gtoATYDxwK9gHXAV6v6jAeeAgw4A1gRO+6CXAYAbdnzPsDrObmMBZbGjrXGfN4CWjtpT2JcKuJtAf5O+MP5JMYEGAO0ARsqtt0JzMmezwHuKMi10/2qJLl8Czgse35HXi5ZW6fvxZLkciNwzQFeV6pxycujqv1HwPWJjEnu8Tfm/pLCGfBpwCZ3/6u7fwT8DphY1Wci8CsPlgNHmNmARgd6IO6+xd1XZ8//DWwEBsaN6pBKYlwqfAPY7O71vFvbIeXufwb+WbV5IrAwe74Q+E7OS2vZrxoqLxd3f8bdP87+uRw4uuGBdUHBuNSiVOPSWR5mZsBk4LcNDaqLOjn+RttfUijAA4F3Kv7dzv5Fq5Y+pWJmg4FhwIqc5hFmts7MnjKzrzU2soPiwDNm9rKZzcxpT21cLqD4YJLKmAD0d/ctEA46wBdy+qQ2NgDTCTMqeQ70XiyLy7Pp9AUFU50pjctoYKu7v1HQXtoxqTr+RttfUijAlrOt+m+naulTGmb2WeAPwFXuvquqeTVhCvRk4KfAYw0O72CMdPc2YBxwmZmNqWpPZlzMrBdwLvD7nOaUxqRWyYwNgJldB3wMLCrocqD3YhncB3wFOAXYQpi+rZbSuEyh87PfUo7JAY6/hS/L2dbtcUmhALcDgyr+fTTwbhf6lIKZ9SQM/iJ3/2N1u7vvcvfd2fMngZ5m1trgMGvi7u9mj9uAxYRpmkrJjAvhILHa3bdWN6Q0Jpmt+6b6s8dtOX2SGRszmwZMAKZ69oFctRrei9G5+1Z33+vuHcDPyI8xiXExs8OAScDDRX3KOCYFx99o+0sKBXglcLyZDcnOUi4AllT1WQJ8L7vq9gxg574phTLJPjP5BbDR3X9c0OeLWT/M7DTCGO1oXJS1MbPDzazPvueEi2U2VHVLYlwyhb/NpzImFZYA07Ln04DHc/rUsl9FZ2bnALOBc939PwV9ankvRld1/cN3yY8xiXEBvgm85u7teY1lHJNOjr/x9pfYV6bV8kW4mvZ1wlVo12XbLgYuzp4bMC9rXw8Mjx1zQR6jCNMWrwBrs6/xVblcDrxKuMpuOXBm7LgLcjk2i3FdFm/K49KbUFD7VmxLYkwIvzRsAfYQfkufAXweeBZ4I3s8Kuv7JeDJitfut1+VMJdNhM/e9u0v86tzKXovljCXX2f7wSuEg/eAso9LXh7Z9l/u2z8q+pZ9TIqOv9H2F92KUkREJIIUpqBFRESajgqwiIhIBCrAIiIiEagAi4iIRKACLCIiEoEKsIiISAQqwCIiIhH8D6VtpUle1YRmAAAAAElFTkSuQmCC\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:                Sun, 16 Aug 2020   Prob (F-statistic):           1.66e-29\n",
      "Time:                        18:00:44   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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
