{
 "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:                Mon, 24 Feb 2020   Prob (F-statistic):           5.44e-29\n",
      "Time:                        22:49:06   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1xV9f/A8de5l+lAxC3uvcVEE81EcWXuHKm5ysz1dZv17furbLhwZZlpZrkqt+VKza3gwAm5ByLgQBQVmffe8/vjI5jlQLhwL/p+Ph48wnPPPefD1Xifz3q/NV3XEUIIIUTWMti6AUIIIcSLSAKwEEIIYQMSgIUQQggbkAAshBBC2IAEYCGEEMIGJAALIYQQNvDUAKxp2nxN065rmhbyt2Nemqbt0zTtqKZpQZqm1c3cZgohhBDPl7T0gH8CWv7j2GRgnK7rXsDH9/8shBBCiDRyeNoJuq7v0jSt1D8PA273v88DRKblZvnz59dLlfrnpYQQQojn06FDh27oul7gUa89NQA/xnBgk6ZpU1C96PqPO1HTtP5Af4ASJUoQFBSUzlsKIYQQ2YumaZce91p6F2ENBEboul4cGAH88LgTdV2fq+u6t67r3gUKPPIhQAghhHjhpDcA9wZW3f9+OSCLsIQQQohnkN4AHAk0uv99E+CsdZojhBBCvBieOgesadovgC+QX9O0cOAT4F3gK03THIAE7s/xpkdycjLh4eEkJCSk9xLZhouLC8WKFcPR0dHWTRFCCGFjaVkF3e0xL9W2RgPCw8PJnTs3pUqVQtM0a1zSLum6TnR0NOHh4ZQuXdrWzRFCCGFjNs+ElZCQQL58+Z7r4AugaRr58uV7IXr6Qgghns7mARh47oNvihfl5xRCCPF0dhGAhRBCiBdNtgvAa45E0GDiNkp/sJ4GE7ex5khEhq9pNBrx8vKiatWq1KxZk2nTpmGxWJ74ntDQUH7++ecM31sIIcSLKb2ZsGxizZEIPlwVTHyyGYCImHg+XBUMQPtanum+rqurK0ePHgXg+vXrdO/endu3bzNu3LjHviclAHfv3j3d9xVCCJE11hyJwH/TaSJj4inq7sqYFhUzFDesIVv1gP03nU4Nvinik834bzpttXsULFiQuXPn8s0336DrOqGhoTRs2JCXXnqJl156iYCAAAA++OADdu/ejZeXF9OnT3/seUIIIWwrpfMWEROPzoPOmzVGUDMiW/WAI2Pin+l4epUpUwaLxcL169cpWLAgW7ZswcXFhbNnz9KtWzeCgoKYOHEiU6ZMYd26dQDExcU98jwhhBC29aTOmy17wdkqABd1dyXiEcG2qLur1e+l6zqgEoUMGTKEo0ePYjQaOXPmzCPPT+t5QgghslZWdd6eVbYagh7ToiKujsaHjrk6GhnToqJV73PhwgWMRiMFCxZk+vTpFCpUiGPHjhEUFERSUtIj35PW84QQQmStx3XSMqPz9iyyVQBuX8uTCR2r4+nuigZ4ursyoWN1qw4hREVFMWDAAIYMGYKmady+fZsiRYpgMBhYtGgRZrMaxsidOzd3795Nfd/jzhNCCGFbWdV5e1bZaggaVBC29ph9fHw8Xl5eJCcn4+DgQM+ePRk5ciQAgwYN4o033mD58uU0btyYnDlzAlCjRg0cHByoWbMmffr0eex5QgghbCslZtjbKmgtZa4zK3h7e+v/XJh08uRJKleunGVtsLUX7ecVQogXmaZph3Rd937Ua9lqCFoIIYR4XkgAFkIIIWxAArAQQghhAxKAhRBCCBuQACyEEEIAZ6PPZun9JAALIYR4oR29epRGPzWi0qxKnLt5Lsvum+32AVtbdHQ0fn5+AFy9ehWj0UiBAgUAOHDgAE5OTrZsnhBCiExgspi4FX+LAjkL4OrgyqWYS0xrPo0iuYpkWRte+ACcL1++1FKEn376Kbly5WL06NEPnaPrOrquYzDIgIEQQmRn8cnx/Hj0R/wD/KlVuBaruq6iYv6KnB96HqPB+PQLWJFdBeDhw+F+LLQaLy+YMePZ33fu3Dnat2/PK6+8wv79+1mzZg01a9YkJiYGgF9//ZU///yTefPmce3aNQYOHEhYWBgGg4GZM2dSr1496/4gQggh0u1W/C1mB81mxr4ZRMVFUa9YPfp49Ul9PauDL9hZALY3J06c4Mcff+S7777DZDI99ryhQ4fy/vvvU69ePUJDQ2ndujUhISFZ2FIhhBBPMmPfDD7b9Rkty7Xkw1c+pGGJhmiaZtM22VUATk9PNTOVLVuWOnXqPPW8P//8k9OnT6f++datW8THx+PqattKG0II8aI6feM0/gH+tK/UntYVWvOfl/9Dx8odqVm4pq2blsquArC9+XtBBYPBwN/zZickJKR+r+u6LNgSQgg7cCDiAJP2TmL1ydU4OzhTvWB1APLnyE/+HPlt3LqHyaqiNDIYDOTNm5ezZ89isVhYvXp16mtNmzZl1qxZqX8+au2JbCGEEE/V97e+vDzvZbZd3MZ/G/6XS8MvMazeMFs367EkAD+DSZMm0bJlS/z8/ChWrFjq8VmzZrF3715q1KhBlSpV+P77723YSiGEeDGYLWZWnFhBgkmNSPqV9mNKsymEDQ/jiyZfUDBnQRu38MmkHGEWe9F+XiGEsLYEUwILji7AP8Cf87fOs7D9QnrW7GnrZj3Sk8oRyhywEEKIbMFkMTElYAoz9s3g2r1r1PWsi38zf9pVamfrpqWLBGAhhBB2LS45jhyOOTBqRlaeXIlXYS/GNhiLbylfm28lyggJwEIIIezS2eizTAmYwoqTKzgz5Az5cuRjR+8d5HTK+fQ3ZwMSgIUQQtiVoMggJu2dxMoTK3EyOtHXqy8mi0qG9LwEX5AALIQQwo5cvHWROt/Xwc3ZjQ9e+YChLw+lcK7Ctm5WppAALIQQwmbMFjOrTq4i5HoI4xqPo3Te0izrtIzmZZuTxyWPrZuXqZ66D1jTtPmapl3XNC3kH8f/o2naaU3T/tI0bXLmNTHzhYeH065dO8qXL0/ZsmUZNmwYSUlJ7Nixg9atW//r/HXr1lGrVi1q1qxJlSpVmDNnjg1aLYQQ2VeCKYG5h+ZSaVYluqzowrITy1L383au2vm5D76QtkQcPwEt/35A07TGQDughq7rVYEp1m9a1tB1nY4dO9K+fXvOnj3LmTNniI2N5aOPPnrk+cnJyfTv35+1a9dy7Ngxjhw5gq+vb9Y2WgghsrEdoTso/VVp3lv3Hu4u7izvvJyQgSG4OLjYumlZ6qlD0Lqu79I0rdQ/Dg8EJuq6nnj/nOtWa9GjglmXLjBoEMTFQatW/369Tx/1deMGdOr08Gs7djzxdtu2bcPFxYW+ffsCYDQamT59OqVLl6Zx48b/Ov/u3buYTCby5csHgLOzMxUrVnz6zyWEEC+wq7FXuRV/i8oFKlMhXwVqFa7FKJ9RNCndJFtvJcqI9KairAA01DRtv6ZpOzVNe2zJIE3T+muaFqRpWlBUVFQ6b5d5/vrrL2rXrv3QMTc3N0qUKMG5c+f+db6Hhwdt27alZMmSdOvWjSVLlmCxWLKquUIIka2cu3mOAesGUGpGKQZvGAxA0dxF2dBjA35l/Owj+GZhRsi/S+8iLAcgL1APqAMs0zStjP6IvJa6rs8F5oJKRfnUKz+px5ojx5Nfz5//qT3eR7Tvkf8AHnccYN68eQQHB/Pnn38yZcoUtmzZwk8//fRM9xVCiOfZsavHGL9nPCtOrMDR4Ejvmr0Z02CMrZv1sMREWLxY1cL9/XcoXTpLb5/eHnA4sEpXDgAWwL7qPKVR1apV+Wd+6jt37nD58mXKli372PdVr16dESNGsGXLFlauXJnZzRRCCLun6zoWXY0Ibr24lT/O/cH79d8ndHgoc9rMoZxHORu38L6YGJg0SQXcfv3AwUFNYWax9AbgNUATAE3TKgBOQNa33gr8/PyIi4tj4cKFAJjNZkaNGkWfPn3IkSPHv86PjY1lx9962UePHqVkyZJZ1VwhhLA7KVWJ6s6ry8Jj6nfpAO8BhA0PY0LTCfa1jzcuDsqVgw8+gGrVYPNmOHwY6jx2JjXTpGUb0i9AIFBR07RwTdPeAeYDZe5vTfoV6P2o4efsQNM0Vq9ezfLlyylfvjwVKlTAxcWF8ePHA7B161aKFSuW+nXkyBEmT55MxYoV8fLy4pNPPpHhZyHEC+nvW4k6L+9MTEIMbs5uAORwzGE/W4n++kv1eEFNZY4fr4Lu5s3QrBnYaB5ayhFmsRft5xVCPL9aLG7B5vOb8S7qzdgGY+lQqQNGg9HWzVJ0HXbtAn9/WL9eBd4zZ8DTM0ub8aRyhOkdghZCCPGCuRp7lY+2fkRMQgwAYxuM5c+ef3Kg3wE6VelkP8H3zBmoV09taz1wAD77DMLCsjz4Po2kohRCCPFE526eY0rAFH46+hPJlmRqF61Nx8odaVK6ia2b9kB8PFy+DBUqQJEiYLHA7NnQuze4utq6dY8kAVgIIcQjJZmT6Lm6Z+pWoj5efRhdf7T9rGYGiI6Gb7+Fr7+GggXh+HHInRsOHrR1y55KArAQQohUuq5zIuoEVQtWxcnohK7rvF//fYbVG2Zfq5lDQ2HqVJg/X61sfv11GDPGZguq0kMCsBBCCMwWMytPrmTS3kkcv3ac80PPUyJPCZZ1Xmbrpj3MYgGDAfbsgTlzoEcPGD0aqla1dcuemSzCEkKIF9jftxJ1XdGV2KRYZr8+m0I5C9m6aQ/oOmzaBE2bwpT7tX+6doWLF+HHH7Nl8AUJwIwYMYIZM2ak/rlFixb069cv9c+jRo1i2rRpVKtW7V/v3bdvHy+//DJeXl5UrlyZTz/9NCuaLIQQVnMt9hqDNwzG3cWdFZ1XcGLQCfq91A9nB2dbNw2Sk2HRIvDygpYt4eRJyJtXveboaHermp/VCx+A69evT0BAAAAWi4UbN27w119/pb4eEBBAgwYNHvne3r17M3fuXI4ePUpISAhdunTJkjYLIUR6Xbl7hbFbxtJ5eWcASrqX5PiA4xzod4A3qrxhP1uJQFW569ULTCbV0714Ed5919atshq7mwP2/cn3X8e6VO3CoDqDiEuOo9WSf5cj7OPVhz5efbgRd4NOyx4uR7ijz44n3q9BgwaMGDECUJWRqlWrxpUrV7h16xY5cuTg5MmT5E154vqH69evU6RIEUCVMaxSpUoafkIhhMh6Z6LP4L/Xn4XHF2KymOhStQtJ5iScjE5ULmAnyYGuXoWZM2HgQCheHIYNg+7d4bXX1Lzvc8buAnBWK1q0KA4ODoSFhREQEICPjw8REREEBgaSJ08eatSogZOT0yPfO2LECCpWrIivry8tW7akd+/euLi8WAWlhRD2b+WJlXRe3hknoxPv1HqHUT6jKOvx+GIzWe7UKbWieeFC1dutWFHt361b19Yty1SSihLo0aMHbdq0YePGjYwcOZKIiAgCAgLIkycP0dHRDBgwgNatWxMSEvKv954/f57Nmzfz66+/omnaQ4UaHsUefl4hxPNN13W2XNiCUTPiV8aPmIQYpgZMZUjdIRTKZUeLqywW6NIFVq4EFxc15DxqlCqW8JyQVJRPkTIPHBwcTLVq1ahXrx6BgYFPnP9NUbZsWQYOHMjWrVs5duwY0dHRWdRqIYR4mMliYmnIUmrPrU2LxS3wD/AHwN3Fnc+bfG4fwddsht271fcGg8pa9fHHcOmSylz1HAXfp5EAjJoHXrduHR4eHhiNRjw8PIiJiSEwMBAfH5/Hvm/9+vWkjCCcPXsWo9GIu7t7VjVbCCFSrTixgorfVOTNlW8SlxzHD21/4Lc3f7N1sx6Ij1f7ditXhldfVRmrQGWwGjdOZbF6wbzwc8AA1atX58aNG3Tv3v2hY7GxseTPn5/Y2FhOnz5NsWLFUl+fPn06K1euZMSIEeTIkQMHBweWLFmC0WhHKwiFEM+1W/G3cDI6kdMpJ7FJsRTIUYApzabQrlI7DJqd9K/u3FELq77+Gq5fB29vWLoUZNGqzAFntRft5xVCWF/4nXCmB05n7uG5jPMdx0ifkVh0Cxoamr2kYkxOVnt1b92CkiXhlVfg/fehUaNslS4yo540Byw9YCGEyCZORp3EP8CfxccXY9EtdK3WlWZlmgHYT4/38GFVg/fcOVUKMG9eOH8eChSwdcvsjgRgIYTIJgZvGMy+8H28V/s9RvqMpHTe0rZukqLrsHmzCrxbt6pqRO+9B4mJanWzBN9HsosArOu6/QybZKKsHO4XQmRvuq6z8dxGpu+bzoL2CyiauyizX5+Nh6sHBXLaWUBbsUJtJypaFCZPhv79IU8eW7fK7tk8ALu4uBAdHU2+fPme6yCs6zrR0dGSqEMI8UTJ5mSW/rWUyXsnE3w9mGJuxTh/8zxFcxelYv6Ktm6ecucOfP895M+vEma0bauSaHTtCo9JXCT+zeYBuFixYoSHhxMVFWXrpmQ6FxeXh1ZSCyHE38Unx1N9dnXO3zpPlQJV+KndT3Sr3g0no50EtchI+OortZ3o9m3o2VMFYGdn9b14JjYPwI6OjpQubSfzGEIIkcWi46L549wf9KjRA1dHV/p49aFmoZq8XuF1+1lYBSpV5IcfqkQanTrBmDFqS5FIN5sHYCGEeBFdirnEtMBpzDsyj/jkeF4p8Qol3Uvyv1f/Z+umKboOO3dCpUpQuDBUq6bmdkeOhDJlbN2654IdPV4JIcTzL+JOBL1W96Lc1+X4NuhbOlXpxPGBxynpXtLWTVNMJli2TBVCaNxYpYcEaNECvvlGgq8VSQ9YCCEyma7r3E68jbuLOy4OLmw+v5khdYYwwmcEJfKUsHXzHpg9W20lungRypeH775T9XhFppAALIQQmcSiW1h3Zh0T90wk2ZLMgX4HyJcjH2EjwuxnYdWdO+Dmpr7fvh0KFVLzvW3bgqTWzVQyBC2EEFaWZE7ixyM/Uu3barT7tR1XYq/Qu2ZvLLoFwD6C79mzMGCAmt89cUId++knCAiADh0k+GYB6QELIYSVLT6+mHd+f4eahWqypOMSulTtgoPBTn7dBgaqYeY1a9Se3V69IFcu9VqOHLZt2wvGTv5FCCFE9nUt9hoz98+kfL7y9PHqQ/fq3fHM7Unzss3tK8FQTAz4+an0kB99BEOGqCFnYRMSgIUQIp0u3LrAlIAp/Hj0RxJNiQx9eSgALg4utCjXwsatAxISVIaqnTth8WJwd4cNG6BOHciZM8uaseZIBP6bThMZE09Rd1fGtKhI+1qeWXZ/eyUBWAgh0uGznZ8xbuc4HAwO9KrRi9H1R9tPqsjoaPj2W7Vt6Pp1qF1blQX08ABf3yxtypojEXy4Kpj4ZDMAETHxfLgqGOCFD8KyCEsIIdJA13W2XthK1D2VNrdO0TqM9hnNxWEX+b7t9/YTfPfsgRIl4OOPVaaq7dvh4EEVfG3Af9Pp1OCbIj7ZjP+m0zZpz5MEXg5kwu4JBF4OzJL7SQ9YCCGewGwxs+rkKibtncShK4f4vPHn/O/V//Fa+dd4rfxrtm6ecvCg6uE2b656u336wKBBULWqrVtGZEz8Mx23lYCwAJosbILJYsLJ6MTWXlvxKe6TqfeUACyEEI/x/aHvmRwwmXM3z1HOoxxzWs+hV007SUxhscDGjWpF886d8NJLKgC7usKsWbZuXaqi7q5EPCLYFnV3tdo9rDHHPGrzKBLNiYDaRrYjdEemB+CnDkFrmjZf07TrmqaFPOK10Zqm6Zqm5c+c5gkhRNaKT34QLDac24C7izvLOy/n1OBT9K/dHxcHOygpun69ys3cujVcuKASZ2zfbutWPdKYFhVxdXx4T7Gro5ExLawzZJ8yxxwRE4/OgznmNUcinvi+u4l3mR44nQu3LgDwbu13cTI4YdSMOBmd8C3la5X2PUlaesA/Ad8AC/9+UNO04kAzIMz6zRJCiKwVcSeC6fumM+/wPA68e4AK+SqwsP1Ccjnlso+tRLdugcGgCt0nJqo9vIsXQ5cu4Oho69Y9VkpPNLNWQT9pjvmf9wi8HMi6M+sIvxPO72d+JyYhBh2dkT4jebvW21TOX5kdoTvwLeWb6b1fSEMA1nV9l6ZppR7x0nTgfeA3K7dJCCGyzMmok/gH+LP4+GLMupmuVbti1FSPLbdzbhu3DggNhRkzYN48GDUKxo2D9u1Vtip7eDBIg/a1PDNtxXNa55gDwgJotKARJosJAN+SvkxsOpGXi72ceo5PcZ8sCbwp0jUHrGlaWyBC1/VjT3sy1DStP9AfoEQJO0o6LoR44d1JvIP3995YdAv9a/dnlM8oSue1k/rkhw+r+d3ly1Wg7dZN1eEF1RMWwNPnmEOuh1CtYDV2XtqJ2aJ6ykbNSPOyzR8KvrbwzH+LmqblAD4CPk7L+bquz9V13VvXde8CBQo86+2EEMJqdF1nw9kNDNs4DAA3Zzd+feNXwoaH8U2rb2wffHX9wfeffabmekeMUPO8CxdC9eq2a5udetQcs4ujRmOvcBr91Ijqs6uzP3w/vqV8cXFw+fccb0KCyoGdlJTlbU9PD7gsUBpI6f0WAw5rmlZX1/Wr1mycEEJYQ7I5mV9DfmVywGRCrodQ3K04H736EQVzFqRNxTa2bp6a0/35Z5g+HVauVKUAZ85U87158mTqrbN7lqqUtn68cRWX4w6Q00UnLkcQX+4/TXG34kxvMZ0qBaqQ2zk3W3ttfTDH61xWDefPmgVRUZA3L7Rrl6Vtf+YArOt6MFAw5c+apoUC3rqu37Biu4QQwiqOXT1G21/bEnY7jKoFqrKw/ULerPYmjkY7WLgUE6Nq7s6cCVeuqB5udLQKwFkwZZddslQ97SGhYL5LnNPHkuSURIzZTBmHMizqsIiuVbs+9PfsU9wHH48aalRh4UL14NO6tZpbb9Qoy3+upwZgTdN+AXyB/JqmhQOf6Lr+Q2Y3TAgh0ivqXhSXbl/Cu6g35TzKUaNQDWa1mkWr8q0waHYyf5qQoALtjRvQtKkaBm3WLEsXVj3LCmJbedJDwsvlDMzcP5Ofg38myZyEWTdj1Iy8U+sd3qrx1oOL6DpcugSlSqmKT0eOQO/eKhBXqmSDn0pJyyrobk95vZTVWiOEEBlw8dZFpgZOZf6R+ZTIU4KTg0+S0ykna7uttXXTlMOHYe1a+OQTVZFo8mSoVQu8vGzSnOyQpepRDwl3TJd4b93XxGh/YrKY8C3pS1RcFEnmJJyMTjQu1VidmJysFrFNnQrnz8Ply5A7N+zbZxf1jiUTlhAi2zsRdYIvdn3Bsr+WYdAM9KzRk9H1R9vH/l1dhz/+gClTYNs2FQDeeQeKFYO+fW3atKzIUpVRKQ8DiYaTJBiC0fQc3HKcg2Z2ZECddxjpM5JyHuUIvBz4YH7XrYr6vL/6CsLDVS/X31/tnQa7CL4gAVgIkU3puo7JYsLR6EjI9RDWnlnL8HrDGVFvBJ5u9jF8ysmTKlFGSAh4eqoeb//+mb6wKq3GtKj40PAuWDdLlTUUyePMydifiXFcBFjQcCCXuSWVcrzDt6+/kXqeT3EffIrVU0P4hw7BmDHQuLGaY3/tNbvcuiUBWAiRrZgtZtacWsOkvZNoU6EN/9fo/3ij8hs0K9OMvK55bd08tbAqNFQNK5coAfnywYIF8OabD3pgdiKzs1RlRJI5iZ+DfybceQIxSWdABzTQdROuWiE+alnvwckHDqhh5ty5VcKS2rXh1CmoaD8PEo8iAVgIkS0kmBJYeGwhUwKmcPbmWcp5lKNM3jIAGA1G2wffS5ceZKzy9FS935w5YccO27brKTIzS1V6rTq5iqEbhxJxN4KahWpSv8hIlp+ZhUVPxqA5MrJRB9rXKAxr1qjAu2ePGlX4z38eXMTOgy9IABZCZBP9fu/HkuAl1C5Sm2WdltGxckeMBjuYywsJgfHjYdmyBxmrRo3KNmki7cXa02vZF7GP1uVb4+7iTqX8lZjfbj7NyjRD0zSGXe70cJ7mTz9V+3hLllQPPm+/rXrA2Yim/z3zSibz9vbWg4KCsux+QojsK/JuJDP2zWCg90BK5y3N0atHuRl/k8alGtt+cZXFolbYOjvD0qXw7rvw3nswdCgUL27btmUzp26cYsyWMaw7sw4NDRcHl0fX4r16Fb75Blq0gIYN1YjD/v3QsSM42G9fUtO0Q7quez/qNftttRDihRWbFEvN72pyM/4mVQpUoXTe0ngVts1WnYckJsKSJWrYs3t3+OgjeOMNaNnSbhZWZReBlwOZtHcSv53+DQeDAxoaOvq/a/GGhKjP++ef1UNPzpwqAJcsqb6yMftbFiaEeOEduXKEG3E3WNppKX28+ti6OaoU4IQJKpHDO++o8n9VqqjXHBwk+KZB4OVAxu8eT+DlQABmB81md9huPmn0Cb+9+duj8zT37q2ygy1bpkYZzpyBDz+03Q9hZdIDFkLYnaBINVXVsERDG7fkvnfegdWr1fDnokXg5ydzvM9gZ+hOmi1qRrIlGWejM9t7b8e/mT+zX59NTqecACpP8/k/8T1vwadIHfXGl19Wi6kGDAAPDxv+BJlDArAQwu4EXQmiuFtxCuUqZJsGHDyohj0nTIDSpdVin08/hRo1bNOebOp2wm3mHJrDF7u+INmSDECyJfnhIWaAmzfxWbQdn29mq5zYhWqrHM2DBtmo5VlDArAQwu4MrjOYDpU6ZO1NLRZV/m/KFNi1C9zc1Dxv6dJSBjAdTBYT1WZXI/xOON5FvQm+FozJYnp4iDk2Fj74AH78EeLioHnzBzmxXwASgIUQdqd+8fpZe0OTSSVvOH5crWKeNk0NO7u5ZW07MkFWlBtMSQNZ0r0kJ6NO8lnjz3AwODCp6SQq5qtI7aK1H6SKLNkIH+3+SvEcOdTDTpcuMHLkC/egI9uQhBB25eKti5yOPp1aQD3TREfDunVqoQ/Al19CmTLQqZNaZPUc+GclIVCpJid0rG61IBx4OZDGCxqTaE4EwMnoxLEBx6iU/x9VhkwmNY8+ZQqcPasKI+TMqY7b8TaijHrSNiRZBS2EsCsrTqzgtSWvcS/pXubc4Nw5GDxY9XT79FHBANSWom7dnpvgC08uN2gNYbfDeHPlm6nBV0NjTP0xDwffu3dVUYTy5VVPNzoaPv/8QUGE5zj4Po0EYCGEXb7aJsgAACAASURBVDkYeZDS7qXJlyOfdS98+bLas1uhgkoX2a2b2mNavrx172NHMqPcYKIpkZDrIQAUzlUYdxd3HA2OGDUjLg4uvF7+dXViyujqX3/B8OEqPeeqVXD6tHoAcsnE0Y1s4sV99BBC2KWgyCC8iz5yxO7Zmc0QGal6u3nyqHq8H34IQ4ZAkSLWuYcds1a5wcDLgWw8t5EbcTdYfWo1Rs3IhWEXUoebHyoFeMMFPnxLfd6zZkG9enDsmKwgfwQJwEIIuxEdF83FmIsM8B6QsQvFxakKRNOmqXSRwcFqQdW5c3ZTCzYrWKPc4O+nf+eNZW9gspgAqOtZl/FNxuNoeDBU7+P5Mj7Hb0Lvj2D7dsiVSz3kpJDg+0gSgIUQduPQlUMA1ClaJ30XuHZN9bq+/VbNNdatq+rC6rpKnPECBV/IWLlBi27BoBlYfXJ1avA1akbaV2yPXxm/1PPWHIng5ugPeXvbIq675ef6sI+o9ulocHfPnB/qOSKroIUQdsNkMXEi6gTlPMqRwzFH2t+YEmAXLVKrmtu2hdGjoUEDyVj1DHRdZ0/YHiYHTKZqgapMbDqRgLAA/Bb5kWxOxsnopAoluJSD2bPZWdKLAeecKHA9nFqRp1hfqSGOLs5WXWWd3UkxBiFEtuBgcKBGoTQOV+q62kM6ZQo0aqQCbteuKn1hhQqZ29DnzJ6wPXwX9B1Hrx7lr6i/yJ8jP74lfQGoX6I+23ptU3O8htL4fPETLFwICQmcav428bU6Epa3CGF51Zy66f4q6+wUgOPj4cgROHBA/Xn48Ky5rwRgIYTd+O/W/9KyXEteLfnq408ymWDlSpUq8uBByJ9fZVACcHKS4PuMUvbxmiwmNDRG+Yzis8afPTQC4VPcB5+P58JP/1Vz6r16wYgRTFxw4ZHXzMgq68xmNsOJEyrYpnwFB6vjoJ7fJAALIV4oV2OvMmHPBPLnyP/kAPz222qouUIF+O47FQxcn21V7/MivVmuYhJi+C7oO9pUaMOO0B1YdAsABs1APtd8KvgmJ8Pvv0P79mruvEYN+OQTlZ+5YEEAirpfscoq68x05Qrs2/fg69AhuHd/i3mePGqZwAcfQJ066qto0axrmwRgIYRdOBT5mAVY4eHw9ddqVW3x4moPaadOKlm/4cVNZfDPLFcRMfF8uCoY4LFBOPxOODP2zWDOoTnEJsVi1Iz4lvLF2ehMkjlJ5WnO7w3+/jBzpvrs16+HVq1gxIh/Xc8aq6ytKTFRDSWnBNvAQAgLU685OoKXl3p+q1tXfZUr9+Cf0JojEXRemLkpO/9JArAQwi4ERQahoVGrSC114NgxNcz8yy+qUEK1atCzpxojFE/McvXPwBF4OZCRm0ZyMPIgAF2rdWVM/TF4FfYC7pcCPLMJ340n8anbURVJaNwYZs+Gli0f24aMrLK2hqtXISDgwdehQ5CUpF4rUUJtQR4+XP23Vq3H5/5Iz8OMNUgAFkLYhaArQVQuUJlcRlfV49q4UeUKHjwYhg1TVYlEqqdludJ1nYORBzGZTTRd1JQEUwJGg5GlnZbSsXLHB2+4dk3N8RatC4OrQrt2qjDCSy+lqR3ta3lmScA1m1VSrb17HwTcC/enoJ2dwdsbhg6F+vXVM9qzDCU/y8OMNUkAFkLYXmIiVyNO413OR803VqqkVjb37w9589q6dXbpcVmuiuRxYtXJVUzeO5n9EfvpV6sfSeYkdHR0Xef0jdMqmq1dq0YYTp+GS5fUPPrx42ohmx1ISFALpPbsgd27VcC9c0e9VqiQ2mE2aJAKuC+9pIJwemVGys60kAAshLCdmzdhzhz4+msOXrlCUvAydXzaNNu2KxtImX+NMYeQYAjG2VIZzSGSUMd1vLHsImXyluHbVt9SOX9llgQveTDHezga3q6ksoKVLAn//e+Di9ow+N66pXq3KQE3KOjBcHLVqip1d4MG6qt0aetu77ZWys5nJQFYCJH1oqJURZwffnioELtT1Zq2blm20b6WJ6duHuKjPf/Doidj0BxxdXChkls5vnptEh0rd8RoUJm/tvbaqvbx3nbHp80gtQJp6VLo2NFm1YiuX1eBdudOtZ37+HG1tdvRUQ0nDx8Or7yierj5rFyX459stZhMArAQIuvcvq32fjg4wOLF0LkzjBzJV3Hb2R/xE0tohuStSpuw22GsDZ2EhUTQQNNMDHl5ABP8JqCldA9DQmDaNHzy5sVn6lQV4YLqqjHbLM4QFhGhAm1KwD15Uh3PkUMF2XHj4NVX1bNBVu8qs9ViMgnAQojMZTar/aRTp6re7qFDal43pSA7sHHx+1yJvfIgcIjHCr4WjH+AP7+E/ILFYsGoqV6uk9GJdhXbqQeYLVvU571pk4pmgwerN2sa1K6dJe28elXVZdixQ/03peyym5vq2fburab5X3rJPqads2ox2d9JABZCZI6/VyQ6dw5KlVLjimaz6gHfD766rhMUGUT7Su1t21479fdSfzo6DeY3IKdjTobUGcLwesOJvBv5oBRgcR/4+GM1vF+4MHzxBQwYkPljuKgh5R07HgTcU6fUcTc31bMdMEAFXC+vF64mxmNJABZCZI5ly9Qy1bp11fcdOjxyvvHS7UtEx0dbrwawlaQ3y5Q17Qnbg99CVQjBxcGFLT23MLPlTHrU6IGHqwcAJS258dmzA3I4QXGgRw+1Sql794wtDX6KO3fUcPLWrbBtm0rnCKoSYcOGKuFF48Yq4NpomtnuyccihLCOEydUb/ell1Tg7dZNpRp6SkWioEhVIc2eArCtEjOkiE+OZ8GxBXy09SOSzGopcJI5iV2XdvFhww/VSefPw/Tp8OOParRB19XwcsWK6svKEhJUZqmtW9XXwYNqMMPFRQ0pd+umAm7t2mohlXi6pwZgTdPmA62B67quV7t/zB9oAyQB54G+uq7HZGZDhRB2SNfVeOOUKSpxhouLSkEEqvf1yitPvYRRM/Ky58tUL1g9kxubdrZKzACw8sRKBq4fSFRcFJXzVyY2ORazxay2EJXyVSf16wfz56uuZY8eKnFGdet+fhaLWpm8ebOaUt6zRwVho1HlTP7gA/DzAx+fx2eYEk+Wlh7wT8A3wMK/HdsCfKjruknTtEnAh8BY6zdPCGHXBgyAuXNVcv7PPoOBA1V1omfQoXIHOlTukEkNTJ+sTMwQeDmQNafWUL94fdpVakcxt2J4F/VmbIOxvFryVfaF72PHhW34Rjrh43k/DWf58vDhhyo/dpEiVmtLeLgKtlu2wJ9/qt1iAFWqwHvvqYD76qtqIbs12cNwvy08NQDrur5L07RS/zi2+W9/3Ad0sm6zhBB2KSYGvv9e9bqKFlXjjnXqwFtvpasbpOs6Ft2Sul/VXmRVYoZFxxbR97e+mHUzRs3I7r678Snuw4YeG9QJd+/is+oAPjPmQWgoeNSAFi1grHX6O/fuqXncTZtUTzdl4VShQmprdrNm0LQpeGZiLLT1cL8tWaOUyNvAxse9qGlaf03TgjRNC4pKeZwSQmQvoaGqGk7x4vD++7BunTru66uGQ9M5Bnnu5jnyTMzD2tNrrdZUaxjToiKujg8/FFgzMcPO0J28tuQ1eq3phVl/MNS9I3SH+iY2VgXZ4sXVynFPT1i1SkXDDNB1Nazs768u5eEBr7+unqlKllQzCceOqRJ+ixerrUKZGXzhycP9z7sMLcLSNO0jwAQsedw5uq7PBeYCeHt76xm5nxAii1ksqne7dKmq29a1K4wapUrLWEFQZBD3ku9RIk8Jq1zPWjIjMYPZYk7t6S/9aymHrxzmvdrvsfDYwgdpIvPe/1xdXGDFCtUNHTUqQxWgbt5UvduUXm5kpDperRr85z+qQ92woe3mcW2Vh/khJpNNlmqn+46apvVGLc7y03VdAqsQzwuLBfbvV6trDAa1X3fUKFVqplgxq94qKDIIFwcXqhSoYtXrWoM1EjMEXg5ky4Ut3E28y+pTq1nQfgENSjTgiyZfMK3FNFwcXOhdoyc7Ns/Fd10IPrPehosX1QK2kJB0pYSyWFRN3I0bYcMG9VdpsajcJ82aqeqCzZtnfs82rWyVhxlQXf3Zs9U6ht271dx6FkpXANY0rSVq0VUjXdfjrNskIYRNpCTOmD5dpS06cQIqV1bjk5kk6EoQXoW9cDQ+f/tWNp3bROtfWmOymAConL9y6mserh5qSfFP8/CZNg2fkydVRExJVALPFHxv3VK9240b4Y8/4No1tfPL2xv+9z947TU1VW+PCTBskof58GGYMQN+/VX1flu3huTkzLvfY6RlG9IvgC+QX9O0cOAT1KpnZ2DL/dRx+3RdH5CJ7RRCZJZbt1TQ/fZbiI5Wv6l//TXTewNmi5nDVw7Tp2afTL3P42TmyluLbuHNFW+mBl+DZuCtGm/RoESDBycFBMC776pMFSl5sdOYk1HXVW3c9evVdHxAgOrlenioIeXXXlP/LVjQKj9OpsryPMwxMSr5tIODWsX/n/9kec83hZaVo8fe3t56UFBQlt1PCPEECQlq4u/GDZUmsmlTNdT8yitZkqg/LjmOyXsn07BEQ/zK+GX6/f7unytvQfW6JnSsnu5f/EevHmXB0QVMaT4Fo8HIxD0TGbdzHMnmZJyMTmz1nY/Pwu0qLeT48SqKBgaqof40fN4JCWrL9bp1KvBeuqSOv/QStGqlvurWtc9erk3duaP2TB84AD//rI5t2qTm1d3dM/32mqYd0nX9kVlmJAAL8SLRdZU3cOpUtTonMFD98o+OzpJ8wfaiwcRtj5x39HR3Ze8HTdJ8nYCwAH448gMhUSEciDhALqdcBLwdQPVCKilGYFgAO3b8iO+Gk/gs3avmdgcPVp9/GkRGqoC7bp3KPhUXp6oHNWumVi+3amU/c7l25/x5+PprFXzv3lUPluvXq+TUWehJAVhSUQrxIkhKUsPK06apfSaFCqkkDimFEWwQfM/dPEeBHAXI42LlrA5pYI2Vt7+d+o0OSzugozoxA2oPYLzfePK65k09x2feH/h8Pk8lJ/nkE5Wi8wnjwrqu/nrWrlUFpFL6K6VKqdzKrVurggaSeeopNmxQH5bRCG++CcOGqQlxOyMBWIgXweLF8M47KqXRDz+oRP02/i3ea3UvHAwO7Oq7K8vvnd6Vt3HJcZyIOoF3UW9CroekBl+jZqREnhLkTTLArClqmXGNGmrbVrFi0LPnIxdVrTkSwaT1Z7h4PAfaZU9MoYW5cdUBTVMjpOPHQ5s2ULVqlpfvzV4SE9UDZu7c0LGjekr5+GPo318ljLFTEoCFeB5duKBWeXp5qa5Tt27qF1GLFnbxm9xkMXHk6hEGeg+0yf2fdeVtdFw0sw7O4usDXwNwecRlmpRugutuV7WH1+CI7+/H4PXiargzMZE15nz4b7pGZIwnRb8KfGhhUUwMfD77Jt8vNhB7riF6kgOao4mcZaIZ3N+J/xuUl0KFMv9zyPauX4fvvlMLCK9dg/btVQDOmRM+/dTWrXsqCcBCPE8CA9X84urVag/v+++r466uagOonQiKDCLBlECdonVscv+0rrxdc3INkwMmc/jKYRLNibSu0JqxDcbibHTGp7gPW3ttZceMYfiuOIRPxIrURCVrtEL/Sq84+qczbFyWk4uH3Nm+HUwmDww5E8hZORLX8tdwKXEDg6OFI66uFCqU9nnoF5a/P/zf/6neb6tWagtXBjKFJZuT2R66neZlm1uxkU8mAViI58Xgwaon4O4OY8ao7RV2uEJnxr4ZvL/lfXI45qBhyYY2a8eTEm2YLWYORBzgzZVvkmhOxKgZWdxhMT1q9FDz5ps3Q/Pm+BT3wcftdejSSCUqKV4cAP+J24hPNpN8IxdxZwsRd7YwSVfcOQ9UqKAWmy+I2ItT0Zh/DUhkaQao7MRiUXO79eqpOfVy5aBvXzW/W6lShi696dwm3l37LpfvXObYgGPUKFTDSo1+MgnAQmRXsbGqFmzXrmphT7t26hdR376qKrodORhxkGJuxSiSuwjVC1bn3ZfeZVT9URRzs25mrYzQdZ2tF7cyee9kKuSrgGduz9R9vABhN86pB5zp0+HcOVUuyM9PLa5KvYbK8fDX7yW4d7owppvq78GpyC3cG50iR/mrnP7eF4DdExOJeEQR1yzJAJWdxMaqBDFffaUSxEyZop5gOnRQX+kUcSeCZEsypdxL4enmSVmPssx+fTbVClazYuOfTAKwENlNRITaXjFnjppMdHJSteKaN1dfdiIloE3cM5GtF7cytsFYJjadiF8Zvyzf9/s4gZcD2XZxG2bdzJpTazhy9QiFcxXm9fKvU9ezLk5GJzXHa9Hw/c9UOHH3QaKSRo0A1SEOCFC1ElatgrAwwFAGl+I3casdimv5qzjkTgTUNqcUNskAlZ3ouio6PGcO3L6tVqX98gu88UaGLht8LZgpgVP4OfhnOlXpxC9v/EK1gtXY3nu7lRqedhKAhcguzGa1knnJEjUc16GD6gn4+Ni6Zf+y6uQqxu8ez6ErhyiSqwj+zfzpX7u/rZv1kMDLgfgt9CPBlICOTgm3EsxrM4+3aryFs4Mz3L6t5njPb8V31Ex8yteH71SiEpNZY8cOVS9hzRq1/sfZWT3/jBsHhpLXmLDt2BODa5ZngMoOdB3OnIGKFdViwdOn1cLBESPU0HMG7AzdyYQ9E9h0fhM5HXMyyHsQw+sNt1LD00cCsBD2zGJRY5re3mpPo8mk9pIOGwZlyti6dQ8xWUw4GNSvlBUnVnAn8Q7ft/menjV6qoBmJ27E3WDWgVlExUWRZE5CR8egGXjP+z3eqfW2GlqeOhVOnMDn/Hl8ivvAzuEkO+di2zZY0V+tcYuOVottW7VSnbJWrdQuGKUIbu6WpwZXaxR8eC4kJ8PKlWp4/+BBNcRfpow6loHUXsnmZBwMDmiaxoazGzh69ShfNvmSAd4DVD5uW9N1Pcu+ateurQsh0iAuTtfnztX1ypV1XdN0/cwZW7fosW4n3NYn75msF51aVD9+9biu67p+M+6mbjKbbNyyh124eUEfsn6I7vqFq86n6APXDdRdv3DVjeOMuusXrnrAd//T9Ro1dB10vXBhXf/iCz3h5j193Tpd79NH1/PmVS/lzq3r3bvr+urV6q9JZMDt27o+caKuFyumPtzy5XX9m290PTY2Y5dNuK1PDZiqF59WXN94dmPqsfjkeGu0+pkAQfpjYqL0gIWwJzExav/ut99CVJTax7twoaqWbmeuxV7jq/1f8e3Bb7mdeJumZZpi0S0AD2WDspXAy4HsCN2BbylfFh5byNzDczFqRnrW6Mmo+qOoUqAKPWv0VOdcc8Gn80ioVg3T3Pn8WbA7v652Zk1pNf2YJ49a49apk0oDKZmoMigpSa1duHdPbSVq2FCVBWzVSm2fS6eIOxHM3D+TOYfmcDvxNo1KNsLNWaWeTPmvPZFc0ELYg8RENYl47ZrKO+jnp+Z3fX3tInHGPyWaEik2vRjRcdF0qtKJsQ3GUrtobVs3K9XEP3/jv3s7o+tmDJojzUp2oaZnEYa+PBRPN0+VJ3jGjNQUkclJOkdm7OS7k41YvUYjJkYF3Q4dVJGipk3TXKhIPI6uq4TWM2aopNbbtqnjERFW2S6n6zrlvy7PxZiLdKrSidE+o6njaZt95n8nuaCFsEcphRGmTVPZk3btUjmaL12yyzpyR68eZWnIUsb7jcfZwZlvW31LzcI1qZCvgq2blspkMTF2/Ry+OvR/6FoyaGDRkzl8wciAmkPxDLkEU4fC6tXoDg5cfn0gn78Lq1Zp3Lzpi5ub6ul27ap6uhJ0rSAhQVUhmjEDgoPVv+1Bg9SiQqMx3cFX13W2h25n/pH5zG83HyejE3PbzKWUeynK5LWv9RGPIwFYiKz2z8IIBQuqJBopv5DsKPjqus7OSzuZuGcim85vIrdTbvq91I+yHmXpXLWzrZuXKtGUyNxDc5m2bxqhMaEY9QKoX28WNBwwmqpya/QHsG0xybnzsu2lDxkTOpjgNUXJlQvatlVBt3lzGV62uh9+UIU/atRQ+9a7dVOjPemUbE5m+YnlTAmYwpGrRyiUsxCnb5ymeqHqNCmdvTKISQAWIqv9+KMqBF6lCsybBz162OVv/dCYUN5c8Sb7I/ZTMGdBxjcZz8A6A3F3yfwaqmmRsoe3SekmeBf1ZmrgVIq5FSPuWk9cLHXQzMcpG/Ubd1xfJTz+ZX64VpHz7vWYHtMHy185adMGPumqph0fUSdBpNfRo6q36+enilD06qX+rVthOiX8TjgN5jcg7HYYlfJXYl6befSo0QMXB/v7/yctJAALkdnOnVNZfOrUUb+MevRQi6rspDDC3yWZkzh/8zyVC1SmSK4iOBgcmP36bHrX7I2ro/1EqeV/Laf7yu6YdBOuu13Z2msrB989SIGcBWj/4VJa7lhAtyN/kCfpHh+7tuTz+IZcM1hwf83A3G6qx/tgy5DIMItFFS2ePh127FBFi6vdzyiVOzc0bpzuS0fejeTo1aO0Kt8Kz9yeNCvTjHYV2/F6hdcxaOlfsGUPJAALkRl0HfbuVcPMa9Y8XHM3Vy67KowAcDfxLnMPzWX6vuk4GBw4+5+zODs4s+ftPVa/15ojEelOPnEo8hCTAyaz/K/lqaUAk8xJ7AjdgU9xH+71HsSKxd+jWSys5A2mMYJjBctTuGoIk8d40NPXfkvTZWtduqg9u8WLw+TJ0K8f5M3YSviQ6yFMDZzKkuNLyOWUi8hRkbg4uDCv7TwrNdr2JAALkRn691fDyx4e8N//qjneIkVs3ap/iboXxcz9M/nm4DfEJMTQuFRjGhZ+h0aTd3HldoLVszOtORLxrypBH64KBnjqPYIig6jzfR3cnN3oUb0HK06uINmcjBMO3A1uhN+n0HSbO64M5rdyAwmtqGEqFcZLxY7f/xmsG3wz8iCR7YWFqa1yY8eqQNu/v1ou3rEjODpm6NLHrx3ngz8/YOO5jeRwzMF7td9jhM+IbDvM/CQSgIWwhjt31GKTXr1UT7djR6hVC3r3VumS7Iyu62iaRmB4IF/u/pIOlTswtsFYIq973g+QCcCzBci08N90+qH0jADxyWb8N53+1/V3X9rNNwe+oUjuIsxoOYPaRWrzQ9sf6FSlE84JTvjtLMbJv36g/ekoPghPJrIcJH4yno7dYHhqxsfMWaGdkQeJbG3fPjXMvHKl+nO9eqoGbwZzkJssJu4k3sHD1QOzxcyhK4f4vPHnDPQeSL4c+azQcPskAViIjLh0CWbOhO+/V1uJPDxU0H3tNVu37JGOXT3GpL2TKO9RnnGNx9G6QmtODzlN+XzlAWiwaFuaA2R6PK7U3t+PxybF8r9t/2Pm/pmpaSI7V+lM/eINqBTVhR2vTKPeoVn00a8T7ODFUb9pTF3qQ22frJtSf5YHiedCfLxaVBUYqDZIjxihyl2WKJGhy8YmxfLD4R+Yvm86r5Z8lYUdFlKrSC0uj7iMk/H53wMmAViI9DCZ4K23VDZ+UHtYRoxQOZvtTMpWokl7J/HHuT/I5ZSL9+u/D4BBM6QGX0hbgMyIou6uRDziWikl+FaeWMm7a9/lVsKt1Nc0ND6d+yfnf25A5EUHQplNZFFvLg0dRa2RjanumPUL2TL7c7ILMTGweze0aaOWiVevDt27Q58+GS53eeXuFb4+8DWzg2YTkxBDwxIN6VK1S+rrL0LwBcjeS8iEyEpmM+zfr753cFAp80aOhIsXVYUiOwy+AB9t+4jGCxpz+MphvmzyJWHDw/i/Rv/3yHMfV4vWWjVqx7SoiKvjw8n1jY7X6PmKmt8rn688jUo1YqrvHBxxRbMYMCbp9F8wlwplTHy/0IVc4afwilhPnbFNcLBB8IXM/5xs6tw51bstVkxNpVy/ro7PmaP281qh1rR/gD8T90zEr7Qf+97Zx66+u2hdoXWGr5vdSA9YiKdJKXw/YwaEhqo0hqVKqew+dijRlMiS4CW8UuIVKuSrQOcqnSmRp0SathJldo3alOHZjzeuIjRuE5pTBLH6cfZGdeftpPpcCKyB8edlHFm/gu+KlOBaqdPUu5YH78796PxFEuRwAPJYpS0Z8VzW8j17FkaPhrVr1QNm9+4wfHiGE8Pous6uS7vwD/BnRL0R+JXx4/0G7zO4zmDKepS1UuOzJwnAQjzOjRswZcqDwvc+PjBxouoZ2KE7iXdStxJF3o3k41c/ZlzjcdQqUotaRWql6RpZUaP2ctIq/jINx+JoAR1aefbEacckPPurj7yH+2YWmbuT4FARl55zVDIHO8uU8dzU8k1MVEU/ihVTn/HBg/C//6lUkYULZ+jSJouJVSdXMSVgCgcjD5I/R36u31O96cK5Mnbt54UUYxDin+LiVCKByEgoW1bNgY0YYZeF71N8vvNzpgZO5XbibZqUbsLYBmNpVqYZmp0k+vh7XdaGPzZkT9j9/cUWI3m2jWRcgIn8VQviPuEDWjSz4LB9i0rGnIHKOOIJoqLgu+/UVqJq1WDLFnXcZFK9XytoML8BAZcDKO9RnlE+o+hVs5ddJXPJKlKMQYinsVhg/XpViN1gUEUSihZVlVo87KBw9yNcvn2ZYm7F0DSNqLgompVtxvv137eLCjAppQDretYl+How0wKnMfu1eSSeaE7ylk+gTFswJOFsgfVhU6mvaWiNBsLrAAaVJUxY36lTKjnMokWqSELLlurhMkUGgu+12GvMPzKf0fVH42h0ZHCdwYz2GU3bim0xGoxPv8ALSAKweLHdu6fq7U6frubAiheHYcNUQDYY7DL4Hr5ymEl7J7HixAp29N5Bw5IN+arlV3bT2w28HEiThU1INCWmZqvyTH6VHp1zczsEihRpysyYN4m9+yO+UTnw6TIIhg5Vn72wPl1X/56NRjW/u2iR2q8+bJjK0ZxBp26cYlrgNBYeW0iSOYn6xevTqFQjulfvboXGP98kAIsX27x5aqFJnTqqQpEVMvlkBl3X2XpxquAXLQAAIABJREFUK5P2TuLPC3/i5uzGmPpjUrcQ2UvwBdgeup0Ek0rkgQ4EDeDupqlMqb6ACnPz0aBvBRyODoJd1VTKQjf7K5T+XIiPh8WL1eLBDz5Qc+kDBkDfvqoOcho9LuPXrfhb9PmtD7+f/h0XBxf6ePVhpM9IuypPae8kAIsXy7FjagiuWTO1j7dvX6hdGxo0sLvCCH+XYEqg28puOBgcmNR0Eu/Vfo88LrZfDZziQMQBvj80jzdcZrH9l8ZQzBkMJgwWR+a6JNDXrQSGI9Fw/Qtw+Eht2bLTbVvZ3tWram539my1qs3L60Ee8ty5n6kKxT8zfoXHxDJy1TqgNW29inAr/hYfv/oxg+sOpmBO+ymjmV1IABbPP4sF/vhDBd6tW1VqyBo11GtubvDKK7Zt3yPEJ8ez4NgC1pxaw/ru63F1dGXzW5upXKDyE3PiZlV+4sDLgWwP3Y6z0ZkVwWvZd3UnhkR35v04iLyJPrzRezuVkkby+sYgfEIXqCr3o0apBx2Rudq0gUOHHiwebNQo3Q+XKRm/LCRwz7iVOw5rsGh3mPiHJ+1rebKzz067Gn3JbiQAi+ffm2/C8uXg6am2EfXvn+FKLZnlVvwtZgfN5qv9X3H93nXqetbl2r1rFM1d9KlbibIqP3HKHG/qMHNsAdgzjYa53+G/A/4ipnk4X+2Mp+WqwkQWasWW6aNp1q6h1e4v/sZigQ0b1Fa5JUvUA+XMmWqIuXz5p7//KS7HXOWOwzruOqzHot3ByVIR96Q+XE2wAPY19ZEdSQAWz5+rV9UWi+HDwd1dDTO3b6+qtdjh/G6K4GvB1J9fn9ikWFqUbcHYBmPxLeWb5l9ymZ2f+G7iXdYePMp3G/eQ4JCs8uhZDDR2GcTSTvkosOhV+PQY3S9OIaJwJb5o0g8A16B7TCgRkf32yNqze/dgwQJVZ/rMGbWP98wZNaxvhe1yFt2CQTPg7naLy8m/4Gp+GTdTB5wtVdHQ8HweMn7ZgacGYE3T5gOtgeu6rle7f8wDWAqUAkKBLrqu33rcNYTIEsHBajXzkiWQnKyGmTt2tNvCCAAno05y9uZZ2lZsS5UCVehXqx+9vXrjVdjrma+VWfmJQ29cZdjPM9lwfTYmswXD8tUYezihk4SzwcCXv8+iwNEbULUqX74xhqD8D2c3eq6LFNjC9etQqRLcugV168Ivv8Abb1jl4TLwciD+Af4UyFGAOW3m8Nlr7Rm1Khfm5AeLtrJ9xi87kpZd7j8B/6we/gGwVdf18sDW+38WwjYSElQ5tBo1YOlSePddOH1aBV87FXA5gHa/tqPKt1UYsmEIZosZo8HI9JbT0xV8wbr5iQMvB/L2L+9T9qN2lJ5Zkt9vTsQ50o+BuTfz41elqeI8iXwJXVmz2JGyLlXUHHtwMPPKNSLJ4d+B4LkqUpCF1hyJoMHEbbTp8xWTO4xgzZEIlRpy6FDYu1eVB3zzzQwFX4tu4bdTv/HK/FeoP78+O0J3UMxNZXtrX8uTqR2b4unuigZ4ursyoWN1eZiykqf2gHVd36VpWql/HG4H+N7/fgGwAxhrxXYJ8WTx8ao0WpMm4OIChQrBhAlqftcO9+6m2B++n9FbRrMnbA8erh58/OrHDKk7xCqJCqyRnzg+Hj5fuItJES2xaIngaMHzbju+bOxPz6ZRXP1kHAnTgvl/9s47rMryjeOflwOHISooKg4QtzhzlIIj3Km5V67UUiv3IhvashyAuEdu/ZkrRU3LMjVcIDlw4wbFLQqCzDPe3x+PI8vBOHA4+nyuyysEzvvcnPC93+d+7vv7/frDuTgopfm01buk5nNmUuEqtFeUl7odSdLPpkNX2OG3iMDQIOpcPcU9+3w0qfA2AO2/+cZk63wb/C3f7fkODycPZrwzgw9qfICj9onhQvsaxWXCzSYyewZcRFXVGwCqqt5QFOW5/eeKogwEBgK4Z9E7UiJ5asQiLg6uXIGiRYW4QC4lzZBGki4JJzsndEYdV+5fYXqL6XxY88OnbnRZJSv6xKcjjHyxZBu/xvmht3oArmlgZUSjaBhcNg99fngfDhzA0T4vG6u/g60+jRQbO+Ls88E/SsyvpEmBOdi7lzfbd6P9vRtczVeYCY0+ZG31FjzAJsvl/LtJd5l3aB5NSjXBy82Lvm/0xbOQJ50rdcbaSrYF5STZ/m6rqroAWABCCzq715O8okRHw9dfPznfbdNGWAFmUTA+O3mQ9oCFhxcSeCCQNuXbMLf1XOq71+fisIvZdqPLyG4lNRV+Dkrj+82rOVvQHwqfwsHWjQ4eXdh6K4I0QxpaNPhMXAXa0jBrFl6XXEnU/nc3+6jE/MqYFJiDK1eEDnnFiuDmxnUHZ75v0Ift5b0w/KNCktlyfmRsJIGhgSw5uoQkXRJG1YiXmxelnEtRyrmUqX4KSQbI7F3glqIoRR/ufosCt00ZlEQCiBGL2FghIqCqsGGDUE4aPhzK5161nduJt5kZNpO5B+cSmxLL2yXfpl2Fdo+/bu5dxvqwUGb9EszxX3yIK7IFGkzCVanKuEb/Y6BHfWzmzCe04ACC67ri494Arxqx0KoVaDQ4Td5F4ktKzLJkmUH+/lvMqK9fLxoGt2wBDw9GDJ5lsnL+kN+GMO/QPDSKhp7VejLaazRVClcxRfSSLJDZO8EvQB9g8sP/bjZZRBLJIwm9adPA3V00+Li7w40bwqUolzN+13gWHllIB88OfOr9KXVK1DF3SBgMYlz0m1VbOVKuA1gbsepoy/iKP1G31jZaJrmiBAbC6n5gNOL10Ud4+X4uXlzyyXVkidmEbNsG338PISFifnfkSBg69PGXs/JeG1UjOy7toEmpJmisNJR2Ls0YrzEMqzOM4vnkw1FuIT1jSKsRDVcuiqJcBb5GJN51iqJ8CFwBumRnkJLXhFu3xPnu3LlCQq9GDSEXqapCySeXJt/D1w/jF+LHqLqjqFOiDl82/JJRXqOo4GL+pHTrFixeDLNXn+VGqQCosRQUAyigKGnYlzhDq59S4LvvhELY4MGiwlDq2SVJWWLOIvHx4vfY2hrCw8VD5fTp8MEH/5GIzMx7napPZdWJVQSEBnD6zmk2dttI+4rtGeU1Klt/LEnmkH7AEvPzKMFOmQKff/7kfLdhw1yrz6yqKn9e+hO//X7sjNxJPtt8zG01l57Vemb6mqaSkVRV2LdP9KmtXw+6Fh9BrYVorWxpWa4lf1zchs6QhlajZWefXXhdBfbuFR3kTk6Zjl/yAqKihELVokWwcCF06ybG52xshEtRFknVpzLtwDRmhs3kxoMbVC9SnTHeY+hWuRs2mtwrPvM6IP2AJbkPoxH++EOcffXrBz16wEcfQYcOufp8F0Tybfa/ZuyM3ElRx6ImMUcwhYxkYqLoUZuyaj+XnBbicO5DPvmkAY5NqqDJO44hZbpTePl6Qn8OJjh/Cj4t3sfLzQvcMIl6kuQZhIaKo5QNG4S9ZZcuTywA7Z6v6Z1eknRJONg4YG1lzZLwJVQpXIXl7ZfTtHRTKRNpAcgELMlZkpPFyND06RARIUzvDQ/PuJyccu0OLEmXxLpT63i/+vtYKVZ09OxI9yrd6VWtF7bWtlm+flZkJC9cEFX7xcvSiH/jO/CZCIqKvtZq3usbjJfbUHG++GMtSE7Gq1UrvEaPhkaNshx3ZkjSJbH86HL2R+9n/rvzMzWKlVOmE1nCaBQPl7duga8vDBkiJCNNwNGbRwkICWBn5E4uDruIg40DhwYeIp+ttHa0JGQCluQsLVvC7t3ifPd//4OuXUGrNXdUz+Vu0l1m/z2bWX/P4m7yXUrmL0mjUo0Y9OYgk66TURnJRwZPs2fDtu1pWHnNxHbQdLC59vh7DEYDwVHBYpeblgY9e4pEbAIT9qzQZnUbdkXuAqCTZyc6eHbI0OtzynQiw9y/L0rMK1eKkr6jozgD8PAQH2eRR8ce/iH+7Li0A0etIwNrDiRVn4qDjYNMvhZIeqQoJZLMc+yYKC0nJIi/f/klBAcLu7RevXJt8o1PjWfYtmG4T3fnm93f4O3mzd5+e/Hx8Mn0NR/JCpb67FfqTd4lZAUfkl4Zybg4UbUvXx5at0shPBy+Gm9NqU6L8CpfgcCm/tgrWjRG0OoM+CQ/1MiZPVucPZoh+Z67e44hvw0hLiUOgHENxrGj9w7sre0JjgrO8PVeVC0wC5GRwvijRAkYM0Z0NN9+OJlZpYpJki/AoeuHaLGyBadun2Jyk8lEj4xmaoupONvnTmcvycuRO2CJ6TEaxYhFYCDs2iW6Prt3Bx8faNbM3NG9kNjkWJztnXGwceD3C7/TtXJXxniNoXLhylm67st2bS8bOYmIgFmzhAFOUunV5GkbSN4CkZwaEUUBR0dGxf5F/iWrYNwM6qppBNd2wafxB3g16CEulsPngaqqsj96PwEhAfxy9he0Gi3vln+Xd8q+Q6NSovRdz70ewZeDM3zt7DKdyBQXLkCFCuJ8t1s3UWGoVcskl05ITWDhkYUkpCbwtc/X1C5Wm43dNtKybEuTHHtIzI9MwBLTEhcHdesKMwQL8N8FkSz+ivoLv/1+HLt1jMjhkdhZ23Fy0Em0GtPs0F92xvuskZPRzSpgfa04zcfCn3+CTekQCgz7jCS7vSQC1kZrwqJ30dKzLflt8oqZ0urV8Ro9B6/WrUVSMANJuiSarGjCgasHKGBfgHENxzH4zcEUcSzy1Pf5N/Mnj02eDF/frHrTej0EBYmu5k8/hbJlRXdz+/bi990EXE+4zsywmcw/NJ/7qfd5p+w7qKqKoii0r9jeJGtIcgcyAUuyzvXrQkygc2fRRNWkiZCN7Nw5V/vvGowGNkRswG+/H4dvHKZIniIMqzMMg1EkSlMlX0jfru1RIr5/H5YuhU+7wcWL4r4+9PvjzNLXI97aHkWvoKKiGg0c/WogLde+K8qcERFmk+ZMTEskJDqEZmWa4WDjQJVCVehVtRd93+hLHu2zk2xmXZ/MIgby6Hx35kwhGVmlihiVs7YWs9MmYvnR5QzYMgCDaqCTZyfGeI/hreJvmez6ktyFTMCSzBMeLkYs1qwRs4zNmkH+/DBnjrkjSxf7o/fTbX03yhUox4J3F9C7em/srJ8/GpKVztv07NrOn4cv5oXyy7Fg0i54U7bORbqNusH/BnyJjU016p1YRaGjF3j3/DekoaI1qviUbiLmSR0czJJ8bz24xZyDc5hzcA7xqfFEj4zG1dGVhW0Xpuv1q06swqga6VWtV7rXzHExkA0boG9fePAA3n5bnAW0bm2S+V1VVdlzeQ8FHQpSpXAV6pSow8BaAxlZdyRlCpR5+QUkFo0U4pBknOPHhVpScLDYeX3wgfAnLZO7bxj3ku8x9+BcjKqRr97+6nFX6SO5vhfx7zNcELuu9HqjPu/1EztUJX9ccaZPhy1HQ+H9JmCdAor4d1m3RF329dsn4gsKgk6d2OrpzIw3yxBXqitftuuRocRjqvGda/HX+G73dyw/tpw0QxptK7RljPcY6rnVy9D8afP/Nefmg5sc/+R4hmPINlRVzO86OgqP6fPn4dtvxY63Zk2TLKE36gmKCMI/xJ9D1w/Rp3oflrVfZpJrS3IXUohDknUSE+HePXBzE5KFUVHg7y/MEXLp7O4jLsddZtqBaSw6sohEXSJdK3d9fKbWvEzzdF0jK3O68N9dm2uePNTSV+frvs6cOAGFCkHVT+ZwwurJLrlfhfdYfLosyrz5MHgwm0vUZGenL/m19FvCHSeFDI3fZHV8R1VVHqQ9IK9tXgyqgVUnV9H3jb6MrDsy07KbPh4+fLnrS2KSYnBxcMnUNUzGo/PdwEAICxONg6tWQblyYrTIRCw7uozvdn9HZFwkZQuUZV7refSp3sdk15dYDjIBS17M1atihGXBAqhXTzi1lCkjDifN1OSTEX489CODfxuMoij0qNqDMV5jqFqkaoavY4rO2/Y1ilO3aHHmzYN5s+DAHShfL4Lvf4TR73uy5WI7um9YA6qK1qgw4Iv1KFEG0cQG+O26xLWyTytWZeQhILMPEQajgU1nNuEf4k9+u/z80esP3PO7c3P0zeee76aXR2Ndu6N206lSpyxdK0ssWgQTJojz3bJlxTFKH9MlxTuJdyjoUBArxYoL9y7g6ujK1OZTaVuh7UurL5JXl9x/B5WYhyNHhDxkqVJip9u4sdBpfkQuTb6qqvJX5F9E3IkAxKjL8DrDuTTsEsvbL89U8oX0z+k+jxMnRKXezTuU7/6ahOu7P+I9qx3nmlUi3Gk8dnbQpXIX9hr7MGGHkZ0/afB6ZwCcOQPz5wNZfwjI6OuTdEnMPTiXCrMr0PnnzsQkxdCuQjseHVtlNfkCvFnsTRxsHDI1D5xloqKeqLBdvSp+1zdvFu/5oEGi0pNFzt89z8dbP8Z9ujtbz20F4Bufbwj5MIQOnh1k8n3NkTtgyRMe3Yw0Gvj1V9i6VdijDR36XHec3ILBaCAoIgi/ED8OXT/EgJoDWNBmAVUKV2Fqi6kvff3LzkYz03mrqrB9O0ydKsaItGVDML7fGJRUTgD5EvLxTf3xDL5SROy83N3xatwHL+tS8PHH4PJ0STar4zcZff3cg3Px/dOXOsXrMKXpFNpXbG/yhGGjsaG+e32i46NNet0XEhoqysxBQUKpqkMHGD8evvnGdEtEh+If4s+mM5uw0djwfrX3qVRIiKCY2w9aknuQTVgSoVK1ZAnMmCEcibp0EbZpIFR9cjnLjy5nwp4JXIy9SLkC5fD19n1pR/M/SW+DVXobmFJSxNFhYCCcOgWuxVMZMsialFpT+OHAOFRUrFD4StOYr+dFiDGuH36AL74wSZyZff25u+cIDA2kSakmdKnchdjkWE7dOZXhxqqMkqpPzX5hCYPhyfnugQOib+Hjj8XDZbFipl3KaKDsrLLcT7nPoDcHMeStIbg6mmc8TGJ+ZBOW5NlcvvzEIi0+XpzxFiokvpbLE++95Hs42zmjKAqn7pzCxcEF/2b+mTpTS+/Z6D8FM55FTIywAJw9WygRVqkdR9cZP7JXN4NKrebg6tiIqYfsSNOloNWrNF+0Eyo1FYa9LVq8NM6sjt887/WFCkbRYe0QNp/ZjFajpZSTqHY42ztT371+uq6dFbI1+RoMT8aFPv9cKILNni3Od00kEZmiT2Hl8ZWsPL6S33v9jp21HZu6baJMgTKZMpqQvD7IHfDriqoKMYGzZ8WOd+RIeCv3D/xHxUURGBrI4vDFbOi6gXfKvkOaIQ0bK5tM79JKffYrz/pXoACRk1u/9PUXL4px6CVLILlgKKVb/YJHtSscjN9CQloCzUo3Y0LJvtRp2IPQ6FCCA4fik1gYr8GToHr1TMVsKvpu6svyY8spYF+AQbXFbu3filU5Qa+gXng4efB94+9Nc8FHD5ebNsHJk2BvL8583dxMMr8LQrZ03qF5zAybya3EW9RwrcHPXX6W87uSp5A7YIkYsdiwQWSJoCDRYLJokRCQd3Mzd3QvJfxGOP4h/qw7tQ4rxYqe1XpS2rk0kHXFqsyerYaFQUCAeDutraH5B6HsKN6ES4ZkLsVAs1JNmWLVnBpzN0JoTwivhNcbXngFHsxxbeZHJOuS+d/x/9G9Snfy2ualbYW21C5Wm35v9DNJU1VmiUmK4ejNo1lPwGFhosy8YYP4e9euorpjby9ciUzExXsXqT6/Oom6RFqUaYGvty+NSzWWHrySDCET8KtOXJxwwZk1C6KjxYhFZKTY/VqICbvBaKDtmrbcT7nPyLojGV53OCXymcZXFTLWYGU0it60gADhOJcvv0r3sftIqjyfN0pUYNvuNAA0KDRaf5gaW3ZA6dLi/S9bVlzEDDfpmKQY5vw9h9kHZxOTFIOdtR3vV3+fjp4dczyWZ+Hj4cPnOz/nduJtCucpnLmLhIcLHfL8+YVoxtChJn24DL8RzsnbJ+ldvTelnUszsu5IulTuQrUi1Uy2huT1QibgV5krV4T9XGKicCKaM0dI6OXSEaJH6I16fj71MytPrGRjt41oNVo2dN1A+YLlcbIzvehHes5WU1OFFkNAgJhScXM30s9vMyfz+/HTjQO4XHehTcVGaDVa0gxpaNMM+OhLwIZF0K6dycqeGUVn0DH89+EsO7qMZH0y75Z/F19vXxq4NzBLPM/j0Tzwnst76Fypc/pelJAgRLPj42HcOHjjDeEx3a4d5M1rkrhUVWX7xe34h/izM3InRR2L8l6V97DR2DCh8QSTrCF5fZEJ+FVCVYXZ/ZkzosPT3R18fcUN6Y3MCd/nJIlpiSwJX0LggUCi4qKo6FKRK/evULZA2WwXpH9eg9X9+/DjjzB9Oty4AeUahdL+k185rF/B0oRoSmlLMafa5/TdfJnUZbOZ1mMK0WmHqGosy60fPgAzGcRHxUXh4eSBjcaGC/cu0L1Kd0Z7j348CpPbqFW0Fnls8hAcFfzyBBwdLSoKCxaI/0HNmonffUURHtMmYv+V/Qz6bRDHbx2nWN5i+DX1Y2Ctgdhocq/BiMSykAn4VSAtTRgiTJsGR4+KstsHHwiz+6+/Nnd06eJS7CXeXPgm95LvUc+tHjPemcG75d/FSjHPbv3GDTGVNW+e2GA1bBFLu0nLWX71Cy7FiTLzd8V78/naa1jvmITe3oH1VZqReq84+W1Kc4WMyUSaAqNq5Ndzv+If4k/YtTCihkdRNG9Rfu/1u9nex/Rio7FhYK2BeDh5vPgbFy8WD5eqKty2Ro6EOnVMFkdCagIJaQkUy1uM/Hb5MapGlrZbSo+qPUzqjiWRgOyCtnz++AP69RMZo1IlGDFC7ALsc8AbNYtcuHeBE7dO0MGzA6qqMvKPkXSp1IV67vXMFtPZs6LMvGKF6Ftr9V40+d+ZzuarC0jVp2JUjRhUAxqsmLDDyOeXisGwYbyTVJEzqf99ni3uZM/+zxpna8wp+hR+Ov4TAaEBnIk5g3t+d0bWHcmAmgPM2lhlEgwGceheqpQwRoiIEM2Dw4ZByZImW+afHrwty7VkdafVAI81wyWSzCK7oF81Tp8W57gVKz65MS1dCs2bm627NiP8fe1v/Pb7ERQRhIuDC63Lt0ar0TL9nenZtubLRDQOHYLJk0VHs1YLnQdEkvzWN2y5vAo1UuW9ch1ofs7Ix8oW0hTRee3Tdyz0/hy0Ws5+9usz182IVnRmuZ5wnYFbB1KtSDV+6vgTXSp1sdgyaYo+hQdpD3BR7WHZMlH7v3ABPvpISHJ6egppMRNx+s5pAkICWHl85WMP3lF1Rz3+uky+kuxEJmBL4ZGu4bRpYtfbtSusXQvly8Pvv5s7unRx+PphRm8fze7Lu8lvm5/P6n/G0LeGZntp73kuQKoK+eOKM2kS7NgBDhVDqP/ln4zt3JzSxZx4a1EQg8v2YOQ+AyUnB0FyMuUGtSH4PS98PHzwcnvSRZ5VmciMcDnuMtMPTOdm4k1Wd1pNaefShH8UTtXCVS06YRhVI+7T3OmSWoY5gWchNlbMpq9dCx1N1639qOqnKArLji5jzck10oNXYhZkArYEVq+G778XO19XV+Ha8tFH5o4qXaQZ0khITaCgQ0EUReFS7CUCmwfSv2Z/8tqaplP1Zfxb6UpV4e4pF95f5siDq1CkqIEW309hu348ezFyaOsUdr6/kxtpQ3HsOVlsiXv3hpEj8apUiWcNb2VGKzqj/HMWWlEUelbticFoQGOlsfxRmJMnsapUiVrFahF8/hA0aiRGiby9TVbVeaQX7h/iz3eNvuOdsu8wtt5YPq33qfmtECWvJTIB51Zu3IDChcX4yunTYGsrDia7dhUf53Lup9xnweEFTA+bTtPSTVnefjk1i9YkcnhkjjvAPCoDqwaFxNPFiA8rg+5uXjQF79IzcAFhmgD+iD3/+PvTDGkERwXjVb8JjLcRzjhFXqwOlVWZyJex/Ohy+m7uS15tXkbUHcHwOsNxy5/7BVReiNEoqjeBgbBzJ2zZgk9JHz678Du3V8zL/Dzwv0jSJbHs6DICQwO5GHuRMs5l0Bv1ABR0KGiSNSSSzCCbsHIbR46IMvPatfDzz2KEKC0NbGws4nz3Wvw1ZoTN4MfDPxKfGk+TUk0YW28szco0M1tMdScEc36vC3FX7mAsdABN3Js4lyqI3ZszuGpcTq3Cb9AhzpUfkv4gTVHRamzY+cHup0rMOY3OoGPtqbW4OrrStHRT7ibdZXH4YgbWGpgts9A5ik4nzncDA8XIXDHRyMbAgfyddJ46i+qwtvNaulbumuWlVFWl1oJahN8Mp07xOvh6+2aLq5NE8jxe1ISFqqo59qdWrVqq5BkYDKq6caOqNmyoqqCqefKo6tChqnrpkrkjyzAjto1Qrb61Ut9b/556+Pphs8YSH6+qfn6q6lRQr1IiRGWcncrXqHxtrZYcF6guDTmi7vr8PdWY11FVQQ1p84Y6cUFvNeTyPvPFnBKvTg2ZqroFuql8g9pjQw+zxWJyUlPFf3U6VS1ZUlVr1FDVlSuffF5VVZ1BpzpOdFQ/2fpJppe5cPeCOvbPsWqqXlx3Y8RGdU/UHtVoNGYleokkUwCH1OfkRLkDNidGo+hmNhhEM5VeL+Tz+vcXdmm5HFVV2XN5D/4h/ozyGkXjUo25+eAmybpkSjmbzz84NlboNEyfLj6u2/YEUW/14Kb+5MPAoVelz/hf10nibNdgEOeNtZ/9kJpTzAybyVd/fcX91Pu8XfJtfL19aVmuZa6f4X0pp0+L3e6OHWLOy9ZWHLG4uj6zqrP6xGoqulSkRtEaGVrm4LWD+If4syFiAxpFw199/jLrSJtEAnIMKfdx+bLIEL/8AsePg52d6Gz28BCq/rkcg9HAxjMb8dvvx8HrBynkUIjbibcBssX3NL0+vLduier93LlCpbBtW7Bq+xGbri7ATrXDGitUoxGtAQa5VhEvWr7crNKcEXciKOlUEgcbBxxsHGhepjljvMdku/JXtqOqsGuXGBnatk23yPsEAAAgAElEQVT8jvftK2RRbW2haNHnvrR71e4ZWupu0l06rev0uLve19uXYXWGUSyvaX1+JRJTk/vv9q8KqgqhoSJDBAWJJ//OnYVZgqvrE6F+C8BnuQ/7ruyjbIGyzGs9jz7V+2Bvkz3CH88bIYInjU/XroG/P8zbEkpaiV1U7W3Fkv6jqF3DlqV/16R2bFs+WX6Ks3EXCX4jPz4N38frjTZiATMkX1VV2XdlH34hfmw9t5XZLWcz+K3B9K/Zn/41++d4PNlCaCg0bSqa1yZMEOpVLunrNNYZdHyz/Se2hqcSH1/imQ9daYY0jt86Tu1itSlgXwBHrSNTm09lQM0BOdZdL5FkFVmCzikOHYI33xSl5YEDYcgQi7ABBOGks+zoMkbUHYG1lTU/Hf8JO2u7HGlmqTd51zPna4s72bOqe2OmTBHqhPoSwdC7OUZFB8DXDcfzTaPv4N498T5XqgSjR4uHHjNVGVRVfTwGE3YtDBcHF4a8OYTBbw22/DGYe/eEUIbBAOPHiwfODRugTZsMd+1vOHyZrls8cTD4UFA3GBAjXZM6VqWRpyM/Hv6RGWEzSEhN4Oqoq+SzzZcdP5FEYhJkCdoc3LsnbABTU+Grr6BWLfjpJ9HVnMcy5AEv3rtIYGggS48uJVmfTK2itWhUqhE9q/XMsRiepSSli3Xg+LYylB0PaHTUGOrHsfwTSTOK5KuooJu3GHy+hQIF4MQJoRhmpi5yo2rESrFCURSmhk7lTtId5rSaQ983+uJg42CWmEzGhQuiqrNsGSQlCcGMR8YIndPpavQvAv+8iK2xMilWJx5/7oEuhsG/jiLh920kpCXQpFQTPq33KXm1crcrsVyylIAVRRkJ9AdU4ATQT1XVFFMEZrGcPStU/JcvFzektm2f3JB69DB3dOkiLiWOAVsGEBQRhLWVNb2q9jKbk84/FaZ0d/Nw/0BZEk8VA7sEhnwMvr7WNNqwhny3CnJfm4wRFa0RtLEe/HLgIm29ygo/XjNwN+kucw/OZVH4Ig4OOEjhPIVZ33U9RfIUeTXGYGbNguHDxYhcz57CGKFq1Sxf9npcMrbWVUm2OYSeGKxxwaDc44Z+Pe95dsXX25eaRWua4AeQSMxLpg/AFEUpDgwDaquqWgXQAO+ZKjCLZOZMoc+8ZAm8955osNq82SLmd42qkfN3hRhFPtt8XIu/xqfenxI1PIrF7RZnKfluCr9Gvcm7KPXZr9SbvItN4dfS/VrfFhWwup+PmyEGrp8JIdFmB9r+nXD4oiTf+cXi7q7Q5dC73JkSzaY1eWge+QYlkr5haa0vmLL7SqZjzgqRsZEM/W0o7tPd+Sr4K6oUrkJ8ajwAxfIWs9zkq9eL+fSTD7vJGzWCL78UTYVLlpgk+YJ46LIziGvdsv0CAK1ahtq2a1ndabVMvpJXhqyWoK0Be0VRdIADcD3rIVkQycmirFyrFtSoAU2awDffwCefCBUrCyBVn8qqE6sICA3gRsINroy8gqPWkf0f7DeJrnB6mqiex+nTsNavOJH7o6BPV9CkggJY2fIRb2Jcsxr6DWJj0bfQtRrJVs8GpFo/0ZXOCSOEf3PzwU3Kzy6PgkLPaj0Z4zWGyoUr53gcJiU+XjgQzZgBV64Ix61p06BKFfHHxPi2qMBnQckk6Jtio5YAxBnwl++YTxhFIskOMp2AVVW9pihKAHAFSAa2q6q6/d/fpyjKQGAggLu7e2aXy13cuCFmXebPh5gYYXpfowZUriz+WAD3U+4/bma5nnCdakWqMbPlTGw1omHGVKL+/9ZhBkjWGfD/4+xzE/CpU6Jxdt06cVz+5ogNHLROFXGpMHafke+274MPykM/cHHJz4aqTf5znewwQvg3qqqy/eJ2/r72N+PfHo+roysL2yykWelmFM+XMz7A2cp334lRovh4ePttmD0bWrfO1iWfyHp+ni2ynhJJbiHTXdCKojgDG4BuQBzwM7BeVdWVz3vNK9EFPWKESL56Pbz7rjj38vGxiDIzPPE3DbsaRt3FdWlauim+3r40K90sW5x0Sn32K8/6DVOAyMlP38hPnhT3+583JWP71jJqNr7CluGTOJcUis+SBhgMBrRG2Hm7JV6DJ0H16sB/d9nwpGs2u27aOoOONSfXEBAawPFbx3HL58bpwadx1Dpmy3o5ytGj4r1VFPjiC4iKyhVCJRKJJZJdXdBNgUhVVe88XCQI8Aaem4AtEoNB2AC2aCFmRgsWFE5Ew4ZBuXLmji7dHL91HP8Qf/Jq8zK39VzqlKhDxOAIKrpUzNZ102PTd+oUDPMLZVf0b1g738Thi80kKXewcqhI/rQheLl5EVzJn+Bjm/HpNAqvGm2fulZ2GyH8m/1X9tN9Q3ei46OpXKgyy9oto3vV7tluq5itGI2wZYvY7e7dK0wSWrSAH36wmIdLicTSyEoCvgLUVRTFAVGCbgJY+Pb2H8THC5P7mTPh0qUnN6Tx480dWbpRVZWdkTvxD/Fn+8Xt5LHJw5C3hjz+enYnX3ixTV9EBHz7LawNCYE+jaBUGnoFKuvdmLnDlQYHzqCoP8Gnn+LVZSReXUY+d532NYpna4ny5oOb3E26S+XClSlboCwVXSoyr/U8y5eKTEsTg9TTpsH58+DuLmQjvR6et8rkK5FkG1k5Aw5TFGU9cATQA+HAAlMFZjYSEsTc7uLF4mNvb5gyRTRYWRgT905k3F/jcHV0ZWLjiXxc+2Oc7Z0zdI30ykA+j2ftTntUqMzPAUVYtesYWlsDDg2DSNLoQQErI3QLjqZyUi2U9bOhffsMxWtqzsacJSAkgBXHV1CneB329NtDEccibO/9n3YHy0KnE+NDIHa5xYrBmjXQqZNFyKFKJK8CUgkLxJzu9etQvLg4261QAerWFee9b75p7ujSTUJqAouOLMLbzZs6JepwKfYSf0X+Ra9qvbC1zriHsKnPVs+fh+8mqPwUsgulvh/GUtspkloLhfe4pf0SRU3D2mhF3ZufoC/Tkf2fNc7wGqbi0PVD/LD3Bzaf2YyttS393ujHKK9RlC1gOZKhz+TkSbHb3b1btJlrtS80RpBIJFlDKmE9j7Q0Mdc4fbpIwJcvixvS6dMWYXr/iBsJN5gZNpP5h+cTlxLH5/U/p06JOpR2Lk1p58yLUGSmg/lZREbCkCmhbLs9H7XEAeh9jkIaZ0ZeLMvANYfp3GsAFPqBVOU4tmo1Lhf0RDHDCJFRNWJUjVhbWRMaHcqey3sY33A8g98abDJzeLOgqsKJaOpUYfphby+MEZKSxO/7C4wRJBJJ9vF6JuCYGDFCNGcO3LwpxDO++UbcqMCiku+nf37KjLAZ6I16Onp2ZIzXGOqUqGOSaz9vjja987VXr8LX3yexbPthjD1agGsKiqIy9nRBvg66i11he+b69CcmjxO2xhLY4vn4tTkxQvSIVH0qK4+vJCA0gFF1RzGg1gD61+zPBzU+II/WMmRDX8iePdC8udjlfv+9MEYoWNDcUUkkrz2vVwI2GECjgWPHRDNVixai0ap5c7Na0mUEVVXZH70frxJeaKw0uDi40L9Gf0Z5jaJMgTImXSs9HczP4uZN+GryXZacmIOh1izKdPMhSpuGQVWxMkA+1Qa7JSugWzeKnbqDLugEPKNJK7uJS4njx0NiFvrGgxvUcK3xeHY3u9ydcoTYWPGAaWUFY8dCw4ZiqLptW4t6uJRIXnVe/QRsNIoO5unThWpPYCA0biw0m8uXN3d06UZv1BMUEURASAAHrx9kQ9cNdPTsyKf1Ps22NV/Uwfws7t6FL/2jWHQqEEO1xdAwiZYPStL58kmGVNaSZkhDq7XBJ2ADuHsDOT9C9E86rO1AcFQwzUo3Y0WHFTQp1SRbZqFzjIsXxe/5kiWivNyli/i8ojz5WCKR5Bpe3QT84AGsWCHk886dEw1WbR56wCqKxSTfNEMaPx76kWkHphEZF0m5AuWY33o+Lcu2zPa105sc79+HkYGh/LQ/mDTPpWhqRtLrdlE+35BEpfs3oHdvPN+bR/D1UHw8fPBy8/rPOjmRcE/dPsX0A9OZ3HQyBR0KMrHxROys7ahRtEa2r53tzJghRGGsrYXpx6hRUK2auaOSSCQv4NXtgv7oI1iwQHQxjxwprNEejV1YAGmGNLQaLQajgQqzK1A4T2F8vX1pW6FtrhHzT0xUGT5jJ8siJ2AochCs09AqVqxbraNdTEEYNAgGDxam7GZCVVX2XtmL334/fj3/Kw42DgR1DaJF2RZmi8kkGAywcaOo6lSsKNSr1q0TPtPFipk7OolE8pAXdUG/GglYVeHAAVF+GzsWatYUu947d8QcrwWVFc/EnCEwNJBtF7ZxdshZHGwciEmKyVWG7UkpegbPXc/KS37oC4Vjq9Ois9FhREWjaJjg2IbPP/4JHMzrdZusS6bxisYcuHoAFwcXhr01jEFvDqKggwU3ID14IPoWpk0T7eWjRonuZolEkit5dceQdDpYv14k3r//BicnISRQs6YoMVtImVlVVfZd2UdAaAC/nP0FO2s7+lTvQ7IuGQcbB7Mn39DoUIKjgmng7sPxnRUZfroW+ryRFLNx5utttlS4kUrLvhrSNKDVaPHp8qnZkm+KPoV9V/bRtHRT7G3sqeFag/ervU+fN/rgYGPeB4Is88MPEBAAcXFCqSogANq1M3dUEokkk1huAjYahQPRqVMi0c6ZA336CPscC+PIjSM0XNaQgvYF+arhV7lq7jQ0OpTGKxqTqk8DvS3qsp00ruPKkFOXaXv+PpouXQke2puyly8SnXQIN4fa3IpxB7ecjTM2OZZ5h+YxI2wGMUkxRA2Pwi2/G3Nbz83ZQExNRIQoMSuK6HJr2hRGjxZCMRKJxKKx7BL0okXivOuddyxmjAggSZfEsqPLiE2O5cuGX6KqKutOraNNhTa5apd26V4k7y7tQUTCAWFfZNTQs+gE/le0Msqe3TB8OJvuanLcieif3E68zZR9U1hwZAEP0h7QokwLxtYbi4+Hj+V2NKuqMACZOhX+/BN27hSd+6pqUccpEonkdTgDthBuJ95m9t+zmXtwLneT79LIoxE73t+R68T8L8Ve4qOfv2DH9Z9RjGClGAHQWtmw84PdT3Ux15u865mzwsWd7LNVSjJVn4qttS3X4q9Rfnb5xyIk1V2rZ9ua2Y5OBytXilG5kyeFQtWwYaKh0DljGt4SiSR38OqeAVsQK4+vpP8v/UkzpNG2Qlt8vX3xdvPOFbu00OhQ/or6i9rFauOW1pyR38I+j80MPaLlswMpRNapyJ721fFpPfg/I0RZVcvKCKqqsufyHvxC/EjVp7Lj/R0Uz1eca6Ou4WTnZPL1coxHAjFGI3z+ORQuDMuWQffuQipSIpG8ksgEnE08UqxytnOmcuHK1C5Wmz7V+zDKaxQVXLJf5Sm9TNwRxLj93VBVPahWKMv24hjrTZhbUzw9rLDZOppiDRpQ7zkPCplVy8oIBqOBzWc347ffj7BrYRRyKMSwOsMwqkasFCvLTb4XL4pu5r/+EupstrYQFiYsAXPBg5lEIsleclft8xXAYDSw/vR6vBZ70WBpA/xC/ADhvftjmx9zTfJNTEuk/4bvGL+3DyoPrQAxUqTOUmZvukG1Yxux+XWzkDF8QTLwbVEBe5un55JNLSW58MhCOq3rRExSDPNaz+PyiMuMazgu15Xu001IiOjWL1dOzKq/9ZawvgQoWVImX4nkNUHugE3IsqPLmLBnApdiL1HGuQxzWs2hT/U+5g7rmQTumcHik19T7rYtl13AoIAVGhzc3Pnx7wjeb5w+h5zskJKMTY5l/qH5lCtYjs6VOtOzak8K2heko2fHXCNCkml27wYfH3Gm+9lnUjhDInmNkQk4i9xOvI2LgwtWihXn756ncJ7C+Dfzp12FdrkmWYRGh7IhYgORsZF0r9yT2NCOrJ3QkV18TZFbJZnQrDY7qilY8wYGG88Mn9+aSkryavxVpoVOe9zRPPjNwXSu1Jm8tnnpUtlCtYwTE4VwhtEoGqoaNIDFi6FrV3B0NHd0EonEjMgu6EzySLFqxbEV/NzlZ9pUaIPOoMNGk7vkLpeGL2XAlgEYVAOo0PaoO79svoyXF7gWX8uR0nn+U/LM7g7mZzFh9wQm7JmAUTXyXpX38PX2teyO5hs3YNYs4UoUGwutWsGvv5o7KolEksPILmgT8UhXOCAkgC3ntmBnbUffN/pSqVAlgFyXfAduGcjCIwtBBRTQGCHfgwJsXpdKm862bD5an4gMuB2ZkkdNalULVyW/XX4qulTk49ofM8prFB5OHtm+frYyf77Y7er10L49jBkjJFElEonkH1hoF4t5MKgG3t/4PqFXQ/n67a+5MuIK89+db3If3syiM+hYfWI1iWmJANQ8ZcfgMLDTg5VBQaOxY+CiubTtYouiiNLxpI5VKe5kj4LY+Wa3gIZRNbLpzCbqLalHg6UNWHRkEQBdKndhZsuZlpl8VRV27BBdzSCkUAcMEHrkQUEy+UokkmciS9Av4EHaA5aGL2X1ydXs6rMLO2s7Ttw6QZkCZXKNYlVodCjbL24nLiWOjRFBXI6/wiynQZyMmMOahQn0t1lOylhPCr39N83L/dcKMKdQVZUl4UvwD/Hn7N2zeDh5MMZrDP1q9Ms172WG0elg7VqhyXzsGAwfLnTJJRKJ5CGyBJ1BbiTcYNbfs5h3aB5xKXF4u3lz88FNPJw8qFqkqsnW2RR+LUvdw3uv7KXx8sbojXoAqt/RMPtPuHfpAYuBTwbnZez4IRQqBNDEZHFnhEe2ioqisObUGuxt7FndaTWdK3XG2sqCf/1mzYIpU+DaNahUSTRW9ehh7qgkEokFYcF3wOzhbMxZqs2vhs6go4NnB0Z7jcbbzfQlxE3h157SUL4Wl8znQScAXpqE7yXfo4B9AfZd3vc4+VoZofqpivxwbiElOnsRMQnKln2Y5BebbkQovdxIuMGMsBksDl/MkYFHcMvvxrrO63Cyc8oV6l+Z4vp1IQ+pKMIkoUIFWLgQWrSwKC1yiUSSO3jtE7CqqvwV9ReXYi/Rv2Z/yhcsz9dvf03Xyl0pW6Bstq3r/8fZp5qfAJJ1Bvz/OPvcBHnw2kH8Q/z59cwvXPrgGG+X9MFW1aJTdagGW8IdFrIg1OuxUU5WknxmORtzloCQAFYcX4HeqKdzpc6PHxKc7S1Uz/jwYWGMsG4dBAdD/fowcyZYv/b/fCQSSRZ4be8gOoOOn0//TEBIAOE3wyntXJp+b/RDY6XhiwZfZPv66dFQ3hR+jfG/BRGZshnF+jIPuED+NCuGhRm5ff9nxh0ZR+r5YFxqBzO6sw9jJ3g9NVGUmSSfFWKSYqg6ryoaKw0f1viQ0V6jc02DWoYxGuG330TiDQ6GvHlh5EgoXVp8XSZfiUSSRV7Lu8gfF/5gwJYBRMdHU9GlIgvbLKRXtV45KpzxMg3lTeHXGBG0jsuasaDRgQpDwmB0ZHl+d/mCWlO7kd8FZn/jxcCBXtg8YwIqu40SVFVl24Vt7L28l0lNJ+Hi4MJPHX/ibY+3c42fcaZJTYV+/cDeXiTh/v0hXz5zRyWRSF4hXpsEfDX+KnqjHg8nD4rnK06ZAmWY23ourcq1MoumsG+LCs/00R3SpATTD0xn+m+/c18tAhiETrMRdjh0Zu6F1dhEWTNqrDDOyZ//+Wtkl1GCzqBj7am1+O3348TtE7jlc2Ns/bE42TlZrmJVTAzMmwd//CHkIu3tYdcuqFiRZz7dSCQSSRZ55RPw0ZtHmRo6lTUn19ClUhdWdVpFlcJV+KvPX2aN61EJ+KttQUQnHcLVvjxVSt1i8PZlxOoTaBgNd/OO5X4Ba1SjAaNBy5m/R5Gnwi1O/locD4+Xr/G8JJ8VoY1D1w/RaV0nrty/QuVClVnefjnvVXkPrcZCbfPOnxeORMuWQXIytGwJ9+5BoUJQ1XQd7xKJRPJvXtkEvPPSTibtm8TOyJ04ah0Z8uYQhtcdbu6wnqKIyxUuqGNJtUklTm/kzDnoGAG+J/NxpkhbJudtgtWuphicD2KTVIOCjQ2U9jyLh0f6zm9NZZQQkxTD9YTrVCtSjbIFyuLp4smcVnPMVj0wGSEhoqHKxgZ69YJRo6ByZXNHJZFIXhNeqQScqk99PHP6+4XfiYiJYErTKQysNTDXecb+fe1vVhxfQZohDSNGFBVGRjgx1Wcip337MHO4FWf/ssPaOZFC5Z2wL3cLB+0DfFtkbFeWFaOEqLgopoZMZXH4Yiq6VOTwwMM42Tnxe6/fM3U9s2MwwObNEB8PfftCnTowaRL06QOuruaOTiKRvGa8EkpY95LvMf/QfGb9PYul7ZbyTtl3iE+Nx87aLleVRh81LfntmsDumweok+jM8fwpQqxCsWZ9251smV+PBQtEv0+HD+KIcA7n5oOkHJ3hPXX7FBP3TWTtybVYKVb0qtaLMd5jHmteWxyJiaLEPG2akIusUwdCQ6XvrkQiyXZeWSWsS7GXmH5gOovDF5OkS6J5meYUsC8AQD7b3NGxGhodSnBUMFaKFSv/XsjJhIuUuA+BYQr9PZpw8uOh7Liyn1sHfOje0IvERBg8GL76ClxcnIBGORKnqqoYVSMaKw1h18L45ewvjKg7ghF1R1AiX4kciSFbWL1aeO7euycS7+TJ0KGDTL4SicTsWGwCNqpGGi9vzPWE6/So2oNRXqOoVqSaucN6ih2XdtB2dVvSDGlYqeB2z8CKA3a8V/8jbFaPQnVz53oQLPVtSGQkvPsu+PuLxtucwmA0sOnMJvxC/OhepTsj6o6gV7VedKjYwXKFMyIixNxuiRLg5gYNGz5xJJKJVyKR5BIsNgFbKVas6LCCsgXKUixvMXOH8xQ3H9xk1v5AAg9MJxU9KiooGj50bUnv7ashf37Cw2FEb9izRzTb/vknNG2aczGm6lNZcWwFAaEBnLt7jtLOpXF1FOegWo0WrX3uKd2nC1UVb2ZAAGzdKna9s2aJJqv69c0dnUQikfyHLCVgRVGcgEVAFYTr7AeqqoaaIrD00LBkw0y/NqtGCM/i3N1zTN31PctPrSINAz6REOJhhV5jhVajpVGv8dxKyc+Xo2HJEihYUFjH9u8PmpzTAAGgR1APgiKCqFm0Jms7r6WTZ6ccFSIxKRs2iNLyoUNifOjbb+GTT8wdlUQikbyQrO6AZwC/q6raWVEULWARvnKm0EjeFH7t8Qyvm0NtvmvZkSlbmhKui6bfURhl35hyg78itJQNwZd3413chz2rvGjxA6SkiImXcePAKYeas28k3GBm2EyG1x2Oq6Mrvt6+fFL7E5qUamKZ5ghJSeDw8Nft999FZ/OPP0Lv3kJEQyKRSHI5me6CVhQlH3AMKK2m8yK5xQ+43uRdz1SIKu5kz/7PGr/09ZvCrzE8aA1XNJ8BehRscDdMYZLdHRpfv0aR4V9C+fKAqIxu3Ai+vnDpErRpI5QNy5Uz9U/1bM7fPY9/iD/Ljy1Hb9Szov0KelbrmTOLZwfXrwsjhPnzhWpVnToi+To6SkciiUSS68iuLujSwB1gqaIo1YHDwHBVVRP/tfhAYCCAu7t7FpYzHVnRSE4zpDF2yzRijPPAWrj8oOq4rx5ltqYP3ec8SeDHj8OIEfDXX0LfYft2aNbMJD/CSzEYDfQM6sm6U+vQarT0e6MfY7zHZKvDU7Zy8qR4cvnpJzHP26nTE21mqdEskUgskKwkYGugJjBUVdUwRVFmAJ8B4//5TaqqLgAWgNgBZ2G9x2T1/DazGskJyfep7FeSaO5T+gFE5we9lQLYYGes+jiB370rxojmzxcl5tmz4aOPst9AR1VVjt48So2iNdBYachnm4/P6n/G8DrDKeJYJHsXz05SUkQnc2oqfPyxeKp55EokkUgkFkpWUsJV4KqqqmEP/74ekYCzFVOc36ZHI/nRGe/lpL24aBSmtg+kfY3i9L3sTMELLhwp0ovtDnlI0ZzEzlgVW6MnRfM6MHu2SL7x8TBokOgHKlDAhG/AM9Ab9Ww4vYEp+6cQfjOcU4NOUalQJRa0WZC9C2cXOh38/DNs2gRr1oCdHQQFQbVq2f9mSiQSSQ6R6QSsqupNRVGiFUWpoKrqWaAJcNp0oT0bU3jcvkwjeVP4NQYFzeamxg/V2kg8MG5lCWAU3805zaYz95gfdAJbnQFbvVCHMl4txLWwGgy9AI0bw4wZUKWKyX7sZ5KiT2Fp+FICQgO4FHuJCgUrsKjNIso4W6gHb3w8LFoE06dDdLQYiL5+Xczz+viYOzqJRCIxKVktig4FfnrYAX0J6Jf1kF6MqTxun6eRfCnmPGNXteVGnjPiEwooRrinufowyTd+KoFfjoLkfVW4d6owHh5iIia7hZZUVUVRFNYdOs+QP0ZibSxJBe23TPTuR8eabtm3cHZy7JgoM8fHi2Q7b55wJpKNVRKJ5BUlSwlYVdWjwDO7u7KL7PC4NapG7ibdpVCeQtjfiSNBPUOXU3Zs9NSht1JBsUa1fvOpJN+8QnGObCyO3zIxw/v99zB6tKiWZhfXE64zLXQax24d45PKS5j06zVc9bOxVl1JSVX4cuMprBSrHNGLNgnHj8Ply6I1vHJl4UjUrx/UztFfKYlEIjELFre98G1RAXubpwUjMutxuztqN52Xt6b0N850/k7Ui4t6vkmrpOkcKrUWF/1knPS9KJL2A7ZGT4o52aOqsH49eHrChAnQsSOcPQtffpl9yfdszFn6/9KfUjNKEXggEBcHF6b8cZxknQEbtSgKYrv9qBSfq1FVIfvVogVUry4GolVVdKjNmSOTr0QieW2wOClKU3jcxqfG88W6j5hzcY0oMSvQN6E4aloailbLuwM6sz/oBEadJ7ZGT0Ak+W5lK9O0KezaJfqB/vc/UTV9HqZQ2/r9wu+0+qkVtta29K/Rn9HeoyntXJpSn/36zO/PaCk+R9m1SyTcY8eE/d/EiaI93BKFQCQSiSSLWFwChuLYyBEAABGxSURBVKx53AIsWzSEOTFrHv/dykqDbffeKFrt4+vDkyRf2M6RIudrMWqSI/nypW+sKLPd2qqq8uelP0nVp9KmQhveLvk23/p8y0e1P6JwnsKPvy87SvHZQny86GouWFD8XacTOpw9eoCtrXljk0gkEjNicSXozHDm+nH6BzRk5cpPAfig0/cs1nTC3sYejaJBq9Hi4+Hz1Gva1yjO3k8b81WF1kTNe5utqx358EM4d07YBb5spvdF3drPQm/Us/bkWmotqEWLlS3wC/EDwN7GnvFvj38q+YJpS/HZwtWrQv7LzU0ckAM0agQnTohzXpl8JRLJa45F7oDTy4LdgUwPnkwEd7DTQbnLD6AXOLq688G49Xg+9Or18fDBy83rqdcePiwMdQ4cAC8v+O03qFUr/WtnpFt785nNjN4+mouxF6lQsAKL2y6mZ9UXy0WaohSfLRw/LhyJVq8WZ7tduwp9ZhClZllulkgkEuAVTsBtJ1RiizECVLBWFTZUm0CrLl889T1ebl7/Sbz37omGqh9/hMKFYdkykT8yOg3zshJxXEocCgr57fJjVI0UdCiIfzN/2lVsh5WSvsWyWoo3Gar6JLFOnSpEM4YMgeHDwcPDrKFJJBJJbuWVKUGn6lJYuvZzYmOiAXC2d0ZRAQVUjRXHilq9cPdlNAoNiPLlYeFCGDZMdDf36ZO5UdTnlYg/fNuJsX+OxX2aO4GhgQC0r9ieAx8eoINnh3Qn31yBTgcrV0KNGnDkiPjcpElCRGPaNJl8JRKJ5AVY7A449GH5uHaRGoTv+B8zon/mup0O3dXLDBy9io+7BfDziiakGdKeecb7Tw4dEue6f/8tvNvnzBFdzlnh3yVi53x3cSn2B5/sXI/eqKdr5a508OwAYHl2gPHx4ill+nRx1uvpKT4HUKyYeWOTSCQSCyHTdoSZwVR2hKHRoTRZ0YQUfTLqw11usxsOfOrZnyb9J6LkyfP4+553xgvCNOGLL0QuKVIE/P2hZ8/sOabssLYD285ve+xKVKaAhcpFGgzCCOHKFaFY5esL77wjFaskEonkGWSXHaHZCI4KJs2QxsPcy2CX1swa/8tTSUDM4CZzPa4aW52S8W1x7fGu1GgUkzBjx8L9+8Jc55tvTOdqp6oqOyN34rffjzmt5lCuYDmmNp/KvNbzcHV0Nc0iOcnx48IU4YcfhOzXpEmiVi9FMyQSiSTTWGQC9vHwQavRivKytZYebb/8T/J93gxuSYozaJDobm7QQJSbq1Y1TVwGo4GNZzYyed9kDt84jKujK5diL1GuYDlKO1uYfZ6qws6doqP5jz/AwUEciFeoIGZ4JRKJRJIlLDIBe7l5sfP9nc8tLz9rBjfxgcIng43cDgMXF1ixQkgPm6rcbDAaqPFjDU7cPkG5AuVY8O4CelfvjZ11NopDZxeRkUJj8+hRUZv/4QfhwyutACUSicRkWGQChmePED3in7O2qgqJp4sT+1dFjEm2DB4kdCGcnLIeQ3xqPL+c/YVe1XqhsdLwfvX38XDyoEPFDmisNC+/QG4iIUGojNSqBcWLC+WqRYvEoXh2OkxIJBLJa4rFJuAX8WgGN+2OI/f+rEJqdEG0RWOp2O8Es2e9meXr3068zYwDM5hzcA73U+9TvUh1qhapyhjvMSaIPoe5fh1mzoT588HREaKiQKuFHTvMHZlEIpG80rySCXhow4oM9k3h3gEPrLR6CrQ4TqHa1/m2U9YOe+8m3eXr4K9ZHL6YVH0qHT07MrbeWKoWMdEhck5y/rxoplq5UnQ2d+woOppfprEpkUgkEpPwyt1tt2yBz4YW495lKFzrOrbep3ArpsG3RdVMq0Y9SHuAo9YRO2s7giKC6FGlB5/W+5QKLrlEdzm9qKoQz9Bq4eJF0dk8cCCMHAllLHQsSiKRSCyUVyYBX7kilA83bRLe7nv2QIMGxYDMC0Psu7KPyfsmczH2Iic/OUkebR4uDruIvU0ucxx6GXo9bNggOpobNQI/P+HHGx39xKVIIpFIJDmKxasn6HQir3h6immZyZOFKmKDBpm7nlE1svXcVuovqU+DpQ0IuxZGz6o90Rl1AJaVfBMTYdYsMbP73ntCrerRzJWiyOQrkUgkZsSid8D794vpmJMnoU0b0UuUVfnhree20m5NO9zzuzOr5Sw+qPEBDjYOJok3xxk+HBYvBm9vCAyEtm2lYpVEIpHkEixSihJg4kThWuTmJjZ57dpl7jpJuiSWhi/FRmPDwFoD0Rv1BEUE0aFiB2w0NiaJNcc4d064EQ0ZIna6Z84Ieydvb3NHJpFIJK8lL5KitNjt0CMZ4tOnM5d8Y5Nj+WHPD3hM92DItiH8dv43AKytrOlauatlJd/9+6F9e6hYEZYvf+JMVLGiTL4SiUSSS7HYErS3d+Zzy7Kjyxi6bSgP0h7QsmxLPq//OfXd65s2wJxAVUUz1Z9/CpWqcePE7rdwYXNHJpFIJJKXYLEJOKOcv3seR60jRfMWpbRzadqUb8PYemOp7lrd3KFljORk2LwZunUTjVRNmoiz3X794KELlEQikUhyPxZ7Bpxejtw4wuR9k1l/ej1D3hrCzJYzc3R9kxETA3PnwuzZcOcOhISA17OlOCUSiUSSO3jl7AjTw+6o3UzcN5HtF7eTzzYfY+uNZXjd4eYOK+PExsL48cI/MTkZWrUSh99165o7MolEIpFkgVcqAauqivLQ3mjZsWUcvXmUSU0m8UntT8hvl9/M0WWQe/fEua69vZD36tYNxowRKiMSiUQisXheiRK0zqBj1YlVTNk/hRUdVlC7WG1ikmLIY5PHsoQzjEbYtg38/eHyZaHXbG0Nqalga2vu6CQSiUSSQV7JMSQQM7wzw2ZSZmYZ+m7ui43GhmSdsCJ0cXCxnOSbmgpLl4rZ3XffhQsX+H97dx8jVXXGcfz7sEIaKPWlUEoRChaiaSPiBl+QlxJbjVACLUYjAUsCCfEdE00FTXyNGmlatQ2KtiWlLbZaW6jBl0i02pgoAQEBgxWoELdSoLRCsTEo+/SPc6mT8d5l2B3m3DP+PslmZu85s/s8nLn32Xvmcg9XXx1uIQkqviIiTSjZKehD7YcY/vBwtv17G2MGjWHRpEVMGDrh/1PQSVm5EmbNCgV4yZJw28gePWJHJSIix1CyBbilWwu3ffM2hpw4JL3/w9vWBg88AH36wLx54cKqF14ICyWk+AeEiIgctWQLMMDlZ1weO4Sjs3FjWDniscfCTTTmzAnbu3WD88+PG5uIiDRU0p8BJ+Wee2D4cHjySbjqqvA570MPxY5KREQiSfoMuNQOr8Hb2grDhsEFF4SrnK+8UssAiohI18+AzazFzNaZ2Yp6BJS86jV4Fy8O2886K9yrWcVXRESozxT0XGBzHX5O+u69FwYNguuug/79YdkyuPvu2FGJiEgJdakAm9nJwHeAn9cnnATt2BEuqIJwdfPYsfDKK58sEdhNH7OLiMindbU6PAD8AGgv6mBmc8xsjZmt2bNnTxd/XYmsWgUXXwxDhoRiC2HqeflyGD06bmwiIlJ6nS7AZjYJ2O3ur3fUz90fdfeR7j6yb9++nf115dDeDitWwLhxYTGEF1+E+fPDRVags10REalZV66CHg1MNrOJwOeAL5jZb9x9Rn1CK6EPPwzr7vbsCfffD7NnQ+/esaMSEZEEdboAu/t8YD6AmY0Hbmy64rtvHzzyCDz9dDjb7dkzPJ52GnTvHjs6ERFJmOZM87S1hTV3Bw6Em24KxXbv3tB2+ukqviIi0mV1uRGHu78EvFSPnxXd6tVw3nnhyuZLLw1r8La2xo5KRESajO6E5Q4vvwy7doVF71tb4eabw2e9gwfHjk5ERJrUZ3cK+tChcF/mc84JqxDddVcoxi0tcMcdKr4iInJMfTYL8HPPwamnwiWXwPvvw6JFYepZSwGKiEiDfHamoPfuDQsk9OsHvXqFtXgXLIApU8JZr4iISAM1/xnwO+/AtdeGezTfeWfYNmYMvPoqTJ2q4isiIlE07xnwunXhDPeJJ0KRnT49rMMLmmoWEZHomqsAu39SXBcuDDfQuOEGmDsXBgyIG5uIiEiF5piC/ugjWLoUzjwzLJIAYRnAd98NZ8EqviIiUjJpF+ADB+DBB2HoUJgxAw4ehA8+CG39+sHxx8eNT0REpEC6U9Dt7TBiBGzbFtbgXbgQJk7UikQiIpKEdAtwt27h5hmDB8OoUbGjEREROSrpFmCAadNiRyAiItIpmq8VERGJQAVYREQkAhVgERGRCFSARUREIlABFhERiUAFWEREJAIVYBERkQhUgEVERCJQARYREYlABVhERCQCFWAREZEIVIBFREQiUAEWERGJwNy9cb/MbA+wo44/sg/wzzr+vJiUS/k0Sx6gXMqqWXJpljyg/rl81d375jU0tADXm5mtcfeRseOoB+VSPs2SByiXsmqWXJolD2hsLpqCFhERiUAFWEREJILUC/CjsQOoI+VSPs2SByiXsmqWXJolD2hgLkl/BiwiIpKq1M+ARUREkqQCLCIiEkESBdjMLjKzv5rZVjObl9NuZvaTrH2DmbXGiPNIzGygmf3ZzDab2ZtmNjenz3gz22dm67OvW2PEWgsz225mG7M41+S0l35czOzUin/r9Wa238yur+pT2jExs8VmttvMNlVsO8nMVprZluzxxILXdrhfNVpBLj80s7ey988yMzuh4LUdvhcbrSCX283s7xXvo4kFry3NuBTk8XhFDtvNbH3Ba8s2JrnH36j7i7uX+gtoAbYBpwA9gDeAr1f1mQg8CxhwLrAqdtwFufQHWrPnvYG3c3IZD6yIHWuN+WwH+nTQnsS4VMTbAvyD8B/nkxgTYBzQCmyq2LYAmJc9nwfcV5Brh/tVSXK5EDgue35fXi5ZW4fvxZLkcjtw4xFeV6pxycujqv1HwK2JjEnu8Tfm/pLCGfDZwFZ3/5u7HwR+B0yp6jMF+JUHrwEnmFn/Rgd6JO6+093XZs//A2wGBsSN6phKYlwqfAvY5u71vFvbMeXufwH+VbV5CrAke74E+G7OS2vZrxoqLxd3f97dP86+fQ04ueGBdULBuNSiVOPSUR5mZsClwG8bGlQndXD8jba/pFCABwDvVnzfxqeLVi19SsXMBgNnAqtymkeZ2Rtm9qyZfaOhgR0dB543s9fNbE5Oe2rjchnFB5NUxgSgn7vvhHDQAb6U0ye1sQGYRZhRyXOk92JZXJNNpy8umOpMaVzGArvcfUtBe2nHpOr4G21/SaEAW8626v87VUuf0jCzzwN/AK539/1VzWsJU6BnAD8Fljc6vqMw2t1bgQnA1WY2rqo9mXExsx7AZOD3Oc0pjUmtkhkbADO7BfgYWFrQ5UjvxTJ4GPgaMALYSZi+rZbSuEyj47PfUo7JEY6/hS/L2dblcUmhALcBAyu+Pxl4rxN9SsHMuhMGf6m7/7G63d33u/uB7PkzQHcz69PgMGvi7u9lj7uBZYRpmkrJjAvhILHW3XdVN6Q0Jpldh6f6s8fdOX2SGRszmwlMAqZ79oFctRrei9G5+y53P+Tu7cDPyI8xiXExs+OAqcDjRX3KOCYFx99o+0sKBXg1MMzMhmRnKZcBT1X1eQr4fnbV7bnAvsNTCmWSfWbyC2Czu/+4oM+Xs36Y2dmEMdrbuChrY2a9zKz34eeEi2U2VXVLYlwyhX/NpzImFZ4CZmbPZwJ/yulTy34VnZldBNwETHb3/xb0qeW9GF3V9Q/fIz/GJMYF+Dbwlru35TWWcUw6OP7G219iX5lWyxfhatq3CVeh3ZJtuwK4IntuwMKsfSMwMnbMBXmMIUxbbADWZ18Tq3K5BniTcJXda8B5seMuyOWULMY3snhTHpeehIJ6fMW2JMaE8EfDTuAjwl/ps4EvAi8AW7LHk7K+XwGeqXjtp/arEuaylfDZ2+H9ZVF1LkXvxRLm8utsP9hAOHj3L/u45OWRbf/l4f2jom/Zx6To+Bttf9GtKEVERCJIYQpaRESk6agAi4iIRKACLCIiEoEKsIiISAQqwCIiIhGoAIuIiESgAiwiIhLB/wDuhq8H4yiwogAAAABJRU5ErkJggg==\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:                Mon, 24 Feb 2020   Prob (F-statistic):           1.66e-29\n",
      "Time:                        22:49:06   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.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
