{
 "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:                Fri, 10 Jul 2020   Prob (F-statistic):           5.44e-29\n",
      "Time:                        05:46:36   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3xN9//A8dfJzUZEbLH3FoQKVSE1qvYsarWq1rd2tV/9tdX2a8WqtpSq1mpra62iiJUYsdVWQRIkQhCZ997z++MjUa0Ryc29N7yfj0cejXPPPeeTS/M+n/V+a7quI4QQQgjrcrB1A4QQQogXkQRgIYQQwgYkAAshhBA2IAFYCCGEsAEJwEIIIYQNSAAWQgghbOCpAVjTtPmapkVpmnbib8d8NE3bq2naEU3TQjVNq5u1zRRCCCGeL+npAf8ItPjHscnAOF3XfYCP7/9ZCCGEEOnk+LQTdF3fqWlayX8eBjzuf58biEzPzfLly6eXLPnPSwkhhBDPp4MHD97QdT3/o157agB+jGHAJk3TpqB60fUfd6Kmaf2B/gDFixcnNDQ0g7cUQgghshdN0y497rWMLsIaCAzXdb0YMBz4/nEn6ro+V9d1X13XffPnf+RDgBBCCPHCyWgA7g2suv/9ckAWYQkhhBDPIKMBOBJodP/7JsA5yzRHCCGEeDE8dQ5Y07SfAX8gn6Zp4cAnwDvAl5qmOQKJ3J/jzYiUlBTCw8NJTEzM6CWyDVdXV4oWLYqTk5OtmyKEEMLG0rMKuttjXqptiQaEh4eTK1cuSpYsiaZplrikXdJ1nZiYGMLDwylVqpStmyOEEMLGbJ4JKzExkbx58z7XwRdA0zTy5s37QvT0hRBCPJ3NAzDw3AffVC/KzymEEOLp7CIACyGEEC+abBeA1xyOoMHEbZT6YD0NJm5jzeGITF/TYDDg4+NDlSpVqFGjBtOmTcNsNj/xPWFhYfz000+ZvrcQQogXU0YzYdnEmsMRfLjqOAkpJgAiYhP4cNVxANrV9M7wdd3c3Dhy5AgAUVFRdO/endu3bzNu3LjHvic1AHfv3j3D9xVCCGEdaw5HELjpDJGxCRTxdGN08wqZihuWkK16wIGbzqQF31QJKSYCN52x2D0KFCjA3Llz+frrr9F1nbCwMBo2bEitWrWoVasWwcHBAHzwwQfs2rULHx8fpk+f/tjzhBBC2FZq5y0iNgGdB503S4ygZka26gFHxiY80/GMKl26NGazmaioKAoUKMCWLVtwdXXl3LlzdOvWjdDQUCZOnMiUKVNYt24dAPHx8Y88TwghhG09qfNmy15wtgrARTzdiHhEsC3i6Wbxe+m6DqhEIUOGDOHIkSMYDAbOnj37yPPTe54QQgjrslbn7VllqyHo0c0r4OZkeOiYm5OB0c0rWPQ+f/31FwaDgQIFCjB9+nQKFizI0aNHCQ0NJTk5+ZHvSe95QgghrOtxnbSs6Lw9i2wVgNvV9GZCh2p4e7qhAd6ebkzoUM2iQwjR0dEMGDCAIUOGoGkat2/fpnDhwjg4OLBo0SJMJjWMkStXLu7evZv2vsedJ4QQwras1Xl7VtlqCBpUELb0mH1CQgI+Pj6kpKTg6OhIz549GTFiBACDBg2iY8eOLF++nMaNG5MjRw4AqlevjqOjIzVq1KBPnz6PPU8IIYRtpcYMe1sFraXOdVqDr6+v/s+FSadOnaJSpUpWa4OtvWg/rxBCvMg0TTuo67rvo17LVkPQQgghxPNCArAQQghhAxKAhRBCCBuQACyEEELYgARgIYQQAjgXc86q95MALIQQ4oV25NoRGv3YiIrfVOT8zfNWu2+22wdsaTExMQQEBABw7do1DAYD+fPnB2D//v04OzvbsnlCCCGygNFs5FbCLfLnyI+boxuXYi8xrdk0CucsbLU2vPABOG/evGmlCD/99FNy5szJqFGj0l43Go04Or7wH5MQQjwXElIS+OHIDwQGB1KzUE1WdV1FhXwVuPDeBQwOhqdfwILsKrIMGwb3Y6HF+PjAjBnP9p4+ffrg5eXF4cOHqVWrFrly5XooMFetWpV169ZRsmRJFi9ezMyZM0lOTuall15i1qxZGAzW/UsUQgjxZLcSbjE7dDYz9s4gOj6aekXr0cenT9rr1g6+IHPAj3X27Fn++OMPpk6d+thzTp06xdKlS9mzZ09aFaQlS5ZYsZVCCCHSY8beGYzdNpbaRWqzo88Ogt8Kpk2FNjZtk131gJ+1p5qVOnfu/NSe7NatWzl48CB16tQBVE7pAgUKWKN5QgghnuDMjTMEBgfSrmI7WpVvxX9e+g8dKnWgRqEatm5aGrsKwPbk78UUHB0dMZvNaX9OTEwEVM3g3r17M2HCBKu3TwghxL/tj9jPpD2TWH1qNS6OLlQrUA2AfO75yOeez8ate5gMQadDyZIlOXToEACHDh3i4sWLAAQEBLBixQqioqIAuHnzJpcuXbJZO4UQ4kXW99e+vDTvJbZd3MZ/G/6XS8MuMbTeUFs367GkB5wOHTt2ZOHChfj4+FCnTh3Kly8PQOXKlfniiy9o1qwZZrMZJycnvvnmG0qUKGHjFgshxPPPZDax+vRqWpVvhaujKwGlAqiavyr9a/cnl0suWzfvqaQcoZW9aD+vEEJYWqIxkQVHFhAYHMiFWxdY2G4hPWv0tHWzHulJ5QilByyEECJbMJqNTAmewoy9M7h+7zp1vesS2DSQthXb2rppGSIBWAghhF2LT4nH3ckdg2Zg5amV+BTyYUyDMfiX9EfTNFs3L8MkAAshhLBL52LOMSV4CitOreDskLPkdc9LUO8gcjjnePqbswEJwEIIIexKaGQok/ZMYuXJlTgbnOnr0xej2Qjw3ARfkAAshBDCjly8dZE639XBw8WDD17+gPdeeo9COQvZullZQgKwEEIImzGZTaw6tYoTUScY13gcpfKUYlmnZTQr04zcrrlt3bws9dREHJqmzdc0LUrTtBP/OP4fTdPOaJr2p6Zpk7OuiVkvPDyctm3bUq5cOcqUKcPQoUNJTk4mKCiIVq1a/ev8devWUbNmTWrUqEHlypWZM2eODVothBDZV6IxkbkH51Lxm4p0WdGFZSeXkWhUWQY7V+n83AdfSF8mrB+BFn8/oGlaY6AtUF3X9SrAFMs3zTp0XadDhw60a9eOc+fOcfbsWeLi4hg7duwjz09JSaF///6sXbuWo0ePcvjwYfz9/a3baCGEyMaCwoIo9WUp3l33Lp6unizvvJwTA0/g6uhq66ZZ1VOHoHVd36lpWsl/HB4ITNR1Pen+OVEWa9GjglmXLjBoEMTHQ8uW/369Tx/1deMGdOr08GtBQU+83bZt23B1daVv374AGAwGpk+fTqlSpWjcuPG/zr979y5Go5G8efMC4OLiQoUKFZ7+cwkhxAvsWtw1biXcolL+SpTPW56ahWoy0m8kTUo1ydZbiTIjo7mgywMNNU3bp2naDk3T6jzuRE3T+muaFqppWmh0dHQGb5d1/vzzT2rXrv3QMQ8PD4oXL8758+f/db6Xlxdt2rShRIkSdOvWjSVLljxUqEEIIcQD52+eZ8C6AZScUZLBGwYDUCRXETb02EBA6QD7CL5WzAj5dxldhOUI5AHqAXWAZZqmldYfkddS1/W5wFxQqSifeuUn9Vjd3Z/8er58T+3xPqJ9j/wH8LjjAPPmzeP48eP88ccfTJkyhS1btvDjjz8+032FEOJ5dvTaUcbvHs+KkytwcnCid43ejG4w2tbNelhSEixerGrh/vYblCpl1dtntAccDqzSlf2AGbCvOk/pVKVKFf6Zn/rOnTtcuXKFMmXKPPZ91apVY/jw4WzZsoWVK1dmdTOFEMLu6bqOWVcjglsvbuX387/zfv33CRsWxpzWcyjrVdbGLbwvNhYmTVIBt18/cHRUU5hWltEAvAZoAqBpWnnAGbB+6y0gICCA+Ph4Fi5cCIDJZGLkyJH06dMHd3f3f50fFxdH0N962UeOHJHqR0KIF5rJbGLFyRXUnVeXhUfV79IBvgO4POwyE16dYF/7eOPjoWxZ+OADqFoVNm+GQ4egzmNnUrNMerYh/QyEABU0TQvXNO1tYD5Q+v7WpF+A3o8afs4ONE1j9erVLF++nHLlylG+fHlcXV0ZP348AFu3bqVo0aJpX4cPH2by5MlUqFABHx8fPvnkExl+FkK8kP6+lajz8s7EJsbi4eIBgLuTu/1sJfrzT9XjBTWVOX68CrqbN0PTpmCjeWgpR2hlL9rPK4R4fjVf3JzNFzbjW8SXMQ3G0L5iewwOBls3S9F12LkTAgNh/XoVeM+eBW9vqzbjSeUIMzoELYQQ4gVzLe4aY7eOJTYxFoAxDcbwR88/2N9vP50qd7Kf4Hv2LNSrp7a17t8Pn30Gly9bPfg+jaSiFEII8UTnb55nSvAUfjzyIynmFGoXqU2HSh1oUqqJrZv2QEICXLkC5ctD4cJgNsPs2dC7N7i52bp1jyQBWAghxCMlm5Lpubpn2laiPj59GFV/lP2sZgaIiYFZs+Crr6BAATh2DHLlggMHbN2yp5IALIQQIo2u65yMPkmVAlVwNjij6zrv13+fofWG2tdq5rAwmDoV5s9XK5tffx1Gj7bZgqqMkAAshBACk9nEylMrmbRnEseuH+PCexconrs4yzovs3XTHmY2g4MD7N4Nc+ZAjx4wahRUqWLrlj0zWYQlhBAvsL9vJeq6oitxyXHMfn02BXMUtHXTHtB12LQJXn0Vptyv/dO1K1y8CD/8kC2DL0gAZvjw4cyYMSPtz82bN6dfv35pfx45ciTTpk2jatWq/3rv3r17eemll/Dx8aFSpUp8+umn1miyEEJYzPW46wzeMBhPV09WdF7ByUEn6VerHy6OLrZuGqSkwKJF4OMDLVrAqVOQJ496zcnJ7lY1P6sXPgDXr1+f4OBgAMxmMzdu3ODPP/9Mez04OJgGDRo88r29e/dm7ty5HDlyhBMnTtClSxertFkIITLq6t2rjNkyhs7LOwNQwrMExwYcY3+//XSs3NF+thKBqnLXqxcYjaqne/EivPOOrVtlMXY3B+z/o/+/jnWp0oVBdQYRnxJPyyX/LkfYx6cPfXz6cCP+Bp2WPVyOMKhP0BPv16BBA4YPHw6oykhVq1bl6tWr3Lp1C3d3d06dOkWe1Ceuf4iKiqJw4cKAKmNYuXLldPyEQghhfWdjzhK4J5CFxxZiNBvpUqULyaZknA3OVMpvJ8mBrl2DmTNh4EAoVgyGDoXu3eG119S873PG7gKwtRUpUgRHR0cuX75McHAwfn5+REREEBISQu7cualevTrOzs6PfO/w4cOpUKEC/v7+tGjRgt69e+Pq+mIVlBZC2L+VJ1fSeXlnnA3OvF3zbUb6jaSM1+OLzVjd6dNqRfPChaq3W6GC2r9bt66tW5alJBUl0KNHD1q3bs3GjRsZMWIEERERBAcHkzt3bmJiYhgwYACtWrXixIkT/3rvhQsX2Lx5M7/88guapj1UqOFR7OHnFUI833RdZ8tfWzBoBgJKBxCbGMvU4KkMqTuEgjntaHGV2QxdusDKleDqqoacR45UxRKeE5KK8ilS54GPHz9O1apVqVevHiEhIU+c/01VpkwZBg4cyNatWzl69CgxMTFWarUQQjzMaDay9MRSas+tTfPFzQkMDgTA09WTz5t8bh/B12SCXbvU9w4OKmvVxx/DpUsqc9VzFHyfRgIwah543bp1eHl5YTAY8PLyIjY2lpCQEPz8/B77vvXr15M6gnDu3DkMBgOenp7WarYQQqRZcXIFFb6uwBsr3yA+JZ7v23zPr2/8autmPZCQoPbtVqoEr7yiMlaBymA1bpzKYvWCeeHngAGqVavGjRs36N69+0PH4uLiyJcvH3FxcZw5c4aiRYumvT59+nRWrlzJ8OHDcXd3x9HRkSVLlmAw2NEKQiHEc+1Wwi2cDc7kcM5BXHIc+d3zM6XpFNpWbIuDZif9qzt31MKqr76CqCjw9YWlS0EWrcocsLW9aD+vEMLywu+EMz1kOnMPzWWc/zhG+I3ArJvR0NDsJRVjSoraq3vrFpQoAS+/DO+/D40aZat0kZn1pDlg6QELIUQ2cSr6FIHBgSw+thizbqZr1a40Ld0UwH56vIcOqRq858+rUoB58sCFC5A/v61bZnckAAshRDYxeMNg9obv5d3a7zLCbwSl8pSydZMUXYfNm1Xg3bpVVSN6911ISlKrmyX4PpJdBGBd1+1n2CQLWXO4XwiRvem6zsbzG5m+dzoL2i2gSK4izH59Nl5uXuTPYWcBbcUKtZ2oSBGYPBn694fcuW3dKrtn8wDs6upKTEwMefPmfa6DsK7rxMTESKIOIcQTpZhSWPrnUibvmczxqOMU9SjKhZsXKJKrCBXyVbB185Q7d+C77yBfPpUwo00blUSja1d4TOIi8W82D8BFixYlPDyc6OhoWzcly7m6uj60kloIIf4uISWBarOrceHWBSrnr8yPbX+kW7VuOBvsJKhFRsKXX6rtRLdvQ8+eKgC7uKjvxTOxeQB2cnKiVCk7mccQQggri4mP4ffzv9Ojeg/cnNzo49OHGgVr8Hr51+1nYRWoVJEffqgSaXTqBKNHqy1FIsNsHoCFEOJFdCn2EtNCpjHv8DwSUhJ4ufjLlPAswUevfGTrpim6Djt2QMWKUKgQVK2q5nZHjIDSpW3duueCHT1eCSHE8y/iTgS9Vvei7FdlmRU6i06VO3Fs4DFKeJawddMUoxGWLVOFEBo3VukhAZo3h6+/luBrQdIDFkKILKbrOreTbuPp6omroyubL2xmSJ0hDPcbTvHcxW3dvAdmz1ZbiS5ehHLl4NtvVT1ekSUkAAshRBYx62bWnV3HxN0TSTGnsL/ffvK65+Xy8Mv2s7Dqzh3w8FDfb98OBQuq+d42bUBS62YpGYIWQggLSzYl88PhH6g6qyptf2nL1bir9K7RG7NuBrCP4HvuHAwYoOZ3T55Ux378EYKDoX17Cb5WID1gIYSwsMXHFvP2b29To2ANlnRYQpcqXXB0sJNftyEhaph5zRq1Z7dXL8iZU73m7m7btr1g7ORfhBBCZF/X464zc99MyuUtRx+fPnSv1h3vXN40K9PMvhIMxcZCQIBKDzl2LAwZooachU1IABZCiAz669ZfTAmewg9HfiDJmMR7L70HgKujK83LNrdx64DERJWhascOWLwYPD1hwwaoUwdy5LBaM9YcjiBw0xkiYxMo4unG6OYVaFfT22r3t1cSgIUQIgM+2/EZ43aMw9HBkV7VezGq/ij7SRUZEwOzZqltQ1FRULu2Kgvo5QX+/lZtyprDEXy46jgJKSYAImIT+HDVcYAXPgjLIiwhhEgHXdfZ+tdWou+ptLl1itRhlN8oLg69yHdtvrOf4Lt7NxQvDh9/rDJVbd8OBw6o4GsDgZvOpAXfVAkpJgI3nbFJe54k5EoIE3ZNIORKiFXuJz1gIYR4ApPZxKpTq5i0ZxIHrx7k88af89ErH/Faudd4rdxrtm6ecuCA6uE2a6Z6u336wKBBUKWKrVtGZGzCMx23leDLwTRZ2ASj2YizwZmtvbbiV8wvS+8pAVgIIR7ju4PfMTl4MudvnqesV1nmtJpDrxp2kpjCbIaNG9WK5h07oFYtFYDd3OCbb2zdujRFPN2IeESwLeLpZrF7WGKOeeTmkSSZkgC1jSwoLCjLA/BTh6A1TZuvaVqUpmknHvHaKE3TdE3T8mVN84QQwroSUh4Eiw3nN+Dp6snyzss5Pfg0/Wv3x9XRDkqKrl+vcjO3agV//aUSZ2zfbutWPdLo5hVwc3p4T7Gbk4HRzS0zZJ86xxwRm4DOgznmNYcjnvi+u0l3mR4ynb9u/QXAO7XfwdnBGYNmwNngjH9Jf4u070nS0wP+EfgaWPj3g5qmFQOaApct3ywhhLCuiDsRTN87nXmH5rH/nf2Uz1uehe0WktM5p31sJbp1CxwcVKH7pCS1h3fxYujSBZycbN26x0rtiWbVKugnzTH/8x4hV0JYd3Yd4XfC+e3sb8QmxqKjM8JvBG/VfItK+SoRFBaEf0n/LO/9QjoCsK7rOzVNK/mIl6YD7wO/WrhNQghhNaeiTxEYHMjiY4sx6Sa6VumKQVM9tlwuuWzcOiAsDGbMgHnzYORIGDcO2rVT2ars4cEgHdrV9M6yFc/pnWMOvhxMowWNMJqNAPiX8GfiqxN5qehLaef4FfOzSuBNlaE5YE3T2gARuq4ffdqToaZp/YH+AMWL21HScSHEC+9O0h18v/PFrJvpX7s/I/1GUiqPndQnP3RIze8uX64Cbbduqg4vqJ6wAJ4+x3wi6gRVC1Rlx6UdmMyqp2zQDDQr0+yh4GsLz/y3qGmaOzAW+Dg95+u6PlfXdV9d133z58//rLcTQgiL0XWdDec2MHTjUAA8XDz4peMvXB52ma9bfm374KvrD77/7DM11zt8uJrnXbgQqlWzXdvs1KPmmF2dNBr7hNPox0ZUm12NfeH78C/pj6uj67/neBMTVQ7s5GSrtz0jPeAyQCkgtfdbFDikaVpdXdevWbJxQghhCSmmFH458QuTgydzIuoExTyKMfaVsRTIUYDWFVrbunlqTvenn2D6dFi5UpUCnDlTzffmzp2lt87uWapS2/rxxlVcid9PDledePdQ/rfvDMU8ijG9+XQq569MLpdcbO219cEcr0sZNZz/zTcQHQ158kDbtlZt+zMHYF3XjwMFUv+saVoY4Kvr+g0LtksIISzi6LWjtPmlDZdvX6ZK/iosbLeQN6q+gZPBDhYuxcaqmrszZ8LVq6qHGxOjArAVpuyyS5aqpz0kFMh7ifP6GJKdk4k1mSjtWJpF7RfRtUrXh/6e/Yr54edVXY0qLFyoHnxatVJz640aWf3nemoA1jTtZ8AfyKdpWjjwia7r32d1w4QQIqOi70Vz6fYlfIv4UtarLNULVueblt/QslxLHDQ7mT9NTFSB9sYNePVVNQzatKlVF1Y9ywpiW3nSQ8JLZR2YuW8mPx3/iWRTMibdhEEz8HbNt3mz+psPLqLrcOkSlCypKj4dPgy9e6tAXLGiDX4qJT2roLs95fWSFmuNEEJkwsVbF5kaMpX5h+dTPHdxTg0+RQ7nHKztttbWTVMOHYK1a+GTT1RFosmToWZN8PGxSXOyQ5aqRz0k3DFe4t11XxGr/YHRbMS/hD/R8dEkm5JxNjjTuGRjdWJKilrENnUqXLgAV65Arlywd69d1DuWTFhCiGzvZPRJvtj5Bcv+XIaD5kDP6j0ZVX+Ufezf1XX4/XeYMgW2bVMB4O23oWhR6NvXpk2zRpaqzEp9GEhyOEWiw3E03Z1bTnPQTE4MqPM2I/xGUNarLCFXQh7M73pUVp/3l19CeLjq5QYGqr3TYBfBFyQACyGyKV3XMZqNOBmcOBF1grVn1zKs3jCG1xuOt4d9DJ9y6pRKlHHiBHh7qx5v//5ZvrAqvUY3r/DQ8C5YNkuVJRTO7cKpuJ+IdVoEmNFwJKepBRXd32bW6x3TzvMr5odf0XpqCP/gQRg9Gho3VnPsr71ml1u3JAALIbIVk9nEmtNrmLRnEq3Lt+b/Gv0fHSt1pGnppuRxy2Pr5qmFVWFhali5eHHImxcWLIA33njQA7MTWZ2lKjOSTcn8dPwnwl0mEJt8FnRAA1034qYVZGyLeg9O3r9fDTPnyqUSltSuDadPQwX7eZB4FAnAQohsIdGYyMKjC5kSPIVzN89R1qsspfOUBsDgYLB98L106UHGKm9v1fvNkQOCgmzbrqfIyixVGbXq1Cre2/geEXcjqFGwBvULj2D52W8w6yk4aE6MaNSedtULwZo1KvDu3q1GFf7znwcXsfPgCxKAhRDZRL/f+rHk+BJqF67Nsk7L6FCpAwYHO5jLO3ECxo+HZcseZKwaOTLbpIm0F2vPrGVvxF5alWuFp6snFfNVZH7b+TQt3RRN0xh6pdPDeZo//VTt4y1RQj34vPWW6gFnI5r+98wrWczX11cPDQ212v2EENlX5N1IZuydwUDfgZTKU4oj145wM+EmjUs2tv3iKrNZrbB1cYGlS+Gdd+Ddd+G996BYMdu2LZs5feM0o7eMZt3ZdWhouDq6ProW77Vr8PXX0Lw5NGyoRhz27YMOHcDRfvuSmqYd1HXd91Gv2W+rhRAvrLjkOGp8W4ObCTepnL8ypfKUwqeQbbbqPCQpCZYsUcOe3bvD2LHQsSO0aGE3C6uyi5ArIUzaM4lfz/yKo4MjGho6+r9r8Z44oT7vn35SDz05cqgAXKKE+srG7G9ZmBDihXf46mFuxN9gaael9PHpY+vmqFKAEyaoRA5vv63K/1WurF5zdJTgmw4hV0IYv2s8IVdCAJgdOptdl3fxSaNP+PWNXx+dp7l3b5UdbNkyNcpw9ix8+KHtfggLkx6wEMLuhEaqqaqGxRvauCX3vf02rF6thj8XLYKAAJnjfQY7wnbQdFFTUswpuBhc2N57O4FNA5n9+mxyOOcAUHmaL/yB/wUzfoXrqDe+9JJaTDVgAHh52fAnyBoSgIUQdif0aijFPIpRMGdB2zTgwAE17DlhApQqpRb7fPopVK9um/ZkU7cTbzPn4By+2PkFKeYUAFLMKQ8PMQPcvInfou34fT1b5cQuWFvlaB40yEYttw4JwEIIuzO4zmDaV2xv3Zuazar835QpsHMneHioed5SpaQMYAYYzUaqzq5K+J1wfIv4cvz6cYxm48NDzHFx8MEH8MMPEB8PzZo9yIn9ApAALISwO/WL1bfuDY1Glbzh2DG1innaNDXs7OFh3XZkAWuUG0xNA1nCswSnok/xWePPcHRwZNKrk6iQtwK1i9R+kCqyRCP8tPsrxd3d1cNOly4wYsQL96Aj25CEEHbl4q2LnIk5k1ZAPcvExMC6dWqhD8D//gelS0OnTmqR1XPgn5WEQKWanNChmsWCcMiVEBovaEySKQkAZ4MzRwccpWK+f1QZMhrVPPqUKXDunCqMkCOHOm7H24gy60nbkGQVtBDCrqw4uYLXlrzGveR7WXOD8+dh8GDV0+3TRwUDUFuKunV7boIvPLncoCVcvn2ZN1a+kRZ8NTRG1x/9cPC9e1cVRShXTvV0Y2Lg888fFER4joPv00gAFkLYlQORBx7GUNMAACAASURBVCjlWYq87nkte+ErV9Se3fLlVbrIbt3UHtNy5Sx7HzuSFeUGk4xJnIg6AUChnIXwdPXEycEJg2bA1dGV18u9rk5MHV39808YNkyl51y1Cs6cUQ9Arlk4upFNvLiPHkIIuxQaGYpvkUeO2D07kwkiI1VvN3duVY/3ww9hyBAoXNgy97Bjlio3GHIlhI3nN3Ij/garT6/GoBn4a+hfacPND5UCvOEKH76pPu9vvoF69eDoUVlB/ggSgIUQdiMmPoaLsRcZ4DsgcxeKj1cViKZNU+kijx9XC6rOn7ebWrDWYIlyg7+d+Y2OyzpiNBsBqOtdl/FNxuPk8GCo3s/7JfyO3YTeY2H7dsiZUz3kpJLg+0gSgIUQduPg1YMA1ClSJ2MXuH5d9bpmzVJzjXXrqrqwuq4SZ7xAwRcyV27QrJtx0BxYfWp1WvA1aAbaVWhHQOmAtPPWHI7g5qgPeWvbIqI88hE1dCxVPx0Fnp5Z80M9R2QVtBDCbhjNRk5Gn6SsV1ncndzT/8bUALtokVrV3KYNjBoFDRpIxqpnoOs6uy/vZnLwZKrkr8LEVycSfDmYgEUBpJhScDY4q0IJrmVh9mx2lPBhwHln8keFUzPyNOsrNsTJ1cWiq6yzOynGIITIFhwdHKleMJ3Dlbqu9pBOmQKNGqmA27WrSl9YvnzWNvQ5s/vybr4N/ZYj147wZ/Sf5HPPh38JfwDqF6/Ptl7b1ByvQyn8vvgRFi6ExERON3uLhJoduJynMJfzqDl14/1V1tkpACckwOHDsH+/+vOwYda5rwRgIYTd+O/W/9KibAteKfHK408yGmHlSpUq8sAByJdPZVACcHaW4PuMUvfxGs1GNDRG+o3ks8afPTQC4VfMD7+P58KP/1Vz6r16wfDhTFzw1yOvmZlV1lnNZIKTJ1WwTf06flwdB/X8JgFYCPFCuRZ3jQm7J5DPPd+TA/Bbb6mh5vLl4dtvVTBwe7ZVvc+LjGa5ik2M5dvQb2ldvjVBYUGYdTMADpoDed3yquCbkgK//Qbt2qm58+rV4ZNPVH7mAgUAKOJ51SKrrLPS1auwd++Dr4MH4d79Lea5c6tlAh98AHXqqK8iRazXNgnAQgi7cDDyMQuwwsPhq6/UqtpixdQe0k6dVLJ+hxc3lcE/s1xFxCbw4arjAI8NwuF3wpmxdwZzDs4hLjkOg2bAv6Q/LgYXkk3JKk9zPl8IDISZM9Vnv349tGwJw4f/63qWWGVtSUlJaig5NdiGhMDly+o1Jyfw8VHPb3Xrqq+yZR/8E1pzOILOC7M2Zec/SQAWQtiF0MhQNDRqFq6pDhw9qoaZf/5ZFUqoWhV69lRjhOKJWa7+GThCroQwYtMIDkQeAKBr1a6Mrj8an0I+wP1SgGc34b/xFH51O6giCY0bw+zZ0KLFY9uQmVXWlnDtGgQHP/g6eBCSk9VrxYurLcjDhqn/1qz5+NwfGXmYsQQJwEIIuxB6NZRK+SuR0+CmelwbN6pcwYMHw9ChqiqRSPO0LFe6rnMg8gBGk5FXF71KojERg4OBpZ2W0qFShwdvuH5dzfEWqQuDq0DbtqowQq1a6WpHu5reVgm4JpNKqrVnz4OA+9f9KWgXF/D1hffeg/r11TPaswwlP8vDjCVJABZC2F5SEtcizuBb1k/NN1asqFY29+8PefLYunV26XFZrgrndmbVqVVM3jOZfRH76FezH8mmZHR0dF3nzI0zKpqtXatGGM6cgUuX1Dz6sWNqIZsdSExUC6R274Zdu1TAvXNHvVawoNphNmiQCri1aqkgnFFZkbIzPSQACyFs5+ZNmDMHvvqKA1evknx8mTo+bZpt25UNpM6/xppOkOhwHBdzJTTHSMKc1tFx2UVK5ynNrJazqJSvEkuOL3kwx3soBt6qqLKClSgB//3vg4vaMPjeuqV6t6kBNzT0wXBylSoqdXeDBuqrVCnLbu+2VMrOZyUBWAhhfdHRqiLO998/VIjduUoNW7cs22hX05vTNw8ydvdHmPUUHDQn3BxdqehRli9fm0SHSh0wOKjMX1t7bVX7eG974td6kFqBtHQpdOhgs2pEUVEq0O7YobZzHzumtnY7Oanh5GHD4OWXVQ83r4XrcvyTrRaTSQAWQljP7dtq74ejIyxeDJ07w4gRfBm/nX0RP7KEpkjeqvS5fPsya8MmYSYJNNA0I0NeGsCEgAloqd3DEydg2jT88uTBb+pUFeFC66oxWytnCIuIUIE2NeCeOqWOu7urIDtuHLzyino2sPauMlstJpMALITIWiaT2k86darq7R48qOZ1UwuyAxsXv8/VuKsPAod4rOPXjxMYHMjPJ37GbDZj0FQv19ngTNsKbdUDzJYt6vPetElFs8GD1Zs1DWrXtko7r11TdRmCgtR/U8sue3ionm3v3mqav1Yt+5h2ttZisr+TACyEyBp/r0h0/jyULKnGFU0m1QO+H3x1XSc0MpR2FdvZtr126u+l/nR0GsxvQA6nHAypM4Rh9YYReTfyQSnAYn7w8cdqeL9QIfjiCxgwIOvHcFFDykFBDwLu6dPquIeH6tkOGKACro/PC1cT47EkAAshssayZWqZat266vv27R8533jp9iViEmIsVwPYQjKaZcqSdl/eTcBCVQjB1dGVLT23MLPFTHpU74GXmxcAJcy58NsdBO7OUAzo0UOtUurePXNLg5/izh01nLx1K2zbptI5gqpE2LChSnjRuLEKuDaaZrZ78rEIISzj5EnV261VSwXebt1UqqGnVCQKjVQV0uwpANsqMUOqhJQEFhxdwNitY0k2qaXAyaZkdl7ayYcNP1QnXbgA06fDDz+o0QZdV8PLFSqoLwtLTFSZpbZuVV8HDqjBDFdXNaTcrZsKuLVrq4VU4umeGoA1TZsPtAKidF2vev9YINAaSAYuAH11XY/NyoYKIeyQrqvxxilTVOIMV1eVgghU7+vll596CYNm4CXvl6hWoFoWNzb9bJWYAWDlyZUMXD+Q6PhoKuWrRFxKHCazSW0hKumvTurXD+bPV13LHj1U4oxqlv38zGa1MnnzZjWlvHu3CsIGg8qZ/MEHEBAAfn6PzzAlniw9PeAfga+BhX87tgX4UNd1o6Zpk4APgTGWb54Qwq4NGABz56rk/J99BgMHqupEz6B9pfa0r9Q+ixqYMdZMzBByJYQ1p9dQv1h92lZsS1GPovgW8WVMgzG8UuIV9obvJeivbfhHOuPnfT8NZ7ly8OGHKj924cIWa0t4uAq2W7bAH3+o3WIAlSvDu++qgPvKK2ohuyXZw3C/LTw1AOu6vlPTtJL/OLb5b3/cC3SybLOEEHYpNha++071uooUUeOOderAm29mqBuk6zpm3Zy2X9VeWCsxw6Kji+j7a19MugmDZmBX3134FfNjQ48N6oS7d/FbtR+/GfMgLAy8qkPz5jDGMv2de/fUPO6mTaqnm7pwqmBBtTW7aVN49VXwzsJYaOvhfluyRCmRt4CNj3tR07T+mqaFapoWGp36OCWEyF7CwlQ1nGLF4P33Yd06ddzfXw2HZnAM8vzN8+SemJu1Z9ZarKmWMLp5BdycHn4osGRihh1hO3htyWv0WtMLk/5gqDsoLEh9ExengmyxYmrluLc3rFqlomEm6LoaVg4MVJfy8oLXX1fPVCVKqJmEo0dVCb/Fi9VWoawMvvDk4f7nXaYWYWmaNhYwAksed46u63OBuQC+vr56Zu4nhLAys1n1bpcuVXXbunaFkSNVaRkLCI0M5V7KPYrnLm6R61lKViRmMJlNaT39pX8u5dDVQ7xb+10WHl34IE1knvufq6srrFihuqEjR2aqAtTNm6p3m9rLjYxUx6tWhf/8R3WoGza03TyurfIwP8RotMlS7QzfUdO03qjFWQG6rktgFeJ5YTbDvn1qdY2Dg9qvO3KkKjVTtKhFbxUaGYqroyuV81e26HUtwRKJGUKuhLDlry3cTbrL6tOrWdBuAQ2KN+CLJl8wrfk0XB1d6V29J0Gb5+K/7gR+37wFFy+qBWwnTmQoJZTZrGribtwIGzaov0qzWeU+adpUVRds1izre7bpZas8zIDq6s+erdYx7Nql5tatKEMBWNO0FqhFV410XY+3bJOEEDaRmjhj+nSVtujkSahUSY1PZpHQq6H4FPLByfD87VvZdH4TrX5uhdFsBKBSvkppr3m5eaklxT/Ow2/aNPxOnVIRMTVRCTxT8L11S/VuN26E33+H69fVzi9fX/joI3jtNTVVb48JMGySh/nQIZgxA375RfV+W7WClJSsu99jpGcb0s+AP5BP07Rw4BPUqmcXYMv91HF7dV0fkIXtFEJklVu3VNCdNQtiYtRv6l9+yfLegMls4tDVQ/Sp0SdL7/M4Wbny1qybeWPFG2nB10Fz4M3qb9KgeIMHJwUHwzvvqEwVqXmx05mTUddVbdz169V0fHCw6uV6eakh5ddeU/8tUMAiP06Wsnoe5thYlXza0VGt4v/Pf6ze802lWXP02NfXVw8NDbXa/YQQT5CYqCb+btxQaSJffVUNNb/8slUS9cenxDN5z2QaFm9IQOmALL/f3/1z5S2oXteEDtUy/Iv/yLUjLDiygCnNpmBwMDBx90TG7RhHiikFZ4MzW/3n47dwu0oLOX68iqIhIWqoPx2fd2Ki2nK9bp0KvJcuqeO1akHLluqrbl377OXa1J07as/0/v3w00/q2KZNal7d0zPLb69p2kFd1x+ZZUYCsBAvEl1XeQOnTlWrc0JC1C//mBir5Au2Fw0mbnvkvKO3pxt7PmiS7usEXw7m+8PfcyL6BPsj9pPTOSfBbwVTraBKihFyOZigoB/w33AKv6V71Nzu4MHq80+HyEgVcNetU9mn4uNV9aCmTdXq5ZYt7Wcu1+5cuABffaWC79276sFy/XqVnNqKnhSAJRWlEC+C5GQ1rDxtmtpnUrCgSuKQWhjBBsH3/M3z5HfPT25XC2d1SAdLrLz99fSvtF/aHh3ViRlQewDjA8aTxy1P2jl+837H7/N5KjnJJ5+oFJ1PGBfWdfXXs3atKiCV2l8pWVLlVm7VShU0kMxTT7Fhg/qwDAZ44w0YOlRNiNsZCcBCvAgWL4a331Ypjb7/XiXqt/Fv8V6re+Ho4MjOvjutfu+MrryNT4nnZPRJfIv4ciLqRFrwNWgGiucuTp5kB/hmilpmXL262rZVtCj07PnIRVVrDkcwaf1ZLh5zR7vijTGsEDeuOaJpaoR0/Hho3RqqVLF6+d7sJSlJPWDmygUdOqinlI8/hv79VcIYOyUBWIjn0V9/qVWePj6q69Stm/pF1Ly5XfwmN5qNHL52mIG+A21y/2ddeRsTH8M3B77hq/1fAXBl+BWalGqC2y43tYfXwQn/347C68XUcGdSEmtMeQncdJ3IWG+KfBny0MKi2Fj4fPZNvlvsQNz5hujJjmhORnKUjmFwf2f+b1AeChbM+s8h24uKgm+/VQsIr1+Hdu1UAM6RAz791NateyoJwEI8T0JC1Pzi6tVqD+/776vjbm5qA6idCI0MJdGYSJ0idWxy//SuvF1zag2Tgydz6OohkkxJtCrfijENxuBicMGvmB9be20laMZQ/FccxC9iRVqikjVawX+lVxz141k2LsvBxYOebN8ORqMXDjkSyVEpErdy13EtfgMHJzOH3dwoWDD989AvrMBA+L//U73fli3VFq5MZApLMaWwPWw7zco0s2Ajn0wCsBDPi8GDVU/A0xNGj1bbK+xwhc6MvTN4f8v7uDu507BEQ5u140mJNkxmE/sj9vPGyjdIMiVh0Awsbr+YHtV7qHnzzZuhWTP8ivnh5/E6dGmkEpUUKwZA4MRtJKSYSLmRk/hzBYk/V4jkq55cAMqXV4vNF0TswblI7L8GJKyaASo7MZvV3G69empOvWxZ6NtXze9WrJipS286v4l31r7DlTtXODrgKNULVrdQo59MArAQ2VVcnKoF27WrWtjTtq36RdS3r6qKbkcORBygqEdRCucqTLUC1Xin1juMrD+Soh6WzayVGbqus/XiVibvmUz5vOXxzuWdto8X4PKN8+oBZ/p0OH9elQsKCFCLq9KuoXI8/Plbce6dKYTxpvp7cC58C89Gp3Evd40z3/kDsGtiEhGPKOJqlQxQ2UlcnEoQ8+WXKkHMlCnqCaZ9e/WVQRF3Ikgxp1DSsyTeHt6U8SrD7NdnU7VAVQs2/skkAAuR3UREqO0Vc+aoyURnZ1Urrlkz9WUnUgPaxN0T2XpxK2MajGHiqxMJKB1g9X2/jxNyJYRtF7dh0k2sOb2Gw9cOUyhnIV4v9zp1vevibHBWc7xmDf//TIWTdx8kKmnUCFAd4uBgVSth1Sq4fBlwKI1rsZt41A7Drdw1HHMlAWqbUyqbZIDKTnRdFR2eMwdu31ar0n7+GTp2zNRlj18/zpSQKfx0/Cc6Ve7Ezx1/pmqBqmzvvd1CDU8/CcBCZBcmk1rJvGSJGo5r3171BPz8bN2yf1l1ahXjd43n4NWDFM5ZmMCmgfSv3d/WzXpIyJUQAhYGkGhMREenuEdx5rWex5vV38TF0QVu31ZzvBe24j9yJn7l6sO3KlGJ0aQRFKTqJaxZo9b/uLio559x48ChxHUmbDv6xOBq9QxQ2YGuw9mzUKGCWix45oxaODh8uBp6zoQdYTuYsHsCmy5sIodTDgb5DmJYvWEWanjGSAAWwp6ZzWpM09dX7Wk0GtVe0qFDoXRpW7fuIUazEUcH9StlxckV3Em6w3etv6Nn9Z4qoNmJG/E3+Gb/N0THR5NsSkZHx0Fz4F3fd3m75ltqaHnqVDh5Er8LF/Ar5gc7hpHikpNt22BFf7XGLSZGLbZt2VJ1ylq2VLtglMJ4eJqfGlwtUfDhuZCSAitXquH9AwfUEH/p0upYJlJ7pZhScHRwRNM0NpzbwJFrR/hfk/8xwHeAysdta7quW+2rdu3auhAiHeLjdX3uXF2vVEnXNU3Xz561dYse63bibX3y7sl6kalF9GPXjum6rus342/qRpPRxi172F83/9KHrB+iu33hpvMp+sB1A3W3L9x0wziD7vaFmx787Ue6Xr26roOuFyqk6198oSfevKevW6frffroep486qVcuXS9e3ddX71a/TWJTLh9W9cnTtT1okXVh1uunK5//bWux8Vl7rKJt/WpwVP1YtOK6RvPbUw7lpCSYIlWPxMgVH9MTJQesBD2JDZW7d+dNQuio9U+3oULVbV0O3M97jpf7vuSWQdmcTvpNq+WfhWzbgZ4KBuUrYRcCSEoLAj/kv4sPLqQuYfmYtAM9Kzek5H1R1I5f2V6Vu+pzrnuil/nEVC1Ksa58/mjQHd+We3CmlJq+jF3brXGrVMnlQZSMlFlUnKyWrtw757aStSwoSoL2LKl2j6XQRF3Ipi5byZzDs7hdtJtGpVohIeLSj2Z+l97IrmghbAHSUlqEvH6dZV3MCBAze/6+9tF4ox/SjImUXR6UWLiY+hUuRNjGoyhdpHatm5Wmol//Mp/93RG1004aE40LdGFGt6Fee+l9/D28FZ5gmfMSEsRmZKsc3jGDr491YjVazRiY1XQbd9eFSl69dV0FyoSj6PrKqH1jBkqqfW2bep4RIRFtsvpuk65r8pxMfYinSp3YpTfKOp422af+d9JLmgh7FFqYYRp01T2pJ07VY7mS5fsso7ckWtHWHpiKeMDxuPi6MKslrOoUagG5fOWt3XT0hjNRsasn8OXB/8PXUsBDcx6Cof+MjCgxnt4n7gEU9+D1avRHR258vpAPn8HVq3SuHnTHw8P1dPt2lX1dCXoWkBioqpCNGMGHD+u/m0PGqQWFRoMGQ6+uq6zPWw78w/PZ37b+TgbnJnbei4lPUtSOo99rY94HAnAQljbPwsjFCigkmik/kKyo+Cr6zo7Lu1g4u6JbLqwiVzOuehXqx9lvMrQuUpnWzcvTZIxibkH5zJt7zTCYsMw6PlRv97MaDhiMFbh1qgPYNtiUnLlYVutDxkdNpjja4qQMye0aaOCbrNmMrxscd9/rwp/VK+u9q1366ZGezIoxZTC8pPLmRI8hcPXDlMwR0HO3DhDtYLVaFIqe2UQkwAshLX98IMqBF65MsybBz162OVv/bDYMN5Y8Qb7IvZRIEcBxjcZz8A6A/F0zfoaqumRuoe3Sakm+BbxZWrIVIp6FCX+ek9czXXQTMcoE/0rd9xeITzhJb6/XoELnvWYHtsH8585aN0aPumqph0fUSdBZNSRI6q3GxCgilD06qX+rVtgOiX8TjgN5jfg8u3LVMxXkXmt59Gjeg9cHe3v/5/0kAAsRFY7f15l8alTR/0y6tFDLaqyk8IIf5dsSubCzQtUyl+JwjkL4+jgyOzXZ9O7Rm/cnOwnSi3/czndV3bHqBtx2+XG1l5bOfDOAfLnyE+7D5fSImgB3Q7/Tu7ke3zs1oLPExpy3cGM52sOzO2merwPtgyJTDObVdHi6dMhKEgVLa56P6NUrlzQuHGGLx15N5Ij147QslxLvHN507R0U9pWaMvr5V/HQcv4gi17IAFYiKyg67BnjxpmXrPm4Zq7OXPaVWEEgLtJd5l7cC7T907H0cGRc/85h4ujC7vf2m3xe605HJHh5BMHIw8yOXgyy/9cnlYKMNmUTFBYEH7F/LjXexArFn+HZjazko5MYzhHC5SjUJUTTB7tRU9/+y1Nl6116aL27BYrBpMnQ79+kCdzK+FPRJ1gashUlhxbQk7nnESOjMTV0ZV5beZZqNG2JwFYiKzQv78aXvbygv/+V83xFi5s61b9S/S9aGbum8nXB74mNjGWxiUb07DQ2zSavJOrtxMtnp1pzeGIf1UJ+nDVcYCn3iM0MpQ639XBw8WDHtV6sOLUClJMKTjjyN3jjQj4FF7d5okbg/m17EDCKmgYS16mVtFj938GywbfzDxIZHuXL6utcmPGqEDbv79aLt6hAzg5ZerSx64f44M/PmDj+Y24O7nzbu13Ge43PNsOMz+JBGAhLOHOHbXYpFcv1dPt0AFq1oTevVW6JDuj6zqaphESHsL/dv2P9pXaM6bBGCKjvO8HyETg2QJkegRuOvNQekaAhBQTgZvO/Ov6uy7t4uv9X1M4V2FmtJhB7cK1+b7N93Sq3AmXRGcCdhTl1J/f0+5MNB+EpxBZFpI+GU+HbjAsLeNj1qzQzsyDRLa2d68aZl65Uv25Xj1VgzeTOciNZiN3ku7g5eaFyWzi4NWDfN74cwb6DiSve14LNNw+SQAWIjMuXYKZM+G779RWIi8vFXRfe83WLXuko9eOMmnPJMp5lWNc43G0Kt+KM0POUC5vOQAaLNqW7gCZEY8rtff343HJcXy07SNm7puZliayc+XO1C/WgIrRXQh6eRr1Dn5DHz2K444+HAmYxtSlftT2s96U+rM8SDwXEhLUoqqQELVBevhwVe6yePFMXTYuOY7vD33P9L3TeaXEKyxsv5CahWtyZfgVnA3P/x4wCcBCZITRCG++qbLxg9rDMny4ytlsZ1K3Ek3aM4nfz/9OTuecvF//fQAcNIe04AvpC5CZUcTTjYhHXCu1BN/Kkyt5Z+073Eq8lfaahsanc//gwk8NiLzoSBiziSziy6X3RlJzRGOqOVl/IVtWf052ITYWdu2C1q3VMvFq1aB7d+jTJ9PlLq/evcpX+79iduhsYhNjaVi8IV2qdEl7/UUIvgDZewmZENZkMsG+fep7R0eVMm/ECLh4UVUossPgCzB221gaL2jMoauH+F+T/3F52GX+r9H/PfLcx9WitVSN2tHNK+Dm9HByfYPTdXq+rOb3yuUtR6OSjZjqPwcn3NDMDhiSdfovmEv50ka+W+hKzvDT+ESsp86YJjjaIPhC1n9ONnX+vOrdFi2qplKiotTxOXPUfl4L1JoODA5k4u6JBJQKYO/be9nZdyetyrfK9HWzG+kBC/E0qYXvZ8yAsDCVxrBkSZXdxw4lGZNYcnwJLxd/mfJ5y9O5cmeK5y6erq1EWV2jNnV49uONqwiL34TmHEGcfow90d15K7k+f4VUx/DTMg6vX8G3hYtzveQZ6l3PjW/nfnT+IhncHYHcFmlLZjyXtXzPnYNRo2DtWvWA2b07DBuW6cQwuq6z89JOAoMDGV5vOAGlA3i/wfsMrjOYMl5lLNT47EkCsBCPc+MGTJnyoPC9nx9MnKh6BnboTtKdtK1EkXcj+fiVjxnXeBw1C9ekZuGa6bqGNWrUXklexZ/GYZidzKBDS++eOAdNwru/+sh7eG5mkak7iY4VcO05RyVzsLNMGc9NLd+kJFX0o2hR9RkfOAAffaRSRRYqlKlLG81GVp1axZTgKRyIPEA+93xE3VO96UI5M3ft54UUYxDin+LjVSKByEgoU0bNgQ0fbpeF71N9vuNzpoZM5XbSbZqUasKYBmNoWropmp0k+vh7XdaGPzRk9+X7+4vNBnJvG8G4YCP5qhTAc8IHNG9qxnH7FpWMOROVccQTREfDt9+qrURVq8KWLeq40ah6vxbQYH4Dgq8EU86rHCP9RtKrRi+7SuZiLVKMQYinMZth/XpViN3BQRVJKFJEVWrxsoPC3Y9w5fYVinoURdM0ouOjaVqmKe/Xf98uKsCklgKs612X41HHmRYyjdmvzSPpZDNStnwCpduAQzIuZlh/eSr1NQ2t0UB4HcBBZQkTlnf6tEoOs2iRKpLQooV6uEyVieB7Pe468w/PZ1T9UTgZnBhcZzCj/EbRpkIbDA6Gp1/gBSQBWLzY7t1T9XanT1dzYMWKwdChKiA7ONhl8D109RCT9kxixckVBPUOomGJhnzZ4ku76e2GXAmhycImJBmT0rJVeae8Qo/Oubh9AgoXfpWZsW8Qd/cH/KPd8esyCN57T332wvJ0Xf17NhjU/O6iRWq/+tChKkdzJp2+cZppIdNYeHQhyaZk6herT6OSjeherbsFGv98kwAsXmzz5qmFJnXqqApFFsjkkxV0XWfrxa1MGWT5PwAAIABJREFU2jOJP/76Aw8XD0bXH522hchegi/A9rDtJBpVIg90IHQAdzdNZUq1BZSfm5cGfcvjeGQQ7KyqUhZ62F+h9OdCQgIsXqwWD37wgZpLHzAA+vZVdZDT6XEZv24l3KLPr3347cxvuDq60senDyP8RthVeUp7JwFYvFiOHlVDcE2bqn28fftC7drQoIHdFUb4u0RjIt1WdsPRwZFJr07i3drvktvV9quBU+2P2M93B+fR0fUbtv/cGIq6gIMRB7MTc10T6etRHIfDMRD1BTiOVVu27HTbVrZ37Zqa2509W61q8/F5kIc8V65nqkLxz4xf4bFxjFi1DmhFG5/C3Eq4xcevfMzguoMpkMN+ymhmFxKAxfPPbIbff1eBd+tWlRqyenX1mocHvPyybdv3CAkpCSw4uoA1p9ewvvt63Jzc2PzmZirlr/TEnLjWyk8cciWE7WHbcTG4sOL4WvZe24FDkifzfhhEniQ/OvbeTsXkEby+MRS/sAWqyv3IkepBR2St1q3h4MEHiwcbNcrww2Vqxi8zidwzbOWO4xrM2h0m/u5Nu5re7Oizw65GX7IbCcDi+ffGG7B8OXh7q21E/ftnulJLVrmVcIvZobP5ct+XRN2Loq53Xa7fu06RXEWeupXIWvmJU+d404aZ4/LD7mk0zPU2/x3wJ7HNwvlyRwItVhUismBLtkwfRdO2DS12f/E3ZjNs2KC2yi1Zoh4oZ85UQ8zlyj39/U9xJfYadxzXcddxPWbtDs7mCngm9+Faohmwr6mP7EgCsHj+XLumtlgMGwaenmqYuV07Va3FDud3Ux2/fpz68+sTlxxH8zLNGdNgDP4l/dP9Sy6r8xPfTbrL2gNH+HbjbhIdU1QePbMDjV0HsbRTXvIvegU+PUr3i1OIKFSRL5r0A8At9B4Tikdkvz2y9uzePViwQNWZPntW7eM9e1YN61tgu5xZN+OgOeDpcYsrKT/jZnoJD2N7XMxV0NDwfh4yftmBpwZgTdPmA62AKF3Xq94/5gUsBUoCYUAXXddvPe4aQljF8eNqNfOSJZCSooaZO3Sw28IIAKeiT3Hu5jnaVGhD5fyV6VezH719euNTyOeZr5VV+YnDblxj6E8z2RA1G6PJjMPy1Rh6OKOTjIuDA//77RvyH7kBVarwv46jCc33cHaj57pIgS1ERUHFinDrFtStCz//DB07WuThMuRKCIHBgeR3z8+c1nP47LV2jFyVE1PKg0Vb2T7jlx1Jzy73H4F/Vg//ANiq63o5YOv9PwthG4mJqhxa9eqwdCm88w6cOaOCr50KvhJM21/aUnlWZYZsGILJbMLgYGB6i+kZCr5g2fzEIVdCeOvn9ykzti2lZpbgt5sTcYkMYGCuzfzwZSkqu0wib2JX1ix2ooxrZTXHfvw488o2Itnx34HguSpSYEVrDkfQYOI2Wvf5ksnth7PmcIRKDfnee7BnjyoP+MYbmQq+Zt3Mr6d/5eX5L1N/fn2CwoIo6qGyvbWr6c3UDq/i7emGBnh7ujGhQzV5mLKQp/aAdV3fqWlayX8cbgv43/9+ARAEjLFgu4R4soQEVRqtSRNwdYWCBWHCBDW/a4d7d1PtC9/HqC2j2H15N15uXnz8yscMqTvEIokKLJGfOCEBPl+4k0kRLTBrSfw/e+cdVmX5xvHPy4HDEBUVFQeIW5w5SsER7tTcK/dIrdyLbGhWViog7pFbf+ZKUdOyzIEL3Lhxg+IWBUHmGe/vj8eR5WAcOBx9PtflFQLnfW5O+N7vcz/3/f1iY6RIXGt+rO9Hj0b3uD3+O5ICTjH+4zk4KCX4vPmHJOfKw8QCFWmjKK91O5Kkno1HrrHddyEBIYHUvH6GB/a5aFj2fQDafPutydb5Lug7vt/zPe5O7kz/YDp9q/bFUfvMcKFN1SIy4WYS6T0DLqiq6i0AVVVvKYry0v5zRVEGAAMA3DLoHSmRPDdiERMD165BoUJCXCCbkmJIIUGXgJOdEzqjjmsPrzGt6TQ+rvbxcze6jJIRfeKzYUa+WryV32N80Vs9ApcUsDKiUTQMKpWDXj/2hAMHcLTPyYYqH2CrTyHJxo4Y+1zwjxLzG2lSYA727uXdNp1p8+AW13MVYEL9j1lTpSmPsMlwOf9+wn3mHplLw+IN8XT1pPc7vfHI70GH8h2wtpJtQVlJpr/bqqrOB+aD0ILO7PUkbyiRkTB+/LPz3ZYthRVgBgXjM5NHKY9YcHQBAQcCaFmmJXNazKGOWx0uD72caTe6tOxWkpPh18AUfti0ivP5/KDAGRxsXWnr3pEtd8JIMaSgRYP3TytBWwJmzsTzigvx2v/uZp+UmN8YkwJzcO2a0CEvVw5cXbnpkIcf6vZiWxlPDP+okKS3nB8eHU5ASACLjy8mQZeAUTXi6epJ8TzFKZ6nuKl+CkkaSO9d4I6iKIUe734LAXdNGZREAogRi+hoISKgqrB+vVBOGjYMymRftZ278XeZcXAGcw7PITopmveLvU/rsq2fft3cu4x1B0OY+VsQJ3/zJqbgZqg7ERelEmPr/48B7nWwmT2PkHz9CarlgrdbXTyrRkPz5qDR4DRpJ/GvKTHLkmUaOXRIzKivWycaBjdvBnd3hg+aabJy/uA/BjP3yFw0ioZulbsxynMUFQtUNEX0kgyQ3jvBb0AvYNLj/24yWUQSyRMJvalTwc1NNPi4ucGtW8KlKJszbuc4FhxbQFuPtnzu9Tk1i9Y0d0gYDGJc9NuVWzhWui1YG7FqZ8u4cr9Qq/pWmiW4oAQEwKo+YDTi+cknePp8KV5c7Nl1ZInZhGzdCj/8AMHBYn53xAgYMuTplzPyXhtVI9uvbKdh8YZorDSUyFOC0Z6jGVpzKEVyyYej7EJqxpBWIRqunBVFuQ6MRyTetYqifAxcAzpmZpCSt4Q7d8T57pw5QkKvalUhF6mqQsknmybfozeP4hvsy8haI6lZtCZf1/uakZ4jKets/qR05w4sWgSzVp3nVnF/qLoEFAMooCgp2Bc9R/NfkuD774VC2KBBosJQ/MUlSVliziCxseL32NoaQkPFQ+W0adC3738kItPzXifrk1l5aiX+If6cvXeWDZ030KZcG0Z6jszUH0uSPqQfsMT8PEmwkyfDl18+O9+tVy/b6jOrqsrfV/7Gd78vO8J3kMs2F3Oaz6Fb5W7pvqapZCRVFfbtE31q69aBruknUH0BWitbmpVuxl+Xt6IzpKDVaNnRayee14G9e0UHuZNTuuOXvIKICKFQtXAhLFgAnTuL8TkbG+FSlEGS9clMPTCVGQdncOvRLaoUrMJor9F0rtAZG032FZ95G5B+wJLsh9EIf/0lzr769IGuXeGTT6Bt22x9vgsi+Tb+X2N2hO+gkGMhk5gjmEJGMj5e9KhNXrmfK04LcLjwMZ99VhfHhhXR5BzL4JJdKLBsHSG/BhGUOwnvpj3xdPUEV0yiniR5ASEh4ihl/Xphb9mx4zMLQLuXa3qnlgRdAg42DlhbWbM4dDEVC1RkWZtlNCrRSMpEWgAyAUuylsREMTI0bRqEhQnTe8PjMy4np2y7A0vQJbD2zFp6VumJlWJFO492dKnYhe6Vu2NrbZvh62dERvLSJVG1X7Q0hdh3vgfvn0BR0VdfxUe9g/B0HSLOF3+uDomJeDZvjueoUVC/fobjTg8JugSWHV/G/sj9zPtwXrpGsbLKdCJDGI3i4fLOHfDxgcGDhWSkCTh++zj+wf7sCN/B5aGXcbBx4MiAI+SyldaOloRMwJKspVkz2L1bnO/+73/QqRNoteaO6qXcT7jPrEOzmHloJvcT71MsdzHqF6/PwHcHmnSdtMpIPjF4mjULtm5LwcpzBrYDp4HNjaffYzAaCIoIErvclBTo1k0kYhOYsGeElqtasjN8JwDtPdrT1qNtml6fVaYTaebhQ1FiXrFClPQdHcUZgLu7+DiDPDn28Av2Y/uV7ThqHRlQbQDJ+mQcbBxk8rVAUiNFKZGknxMnRGk5Lk78/euvIShI2KV1755tk29scixDtw7FbZob3+7+Fi9XL/b22Yu3u3e6r/lEVrD4F79Te9JOISv4mNTKSMbEiKp9mTLQonUSoaHwzThrirdfiGeZsgQ08sNe0aIxglZnwDvxsUbOrFni7NEMyffC/QsM/mMwMUkxAIytO5btPbZjb21PUERQmq/3qmqBWQgPF8YfRYvC6NGio/nu48nMihVNknwBjtw8QtMVTTlz9wyTGk4ickQkU5pOIY999nT2krweuQOWmB6jUYxYBATAzp2i67NLF/D2hsaNzR3dK4lOjCaPfR4cbBz489KfdKrQidGeo6lQoEKGrvu6XdvrRk7CwmDmTGGAk1BiFTlaBZAzbzhnhkeQ19GRkdG7yL14JYydTi01haAazng36Itn3a7iYll8HqiqKvsj9+Mf7M9v539Dq9HyYZkP+aDUB9QvLkrftd1qE3Q1KM3XzizTiXRx6RKULSvOdzt3FhWG6tVNcum45DgWHFtAXHIc473HU6NwDTZ03kCzUs1McuwhMT8yAUtMS0wM1KolzBAswH8XRLLYFbEL3/2+nLhzgvBh4dhZ23F64Gm0GtPs0F93xvuikZNRjctifaMITcbA33+DTYlg8g79ggS7vcQD1kZrDkbupJlHK3Lb5BQzpVWq4DlqNp4tWoikYAYSdAk0XN6QA9cPkNc+L2PrjWXQu4Mo6Fjwue/za+xHDpscab6+WfWm9XoIDBRdzZ9/DqVKie7mNm3E77sJuBl3kxkHZzDvyDweJj/kg1IfoKoqiqLQplwbk6whyR7IBCzJODdvCjGBDh1EE1XDhkI2skOHbO2/azAaWB+2Ht/9vhy9dZSCOQoytOZQDEaRKE2VfCF1u7YnifjhQ1iyBD7vDJcvi/v6kB9OMlNfm1hrexS9goqKajRw/JsBNFvzoShzhoWZTZozPiWe4MhgGpdsjIONAxXzV6R7pe70fqc3ObQvTrLpdX0yixjIk/PdGTOEZGTFimJUztpazE6biGXHl9F/c38MqoH2Hu0Z7TWa94q8Z7LrS7IXMgFL0k9oqBixWL1azDI2bgy5c8Ps2eaOLFXsj9xP53WdKZ23NPM/nE+PKj2ws375aEhGOm9Ts2u7eBG+mhvCbyeCSLnkRamal+k88hb/6/81NjaVqX1qJfmPX+LDi9+SgorWqOJdoqGYJ3VwMEvyvfPoDrMPz2b24dnEJscSOSISF0cXFrRakKrXrzy1EqNqpHvl7qleM8vFQNavh9694dEjeP99cRbQooVJ5ndVVWXP1T3kc8hHxQIVqVm0JgOqD2BErRGUzFvy9ReQWDRSiEOSdk6eFGpJQUFi59W3r/AnLZm9bxgPEh8w5/AcjKqRb97/5mlX6RO5vlfx7zNcELuu1Hqjvuz1P7WtRO6YIkybBpuPh0DPhmCdBIr4d1mraC329dkn4gsMhPbt2eKRh+nvliSmeCe+bt01TYnHVOM7N2Jv8P3u71l2YhkphhRalW3FaK/R1Hatnab50yb/a8LtR7c5+dnJNMeQaaiqmN91dBQe0xcvwnffiR1vtWomWUJv1BMYFohfsB9Hbh6hV5VeLG2z1CTXlmQvpBCHJOPEx8ODB+DqKiQLIyLAz0+YI2TT2d0nXI25ytQDU1l4bCHxung6Vej09EytSckmqbpGRuZ04b+7NpccOaiur8L43nk4dQry54dKn83mlNWzXXKfsh+x6GwplLnzYNAgNhWtxo72X/N7ifeEO04SaRq/yej4jqqqPEp5RE7bnBhUAytPr6T3O70ZUWtEumU3vd29+Xrn10QlROHs4Jyua5iMJ+e7AQFw8KBoHFy5EkqXFqNFJmLp8aV8v/t7wmPCKZW3FHNbzKVXlV4mu77EcpAJWPJqrl8XIyzz50Pt2sKppWRJcThppiaftPDzkZ8Z9McgFEWha6WujPYcTaWCldJ8HVN03rapWoRahYowdy7MnQkH7kGZ2mH88DOM6unB5sut6bJ+NagqWqNC/6/WoUQYRBMb4LvzCjdKPa9YlZaHgPQ+RBiMBjae24hfsB+57XLzV/e/cMvtxu1Rt196vptanox17Y7YTfvy7TN0rQyxcCFMmCDOd0uVEscovUyXFO/F3yOfQz6sFCsuPbiEi6MLU5pMoVXZVq+tvkjeXLL/HVRiHo4dE/KQxYuLnW6DBkKn+QnZNPmqqsqu8F2E3QsDxKjLsJrDuDL0CsvaLEtX8oXUz+m+jFOnRKXe1SuE73dNxOXDn/Ga2ZoLjcsT6jQOOzvoWKEje429mLDdyI5fNHh+0B/OnYN584CMPwSk9fUJugTmHJ5D2Vll6fBrB6ISomhdtjVPjq0ymnwB3i38Lg42DumaB84wERHPVNiuXxe/65s2ifd84EBR6ckgF+9f5NMtn+I2zY0tF7YA8K33twR/HExbj7Yy+b7lyB2w5BlPbkYaDfz+O2zZIuzRhgx5qTtOdsFgNBAYFohvsC9Hbh6hf7X+zG85n4oFKjKl6ZTXvv51Z6Pp6bxVVdi2DaZMEWNE2lLBGHs2ACWZU0CuuFx8W2ccg64VFDsvNzc8G/TC07o4fPopOD9fks3o+E1aXz/n8Bx8/vahZpGaTG40mTbl2pg8YdhobKjjVofI2EiTXveVhISIMnNgoFCqatsWxo2Db7813RKRIfgF+7Hx3EZsNDb0rNyT8vmFCIq5/aAl2QfZhCURKlWLF8P06cKRqGNHYZsGQtUnm7Ps+DIm7JnA5ejLlM5bGh8vn9d2NP+T1DZYpbaBKSlJHB0GBMCZM+BSJJnBA61Jqj6ZHw+MRUXFCoVvNA0YPzdMjHH9+CN89ZVJ4kzv6y/cv0BASAANizekY4WORCdGc+bemTQ3VqWVZH1y5gtLGAzPzncPHBB9C59+Kh4uCxc27VJGA6VmluJh0kMGvjuQwe8NxsXRPONhEvMjm7AkL+bq1WcWabGx4ow3f37xtWyeeB8kPiCPXR4UReHMvTM4Ozjj19gvXWdqqT0b/adgxouIihIWgLNmCSXCijVi6DT9Z/bqplO++WxcHOsz5YgdKboktHqVJgt3QPlGwrC3adPXxpnR8ZuXvT5/vgjarhnMpnOb0Gq0FHcS1Y489nmo41YnVdfOCJmafA2GZ+NCX34pFMFmzRLnuyaSiEzSJ7Hi5ApWnFzBn93/xM7ajo2dN1Iyb8l0GU1I3h7kDvhtRVWFmMD582LHO2IEvJf9B/4jYiIICAlgUegi1ndazwelPiDFkIKNlU26d2nFv/idF/0rUIDwSS1e+/rLl8U49OLFkJgvhBLNf8O98jUOx24mLiWOxiUaM6FYb2rW60pIZAhBAUPwji+A56CJUKVKumI2Fb039mbZiWXktc/LwBpit/ZvxaqsoHtgd9yd3PmhwQ+mueCTh8uNG+H0abC3F2e+rq4mmd8FIVs698hcZhycwZ34O1R1qcqvHX+V87uS55A7YIkYsVi/XmSJwEDRYLJwoRCQd3U1d3SvJfRWKH7Bfqw9sxYrxYpulbtRIk8JIOOKVek9Wz14EPz9xdtpbQ1N+oawvUhDrhgSuRIFjYs3YrJVE6rO2QAh3SC0PJ7veOIZcDjLtZmfkKhL5H8n/0eXil3IaZuTVmVbUaNwDfq808ckTVXpJSohiuO3j2c8AR88KMrM69eLv3fqJKo79vbClchEXH5wmSrzqhCvi6dpyab4ePnQoHgD6cErSRMyAb/pxMQIF5yZMyEyUoxYhIeL3a+FmLAbjAZarW7Fw6SHjKg1gmG1hlE0l2l8VSFtDVZGo+hN8/cXjnO5cqt0GbOPhArzeKdoWbbuTgFAg0L9dUepunk7lCgh3v9SpcRFzHCTjkqIYvah2cw6PIuohCjsrO3oWaUn7TzaZXksL8Lb3Zsvd3zJ3fi7FMhRIH0XCQ0VOuS5cwvRjCFDTPpwGXorlNN3T9OjSg9K5CnBiFoj6FihI5ULVjbZGpK3C5mA32SuXRP2c/Hxwolo9mwhoZdNR4ieoDfq+fXMr6w4tYINnTeg1WhZ32k9ZfKVwcnO9KIfqTlbTU4WWgz+/mJKxdXNSB/fTZzO7csvtw7gfNOZluXqo9VoSTGkoE0x4K0vCusXQuvWJit7phWdQcewP4ex9PhSEvWJfFjmQ3y8fKjrVtcs8byMJ/PAe67uoUP5Dql7UVycEM2OjYWxY+Gdd4THdOvWkDOnSeJSVZVtl7fhF+zHjvAdFHIsxEcVP8JGY8OEBhNMsobk7UUm4DcJVRVm9+fOiQ5PNzfw8RE3pHfSJ3yflcSnxLM4dDEBBwKIiImgnHM5rj28Rqm8pTJdkP5lDVYPH8LPP8O0aXDrFpSuH0Kbz37nqH45S+IiKa4tzuzKX9J701WSl85iatfJRKYcoZKxFHd+7AtmMoiPiInA3ckdG40Nlx5cokvFLozyGvV0FCa7Ub1QdXLY5CAoIuj1CTgyUlQU5s8X/4MaNxa/+4oiPKZNxP5r+xn4x0BO3jlJ4ZyF8W3ky4DqA7DRZF+DEYllIRPwm0BKijBEmDoVjh8XZbe+fYXZ/fjx5o4uVVyJvsK7C97lQeIDarvWZvoH0/mwzIdYKebZrd+6Jaay5s4VG6x6TaNpPXEZy65/xZUYUWb+vkgPvlxzA+vtE9HbO7CuYmOSHxQht00JrpE2mUhTYFSN/H7hd/yC/Th44yARwyIolLMQf3b/02zvY2qx0dgwoPoA3J3cX/2NixaJh0tVFW5bI0ZAzZomiyMuOY64lDgK5yxMbrvcGFUjS1ovoWulriZ1x5JIQHZBWz5//QV9+oiMUb48DB8udgH2WeCNmkEuPbjEqTunaOvRFlVVGfHXCDqW70htt9pmi+n8eVFmXr5c9K01/yiS3B9MY9P1+STrkzGqRgyqAQ1WTNhu5MsrhWHoUD5IKMe55P8+zxZxsmf/Fw0yNeYkfRK/nPwF/xB/zkWdwy23GyNqjaB/tf5mbawyCQaDOHQvXlwYI4SFiebBoUOhWDGTLfNPD95mpZuxqv0qgKea4RJJepFd0G8aZ8+Kc9xy5Z7dmJYsgSZNzNZdmxYO3TiE735fAsMCcXZwpkWZFmg1WqZ9MC3T1nydiMaRIzBpkuho1mqhQ/9wEt/7ls1XV6KGq3xUui1NLhj5VNlMiiI6r717j4EeX4JWy/kvfn/humnRik4vN+NuMmDLACoXrMwv7X6hY/mOFlsmTdIn8SjlEc6qPSxdKmr/ly7BJ58ISU4PDyEtZiLO3juLf7A/K06ueOrBO7LWyKdfl8lXkpnIBGwpPNE1nDpV7Ho7dYI1a6BMGfjzT3NHlyqO3jzKqG2j2H11N7ltc/NFnS8Y8t6QTC/tvcwFSFUhd0wRJk6E7dvBoVwwdb7+mzEdmlCisBPvLQxkUKmujNhnoNikQEhMpPTAlgR95Im3uzeers+6yDMqE5kWrsZcZdqBadyOv82q9qsokacEoZ+EUqlAJYtOGEbViNtUNzoml2R2wHmIjhaz6WvWQDvTdWs/qfopisLS40tZfXq19OCVmAWZgC2BVavghx/EztfFRbi2fPKJuaNKFSmGFOKS48jnkA9FUbgSfYWAJgH0q9aPnLam6VR9Hf9WulJVuH/GmZ5LHXl0HQoWMtD0h8ls049jL0aObJnMjp47uJUyBMduk8SWuEcPGDECz/LledHwVnq0otPKP2ehFUWhW6VuGIwGNFYayx+FOX0aq/LlqV64OkEXj0D9+mKUyMvLZFWdJ3rhfsF+fF//ez4o9QFjao/h89qfm98KUfJWIhNwduXWLShQQIyvnD0LtrbiYLJTJ/FxNudh0kPmH53PtIPTaFSiEcvaLKNaoWqEDwvPcgeYJ2Vg1aAQf7YwsQdLorufE02++3QLmM9BjT9/RV98+v0phhSCIoLwrNMQxtkIZ5yCr1aHyqhM5OtYdnwZvTf1Jqc2J8NrDWdYzWG45s7+AiqvxGgU1ZuAANixAzZvxruYN19c+pO7y+emfx74XyToElh6fCkBIQFcjr5MyTwl0Rv1AORzyGeSNSSS9CCbsLIbx46JMvOaNfDrr2KEKCUFbGws4nz3RuwNph+czs9HfyY2OZaGxRsypvYYGpdsbLaYak0I4uJeZ2Ku3cOY/wCamHfJUzwfdu9O57pxGdULvEPbGBd+TPiLFEVFq7FhR9/dz5WYsxqdQceaM2twcXShUYlG3E+4z6LQRQyoPiBTZqGzFJ1OnO8GBIiRucKikY0BAziUcJGaC2uypsMaOlXolOGlVFWl+vzqhN4OpWaRmvh4+WSKq5NE8jJe1YSFqqpZ9qd69eqq5AUYDKq6YYOq1qunqqCqOXKo6pAhqnrlirkjSzPDtw5Xrb6zUj9a95F69OZRs8YSG6uqvr6q6pRPr1I0WGWsncp4VMZbq8XGBqhLgo+pO7/8SDXmdFRVUINbvqP+NL+HGnx1n/liTopVpwRPUV0DXFW+Re26vqvZYjE5ycnivzqdqhYrpqpVq6rqihXPPq+qqs6gUx1/clQ/2/JZupe5dP+SOubvMWqyXlx3Q9gGdU/EHtVoNGYkeokkXQBH1JfkRLkDNidGo+hmNhhEM5VeL+Tz+vUTdmnZHFVV2XN1D37Bfoz0HEmD4g24/eg2ibpEiucxn39wdLTQaZg2TXxcq9UpIt7rym396ceBQ/fyX/C/ThPF2a7BIM4ba7z4ITWrmHFwBt/s+oaHyQ95v9j7+Hj50Kx0s2w/w/tazp4Vu93t28Wcl62tOGJxcXlhVWfVqVWUcy5H1UJV07TM4RuH8Qv2Y33YejSKhl29dpl1pE0iATmGlP24elVkiN9+g5Mnwc5OdDa7uwtV/2yOwWhgw7kN+O735fDNw+R3yM/d+LsAmeJ7mlof3jt3RPV+zhyhUtiqFVi1+oSN1+djp9phjRWq0YjWAANdKooXLVtmVmnOsHthFHMqhoONAw42DjQp2YTRXqMzXfkr01FV2LlTjAxt3Sp+x3sKPrijAAAgAElEQVT3FrKotrZQqNBLX9qlUpc0LXU/4T7t17Z/2l3v4+XD0JpDKZzTtD6/Eompyf53+zcFVYWQEJEhAgPFk3+HDsIswcXlmVC/BeC9zJt91/ZRKm8p5raYS68qvbC3yRzhj5eNEMGzxqcbN8DPD+ZuDiGl6E4q9bBicb+R1Khqy5JD1agR3YrPlp3hfMxlgt7JjXe9nni+01IsYIbkq6oq+67twzfYly0XtjCr2SwGvTeIftX60a9avyyPJ1MICYFGjUTz2oQJQr3KOXWdxjqDjm+3/cKW0GRiY4u+8KErxZDCyTsnqVG4Bnnt8+KodWRKkyn0r9Y/y7rrJZKMIkvQWcWRI/Duu6K0PGAADB5sETaAIJx0lh5fyvBaw7G2suaXk79gZ22XJc0stSftfOF8bREne1Z2acDkyUKdUF80CHo0wajoABhfbxzf1v8eHjwQ73P58jBqlHjoMVOVQVXVp2MwB28cxNnBmcHvDmbQe4MsfwzmwQMhlGEwwLhx4oFz/Xpo2TLNXfvrj16l02YPHAze5NMNAsRI18R2lajv4cjPR39m+sHpxCXHcX3kdXLZ5sqMn0giMQmyBG0OHjwQNoDJyfDNN1C9Ovzyi+hqzmEZ8oCXH1wmICSAJceXkKhPpHqh6tQvXp9ulbtlWQwvUpLSRTtwcmtJSo0DNDqqDvHlRO6fSDGK5KuooJu7CLy/g7x54dQpoRhmpi5yo2rESrFCURSmhEzhXsI9ZjefTe93euNg42CWmEzGpUuiqrN0KSQkCMGMJ8YIHVLpavQvAv6+jK2xAklWp55+7pEuikG/jyTuz63EpcTRsHhDPq/9OTm1crcrsVwylIAVRRkB9ANU4BTQR1XVJFMEZrGcPy9U/JctEzekVq2e3ZC6djV3dKkiJimG/pv7ExgWiLWVNd0rdTebk84/FaZ093Pw8EAp4s8UBrs4Bn8KPj7W1F+/mlx38vFQm4gRFa0RtNHu/HbgMq08Swk/XjNwP+E+cw7PYWHoQg73P0yBHAVY12kdBXMUfDPGYGbOhGHDxIhct27CGKFSpQxf9mZMIrbWlUi0OYKeKKxxxqA84JZ+HR95dMLHy4dqhaqZ4AeQSMxLug/AFEUpAgwFaqiqWhHQAB+ZKjCLZMYMoc+8eDF89JFosNq0ySLmd42qkYv3hRhFLttc3Ii9wedenxMxLIJFrRdlKPluDL1B7Uk7Kf7F79SetJONoTdS/VqfpmWxepiL28EGbp4LJt5mO9p+7XH4qhjf+0bj5qbQ8ciH3JscycbVOWgS/g5FE75lSfWvmLz7Wrpjzgjh0eEM+WMIbtPc+CboGyoWqEhsciwAhXMWttzkq9eL+fTTj7vJ69eHr78WTYWLF5sk+YJ46LIziGvdsf0KAK1akhq2a1jVfpVMvpI3hoyWoK0Be0VRdIADcDPjIVkQiYmirFy9OlStCg0bwrffwmefCRUrCyBZn8zKUyvxD/HnVtwtro24hqPWkf1995tEVzg1TVQv4+xZWONbhPD9EdCrE2iSQQGsbPmEdzGuXgV9BrKh0Hvomo9gi0ddkq2f6UpnhRHCv7n96DZlZpVBQaFb5W6M9hxNhQIVsjwOkxIbKxyIpk+Ha9eE49bUqVCxovhjYnyaluWLwETi9I2wUYsC4gz46w/MJ4wikWQG6U7AqqreUBTFH7gGJALbVFXd9u/vUxRlADAAwM3NLb3LZS9u3RKzLvPmQVSUML2vWhUqVBB/LICHSQ+fNrPcjLtJ5YKVmdFsBrYa0TBjKlH/f+swAyTqDPj9df6lCfjMGdE4u3atOC5/d/h6Dlsni7hUGLPPyPfb9kHfMtAHnJ1zs75Sw/9cJzOMEP6Nqqpsu7yNQzcOMe79cbg4urCg5QIal2hMkVxZ4wOcqXz/vRglio2F99+HWbOgRYtMXfKZrOeXmSLrKZFkF9LdBa0oSh5gPdAZiAF+BdapqrriZa95I7qghw8XyVevhw8/FOde3t4WUWaGZ/6mB68fpNaiWjQq0QgfLx8al2icKU46xb/4nRf9hilA+KTnb+SnT4v7/a8bE7F9bynVGlxj87CJXEgIwXtxXQwGA1oj7LjbDM9BE6FKFeC/u2x41jWbWTdtnUHH6tOr8Q/x5+Sdk7jmcuXsoLM4ah0zZb0s5fhx8d4qCnz1FUREZAuhEonEEsmsLuhGQLiqqvceLxIIeAEvTcAWicEgbACbNhUzo/nyCSeioUOhdGlzR5dqTt45iV+wHzm1OZnTYg41i9YkbFAY5ZzLZeq6qbHpO3MGhvqGsDPyD6zz3Mbhq00kKPewcihH7pTBeLp6ElTej6ATm/BuPxLPqq2eu1ZmGyH8m/3X9tNlfRciYyOpkL8CS1svpUulLpluq5ipGI2webPY7e7dK0wSmjaFH3+0mIdLicTSyEgCvgbUUhTFAVGCbghY+Pb2H8TGCpP7GTPgypVnN6Rx48wdWapRVZUd4TvwC/Zj2+Vt5LDJweD3Bj/9emYnX3i1TV9YGHz3HawJDoZe9aF4CnoFKuhdmbHdhboHzqGov8Dnn+PZcQSeHUe8dJ02VYtkaony9qPb3E+4T4UCFSiVtxTlnMsxt8Vcy5eKTEkRg9RTp8LFi+DmJmQjPR+ft8rkK5FkGhk5Az6oKMo64BigB0KB+aYKzGzExYm53UWLxMdeXjB5smiwsjB+2vsTY3eNxcXRhZ8a/MSnNT4lj32eNF0jtTKQL+NFu9OuZSvwq39BVu48gdbWgEO9QBI0elDAygidgyKpkFAdZd0saNMmTfGamvNR5/EP9mf5yeXULFKTPX32UNCxINt6/KfdwbLQ6cT4EIhdbuHCsHo1tG9vEXKoEsmbgFTCAjGne/MmFCkiznbLloVatcR577vvmju6VBOXHMfCYwvxcvWiZtGaXIm+wq7wXXSv3B1b67R7CJv6bPXiRfh+gsovwTtR6vhiLL6NgsnVUfiIO9qvUdQUrI1W1Lr9GfqS7dj/RYM0r2Eqjtw8wo97f2TTuU3YWtvS550+jPQcSam8liMZ+kJOnxa73d27RZu5VvtKYwSJRJIxpBLWy0hJEXON06aJBHz1qrghnT1rEab3T7gVd4sZB2cw7+g8YpJi+LLOl9QsWpMSeUpQIk/6RSjS08H8IsLDYfDkELbenYda9AD0uEB+TR5GXC7FgNVH6dC9P+T/kWTlJLZqZa7m80AxwwiRUTViVI1YW1kTEhnCnqt7GFdvHIPeG2Qyc3izoKrCiWjKFGH6YW8vjBESEsTv+yuMESQSSebxdibgqCgxQjR7Nty+LcQzvv1W3KjAopLv539/zvSD09Eb9bTzaMdoz9HULFrTJNd+2Rxtaudrr1+H8T8ksHTbUYxdm4JLEoqiMuZsPsYH3seugD1zvPsRlcMJW2NRbPF4+tqsGCF6QrI+mRUnV+Af4s/IWiPpX70//ar1o2/VvuTQWoZs6CvZsweaNBG73B9+EMYI+fKZOyqJ5K3n7UrABgNoNHDihGimatpUNFo1aWJWS7q0oKoq+yP341nUE42VBmcHZ/pV7cdIz5GUzFvSpGulpoP5Rdy+Dd9Mus/iU7MxVJ9Jyc7eRGhTMKgqVgbIpdpgt3g5dO5M4TP30AWeghc0aWU2MUkx/HxEzELfenSLqi5Vn87uZpa7U5YQHS0eMK2sYMwYqFdPDFW3amVRD5cSyZvOm5+AjUbRwTxtmlDtCQiABg2EZnOZMuaOLtXojXoCwwLxD/bn8M3DrO+0nnYe7fi89ueZtuarOphfxP378LVfBAvPBGCovAjqJdDsUTE6XD3N4ApaUgwpaLU2ePuvBzcvIOtHiP5J2zVtCYoIonGJxixvu5yGxRtmyix0lnH5svg9X7xYlJc7dhSfV5RnH0skkmzDm5uAHz2C5cuFfN6FC6LBquVjD1hFsZjkm2JI4ecjPzP1wFTCY8Ipnbc081rMo1mpZpm+dmqT48OHMCIghF/2B5HisQRNtXC63y3El+sTKP/wFvTogcdHcwm6GYK3uzeerp7/WScrEu6Zu2eYdmAakxpNIp9DPn5q8BN21nZULVQ109fOdKZPF6Iw1tbC9GPkSKhc2dxRSSSSV/DmdkF/8gnMny+6mEeMENZoT8YuLIAUQwpajRaD0UDZWWUpkKMAPl4+tCrbKtuI+cfHqwybvoOl4RMwFDwM1iloFSvWrtLROiofDBwIgwYJU3Yzoaoqe6/txXe/L79f/B0HGwcCOwXStFRTs8VkEgwG2LBBVHXKlRPqVWvXCp/pwoXNHZ1EInnMq7qg34wErKpw4IAov40ZA9WqiV3vvXtijteCyornos4REBLA1ktbOT/4PA42DkQlRGUrw/aEJD2D5qxjxRVf9PlDsdVp0dnoMKKiUTRMcGzJl5/+Ag7m9bpN1CXSYHkDDlw/gLODM0PfG8rAdweSz8GCG5AePRJ9C1OnivbykSNFd7NEIsmWvLljSDodrFsnEu+hQ+DkJIQEqlUTJWYLKTOrqsq+a/vwD/Hnt/O/YWdtR68qvUjUJeJg42D25BsSGUJQRBB13bw5uaMcw85WR58znMI2eRi/1Zayt5Jp1ltDiga0Gi3eHT83W/JN0iex79o+GpVohL2NPVVdqtKzck96vdMLBxvzPhBkmB9/BH9/iIkRSlX+/tC6tbmjkkgk6cRyE7DRKByIzpwRiXb2bOjVS9jnWBjHbh2j3tJ65LPPxzf1vslWc6chkSE0WN6AZH0K6G1Rl+6gQU0XBp+5SquLD9F07ETQkB6UunqZyIQjuDrU4E6UG7hmbZzRidHMPTKX6QenE5UQRcSwCFxzuzKnxZysDcTUhIWJErOiiC63Ro1g1CghFCORSCwayy5BL1wozrs++MBixogAEnQJLD2+lOjEaL6u9zWqqrL2zFpalm2ZrXZpVx6E8+GSroTFHRD2RUYN3QpN4H+FKqDs2Q3DhrHxvibLnYj+yd34u0zeN5n5x+bzKOURTUs2ZUztMXi7e1tuR7OqCgOQKVPg779hxw7Rua+qFnWcIpFI3oYzYAvhbvxdZh2axZzDc7ifeJ/67vXZ3nN7thPzvxJ9hU9+/YrtN39FMYKVYgRAa2XDjr67n+tirj1p5wtnhYs42WeqlGSyPhlba1tuxN6gzKwyT0VIqrhUybQ1Mx2dDlasEKNyp08LhaqhQ0VDYZ60aXhLJJLswZt7BmxBrDi5gn6/9SPFkEKrsq3w8fLBy9UrW+zSQiJD2BWxixqFa+Ca0oQR38E+900MOabliwNJhNcsx542VfBuMeg/I0QZVctKC6qqsufqHnyDfUnWJ7O953aK5CrCjZE3cLJzMvl6WcYTgRijEb78EgoUgKVLoUsXIRUpkUjeSGQCziSeKFblsctDhQIVqFG4Br2q9GKk50jKOme+ylNq+Wl7IGP3d0ZV9aBaoSzdi2O0FwddG+HhboXNllEUrluX2i95UEivWlZaMBgNbDq/Cd/9vhy8cZD8DvkZWnMoRtWIlWJlucn38mXRzbxrl1Bns7WFgweFJWA2eDCTSCSZS/aqfb4BGIwG1p1dh+ciT+ouqYtvsC8gvHd/bvlztkm+8Snx9Fv/PeP29kLlsRUgRgrWXMKsjbeofGIDNr9vEjKGr0gGPk3LYm/z/FyyqaUkFxxbQPu17YlKiGJui7lcHX6VsfXGZrvSfaoJDhbd+qVLi1n1994T1pcAxYrJ5CuRvCXIHbAJWXp8KRP2TOBK9BVK5inJ7Oaz6VWll7nDeiEBe6az6PR4St+15aozGBSwQoODqxs/HwqjZ4PUOeRkhpRkdGI0847Mo3S+0nQo34FulbqRzz4f7TzaZRsRknSzezd4e4sz3S++kMIZEslbjEzAGeRu/F2cHZyxUqy4eP8iBXIUwK+xH63Lts42ySIkMoT1YesJjw6nS4VuRIe0Y82EduxkPAXvFGNC4xpsr6xgzTsYbDzSfH5rKinJ67HXmRoy9WlH86B3B9GhfAdy2uakYwUL1TKOjxfCGUajaKiqWxcWLYJOncDR0dzRSSQSMyK7oNPJE8Wq5SeW82vHX2lZtiU6gw4bTfaSu1wSuoT+m/tjUA2gQqvjbvy26SqenuBSZA3HSuT4T8kzszuYX8SE3ROYsGcCRtXIRxU/wsfLx7I7mm/dgpkzhStRdDQ0bw6//27uqCQSSRYju6BNxBNdYf9gfzZf2IydtR293+lN+fzlAbJd8h2weQALji0AFVBAY4Rcj/KyaW0yLTvYsul4HcLS4HZkSp40qVUqUIncdrkp51yOT2t8ykjPkbg7uWf6+pnKvHlit6vXQ5s2MHq0kESVSCSSf2ChXSzmwaAa6LmhJyHXQxj//niuDb/GvA/nmdyHN73oDDpWnVpFfEo8ANXO2DHoINjpwcqgoNHYMWDhHFp1tEVRROl4YrtKFHGyR0HsfDNbQMOoGtl4biO1F9em7pK6LDy2EICOFToyo9kMy0y+qgrbt4uuZhBSqP37Cz3ywECZfCUSyQuRJehX8CjlEUtCl7Dq9Cp29tqJnbUdp+6comTektlGsSokMoRtl7cRkxTDhrBArsZeY6bTQE6HzWb1gjj62SwjaYwH+d8/RJPS/7UCzCpUVWVx6GL8gv04f/887k7ujPYcTZ+qfbLNe5lmdDpYs0ZoMp84AcOGCV1yiUQieYwsQaeRW3G3mHloJnOPzCUmKQYvVy9uP7qNu5M7lQpWMtk6G0NvZKh7eO+1vTRY1gC9UQ9AlXsaZv0ND648YhHw2aCcjBk3mPz5ARqaLO608MRWUVEUVp9Zjb2NPavar6JD+Q5YW1nwr9/MmTB5Mty4AeXLi8aqrl3NHZVEIrEgLPgOmDmcjzpP5XmV0Rl0tPVoyyjPUXi5mr6EuDH0xnMayjdiEvky8BTAa5Pwg8QH5LXPy76r+54mXysjVDlTjh8vLKBoB0/CJkKpUo+T/CLTjQillltxt5h+cDqLQhdxbMAxXHO7srbDWpzsnLKF+le6uHlTyEMqijBJKFsWFiyApk0tSotcIpFkD976BKyqKrsidnEl+gr9qvWjTL4yjH9/PJ0qdKJU3lKZtq7fX+efa34CSNQZ8Pvr/EsT5OEbh/EL9uP3c79xpe8J3i/mja2qRafqUA22hDosYH6I51OjnIwk+fRyPuo8/sH+LD+5HL1RT4fyHZ4+JOSxt1A946NHhTHC2rUQFAR16sCMGWD91v/zkUgkGeCtvYPoDDp+Pfsr/sH+hN4OpUSeEvR5pw8aKw1f1f0q09dPjYbyxtAbjPsjkPCkTSjWV3nEJXKnWDH0oJG7D39l7LGxJF8MwrlGEKM6eDNmgudzE0XpSfIZISohikpzK6Gx0vBx1Y8Z5Tkq2zSopRmjEf74QyTeoCDImRNGjIASJcTXZfKVSCQZ5K28i/x16S/6b+5PZGwk5ZzLsaDlArpX7p6lwhmv01DeGHqD4YFruaoZAxodqDD4IIwKL8Ofzl9RfUpncjvDrG89GTDAE5sXTEBltlGCqqpsvbSVvVf3MrHRRJwdnPml3S+87/5+tvEzTjfJydCnD9jbiyTcrx/kymXuqCQSyRvEW5OAr8deR2/U4+7kTpFcRSiZtyRzWsyheenmZtEU9mla9oU+uoMbFmXagWlM++NPHqoFAYPQaTbCdocOzLm0CpsIa0aOEcY5uXO/fI3MMkrQGXSsObMG3/2+nLp7CtdcroypMwYnOyfLVayKioK5c+Gvv4RcpL097NwJ5crxwqcbiUQiySBvfAI+fvs4U0KmsPr0ajqW78jK9iupWKAiu3rtMmtcT0rA32wNJDLhCC72ZahY/A6Dti0lWh9HvUi4n3MMD/NaoxoNGA1azh0aSY6ydzj9exHc3V+/xsuSfEaENo7cPEL7te259vAaFfJXYFmbZXxU8SO0Ggu1zbt4UTgSLV0KiYnQrBk8eAD580Ml03W8SyQSyb95YxPwjis7mLhvIjvCd+CodWTwu4MZVmuYucN6joLO17ikjiHZJpkYvZFzF6BdGPiczsW5gq2YlLMhVjsbYchzGJuEquRrYKCEx3nc3VN3fmsqo4SohChuxt2kcsHKlMpbCg9nD2Y3n2226oHJCA4WDVU2NtC9O4wcCRUqmDsqiUTylvBGJeBkffLTmdM/L/1JWFQYkxtNZkD1AdnOM/bQjUMsP7mcFEMKRowoKowIc2KK90+c9enFjGFWnN9lh3WeePKXccK+9B0ctI/waZq2XVlGjBIiYiKYEjyFRaGLKOdcjqMDjuJk58Sf3f9M1/XMjsEAmzZBbCz07g01a8LEidCrF7i4mDs6iUTylvFGKGE9SHzAvCPzmHloJktaL+GDUh8QmxyLnbVdtiqNPmla8t05gd23D1AzPg8ncycJsQrFmnWtdrB5Xm3mzxf9Pm37xhCWJ5TbjxKydIb3zN0z/LTvJ9acXoOVYkX3yt0Z7TX6qea1xREfL0rMU6cKuciaNSEkRPruSiSSTOeNVcK6En2FaQemsSh0EQm6BJqUbEJe+7wA5LLNHh2rIZEhBEUEYaVYseLQAk7HXaboQwg4qNDPvSGnPx3C9mv7uXPAmy71PImPh0GD4JtvwNnZCaifJXGqqopRNaKx0nDwxkF+O/8bw2sNZ3it4RTNVTRLYsgUVq0SnrsPHojEO2kStG0rk69EIjE7FpuAjaqRBssacDPuJl0rdWWk50gqF6xs7rCeY/uV7bRa1YoUQwpWKrg+MLD8gB0f1fkEm1UjUV3duBkIS3zqER4OH34Ifn6i8TarMBgNbDy3Ed9gX7pU7MLwWsPpXrk7bcu1tVzhjLAwMbdbtCi4ukK9es8ciWTilUgk2QSLTcBWihXL2y6nVN5SFM5Z2NzhPMftR7eZuT+AgAPTSEaPigqKho9dmtFj2yrInZvQUBjeA/bsEc22f/8NjRplXYzJ+mSWn1iOf4g/F+5foESeErg4inNQrUaL1j77lO5ThaqKN9PfH7ZsEbvemTNFk1WdOuaOTiKRSP5DhhKwoihOwEKgIsJ1tq+qqiGmCCw11CtWL92vzagRwou4cP8CU3b+wLIzK0nBgHc4BLtboddYodVoqd99HHeScvP1KFi8GPLlE9ax/fqBJus0QADoGtiVwLBAqhWqxpoOa2jv0T5LhUhMyvr1orR85IgYH/ruO/jsM3NHJZFIJK8kozvg6cCfqqp2UBRFC1iEr5wpNJI3ht54OsPr6lCD75u1Y/LmRoTqIulzHEbaN6D0oG8IKW5D0NXdeBXxZs9KT5r+CElJYuJl7FhwyqLm7Ftxt5hxcAbDag3DxdEFHy8fPqvxGQ2LN7RMc4SEBHB4/Ov255+is/nnn6FHDyGiIZFIJNmcdHdBK4qSCzgBlFBTeZHs4gdce9LOFypEFXGyZ/8XDV77+o2hNxgWuJprmi8APQo2uBkmM9HuHg1u3qDgsK+hTBlAVEY3bAAfH7hyBVq2FMqGpUub+qd6MRfvX8Qv2I9lJ5ahN+pZ3mY53Sp3y5rFM4ObN4URwrx5QrWqZk2RfB0dpSORRCLJdmRWF3QJ4B6wRFGUKsBRYJiqqvH/WnwAMADAzc0tA8uZjoxoJKcYUhizeSpRxrlgLVx+UHU8VI8zS9OLLrOfJfCTJ2H4cNi1S+g7bNsGjRub5Ed4LQajgW6B3Vh7Zi1ajZY+7/RhtNfoTHV4ylROnxZPLr/8IuZ527d/ps0sNZolEokFkpEEbA1UA4aoqnpQUZTpwBfAuH9+k6qq84H5IHbAGVjvKRk9v02vRnJc4kMq+BYjkoeUeASRuUFvpQA22BkrPU3g9++LMaJ580SJedYs+OSTzDfQUVWV47ePU7VQVTRWGnLZ5uKLOl8wrOYwCjoWzNzFM5OkJNHJnJwMn34qnmqeuBJJJBKJhZKRlHAduK6q6sHHf1+HSMCZiinOb1OjkfzkjPdqwl6cNQpT2gTQpmoRel/NQ75Lzhwr2J1tDjlI0pzGzlgJW6MHhXI6MGuWSL6xsTBwoOgHypvXhG/AC9Ab9aw/u57J+ycTejuUMwPPUD5/eea3nJ+5C2cWOh38+its3AirV4OdHQQGQuXKmf9mSiQSSRaR7gSsquptRVEiFUUpq6rqeaAhcNZ0ob0YU3jcvk4jeWPoDQYGzuK2xhfV2kgsMHZFUWAk388+y8ZzD5gXeApbnQFbvVCHMl7Pz42DVRlyCRo0gOnToWJFk/3YLyRJn8SS0CX4h/hzJfoKZfOVZWHLhZTMY6EevLGxsHAhTJsGkZFiIPrmTTHP6+1t7ugkEonEpGS0KDoE+OVxB/QVoE/GQ3o1pvK4fZlG8pWoi4xZ2YpbOc6JTyigGOGB5vrjJN/guQR+NQIS91XkwZkCuLuLiZjMFlpSVRVFUVh75CKD/xqBtbEYZbXf8ZNXH9pVc828hTOTEydEmTk2ViTbuXOFM5FsrJJIJG8oGUrAqqoeB17Y3ZVZZIbHrVE1cj/hPvlz5Mf+Xgxx6jk6nrFjg4cOvZUKijWq9bvPJfkmZYtwbEMRfJeKGd4ffoBRo0S1NLO4GXeTqSFTOXHnBJ9VWMzE32/gop+FtepCUrLC1xvOYKVYZYletEk4eRKuXhWt4RUqCEeiPn2gRpb+SkkkEolZsLjthU/TstjbPC8YkV6P290Ru+mwrAUlvs1Dh+9FvbiQx7s0T5jGkeJrcNZPwknfnYIpP2Jr9KCwkz2qCuvWgYcHTJgA7drB+fPw9deZl3zPR52n32/9KD69OAEHAnB2cGbyXydJ1BmwUQuhILbbT0rx2RpVFbJfTZtClSpiIFpVRYfa7Nky+UokkrcGi5OiNIXHbWxyLF+t/YTZl1eLErMCveOKoKakoGi1fNi/A/sDT2HUeWBr9ABEku9cqgKNGsHOnaIf6H//E1XTl2EKta0/L/1J81+aY2ttS7+q/RjlNYoSeUpQ/IvfX/j9aS3FZyk7d4qEe+KEsP/76SfRHm6JQiASiUSSQSwuAUPGPKeVeOgAABGySURBVG4Bli4czOyo1U//bmWlwbZLDxSt9un14VmSL2DnSMGL1Rk50ZFcuVI3VpTebm1VVfn7yt8k65NpWbYl7xd7n++8v+OTGp9QIEeBp9+XGaX4TCE2VnQ158sn/q7TCR3Orl3B1ta8sUkkEokZsbgSdHo4d/Mk/fzrsWLF5wD0bf8DizTtsbexR6No0Gq0eLt7P/eaNlWLsPfzBnxTtgURc99nyypHPv4YLlwQdoGvm+l9Vbf2i9Ab9aw5vYbq86vTdEVTfIN9AbC3sWfc++OeS75g2lJ8pnD9upD/cnUVB+QA9evDqVPinFcmX4lE8pZjkTvg1DJ/dwDTgiYRxj3sdFD66iPoDo4ubvQduw6Px1693u7eeLp6Pvfao0eFoc6BA+DpCX/8AdWrp37ttHRrbzq3iVHbRnE5+jJl85VlUatFdKv0arlIU5TiM4WTJ4Uj0apV4my3Uyehzwyi1CzLzRKJRAK8wQm41YTybDaGgQrWqsL6yhNo3vGr577H09XzP4n3wQPRUPXzz1CgACxdKvJHWqdhXlcijkmKQUEht11ujKqRfA758GvsR+tyrbFSUrdYRkvxJkNVnyXWKVOEaMbgwTBsGLi7mzU0iUQiya68MSXoZF0SS9Z8SXRUJAB57POgqIACqsaKE4WsXrn7MhqFBkSZMrBgAQwdKrqbe/VK3yjqy0rEH7/vxJi/x+A21Y2AkAAA2pRrw4GPD9DWo22qk2+2QKeDFSugalU4dkx8buJEIaIxdapMvhKJRPIKLHYHHPK4fFyjYFVCt/+P6ZG/ctNOh+76VQaMWsmnnf35dXlDUgwpLzzj/SdHjohz3UOHhHf77Nmiyzkj/LtEnCfXfZwL/8VnO9ahN+rpVKETbT3aAlieHWBsrHhKmTZNnPV6eIjPARQubN7YJBKJxEJItx1hejCVHWFIZAgNlzckSZ+I+niX2/iWA5979KNhv59QcuR4+n0vO+MFYZrw1VcilxQsCH5+0K1b5hxTtl3Tlq0Xtz51JSqZ10LlIg0GYYRw7ZpQrPLxgQ8+kIpVEolE8gIyy47QbARFBJFiSOFx7mWQcwtmjvvtuSQgZnATuRlTmS1Oifg0vfF0V2o0ikmYMWPg4UNhrvPtt6ZztVNVlR3hO/Dd78vs5rMpna80U5pMYW6Lubg4uphmkazk5ElhivDjj0L2a+JEUauXohkSiUSSbiwyAXu7e6PVaEV52VpL11Zf/yf5vmwGtxhFGDhQdDfXrSvKzZUqmSYug9HAhnMbmLRvEkdvHcXF0YUr0Vcona80JfJYmH2eqsKOHaKj+a+/wMFBHIiXLStmeCUSiUSSISwyAXu6erKj546XlpdfNIMb/0jhs0FG7h4EZ2dYvlxID5uq3GwwGqj6c1VO3T1F6bylmf/hfHpU6YGddSaKQ2cW4eFCY/P4cVGb//FH4cMrrQAlEonEZFhkAoYXjxA94Z+ztqoK8WeLEL2rHMYEWwYNFLoQTk4ZjyE2OZbfzv9G98rd0Vhp6FmlJ+5O7rQt1xaNleb1F8hOxMUJlZHq1aFIEaFctXChOBTPTIcJiUQieUux2AT8Kp7M4Kbcc+TB3xVJjsyHtlA05fqcYtbMdzN8/bvxd5l+YDqzD8/mYfJDqhSsQqWClRjtNdoE0WcxN2/CjBkwbx44OkJEBGi1sH27uSOTSCSSN5o3MgEPqVeOQT5JPDjgjpVWT96mJ8lf4ybftc/YYe/9hPuMDxrPotBFJOuTaefRjjG1x1CpoIkOkbOSixdFM9WKFaKzuV070dH8Oo1NiUQikZiEN+5uu3kzfDGkMA+uQoHqN7H1OoNrYQ0+TSulWzXqUcojHLWO2FnbERgWSNeKXfm89ueUdc4musupRVWFeIZWC5cvi87mAQNgxAgoaaFjURKJRGKhvDEJ+No1oXy4caPwdt+zB+rWLQykXxhi37V9TNo3icvRlzn92WlyaHNweehl7G2ymePQ69DrYf160dFcvz74+go/3sjIZy5FEolEIslSLF49QacTecXDQ0zLTJokVBHr1k3f9YyqkS0XtlBncR3qLqnLwRsH6VapGzqjDsCykm98PMycKWZ2P/pIqFU9mblSFJl8JRKJxIxY9A54/34xHXP6NLRsKXqJMio/vOXCFlqvbo1bbjdmNptJ36p9cbBxMEm8Wc6wYbBoEXh5QUAAtGolFaskEokkm2CRUpQAP/0kXItcXcUmr3Xr9F0nQZfAktAl2GhsGFB9AHqjnsCwQNqWa4uNxsYksWYZFy4IN6LBg8VO99w5Ye/k5WXuyCQSieSt5FVSlBa7HXoiQ3z2bPqSb3RiND/u+RH3ae4M3jqYPy7+AYC1lTWdKnSyrOS7fz+0aQPlysGyZc+cicqVk8lXIpFIsikWW4L28kp/bll6fClDtg7hUcojmpVqxpd1vqSOWx3TBpgVqKpopvr7b6FSNXas2P0WKGDuyCQSiUTyGiw2AaeVi/cv4qh1pFDOQpTIU4KWZVoypvYYqrhUMXdoaSMxETZtgs6dRSNVw4bibLdPH3jsAiWRSCSS7I/FngGnlmO3jjFp3yTWnV3H4PcGM6PZjCxd32RERcGcOTBrFty7B8HB4PliKU6JRCKRZA/eODvC1LA7Yjc/7fuJbZe3kcs2F2Nqj2FYrWHmDivtREfDuHHCPzExEZo3F4fftWqZOzKJRCKRZIA3KgGrqory2N5o6YmlHL99nIkNJ/JZjc/IbZfbzNGlkQcPxLmuvb2Q9+rcGUaPFiojEolEIrF43ogStM6gY+WplUzeP5nlbZdTo3ANohKiyGGTw7KEM4xG2LoV/Pzg6lWh12xtDcnJYGtr7ugkEolEkkbeyDEkEDO8Mw7OoOSMkvTe1Bub/7d39zFSVWccx78PK6SBUtQupRSwYDWYNiJs8IUXKbG1UUqkxWgkYE00ISoqJpoCmvgaNdK0ahsUbYulLbZaW6zBl2iw2pgIAREFgxWoGLZSoLRCsVWRffrHudTJeO8y7A5z7hl/n2Qzwz1ndp+HM/c+e8/cvaelJ//dF5YibO3dmk7x/eADePDB8Le7kyfDpk0wa1a4hSSo+IqINKFkp6D3d+xnxH0j2PyvzYw/ZjwLJy/k7OPO/v8UdFKefRYuvjgU4MWLw20je/WKHZWIiBxGyRbglh4t3Pj1Gxl21LD0/oa3vR3uvhtaW2Hu3HBh1fLlYaGEFH+BEBGRQ5ZsAQa48KQLY4dwaNatCytHPPRQuInGzJlhe48ecMYZcWMTEZGGSvoz4KTcfjuMGAGPPgqXXx4+57333thRiYhIJEmfAZfagTV429rg+OPhzDPDVc6XXaZlAEVEpPtnwGbWYmavmNmyegSUvOo1eBctCttPPjncq1nFV0REqM8U9GxgQx2+T/ruuAOOOQauugoGDoSlS+G222JHJSIiJdStAmxmg4FvAz+rTzgJevvtcEEVhKubTz8dXnzx4yUCe+hjdhER+aTuVoe7ge8DHUUdzGymma02s9U7d+7s5o8rkZUr4dxzYdiwUGwhTD0/9hiMGxc3NhERKb0uF2AzmwzscPeXO+vn7g+4+2h3H92/f/+u/rhy6OiAZctgwoSwGMJzz8G8eeEiK9DZroiI1Kw7V0GPA84xs0nAZ4DPmdmv3X1GfUIrofffD+vu9u4Nd90Fl1wCffvGjkpERBLU5QLs7vOAeQBmNhG4tumK7+7dcP/98MQT4Wy3d+/weMIJ0LNn7OhERCRhmjPN094e1twdMgTmzAnFdteu0HbiiSq+IiLSbXW5EYe7Pw88X4/vFd2qVTB2bLiy+fzzwxq8bW2xoxIRkSajO2G5wwsvwPbtYdH7tja47rrwWe/QobGjExGRJvXpnYLevz/cl/nUU8MqRLfeGopxSwvcfLOKr4iIHFafzgL89NMwfDicdx68+y4sXBimnrUUoIiINMinZwp6166wQMKAAdCnT1iLd/58mDIlnPWKiIg0UPOfAb/1Flx5ZbhH8y23hG3jx8NLL8HUqSq+IiISRfOeAb/ySjjDfeSRUGSnTw/r8IKmmkVEJLrmKsDuHxfXBQvCDTSuuQZmz4ZBg+LGJiIiUqE5pqD37YMlS2DUqLBIAoRlALduDWfBKr4iIlIyaRfgvXvhnnvguONgxgz48EN4773QNmAA9OsXNz4REZEC6U5Bd3TAyJGweXNYg3fBApg0SSsSiYhIEtItwD16hJtnDB0KY8bEjkZEROSQpFuAAaZNix2BiIhIl2i+VkREJAIVYBERkQhUgEVERCJQARYREYlABVhERCQCFWAREZEIVIBFREQiUAEWERGJQAVYREQkAhVgERGRCFSARUREIlABFhERiUAFWEREJAJz98b9MLOdwNt1/JatwD/q+P1iUi7l0yx5gHIpq2bJpVnygPrn8mV375/X0NACXG9mttrdR8eOox6US/k0Sx6gXMqqWXJpljygsbloClpERCQCFWAREZEIUi/AD8QOoI6US/k0Sx6gXMqqWXJpljyggbkk/RmwiIhIqlI/AxYREUmSCrCIiEgESRRgMzvLzP5iZpvMbG5Ou5nZj7P218ysLUacB2NmQ8zsT2a2wcxeN7PZOX0mmtluM1ubfd0QI9ZamNkWM1uXxbk6p73042Jmwyv+r9ea2R4zu7qqT2nHxMwWmdkOM1tfse1oM3vWzDZmj0cVvLbT/arRCnL5gZm9kb1/lprZkQWv7fS92GgFudxkZn+reB9NKnhtacalII+HK3LYYmZrC15btjHJPf5G3V/cvdRfQAuwGTgW6AW8Cny1qs8k4CnAgNOAlbHjLshlINCWPe8LvJmTy0RgWexYa8xnC9DaSXsS41IRbwvwd8IfzicxJsAEoA1YX7FtPjA3ez4XuLMg1073q5Lk8i3giOz5nXm5ZG2dvhdLkstNwLUHeV2pxiUvj6r2HwI3JDImucffmPtLCmfApwCb3P2v7v4h8FtgSlWfKcAvPVgBHGlmAxsd6MG4+zZ3X5M9/zewARgUN6rDKolxqfANYLO71/NubYeVu/8Z+GfV5inA4uz5YuA7OS+tZb9qqLxc3P0Zd/8o++cKYHDDA+uCgnGpRanGpbM8zMyA84HfNDSoLurk+Bttf0mhAA8Ctlb8u51PFq1a+pSKmQ0FRgErc5rHmNmrZvaUmX2toYEdGgeeMbOXzWxmTntq43IBxQeTVMYEYIC7b4Nw0AG+kNMntbEBuJgwo5LnYO/Fsrgim05fVDDVmdK4nA5sd/eNBe2lHZOq42+0/SWFAmw526r/dqqWPqVhZp8Ffg9c7e57qprXEKZATwJ+AjzW6PgOwTh3bwPOBmaZ2YSq9mTGxcx6AecAv8tpTmlMapXM2ACY2fXAR8CSgi4Hey+WwX3AV4CRwDbC9G21lMZlGp2f/ZZyTA5y/C18Wc62bo9LCgW4HRhS8e/BwDtd6FMKZtaTMPhL3P0P1e3uvsfd92bPnwR6mllrg8Osibu/kz3uAJYSpmkqJTMuhIPEGnffXt2Q0phkth+Y6s8ed+T0SWZszOwiYDIw3bMP5KrV8F6Mzt23u/t+d+8Afkp+jEmMi5kdAUwFHi7qU8YxKTj+RttfUijAq4DjzWxYdpZyAfB4VZ/Hge9lV92eBuw+MKVQJtlnJj8HNrj7jwr6fDHrh5mdQhijXY2LsjZm1sfM+h54TrhYZn1VtyTGJVP423wqY1LhceCi7PlFwB9z+tSyX0VnZmcBc4Bz3P0/BX1qeS9GV3X9w3fJjzGJcQG+Cbzh7u15jWUck06Ov/H2l9hXptXyRbia9k3CVWjXZ9suBS7NnhuwIGtfB4yOHXNBHuMJ0xavAWuzr0lVuVwBvE64ym4FMDZ23AW5HJvF+GoWb8rj0ptQUPtVbEtiTAi/NGwD9hF+S78E+DywHNiYPR6d9f0S8GTFaz+xX5Uwl02Ez94O7C8Lq3Mpei+WMJdfZfvBa4SD98Cyj0teHtn2XxzYPyr6ln1Mio6/0fYX3YpSREQkghSmoEVERJqOCrCIiEgEKsAiIiIRqACLiIhEoAIsIiISgQqwiIhIBCrAIiIiEfwPpW2lSeUf0RUAAAAASUVORK5CYII=\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:                Fri, 10 Jul 2020   Prob (F-statistic):           1.66e-29\n",
      "Time:                        05:46:36   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.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
