{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heterogenous Wage Effects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.011184,
     "end_time": "2021-02-28T17:17:25.882205",
     "exception": false,
     "start_time": "2021-02-28T17:17:25.871021",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Application: Heterogeneous Effect of Gender on Wage Using Double Lasso\n",
    "\n",
    " We use US census data from the year 2012 to analyse the effect of gender and interaction effects of other variables with gender on wage jointly. The dependent variable is the logarithm of the wage, the target variable is *female* (in combination with other variables). All other variables denote some other socio-economic characteristics, e.g. marital status, education, and experience.  For a detailed description of the variables we refer to the help page.\n",
    "\n",
    "\n",
    "\n",
    "This analysis allows a closer look how discrimination according to gender is related to other socio-economic variables.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hdmpy\n",
    "import pyreadr\n",
    "import os\n",
    "from urllib.request import urlopen\n",
    "import patsy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import figure\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import scipy.linalg as sci_lag\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>lnw</th>\n",
       "      <th>female</th>\n",
       "      <th>widowed</th>\n",
       "      <th>divorced</th>\n",
       "      <th>separated</th>\n",
       "      <th>nevermarried</th>\n",
       "      <th>hsd08</th>\n",
       "      <th>hsd911</th>\n",
       "      <th>hsg</th>\n",
       "      <th>cg</th>\n",
       "      <th>ad</th>\n",
       "      <th>mw</th>\n",
       "      <th>so</th>\n",
       "      <th>we</th>\n",
       "      <th>exp1</th>\n",
       "      <th>exp2</th>\n",
       "      <th>exp3</th>\n",
       "      <th>exp4</th>\n",
       "      <th>weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>29217.0</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "      <td>29217.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>2.797007</td>\n",
       "      <td>0.428757</td>\n",
       "      <td>0.007975</td>\n",
       "      <td>0.113393</td>\n",
       "      <td>0.016600</td>\n",
       "      <td>0.156347</td>\n",
       "      <td>0.004107</td>\n",
       "      <td>0.022179</td>\n",
       "      <td>0.247288</td>\n",
       "      <td>0.283431</td>\n",
       "      <td>0.155800</td>\n",
       "      <td>0.291645</td>\n",
       "      <td>0.282849</td>\n",
       "      <td>0.199644</td>\n",
       "      <td>18.756939</td>\n",
       "      <td>4.286811</td>\n",
       "      <td>10.875998</td>\n",
       "      <td>29.408779</td>\n",
       "      <td>1513.842566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.662406</td>\n",
       "      <td>0.494907</td>\n",
       "      <td>0.088947</td>\n",
       "      <td>0.317078</td>\n",
       "      <td>0.127769</td>\n",
       "      <td>0.363191</td>\n",
       "      <td>0.063957</td>\n",
       "      <td>0.147267</td>\n",
       "      <td>0.431443</td>\n",
       "      <td>0.450671</td>\n",
       "      <td>0.362672</td>\n",
       "      <td>0.454528</td>\n",
       "      <td>0.450391</td>\n",
       "      <td>0.399740</td>\n",
       "      <td>8.767040</td>\n",
       "      <td>3.321506</td>\n",
       "      <td>11.121864</td>\n",
       "      <td>36.569919</td>\n",
       "      <td>1009.811610</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>-7.469874</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>106.790000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>2.408296</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>11.500000</td>\n",
       "      <td>1.322500</td>\n",
       "      <td>1.520875</td>\n",
       "      <td>1.749006</td>\n",
       "      <td>654.240000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>2.774540</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>19.000000</td>\n",
       "      <td>3.610000</td>\n",
       "      <td>6.859000</td>\n",
       "      <td>13.032100</td>\n",
       "      <td>1472.100000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>3.181569</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>26.000000</td>\n",
       "      <td>6.760000</td>\n",
       "      <td>17.576000</td>\n",
       "      <td>45.697600</td>\n",
       "      <td>1966.630000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>2012.0</td>\n",
       "      <td>5.970942</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>43.500000</td>\n",
       "      <td>18.922500</td>\n",
       "      <td>82.312875</td>\n",
       "      <td>358.061006</td>\n",
       "      <td>6444.150000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          year           lnw        female       widowed      divorced  \\\n",
       "count  29217.0  29217.000000  29217.000000  29217.000000  29217.000000   \n",
       "mean    2012.0      2.797007      0.428757      0.007975      0.113393   \n",
       "std        0.0      0.662406      0.494907      0.088947      0.317078   \n",
       "min     2012.0     -7.469874      0.000000      0.000000      0.000000   \n",
       "25%     2012.0      2.408296      0.000000      0.000000      0.000000   \n",
       "50%     2012.0      2.774540      0.000000      0.000000      0.000000   \n",
       "75%     2012.0      3.181569      1.000000      0.000000      0.000000   \n",
       "max     2012.0      5.970942      1.000000      1.000000      1.000000   \n",
       "\n",
       "          separated  nevermarried         hsd08        hsd911           hsg  \\\n",
       "count  29217.000000  29217.000000  29217.000000  29217.000000  29217.000000   \n",
       "mean       0.016600      0.156347      0.004107      0.022179      0.247288   \n",
       "std        0.127769      0.363191      0.063957      0.147267      0.431443   \n",
       "min        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "25%        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "50%        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "75%        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "max        1.000000      1.000000      1.000000      1.000000      1.000000   \n",
       "\n",
       "                 cg            ad            mw            so            we  \\\n",
       "count  29217.000000  29217.000000  29217.000000  29217.000000  29217.000000   \n",
       "mean       0.283431      0.155800      0.291645      0.282849      0.199644   \n",
       "std        0.450671      0.362672      0.454528      0.450391      0.399740   \n",
       "min        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "25%        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "50%        0.000000      0.000000      0.000000      0.000000      0.000000   \n",
       "75%        1.000000      0.000000      1.000000      1.000000      0.000000   \n",
       "max        1.000000      1.000000      1.000000      1.000000      1.000000   \n",
       "\n",
       "               exp1          exp2          exp3          exp4        weight  \n",
       "count  29217.000000  29217.000000  29217.000000  29217.000000  29217.000000  \n",
       "mean      18.756939      4.286811     10.875998     29.408779   1513.842566  \n",
       "std        8.767040      3.321506     11.121864     36.569919   1009.811610  \n",
       "min        0.000000      0.000000      0.000000      0.000000    106.790000  \n",
       "25%       11.500000      1.322500      1.520875      1.749006    654.240000  \n",
       "50%       19.000000      3.610000      6.859000     13.032100   1472.100000  \n",
       "75%       26.000000      6.760000     17.576000     45.697600   1966.630000  \n",
       "max       43.500000     18.922500     82.312875    358.061006   6444.150000  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "link=\"https://raw.githubusercontent.com/d2cml-ai/14.388_py/main/data/cps2012.RData\"\n",
    "response = urlopen(link)\n",
    "content = response.read()\n",
    "fhandle = open( 'cps2012.Rdata', 'wb')\n",
    "fhandle.write(content)\n",
    "fhandle.close()\n",
    "result = pyreadr.read_r(\"cps2012.Rdata\")\n",
    "os.remove(\"cps2012.Rdata\")\n",
    "\n",
    "# Extracting the data frame from rdata_read\n",
    "cps2012 = result[ 'data' ]\n",
    "cps2012.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "136"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_basic =  '''lnw ~ -1 + female + female:(widowed + divorced + separated + nevermarried +\n",
    "hsd08 + hsd911 + hsg + cg + ad + mw + so + we + exp1 + exp2 + exp3) + +(widowed +\n",
    "divorced + separated + nevermarried + hsd08 + hsd911 + hsg + cg + ad + mw + so +\n",
    "we + exp1 + exp2 + exp3) ** 2'''\n",
    "\n",
    "y, X = patsy.dmatrices(formula_basic, cps2012, return_type='dataframe')\n",
    "X.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the same number of covariables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "variance_cols = X.var().to_numpy()\n",
    "X = X.iloc[ : ,  np.where( variance_cols != 0   )[0] ]\n",
    "\n",
    "def demean(x):\n",
    "    dif = x - np.mean( x )\n",
    "    return dif \n",
    "\n",
    "X = X.apply( demean, axis = 0 )\n",
    "\n",
    "index_gender = np.where( X.columns.str.contains('female'))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter estimates for the target parameters, i.e. all coefficients related to gender (i.e. by interaction with other variables) are calculated and summarized by the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "effect_female = hdmpy.rlassoEffects( x = X , y = y , index = index_gender )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                     Estimate.  Std. Error   t value      Pr(>|t|)\n",
      "female               -0.154593    0.048599 -3.180992  1.467717e-03\n",
      "female:widowed        0.136096    0.094115  1.446056  1.481614e-01\n",
      "female:divorced       0.136939    0.021680  6.316402  2.677227e-10\n",
      "female:separated      0.023303    0.048849  0.477034  6.333381e-01\n",
      "female:nevermarried   0.186853    0.020094  9.299104  1.416347e-20\n",
      "female:hsd08          0.027810    0.140357  0.198140  8.429356e-01\n",
      "female:hsd911        -0.119335    0.052223 -2.285101  2.230692e-02\n",
      "female:hsg           -0.012885    0.018229 -0.706850  4.796595e-01\n",
      "female:cg             0.010139    0.018079  0.560797  5.749360e-01\n",
      "female:ad            -0.030464    0.022807 -1.335714  1.816429e-01\n",
      "female:mw            -0.001063    0.018715 -0.056822  9.546869e-01\n",
      "female:so            -0.008183    0.019027 -0.430096  6.671256e-01\n",
      "female:we            -0.004226    0.021422 -0.197281  8.436076e-01\n",
      "female:exp1           0.004733    0.007483  0.632488  5.270682e-01\n",
      "female:exp2          -0.159519    0.043479 -3.668908  2.435892e-04\n",
      "female:exp3           0.026512    0.008244  3.215830  1.300681e-03\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'\\\\begin{tabular}{lrrrr}\\n\\\\toprule\\n{} &  Estimate. &  Std. Error &  t value &  Pr(>|t|) \\\\\\\\\\n\\\\midrule\\nfemale              &     -0.155 &       0.049 &   -3.181 &     0.001 \\\\\\\\\\nfemale:widowed      &      0.136 &       0.094 &    1.446 &     0.148 \\\\\\\\\\nfemale:divorced     &      0.137 &       0.022 &    6.316 &     0.000 \\\\\\\\\\nfemale:separated    &      0.023 &       0.049 &    0.477 &     0.633 \\\\\\\\\\nfemale:nevermarried &      0.187 &       0.020 &    9.299 &     0.000 \\\\\\\\\\nfemale:hsd08        &      0.028 &       0.140 &    0.198 &     0.843 \\\\\\\\\\nfemale:hsd911       &     -0.119 &       0.052 &   -2.285 &     0.022 \\\\\\\\\\nfemale:hsg          &     -0.013 &       0.018 &   -0.707 &     0.480 \\\\\\\\\\nfemale:cg           &      0.010 &       0.018 &    0.561 &     0.575 \\\\\\\\\\nfemale:ad           &     -0.030 &       0.023 &   -1.336 &     0.182 \\\\\\\\\\nfemale:mw           &     -0.001 &       0.019 &   -0.057 &     0.955 \\\\\\\\\\nfemale:so           &     -0.008 &       0.019 &   -0.430 &     0.667 \\\\\\\\\\nfemale:we           &     -0.004 &       0.021 &   -0.197 &     0.844 \\\\\\\\\\nfemale:exp1         &      0.005 &       0.007 &    0.632 &     0.527 \\\\\\\\\\nfemale:exp2         &     -0.160 &       0.043 &   -3.669 &     0.000 \\\\\\\\\\nfemale:exp3         &      0.027 &       0.008 &    3.216 &     0.001 \\\\\\\\\\n\\\\bottomrule\\n\\\\end{tabular}\\n'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_coeff = pd.concat( [ effect_female.res.get( 'coefficients' ).rename(columns = { 0 : \"Estimate.\" }) , \\\n",
    "             effect_female.res.get( 'se' ).rename( columns = { 0 : \"Std. Error\" } ) , \\\n",
    "             effect_female.res.get( 't' ).rename( columns = { 0 : \"t value\" } ) , \\\n",
    "             effect_female.res.get( 'pval' ).rename( columns = { 0 : \"Pr(>|t|)\" } ) ] ,\\\n",
    "             axis = 1 )\n",
    "\n",
    "print( result_coeff )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_value = stats.t.ppf(1-0.05, 29217 - 116 -1 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "pointwise_CI = pd.DataFrame({ '5%' : result_coeff.iloc[ : , 0 ] \\\n",
    "                                     - result_coeff.iloc[ : , 1 ] * t_value ,\\\n",
    "                              '95%' : result_coeff.iloc[ : , 0 ] \\\n",
    "                             + result_coeff.iloc[ : , 1 ] * t_value })\n",
    "pointwise_CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAH4CAYAAADNU5vyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAABHVklEQVR4nO3de5heZ13v//eXtDMjlgyMvxZKUyjQEqUCobQBuqFaiopy2kAF0cpBUFRUtkEOce/9Q06OqFROF8i5QDeHDfyoVVSk0gKlQFp6pGoAkZZEyimQaYWZ0OH7+2OtJ51M5vBMsjL3uifv13XNlTxrnnQ+13TNmvVd931/78hMJEmSJKlLtysdQJIkSdLaY6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6d0TpACsxOjqaRx99dOkYkiRJkoCdO3fuyczRhT5XVaFx9NFHs2PHjtIxJEmSJAER8a3FPufUKUmSJEmds9CQJEmS1DkLDUmSJEmds9CQJEmS1DkLDUmSJEmd67TQiIiTIuKyiPhiRFweEScv8J6HRMTV7cf1EfGmiFiwJZYkSZKkOnU9ovEm4M2ZeW/glcB5C7znGuC0zNwE3Bc4BvjdjnNIkiRJKqizQiMijgFOBc5vD30IOD4iTpz7vsz8fmb+sH05AvwYkF3lkCRJklRelyMaxwNfz8xbATIzgRuBu81/Y0ScEBHXAN8GdgNv6DCHJEmSpMKKLAbPzK9m5v2BuwCjwBMWel9EbImIHYOPW265ZVVzSpIkSTowXRYaXwOOjYgjACIiaEYzblzsH2TmLcD7gF9b5PPnZuaGwcdRRx3VYVxJkiRJh0pnhUZmfhO4EjinPfREYEdmfnnu+yLixIg4sv37CPB44NquckiSJEkqr+upU88Gnh0RXwReBDwDICLeGhGPbd/zcOCqdo3GVcA3gJd1nEOSJElSQdGs2a7Dhg0bcseOHaVjSJIkaY6ZmRlmZmb2Oz46OsroqNulrWURsTMzNyz0OXcGlyRJ0kGZnJxkfHx8v4/JycnS0VSQIxqSJEk6KIMRjc2bNwOwbds2wBGNw8FSIxpHrHYYSZIkrS2DgmLdunUArF+/vnAi9YFTpyRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUucsNCRJkiR1zkJDkiRJUueOKB1AkiRJt5mZmWFmZma/46Ojo4yOjhZIJB0YRzQkSZJ6ZHJykvHx8f0+JicnS0eTViQys3SGoW3YsCF37NhROoYkSdIhMxjR2Lx5MwDbtm0D6hjROPnkkwG4/vrrCyfRaomInZm5YaHPOXVKkiSpRwYFxbp16wBYv3594UTSgXHqlCRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOHVE6gCRJkqTlzczMMDMzs9/x0dFRRkdHCyRamiMakiRJUgUmJycZHx/f72NycrJ0tAU5oiFJkiRVYOvWrWzZsoXNmzcDsG3bNoBejmaAhYYkSZJUhcEUqXXr1gGwfv36womW5tQpSZIkSZ2z0JAkSZLUOQsNSZIkSZ2z0JAkSZLUOQsNSZIkSZ3rtNCIiJMi4rKI+GJEXB4RJy/wnodHxLaI+JeIuD4i/jwiLHgkSZKkNaTrG/w3AW/OzHsDrwTOW+A93wV+JTPvAzwQOB14asc5JEmSJBXUWaEREccApwLnt4c+BBwfESfOfV9mXpWZX2n/Pg1cDZzQVQ5JkiRJ5XU5onE88PXMvBUgMxO4EbjbYv8gIu4CnA38XYc5JEmSJBVWbG1ERKwH/hb488y8YpH3bImIHYOPW265ZXVDSpIkSTogXRYaXwOOjYgjACIiaEYzbpz/xoi4A/CPwN9k5rmL/Qcz89zM3DD4OOqoozqMK0mSJOlQ6azQyMxvAlcC57SHngjsyMwvz31fRBxFU2T8Y2a+vKuvL0mSJKk/up469Wzg2RHxReBFwDMAIuKtEfHY9j3PBTYDT4iIq9uP/9lxDkmSJEkFHdHlfywztwMPWeD4s+b8/RXAK7r8upIkSZL6xY3yJEmSJHXOQkOSJEmqxPT0NLOzs8zOzjI1NcX09HTpSIvqdOqUJEmSpENjenqa4447jl27dgEwPj7OxMQEO3fuZGxsrHC6/TmiIUmSJFVgz549e4uMgV27drFnz55CiZZmoSFJkiSpcxYakiRJkjpnoSFJkiRVYGRkhImJiX2OTUxMMDIyUijR0lwMLkmSJFVgbGyMnTt3smnTJgC2bdvGyMhILxeCg4WGJEmSVI2xsTHWrVsHwPr16wunWZpTpyRJkiR1zkJDkiRJUuecOiVJkqTDzszMDDMzM/sdHx0dZXR0tECitccRDUmSJB12JicnGR8f3+9jcnKydLQ1wxENSZIkHXa2bt3Kli1b2Lx5M9B0cAIczeiQhYYkSZIOO4MpUrV0cKqRU6ckSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLn3LBPkiRJqsDMzAwzMzPMzs4CMDU1Bdy2+WDfOKIhSZIkVWBycpLx8XG2b9/O9u3bGR8fZ3x8nMnJydLRFuSIhiRJklSBrVu3smXLlv2O93E0Ayw0JEmSpCr0dYrUYpw6JUmSJKlzFhqSJEmSOmehIUmSJKlzFhqSJEmSOmehIUmSJKlzFhqSJEk6aNPT08zOzjI7O8vU1BTT09OlIy2rxsw1sb2tJEmSDsr09DTHHXccu3btAmB8fJyJiQl27tzJ2NhY4XQLqzFzbRzRkCRJ0kHZs2fP3hv2gV27drFnz55CiZZXY+baWGhIkiRJ6pyFhiRJkqTOWWhIkiTpoIyMjDAxMbHPsYmJCUZGRgolWl6NmWvjYnBJkiQdlLGxMXbu3MmmTZsA2LZtGyMjI71eVF1j5tpYaEiSJOmgjY2NsW7dOgDWr19fOM1wasxcE6dOSZIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkznVaaETESRFxWUR8MSIuj4iTF3jPCRFxSUTsjoiru/z6kiRJkvrhiI7/e28C3pyZ50XE2cB5wGnz3jMF/C9gHHhFx19fkiQdIjMzM8zMzOx3fHR0lNHR0QKJJPVZZyMaEXEMcCpwfnvoQ8DxEXHi3Pdl5q7MvBT4r66+tiRJOvQmJycZHx/f72NycrJ0NGnFZmZmmJqaYnZ2ltnZWaamppiamlqwmNaB6XLq1PHA1zPzVoDMTOBG4G4dfg1JklTI1q1b2b17Nxs3bmTjxo3s3r2b3bt3s3Xr1tLRpBUbFM7bt29n+/btFs6HQNdTpzoVEVuALYPX4+PjBdNIknR4G0yRWrduHQDr168vnEg6cFu3bmXLli37HXcaYHe6LDS+BhwbEUdk5q0RETSjGTce6H8wM88Fzh283rBhQx58TEmSJB3uXFt06HU2dSozvwlcCZzTHnoisCMzv9zV15AkSZJUh6730Xg28OyI+CLwIuAZABHx1oh4bPv320fEDuADwH0iYkdEOBlOkiRJWkM6XaORmduBhyxw/Flz/v59YEOXX1eSJElSv7gzuCRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6pyFhiRJkqTOWWhIkiRJ6twRpQNIknQ4mpmZYWZmZr/jo6OjjI6OFki0Nvl9lspxREOSpAImJycZHx/f72NycrJ0tDXF77NUTmRm6QxD27BhQ+7YsaN0DEmSDtrgSfvmzZsB2LZtG1DHk/aTTz4ZgOuvv75wkuX5fV5dNWbWwYmInZm5YaHPOXVKkqQCBje669atA2D9+vWFE61Nfp+lcpw6JUmSJKlzjmhI0iHiIlRJB2p6eprZ2VkApqamGBkZYWxsrHAqaWUc0ZCkQ6TGRagzMzNMTU3t97FQwSTVYnDTPjs7y9TUFNPT06UjLWl6eprjjjuO7du3s337dsbHxznuuON6n1uaz0JDkg6RrVu3snv3bjZu3MjGjRvZvXs3u3fvZuvWraWjLarG4khaSo037Xv27GHXrl37HNu1axd79uwplEg6ME6dkqRDpMZFqFu3bmXLli0LduiRoL4pPUvdtPc5t7QWWGhIkvaqsTjS6hmMDgxu3MfHx5mYmGDnzp3etEvaj4WGpCq4sFoqz9GB1TEyMsLExMQ+3+uJiQlGRkYKppJWzjUakqrg2gFJB2Jw0z5X32/ax8bG2Llz5z7ruxw1Uo0c0ZBUBdcOSDoQg5v2TZs2Ac21o+/rSqDJ7RRG1c5CQ1IVXDsglVfrlB5v2qUyLDSkw5DrHbSWeD6vnlpHBySV4RoN6TDkegetJZ7Pq2swOrBu3TrWr19vkSFpUY5oSIch1ztoLan5fK5tTwpJWgkLDekw5HoHrSW1ns/uSbE6BlPr5hZ04NQ6aTU4dUqSpAKW2pNC3RlMrdu+fTvbt293ap20ihzRkCRJa9Zgat18jmZIh56FhnSQ7HgjledaBy3Ga7FUjlOnpINkxxuprMFah7lTY4477jimp6dLR1tSjTtWS9JKdFpoRMRJEXFZRHwxIi6PiJMXed8zI+JLEfHvEfGWiDiyyxzSatq6dSu7d+9m48aNbNy4kd27d7N79262bt1aOpp0WKh1rcNgT4q51w4XgktaS7oe0XgT8ObMvDfwSuC8+W+IiHsALwMeBpwI3Bn4rY5zSKtmdHSU9evX79NXfv369Q7Vq1qDaUizs7NMTU31fmSgZu5JIWkt66zQiIhjgFOB89tDHwKOj4gT5731bODCzLwpMxP4a+ApXeWQJB24WqchSZL6p8sRjeOBr2fmrQBtEXEjcLd577sbcMOc119d4D2StCbUNjpQ4zQk1zpIUj/1uutURGwB9vakGx8fL5hGklbGDdlWx2Ctw6ZNm4BmZ3C7TklSeV2OaHwNODYijgCIiKAZqbhx3vtuBO4+5/UJC7wHgMw8NzM3DD6OOuqoDuNK0qFV4+hArVzrIEn901mhkZnfBK4EzmkPPRHYkZlfnvfWDwGPjYi7tMXIbwPv6yqHpLWrtmlINXIakiSpK11PnXo2cF5E/DEwBTwDICLeSrMA/MLM/EpEvBj4dPtvLqHpViVJi3Ia0uqocRrSYNPMuRv2gRu1SVJpnba3zcztmfmQzLx3Zp6amde1x5+VmRfOed9bMvNe7cczM/OHXeaQtPbUOA2p1tGB2qYhDTbNnNspy00zJam8Xi8Gl3ToDKYhQfMEuO9PrWtU4+hAjbZu3cqWLVv2O+5ohiSVZaEhHYachrR6BqMDAOvXry+cZm1yipQk9VPXO4NLqoDTkCRJ0qHmiIakKjgNSZKkulhoSKqG05AkSaqHU6ekw5DTkCRJXZqZmWFqamqfvY6mpqaYmZkpHU0FWWhIh6HBNKSNGzeyceNGdu/e7UJwaZV5Y6a1xDbTWoiFhnSYqm2vBGmt8cZMa8nWrVvZvXv3fh9bt24tHU0FuUZDkqQC3P9Da4ltprUQCw1JkgrwxkzSWufUKUmSJEmds9CQOjA9Pb3Pgs7p6enSkSRJkopy6pR0kKanpznuuOP27rQ9Pj7OxMSEXZwkSdJhzREN6SDt2bNnb5ExsGvXLvbs2VMokSRJUnmOaEiS9pqZmWFmZobZ2VkApqamABcuS5JWzhEN9cpgA6v5H25gJa0O93aQJHXFQkO9MrjJmf/R55uckZERJiYm9jk2MTHByMhIoUTSgXPTLUlSV5w6pV4ZbGC1efNmALZt2wb0ewOrsbExdu7cyaZNm4Am88jIiAvBVSWnSEmSumKhoV4Z3OSsW7cOgPXr1xdONJyxsbHqMtemxrUDNWaWJKkrTp2SVIUa1w7UmFmSpK44oiGpCoNpdfP1eWSgxsySJHXFQkM6DNU4pafP2RZTY2ZJkrri1CnpMOSUHkmSdKg5oiEdhpzSI0mSDjULDekw5JQeSZJ0qDl1SpIkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc7F4JIkaSg17sEjqRxHNCRJ0lDcg0fSSjiiIUmShuIePJJWwkJDkiQNxSlSklbCqVOSJEmSOmehIUmSJKlzFhqSJEmSOucaDekg2e5RkiRpf45oSAfJdo+SJEn7c0RDOki2e5QkSdqfhYZ0kJwiJUmStD+nTkmSJEnqnCMaa9hgkfJ8PoGXJEnSoeaIxho2WKQ8/8NFypIkSTrUHNFYwwaLlDdv3gzAtm3bABcpS5Ik6dCz0FjDBlOk1q1bB8D69esLJ5IkSdLhwqlTkiRJkjpnoSFJkiSpcxYakiRJkjpnoSFJkiSpcy4GH5J7UkiSJEnDc0RjSO5JIUmSJA3PEY0huSeFJEmSNDwLjSG5J4UkSZI0PKdOSZIkSepcJ4VGRNwuIl4XEf8eEV+OiN9b4r2vjYivRkRGxKYuvr4kSZKkfulqROMc4D7AvYHNwPMj4uRF3vtB4KHADR19bUmSJEk901Wh8WTgLZk5m5m7gPcDT1nojZn5yczc0dHXlSRJktRDXRUad2PfEYqvtscOSkRsiYgdg49bbrnlYP+TkiRJklbBUF2nIuIzwEmLfPoB3cXZV2aeC5w7eL1hw4Y8VF9LkiRJUneGKjQy8yFLfT4ibgTuDnymPXQCcONBJZMkSZJUra6mTn0A+M2IWBcREzRrNt7f0X9bkiRJUmW6KjTeDfwb8CXgcuDczLwOICIeGxFvHbwxIt4UETuADcBHI+LLHWWQJEmS1BOd7AyembPAcxb53IXAhXNeP7uLrylJkiSpv9wZXJIkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ70zPT3N7Owss7OzTE1NMT09XTqSJEmSVqiTrlNSV6anpznuuOPYtWsXAOPj40xMTLBz507GxsYKp5MkSdKwHNFQr+zZs2dvkTGwa9cu9uzZUyiRJEmSDoSFhiRJkqTOWWhIkiRJ6pyFhnplZGSEiYmJfY5NTEwwMjJSKJEkSZIOhIvB1StjY2Ps3LmTTZs2AbBt2zZGRkZcCC5JklQZCw31ztjYGOvWrQNg/fr1hdNIkiTpQDh1SpIkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CY42bnp5mdnaW2dlZpqammJ6eLh1JkiRJh4EjSgfQoTM9Pc1xxx3Hrl27ABgfH2diYoKdO3cyNjZWOJ0kSZLWMkc01rA9e/bsLTIGdu3axZ49ewolkiRJ0uHCQmMFnIYkSZIkDcepU0NyGpIkSZI0PEc0hlTjNKSRkREmJib2OTYxMcHIyEihRJIkSTpcOKKxho2NjbFz5042bdoEwLZt2xgZGXEERpIkSYechcYaNzY2xrp16wBYv3594TSSJEk6XDh1akhOQ5IkSZKG54jGkJyGJEmSJA3PQmMFnIYkSZIkDcepU5IkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXMWGpIkSZI6Z6EhSZIkqXOdFBoRcbuIeF1E/HtEfDkifm+R941FxAUR8cWIuCYiPhYRJ3aRQZIkSVJ/dDWicQ5wH+DewGbg+RFx8iLvfTOwMTPvD/wN8NaOMkiSJEnqia4KjScDb8nM2czcBbwfeMr8N2XmdGb+fWZme+izwAkdZZAkSZLUE10VGncDbpjz+qvtseU8l2ZUQ5IkSdIacsQwb4qIzwAnLfLpBxzIF46IPwZOBM5a4j1bgC2D1+Pj4wfypSRJkiStsqEKjcx8yFKfj4gbgbsDn2kPnQDcuMT7/wh4AvCIzPz+El/3XODcwesNGzbkYu+VJElaC2ZmZpiZmWF2dhaAqakpAEZHRxkdHS0ZTVqRrqZOfQD4zYhYFxETNGs23r/QG9tRiqcAP5eZ3+vo62uNmJmZYWpqitnZWWZnZ5mammJqaoqZmZnS0SRJWhWTk5OMj4+zfft2tm/fzvj4OOPj40xOTpaOJq1IV4XGu4F/A74EXA6cm5nXAUTEYyPire3fNwCvAu4IXBwRV0fE5zrKoDXAi6sk6XC3detWdu/evd/H1q1bS0eTViRuawDVfxs2bMgdO3YUzXDyyU3X3uuvv75ojpWoKfNguHg+h4slSZL6JyJ2ZuaGhT431BoNabVYUEiSJK0NXU2dkiRJkqS9LDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnLDQkSZIkdc5CQ5IkSVLnjigdoBYzMzPMzMwwOzsLwNTUFACjo6OMjo6WjCZJkiT1jiMaQ5qcnGR8fJzt27ezfft2xsfHGR8fZ3JysnQ0SZIkqXciM0tnGNqGDRtyx44dRb72YERjvj6PaAwyb968GYBt27YB/c4sSZKkekTEzszcsNDnHNEY0ujoKOvXr9/vo8837I7CSJIkqRRHNNawGkdhJEmSVI+lRjRcDL6GWVBIkiSpFKdOSZIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSeqchYYkSZKkzlloSJIkSepcZGbpDEOLiBngW6VzAEcBt5QOsUJmXh1mXh1mXh1mXj015jbz6jDz6jDzgTs6M0cX+kRVhUZfRMSOzNxQOsdKmHl1mHl1mHl1mHn11JjbzKvDzKvDzIeGU6ckSZIkdc5CQ5IkSVLnLDQOzLmlAxwAM68OM68OM68OM6+eGnObeXWYeXWY+RBwjYYkSZKkzjmiIUmSJKlzFhqSJEmSOmehIUmSJKlzFhqSJEmSOndE6QA6tCJiNDNnSudYayLijKU+n5mfXK0sK7VI9u8BX8zM6VWOsyIRcVeAzPzP0lmGERFHAM8FTszM34mIewF3z8yPF462qIh4K3AR8PHM/GbpPJKketl1ahntTdkJwEVzb24i4mmZ+c5iwZYREfcD3gPcMTM3RMQDgSdn5gsKR1uRiHh0Zv5d6RzzRcTl7V/XAZuArwAJ3Au4OjNPKRRtWRFxNXBf9s28HRgHzsnMi8ulW1hE/BTwQeCu7aEdwC9n5r+VS7W8iPhrmnPkoZn5UxFxR5pryallky0uIn4d+Dng4cAumqLjosz8+6LBVigifjsz/7p0jsVExP+7wOHvAZ/JzMsX+FwxEfEOmmvFgjLzN1YxztAi4nFAZuaFEfFQ4JeBazPzbYWjLSoiHgl8PTOviYizgJ8FvpCZ7y+bbGGLnMd7ZeZLVyvLsCLiqUt9PjPftVpZhhURP2Lpn8F1qxhnRZw6tYSI2AK8jebidE1EPGHOp59bJtXQXgv8NvCt9vWVwKPKxTlgbygdYCGZeVpmngZcDfxCZp6YmScBP0/zve6zzwNnZeZJmXlvmhvKzwKPA/6iaLLFvQF4RWbeKTPvBLwCeGPhTMN4cGb+JjANkJnfA44smmgZmfnuzHwqcDzN+fB44G/Lpjogf1w6wDJ+Cvgdmu/zBprr9c8C/yci/qBgroVcQXPd2AM8mOYhxb8Dm4FejphHxMtozoGXRMSraK4Z/wk8LSL+d9Fwi4iIvwD+HHhvRLwIeA0wCvxRRPxp0XCLu0P78VPA7wF3ozmnnwP8ZMFcS3lM+/FrwJuBpwK/DrwJ+NWCuZZyB5qHgS8G/gy4e/sxCSxZ7JXmiMYSIuJamieRUxFxH+AC4KWZeX5EXJWZDyibcHERcUVmnjo3Z18zR8RiG84E8BuZOb6aeVYiIq7NzPvNO3ZNZt6/VKblLJRvcKyv2SPi6szctNyxvomIz2bmgwc/exGxjmbE676lsy0mIp4FPAI4BbgO+GfgY5n5paLBFhAR/99inwJ+PjN/fDXzrERE/CPwtMz8Rvv6zsC7gV8BPpWZJ5fMt5CI+CTw6Mycal+vB/4uM5ecSlpCRHyBZrT59sBNwPGZ+Z2IuAPw6fnX7T6IiH+hyfzjNKO2d8/Mb0fEjwPb+nhODETEPwFPH8z8iIhjgfMy8xfKJltcRHwY+N+Z+YX29ck093hPLJtscRHx+cx84HLH+sQ1GssYXFAz818i4uHAx9qbhb5XaLdGxJG0OSPieGC2bKRF/S7NU5yF8vX9+zwbEWcOphtFxM8APyqcaTk/iogzButI2umBg8x9/X7PRsR9MvNfANrCv6/n81zXRsQ5wO0i4kTghcAlZSMt6w3ANmAL8M+Z+YPCeZbyC8D/oHnSPlcAD1v1NCuzYVBkAGTmNyLirpm5KyJ+WDLYEo4e/E6E5vdjRBxdMtAS9mTmrcBURHw5M78DkJk3R0Rfrx0zmbkH2BMR38vMbwNk5n9FxPxzvG/uOnd6eWZ+PSKOKxloCCcOigyAzLw+Ik4qGWgId4iIYwbr5yLiGJrRjt6y0Fja7Nz/oZm5o50zeRHN0GCfvZ5mBOboiHg5cA7Q1/UZXwA+kJnXzf9E+3S1z54DvG/OjcERwJML5hnG3MxBk/lXIuIo4K+KJlvcHwOfbEcZA/hpmmHvvtsCvAq4C/Bpmp/JF5YMNISfAM6kmQb4soj4Hs0ajVcUTbWwq4GrMvOK+Z9op8702c6IeDHw9vb1M4D/7PmDrGsi4jyaKcXQZL6mXJwlzZ2z/uzBXyIigJHVjzOU70bE79FMkfl2RLwQeCfwSOC/iiZb3o6IeAnw1vb1M2lGZfpsKiKeTvM9BngacEu5OEN5Fc3P4WDN3COBPykXZ3lOnVpCRDwZuDEzPzPv+LHASzLzt8okG05EnE4z7z6ACzPz0sKRFhQRv0TT8ejLC3zu4X3u0APQjhwN5qL+W2b29WnkXvMyb2+fovVa++T0Qe3Lzw6e9ql77fnxUJopVL8KrM/Mnyiban8RcSrwnwt1IYuIe2fmFwvEGkpE3IVmLd1ZNIXFx2lGZ3YBJy304KW09mHEi2nWdUHz0O1lmdm7m7OI+O80U/7+a97xjTSNUfq4SPlEmnVRPwL+kOah0O/SrId5Wmb2tahb6Hy+CPgfmXlT0WBLaM+Fd9NMV0vgKprv8/aSuZYTET9N8zAImu6A15fMsxwLDakDETFKs2gPuG3KXR9FxGnA9Zn5/Yh4Es2CznMXulnri3bq3zcyc09E/DfgAcA7M/PmwtGWFBH/DHyMZp3DFVnBBTciPgo8kGa04CIqyi5JK9Wu26Hvv08G2ofdGzPzkmhaqN+uzw8LLTSGEBF3olnlP7dS/+PM/G7RYAtoFzct1QLtCYt9ro+ip+1tByLiwcA7gHvPPd7nVnMRcQ3NQt97An9P0zb2lJ4v2rsSOJ1mWs9ngUuBIzLzl4sGW0ZEPIymVewjgBNpcl+Umb3spgYQET9Hsxi513uqzFXTNXqgnWby2sHagYj4f4DnZOZLyiZbXFvwv5FmfcmmiNgEnJmZfZ1yWd250eZ9PE33JoAbgQ/3Ne9ARIzTdPa6e2Y+pl1Hd//MfG/haIuK2/Y5uldm/m7Usc/R2TTTp36UmfeIiPsDk5n5S4WjLcpCYwgR8Tc0cw3f3B56FnC3zHxcuVQLi4inLfX57PHeHwuJiBsz827Lv7OMiPgc8AfAXwNntH+fzsxXFQ22hIi4MjNPiaaF5hGZeW70tCPZwJzMvwUck5kv72uHrIW0v4QfTzPt5NjMHCscaUltN6ETmLOOLzN727a5pmv0QCzcSe3K7PcePH9Psz/T87PpUncEzRqZPndRq+bciIgn0jRjuAT4anv4BOBnaIrQDxUJNoSIeB/Nestfycyfjogfo9kTZlPZZIuLOvc5+jzN+rmL8raOotf3uSOZi8GHc+95F6Xfj4h/LZZmCbUVErBse9vetrZtHZmZn4uII9ph11dEs5lfbwsNYDSaVpqP4baFyb0dgWmNttPTfg54deEsQ2sbMZxFM63uYpo5158oGmoZEfGHwEtp9uAZdOdJ5o3a9Uw11+g5FtrHqq+LlAeOyaa9+/MAMvPWiLi1dKhl1HRuvAJ4UGZ+de7BiLgH8A9AbwsNmu/zr7TFEpn5g3bhfZ89uB2ZuwqafY7a9Wl9NptNm+a5x3o7bQosNIb1nxFxdGZ+C/YuSt1ZONOy2vn3m4C9T08zc0uxQIurub3tYOH3dyLiFOBrQF/bPQ78Fc1O4Bdl5pXtcHGvh+WB99L0wv8icFk7R/X7ZSMN5Vk0m5u9hWZh6n4ND3ro92nm//Z2zc4CarxGb4+IF9A8lAjgeUCvd7qnaZu+9w6nnebT95vJms6NdfOLDIDM/I929KjP9rnZbUc0+n5u7DM9tO341veNrG9uHxQOti44i6aBRG/1/cTti+8C10XER9rXvwR8avAkvo837xHxWuAeNIs630uzu/nHioZaXM3tbd8XET8B/CnNk+ojgf9VNtLSMvOt3NaCEOA/aEYKequdKvV6YCozMyJuBs4unWs5mXmXiLgfzRqN10TECcBl2ewW3lc7KysyoMJrNM3c8POBl9PcNHySZnfiPvsAze7J69tr82+z77Wkj2o6Ny6PiLfTTMW9oT12d5rv834tnHvm4oj4n8BYRDyCpmvWYhtq9kWN+xy9kGZ0654RcSnNfd6jykZamms0hhBNr/NF9XHxXkRcB9yfZv7s/dvWc+/s44LfqLy97UA75DpWQ+eKiHgQcC/2nYP/rnKJ1q6I2MBtC8LPAm7q+bzlXwR+Efg75jzxy3aDxz6q8Ro9EM2uz8xvw9pXEfEU4L/TPK2+IDPfUzbR0mo6N9pRgD+i2YtpsDbxBpqGHX+Rmb0dxW1HXJ7PnHMDeGVm9nVzxEG75lfRZIYm85a+/yy2a/5Op/k+X5aZ3yubaGkWGmtURFyemadFxNXAaZn5w4i4rs+L9rQ6IuKNNDsqX82cOfiZ+aRioYYwf6Fs3xfOAkTEdpp59/88+Mh2A9C+ioiX0jQ1+Ar7nh+by6VSH0TEIzPzH5c7JkkDFhpDiIg3Ab+X7UZs7bzUd2XmY8omW1xEfBx4NM3ah6Np5rc/ODMftOQ/LKSdc/9kmg4btwLXA+/JzJmSuYZR2w1wRHwJuG9N7UuhOUcy8+uLve6jiDgpM79UOsdKRMRXgU19f0o2V43XaKjy2rFfvgoyV3NuRMQZS32+56OKX+a2fXc+nm3b5j6LilpMR8QnMvNnIuK77Lt2NWgeBE0UirYs12gM54fA5yLil4E70+wk2ds++K2n0DyNfD7NIsM70tM57e2i9T8HrqEZDrwIOBl4cUT8Ymb2tUPIwPz5kb2eLwl8Heh9ATdfZn49IkYHxWffiwyA2oqM1g01FRmtGq/RUMm1IyLuDfwkMB4Rj53zqXHg9mVSDa2mc2PQrfAImqnPX6G5qbwXzQh0bws6mp2qf46mjfdrIuImmoYjLygba0mPy8y9U+sy89sR8Tigd4UG8Cvtn5tKhjgQjmgMqb1IvQn4L+BJmfmZwpHWjHY9yZntD/k9gVdl5uMj4heAF2TmWYUjLmvuDXBfzblBeBjNZn3vZ985+BeWyDWMdkH1e4A7ZuaGiHgg8OSe/xIDqnxq/Zc088M/SCXnB9R7ja7k2vE04OnAqey7KHkKeHNmfmShf9cXtZ0b7YLw92bmx9rXj6DZn6LvzVGIiPvSFBy/T7NP0/GFIy0qIq7NzPvNO/aFzPzpUpmW0nbF+mhmPqJ0lpWw0BhCNJu4vJNmrvU9gXdn5suLhlpGNLty/gnNbsRzF/zeb7F/U0rM27hq7o1YRPxrZv5UsXDLqOkGOCIuXuLTmZkPX7UwKxQRl9B083pdZj6gbbH5hezxJkUDtU35WuQ86fv5cUfqu0ZXc+0YiIhnZubbSudYiUrPjYVugHu9QWlEnA+cRtOi+SKa9Wj/UjbV0iLiA8Bg36tBi+nTMrOXsz8AIuIymg0Gf1Q6y7AsNIbQzj18S2a+MiLuALwNuFNm9rYlaERcC7wL2Mac/Sky89PFQi0iIv6JpgXvPwDn0KwlObu9mfxiZp5UNOASar4BrklEXJGZp8acHcyj57uZz1XDU+uaVXqNvoQKrx1td717sO/+TNeWS7S0Ss+Nq2i6H13cvv4Z4NV9vt5FxDaa9u7/RNNK/1N9v+ZFxF1pWkw/lDktpjPzpqLBlhARrwFOosl9y+B4n0ecXaMxnGdk5qcAsmld+qSIeE7hTMuZzcy/LB1iSL8LnAe8lmZY/mnt8aOBPyuUaVhHZeal0e5hlZkZEb3epbNSt7Y3OINNio5n4Q0ee2XuU2ugiqfWlarxGl3dtSMiHk2z+eSdaKYh3Ymm/eo9SuZaRo3nxnNo9mj6Ic2T9nU0zVJ6KzM3t6NHD6dZp/GGiLihrwVdOw3p7Mx8eNTVYnow0jV3L6YEelto9H0HxF7IzE9FxBMj4o8BIuI4msq3zy5eroNFX2TmlzPzoZl5h8w8MzNvbI9/s4Jh+lpvgK9c6nUPvZ6mx/nREfFy4FM0DQT67rU0m219q319JT1d8DtXbedHpdfoGq8dLwMeDPxrZv4E8FSatTy9VeO5kZmX0SwA/+/A44ATM/OzRUMto71x/2maReybgAmaYrSXstnf42nt3/+rhiKj/R5/uL1PmvvR22mt4NSpoUTTV/404F6Zee9oWrF+KDNPLxxtURFxOvBR4GaaBZ2DFmj3LBpsAcsVRNnvln7n0HT4uh/NPOBzaBaw/9+iwZZR27oB2HtOP47mXL4wMy8tHGlZtU75qu38qPQaXd21IyI+n5kPjDl7Mg2Olc62mErPjdOA6zPz+9F0ZdwMnJuZ/1k42qIi4jvAdTTtbS8CtmWPN+sDiIg/A67LzP9TOsuwBr9TSudYCadODedxNG3lroC9bTaPKhtpWe8AnkuTudc/7FTc0i8zz4+Ir9CcIyPAOTXcAGedrWIvAy4rnWOFanxqXeP5Ud01utJrxw/bP3dExOOBr9JMn+qz6s4N4K3AKRFxEvAKmlGjd9BstNpXx2ePdy5fxLNpWja/Dfg+FexJAVwUEb9WU3FkoTGcH2Tm7GAubSsWe3NP3JKZby8dYhiZeRrsben3gpzX0q9ktmHUdgNc07qBiPgw+25OtI/MfMIqxjkQ86d8nQP07vs8V03nxxw1XqOru3bQ7I9wJ5pF7O+jOUeeWzTR8mo8N2bbzL8IvDEzz20XiPfZnoh4Pk1rW2gWhb86M28tmGk5m0oHOADVFUcWGsO5ISIeBmT7dPKPaZ6099lHIuIxmfm3pYOswKmZ+RuDF5l5UUS8aql/UErlN8CDdQOva19fSdOhrI83kheUDnAwKn1qXdP5MVDNNbrWa0c7P3xPZn4X+DxN55saVHNuzDEaEXcGHgO8sD22rmCeYZxLMwvhDTTn97No9uP5g5KhlpKZN0TE7bmt4Li6glGZTaUDrJSFxnD+gGYO7X1pFjddDPxa0UTL+32aqvcHNLtA977qBWYj4sx5Lf362iv6gtIBDkI13W4y852lMxysCp9aV3N+zFHTNfqC0gEORPuE/X8CHyqdZYVqOjcG/grYTrOz9pURcS/gu4UzLedngU3Z7u8QER+heUjRW+26vw8Bg3a2d46IJ2aPN3TMzBtKZ1gpC40hZOY3gEe2lW/U0J2ACqteKmrpV/kNcJXrBtpFkZvYt3//lmKBllDrU+tWdedHTdfoyq8dV0bEQysYldurpnNjIDPfSrNOY+A/uG1KUl8FTSfTH8153fcpaufStLj9NOwtPP6KprNaL0XEMcBLaNazzv1d2Nu1rBYaK1DBkNpe7ZDgscDGzLwkIo6g5+2MM/Oy9snNT7aH/i0zf7jUv+mDmm6AWzWuG3gtTa/+B9Js7vjLNJtC9dUFpQMchOrOj4GartFQ5bXjwcDT2+mAczcL6+1NzkCF58aDaKYizb1Pe1ehOMP4R+CfIuK89vVTaTbh7bMfyzmbGLf3IGNL/YMeeBtwKXAWzU7mzwZ6vX7H9rZDiogr515M57/um4g4m6abU2bmCRFxf2AyM3+pcLRFVdrSb8Eb4Mx8ZtFgy6itVWxEXEfzBOeqzLx/RNwFeGdm9rkLS7VqOz+gymt0ddeOdjrrfjLzE6udZSUqPDfeSNNh6mpuG03MzHxSsVDLiIjb0dz0ntUeugh482AqVR9FxKeBF2fmRe3rs4CXZuZ/K5tscRFxdWZuGrSYjogR4BOZ+ZDS2RZjoTGkCvvKfx74eZo5noP+/ddn5sllky0uIq6haUN4T+DvaVr6ndLnm0lvgFdHRFyemadFxNXAaZn5w7m9/PuswqfWVarwGu21Y5VUeG58CbhvZk6XzrKWRcSpNGs0BsXc7YAnZGZv15ZExLZsdmG/HHgkzdqd7ZnZ2+YMTp0aUtbXV342M78zr6Vf3xd01tjSbzozfxQRGRFHZuZNEXHX0qEWUvm6gZvbOdaXAudHxE00rf16raYpX5WfHzVeo6u5dgxExI/RNBrZxL6Fs+dGt75O08Sl9yLiHSx93fiNxT5XWmZeEREnAhvbQ9srmK79xYj4CeB84HPAFE0XuN6y0BhC1NlX/ua2Pd5gQedZwK6ykZZVY0u/mm6ALygd4CA8heap0/Np5qXeETi7ZKAhncltT62fFxF/QdMBp48uKB3gQFV8ja7l2jHwFpobm9NppuY+HfhkyUDLqenciIjHtn/9HPDBiHg/sHdUIzMvLBJsaVe0f94XOIPme5001+xPlQo1jGja578tM79QOsuwMvOc9q+viYgraDbM/MeCkZbl1KkhRMQlNBsUvS4zHxDNMMEX+jgNKSJOzszr2yHBN9NMQ/oCzVPVR2Xm1SXzLSUingX8Jc10r7PbheFvz8wF5wX3QVsYfY9myHVwA/yazPxawVjqiZqnfNWkpmv0QI3Xjjnzwq/NzPtFxB2Aj2TmGaWzLaamcyMiLl7i05mZD1+1MCsUEZ8EHp2ZU+3r9cDf9fzceDHwNODbNDuvvyczd5dNNZz2+7t3sCAze/sg2RGN4dTUV/7dNOscXk3zNPV0mgWdl2Xm98rFWl6NLf2yaZ048PJiQVaotnUDEXEf4E+AE9n34nq/UpmGVONT6+rOD+q6RgPVXjt+0P55a0T8eGbeHBFHF020vGrOjcw8s3SGg3D0oMgAyMypvp8bmfkS4CURcSZNwfGSiPhYZvZ2n5WIeDLNZqp3ohk5ivbPkZK5lmKhMZya+sqPtSfisTTDmINFGmdERF+HXveKylr61XgDXNO6gTneR3MevJ7+/uwtpLopX5WeHzVdo4E6rx3Aroi4E02zjo9GxLeBHYUzLae6c6NS10TT2vZt7etnANeUizO8zLw4Im6hGV18Mv3e0HES+KXMvGLZd/aEU6eGEBHn0Nww3I9mfvU5wAsy8/8WDbaAdo7nbwMP47a5kwN9H3qtsaXftTQ3wNuY88sr5/Tm7psau91ExFWD7mk6tCo9P6q5Rg9Ueu1Y1zbsCOBXaZ6qvmvuk+y+qfHcgCpb8h4FvBgY3GNcBLwsM29Z/F+VFc3md79OUxQFzfSp8zPzpiX/YUERcVlmnl46x0pYaAwpKusrHxGvycznls6xElFhS78ab4BrXDcQEecCF2RmrxeezlfjU+sazw+o8hpd3bUDICJ+HHgAzQjB1VnBTtu1nRtQX0veGrUjch8EzsvMz5bOM4yIeDpwV5rccxsF3Fgq03KcOjWkzLwMuKx0jmHVVmS0qmnpN8fFEXFGZTfANa4b+CDNVI2baS6uQTPadc+ysZZV45SvGs+P6q7RVHjtiKZ74XuAwSaqx0bEUzJzqUXMxVV4blTTkrf9///eiPiDhT6fma9d7UwrcHxm/mD5t/XKKE1zgz9izswP4JhiiZbhiMYSovK+8rWY09LvYTRdsmpo6QfsfVL2UaCaG+BKu91sB15JMx1w7jST64uFGkKNT61rOj9qvkZXeu24DnhWZn6ufb2Zpj1o70a7Kj839rbkzcy+t+R9SWa+OJr9NObL7OE+GjUXRxFxA/DwzPz30lmG5YjG0i4oHeAw8YfzXv/OnL8n0NtCg2ZO53OZdwPcZ5V2u7klM99eOsQBqO6pdWXnxwWlAxyE6q4dwI8GRQZAZm6LiL5mv6B0gIPwWpq1lq9rX19JMzLau0IjM1/c/vmM0llW4CfbPxd6CNT3p+87aioywBEN6aBExOcz84Glc6xEpesGXgpcnpl/WzrLSlT61Lq686NGlV47/gK4nts2nfx14Kf7+KS9ZhFxRWaeOndEtO+joxHxzzTd6f4ZuCK9uTwkIuJlwO3Zf+bHtcVCLcNCY0gV9pXXKqjxBrjSbjffBcZp+vjPcNsN+0TRYMuoccpXjecH1HeNrvTaMfg5/GF76EhgsMFZb38eKzw3PkszlfhzmXlK25L3w5l5auFoi4qIh9Hse/UImocUl9JsvvuGosEWEBFPXerzmdnnlvr/scDhfj+8stBY3mJ95TPzmUWDrUEVtvSr7ga470/GFhIRd1/oeGbesNpZVqLSp9Y1nh/VXaMrvXYs+HM40Mefx0rPjSpb8gJExDjweJpWt8dm5tgy/2TVRcQH2r+uB36GpihK4KHAJzLzkaWyLSUi1gGPz8wPls6yEhYaQ6ixr3ytamvpV+MNcMWtYo8FNmbmJRFxBHC7zOzlDr8DlT61ru78qPEaXeO1A+r7Oazx3ID6WvJGxMuBs2i6Il1MM4XqE9nj9sdtw4D/nZlfaF+fDLw0M59YNtni+v7wdSEuBh/OdGb+KCIyIo7MzJsi4q6lQ61FtbT0G8jMGxb6xVs61zKqaxUbEWcDr6J56nQCcDLtDqkFYw3j94HxiKjmqTUVnh9UeI2u8dpR6c9hdecGVNmS91nAV4C30IwYfblwnmGcOCgyoJnSGhEnlQw0hCsj4qF9LzznstAYTpV95Wszt6Uf0OuWfgOV/uKtsdvNVuAUmt1mycxrlpvG0RObSgc4ADWeH9Vdoyu9dtT4c1jNuVFzS97MvEv7O/wRwGsi4gTgssz8zbLJljQVzQZ4g+YGTwN6u5N568HA0yPiK8zJ2udRDguN4TyF5hfu87mtr/zZJQOtUdW09Jujxl+8NbaKnc3M70TE3GO9na4xUONTa+o8P2q8Rtd47ajx57Cmc+OC0gEO0i7guzT78PwEcFrRNMv7DeDdwJva11fRFBt99pzSAVbKQmMIWVdf+ZodlZmXDn6JZWZGRN9/idX4i/cjEfGYmtYN0DyVvDPt075odijeVTbS8ip9al3d+VHpNbrGa0d1P4c1nRuZ+c7l39VPbYe9EZq1GX8H/GFmfrNsqqVl5nZgc0TcoX19c+FIy8rMT0TEkcDdspL9NCw0hmBf+VVza/sDNPgldjz9n7pR3S9eKlo3EBEnt61gXwj8A3DPiLiUpovMo4qGG06NT62rOT8GKr1G13jtqO7nsNJzo7qWvMCjM/NLpUMciBoKjIGI+FmaKea3AneLiNOA52bmOSVzLcVCYzjvo5nC83r6f+Nbs9fTDB0f3XawOIeeTpuq/AZ4U+kAK/Bumhv1VwNnAqfT3PhelpnfKxdraDU+td5UOsABqOYaXfO1IzOviIjafg6rOTcGFmvJWzTUMmotMuZ3caqgq9Of0eyx8kGAzLw8InrdjtxCYzizmfmXpUOsdZl5frvA6XE0Q7Dn9LizQrU3wJWtGxiLiCcDxwJn0HyPAc6ICDLzwnLRhlLdU+vKzo+Bmq7R1V472lHmb2TmP0TEfwPOiYh39vyJcE3nxsCZ3NaS93nR7Mje+2lVFd60w/7Ffa+LfWBdZv57TQ+vLDSGc3FEnFFTX/laVdTSr9ob4MrWDbyIpkHAMcD8aQMJ9PL7XPNT68rOj4GartHVXjuAvwFOj4jjaEYKLqXZ8OyXi6ZaWk3nxkCVLXmp76a9upb6wHREHMVtD6/uS7PpZ2+5Yd8Qotk456NATX3lq1FjS7+IeCzNDfDDaNqAzpWZ+fDVTzWciPg88PPARdnuAB0R12fmyWWTLS4iXpOZzy2dY1iDJ3ltcfEo6npqXeP5Uc01uvJrx+C8/i3gmMx8eURck5n3L51tMTWdGwMR8XHg0cCfA0cDNwEPzswHFQ02hLk37X03t6V+ZtbSUv/nadYc3YtmOt0jgF/NzI+XzLUURzSGU2Nf+ZpcUDrASrVPHS+s7Qa4Vd26gQq/xzU/ta7u/KCia3Tl147RiBgFfo5m6lcNqjk35qipJS9Q5z5YVNhSPzP/KSK+BDyS5vfKi/vefcpCYzg19pWvRs0t/Sq8UYAK1w1UqMopX60az4/qrtGVXjveS/N0/YvAZe1anl5ufjdHjedGNS1556jupp0KW+pHxGOAj2TmG0tnGZaFxnCq6ytfqwpb+lWj5nUDtanxqXXl54fX6FXQTpV6PTDV3pTdTM+ftFPhuVFpS97qbtqps6X+FuCvI+L/AO/IzH8tHWg5rtEYQkR8FxinWXBTRV/5Gi3W0i8zn1k02BpR87oBHXo1nx9eo7WYGs+NiLiWZjRgG3NufDPz08VCLSMiPkuz7uhz7XXkeODDmXlq4WiLiohzaKap3Y+mq9c5wAsy8/8WDbaMiLgn8NT24xvA2zPzLWVTLc5CYwiLbbCVmTesdpa1LCKu47aWfvePiLsA78zMXygcbU2IiH8BXgL8KfA/uG3dAEDf1w3oEKv5/PAavXpqa2Fa47kREVcNGjHUouKb9tNpWuoHcGGPW+rvJyLGgL8Cfisz15XOsxgLjSEt1Fc+M/s+LFiViLg8M0+LiKuB0zLzhxFxXWbet3S2taDmbjc69Go/P7xGr46IOHZuC9D5r/uotnMjIs4FLqisJW/VN+01iYhTgGcATwIuB87LzA+WTbU4C40hzO0rn5knRMT9gcnM7HNf+erU3NKvJjWtG9Dqq/H88Bq9uiprYVrduVFjS96a1NhSf6CdVjcCnEcz46PXRT5YaAylxr7yNWo73XyPZhfiQUu/12Tm1wrGktRzXqNXR6X7DlR3bkTEduCVzGvJ2zZr6JUab9oj4mlLfb7PnTAj4vRsNjauhl2nhlNjX/nqVNrST1J5XqNXR40tTGs8N2pqyXtB6QAr1edCYgjbIuJ5wImZ+TsRcS/g7umGfdWrsa98dSpt6SepPK/Rq6PGFqY1nhvVtOSt/Ka9xpb6rwfWAQ9tX38HeD/Q2+5eFhpLqLyvfI3eR/N07PX0v5e1pMK8Rq+6avYdqPzc+H1gPCKqackL9d20L9ZSv2io5T04MzdFxFUAmfm99meytyw0lvZu4BTg1cCZVNRXvlKzmfmXpUNIqobX6NX1epqpMkdHxMtpW5gWTbS4ms+NTaUDrFSlN+1ncltL/edFxF/QtObts+m5LyJiHc261t7qdbgeGIuIJwPHAmcAR9IUZ2e0rSDVrYsj4ozSISRVw2v0KsrM84FX0CwIHwHO6fE+CdWeG+0eH3uAe7R/3wn0vbvQmTStbb+Vmc8DNgMbykZa1nRm/gjIiDgyM28C7lo61DKubfcsuV1EnAj8NXBJ2UhLc0RjaS+iWfh2DM2273Ml0NsNrCr1QeCjEWFLP0nD8Bq9ytqONzV0van23Jjbkhc4ATgZmAR625KX9qY9IvbetEdE32/ab46I2wOXAudHxE3A9wtnWs4WmnPjLsCnaUYYX1Qy0HJsbzuEGvvK16imln6S+sNr9KFVYwvTgRrPjUpb8la3D5Yt9VeHhYZ6IyI+n5kPLJ1DknSbmvcdqFFEbMvMzRFx1ZxCY+/f+8ib9tUTEQ8C7sW+3TnfVS7R0pw6pT6ppqWfJB0uLCRWXXUteWvcB6vGlvoR8UbgF4CruW3mR9J07OwlRzTUGxHxXWAcqKqlnyQdLmprYVqTQUveiDgVeDNwT+ALtC15M/PqkvmWUulN+7U0N+jb2He69qeLhVpGRHwJuG9mTi/75p5wREN9sql0AEnSwiptYVqTmlvy1rgPVo0t9b9O8yC2GhYa6o3MvCEijgU2ZuYlEXEEtmCWpL6ocd+BmsxvyRvt8TMigszsbacs6rxpvzgizsjMT5YOsgKfAz4YEe9nzp4afT43LDTUG5W29JOkw0WNLUxrUm1LXuq8aa+xpf6p7Z+/M+dYr88N12ioN2ps6SdJh4saW5jWqNKWvKcDHwWquWm3pf7qcERDfTKbmd+JiLnH9pQKI0nax1Nobsiez20tTM8uGWgtqq3IaL0DeC7zbtp77pbMfHvpEGudhYb6pLqWfpJ0uKixhalWTY037VW21I+IKzPzlMVe941Tp1RczS39JOlwUWMLU62OiHgpcHlNN+21ttSPiGMz8+uLve4bCw0VN6jGI+JS4FHU1dJPkg4LNe47oNVR4017RNx9oeOZecNqZ1mpiBjNzCra3Dp1Sn1Qc0s/STpc1NjCVKtjU+kAK1VjS/2IuB/wHpr1URsi4oHAkzPzBUWDLcERDRUXEY+laen3MJqFZHNlZj589VNJkuaKiHOBCyprYapVstBNe2b2tqHL3Jb6mXlCRNwfmMzM3rbUj4hLgP8FvC4zHxBN95wv9Lk7pyMaKq4dsbiwxpZ+knQYqXHfAa2CSvfB2kqzE/tFAJl5zWLTqXrkqMy8dNCdMzMzInpbzIGFhnrEIkOSeq3GFqZaHTXetNfYUv/WiDiS27pzHk/PfxYtNCRJ0jBqbGGq1VHjTXuNLfVfD1wAHB0RLwfOAXq7PgMsNCRJ0nCq3HdAq6Kam/ZBS33ghcA/APdsu17eg6bzZW9l5vkR8RXgccAIcE5mXlo41pJcDC5JkpZVYwtTHVo17oNlS/3V5YiGJEkaxqbSAdQ776ZZm/Fq4EzquGmvrqV+RHyYdrRoIZn5hFWMsyIWGpIkaVk17jugQ666m3bgRTQt9Y8Btsz7XAJ9zHxB6QAHyqlTkiRpWTXuO6BDq+Z9sGypvzosNCRJ0rIi4vPAzwMXZeYD2mPX93mzMK0Ob9pXT0Q8iWYa49jgWGbOH5npDYc8JUnSMGYz8zvzjvW9halWgUXG6oiI1wK/DjydZprX2TQNGnrLQkOSJA2jmham0hp1Jk1r229l5vOAzcCGspGW5mJwSZK0qJr3HZDWmOnM/FFEZEQcmZk3RcRdS4daioWGJElaSo0tTKW16OaIuD1wKXB+RNwEfL9wpiU5dUqSJC1lfgvTI2keVJ7Rdh2StDqeAswCzweuA35Is06jt+w6JUmSFlVzC1NJZVloSJKkZdnCVCorIu4D/AlwInOWP2Tm/UplWo6FhiRJktRzEXEt8C5gG80UKgAy89PFQi3DQkOSJEnquYi4arBZZi1cDC5JkiT138URcUbpECvhiIYkSZLUcxFxOvBR4GZgmqbNdGbmPYsGW4L7aEiSJEn99w7guTTd32aXeW8vOKIhSZIk9VxEfD4zH1g6x0q4RkOSJEnqv49ExGNKh1gJRzQkSZKknouI7wLjwA+AGW5bozFRNNgSXKMhSZIk9d+m0gFWyhENSZIkqQIRcSywMTMviYgjgNtl5p7SuRbjGg1JkiSp5yLibOCzwHntoZOBC0rlGYaFhiRJktR/W4FTgO8CZOY1wN2LJlqGhYYkSZLUf7OZ+Z15x3o7bQosNCRJkqQa3BwRdwYSICLOAnaVjbQ0F4NLkiRJPRURJ2fm9RFxKvBm4J7AF4B7AI/KzKtL5luKhYYkSZLUUxFxZWaeEhGXAo8CTqfZQ+OyzPxe0XDLcB8NSZIkqb/GIuLJwLHAGTRFBsAZEUFmXlgu2tIc0ZAkSZJ6KiIeC/w28DDginmfzsx8+OqnGo6FhiRJktRzEfGazHxu6RwrYaEhSZIkqXO2t5UkSZLUOQsNSZIkSZ2z0JAkSZLUOQsNSZIkSZ2z0JAkSZLUOQsNSZIkSZ37/wG0FFeD3SaIBQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 960x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "result_coeff = result_coeff.sort_values('Estimate.')\n",
    "\n",
    "x = result_coeff.index\n",
    "\n",
    "coef = result_coeff.iloc[ : , 0 ].to_numpy()\n",
    "\n",
    "sd_error = ( result_coeff.iloc[ : , 1 ] * t_value ).to_numpy()\n",
    "\n",
    "figure(figsize=(12, 6), dpi=80)\n",
    "\n",
    "plt.errorbar( x = x , y = coef , yerr = sd_error , linestyle=\"None\" , color = \"black\", \\\n",
    "              capsize = 3 , marker = \"s\" , markersize = 3 , mfc = \"black\" , mec = \"black\" )\n",
    "plt.xticks(x, x, rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we estimate and plot confident intervals, first \"pointwise\" and then the joint confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-0.1545932731827906,\n",
       " 0.13609550288298192,\n",
       " 0.13693939221511991,\n",
       " 0.023302771994035233,\n",
       " 0.18685348640356858,\n",
       " 0.027810312784543743,\n",
       " -0.11933503921829447,\n",
       " -0.01288483245739874,\n",
       " 0.010138550922289817,\n",
       " -0.0304637472338401,\n",
       " -0.0010634393390639041,\n",
       " -0.008183337595043904,\n",
       " -0.004226127118997856,\n",
       " 0.004732893281720598,\n",
       " -0.15951933512468575,\n",
       " 0.026512194276155938]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "effect_female.res.get( 'coefficients' ).iloc[ :, 0 ].to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def confint_joint_python( rlassomodel , level = 0.95 ):\n",
    "    \n",
    "    coef = rlassomodel.res['coefficients'].to_numpy().reshape( 1 , rlassomodel.res['coefficients'].to_numpy().size )\n",
    "    se = rlassomodel.res['se'].to_numpy().reshape( 1 , rlassomodel.res['se'].to_numpy().size )\n",
    "    \n",
    "    e = rlassomodel.res['residuals']['e']\n",
    "    v = rlassomodel.res['residuals']['v']\n",
    "    \n",
    "    n = e.shape[0]\n",
    "    k = e.shape[1]\n",
    "    \n",
    "    ev = e.to_numpy() * v.to_numpy()\n",
    "    Ev2 = (v ** 2).mean( axis = 0 ).to_numpy()\n",
    "    \n",
    "    \n",
    "    Omegahat = np.zeros( ( k , k ) )\n",
    "    \n",
    "    \n",
    "    for j in range( 0 , k ):\n",
    "        for l in range( 0 , k ):\n",
    "            Omegahat[ j , l ] = ( 1 / ( Ev2[ j ] * Ev2[ l ] ) ) * np.mean(ev[ : , j ] * ev[ : , l ] )\n",
    "            Omegahat[ l , j ] = ( 1 / ( Ev2[ j ] * Ev2[ l ] ) ) * np.mean(ev[ : , j ] * ev[ : , l ] )\n",
    "    \n",
    "    var = np.diagonal( Omegahat )\n",
    "    \n",
    "    B = 500\n",
    "    sim = np.zeros( B )\n",
    "    \n",
    "    for i in range( 0 , B ):\n",
    "        beta_i = mvrnorm( mu =  np.repeat( 0, k ) , Sigma =  Omegahat / n )\n",
    "        sim[ i ] = ( np.abs( beta_i/ ( var ** 0.5 ) ) ).max()\n",
    "    \n",
    "    t_stat = np.quantile( sim , level )\n",
    "    \n",
    "    low_num = int( np.round( ( 1 - level ) / 2, 2) * 100 )\n",
    "    upper_num = int( np.round( ( 1 + level ) / 2, 2) * 100 )\n",
    "    \n",
    "    sd_tstat = t_stat * ( var ** 0.5 )\n",
    "    \n",
    "    low = coef - sd_tstat\n",
    "    upper = coef + sd_tstat\n",
    "    \n",
    "    table = pd.DataFrame( { f'{low_num}%' : low.tolist()[0] , f'{upper_num}%' : upper.tolist()[0] , \\\n",
    "                           'Coef.' : rlassomodel.res.get( 'coefficients' ).iloc[ :, 0 ].to_list() }, \\\n",
    "                         index = rlassomodel.res.get( 'coefficients' ).index )\n",
    "    \n",
    "    return ( table, sd_tstat )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mvrnorm( mu , Sigma , n = 1 ):\n",
    "    p = mu.size\n",
    "    \n",
    "    ev = np.linalg.eig( Sigma )[ 0 ]\n",
    "    vectors = np.linalg.eig( Sigma )[ 1 ]\n",
    "    \n",
    "    X = np.random.normal(0, 1 ,  p * n ).reshape( n , p )\n",
    "    \n",
    "    diagonal = np.diag( np.maximum( ev, 0 ) ** 0.5 )\n",
    "    \n",
    "    Z = mu.reshape( p , 1 ) + vectors @ ( diagonal @ X.transpose())\n",
    "    \n",
    "    resultado = Z.transpose()\n",
    "    \n",
    "    return( resultado )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we compare the pointwise confidence intervals to joint confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "joint_CI =  confint_joint_python( effect_female, level = 0.9 )[0]\n",
    "size_err =  confint_joint_python( effect_female, level = 0.9 )[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>5%</th>\n",
       "      <th>95%</th>\n",
       "      <th>Coef.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>female</th>\n",
       "      <td>-0.287258</td>\n",
       "      <td>-0.021928</td>\n",
       "      <td>-0.154593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:widowed</th>\n",
       "      <td>-0.120818</td>\n",
       "      <td>0.393009</td>\n",
       "      <td>0.136096</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:divorced</th>\n",
       "      <td>0.077758</td>\n",
       "      <td>0.196121</td>\n",
       "      <td>0.136939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:separated</th>\n",
       "      <td>-0.110045</td>\n",
       "      <td>0.156651</td>\n",
       "      <td>0.023303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:nevermarried</th>\n",
       "      <td>0.132002</td>\n",
       "      <td>0.241705</td>\n",
       "      <td>0.186853</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:hsd08</th>\n",
       "      <td>-0.355334</td>\n",
       "      <td>0.410954</td>\n",
       "      <td>0.027810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:hsd911</th>\n",
       "      <td>-0.261893</td>\n",
       "      <td>0.023223</td>\n",
       "      <td>-0.119335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:hsg</th>\n",
       "      <td>-0.062645</td>\n",
       "      <td>0.036875</td>\n",
       "      <td>-0.012885</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:cg</th>\n",
       "      <td>-0.039213</td>\n",
       "      <td>0.059490</td>\n",
       "      <td>0.010139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:ad</th>\n",
       "      <td>-0.092722</td>\n",
       "      <td>0.031795</td>\n",
       "      <td>-0.030464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:mw</th>\n",
       "      <td>-0.052152</td>\n",
       "      <td>0.050025</td>\n",
       "      <td>-0.001063</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:so</th>\n",
       "      <td>-0.060122</td>\n",
       "      <td>0.043756</td>\n",
       "      <td>-0.008183</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:we</th>\n",
       "      <td>-0.062703</td>\n",
       "      <td>0.054251</td>\n",
       "      <td>-0.004226</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:exp1</th>\n",
       "      <td>-0.015694</td>\n",
       "      <td>0.025160</td>\n",
       "      <td>0.004733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:exp2</th>\n",
       "      <td>-0.278207</td>\n",
       "      <td>-0.040832</td>\n",
       "      <td>-0.159519</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female:exp3</th>\n",
       "      <td>0.004007</td>\n",
       "      <td>0.049017</td>\n",
       "      <td>0.026512</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           5%       95%     Coef.\n",
       "female              -0.287258 -0.021928 -0.154593\n",
       "female:widowed      -0.120818  0.393009  0.136096\n",
       "female:divorced      0.077758  0.196121  0.136939\n",
       "female:separated    -0.110045  0.156651  0.023303\n",
       "female:nevermarried  0.132002  0.241705  0.186853\n",
       "female:hsd08        -0.355334  0.410954  0.027810\n",
       "female:hsd911       -0.261893  0.023223 -0.119335\n",
       "female:hsg          -0.062645  0.036875 -0.012885\n",
       "female:cg           -0.039213  0.059490  0.010139\n",
       "female:ad           -0.092722  0.031795 -0.030464\n",
       "female:mw           -0.052152  0.050025 -0.001063\n",
       "female:so           -0.060122  0.043756 -0.008183\n",
       "female:we           -0.062703  0.054251 -0.004226\n",
       "female:exp1         -0.015694  0.025160  0.004733\n",
       "female:exp2         -0.278207 -0.040832 -0.159519\n",
       "female:exp3          0.004007  0.049017  0.026512"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "joint_CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAH4CAYAAADNU5vyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAABFxklEQVR4nO3de5xlVXnn/89jQVWNwSqt/CAiDaKCTESlRWgNUQxgEpN4SZRozBAviRPNxTBp46WdmR/iJRVjJGqceFdUoibRkZCYxIiiEYk2yE0waSRGsDugxgYKolUtxTN/7H26q6vr2r3r7L1Ofd6vV73os8/B/ubkUGc9e631rMhMJEmSJKlJ92o7gCRJkqTBY6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIad1DbAVZjZGQkDz300LZjSJIkSQJ27NixKzNHFnquqELj0EMPZfv27W3HkCRJkgRExHcWe86lU5IkSZIaZ6EhSZIkqXEWGpIkSZIa12ihERHHRsRlEXFDRFweEccv8Jp7RcR5EfHViLg2Ii6JiGOazCFJkiSpXU3PaLwDeGdmPhR4PXD+Aq95KvDjwAmZ+Ujg08DvN5xDkiRJUosaKzQi4jDgJOCC+tLHgCMXmK1IYAQYjYgAxgBbSUmSJEkDpMn2tkcCt2Tm3QCZmRFxM3AUcOOc1/01cBpwK3AnsAN4QoM5JEmSJLWsjc3gJwEPB44AHkC1dOrtC70wIjZHxPbez1133dXHmJIkSZL2V5OFxjeBwyPiIIB6WdRRwM3zXvcc4DOZeXtm3gO8n2qGYx+ZeV5mbuj9HHLIIQ3GlSRJkrRWGis0MvPbwJXAWfWlZwDbM/PGeS/9OnB6RAzXj58MXNdUDkmSJEnta3KPBsALgfMj4pXAFPB8gIh4N3BRZl4E/B/gR4FrIuIHVHs1XtRwDkmSJEktisxsO8OKbdiwIbdvt0GVJEmS1AURsSMzNyz0nCeDS5IkSWpc00unJGlNzMzMMDMzs8/1kZERRkZGWkgkSZKW4oyGpCJMTk4yPj6+z8/k5GTb0SRJ0gLcoyGpCL0ZjU2bNgGwdetWwBkNSZLatNQeDZdOSSpCr6AYGhoCYGxsrOVEkiRpKS6dkiRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktS4g9oOoLUzMzPDzMzMPtdHRkYYGRlpIZEkSZLWC2c0Btjk5CTj4+P7/ExOTrYdTZIkSQMuMrPtDCu2YcOG3L59e9sxitGb0di0aRMAW7duBZzRUNmOP/54AK6//vqWk0iSpIjYkZkbFnrOpVMDrFdQDA0NATA2NtZyIkmSJK0XLp2SJEmS1DgLDUmSJEmNs9CQJEmS1DgLDUmSJEmNs9CQJEmS1DgLDUmSJEmNs9CQJEmS1DgLDUmSJEmNa7TQiIhjI+KyiLghIi6PiOMXed0jIuKzEfHP9c/Tm8whSZIkqV1Nnwz+DuCdmXl+RJwJnA+cPPcFEXFv4K+A52TmpRExBEw0nEOSJElSixqb0YiIw4CTgAvqSx8DjoyIY+a99JeBL2bmpQCZOZuZ32kqhyRJkqT2Nbl06kjglsy8GyAzE7gZOGre6x4GzETE30TE1RHxgYg4dKH/wYjYHBHbez933XVXg3ElSZIkrZU2NoMfBDwReCHwKGAH8LaFXpiZ52Xmht7PIYcc0seYkiRJkvZXk4XGN4HDI+IggIgIqtmMm+e97mbgkszcUc96XAA8tsEckiRJklrWWKGRmd8GrgTOqi89A9iemTfOe+lfACdHxFj9+GeBa5rKIUmSJKl9TXedeiFwfkS8EpgCng8QEe8GLsrMizLz5oj4feCyiLiHaunUrzecQ5IkSVKLGi00MnMb8GMLXH/BvMcfBD7Y5N8tSZIkqTs8GVySJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkOSJElS4yw0JEmSJDXOQkNSMaanp5mdnWV2dpapqSmmp6fbjiRJkhZxUNsBJGklpqenOeKII9i5cycA4+PjTExMsGPHDkZHR1tOJ0mS5nNGQ1IRdu3atbvI6Nm5cye7du1qKZEkSVpKo4VGRBwbEZdFxA0RcXlEHL/EayMiPhMRtzeZQZIkSVL7mp7ReAfwzsx8KPB64PwlXvu7wL82/PdLkiRJ6oDGCo2IOAw4CbigvvQx4MiIOGaB1x4P/DzwB039/ZIG2/DwMBMTE3tdm5iYYHh4uKVEkiRpKU3OaBwJ3JKZdwNkZgI3A0fNfVFEHAy8C3ghMNvg368F2KVHg2J0dJQdO3Zw3HHHcdxxx3HHHXe4EVySpA5ro+vUOcD/zcx/joijl3phRGwGNvcej4+Pr3G0wWKXHg2a0dFRhoaGABgbG2s5jSRJWkqTMxrfBA6PiIOg2uxNNZtx87zXPQF4cUR8A7gUGIuIb0TEofP/BzPzvMzc0Ps55JBDGow7+OzSI0mSpLY0Vmhk5reBK4Gz6kvPALZn5o3zXvf4zHxgZh4NPA6YysyjM/M7TWWRJEmS1K6mu069EHhhRNwAvAJ4PkBEvDsintrw3yVJkiSpoxrdo5GZ24AfW+D6CxZ5/TeA+zaZQXv0uvTMXT5llx5JkiT1QxubwdUnvS49GzduBGDr1q0MDw+7EVySJElrzkJjwNmlR5IkSW1oeo+GJEmSJFloSJIkSWqehYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxh3UdgBJ/TczM8PMzMw+10dGRhgZGWkhkSRJGjTOaEjr0OTkJOPj4/v8TE5Oth1NkiQNCGc0pHVoy5YtbN68mU2bNgGwdetWAGczJElSYyw0pHWot0RqaGgIgLGxsZYTSZKkQWOhIR0g9ztIkiTtyz0a0gFyv4MkSdK+nNGQDpD7HSRJkvZloSEdIPc7SJIk7culU5IkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIa5zka6pSZmRlmZmb2ud47q0KSJEllcEZDnTI5Ocn4+Pg+P5OTk21HkyRJ0io4o6FO2bJlC5s3b2bTpk0AbN26FcDZjDUwPT3N7OwsAFNTUwwPDzM6OtpyKkmSNCgsNNQpvSVSQ0NDAIyNjbWcaDBNT09zxBFHsHPnTgDGx8eZmJhgx44dFhuSJKkRLp2S1qFdu3btLjJ6du7cya5du1pKJEmSBo2FhiRJkqTGWWhIkiRJapyFhrQODQ8PMzExsde1iYkJhoeHW0okSZIGjZvBpXVodHSUHTt2sHHjRqDq7mXXKUmS1CQLDWmdGh0dtbuXJElaMy6dkiRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjbPQkCRJktQ4Cw1JkiRJjWu00IiIYyPisoi4ISIuj4jjF3jN6RGxNSK+GhHXR8QfRoQFjyRJkjRAmh7gvwN4Z2Y+FHg9cP4Cr7kN+KXMfBjwaOAU4DkN55AkSZLUosYKjYg4DDgJuKC+9DHgyIg4Zu7rMvOqzPx6/edp4Grg6KZyqHzT09PMzs4yOzvL1NQU09PTbUeSJEnSKjV5MviRwC2ZeTdAZmZE3AwcBdy40L8QEfcHzgSe3GAOFWx6epojjjiCnTt3AjA+Ps7ExAQ7duxgdHS05XSSJElaqdb2RkTEGPDXwB9m5hWLvGZzRGzv/dx11139Dam+27Vr1+4io2fnzp3s2rWrpUSSJEnaH00WGt8EDo+IgwAiIqhmM26e/8KIuA/w98BfZeZ5i/0PZuZ5mbmh93PIIYc0GFeSJEnSWmms0MjMbwNXAmfVl54BbM/MvZZNRcQhVEXG32fma5v6+yWt3MzMDFNTU3vthZmammJmZqbtaJIkaUA0vXTqhcALI+IG4BXA8wEi4t0R8dT6NWcDm4CnR8TV9c//bDiHCjU8PMzExMRe1yYmJhgeHm4p0WCanJxkfHycbdu2sW3bNsbHxxkfH2dycrLtaJIkaUBEZradYcU2bNiQ27dvbztGcY4/vjrO5Prrr285ycpMT0+zceNGALZu3crw8HARG8FLep9nZmYWnL0YGRlhZGSkhUQrV9L7LEnSoIuIHZm5YaHnmuw6JTVidHSUoaEhAMbGxlpOM5hKKCgkSVLZPJFbkiRJUuMsNCRJkiQ1zkJDkiRJUuMsNCRJkiQ1zkJDkiRJUuMsNCRJkiQ1zkJDasD09PRep2xPT0+3HUmSJKlVnqMhHaDp6WmOOOIIdu7cCcD4+DgTExPs2LGjiIMGtXZKPhhRkqQD5YyGdIB27dq1u8jo2blzJ7t27WopkbpicnKS8fHxfX4mJyfbjiZJ0ppzRkOS1siWLVvYvHkzmzZtAmDr1q0AzmZIktYFCw1JWiO9JVJDQ0MAjI2NtZxIkqT+cemUdICGh4eZmJjY69rExATDw8MtJZIkSWqfMxoDrLcRdXZ2FoCpqSnAjahNGx0dZceOHWzcuBGolscMDw+7EVySJK1rzmgMsN5G1G3btrFt2zY3oq6h0dFRhoaGGBoaYmxszCJDkiSte85oDLDeRtT5ujyb4SyMJEnSYHBGY4CNjIwwNja2z0+XB+zOwkiSJA0GZzTUKSXOwkiSJGlfFhrqFJdISZIkDQaXTkmSJElqnIWGJEmSpMZZaEiSJElqnIWGJEmSpMZZaEiSJElqnIWGJEmSpMbZ3laSJKlDZmZmmJmZ2ee6LeBVGmc0JElFm5mZYWpqap+fhQZqUgkmJycZHx/f52dycrLtaNKqWGiskF9kktRNDso0aLZs2cIdd9zBcccdx3HHHccdd9zBHXfcwZYtW9qOJq2KhcYK+UUmSd3koEyDZmRkhLGxMYaGhhgaGmJsbIyxsTGXTak47tFYoS1btrB582Y2bdoEwNatWwH8j16SWtZbtz40NATA2NhYy4kGl3sHJK2GMxor5N0FSdJ65+y+pNVwRkOSJK1IibP7zsJI7XFGQ5KkFpTYZKTE2X1nYaT2WGhIktQCB8D9YbMAqT0unZIk7eYyk/4pcRlSiWwWILXHGQ1J0m7eZe+fEpchSdJqOKMhSdrNu+ySpKZYaEiSdnOZiSSpKRYakiRJWnfck7b23KMhSZKkdafEPWmltcW20JAkSdK6U2Lr49KKI5dOSZKKNz09zezsLABTU1MMDw8zOjracipJXVbinrTSGnZYaEiSijY9Pc0RRxzBzp07ARgfH2diYoIdO3ZYbEgaKKUVRy6dkiQVbdeuXbuLjJ6dO3eya9eulhJJB643Szc7O8vU1BTT09NtR5JWzRkNSZK0Yi5TW3vO0mlQWGhIkqQVcQDcH0vN0vk+q6Ri30JDklS04eFhJiYm9hqYTUxMMDw83GKqweQAWGpXacW+hYYkqWijo6Ps2LGDjRs3AlUXli7f4VP/lXQHWFpKacW+hYYkqXijo6PFdGFRf5V2BxjKnKXzlG0txK5TkorQOw11bheWLp+GWjK73WgxvQHwXF0fAJfYlaw3Szf3ILkuF0ZQ3kFy6g8LDUlF6H2Jbdu2jW3btvkltkZ6d3/nvs9HHHGExYaAMgfAperN0g0NDTE2Ntb597jEU7ZLVFqx79IpSUXonYY6n1PyzSpt/W/pStw74DI1LaS0g+RKVdqeNAsNSUVwna8GTYl7B0pU4n4HaSklFfsWGpIktcDZo/4o7Q5wyUqcodPastCQDlCv08bcX67gHXhVSvvi9e6vBlFJd4BL5Qxdf5Q25nAzuHSA3KSsxZS4sdrNvpL2R4ndvUpU2pjDGQ3pALlJWYspdWmMd3/7w9kjSatV2pjDQkM6QF2drpTUbe4dkLRapY05XDolSVJLSjsrQVpMaec7qD8sNCRpjfjFK2m9cH+XFuLSKUlaIy6N6Y/SurBIg8r9XZqv0RmNiDg2Ii6LiBsi4vKIOH6R1/1aRHwtIv41It4VEQc3mUOSusKlMWuvtC4skrReND2j8Q7gnZl5fkScCZwPnDz3BRHxIOA1wInAt4C/An4d+D8NZ5EkrQOldWGRBpEzi1pIYzMaEXEYcBJwQX3pY8CREXHMvJeeCVyUmbdmZgJvB57dVA5J0voyMjLC2NjYPj8ObqT+cWZRC2lyRuNI4JbMvBsgMzMibgaOAm6c87qjgJvmPP5GfU2SJHWYd621GGcWtZBOd52KiM0Rsb33c9ddd7UdSZKkdcu71v0xMzPD1NQUs7OzzM7OMjU1xdTUFDMzM21HW5Qzi1pIk4XGN4HDI+IggIgIqpmKm+e97mbggXMeH73AawDIzPMyc0Pv55BDDmkwriRJWo0tW7Zwxx137POzZcuWtqMtqsRBuwWdBkVjhUZmfhu4EjirvvQMYHtm3jjvpR8DnhoR96+LkRcBH2kqhyRJWhsl3rUucdBeYkEnLaTprlMvBM6PiFcCU8DzASLi3VQbwC/KzK9HxDnAF+p/57NU3aokSS1zDb4GTYl7B/zvTYOi0UIjM7cBP7bA9RfMe/wu4F1N/t2SpAM3OTnJueeeu/vx+Pg4AOeccw6vetWrWkol7T8H7VJ7PBl8Faanp/e6y+cJv5IGTYl3fyVJ3WShsULT09McccQR7Ny5E6ju8k1MTLBjxw6LDUkDw7u/kqSmdLq9bZfs2rVrd5HRs3PnTnbt2tVSIkmSJKm7LDQkSZIkNc6lU5IktcAOX5IGnTMaKzQ8PMzExMRe1yYmJhgeHm4pkSSpZCWe7yBJq+GMxgqNjo6yY8cONm7cCMDWrVvtOiVJ2m92+JI06Cw0VmF0dJShoSEAxsbGWk4jSSqZS6Sk9nl0wdqy0JAkSdK649EFa889GpIkSVp3PLpg7VloSJIkSWqchYYkSZKkxlloSJIkad3x6IK152ZwSZIkrTseXbD2LDQkSZK0Lnl0wdpy6ZQkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxlloSJIkSWqchYYkSZKkxh3UdgBJGlQzMzPMzMwwOzsLwNTUFAAjIyOMjIy0GU2SpDXnjIYkrZHJyUnGx8fZtm0b27ZtY3x8nPHxcSYnJ9uOJknSmnNGQ5LWyJYtW9i8efM+153NkCStBxYakrRGXCIlSVrPXDolSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXGeDL5CMzMzzMzMMDs7C8DU1BTgyb+SJEnSQpzRWKHJyUnGx8fZtm0b27ZtY3x8nPHxcSYnJ9uOJkmSJHVOZGbbGVZsw4YNuX379lb+7t6MxnzOaEiSJJWnN7bbtGkTAFu3bgUc261WROzIzA0LPeeMxgqNjIwwNja2z48fREmSpPK4WmXtOaMhSZKkdcfVKs1YakbDzeCSJEladywo1p5LpyRJkiQ1zkJDkiRJUuMsNCRJkiQ1zkJDkiRJUuMsNCRJkiQ1zkJDkiRJUuMsNCRJkiQ1zkJDkiRJUuMsNCRJkiQ1rpFCIyLuFRF/EhH/GhE3RsRvL/K60Yi4MCJuiIhrIuJTEXFMExkkSZIkdUdTMxpnAQ8DHgpsAl4aEccv8tp3Asdl5gnAXwHvbiiDJEmSpI5oqtB4FvCuzJzNzJ3AnwPPnv+izJzOzL/NzKwvfRE4uqEMkiRJkjqiqULjKOCmOY+/UV9bztlUsxqSJEmSBshBK3lRRPwTcOwiTz9qf/7iiHglcAxwxhKv2QxsnnNpNiJu3Z+/r2GHAHe1HWKVzNwfZu4PM/eHmfunxNxm7g8z94eZ99+hiz2xokIjM39sqecj4mbggcA/1ZeOBm5e4vW/BzwdeGJmfm+Jv/c84LyVZOyniNiemRvazrEaZu4PM/eHmfvDzP1TYm4z94eZ+8PMa6OppVN/Cfz3iBiKiAmqPRt/vtAL61mKZwM/mZm3N/T3S5IkSeqQpgqNDwL/AnwNuBw4LzO/AhART42Id9d/3gC8EbgvcElEXB0RX2oogyRJkqSOWNHSqeVk5izwW4s8dxFwUf3n7UA08Xe2rHPLuVbAzP1h5v4wc3+YuX9KzG3m/jBzf5h5DcSeTrOSJEmS1Iymlk5JkiRJ0m4WGpIkSZIaZ6EhSZIkqXEWGpIkSZIa10jXqfUiIkYyc6btHJK0ViLi1AUu3w7ckJnTfY4zcBZ5f3fLzH/sV5bVqlvVXwx8JjO/3Xae1YiIBwBk5r+3nWU5EXEQcDZwTGb+RkQ8BHhgZn6m5WjSqtl1agUi4pHAh4D7ZuaGiHg08KzMfFnL0fYREf//Us9n5qv7lWWlIuI5Sz2fmR/oV5aVioj3AYv+x5OZv9rHOCsSEfewdOahPsZZsYh4EnBLZl4TEWcAPwFcl5kLHgradRHx5Mz8m7ZzLCYirgYeAXyd6vPyEGAbMA6clZmXtJduYYv83rsd+KfMvLzPcZYUEb08Q8BG9n6fr87ME1uKtqyI+BXgJ4HTgZ1URcfFmfm3rQZbQkT8KPBR4AH1pe3AL2bmv7SXamkR8Xaqz8fjMvNHI+K+VO/zSe0mW52IeFFmvr3tHAupC/6jqd7Xf59z/bmZ+f7Wgi0hIp4GZGZeFBGPA34RuDYz39NytCW5dGpl3gK8CPhO/fhK4Ofai7Ok+9Q/Pwr8NnAUcCTVOSf/tcVcS3lK/fPfgHcCzwF+BXgH8Mst5lrKFcCXgV3AY6kGC/8KbAK6Out1H6rB4jnAHwAPrH8mgSUL1LZExBuAPwQ+HBGvAN4MjAC/FxG/32q4/fenbQdYxpeBMzLz2Mx8KNWg8ovA04A3tJpscT8K/AbV77oNVL+vfwL4s4j4nRZz7SMzT87Mk4GrgZ/OzGMy81jgp6i+WzorMz+Ymc+hep/fAPwC8NftplrWnwKvy8z7Zeb9gNcBb2s503Iem5n/HZgGyMzbgYNbTbR/Xtl2gIVExGbgPVQD9Wsi4ulznj67nVRLi4jXUL2f50bEG6k+x/8OPDci/ner4ZbhjMYKRMQVmXlSRFyVmY+qr+3+cxdFxD8Az+tV6hFxOHB+Zv50u8kWFxEfB/53Zl5XPz4eeHVmPqPdZIuLiH8EnpyZU/XjMeBvMnPJ5RFtiogvZ+ajl7vWBRHxVaq7vj9EdSfygZn5HxHxQ8DWzDy+zXyLiYjFDlEK4Fczc7yfeVYjIq7JzBMWurbQc10QEX8PPDczv1U//hHgg8AvAZ/v4uckIq7NzEfOu9bJ97cnIl4APBE4EfgK8GngU5n5tVaDLSEirs7Mjctd65KI+GJmPrY3zoiIIarZrke0nW2+iPi/iz0F/FRm/lA/86xERFxLNVs0FREPAy6kGmtc0NWxXURcR/VdeG/gVuDIzPxuRNwH+ML83yVd4h6Nlbk7Ig6mXnYSEUcCs+1GWtYD5k4HZuYtEXFEm4FW4JhekQGQmddHxLFtBlqBQ3tFBkD9i+vQNgOtwH0i4rDeGuuIOIxqtqOLZjJzF7ArIm7PzP8AyMz/jIhdLWdbym9SzcQs9Hui63d37omIU3t7BeolBvfUz3U1+4ZekQGQmd+KiAdk5s6I+EGbwZYwGxGn9ZaiRcQT2PM+d9WfAluBzcCnM/P7LedZidmIeFhmfhWgHlh2/fv72og4C7hXRBwDvBz4bLuRFvXTwP+gmt2fK4DH9z3NCvW+tzPzqxFxOvCpuqDr6u+4XZl5NzAVETdm5ncBMvPOiOj059lCY2XeSlXxHhoRrwXOAjq3P2Oe7RFxLvDu+vGvUd0R7rKpiHge0Fsf+VzgrvbirMg1EXE+1TQswPOBa9qLsyJvpMrdW1f9JOBV7cVZ0m0R8dtUS77+IyJeTvX5eBLwn60mW9p1wF9m5lfmP1HfFe6y3wI+Ug/Qg+p74pci4hDgj1tNtrgdEXEO8N768fOBf+/4wGHu+wzV+/ysFvOsxA8Dp1Et83pNRNxOtcb9da2mWtorgX+s72IH8HCqZbpdtpnq9/T9gS9QjT9e3magJVwNXJWZV8x/ol7u00Wzc2+2Zeb2ev/fxVTLArto7h7KF/b+EBEBDPc/zsq5dGqFIuIUqjXKAVyUmZe2HGlJEXF/qr0lZ1B90V4M/I/MvLXVYEuIiOOoljtspMp8FdVyiG1t5lpKPfg6h2odO1Tv82sys9MFUkQ8nGrAAFUHmevbzLOY+m7eG6ju9P4u1eDsN6n2wzw3MztZ1EXEz1J1abpxgedO73r3mHoGt7ena1s9q9RZC/y++wzVXdadwLELFXxdMO99/pfM7Orsy2515sdRLaH6ZWAsM3+43VRLq2eZH1M//GJvZlQHLiJOAv59oW5eEfHQzLyhhVhLiohnATdn5j/Nu344cG5m/no7yRYXET9PtUzxP+ddP46qOVHnGv30WGioc+o1h2TmnW1nGVT1L9TjMvOzUbVSvFfXB5Pqj4g4Gbg+M78XEc+kanBw3kIDCR24iBihanAA7FnS0UUR8Ung0VR3sS+m2qNxRXZ4IFEvdf5WZu6KiB8HHgW8v8vfLxHxaeBTFPD+Ssux0FhCvTl5qXagT1/subZFxDhVV4IHZuZT6nWpJ2Tmh1uOtqjY0zv8IZn5m1FA7/D6S+xtVGvEN0bERuC0zOzqEhMi4kyqafl7MvNBEXECMJmZP9tytAVFxP2outscVV+6Gfh4Zt7WXqr9F91vb3sN1WbfBwN/S9Ua9MSON5I4F3hLb91yRPx/wG9l5rntJltcRDwWeB/w0LnXs6NtpgEi4iepNtcXc55KRFwJnEK17OuLwKXAQZn5i60GW0JEPJ6qjfATgWOoMl+cmZ3tWFf/nv4D9l5F8cou/542c39YaCwhIp671PPZ0V7LABHxEap14r+UmQ+PiP9C1VN+Y7vJFhcF9g6v9zl8CHhp3ZXnIKr1qp3rDtITEV+mWmN9ce7ponZ9RzvzPINqA+pngW/Ul48GnkA1kPxYK8EOQETcnJlHLf/KdkTElZl5YlRtYQ/KzPO62omlZ5HOQldmt8+k+BLwO8DbgVPrP09n5htbDbaMqDrrHc2cPZ6Z2dm2vHM+z78OHJaZr42Od/fqqW8Y/gLV8tzDM3O05UiLioi/otoH+s760guAozLzae2lWpqZ+8PN4EvociGxAg/NzF+qB2pk5vfrTUNd9th6VuAqqHqH1+uBu+ywuiXeSwAy8+6IuLvtUMuYrdvizb3W1WVTrwMek5nfmHsxIh4E/B3QyUIjlm5v29nWtrWRqNrDPoU9G1A7e5e9ttCZUJ3eIAkcnJlfioiD6mU8r4vqML/OFhoR8bvAq6nOlOp1uknmzcp0zEi9PO0ngTe1nGVF6qYzZ1AtqbuEam/a51oNtbyHzhvsvjgi/rm1NCtj5j6w0Fiheq3yRmD3HYXM3NxaoOXtNXCsZzS6XmjsNR1fd4zp+qGSd88t4Oppza6/z3fWA8leu+YzqDbNdtHQ/CIDIDP/rZ496qqS29v+MdVJ4Bdn5pX1EsbOTsvXtkXEy6gG6QG8BOjsyc+13sbv70bEicA3ga63xn4x1d6ukvbrfJjq3IEbgMvq/WnfazfSsl5AdQjsu6g2AO/TVKKD/j0iDs3M78DuDfg7Ws60HDP3QZe/qDsjIt4CPIhqE9yHqU6T/FSroZZ3SUT8T2A0Ip5I1bFnsYN1uqKk3uE9f0l1gvlYVG1LX8SelsJd9XKq2YAHR8SlVJ/trp50f3lEvJdqeclN9bUHUr3P+7RT7JBi29tm5rvZ+zP8b1R3g7vsbOAC4LVUhdw/Ar/SaqLlfSQifhj4faq71QcD/6vdSMvaUViRQb1U6q3AVGZmRNwJnNl2rqVk5v0j4pFUezTeHBFHA5dldVp4V90GfCUiPlE//lng873Z3Y7emDVzH7hHYwUi4ivACVRr70+oWym+v+ObIw8CXgr8PNUdvguB12dmZw92iapV7BupMkOVefP8dm5dExHPZs77nJkfajfR8uq1v6dQZb4sM29vN9HC6pm436M6X6C3r+Emqg3Kb8jMTt6ZjPLb2z4GeAh7r8P/QHuJViaqE+Pp+u+M+eoloqNd7oQEEBE/A/wM8DfMmYHO+nBHNSciNrBnQ/gZwK0d32N5zlLPd7Exg5n7w0JjBSLi8sw8OSKuBk7OzB9ExFe6vOFX/RERT8rMv1/umlSKiHgb1Wm/VzNnHX5mPrO1UOqEiHg11ab1r7P3Z2NTe6mWN78xQAGNArZR7TH6dO8n68PlpNJYaKxARHwGeDLVmutDqdZ7PjYzH7Pkv9iiiLiRPX3OP9Nr+9hlhbao3OcLq6tfYhHxucx8QkTcxt77BIJqsDDRUrRFRcSpSz3f5Tup9VrwZ1F16LkbuB74UGbOtJlrORHxNeARJbUwhfIGk1Be5oj4BrCxqzOgi4mIwzPzlsUed01EHJuZX2s7x2pExDuA38760Ml6v+IHMvMp7SZbnJn7wz0aK/Nsqrs3L6XaZHhfOr7Gk+rU55+kao335oi4lWpz58vajbWkp2Xm7mnBzPyPiHga0LlCIyIeSnWi73hEPHXOU+PAvdtJtaxfqv+5sc0Qq9TrwHMQ1fLFr1MVSQ+huuPeyUFZ3TziD4FrqJaoXQwcD5wTET+TmV3uEnIL0OliaBHz9xl1dd/RXKVlvqm0IgMgM2+JiJFekd/lIgOgtCKj9gPgSxHxi8CPAB+kak3eZWbuA2c0BlxEPIKq4HgxVU/8I1uOtKiIuDYzHznv2nWZ+fC2Mi0mqjNWngecxN6bkqeAd2bmJxb699pWd/L6ZGY+se0sq1FvCP9wZn6qfvxEqjNiOrmxut7XdVpdLD8YeGNm/kJE/DTwssw8o+WI+5hTMD+e6rC+P2fvdfgXtZFrNeYOJktRUuaI+COqvVIfpZDPRr2p+kPAfTNzQ0Q8GnhWx2+6FTfbBVAPft8B/CfwzMz8p5YjLcvMa89CYwWiOlX7VVQndM7dHPnIxf6dtkXEBcDJVC0eL6Za4/nVdlMtLSL+Euj1ke+1qDw5Mzs7exQRv5aZ72k7x2pExGVUhyLe03aWlVqkCO3soVsx7wC5uYOEiPjnzPzR1sItIiIuWeLpzMzT+xZmlUocTBaaeaHPSNc/G5+l6ub1J5n5qIgI4Lrs4AGlcxW43Ou+wPup9pY8GPhgZr621VDLMHN/WGisQERcC3wA2MqcvviZ+YXWQi0jIrZStUv8B6pWvJ/v+l2ziHgAVYvKxzGnRWVm3tpqsGXUHWMexN5nrFzbXqKlRcSbgWOp3uu7etc7flfyKqoOZJfUj58AvCk7elp1RPwDVSvsvwPOotrTdWY9yLkhM49tNeCAKXEwWWLmEkXEFZl5Usw53T46ftJ9T2GzXTcC78rM10fEfYD3APfLzM62xjZzf7hHY2VmM/OP2g6xGpm5qa58T6fap/GnEXFTVz+M9ZKeMzPz9CioRWVEPJnqUKX7UU1j3o+q/eqD2sy1jN7MwNye7Al0ttCgOhn3IxHxA6rZriGqjdZd9ZvA+cBbqJbWPbe+fijwBy1lGmSHZOalUZ+dmZkZEV097b6nxMwluru+GdQ7oPRIFj5IszPmznYBRcx2Ac/PzM8DZNWm+ZkR8VstZ1qOmfug66cud8Uly3W/6Zp64P5wqg20G4EJqoFwJ2V1vsdz6z//ZwlFRu01wGOBf87MHwaeQ7V+uZPqz8XHM/O0eT+dXfoAkJmXUW0A/3ngacAxmfnFVkMtITNvzMzHZeZ96vf35vr6t0tYahcRVy71uIOKG0xSZuYSPxtvpTqT6dCIeC3weapGDV32FqpDSb9TP76SjjcKyMzPR8QzIuKVABFxBNWqhM4yc3+4dGoFIuIU4JPAnVQb4HrtQB/carAlRMR3ga9Qtbe9GNiaHT6sDyAi/gD4Smb+WdtZVioivpyZj44556r0rrWdbTG9pQRt51iNiDgZuD4zv1d3dNoEnJcdPaV4uRsT2eG2vFDk+vCzqLoDPpJq/fJZVJvu/6LVYEsoMTOU99mA3d/hT6P67r4oMy9tOdKSSlzuFdUZKycDD8nMh0bV3vtjmXlKy9EWZeb+cOnUyrwPOJtqCUSnB+tzHJkdPTV5CS+kahf7HuB7dPh8hzl+UP9ze0T8AvANquVTXXZxRPy3kgo64N3AiRFxLPA6qlmj91EdLNdFRbbl7cny2oFeEBFfpxpMDgNndX0wWWJmKO+zAbtnRC9rO8cqlDjb9TSq32tXwO7PySHtRlqWmfvAQmNl7srM97YdYpV2RcRLqVrbQrUp/E2ZeXeLmZazse0A++HNUR2Y87+Aj1CtqT271UTLK7Ggm83M2Yj4GeBtmXlevUG8kzLzZNjdlvdlOa8tb5vZllPo+vASB5PFZS7psxERH2fvg0n3kplP72Oc1Zq/3OssoHPv8Tzfr39Hz70Wi724I8zcBxYaK/OJiHhKZv5120FW4Tyqu6d/SvXL9gVU/c9/p81QS8nMmyLi3uwpOK7u8qxMvd9hV2beBnyZqpNTCTa2HWA/jETEjwBPAV5eXxtqMc9KnZSZv9p7kJkXR8Qbl/oXOqC3PvxP6sdXUnXd69xAp8TBZImZ5yjms0E1UC9SobNdN0XE44GsZ2NeSTV722Vm7gMLjZV5MdUd4O9TnZhbwh3gnwA2Zn1WQkR8gupLobPqdbQfA3rtbH8kIp6RHT2Mpr6r8D+pMhcjM29qO8N++GNgG9Xp9ldGxEOA21rOtBKzEXFa7t2Wt+vnl5TUDenCtgPshwvbDnAAivlsZOb7285wIEqb7aK6ifl+4BFUjWcuAf5bq4mWZ+Y+sNBYmY1tB9gPQdVV7J45jzs9vUY1C3Nm1ueT1IXHH1N1deqqKyPicQXcbdotIg4DzqXaOzD37I/O7hvIzHdT7dPo+Tf2LAvsstLa8kJB68NLHEyWmHmOYj4bc9UNJDay9++7za0FWkTJs12Z+S3gSfWqhCihc6SZ+8NCYwXqJT2HA8dl5mcj4iC63xr474F/iIjz68fPoTo8rMv+S845BDEzL4uI0aX+hQ54LPC8epp77uF3nR20Ux3wcylwBtXp6y8EOrvfoSciHkO1HHDu760PtBRnRerP8EOA/1pf+pfM/MFS/04HlLg+vJjB5FwFZi7usxERb6E61+jRVIdo/iLVIbZddGHbAQ5Ul5c7L8bMa8v2tisQEWdSdZHJzDw6Ik4AJjPzZ1uOtqiIuBfVAPKM+tLFwDt7S6m6KCK+AJyTmRfXj88AXp2ZP95ussXVS2H2kZmf63eWlYqIqzNzY68lb0QMA5/LzB9rO9tiIuJtVB2mrmbPHdTMzGe2FmoFSmvL21NgO9AFB5OZ+WutBltCiZmhyM/GV6hmb6/KzBMi4v7A+zOzqx3rihURV869yTb/cReZee1ZaKxARHwZ+Cmq9eG9ntbXZ+bx7SYbLBFxEtV+h95A8l7A0zOz03tLShMRW7M6Of5y4ElUex22ZWZnN7NHxNeAR2TmdNtZViMirqFqRfhg4G+p2vKe6CCnWSUOJkvMXKKIuDwzT46Iq4GTM/MHMefco64qcLar1DNWzLzGXDq1MrOZ+d157cQ6uQEuIt7H0ms8f3Wx59qWmVdExDHAcfWlbV1fZhIR/4WqWcBG9v5C6OxaWuCGiPhh4ALgS8AUVdesLruFqhFDaYppy1vy+nBgOjPviYiMiIMz89aIeEDboZZRTObCPxt31uvZLwUuiIhbqdp6d1Zhy712yzLPWDHzGrPQWJk769aavQ1wZwA72420qCvqfz4COJWq53lSnUD7+bZCrUTd9vM9mXld21lW4V1UA/VTqJbXPQ/o9KnPmXlW/cc3R8QVVAcM/n2LkRYVEU+t//gl4KMR8efA7lmNzLyolWArV1Jb3gvbDnAAihtMUlbmC9sOcACeTTVL/lKqPWn3Bc5sM9AKnMae2a6XRMQbqDoNdVYUdMZKj5n7w6VTS4iI4zPz+npJzzuplj9cR3Wn4ecy8+o28y0lIv4ReHJmTtWPx4C/ycxT2022uIg4B3gu8B9Upz5/KDPvaDfV0ubsc7g2Mx8ZEfcBPtHl97mn/kzsvtmQmZ0rniPikiWezsw8vW9h9kNEvAD4I6pll2fWG8Pfm5kL7u3R/qmLudupllv2BpNvzsxvthhrSSVmVn+UuNwrIj5LdXDtn2Tmo6JaAnJdl5eYm7k/nNFY2gep1le/ieoOwylUG+Auy8zb24u1Iof2igyAzJyKiEPbDLSczDwXODciTqMqOM6NiE9lZpd7RH+//ufdEfFDmXln19/niHgW1YFb96Oa7Yr6n8Nt5lpIZp7WdoYDkYW25S1tfXhWLR97XttakFUoMTOU99mIiIcBrwKOYe8bK49sK9MKlDTb1VPMGStzmLkPLDSWNloPyg6nWobU26RxakR0fdnGNVG1tn1P/fj5wDXtxVm5zLwkIu6iutP3LLp9GM3OiLgf1UbfT0bEfwDbW860nEngZzPzimVfqQMWhbXlLXF9eImDyUIzF/fZAD5C9d/bWyngzI9aicu9Sjxjxcx94NKpJdTrw18EPJ49ex96Or1sIyIOAc4BehkvBl6TmXct/m+1K6qD5H6FqigKquVTF2TmrUv+iy2KiKF6s28Av0w1S/CBubNJXRMRl2XmKW3nWK0orKUflNmWt8RuSBFxLdVgcitzvnRzzrk8XVNo5hI/G1f1ukVq7UTEWVQF0iOp9pOcBbwsM/+i1WBLMHN/WGisQES8OTPPbjvHoKtnAz4KnJ+ZX2w7z0pFxA8Bj6K6w3B1dvykzoh4HvAAqvd67sbqm9vKtBJRWEs/KLMtb6Hrw4sbTBaaucTPxnnAhZnZ6SYdc5U42wXlnbECZu4Hl06tQElFRkQ8OzM/HBG/s9DzmfmWfmdahSMz8/vLv6w76g5kHwJ6B7AdXv//YKlNzG0bodpM9nvMucsOHNZaohUoraVfrcS2vCWuD78kIk4taTBJmZlL/Gx8lGpZ651UN1aCalbxwe3GWlKJy73IzMuAy9rOsRpmXnvOaAyYiDg3M8+J6jyN+TI7eI5GycVRvZTgBZn5pfrxJqoWvV2+w3cTcHpm/mvbWVZqbku/zOx8S7/Y05b38VTd6oppy1tiN6T6Dt8ngWIGk4VmLvGzsQ14PdXy57lL1K5vLdQySprtigLPWDFzfzmjMWAy85z6n89vO8sq/Nf6nwv9Yu16JXxPr8gAyMytEdH1O1DbSyoyam+h2i/1J/XjK6nu+HWy0AB+d97j35jz5wQ6W2gU2g3pfcDZzBtMdlxxmQv9bNyVme9tO8QqlTTbdWHbAfbDhW0H2A8Xth1gfzmjMaAi4tNU3UA+DVyR/j96TUR1kNL17DlM6VeAh3f1TjtARLwGuDf73mW/trVQy4iIKzLzpLl3+kq661eSEteHR8SXM/PRbedYjUIzl/jZeDVweWb+ddtZVqrE2S5pMRYaAyoiHk/Vr/+JVF8Kl1IdGvanrQZbQEQ8Z6nnM7PLrUBvA8aBH9SXDgZ6hwxmZk60EmwJEfFvC1zu9JdYRHyRahnSlzLzxLql38cz86SWow2cQrshlTiYLDFziZ+N3u/o71Ptl+oN2jv3u7mnxOVeUN4ZK2DmfrDQGHARMQ78AlWr28Mzc3SZf6XvIuIv6z+OAU+gKooSeBzwucx8UlvZlhMRD1zq+cy8qV9ZViIihoBfyMyPtp1lNUps6ddTWlveEmeKCh1Mlpi5xM/Ggr+ju/a7ea5CZ7sWPGMlM3+t1WBLMHN/WGgMqIh4LXAGVYehS6iWUH2uy61X681O/zszr6sfHw+8OjOf0W6ypUXE4cBxmfnZiDgIuFdmdvakzq4PdBdTWku/ntLa8hbaDrTEwWSJmYv7bECRv6NLnO0q8YwVM/eBm8EH1wuArwPvoqp2b2w5z0oc0ysyoJomjohj2wy0nIg4E3gj1QzM0cDx1CdvtxhrOVdGxONKGaj3lNbSryfLa8tbXDvQzLxpocFk27mWUmJmCvxsFPo7+sXAeEQUM9sFTGfmPRGREXFwZt4aEQ9oO9QyzNwHFhoDKjPvX7cEfSLw5og4GrgsM/97u8mWNBXVYXK9jdXPBTp7knltC3Ai1cnrZOY1yy2n6oDHAs+LiK8z5/3t4ixHyS39YO+2vEDn2/LWiuuGVOJgssTMFPjZoMzf0RvbDrAfSjxjxcx9YKEx2HYCt1H1Pf9h4ORW0yzvV4EPAu+oH19FVWx02Wxmfjci5l7r7JR87bfaDrAKF7Yd4ACV1pYXymwHWuJgssTMJX42ivsdXehs17Opis+XsueMlTPbDLQCZu4DC40BVXetGKbam/E3wO9m5rfbTbW0zNwGbIqI+9SP72w50krcWR9ilbD7pPCd7UZaWmZ+LiIOBo7Kjp+nkZnvX/5VnXZIZl7aG+RkZkZEpwc5wCci4iklrQ+nwMEkZWYu8bNR3O/oEme7ssAzVszcHxYag+vJmfm1tkPsj0IKjJ6XA38HPDgiLqXqBvFz7UZaWkT8BNVynruBoyLiZODszDyrzVzLKa2lX+3uuqjrDXKOpPtLTkpcH17cYJIyMxfz2YiI4+t2sMX9jqbA2a4o84wVM/eBhcaAKrXIKK0VaGZeERGnAadQfelelpm3t5tqWX9AdSbFRwEy8/KI6HTLysVa+rUaamXeSrX869C6E9xZdHvZFBS0PrzEwWSJmefY2HaAVfgg1WD9TUBpv6NLnO36CNWy0LfS/ZspPWbuAwuNAVbaoL02/4u201+89R3qb2Xm30XEjwNnRcT7Oz4rM5SZ/1rYl9hp7Gnp95KoTmTv/LKqzLyg3nT/NKqljGd1vdtXYevDSxxMlpgZKO6zMRoRzwIOB06leo8BTo0IMvOi9qItq8TZrtnM/KO2Q6ySmfvAQmOwFTVohyJbgf4VcEpEHEF1p+FSqkMHf7HVVEubjohD2PMl9giqQ8O6rLiWfj2lteUtbH14iYPJEjMDxX02XkHViOEwYP4SywQ69z4XPtt1SUScWtgZK2buAw/sG3BzB+1dN7cVaGYW0Qq0N0sUEb8OHJaZr42IazLzhLazLSYifopqjedDqJYfPRH45cz8TJu5lhIRnwGeDPwhcChwK/DYzHxMq8EWUXJb3oj4MvBTwMVZnwIdEddn5vHtJttXRDyVajD5eKqWq3NlZp7e/1RLKzFzT0mfjZ6IeHNmnt12jpWY831yKVVhUcxsV1QHqn4SKOmMFTP3gTMaA6rQ/v0ltgIdiYgR4CeplkJ0Xmb+Q0R8DXgS1S+pc7refYryWvpd2HaAA1DM+vD67v9FJQ0mS8w8RzGfjZ7C3uNiZ7so84wVM/eBhcbgKnHQXmIr0A9T3V2/AbisXr/c6cNzIuIpwCcy821tZ1mp0lr6Fd6Wt7j14YUNJoEyM1PgZ6MwxS33mqPEM1bM3AcWGoOrxEF7ca1A66VSbwWm6vf4Trp9px2qL7C3R8SfAe/LzH9uO9BySmzp11NKW97C14drDfnZ6I/CZ7tKPGPFzH3gHo0BFRFfpFoD/KV6zeeRwMcz86SWoy0qIs6iWiLzSKqOQmcBL8vMv2g12ACKiAcDz6l/vgW8NzPf1W6qxUXEtVQzcluZU3xm5hdaC7UCi7XlzcxfazXYAkpeH6615WdDy4mI24BxqsYinT5jpcfM/WGhMaBKHbTXG52eRvUfz0VdbwUKxbYRBiAiRoE/Bn49M4fazrOYiLiqt/m0JBHxFfa05T0hIu4PvD8zf7rlaPuIiK8C5wK/D/wP9qwPB+j6+nCtIT8bWk4scqBgZt7U7ywrZeb+sNAYYCUO2ksUEYfPbcM7/3EXRcSJwPOBZwKXA+dn5kfbTbW4iDgPuLCkln4AEXF5Zp4cEVcDJ2fmDyLiK5n5iLazzVdyNyStLT8bWomFzljJzE4v2Tbz2rPQUOtKbgXaU1gb4WupDo87n+rueqeLIiizpR+U15YXymoHqv7ys6HFzD1jJTOPjogTgMnM7OIZK4CZ+8VCY8CUOGiPiOcu9XyXO/gUevbHKVkdIleMiNgGvJ55Lf3qDaqdVXfouZ3q9OReW943Z+Y3W4wlSY0q9IwVM/eBXacGz4VtB1itLhcSK1BiG+GtEfES4JjM/I2IeAjwwOzwgX0U2NIPymvLK0n7qbgzVjBzX1hoDJjCB+3FtAKdo8Q2wm8FhoDH1Y+/C/w50NmOZBTY0g/KbssrSatQ4hkrZu4DC40BVtqgfbFWoK2GWl5xZ39Q7RHYGBFXAWTm7fX/DV32YmA8Iopp6Vf7CNUM11vp/udCklalxDNWzNxfFhoDqtBB+2nsaQX6koh4A1Vr3i57K9VytUMj4rXUbYRbTbS86bkPImKIag9Bl21sO8B+ms3MP2o7hCStkQ8CJwJvovoOL+GMFTP3UdcHF9p/p1G1tv1OZr4E2ARsaDfSsqYz8x4gI+LgzLwVeEDboZaSmRcAr6PaED4MnNX1s0qAa+tzVu4VEccAbwc+226kpdU9wncBD6r/vAPofLcs4JKIOLXtEJK0RkYj4lnA4cCpwMFUN7FPrdsid5GZ+8gZjcE1nZn3RMTuQXtEdHrQTrX28N7ApcAFEXEr8L2WMy2r7uBUUhenzVTt8e4PfIFqRuYVbQZaztyWfsDRwPHAJNDZln61jwKfjIii2vJK0gq9gqohymFU3y1zJdDFwxzN3Ee2tx1QhfbvL6YVaIlthEtWYks/KLctryStRolnrJi5Pyw0BlRJg/YSlXz2B0BEPAZ4CHt3QvpAe4mWFhFbM3NTRFw1p9DY/eeuiogvZ+aj284hSVIbXDo1oErs319SK9CuFxJLiYi3AT8NXM2eu+xJ1R2pq4pr6Vcrsi2vJElNcEZjQJU0aO+JiGupBrtb2XuZyRdaC7UCBbYR/hrwiMycXvbFLeu19IuIk4B3Ag8GrqNu6ZeZV7eZbzkRcRswDpTWlleSpAPmjMbgKrF/f3GtQAttI3wL1aC3BMW29KttbDuAJEltsdAYXMUN2qlbgWbmP7YdZBVKPPvjS8BHI+LPmXOmRmZ2sWvF/JZ+UV8/NSK6mnm3zLwpIg4HjsvMz0bEQdhWXJK0TlhoDK4SB+0ltgItsY3wSfU/f2POta62xyu2pR8U3ZZXkqQD5h6NARURpwCfBIoZtJfYCrTENsIlKrGlH5TblleSpCY4ozG43geczbxBe8fdlZnvbTvEKj2b6v19KXvaCJ/ZZqBBVGKRUZvNzO9GxNxru9oKI0lSP1loDK4SB+3FtQItsY0wQERcmZknLvZYjSm1La8kSQfMpVMDKiJeDVxe0qC9xFagJbYRBoiIwzPzlsUe68CU3pZXkqQmWGgMqEIH7Q9c6Hpm3tTvLCtV6tkfABExkpmltLktSm+GKCIuBX6O8trySpJ0wFw6Nbg2th1gtQptBVpcG+GIeCTwIar9JBsi4tHAszLzZa0GGyxFt+WVJKkJzmgMsIUG7ZnZ2Y2oc1uBZubREXECMJmZnW0FGhHnAReW1EY4Ij4L/C/gTzLzUVHtVL7OTkjNiYinUrXlfTxVQ4a5MjNP738qSZL6yxmNAVVo//4tVKdAXwyQmdcstpyqQ0o8++OQzLy01wkpMzMiOluAlqiesbio1La8kiQ1wUJjcJU4aC+xFWiJbYTvjoiD2dMJ6UjKyV4UiwxJ0npmoTG4Shy0l9gKtMQ2wm8FLgQOjYjXAmcB7s+QJEmNstAYXMUM2nutQIGXA38HPLju1vMgqo49XVbi2R8XRMTXgacBw8BZmXlpy7EkSdKAcTP4gCmxf3/JrUBLbCMsSZLUD85oDJ4PUu3NeBNwGmUM2ktuBbqx7QArFREfp57hWkhmPr2PcSRJ0oCz0Bg8JQ7aX0HVCvQwYPO85xLoYmaguLM/Lmw7gCRJWj9cOjVgSu7fX2Ir0BLP/pAkSeoHC40BVeKgvUQR8WXgp4CLM/NR9bXru374XUQ8k2rZ12jvWmbOn02SJEnab11d4qEDZJHRN7OZ+d151zrdRjgi3gL8CvA8qqVpZ1JtaJckSWqMhYZ0YIppIzzHaVStbb+TmS8BNgEb2o0kSZIGjZvBpf1Q+Nkf05l5T0RkRBycmbdGxAPaDiVJkgaLhYa0f0psI9xzZ0TcG7gUuCAibgW+13ImSZI0YFw6Je2f+W2ED6Yq3E+tO3912bOBWeClwFeAH1Dt05AkSWqMXaek/VByG2FJkqR+sNCQDkCJbYQj4mHAq4BjmLN8MjMf2VYmSZI0eCw0pHUmIq4FPgBspVpCBUBmfqG1UJIkaeBYaEjrTERc1TtcUJIkaa24GVxafy6JiFPbDiFJkgabMxrSOhMRpwCfBO4Epqna8mZmPrjVYJIkaaB4joa0/rwPOJuqW9bsMq+VJEnaL85oSOtMRHw5Mx/ddg5JkjTY3KMhrT+fiIintB1CkiQNNmc0pHUmIm4DxoHvAzPs2aMx0WowSZI0UNyjIa0/G9sOIEmSBp8zGtI6FBGHA8dl5mcj4iDgXpm5q+1ckiRpcLhHQ1pnIuJM4IvA+fWl44EL28ojSZIGk4WGtP5sAU4EbgPIzGuAB7aaSJIkDRwLDWn9mc3M78675rIpSZLUKAsNaf25MyJ+BEiAiDgD2NluJEmSNGjcDC6tExFxfGZeHxEnAe8EHgxcBzwI+LnMvLrNfJIkabBYaEjrRERcmZknRsSlwM8Bp1CdoXFZZt7eajhJkjRwPEdDWj9GI+JZwOHAqVRFBsCpEUFmXtReNEmSNGic0ZDWiYh4KvAi4PHAFfOezsw8vf+pJEnSoLLQkNaZiHhzZp7ddg5JkjTYLDQkSZIkNc72tpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXEWGpIkSZIaZ6EhSZIkqXH/Dy3mbbhDtxlbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 960x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = joint_CI.index\n",
    "\n",
    "coef = joint_CI.iloc[ : , 1 ].to_numpy()\n",
    "\n",
    "sd_error = size_err\n",
    "\n",
    "figure(figsize=(12, 6), dpi=80)\n",
    "\n",
    "plt.errorbar( x = x , y = coef , yerr = sd_error , linestyle=\"None\" , color = \"black\", \\\n",
    "              capsize = 3 , marker = \"s\" , markersize = 3 , mfc = \"black\" , mec = \"black\" )\n",
    "plt.xticks(x, x, rotation=90)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}