{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"id": "b2ebc3f6-505f-4c2c-9463-5d34b01c6a86",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import scipy as sp\n",
"import scipy.stats\n",
"import statsmodels.api as sm\n",
"from sklearn.linear_model import LogisticRegression"
]
},
{
"cell_type": "markdown",
"id": "3ffce874-3d61-4dad-a179-aa1f2fb7db18",
"metadata": {},
"source": [
"# What Is Qualitative Response Model?"
]
},
{
"cell_type": "markdown",
"id": "76c0fa22-cd7e-4727-9fb2-df9ece1cf899",
"metadata": {},
"source": [
"So far we have only discussed the model that takes $Y$ as a quantitative variable by default, however there are situations we want $Y$ to be qualitative variable, for example, $Y = 1$ represents a family owns house, $Y=0$ means not owning a flat, independent variables can be quantitative variables such as their income, ages, education and etc.\n",
"\n",
"Specifically, if the response variables only takes two values, such as $1$ or $0$, we call them **binary variable**, the regression model with binary variable as dependent variable is called **binary response regression model**, similarly there is also **trichotomous response regression model** or more generally **polychotomous response regression model**, but names do not matter, they are all **qualitative response model** as in the question of this section.\n",
"\n",
"For easy demonstration, we will mostly deal with _binary response model_, there are four approaches to develop this type of model: **linear probability model**, **logit model**, **probit model** and **tobit model**. \n",
"\n",
"We will go through all of them in this chapter."
]
},
{
"cell_type": "markdown",
"id": "e6df86b8-ba84-4bea-b091-d680f122b7f7",
"metadata": {},
"source": [
"# The Linear Probability Model (LPM)"
]
},
{
"cell_type": "markdown",
"id": "db69ebdf-e008-48bd-98fd-4811f73285a4",
"metadata": {},
"source": [
"Consider a simple regression model\n",
"$$\n",
"Y_{i}=\\beta_{1}+\\beta_{2} X_{i}+u_{i}\n",
"$$\n",
"where $X$ represents family income, $Y=1$ represents if the family owns a flat and $0$ the opposite. There are only two outcomes, either own a flat or not, so $Y_i$ follows **Bernoulli Distribution**.\n",
"\n",
"It would be fast to recognize the $u_i$ can't be normally distributed, since\n",
"$$\n",
"u_{i}=Y_{i}-\\beta_{1}-\\beta_{2} X_{i}\n",
"$$\n",
"If we denote the probability $Y_i=1$ as $P_i$ and $Y_i=0$ as $1-P_i$, then\n",
"$$\n",
"\\begin{array}{lcc} \n",
"& u_{i} & \\text { Probability } \\\\\n",
"\\text { When } Y_{i}=1 & 1-\\beta_{1}-\\beta_{2} X_{i} & P_{i} \\\\\n",
"\\text { When } Y_{i}=0 & -\\beta_{1}-\\beta_{2} X_{i} & \\left(1-P_{i}\\right)\n",
"\\end{array}\n",
"$$\n",
"it shows that $u_i$ also follows Bernoulli distribution."
]
},
{
"cell_type": "markdown",
"id": "10aecef9-9e8b-4d9c-a3ff-3e055cd5293a",
"metadata": {},
"source": [
"Recall the fact from statistics course, that the variance of Bernoulli distribution is $P_i(1-P_i)$, what's more\n",
"$$\n",
"E\\left(Y_{i} \\mid X_{i}\\right)=\\beta_{1}+\\beta_{2} X_{i}=0\\left(1-P_{i}\\right)+1\\left(P_{i}\\right)=P_{i}\n",
"$$\n",
"The $P_i$ is a function of $X_i$, hence $u_i$ is heteroskedastic. "
]
},
{
"cell_type": "markdown",
"id": "5f55d5c0-9b6a-4472-b4af-37024a001aa7",
"metadata": {},
"source": [
"## Home Owners' Example"
]
},
{
"cell_type": "markdown",
"id": "69322a7b-7474-49f7-a3d1-c537a52d969c",
"metadata": {},
"source": [
"Here is an excerpt data from a real estate survey in Manchester, UK. The family income is measured in thousand pounds."
]
},
{
"cell_type": "code",
"execution_count": 84,
"id": "445a0ab0-9f82-4655-a3b2-ed9ac4bc6372",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.scatter(df[\"Annual_Family_Income\"], df[\"Own_House\"])\n",
"ax.plot(df[\"Annual_Family_Income\"], model.fittedvalues, color=\"tomato\")\n",
"ax.grid()\n",
"ax.set_xlabel(\"Income\")\n",
"ax.set_ylabel(\"Owning a flat = 1, Not Owning = 0\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "581c58bb-6043-4cc6-8158-cffbc5e94d2a",
"metadata": {},
"source": [
"Let's interpret the results, the constant term is $-0.1981$, but probability can never be negative, hence we ignore constant term, or simply treat it as $0$. The slope term $0.0162$ indicates that in average with every $1000$ pounds increment in family income, the probability of owning a house will increase $1.6\\%$.\n",
"\n",
"However, you could doubt its reliability that assuming income and probability is linear relation, actually it's more rational to assume there is some threshold to afford a house, if the house cost $500000$ pounds, earning $20000$ shouldn't be any obvious improvement than earning $10000$. Similarly, family with income of $200000$ wouldn't be much different than family with income of $190000$. We will talk about how to model this assumption in next section. So far we will stick with LPM.\n",
"\n",
"What if a family has $39000$ pounds of annual income? Our model predicts the probability they own a house is $-0.1981 + 39\\times0.0162=.43$. But how about family income of $11000$ pounds? The models says $-0.1981 + 11\\times0.0162=-0.0199$, it predicts a negative probability and apparent it has no sense. Similarly, if a family has $99000$ pounds income, the model predicts a probability $-0.1981 + 99\\times0.0162=1.4057$ which is more than $1$.\n",
"\n",
"If the vertical axis represents probabilities, i.e. $0 \\leq E\\left(Y_{i} \\mid X_{i}\\right) \\leq 1 $, it would be a natural defects of LPM to predict out of this range."
]
},
{
"cell_type": "markdown",
"id": "d320a657-2c4c-4a48-afcc-c3183b149d06",
"metadata": {},
"source": [
"However, we still haven't address the heteroskedasticity issue of LPM, which means the estimated results above are actually all _invalid_! For fast remedy, you can invoke the White's robust standard error."
]
},
{
"cell_type": "code",
"execution_count": 87,
"id": "3838467c-88c9-4aa5-9373-d2fe5df51dff",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"const 0.086859\n",
"Annual_Family_Income 0.001809\n",
"dtype: float64"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.HC0_se"
]
},
{
"cell_type": "markdown",
"id": "2fe93e84-b8b4-4386-a647-0e05c3dbe5e3",
"metadata": {},
"source": [
"To handle the issue more seriously, we use WLS to obtain more efficient estimates. Because variance of disturbance term is \n",
"$$\n",
"\\operatorname{var}\\left(u_{i}\\right)=P_{i}\\left(1-P_{i}\\right)=\\sqrt{E\\left(Y_{i} \\mid X_{i}\\right)\\left[1-E\\left(Y_{i} \\mid X_{i}\\right)\\right]}\n",
"$$\n",
"So we estimate the disturbance term by using fitted values, note that these are different than residuals!\n",
"$$\n",
"w_i=\\hat{Y}_{i}\\left(1-\\hat{Y}_{i}\\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "a1b66e3a-d9b8-4621-aab9-e8ce727a7c23",
"metadata": {},
"source": [
"Add fitted value into the data frame."
]
},
{
"cell_type": "code",
"execution_count": 71,
"id": "3439f127-c5be-4803-bf14-ef703b3f016b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" Annual_Family_Income Own_House fitted_value weight\n",
"0 56 1 0.710644 0.453463\n",
"2 67 0 0.889155 0.313940\n",
"3 34 1 0.353621 0.478093\n",
"4 23 0 0.175110 0.380061\n",
"5 19 0 0.110197 0.313135"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.head()"
]
},
{
"cell_type": "markdown",
"id": "e3c62153-8caa-4096-96bf-d3c217c423d9",
"metadata": {},
"source": [
"Divide weight onto every terms in the model\n",
"$$\n",
"\\frac{Y_{i}}{\\sqrt{w_{i}}}=\\frac{\\beta_{1}}{\\sqrt{w_{i}}}+\\beta_{2} \\frac{X_{i}}{\\sqrt{w_{i}}}+\\frac{u_{i}}{\\sqrt{w_{i}}}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "42993b97-86bf-414b-ba27-23b3e0c2511c",
"metadata": {},
"outputs": [],
"source": [
"df[\"Annual_Family_Income_WLS\"] = df[\"Annual_Family_Income\"] / df[\"weight\"]\n",
"df[\"Own_House_WLS\"] = df[\"Own_House\"] / df[\"weight\"]"
]
},
{
"cell_type": "markdown",
"id": "692b307d-132f-477b-a98b-42df52903626",
"metadata": {},
"source": [
"Here's the more reliable result."
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "2e401f00-aff8-496d-b279-b668be53f8f4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Own_House_WLS R-squared: 0.765\n",
"Model: OLS Adj. R-squared: 0.757\n",
"Method: Least Squares F-statistic: 94.47\n",
"Date: Thu, 02 Sep 2021 Prob (F-statistic): 1.25e-10\n",
"Time: 14:17:57 Log-Likelihood: -40.040\n",
"No. Observations: 31 AIC: 84.08\n",
"Df Residuals: 29 BIC: 86.95\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"============================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"--------------------------------------------------------------------------------------------\n",
"const -0.3322 0.245 -1.358 0.185 -0.833 0.168\n",
"Annual_Family_Income_WLS 0.0146 0.002 9.720 0.000 0.012 0.018\n",
"==============================================================================\n",
"Omnibus: 5.495 Durbin-Watson: 2.649\n",
"Prob(Omnibus): 0.064 Jarque-Bera (JB): 3.808\n",
"Skew: -0.720 Prob(JB): 0.149\n",
"Kurtosis: 3.934 Cond. No. 244.\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"X = df[\"Annual_Family_Income_WLS\"]\n",
"Y = df[\"Own_House_WLS\"]\n",
"\n",
"X = sm.add_constant(X) # adding a constant\n",
"\n",
"model = sm.OLS(Y, X).fit()\n",
"\n",
"print_model = model.summary()\n",
"print(print_model)"
]
},
{
"cell_type": "markdown",
"id": "ec653d5c-4c76-4cb0-8ac6-59000227bdfc",
"metadata": {},
"source": [
"The slope estimates of WLS is similar but higher significance due to smaller standard error and $R^2$ also much higher than previous model."
]
},
{
"cell_type": "markdown",
"id": "2575796b-11f4-4596-bb76-fbda483b45cc",
"metadata": {},
"source": [
"Though we have fixed heteroskedasticity to some extent, the model is not attractive in either breaking probabilistic range or its unrealistic linear shape. In Next section we will discuss about alternatives to LPM that fulfill both probabilistic range and nonlinear shape."
]
},
{
"cell_type": "markdown",
"id": "e5d019cb-9af3-452a-a812-0077bd630f6b",
"metadata": {},
"source": [
"# The Logit Model"
]
},
{
"cell_type": "markdown",
"id": "e8a726a4-d3a4-4297-a780-b2c0fdb5dc65",
"metadata": {},
"source": [
"The **logit model** hypothesizes a sigmoid shape function \n",
"$$\n",
"P_i = \\frac{1}{1+e^{-(\\beta_1+\\beta_2X_i)}}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "044be4dd-feb1-484e-84c1-aa52d2d27ddf",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAGvCAYAAABRkLPAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4+UlEQVR4nO3deXxU9d3+/+s9k8mesBOWsBM2EVQQcI87WK1d1LpUK9a6tFrb21837+79trba9u7iVm632lppa7VFb5S6EEWtbCqyQwhb2CGE7MvMfH5/TMCIICHM5EzmvJ4P5zFny5wLP5JcfnLmjDnnBAAAAPhNwOsAAAAAgBcowgAAAPAlijAAAAB8iSIMAAAAX6IIAwAAwJcowgAAAPAlijAAAAB8Kc3rAACAzsXM5kjKaFntL6nSOXeyh5EAoF2MD9QAALSHmRVIeknSNc65pV7nAYCjxaURAJBEzOxuM/ua1zmOxMyyJT0r6Rv7S7CZLTCz47xNBgBtRxEGAElm1s3MnJnVmFmdmW091kJqZreZ2SIzazSzxw9zTD8zK29Z7iXpOkl/OJrX6GhmFpD0Z0lPOOfmtNr1S0k/9iYVABw9rhEGgJgTJO1yzvWWJDO7VNI/zexp51x5O19zq6T/J+lCSVmHOeYiSS+2LF8vabZzrv4oX+NDzOyHkuSc++HRBm75+j6Snj7Ersucc9sl3SOpzDn30EH7Z0l6yMz6Oue2tefcANCRKMIAEHOCpIWt1ue3PKe39wWdc89IkplNlFR4mMMuUmx2VZKmSXq0Ha/RLmb2JUlfV+wNb29J+oJzbmdL2T39MF9zs6Shki47eJ9zrsHMFku6QNIf45kVABKBSyMAIOZESQskycy6SvqppMWS1u8/wMyeN7PKwzyeP9oTmllI0pmKveFMko6XtPoY/xxtPfddkm6R9ElJvSRtUWzm+eO+po+k+xQrzq+aWYmZzTzosJWSxsc/MQDEHzPCABBzgqTPmNkdkvZKmivpEtfq1jrOuYvjfM4zJS1xzlW3rHeVVH34w+PDzHpL+q6kE5xzpS3bHpH0wMd9XctMcegIL18tqW88cgJAolGEAfiemWVIGi1pyDFcD9weF0ma3Wp9r6S89rxQy4z0/ssZMlu2fa1l/Y2DSvy5il3yscDMDryEpHfbc+6D5EmqjMPrAEDCUYQBQBorqfZIJdjMXpB0xmF2z3POTTvK814k6dOt1t+XNEIfvla5TVoX3Ta8Wa67pGedc5cf7XnaYLQ+uOYZAJIa1wgDQOz64OVHOsg5N805l3uYx0dKsJmlmVmmpKCkoJllmllay74hkjKcc6tafclsSWe19TWOwTuSzjazk1rOkW9ml1qr6eH2aJlZn6APrnkGgKRGEQaA2PXByxLwut+VVC/p25I+37L83ZZ9n9CHL4uQpCckXWRmrW+T9nGv0S7Ouf8odr/ff5hZjaQVkqa2vh66nT4pqcQ5t/UYXwcAOgQfsQwAHjCz2ZLuc87NPmj7zyTtdM79xpNgx8DM5kv6onMuEf9TAQBxRxEGAA+Y2Tcl/f6gD88AAHQgijAAAAB8iWuEAQAA4EsUYQAAAPgSRRgAAAC+5NkHavTs2dMNHjzYk3PX1tYqJyfHk3OjYzDG/sA4+wPj7A+Mc+rzcowXL1682znX6+DtnhXhwYMHa9GiRZ6cu6SkRMXFxZ6cGx2DMfYHxtkfGGd/YJxTn5djbGYbD7WdSyMAAADgSxRhAAAA+BJFGAAAAL5EEQYAAIAvUYQBAADgSxRhAAAA+BJFGAAAAL5EEQYAAIAvUYQBAADgSxRhAAAA+BJFGAAAAL5EEQYAAIAvHbEIm9mjZrbTzJYdZr+Z2e/MrNTM3jezk+IfEwAAAIivtswIPy5p6sfsnyapqOVxk6QHjz0WAAAAkFhpRzrAOfe6mQ3+mEMulfSEc85JetvMuppZX+fctniFBAAAQNs55xR1rZ7l5JwUda2eJblobN/+Y51i+xT7R84dtN/tf/0PXtPpg6/dv18f2hc7dnN1VM2RqELB5Lky94hFuA36S9rcar28ZdtHirCZ3aTYrLEKCgpUUlISh9MfvZqaGs/OjY7BGPsD4+wPjLM/JHKcnXOKOKk5GnuEo07hqNQckZr3L+/f7qRwVIpEpbBzikSlyP5tLnZMxO3fL0Va1qNOrZ6dolF9ZHvUuVbLaimXUlQf7I8evF+tnz8ol/u/9pDLCfm3GB/ZaSXqkZVaRdgOse2QY+CcmyFphiRNnDjRFRcXx+H0R6+kpERenRsdgzH2B8bZHxhnf3h17lydNOU01TaGVdsYVnVDWLWNEdU0hlXTsq2mMayG5ojqmyKqb449GpujseWWbQ3NrZ6bImoIR9UUjsY9bzBgSguYQsGA0oKx5di2wIF9wf3bgqZgIKD0gClopkBASgsEFAiYAqaWbbF9wcD+ZR3YFmjZH7DYeQNmMlNs+/7nwAfLZibTB/vNJLPWXyeZrNX2WJkz+2BbbD123IHXUOzA/ccGWhqgtXo9teyT9KHXMJOWL1+ui849S1npwbiPR3vFowiXSxrQar1Q0tY4vC4AAOhEwpGo9tU3q7K+WZV1zdpX36TKuthyZX2z9tU1HdhXWd+s6vpmVbeU3LqmiDTn30c8h5mUFQoqKxRUZiiorPT9ywHlZaapd15Gq21BZaQFlJEWUHrLIyMtGFsO7l9vvS+g9GBQobRYoU1vKbmhYEChoCmt5TnUUmJxdLL3rE6qEizFpwjPknSbmc2UNFnSPq4PBgAgNUSjThV1TdpV3XjgsXP/ck2jdlY1aFdNbL26IXzY1zGT8jND6podUteskLpkp2tAtyzlZaYpJz1Nu7aV6/hRw5WTkaacjDTlZgSVmxFSTkZQuRlpym3ZnpEWODDjCByrIxZhM3tKUrGknmZWLukHkkKS5Jx7SNJsSRdJKpVUJ2l6osICAID4qmsKa3NFvTZX1Kl8b502740tb91Xr13Vjdpd06RI9KNXPOakB9UrL0O98jI0qk+ezhjeU91zMmJFNzukLlkhdc1OV9es2HpeZkjBj5lFLSnZqeIzhibyjwp8RFvuGnHVEfY7SV+JWyIAABA30ahT+d56bdhTq81762Kld2+dyvfWq7yiTntqmz50fGYooAHdstWva5bG9M2Pld3cDPXOzzyw3CsvQzkZ8filMuAt/isGACAFOOe0o6pRq3dUa832aq3Zsf9Ro/rmyIHjQkFT/65ZKuyWrQuOK1Bht2wN6J6twm5ZGtAtWz1z07n0AL5BEQYAoJPZW9sUK7w7qrV6+wfPVa2u0e2Vl6GRBXm6atJAjSjI1ZCeORrQPVsF+Zkfe4kC4CcUYQAAklg06lS6q0aLNuzVoo0VemfjXm3YU3dgf15mmkb1ydMl4/tpZJ88jSiIPbrnpHuYGugcKMIAACSR+qaIlpRXavHGvVq0oULvbKrUvvpmSVL3nHRNGNRNnzt5oMb0y9fIgjwV5GdwKQPQThRhAAA8tKemUfPXV2jRhr1avLFCy7dWKdxyl4bhvXM1bWwfTRjUTRMGddOQnjmUXiCOKMIAAHSwdbtq9NKKHXp5xQ4t3rRXzkkZaQGNH9BVN505VBMGddNJA7upG5c3AAlFEQYAIMEiUad3Nu09UH7LdtdKko7rl687zi3SmSN6aWy/LkpPC3icFPAXijAAAAlQ2xjWvLW79dKKHZq7eqcqapsUCpqmDO2h608brPNGF6hf1yyvYwK+RhEGACBOqhua9fz72/TSih16o3S3msJR5Wem6exRvXX+mAKdOaKX8jNDXscE0IIiDADAMVpavk9Pzt+oWUu2qq4posJuWbpm8kCdP7pAJw/prlCQSx6AZEQRBgCgHeqawpr13lb9ZcEmvV++T5mhgD45vp+unjxI4wu7cHcHoBOgCAMAcBRWba/SX+Zv0rPvbFF1Y1gjCnL1o08ep0+d2F9dsrjsAehMKMIAABxBQ3NEs5du05PzN2nxxr1KTwvoE8f31TWTB2rCoG7M/gKdFEUYAIDD2FxRp8ff2qCnF5drX32zhvbM0Xc/MVqfPamQe/wCKYAiDADAQXZVN+q+V9fqLws2SZIuOK6Prpk8UKcM7cHsL5BCKMIAALSoamjW/75epkfeWK/GcFSfO3mAvnpOkfp0yfQ6GoAEoAgDAHyvoTmiP/1no+4vKVVlXbMuGd9P/3X+CA3pmeN1NAAJRBEGAPhWOBLV04vL9dtX1mrbvgadNaKXvnHhSI3t38XraAA6AEUYAOA7zjm9sGy7fvnv1SrbVasTB3bVr684QacM6+F1NAAdiCIMAPCVN9bu1j1zVun98n0q6p2rGddO0PljCngTHOBDFGEAgC9srqjTXc8u1by1u9W/a5Z+efl4ffrE/goGKMCAX1GEAQApzTmnvy3arB8/t0Jmpu9dPEafnzJQGWlBr6MB8BhFGACQsnZWN+g7/1iqV1bt1JSh3fXLy8ersFu217EAJAmKMAAgJb2wdJvuenapapsi+t7FYzT91MEKcBkEgFYowgCAlLKvvlk/nLVcz767Rcf376JfXzFeRQV5XscCkIQowgCAlDFv7S598+n3tbO6UXecW6TbzhmuUDDgdSwASYoiDADo9OqbIvr5Cyv1x/9s1LBeOXrm1lM1fkBXr2MBSHIUYQBAp/bupr26829LVLa7VjecNkTfnDpSmSHuCAHgyCjCAIBOqTkS1e9eWav755aqb5cs/eXGyTp1eE+vYwHoRCjCAIBOp6qhWbf+ebHeLN2jyyYU6vuXjFF+ZsjrWAA6GYowAKBT2VpZr+mPLdS6XTW697JxunziAK8jAeikKMIAgE5jxdYq3fD4QtU0hvX49Ek6vYhLIQC0H0UYANApzFu7S7f++R3lZqTp77ecotF9872OBKCTowgDAJLe3xdt1neeWarhvXP12PST1bdLlteRAKQAijAAIGk55/TP0ib9s/R9nT68px74/Em8KQ5A3FCEAQBJqTkS1V3PLNU/S5v12ZMKdfdnjld6Gp8SByB+KMIAgKRT3dCsLz/5juat3a1Lh4X0y8vHycy8jgUgxVCEAQBJZfu+Bk1/fKHW7KjWPZ8dp9616yjBABKC3zEBAJLGqu1V+vQDb2rTnlo9ev3JuuJk7hEMIHGYEQYAJIW3Snfr5j8tVlZ6UH+75RQd16+L15EApDiKMADAcws3VGj64ws1qEe2Hps+Sf27cns0AIlHEQYAeGrV9tinxfXvmqWnvjRFPXIzvI4EwCe4RhgA4JnNFXW67pEFyk4P6okvTqIEA+hQFGEAgCd21zTqukcXqKE5oidumKzCbtleRwLgM1waAQDocDWNYU1/bKG2VtbryRsna2SfPK8jAfAhijAAoEM1hiO6+U+LtGJblWZcO0ETB3f3OhIAn+LSCABAh4lEnf7rb0v0Zuke/eKz43Tu6AKvIwHwMYowAKBDOOf0o+eW6//e36bvTBulyyYUeh0JgM9RhAEAHeJ3r5Tqif9s1E1nDtXNZw3zOg4AUIQBAIn357c36n9eXqPPnNRf3546yus4ACCJIgwASLDZS7fpe/9apnNG9dYvPjtOgYB5HQkAJFGEAQAJ9Fbpbn1t5ns6aWA33X/1SQoF+bEDIHnwHQkAkBDLtuzTTX9arME9s/XIFyYqKz3odSQA+BCKMAAg7jZX1On6xxaoS1ZIT9wwWV2z072OBAAfQREGAMRVUziq2556V43NUf3xhknq0yXT60gAcEh8shwAIK7ueXGVlmyu1IPXnKThvXO9jgMAh8WMMAAgbl5esUMPv7Fe150ySNOO7+t1HAD4WBRhAEBcbKms151/X6Lj+uXrrotGex0HAI6IIgwAOGbNkahu/8s7ikSd7r/6JGWGuEMEgOTHNcIAgGP2yzmr9c6mSv3+qhM1uGeO13EAoE3aNCNsZlPNbLWZlZrZtw+xv4uZPWdmS8xsuZlNj39UAEAyenXVDv3h9TJdM3mgLhnfz+s4ANBmRyzCZhaUdL+kaZLGSLrKzMYcdNhXJK1wzo2XVCzpV2bGTSMBIMVt21evO/+2RKP75ut7Fx/8owEAkltbZoQnSSp1zpU555okzZR06UHHOEl5ZmaSciVVSArHNSkAIKmEI1F99al31RSO6v6rT+S6YACdTluKcH9Jm1utl7dsa+0+SaMlbZW0VNIdzrloXBICAJLSr19ao4Ub9upnnzleQ3txv2AAnY855z7+ALPLJV3onLuxZf1aSZOcc7e3OuYySadJ+i9JwyS9JGm8c67qoNe6SdJNklRQUDBh5syZcfyjtF1NTY1yc/mmncoYY39gnL2zdFdYv1rcqDML03TD2IyEnotx9gfGOfV5OcZnn332YufcxIO3t+WuEeWSBrRaL1Rs5re16ZJ+7mKtutTM1ksaJWlB64OcczMkzZCkiRMnuuLi4jb/AeKppKREXp0bHYMx9gfG2Rs7qhp052/naWRBnv5w02nKSk/sJRGMsz8wzqkvGce4LZdGLJRUZGZDWt4Ad6WkWQcds0nSuZJkZgWSRkoqi2dQAID39l8XXNcU0f3XnJjwEgwAiXTEGWHnXNjMbpM0R1JQ0qPOueVmdkvL/ock/UTS42a2VJJJ+pZzbncCcwMAPPDbV9Zq/voK/ery8RreO8/rOABwTNr0gRrOudmSZh+07aFWy1slXRDfaACAZPLG2t26b26pLptQqM9OKPQ6DgAcMz5iGQBwRDurGvS1v76r4b1y9eNLj/M6DgDEBR+xDAD4WM45fePp91XTGNZfvjRF2en86ACQGpgRBgB8rFlLtuq1Nbv0ramjNKKA64IBpA6KMADgsCrrmvTj51Zo/ICuuu6UwV7HAYC44vdbAIDDunv2KlXWN+tPnz5ewYB5HQcA4ooZYQDAIc0v26O/LtqsG08fojH98r2OAwBxRxEGAHxEYzii7zy7VIXdsnTHeUVexwGAhODSCADARzxYsk5lu2r1+PSTuUsEgJTFjDAA4ENKd9bogbnrdMn4fioe2dvrOACQMBRhAMAB0ajTXc8uVWYooO9fPMbrOACQUBRhAMABf1+8WQvWV+iui0arV16G13EAIKEowgAASdLumkb9bPYqTRrcXVdMHOB1HABIOIowAECS9JPnV6iuKayffWasAtwzGIAPUIQBAHptzS79672turV4uIb35mOUAfgDRRgAfK6+KaLv/nOphvbM0ZeLh3kdBwA6DDeHBACf++0ra7W5ol5PfWmKMkNBr+MAQIdhRhgAfGzltir977wyXTGxUKcM6+F1HADoUBRhAPCpSNTpO88sVdeskO66aLTXcQCgw1GEAcCnnpy/Ue9trtT3Lh6jrtnpXscBgA5HEQYAH9q+r0H3vLhaZxT11KUn9PM6DgB4giIMAD70g1nL1ByJ6v99aqzMuGcwAH+iCAOAz8xbu0tzlu/QV88t0qAeOV7HAQDPUIQBwEciUaefzV6lAd2zdOMZQ7yOAwCeoggDgI88++4WrdxWpW9cOEoZadwzGIC/UYQBwCcamiP61b9Xa3xhF118fF+v4wCA5yjCAOATj765Xtv2Neg7F41WIMAb5ACAIgwAPrCnplEPzl2n80b31pShfIIcAEgUYQDwhd+/WqraprC+PW2U11EAIGlQhAEgxW3YXas/v71Rnzt5oIb3zvM6DgAkDYowAKS4e+asUnpaQF8/v8jrKACQVCjCAJDCFm/cq9lLt+umM4eqd16m13EAIKlQhAEgRTnndPfsleqVl6EvnTHU6zgAkHQowgCQouYs36FFG/fq6+eNUE5GmtdxACDpUIQBIAU1R6L6xYurNLx3rq6YWOh1HABIShRhAEhBMxds0vrdtfrOtFFKC/KtHgAOhe+OAJBiqhua9ZuX12rykO46Z1Rvr+MAQNLiojEASDF/eK1Me2qb9NgnRsuMj1IGgMNhRhgAUsj2fQ16+I0yfXJ8P40r7Op1HABIahRhAEghv35ptaJR6RsXjvQ6CgAkPYowAKSIVdur9PfF5brulEEa0D3b6zgAkPQowgCQIu6evUp5GWm67ZzhXkcBgE6BIgwAKeCNtbv12ppduv2cInXNTvc6DgB0ChRhAOjkolGnu19Yqf5ds3TtKYO8jgMAnQZFGAA6uX8t2aLlW6v0zakjlRkKeh0HADoNijAAdGLhSFS/fXmtxvTN1yXj+nkdBwA6FYowAHRi/3pvqzbsqdPXzitSIMCHZwDA0aAIA0AnFY5E9ftXY7PB548p8DoOAHQ6FGEA6KRazwbzUcoAcPQowgDQCTEbDADHjiIMAJ3Q/tngO5gNBoB2owgDQCfTejb4AmaDAaDdKMIA0MnMWsJsMADEA0UYADqR2GxwKbPBABAHFGEA6ERmLdmq9btrmQ0GgDigCANAJ8FsMADEF0UYADqJ/bPBXz2X2WAAiAeKMAB0Avtng0czGwwAcUMRBoBO4MC1wecWKRBgNhgA4oEiDABJLhyJ6j5mgwEg7ijCAJDknnt/q8qYDQaAuGtTETazqWa22sxKzezbhzmm2MzeM7PlZvZafGMCgD+FI1H9/pVSjeqTx2wwAMRZ2pEOMLOgpPslnS+pXNJCM5vlnFvR6piukh6QNNU5t8nMeicoLwD4yv7Z4Ic+fxKzwQAQZ22ZEZ4kqdQ5V+aca5I0U9KlBx1ztaRnnHObJMk5tzO+MQHAfz48G9zH6zgAkHLaUoT7S9rcar28ZVtrIyR1M7MSM1tsZtfFKyAA+NX+2eCvnce1wQCQCEe8NELSob77ukO8zgRJ50rKkvQfM3vbObfmQy9kdpOkmySpoKBAJSUlRx04Hmpqajw7NzoGY+wPqTzOkajTL96o14C8gNJ3rVJJyWqvI3kmlccZH2CcU18yjnFbinC5pAGt1gslbT3EMbudc7WSas3sdUnjJX2oCDvnZkiaIUkTJ050xcXF7Yx9bEpKSuTVudExGGN/SOVxfvbdcm2vW6KHPn+izhnb1+s4nkrlccYHGOfUl4xj3JZLIxZKKjKzIWaWLulKSbMOOuZfks4wszQzy5Y0WdLK+EYFAH+IRB3XBgNABzjijLBzLmxmt0maIyko6VHn3HIzu6Vl/0POuZVm9qKk9yVFJT3snFuWyOAAkKqeWxK7NvjBa7hTBAAkUlsujZBzbrak2Qdte+ig9Xsl3Ru/aADgP5Go0+9eXatRffJ04XHMBgNAIvHJcgCQRF5ctl1lu2p12znDmQ0GgASjCANAknDO6YGSUg3tmaNpPn+DHAB0BIowACSJ19bs0vKtVbrlrGEKMhsMAAlHEQaAJPFAyTr17ZKpT5148GcWAQASgSIMAElg0YYKLVhfoRvPGKr0NL41A0BH4LstACSBB0rWqVt2SFdNGnDkgwEAcUERBgCPrdxWpVdX7dT004YoO71Nd7UEAMQBRRgAPPZgyTrlpAf1hVMGex0FAHyFIgwAHtqwu1bPv79Vn58ySF2yQ17HAQBfoQgDgIf+8HqZ0oIBffH0IV5HAQDfoQgDgEd2VDXoH4vLddmEQvXOz/Q6DgD4DkUYADzy8LwyhaNR3XLmMK+jAIAvUYQBwAOVdU16cv4mXTK+nwb2yPY6DgD4EkUYADzwx7c2qq4poluLmQ0GAK9QhAGgg9U2hvXYW+t13ujeGtUn3+s4AOBbFGEA6GBPLdikyrpm3Vo83OsoAOBrFGEA6ECN4YgenrdeU4Z214RB3byOAwC+RhEGgA707DtbtL2qQV9mNhgAPEcRBoAOEok6/eH1Mo3tn68zinp6HQcAfI8iDAAd5IVl27R+d62+UjxcZuZ1HADwPYowAHQA55zun7tOQ3vl6MLj+ngdBwAgijAAdIiSNbu0cluVbjlrmAIBZoMBIBlQhAGgAzw4d536dcnUp07o73UUAEALijAAJNjCDRVasKFCXzpzqNLT+LYLAMmC78gAkGAPzC1V95x0XXnyQK+jAABaoQgDQAKt2Fqluat36YbTBisrPeh1HABAKxRhAEigB19bp9yMNF07ZbDXUQAAB6EIA0CCbNhdq/97f6uumTJQXbJDXscBAByEIgwACTJjXpnSggF98bQhXkcBABwCRRgAEmBnVYOeXlSuyyYUqnd+ptdxAACHQBEGgAR45M31CkejuvnMoV5HAQAcBkUYAOJsX32znnx7kz4xrp8G9cjxOg4A4DAowgAQZ3/6zwbVNIZ161nDvI4CAPgYFGEAiKP6pogee3ODikf20ph++V7HAQB8DIowAMTR3xZt1p7aJn25eLjXUQAAR0ARBoA4aY5ENeP1Mk0Y1E0nD+7mdRwAwBFQhAEgTp5bslVbKuv15eJhMjOv4wAAjoAiDABxEI06PViyTiML8nTOqN5exwEAtAFFGADi4JVVO7V2Z41uZTYYADoNijAAHCPnnB4oKVVhtyxdPK6v13EAAG1EEQaAYzR/fYXe3VSpm88cqrQg31YBoLPgOzYAHKMHStapZ266Lp84wOsoAICjQBEGgGOwbMs+vb5ml244fYgyQ0Gv4wAAjgJFGACOwYOvrVNeRpo+P2WQ11EAAEeJIgwA7bR+d61eWLpNnz9lkPIzQ17HAQAcJYowALTTjNfXKS0Y0PTTBnsdBQDQDhRhAGiHHVUN+sfiLbpiYqF652V6HQcA0A4UYQBoh4fnlSninG4+c5jXUQAA7UQRBoCjVFnXpL/M36SLx/XVgO7ZXscBALQTRRgAjtIT/9mo2qaIbi1mNhgAOjOKMAAchbqmsB57c73OGdVbo/rkex0HAHAMKMIAcBT+unCz9tY168vMBgNAp0cRBoA2agxHNOP1Mp08uJsmDu7udRwAwDGiCANAG/1j8RZt29eg288p8joKACAOKMIA0AbNkageKCnV+AFddUZRT6/jAADigCIMAG3w7LtbVL63XnecO1xm5nUcAEAcUIQB4AjCkagemFuqsf3zdfbI3l7HAQDECUUYAI7gufe3asOeOt1+ThGzwQCQQijCAPAxIlGn+14t1ag+eTp/dIHXcQAAcUQRBoCPMXvpNq3bVavbzylSIMBsMACkEoowABxGtGU2eHjvXE0b28frOACAOGtTETazqWa22sxKzezbH3PcyWYWMbPL4hcRALzx7xXbtXpHtW4/ZzizwQCQgo5YhM0sKOl+SdMkjZF0lZmNOcxxv5A0J94hAaCjOef0u1dKNaRnji4e18/rOACABGjLjPAkSaXOuTLnXJOkmZIuPcRxt0v6h6SdccwHAJ54ZeVOrdhWpa+cPVxBZoMBICWlteGY/pI2t1ovlzS59QFm1l/SpyWdI+nkw72Qmd0k6SZJKigoUElJyVHGjY+amhrPzo2OwRj7Q6LG2Tmnn/6nQb2yTF33rVVJSWncz4G24++zPzDOqS8Zx7gtRfhQUyHuoPXfSPqWcy7ycffYdM7NkDRDkiZOnOiKi4vbljLOSkpK5NW50TEYY39I1DiXrN6p9VUL9fPPHK/zJg2M++vj6PD32R8Y59SXjGPcliJcLmlAq/VCSVsPOmaipJktJbinpIvMLOyc+2c8QgJAR4ldG7xW/btm6TMnFXodBwCQQG0pwgslFZnZEElbJF0p6erWBzjnhuxfNrPHJT1PCQbQGb21bo/e2VSpn3xqrNLTuMMkAKSyIxZh51zYzG5T7G4QQUmPOueWm9ktLfsfSnBGAOgwv31lrQryM3T5BGaDASDVtWVGWM652ZJmH7TtkAXYOXf9sccCgI73dtkeLVhfoR9cMkaZoaDXcQAACcbv/QCgxe9fXaueuRm6ijfIAYAvUIQBQNLijRV6s3SPbj5zKLPBAOATFGEAkPS7V0rVPSdd10xhNhgA/IIiDMD3lmyu1GtrdunGM4YoO71Nb50AAKQAijAA3/v9q2vVJSuk604Z7HUUAEAHoggD8LVlW/bp5ZU79cXThyg3g9lgAPATijAAX7vv1VLlZaTpC6cO9joKAKCDUYQB+NbKbVV6cfl2TT9tsLpkhbyOAwDoYBRhAL71yzmrlZeZphtOH3LkgwEAKYciDMCXFm6o0CurdurW4mHqmp3udRwAgAcowgB8xzmnX7ywSr3zMjT9VGaDAcCvKMIAfOeVlTu1aONe3XFekbLS+RQ5APArijAAX4lEne6ds1pDeuboiokDvI4DAPAQRRiAr/zz3S1avaNad14wQqEg3wIBwM/4KQDANxrDEf36pTU6vn8XXTS2r9dxAAAeowgD8I0n396kLZX1+ubUkQoEzOs4AACPUYQB+EJ1Q7Pum1uq04b30BlFvbyOAwBIAhRhAL7w8Lz1qqht0jcvHOV1FABAkqAIA0h5u2sa9fC8Ml10fB+NH9DV6zgAgCRBEQaQ8u57tVQN4ajuvGCk11EAAEmEIgwgpW2uqNOT8zfqiomFGtYr1+s4AIAkQhEGkNJ+/dIaBcx0x7kjvI4CAEgyFGEAKWvltir9870tuv60werTJdPrOACAJEMRBpCy7p2zWnkZafryWcO9jgIASEIUYQApacH6Cr26aqduKR6mLtkhr+MAAJIQRRhAynHO6RcvrlJBfoamnzrE6zgAgCRFEQaQcl5euVOLN+7VHeeOUFZ60Os4AIAkRREGkFIiUad756zSkJ45unxioddxAABJjCIMIKU8++4WrdlRo//vgpEKBfkWBwA4PH5KAEgZjeGI/uelNTq+fxdNG9vH6zgAgCRHEQaQMv789iZtqazXt6aOUiBgXscBACQ5ijCAlLCrulG/eXmNzijqqdOLenodBwDQCVCEAaSEu2evVENzRD/85HFeRwEAdBIUYQCd3vyyPXrm3S266cyhGtYr1+s4AIBOgiIMoFNrjkT1/X8tV/+uWbrt7CKv4wAAOhGKMIBO7Y9vbdDqHdX6/iVj+PAMAMBRSfM6AAC0196GqP7nrTU6e2QvXTCmwOs4AIBOhhlhAJ3WU6ua1Bx1+uEnj5MZt0sDABwdijCATunN0t1asD2iLxcP06AeOV7HAQB0QhRhAJ1OYzii7/1rmXpnm245a5jXcQAAnRRFGECn88gb61W2q1afH52uzBBvkAMAtA9FGECnsqWyXr9/pVQXHlegcb14vy8AoP0owgA6lR8/t1xOTt+/hE+QAwAcG4owgE5j7uqdmrN8h24/p0j9u2Z5HQcA0MlRhAF0Cg3NEf1w1nIN7ZWjL50x1Os4AIAUwAV2ADqFh15bp4176vTkjZOVnsb/wwMAjh0/TQAkvY17avVAyTpdPK6vThve0+s4AIAUQREGkNScc/rhrOUKBUzf/cQYr+MAAFIIRRhAUntpxQ7NXb1LXz9/hPp0yfQ6DgAghVCEASSt+qaIfvTcCo0syNMXTh3sdRwAQIrhzXIAktZ9c9dqS2W9/nrTFIWC/H87ACC++MkCICmt2FqlGa+X6TMn9tfkoT28jgMASEEUYQBJp6E5ojtmvquu2en67sW8QQ4AkBhcGgEg6dw9e6XW7qzREzdMUvecdK/jAABSFDPCAJLK3FU79cf/bNQXTx+iM0f08joOACCFUYQBJI1d1Y36xtNLNKpPnr5x4Uiv4wAAUhyXRgBICs45ffPpJapuCOsvX5qizFDQ60gAgBTHjDCApPCntzdq7upduuui0RpRkOd1HACAD1CEAXhuzY5q/fT/Vurskb103SmDvI4DAPAJijAATzWGI/rqU+8qNyNN91w2XmbmdSQAgE+0qQib2VQzW21mpWb27UPsv8bM3m95vGVm4+MfFUAquvfF1Vq1vVr3Xj5OvfIyvI4DAPCRIxZhMwtKul/SNEljJF1lZgff4X69pLOcc+Mk/UTSjHgHBZB6Xl+zSw+/sV7XnTJI54wq8DoOAMBn2jIjPElSqXOuzDnXJGmmpEtbH+Cce8s5t7dl9W1JhfGNCSDVVNQ26c6/L9Hw3rm666LRXscBAPhQW4pwf0mbW62Xt2w7nC9KeuFYQgFIbc45fesf72tfXbN+e+UJ3CoNAOAJc859/AFml0u60Dl3Y8v6tZImOeduP8SxZ0t6QNLpzrk9h9h/k6SbJKmgoGDCzJkzj/1P0A41NTXKzc315NzoGIxxcivZ3KzHlzfpypHpmjok1O7XYZz9gXH2B8Y59Xk5xmefffZi59zEg7e35QM1yiUNaLVeKGnrwQeZ2ThJD0uadqgSLEnOuRlquX544sSJrri4uA2nj7+SkhJ5dW50DMY4eZXurNHMV+bp9OE99bMvTFIg0P67RDDO/sA4+wPjnPqScYzbcmnEQklFZjbEzNIlXSlpVusDzGygpGckXeucWxP/mABSQVM4qq/99V1lhoL61RXjj6kEAwBwrI44I+ycC5vZbZLmSApKetQ5t9zMbmnZ/5Ck70vqIemBlnuAhg81/QzA33790hot21Klhz4/QQX5mV7HAQD4XFsujZBzbrak2Qdte6jV8o2SboxvNACp5M3S3frD6+t01aQBmjq2j9dxAADgk+UAJF7pzhrd+ufFGtYrV9+7+ODbkAMA4A2KMICE2lPTqBseX6j0tIAeu/5kZae36RdRAAAkHD+RACRMQ3NENz6xSDuqGjTzpika0D3b60gAABxAEQaQENGo03/97T29t7lSD1x9kk4c2M3rSAAAfAiXRgBIiHvmrNbspdv1nWmjNO34vl7HAQDgIyjCAOLuqQWb9NBr63TN5IH60hlDvY4DAMAhUYQBxNVra3bpu/9cpuKRvfSjTx6nlnuLAwCQdCjCAOJm1fYqfeXJd1TUO1f3XX2S0oJ8iwEAJC9+SgGIix1VDbrhsYXKyQjqseknKzeD9+ICAJIbP6kAHLPaxrC++MeFqqxv1t9uPkV9u2R5HQkAgCNiRhjAMYlEne6Y+a5WbK3S/VefpLH9u3gdCQCANmFGGMAx+cnzK/Tyyp36yaXH6exRvb2OAwBAmzEjDKDdHntzvR5/a4O+ePoQXXvKYK/jAABwVCjCANrlpRU79OPnV+jC4wp010WjvY4DAMBRowgDOGqvrNyh2/7yjsb176LffO5EBQPcKxgA0PlQhAEclX+9t0U3/2mxRvbJ02PTJykrPeh1JAAA2oU3ywFosz/9Z4O+P2u5Jg/prv+9bqLyMkNeRwIAoN0owgCOyDmnB0rW6d45q3Xe6N667+qTlBliJhgA0LlRhAF8LOec7n5hlWa8XqZPn9hf91w2TiE+OhkAkAIowgAOKxJ1uuuZpfrros36wimD9INLjlOAN8YBAFIERRjAITWGI/r6X9/T7KXb9dVzi/T184pkRgkGAKQOijCAj6hrCuvmPy3WvLW79b2Lx+iLpw/xOhIAAHFHEQbwIfvqmjX98QV6b3Ol7rlsnK6YOMDrSAAAJARFGMABO6sbdN0jC1S2q1YPXHOSpo7t63UkAAAShiIMQJK0uaJO1z4yXzurG/Xo9Sfr9KKeXkcCACChKMIAtGzLPt34x0Wqb47ozzdO1kkDu3kdCQCAhONmoICPOef06Bvr9ekH3pSZ9Nebp1CCAQC+wYww4FN7a5v0jaeX6OWVO3Xe6ALde9k4dctJ9zoWAAAdhiIM+ND8sj26Y+Z7qqht0g8uGaPrTx3MPYIBAL5DEQZ8JBJ1uu/VUv32lTUa1CNHz3zhVI3t38XrWAAAeIIiDPjE9n0NumPmu5q/vkKfObG/fvypscrN4FsAAMC/+CkI+MCrq3bozr8tUWM4ql9dPl6fnVDodSQAADxHEQZSWGM4onteXK1H3livMX3z9furT9SwXrlexwIAIClQhIEUtWF3rW5/6l0t3bJP1586WN+eNkqZoaDXsQAASBoUYSDFRKNOT79Trh/NWq60YEAzrp2gC47r43UsAACSDkUYSCGLN1box8+t0JLyfZo0uLt+c+UJ6tc1y+tYAAAkJYowkALK99bp5y+s0vPvb1Of/Ez9z+fG69Lx/RUIcG9gAAAOhyIMdGI1jWE9WFKq/523XgGT7ji3SDefNVTZ6fzVBgDgSPhpCXRC+68DvnfOau2qbtSnTuinb04dxWUQAAAcBYow0Mm8XbZHP3l+hZZvrdKJA7tqxrUTdOLAbl7HAgCg06EIA53Epj11uvuFlXph2Xb165Kp3111oi4Z11dmXAcMAEB7UISBJLezukGPzFuvx97coLSg6c7zR+hLZw7lnsAAABwjijCQpJZt2adH31yv55dsU1MkqssmFOobF45UQX6m19EAAEgJFGEgiUSiTi+t2K5H39igBRsqlJ0e1FWTBugLpw7WUD4aGQCAuKIIA0lgX32z/rZwsx5/a4O2VNarsFuWvvuJ0bp84gB1yQp5HQ8AgJREEQY8VLarRo+/tUFPLy5XXVNEk4Z01/cuHqPzxxQoyIdhAACQUBRhoIM55zRv7W499uZ6zV29S+nBgC4Z30/TTxussf27eB0PAADfoAgDHcA5p/c2V+qFZds1e+k2le+tV8/cdH3tvCJdM3mQeuVleB0RAADfoQgDCRKNOr27uVKzl27TC0u3aeu+BoWCptOH99R/nT9CnxjXVxlp3AINAACvUISBOIpGnRZv2qvZS7fpxWXbtW1fg9KDAZ1R1FN3XjBS540uUJds3vwGAEAyoAgDxygSdVq0oSI287tsu3ZWNyo9LaCzRvTSN6eO1LmjC5SfSfkFACDZUISBo+Sc09qdNZpftkdvl1Xo7bI92lPbpIy0gIpH9tJFx/fVOaN6K4/yCwBAUqMIA0cQjcaK79tle/R22R4tWF+hPbVNkqR+XTJ15oheOmdUb50zqrdyMvgrBQBAZ8FPbeAg0ajT6h3Vertsj+aXVWj++j3aW9csSerfNUtnjeylKUN76JShPVTYLUtm3O8XAIDOiCIMXwtHolq/u1YrtlXFHlurtHTLPlW2FN/Cblk6d3SBJg/prilDe2hA92yPEwMAgHihCMM3qhuatWp7tVZsrdLKluK7enu1GsNRSVJ6MKCiglxdOKaPJg3prslDu6uwG8UXAIBURRFGSnHOaXdNk9bujWjP4nJtrKjT6u1VWrmtWpsq6g4c1y07pDH98nXdKYM0pl++RvfN17BeuQoFAx6mBwAAHYkijE4nEnXatq9em/bUacOeOm2sqNXG3XXaWFGnTXtqVdsUiR04f4nMpCE9cnR8/y763MkDNKZvrPQW5GdwbS8AAD5HEUZSaQpHtbO6QTuqGrWzqkE7qhq0ff9ydYO27WtQeUW9miLRA1+THgyosHuWBvfI0eQh3TW4R7b2bV2nS4qnqLBbttLTmOUFAAAfRRFGQjnnVN0YVmVtsyrqmrS3rkl7a5u0t65Ze2ubtLum8UNld/9tyVoLBU298zJVkJ+hUX3ydP6YAg3ukaNB3bM1qGeO+uRnKhj48OxuSclGDe2V21F/TAAA0Am1qQib2VRJv5UUlPSwc+7nB+23lv0XSaqTdL1z7p04Z4UHnHNqaI6qurFZ1Q1h1TSEY88t67HlsKobYuuVdbHCW1nXpIraZlXWNSkcdYd87YBJ3XMy1KdLhvp1ydQJA7qqT36s8BbkZ7Y8MtQtO12BAJcxAACA+DpiETazoKT7JZ0vqVzSQjOb5Zxb0eqwaZKKWh6TJT3Y8owEiUSdmsJRNYYjagxH1djcavkw2xuaI6prij1iy+FWy7FHfVNE9c2x59qmWPE9XJFtLSsUVG5mmrplh9Q1O11De+ZqwqCQumWnxx456eqWHWp5ji3nZ4YouAAAwDNtmRGeJKnUOVcmSWY2U9KlkloX4UslPeGcc5LeNrOuZtbXObct7omPQUNzRO9trtTKPRGlrd2tqHOKOCfnnKJRfbDsYkUz6pxcy3LEudhzy/YPL0uRaDT27JyiUadw1CkSjao5Ejtu/3o4sn/ZKdxqPRyNqjns1BSJqjkS294ciR5Yb25Z378caUM5/TiZoYCy09OUFQoqKz2o7PSgMkNB9chNV3Z6UFmhNGWnB5WXmabczDTlZYaUl5EWW89oWc+MredkpHG3BQAA0Om0pQj3l7S51Xq5Pjrbe6hj+ktKqiK8p7ZJV854O7aycH7CzhMwKS0YUFrAFAyY0gL2ofVQMHBg+/7nUDCgUDCg/PSQQvvX0wIKBU3pLfti20yhQGw5IxRQRlpAGWnB2HMooPRgQBmhlvWWfektywcKb1qQmVgAAOB7bSnCh2pMB09HtuUYmdlNkm6SpIKCApWUlLTh9PHTFHH61smZamioV3ZWlgImmcWKa0D7l022f1mt9rc6LmD2oa8Nttof+9q2lEynQ/wrOjpOUnPL46DNDS0Pv6qpqenw/77Q8Rhnf2Cc/YFxTn3JOMZtKcLlkga0Wi+UtLUdx8g5N0PSDEmaOHGiKy4uPpqscVNSUiKvzo2OwRj7A+PsD4yzPzDOqS8Zx7gtF3YulFRkZkPMLF3SlZJmHXTMLEnXWcwUSfuS7fpgAAAAoLUjzgg758JmdpukOYrdPu1R59xyM7ulZf9DkmYrduu0UsVunzY9cZEBAACAY9em+wg752YrVnZbb3uo1bKT9JX4RgMAAAASh3teAQAAwJcowgAAAPAlijAAAAB8iSIMAAAAX6IIAwAAwJcowgAAAPAlijAAAAB8iSIMAAAAX6IIAwAAwJcowgAAAPAlijAAAAB8iSIMAAAAXzLnnDcnNtslaaMnJ5d6Strt0bnRMRhjf2Cc/YFx9gfGOfV5OcaDnHO9Dt7oWRH2kpktcs5N9DoHEocx9gfG2R8YZ39gnFNfMo4xl0YAAADAlyjCAAAA8CW/FuEZXgdAwjHG/sA4+wPj7A+Mc+pLujH25TXCAAAAgF9nhAEAAOBzvizCZvYTM3vfzN4zs3+bWT+vMyH+zOxeM1vVMtbPmllXrzMh/szscjNbbmZRM0uqdyPj2JjZVDNbbWalZvZtr/MgMczsUTPbaWbLvM6CxDCzAWY218xWtny/vsPrTPv5sghLutc5N845d4Kk5yV93+M8SIyXJI11zo2TtEbSdzzOg8RYJukzkl73Ogjix8yCku6XNE3SGElXmdkYb1MhQR6XNNXrEEiosKQ7nXOjJU2R9JVk+fvsyyLsnKtqtZojiQulU5Bz7t/OuXDL6tuSCr3Mg8Rwzq10zq32OgfibpKkUudcmXOuSdJMSZd6nAkJ4Jx7XVKF1zmQOM65bc65d1qWqyWtlNTf21QxaV4H8IqZ/VTSdZL2STrb4zhIvBsk/dXrEADarL+kza3WyyVN9igLgDgxs8GSTpQ03+MoklK4CJvZy5L6HGLXfzvn/uWc+29J/21m35F0m6QfdGhAxMWRxrnlmP9W7NcyT3ZkNsRPW8YZKccOsY3f3gGdmJnlSvqHpK8d9Nt5z6RsEXbOndfGQ/8i6f9EEe6UjjTOZvYFSRdLOtdxr8BO6yj+PiN1lEsa0Gq9UNJWj7IAOEZmFlKsBD/pnHvG6zz7+fIaYTMrarX6SUmrvMqCxDGzqZK+JemTzrk6r/MAOCoLJRWZ2RAzS5d0paRZHmcC0A5mZpIekbTSOfdrr/O05ssP1DCzf0gaKSkqaaOkW5xzW7xNhXgzs1JJGZL2tGx62zl3i4eRkABm9mlJv5fUS1KlpPeccxd6GgpxYWYXSfqNpKCkR51zP/U2ERLBzJ6SVCypp6Qdkn7gnHvE01CIKzM7XdI8SUsV616SdJdzbrZ3qWJ8WYQBAAAAX14aAQAAAFCEAQAA4EsUYQAAAPgSRRgAAAC+RBEGAACAL1GEAQAA4EsUYQAAAPgSRRgAAAC+9P8DBYdYXe2lIaMAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"beta1, beta2 = 2, 3\n",
"X = np.linspace(-3, 2)\n",
"P = 1 / (1 + np.exp(-beta1 - beta2 * X))\n",
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.plot(X, P)\n",
"ax.set_title(r\"$P = 1/(1+e^{-Z})$\")\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"id": "465224ba-c927-4a27-a32a-cc52b25a455d",
"metadata": {},
"source": [
"We usually denote $Z_i = \\beta_1+\\beta_2X_i$, ranges from $(-\\infty, +\\infty) $. This model solves the issues of probabilistic range and linear model. \n",
"\n",
"In order to estimate the logit model, we are seeking ways to linearize it. It might look confusing at first sight, here is the steps\n",
"$$\n",
"P_{i}=\\frac{1}{1+e^{-Z_{i}}}=\\frac{e^{Z_i}}{e^{Z_i}}\\frac{1}{1+\\frac{1}{e^{Z_i}}}=\\frac{e^{Z_i}}{1+e^{Z}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "970d9bfd-f897-43e5-aa2b-e9bc2737a11e",
"metadata": {},
"source": [
"And $1-P_i$\n",
"$$\n",
"1-P_{i}=\\frac{1}{1+e^{Z_{i}}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "095498cc-58eb-4939-a549-9e92ccecb257",
"metadata": {},
"source": [
"Combine them, we have the odds ratios\n",
"$$\n",
"\\frac{P_{i}}{1-P_{i}}=\\frac{1+e^{Z_{i}}}{1+e^{-Z_{i}}}=e^{Z_{i}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "fd3a616f-2064-40a0-a61b-6dc271833c1b",
"metadata": {},
"source": [
"Here's the interesting part, take natural log, we get a linearized model and we commonly call the log odds ratios the **logit**.\n",
"$$\n",
"\\underbrace{\\ln{\\bigg(\\frac{P_{i}}{1-P_{i}}\\bigg)}}_{\\text{logit}}= Z_i =\\beta_1+\\beta_2X_i\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "4623c027-c2ea-494a-850a-0c1606cab5bf",
"metadata": {},
"source": [
"To estimate the model, some technical procedures have to be carried out, because logit wouldn't make any sense, if we simply plug in the data $Y_i$, because we don't observe $P_i$, the results are nonsensical.\n",
"$$\n",
"\\ln \\left(\\frac{1}{0}\\right) \\quad \\text{if a family own a house}\\\\\n",
"\\ln \\left(\\frac{0}{1}\\right) \\quad \\text{if a family does not own a house}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "ad79dcb7-6d0f-4118-9fcc-039f08f5e976",
"metadata": {},
"source": [
"One way to circumvent this technical issue, the data can be grouped to compute\n",
"$$\n",
"\\hat{P}_{i}=\\frac{n_{i}}{N_{i}}\n",
"$$\n",
"where $N_i$ is number of families in a certain income level, e.g. $[20000, 29999]$, and $n_i$ is the number of family owning a house in the that level."
]
},
{
"cell_type": "markdown",
"id": "120d1233-e26d-4c5e-a736-59267a7dc16c",
"metadata": {},
"source": [
"However, this is not a preferred method, since we have nonlinear tools. Recall that owning a house is following a Bernoulli distribution whose probability mass function is\n",
"$$\n",
"f_i(Y_i)= P_i^{Y_i}(1-P_i)^{1-Y_i}\n",
"$$\n",
"The joint distribution is the product of Bernoulli PMF due to their independence\n",
"$$\n",
"f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right)=\\prod_{i=1}^{n} f_{i}\\left(Y_{i}\\right)=\\prod_{i=1}^{n} P_{i}^{Y_{i}}\\left(1-P_{i}\\right)^{1-Y_{i}}\n",
"$$\n",
"In order to get its log likelihood function, we take natural log\n",
"$$\n",
"\\begin{aligned}\n",
"\\ln f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right) &=\\sum_{i=1}^{n}\\left[Y_{i} \\ln P_{i}+\\left(1-Y_{i}\\right) \\ln \\left(1-P_{i}\\right)\\right] \\\\\n",
"&=\\sum_{i=1}^{n}\\left[Y_{i} \\ln P_{i}-Y_{i} \\ln \\left(1-P_{i}\\right)+\\ln \\left(1-P_{i}\\right)\\right] \\\\\n",
"&=\\sum_{i=1}^{n}\\left[Y_{i} \\ln \\left(\\frac{P_{i}}{1-P_{i}}\\right)\\right]+\\sum_{i=1}^{n} \\ln \\left(1-P_{i}\\right)\n",
"\\end{aligned}\n",
"$$\n",
"Replace \n",
"$$\n",
"\\ln{\\bigg(\\frac{P_{i}}{1-P_{i}}\\bigg)}=\\beta_1+\\beta_2X_i\\\\\n",
"1-P_{i}=\\frac{1}{1+e^{\\beta_1+\\beta_2X_i}}\n",
"$$\n",
"Finally we have log likelihood function of $\\beta$'s\n",
"$$\n",
"\\ln f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right)=\\sum_{i=1}^{n} Y_{i}\\left(\\beta_{1}+\\beta_{2} X_{i}\\right)-\\sum_{i=1}^{n} \\ln \\left[1+e^{\\left(\\beta_{1}+\\beta_{2} X_{i}\\right)}\\right]\n",
"$$\n",
"Take partial derivative w.r.t to $\\beta_2$\n",
"$$\n",
"\\frac{\\partial LF}{\\partial \\beta_2} = n\\beta_2 -\\sum_{i=1}^n\\frac{e^{\\beta_1+\\beta_2X_i}X_i}{1+e^{\\beta_1+\\beta_2X_i}}=0\n",
"$$\n",
"And stop here, apparently there won't be analytical solutions, i.e. numerical solutions are needed. And this is exactly what Python is doing for us below."
]
},
{
"cell_type": "markdown",
"id": "278e2171-c593-4234-8231-36bc99d4a600",
"metadata": {},
"source": [
"## Single Variable Logit Model "
]
},
{
"cell_type": "markdown",
"id": "01a5d91e-b94a-4a3e-9f58-3eeb15d5ff57",
"metadata": {},
"source": [
"Here we come back to the house owning example.\n"
]
},
{
"cell_type": "code",
"execution_count": 149,
"id": "aa14ab17-f25c-498a-b064-f3eb905daaa5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.257124\n",
" Iterations 8\n"
]
}
],
"source": [
"df = pd.read_excel(\n",
" \"Basic_Econometrics_practice_data.xlsx\", sheet_name=\"HouseOwn_Quali_Resp\"\n",
")\n",
"df.head()\n",
"\n",
"Y = df[[\"Own_House\"]]\n",
"X = df[[\"Annual_Family_Income\"]]\n",
"X = sm.add_constant(X)\n",
"\n",
"log_reg = sm.Logit(Y, X).fit()"
]
},
{
"cell_type": "markdown",
"id": "834a1a62-aa88-4ccd-9d09-5361ca5430a8",
"metadata": {},
"source": [
"Optimization of likelihood function is an iterative process and global optimization has been reached as the notice shows. \n",
"\n",
"Print the estimation results."
]
},
{
"cell_type": "code",
"execution_count": 118,
"id": "713ee3c5-0d95-4339-8607-631964b57c5c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Own_House No. Observations: 38\n",
"Model: Logit Df Residuals: 36\n",
"Method: MLE Df Model: 1\n",
"Date: Wed, 08 Sep 2021 Pseudo R-squ.: 0.6261\n",
"Time: 20:42:19 Log-Likelihood: -9.7707\n",
"converged: True LL-Null: -26.129\n",
"Covariance Type: nonrobust LLR p-value: 1.067e-08\n",
"========================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"----------------------------------------------------------------------------------------\n",
"const -6.7617 2.189 -3.088 0.002 -11.053 -2.471\n",
"Annual_Family_Income 0.1678 0.053 3.176 0.001 0.064 0.271\n",
"========================================================================================\n"
]
}
],
"source": [
"print(log_reg.summary())"
]
},
{
"cell_type": "markdown",
"id": "5212ea16-a24d-43d5-95a4-679bfbc291bc",
"metadata": {},
"source": [
"Here is the estimated model\n",
"$$\n",
"\\hat{P_i} = \\frac{1}{1+e^{6.7617-0.1678 X_i}}\n",
"$$\n",
"And the fitted curve"
]
},
{
"cell_type": "code",
"execution_count": 119,
"id": "43d83892-c8d8-4ebc-a1c2-65b702042b25",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAGpCAYAAACkkgEIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAoUlEQVR4nO3dd5xddZ3/8ddnZtITSCAhkFBCiXQ0ELpIQpEiCCIiCKJYWPYnlnXXtuvuuuu6uotlXRuiYlmVqIB0CDV0MAkBqZEQJJUQSnqbzHx/f5wb7mQyM7mTzJ1zy+v5eNzH/X7POffc98x3Mvnk5Hu/J1JKSJIkSSpNQ94BJEmSpGpiAS1JkiR1gwW0JEmS1A0W0JIkSVI3WEBLkiRJ3dCUd4DuGj58eBozZkwu771y5UoGDRqUy3urdzjG9cFxrg+Oc31wnGtfnmM8ffr0V1NKI9pvr7oCesyYMUybNi2X954yZQoTJkzI5b3VOxzj+uA41wfHuT44zrUvzzGOiJc62u4UDkmSJKkbLKAlSZKkbrCAliRJkrrBAlqSJEnqBgtoSZIkqRssoCVJkqRusICWJEmSusECWpIkSeoGC2hJkiSpGyygJUmSpG6wgJYkSZK6wQJakiRJ6oayFdARcWVEvBIRT3WyPyLifyNiVkT8OSIOLlcWSZIkqac0lfHcvwC+D/yqk/2nAGMLj8OBHxWeJdWI62bM57LJM1mwZDWjhg7gcyftzZnjRucda6ttydfVU9+Lzs5T7u2d+fJ1T3LVo3NpSYnGCM47fBf+48wDe+V7sbnzn7vLcv7pG3dX7M/dlnwftuT7Xe5M1fK+PfUzXwsq7WuutDyliJRS+U4eMQa4KaV0QAf7fgxMSSldVejPBCaklBZ2dc7x48enadOmlSPuZk2ZMoUJEybk8t7qHY5xz7luxny+dO2TrG5ueXPbgD6NfP2sA3P/xbg147wlX1dPfS86O897DxnNNdPnl217Zzm/fN2T/PqROZtsv+CIXTst6sr9c9H2/H9/4Hq+9WRTxfzctbUl34ct+X6XO1OHUoKWFmhZ3+bR0mZby0bb73v2ZX527/O0rl9PEy30Sa0MbITzD92ZQ3fZJjuutQVaW7NHSlk/JWhtZfYLs9hjzJg2+9sel7VfeGUZj8x6ldbWVhpppSEl+jTAHtsPYO5rK0mtrTSmVgLo0wAHjNqGUdv0A1J2nkTHbQr9tm0K70vb7e1fV9jX4TnaPKDjY9p+r7sahw76K9a28PrKtbSt/xoCthvUl0F9mzY+f/vzbPJ2W5Cl3XFrmltYvqaZlCA2bAwY0q+J/k3ZRIl1zevou+PO8OX/7fw9yiQipqeUxrffXs4r0JszGpjbpj+vsK3LAlpSdbhs8syN/iIGWN3cwmWTZ1ZUIdNdW/J19dT3orPzbLgiWa7tneW86tG5m2zbsL2zgq7cPxfV8nO3JTm35Pu9kdYWWLsme6xZDWtXF/trV/P0tTM4Z+VKBqVmBqZm+qf19KOFpl/dBVO3h+Z1sL4ZmpuzdvM6WL9u435zoZ9aS/5evKPw2MQdpb1+D4DH7+vymD0Lj03Mg3EdbZ9d2ntXo8GFxyaW9XKQgv6Fx0YSsLrY7QuwcnlvRSpJnlegbwa+nlJ6oNC/C/h8Sml6B8deDFwMMHLkyEMmTZpUtsxdWbFiBYMHd/hjpxrhGPecJ+cv7XTfgaO37cUkm9qacd6Sr6unvhddnafcOsqZ5/eiM23PP3IALGrzl3DeP3dtbdH3bt4S+rSsZ0DzagY0r6Zf81r6rV9Hv/Vr2W1g0LRuDU3r1mbPzWtoWruGpuas37i+mcaW9eX6cqSyWz14KI+eeUmvv+/EiRM7vALtFI5u8L/3a59j3HOO/sbdzF+yepPto4cO4MEvHpdDoqKtGect+bp66nvR2XkaIza5ctyT2zvLueeXbun0PC98/dRufQ099XPR9vwbpnD05Pl7Stuc/VMzO7auZKeWFezbfx3/cuwoWL4Eli3JnpcvheVLWP3a6wygSorgaIA+faChERobobGpg+es/eeXV7JqPTRHA+tpYH3huU/fPpxwwKjia6IBGgqPDe0IXpo3j93GjNloGw2NbdoNfPfuF3h99XpaCBJBSwQtNBARNKeglaA14s39wwb14z/POjB7PVGYWxBZf6N2YdJB2+OiYePtmzvHm6+Ldse3Pweb7iucpt03v02z/c7gnB8/zKLla0jtXjhym/784ZIjOz5p224H5+x8Xyeva9M+4/sPsGDZ2k0O33Gb/tz4ybcDwYMPPcTRb387DOn9fwRX4hSOG4BLI2IS2YcHl26ueJZUPT530t4dzqf83El755hq623J19VT34vOzlPuOdCd5Tzv8F06nJN73uG7dPtr6Kmfi4r7uWtthWVvwOKF8MpCeP0VeONVrlk3j6VLFrJDywqGpTXF45cBf+j4VAO2NksE9O0P/ftDvwHQr3/hMQD6D2DOylYemLOCZamJVdGHVTSRmvpy2iG7MW7PHaBP36wo7tMXmjppb3g0NpYca3YXc68pYdrNi1OmsNtm/kG824j5XN7Jn50bO5n3X8p7V6MPnH5Eh9/vz552IIzYqdfzXPSu8R3m+ei7DoRttwOgecCgXIrnrpStgI6Iq4AJwPCImAf8K9AHIKV0OXALcCowC1gFXFSuLJJ634Z5nNX2yerN2ZKvq6e+F12dZ/xu25V1e0c2zLvtzqoQ5f65aHt+WM7o3vi5SwmWvg4L58CCOfDKAlj8Mry6MHtuXrfJS3YsPLqruaGJxfTntRjA0ob+DB0+jAPGjoYBg2DQYBgwGAYOgkFD2mwbBP0HZoVtQ+er1+4KDJwxnx+0G5txZf4z2xu/K7bkz06tqrTfzZWWp1RlncJRDk7hUDk5xvXBca4PZRnnpW/A3BeyQnlDwbxwDqxaseXnbGyEodvDsOEwdDgM3Q6GDM2uuA0ZCttseB6aXS3u6r/J65B/nmtfnmNciVM4JEmqTCnBG6/CS7Ngzqzi85LXun+uQUNgxI4wYhQMH5kVysNGFJ6HZ4VyF1eGJVUeC2hJktavhzkvwAvPwKyn4YVnu1cs9x8Io3aFnXaFkaNhh52y+aQjdoKBruwj1RoLaElS/Vm/Hl58Dp6ZATP/DH/9C6zbdCWATfTtBzvvnj1G7ZYVzDvtkl1JdmqFVDcsoCVJtS8lWDgXnnkMnp0BM5+ENau6fk3/gbDbXrDrXrDrnll7x52zJdIk1TULaElSTYrWFnj2cXj84ezx2itdv2D4SNhzP9hrf9hrPxi9m8WypA5ZQEuSakfzOnhqOky/n6MeewjWren82O13gP0Ohn3HwVsOyFbCkKQSWEBLkqpba0s2JePRe+CxB99cUq5P++MGDIR93gb7jcsK5x1GOW9Z0haxgJYkVadF8+H+2+Dhu7IbmXRk2HB425Ew7kh4y4HQtElZLUndZgEtSaoezetgxkNw363w3BMdHzN8JBw2gWlpAOPPer9XmSX1OAtoSVLlW/Ia3H1DVjivWLbp/iFD4dB3wOETYI99IYIVU6ZYPEsqCwtoSVLlmjsb7rgWHp0CLes33hcNcNChcMwpcOCh2S2xJakXWEBLkirPX56Em67K1m1ub7sd4JiT4Oh3wnYjej+bpLpnAS1JqhzPPwXX/xqee3zTfXvuB+88K/tAoOszS8qRBbQkKX8vPJMVzu2vOEcDHHI0nHgW7LlvPtkkqR0LaElSfhYvhKt/BtMf2Hh7QwMcdQK86zwYsVM+2SSpExbQkqTet3ol3DwJ7rwO1jcXt0cDHHkcnPaB7EYnklSBLKAlSb2ntRUemAx//CUsX7LxvsMmwLsvgB13ziOZJJXMAlqS1DvmvQi/+i7Mfm7j7XvsA+dekj1LUhWwgJYklde6tXDTb2Hy1dDSUty+3Qh470eyK8/e8ERSFbGAliSVzwvPwpXfhEXzi9sam+CUc7JHv/75ZZOkLWQBLUnqeeub4YbfwK2/h9Ra3D52f/jgp2DUbvllk6StZAEtSepZ8/8KP/3v7DbcGwwYCO/7OLz9pGyJOkmqYhbQkqSekRI8eDv89ofZvOcN9nkrXPRZ2H5kftkkqQdZQEuStt6a1fDr78Ejdxe39embfUjwuHd71VlSTbGAliRtnQVz4If/Di/PK24btSv8zT/C6DG5xZKkcrGAliRtuccfyeY7r1lV3Hb0O+ED/88VNiTVLAtoSVL3pZTdivv6X2VtgL794IJPwlEn5JtNksrMAlqS1D3N6+Bn34Rp9xW3bb8DXPoV2GWP3GJJUm+xgJYklW7FMvj+V2DWM8Vtex8El/wjDBmaVypJ6lUW0JKk0ix+Gb775Y0/LDjxNHj/JdDkXyeS6oe/8SRJm/fX5+G7/wzLlxS3nXMxvPOs3CJJUl4soCVJXXv+KfjuvxRX2mjqAx/7HIx/R765JCknFtCSpM49/Rj84N+KdxYcOBg++RUYe0CusSQpTxbQkqSOzXgYfvyfsL45628zDD77n7Dz7vnmkqScWUBLkjY1/YGseG5tzfrbjYDPfh123DnfXJJUASygJUkbe/wRuOLrxeJ5xE7wD9+A7Ufmm0uSKoQFtCSp6KlpcPnXoKUl648cDZ/7bxi6fb65JKmCNOQdQJJUIZ59HH7w78U5zyN2gn/4L4tnSWrHAlqSBC/OzO4w2Lwu62+3QzZtY9jwXGNJUiWygJakerdofrbO89o1WX/YcOc8S1IXLKAlqZ4tfR2+80+wYmnWHzQE/u4/YYdR+eaSpApmAS1J9Wr1Sviff4ZXX876ffvBp/4dRu2aby5JqnAW0JJUj1pa4Edfg7kvZP2GBrjkH2HPffPNJUlVwAJakurRpMvhmceK/Qs/AwcdnlscSaomFtCSVG/uvgHuubHYP/18ePs788sjSVXGAlqS6snT07Orzxsceiy8+4L88khSFbKAlqR68fJcuPw/i7foHvMWuOizEJFvLkmqMhbQklQP1qyGH3w1W3kDsrWeL/3XbOUNSVK3WEBLUq1LCX7xHVg4J+v36ZsVz96iW5K2iAW0JNW6O/4I0+4r9j/4KdhtbH55JKnKWUBLUi37y5Nw9U+L/QmnwVEn5JdHkmqABbQk1aqlb2z8ocE99oH3X5xvJkmqARbQklSLWlvhystg2RtZf8i2cMk/ZfOfJUlbxQJakmrR7dfA023uNPixL8B2I/LLI0k1xAJakmrN7Jnwx18U+yefA/sfnFscSao1FtCSVEtWr4SffANaWrL+HvvAmRfmm0mSaowFtCTVkt/8ABYvzNoDBsLHvwhNTflmkqQaYwEtSbVi2v3wyN3F/oWfgRE75hZHkmqVBbQk1YKlr8Ovv1fsH3kCHPqO/PJIUg2zgJakapcS/Oq7sGJZ1t9uBJz3t/lmkqQaZgEtSdXuwdvhiUeL/Ys+CwMH5ZdHkmpcWQvoiDg5ImZGxKyI+GIH+7eNiBsj4omIeDoiLipnHkmqOa8tgqsuL/aPezfsOy6/PJJUB8pWQEdEI/AD4BRgP+C8iNiv3WGfAJ5JKb0VmAB8KyK8TZYklWLD1I21q7P+yNHw3o/km0mS6kA5r0AfBsxKKc1OKa0DJgFntDsmAUMiIoDBwOvA+jJmkqTa8fBdxbsNRsBH/h769c83kyTVgUgplefEEWcDJ6eUPlbofxA4PKV0aZtjhgA3APsAQ4D3p5Ru7uBcFwMXA4wcOfKQSZMmlSXz5qxYsYLBgwfn8t7qHY5xfaiFce6zZiWH3fAT+qxbA8C8vQ9h1qEn5pyqstTCOGvzHOfal+cYT5w4cXpKaXz77eVcXT862Na+Wj8JeBw4DtgTuCMi7k8pLdvoRSldAVwBMH78+DRhwoQeD1uKKVOmkNd7q3c4xvWhJsb5x1+HQvHM9juw8ye/zM79B+SbqcLUxDhrsxzn2leJY1zOKRzzgF3a9HcGFrQ75iLg2pSZBbxIdjVaktSZxx+BqfcW+x/8FFg8S1KvKWcBPRUYGxG7Fz4YeC7ZdI225gDHA0TESGBvYHYZM0lSdVuzCn7z/WL/yBPggE3+d1GSVEZlm8KRUlofEZcCk4FG4MqU0tMRcUlh/+XAV4FfRMSTZFM+vpBSerVcmSSp6t3wG3ij8GtyyLbw/ovzzSNJdaicc6BJKd0C3NJu2+Vt2guAd5YzgyTVjHkvwp1/LPbPuRgGb5NfHkmqU96JUJKqQUrZ1I3W1qz/lgPhiOPyzSRJdcoCWpKqwUN3wvNPZ+3GRrjg0mztZ0lSr7OAlqRKt2I5XP3TYv/E98Ko3fLLI0l1zgJakird9b+E5Uuz9nYj4PQP5JtHkuqcBbQkVbL5f4UpbT6Lfe4l3q5bknJmAS1JlSolmPRjSIUPDu53MIw7Kt9MkiQLaEmqWE88Cs/OyNrRkK357AcHJSl3FtCSVInWN8MfflLsH3sqjB6TWxxJUpEFtCRVortvgEXzs/aAQXDGB/PNI0l6kwW0JFWa5Uvgxt8W+6efn922W5JUESygJanSXP9rWL0ya48cDcednm8eSdJGLKAlqZLM/yvc22bZunMuhqY+ucWRJG3KAlqSKskffrbxsnUHHZZvHknSJpq62hkR+wBnAKOBBCwAbkgpPdsL2SSpvsz8Mzw1NWtHwDkfd9k6SapAnV6BjogvAJOAAP4ETC20r4qIL/ZOPEmqEynBNVcW+0ceDzvvnl8eSVKnuroC/VFg/5RSc9uNEfFt4GngG+UMJkl1ZcZDMPu5rN3Ux2XrJKmCdTUHuhUY1cH2nQr7JEk9oaUFrv1FsT/xNNh+ZG5xJEld6+oK9GeAuyLieWBuYduuwF7ApWXOJUn146E74OXCr9kBA+HUc/PNI0nqUqcFdErptoh4C3AY2YcIA5gHTE0ptfRSPkmqbevWwg2/LvZPOtubpkhShetyFY6UUivwSC9lkaT6c/cN8MarWXubYXDiWfnmkSRtlutAS1JeVq2AW35X7J/+AejXP788kqSSWEBLUl4mX50V0QAjdoJjTsk3jySpJBbQkpSHFcvgzuuL/TMvhKYuZ9VJkipESQV0RPy+7bMkaSvdfg2sXZ21R+0Kh74j3zySpJKVegV6r8Lz2HIFkaS6sXwJ3NXm6vPpF0BDY25xJEnd4xQOSeptk6+BtWuy9ugxcMjbc40jSeoeC2hJ6k3LlmRL123w7gugwV/FklRN/K0tSb1p8h+ym6cA7Lw7jDsq3zySpG4rtYCOsqaQpHqw9A2456Zi36vPklSVSv3NfVm7Z0lSd93W5urzLnt69VmSqlRJBXRK6bdtnyVJ3bT0dbj35mL/jAsg/M89SapG/t+hJPWG268pXn3ebSy89Yh880iStpgFtCSV24plMKXN1efTP+DVZ0mqYhbQklRud12/8brPBx2eaxxJ0tZp2twBEfHZDjYvBaanlB7v8USSVEvWrNr4roOnvt+VNySpypXyW3w8cAkwuvC4GJgA/CQiPl++aJJUA+65GVatyNo7jIJD35FvHknSVtvsFWhge+DglNIKgIj4V+Bq4B3AdOC/yxdPkqrYurVwx7XF/innQENjfnkkST2ilCvQuwLr2vSbgd1SSquBtWVJJUm14IHJsOyNrD1sOBx5fL55JEk9opQr0L8FHomIDZP4TgeuiohBwDNlSyZJ1Wz9erjt6mL/pLOhqU9+eSRJPWazBXRK6asRcStwNNktvS9JKU0r7D6/nOEkqWo9eje8/krWHrItHHNyvnkkST2mlCvQADOABRuOj4hdU0pzypZKkqpZawvc8rti/8T3QL/++eWRJPWoUpax+yTwr8AioIXsKnQCDipvNEmqUjMehkXzs/aAQTDh9HzzSJJ6VClXoD8N7J1Seq3cYSSp6qUEk9vMfZ54OgwclF8eSVKPK2UVjrlkN06RJG3O80/D7OeydlMfOP7d+eaRJPW4Uq5AzwamRMTNtFm2LqX07bKlkqRqNfkPxfaRx8O22+WXRZJUFqUU0HMKj76FhySpIwvmwBOPFvvvfG9+WSRJZVPKMnb/1htBJKnq3X5Nsf22I2CnXfLLIkkqm04L6Ij4n5TSZyLiRrJVNzaSUnJinyRtsOQ1eOTuYv+ks/PLIkkqq66uQP9f4fmbvRFEkqra3TfC+uasvcc+sNf++eaRJJVNpwV0Sml64fne3osjSVVozWqYclOxf9LZEJFfHklSWZVyI5Wjga8AuxWODyCllPYobzRJqhIPTIZVK7L2DqNg3JH55pEklVUpq3D8DPg7YDrZnQglSRu0tMAd1xb7J54FDY355ZEklV0pBfTSlNKtZU8iSdVo2n3w2itZe/C2cPSJ+eaRJJVdKQX0PRFxGXAtG99I5bGypZKkatD+tt3HnQ59++WXR5LUK0opoA8vPI9vsy0Bx/V8HEmqIs89AXNeyNp9+8HE0/PNI0nqFaXcSGVibwSRpKpzxx+L7aNOhCHb5pdFktRrSlmF4wXgEeB+4L6U0jNlTyVJlW7RfPhzm9t2n3BmblEkSb2roYRj9gN+DGwPfDMiZkfEHzfzGkmqbXddX2wfdBjsuHN+WSRJvaqUAroFaC48twKLgFfKGUqSKtqqFfDg7cW+V58lqa6U8iHCZcCTwLeBn6SUXitvJEmqcPdPhrVrsvboMbDvuFzjSJJ6VylXoM8D7gP+HzApIv4tIo4v5eQRcXJEzIyIWRHxxU6OmRARj0fE0xHhbcMlVbaWlo2nb5xwprftlqQ6U8oqHNcD10fEPsApwGeAzwMDunpdRDQCPwBOBOYBUyPihrYfQoyIocAPgZNTSnMiYoct/DokqXfMeAheb3PjlMNdqEiS6s1mr0BHxDWFlTi+CwwGLgSGlXDuw4BZKaXZKaV1wCTgjHbHfAC4NqU0ByCl5NxqSZXtzjafoZ5wqjdOkaQ6VMoc6G8Aj6WUWrp57tHA3Db9eRRvyrLBW4A+ETEFGAJ8N6X0q/YnioiLgYsBRo4cyZQpU7oZpWesWLEit/dW73CM68OWjvOQVxdyyKzsP9FaGxp4pO/2rPPnpWL557k+OM61rxLHuMsCujCl4jTgCxGRgGeAH6aUFpVw7o4mBaYO3v8Q4HiyKSEPR8QjKaW/bPSilK4ArgAYP358mjBhQglv3/OmTJlCXu+t3uEY14ctHuef/NebzYbDJ3LUqaf1XCj1OP881wfHufZV4hh3OoUjIo4Gpha6vwJ+XWg/Wti3OfOAXdr0dwYWdHDMbSmllSmlV8k+rPjWUoJLUq9641WYdl+xf/yZuUWRJOWrqyvQ3wLOTCnNaLPt+sJNVH7MptMx2psKjI2I3YH5wLlkc57buh74fkQ0AX0L5/xON/JLUu+456ZsBQ6AsQfAmLH55pEk5aarAnqbdsUzACmlxyNiyOZOnFJaHxGXApOBRuDKlNLTEXFJYf/lKaVnI+I24M9kN2n5aUrpqS36SiSpXNaugftuKfZPfE9+WSRJueuqgI6IGJZSeqPdxu0obf1oUkq3ALe023Z5u/5lwGWlxZWkHDx6D6xYlrWHj4S3HZFvHklSrroqhL8D3B4Rx0bEkMJjAnArTrOQVC9S2vjGKce9Gxoa88sjScpdp1egU0pXRMQC4KvA/mQraDwD/EdK6cZeyidJ+Xr+KZj/16zdtx+8/aRc40iS8tflMnYppZuAm3opiyRVnrZXn488HgYOzi+LJKkilDSXWZLq0uuLs1t3b3Dcu/PLIkmqGBbQktSZe2+G1tasvfdBMHpMrnEkSZVhswV0YR3nzW6TpJrSvA7uu7XY9+qzJKmglCvQ13Sw7eqeDiJJFWXafbB8adYeNhzedmS+eSRJFaPTDxFGxD5kq29sGxFntdm1DdC/3MEkKVd3t1lsaMJp0OjSdZKkTFercOwNnAYMBU5vs3058PEyZpKkfM2eCS/OzNpNfeAdJ+ebR5JUUbpaB/p64PqIODKl9HAvZpKkfN1zQ7F96LEwZGhuUSRJlaeUOdBzI+KPEfFKRCyKiGsiYueyJ5OkPCxbAlPvK/aP98ODkqSNlVJA/xy4ARgFjAZuLGyTpNpz362wvjlr77EPjHlLvnkkSRWnlAJ6h5TSz1NK6wuPXwAjypxLknpfS0u29vMGLl0nSepAKQX04oi4ICIaC48LgNfKHUySet2Mh+CNV7P2kKFwyNtzjSNJqkylFNAfAc4BXgYWAmcXtklSbbm7zYcHjz0V+vTNL4skqWJ1tYwdACmlOYD/jympts2dDX95Mms3NmYFtCRJHejqRir/0sXrUkrpq2XII0n5uKfNjVPGHZ3dfVCSpA50dQV6ZQfbBgEfBbYHLKAl1YaVy+GRu4t9PzwoSepCVzdS+daGdkQMAT4NXARMAr7V2eskqeo8eDusW5u1d9kDxu6fbx5JUkXrcg50RGwHfBY4H/glcHBK6Y3eCCZJvaK1Be65qdifeDpE5JdHklTxupoDfRlwFnAFcGBKaUWvpZKk3vLUNFi8MGsPHAyHT8w3jySp4nW1jN3fk9198MvAgohYVngsj4hlvRNPksrsrjZL1x1zEvTrn18WSVJV6GoOdClrREtS9Xp5Hjw9PWtHwITT8s0jSaoKFsmS6lfbpesOOgxG7JRfFklS1bCAllSf1qyCh+4o9l26TpJUIgtoSfXp4btg9aqsvePOsO+4fPNIkqqGBbSk+pMS3N1m+saE06DBX4eSpNJ0+2+MiLgzIm6NCD9tI6kqDV30Eiyck3X6DYCjT8w3kCSpqnR5I5VOXAjsBBzRw1kkqVeMnvlYsXPU8TBgUH5hJElVp9sFdEppAbAAmN7zcSSpzF57heHzni/2J56eXxZJUlXaokl/EXFrTweRpF5x7y1ESll7n7fBqN1yjSNJqj5d3cr74M52AW8rSxpJKqfmdXBfm3//H+fVZ0lS93U1hWMqcC9Zwdze0LKkkaRymnYfrFiatbcbAW/1oxySpO7rqoB+FviblNLz7XdExNzyRZKkMmm/dF1jY35ZJElVq6s50F/pYv8nez6KJJXR7Jnw4kwAWhsa4ZiTcg4kSapWnV6BTild3cW+68qSRpLK5Z7i1edXxuzLjkOG5pdFklTVvPWWpNq3fAlMvffN7vy9D8kviySp6llAS6p9998G65uz9h77sHz7nfLNI0mqahbQkmpbSwtMubnY98YpkqStVFIBHRH7tH2WpKrxxCPw+uKsPWRbGH9MvnkkSVWv1CvQv233LEnVoe3Sde84Bfr0zS+LJKkmdHcKR0c3VZGkyrTgJXju8azd0ADHvivXOJKk2uAcaEm1q83SdYw7Krv7oCRJW8kCWlJtWrUSHrqz2PfDg5KkHtLdAjqVJYUk9bSH74S1a7L2qN1g74PyzSNJqhmlFtDR7lmSKldrK9x9Q7F/3Lsh/PUlSeoZpRbQx7R7lqTK9ezjsGh+1h4wEI44Ltc4kqTaUlIBnVJa0fZZkiraPW2uPh/9Tug/IL8skqSa44cIJdWWV1+GJx4t9iecll8WSVJNsoCWVFum3Ayp8Hnn/Q+BHXfON48kqeZstoCOiE+Xsk2ScrduLdx/W7F/3LvzyyJJqlmlXIH+UAfbPtzDOSRp6/1pCqxcnrWH7wgHjs81jiSpNjV1tiMizgM+AOweEW0+kcMQ4LVyB5OkbkkJ7ryu2J/wLmhozC2OJKl2dVpAAw8BC4HhwLfabF8O/LmcoSSp22b+Gea9mLX79oNjTs43jySpZnVaQKeUXgJeAo7svTiStIXaXn0+6kQYNCS3KJKk2lbKhwiPiIipEbEiItZFREtELOuNcJJUksUL4YlHiv3jz8gviySp5pXyIcLvA+cBzwMDgI8B3ytnKEnqlrtv2Hjpup12yTePJKmmdTUH+k0ppVkR0ZhSagF+HhEPlTmXJJVmzSp4YHKxf8KZuUWRJNWHUgroVRHRF3g8Iv6b7IOFg8obS5JK9NCdsHpV1t5x5+wKtCRJZVTKFI4PAo3ApcBKYBfgveUMJUklaW2Fu64v9o8/Axq8waokqbw2+zdNSumllNLqlNKylNK/pZQ+m1KaVcrJI+LkiJgZEbMi4otdHHdo4cOJZ3cnvKQ699RUWDQ/aw8YBEeekG8eSVJd6OpGKk8CqbP9KaWDujpxRDQCPwBOBOYBUyPihpTSMx0c91/A5E3PIklduLPN1edjTob+A/LLIkmqG13NgX4fsHorzn0YMCulNBsgIiYBZwDPtDvuk8A1wKFb8V6S6s38v8Izj2XtaIDjTs81jiSpfnRVQP82pXRwRPxfSumDW3Du0cDcNv15wOFtD4iI0cB7gOPoooCOiIuBiwFGjhzJlClTtiDO1luxYkVu763e4RhXj7c8chujCu3FO+/F0089BzxX0msd5/rgONcHx7n2VeIYd1VA942IDwFHRcRZ7XemlK7dzLmjg23tp4T8D/CFlFJLREeHv/leVwBXAIwfPz5NmDBhM29dHlOmTCGv91bvcIyrxIrl8PvvvNkdcd7HmfCWA0t+ueNcHxzn+uA4175KHOOuCuhLgPOBoUD7/xtNwOYK6HlkK3ZssDOwoN0x44FJheJ5OHBqRKxPKV23mXNLqmf33wrr1mbtXfaEsQfkm0eSVFc6LaBTSg8AD0TEtJTSz7bg3FOBsRGxOzAfOBf4QLv32H1DOyJ+Adxk8SypSy0tcM+Nxf4JZ0AX/4MlSVJP2+yNVLaweCaltD4iLiVbXaMRuDKl9HREXFLYf/mWnFdSnZvxILy+OGsP2RYOm5BrHElS/SnpVt5bKqV0C3BLu20dFs4ppQ+XM4ukGnHHH4vtCadBn775ZZEk1SVv2SWpesx6Bl54Nms39YFjT803jySpLpV0BToihgFjgf4btqWU7itXKEnq0OSri+0jJsLQ7fPLIkmqW5stoCPiY8CnyVbReBw4AniYbO1mSeodi+bD4w8X+ydusrqmJEm9opQpHJ8mu8nJSymlicA4YHFZU0lSe3f8EVJhKfkDD4XRY3KNI0mqX6UU0GtSSmsAIqJfSuk5YO/yxpKkNpYvgQdvL/ZPOju3KJIklTIHel5EDAWuA+6IiDfY9IYoklQ+U26G5nVZe9e9YO+D8s0jSaprpawD/Z5C8ysRcQ+wLXBbWVNJ0gbr1sLdNxT7J73XG6dIknLVrXWgU0r3liuIJHXo4btg+dKsvd0OcMgx+eaRJNU914GWVLlaW+H2a4r9E86EprLe/0mSpM2ygJZUuf78aLZ8HcCAQfCOk/PNI0kSFtCSKtnkNlefjz0V+g/ML4skSQWl3EhlOZDabV4KTAP+PqU0uxzBJNW52c/B809l7cZGOP6MfPNIklRQymTCb5MtW/dbIIBzgR2BmcCVwIRyhZNUx9retvuwiTBseH5ZJElqo5QpHCenlH6cUlqeUlqWUroCODWl9DtgWJnzSapHL8+Fxx4s9t/pbbslSZWjlAK6NSLOiYiGwuOcNvvaT+2QpK132x+Kt+0+4FDYZY9880iS1EYpBfT5wAeBV4BFhfYFETEAuLSM2STVo9cXZ2s/b/Cu9+eXRZKkDpRyJ8LZwOmd7H6gZ+NIqnu3XwMtLVl77P4w9oB880iS1E4pq3CMAD4OjGl7fErpI+WLJakuLV8C991a7J96bm5RJEnqTCmrcFwP3A/cCbSUN46kunbX9bBubdbeZU84YHy+eSRJ6kApBfTAlNIXyp5EUn1bvRLuuqHYP/X9EJFfHkmSOlHKhwhviohTy55EUn2bcktWRAOMHA2HHJ1vHkmSOlFKAf1psiJ6dUQsi4jlEbGs3MEk1ZF1a+GOa4v9k98HDY355ZEkqQulrMIxpDeCSKpjD94By97I2sOGw5HH55tHkqQudFpAR8Q+KaXnIuLgjvanlB4rXyxJdWN9M9z6+2L/pLOhqU9+eSRJ2oyurkB/FrgY+FYH+xJwXFkSSaovD94Br7+StYdsC8ecnG8eSZI2o9MCOqV0ceF5Yu/FkVRX1jfDzZOK/ZPfB/3655dHkqQSlLKMHRFxFJveSOVXZcokqV60v/o84bR880iSVIJS7kT4f8CewOMUb6SSAAtoSVtufTPc0ubq80lne/VZklQVSrkCPR7YL6WUyh1GUh156E54rXD1efC2MPH0fPNIklSiUtaBfgrYsdxBJNWR9c1w81XFvlefJUlVpJQr0MOBZyLiT8DaDRtTSu8uWypJtW2Tq8/OfZYkVY9SCuivlDuEpDqyydzn90L/AfnlkSSpm0opoPcE7k8pPV/uMJLqwEN3wquLsvbgbZz7LEmqOqUU0GOACyJiN2A6cD9ZQf14GXNJqkXN6+DG3xT77/TqsySp+mz2Q4QppX9JKR0HHAA8AHyOrJCWpO6550Z449WsPWQoHOdHKSRJ1aeUdaC/DBwNDAZmAP9AdhVakkq3eiXc8rti/7TzvPosSapKpUzhOAtYD9wM3As8klJaU9ZUkmrP7dfCimVZe/sd4B2n5JtHkqQtVMoUjoOB44E/AScCT0bEA+UOJqmGLF+SFdAbnHEh9OmbWxxJkrZGKVM4DgCOAY4luyvhXJzCIak7bvkdrF2dtUftBkdMzDePJElboZQpHP9FNnXjf4GpKaXm8kaSVFNeewXuuanYf8+HoKExvzySJG2lUgro9wJ7AQloBCygJZXuxl9nN08B2GMfeNuR+eaRJGkrdToHOiKaIuK/yaZs/BL4NTA3Iv47Ivr0VkBJVWzhXHjwzmL/rIsgIr88kiT1gK4+RHgZsB2wR0rpkJTSOLK7Eg4FvtkL2SRVu2t+Bqk1a+93MOzz1nzzSJLUA7oqoE8DPp5SWr5hQ0ppGfC3wKnlDiapyj37ODz+SNaOyK4+S5JUA7oqoFNKKXWwsYVsPrQkday1BX7/k2L/yONhzNj88kiS1IO6KqCfiYgL22+MiAuA58oXSVLVe+hOmPtC1u7bD97z4VzjSJLUk7paheMTwLUR8RFgOtlV50OBAcB7eiGbpGq0ZjX88RfF/klnw7DhucWRJKmndVpAp5TmA4dHxHHA/kAAt6aU7uqtcJKq0G1/gKVvZO2h28PJ78s3jyRJPWyz60CnlO4G7u6FLJKq3euL4fZriv33fBj69c8tjiRJ5dDVHGhJ6p5rfw7r1mbtXffKPjwoSVKNsYCW1DNmz4RH2vxn1TkfhwZ/xUiSao9/u0naeq0t8NvvF/vjjvKmKZKkmmUBLWnr3T8Z/vp81m7qk119liSpRllAS9o6K5Zlc583OOUcGLFTfnkkSSozC2hJW+fan8PK5Vl7+I5ZAS1JUg2zgJa05V6cCfffVuyfe0l250FJkmqYBbSkLdPaAr/5PqSU9Q86DN52RL6ZJEnqBRbQkrbM3Tdu/MHBcy/JN48kSb3EAlpS9732CvzxF8X+u86FHUblFkeSpN5kAS2pe1LKpm6sXZP1R+3qBwclSXWlrAV0RJwcETMjYlZEfLGD/edHxJ8Lj4ciwjsvSJVu2v3w5z8V+xd+JpvCIUlSnShbAR0RjcAPgFOA/YDzImK/doe9CBybUjoI+CpwRbnySOoBK5fDVT8q9iecBnu1/2MtSVJtK+cV6MOAWSml2SmldcAk4Iy2B6SUHkopvVHoPgLsXMY8krbW738Cywp/ZIduD2ddlG8eSZJyEGnDElQ9feKIs4GTU0ofK/Q/CByeUrq0k+P/Adhnw/Ht9l0MXAwwcuTIQyZNmlSWzJuzYsUKBg8enMt7q3c4xp3bbt4sDppy9Zv9p97xHl7dde8cE205x7k+OM71wXGufXmO8cSJE6enlMa3395UxveMDrZ1WK1HxETgo8DbO9qfUrqCwvSO8ePHpwkTJvRQxO6ZMmUKeb23eodj3IkVy+HGNjOsDj2WAy78m/zybCXHuT44zvXBca59lTjG5Syg5wG7tOnvDCxof1BEHAT8FDglpfRaGfNI2lKTfgRLX8/a2wyD8z+Rbx5JknJUzjnQU4GxEbF7RPQFzgVuaHtAROwKXAt8MKX0lzJmkbSlZjwEj9xd7H/wUzB4m/zySJKUs7JdgU4prY+IS4HJQCNwZUrp6Yi4pLD/cuBfgO2BH0YEwPqO5plIysnyJfCr/y32jzwexh2ZWxxJkipBOadwkFK6Bbil3bbL27Q/BmzyoUFJFSAl+MV3siIaslU3vF23JEneiVBSJ6bcBE88Wux/6DMwaEhucSRJqhQW0JI2Nf+v2ZrPG5xwJhx4aF5pJEmqKBbQkjbWvA6u+Eb2DLDz7vDej+SbSZKkCmIBLWljf/hpdgUaoE9fuPiL2bMkSQIsoCW1Ne0+uLvNapPv/xsYtVt+eSRJqkAW0JIyL8+Fn3+n2B93FBx7an55JEmqUBbQkmDtGvjR12Dt6qw/Yie46LOQrc8uSZLasICW6l1K8OvvFec9N/WBv/0yDBycayxJkiqVBbRU7+69BR6+q9g//xOw65755ZEkqcJZQEv17Lkn4KofFvtHnQhvPym/PJIkVQELaKleLX4ZLv8atLRk/V33zK4+O+9ZkqQuWUBL9WjNKvj+V2DFsqy/zTC49CvQr3+eqSRJqgoW0FK9aW2Fn31z4w8NfuKfYbsRucaSJKlaWEBL9eYPP4EZDxX7H/wk7LlffnkkSaoyFtBSPbn9Wrjjj8X+iWfB0e/ML48kSVXIAlqqF9Puy64+b3Dw0fC+j+aXR5KkKmUBLdWDvzwFP70su2kKwF77wcc+Dw2N+eaSJKkKWUBLte6l5+F7/wLrm7P+jjtnK2707ZdrLEmSqpUFtFTL5v8VvvNPsHpV1t9mGHz6P2DwNrnGkiSpmllAS7Vq0QL49j8W13oeOBj+7mswYsd8c0mSVOUsoKVa9Nor8O0vwtLXs36/AfCZr8Eue+SbS5KkGmABLdWaxS/DZZ/LimiAPn3hU/8Ge+ydby5JkmpEU94BJPWgRfPhW1+E1xdn/cYm+MS/wN4H5ZtLkqQaYgEt1YoFc7LiecO0jaY+WfF8wPh8c0mSVGMsoKVa8NLz8D//DMuXZP2+/eDSf4X9Ds41liRJtcgCWqp2zzwGP/gqrF2d9fv1h0/9u9M2JEkqEwtoqZo9eg9c+S1oWZ/1Bw7Oiue99ss3lyRJNcwCWqpGKcHt18IfflLcNmw4fOY/YPSY3GJJklQPLKClatO8Dn7zfXjg9uK2Ubtm6zxvNyK/XJIk1QkLaKmaLF8CP/wqPP90cdvY/eETX4HBQ/JKJUlSXbGAlqrFnBfgh/8Ory4qbjvyBLjwU9nNUiRJUq+wgJYqXUpw/23w2x/C+uZsWwSc/VF453uztiRJ6jUW0FIlW7smm+/80J3Fbf0Hwse/AG89PL9ckiTVMQtoqVLNnQ0/+S9Y8FJx2+gx8Lf/BDvuklssSZLqnQW0VGlaW+HOP8K1vyhO2QA46gQ4/9LsRimSJCk3FtBSJXl9cXZjlOceL27r2w/OvQSOOdn5zpIkVQALaKkStLbCvbfANVfCmlXF7buNhY99HnZyyoYkSZXCAlrK28tz4ZffheefKm6LBjj1HDj9fGjqk182SZK0CQtoKS9r18Ctv4Pbrt54rvPI0fDhz2Y3SJEkSRXHAlrqbSnB1HvhDz+FN14tbm9ogJPPgdM/4I1RJEmqYBbQUm+aMwuuunzj6RoAY8bChZ+BXffMJZYkSSqdBbTUGxbNh+v/L7vynFJx+5Ch8N6L4KgTsyvQkiSp4llAS+X02iK48bfw0B3ZShsbNDbC8WfCaR+AgYNyiydJkrrPAloqh8Uvw+1Xw/2TN/6AIGS34D77Yy5NJ0lSlbKAlnrSS8/D5Kth6v2QWjfet+/b4MwPwZ775hJNkiT1DAtoaWu1tsLT0+H2a+DZxzfdv+e+8J4Pwz5v7e1kkiSpDCygpS219A148Ha471Z49eVN9+87Dk4+G/Y72FtwS5JUQyygpe5obYHnnsjmNj/2ILSs33h/NMChx8BJZ2e34ZYkSTXHAlranJRg3ovwyF3w6BRY8tqmxwwcDEefCMedASN27PWIkiSp91hASx1JCRa8BDMegqn3wfy/dnzcnvvChHfBIcdA3369GlGSJOXDAlraoLWVbRbPh6t/lhXOi+Z3fNyQbeHQY+GYk2GXPXo3oyRJyp0FtOrbkteyFTSeng7PzODgFcs6Pq5vPxh3FBxxXPbhwCb/6EiSVK+sAlRfViyHF56BmX/OiubOpmYA9OsPBx6aFc5vPRz6D+y1mJIkqXJZQKt2pQSvL4bnn4Lnn4ZZT3ddMAPr+g+k78FHw8FHw37jnNcsSZI2YQGt2pBSNh3jpefhpVnZ81+fh2VvdP26xiYYuz/sfwjsfwgPvTCHCRMn9k5mSZJUlSygVX1Wr4QFc7JVMjY85r2Y3dhkcxoaYNe9sqJ537fBWw6C/gOK+2fPLVtsSZJUGyygVZnWN8Ori7I7/L2yABYvhIVzs2L59cWln6ffgGypubH7w177wx77ZHObJUmStpAFtPLRvA7eeDWbdvHGq1lRvHghvLIQFi+A11+F1Nq9c/YbALvumd0BcMzY7Hnk6OyqsyRJUg+xgFbPaW2FVSth+RJYvjR7XrYEViwtFspvvApvvJZt21KNTVlhPGo3GL1b9jxqNxg5Choae+iLkSRJ6pgFtDbVvA5WrciK4VUrsjnHbftvPlbCyuWFgnkJrFgGLS09kyEChg2HETtljx12gh1GZYXyDqNdh1mSJOWmrFVIRJwMfBdoBH6aUvpGu/1R2H8qsAr4cErpsXJmqlopZcXp+ubio3kdNDe329YMzWth7ZrssW5De3WbbWuK7bbb1qyBNauy85RbQwMM3T4rkoduD8NGwIgdC49RMHwk9Olb/hySJEndVLYCOiIagR8AJwLzgKkRcUNK6Zk2h50CjC08Dgd+VHiuLFPvg9UrGfWXZ6FlaVbIthYeLS3Z1IWWtv127ZbW0va1NBcL4uZ2hfL65qyIrnQDBsLgbWGboTBkaHbb6yFDs/52I7KCedhwrnthFZfd8TwLlqxmVAzgc0fszZnjRpc93peve5KrHp1LS0o0RnDe4bvwH2ce2Onx182Yz2WTZ2Y5hw7gcydtPmdnr6m07b2hu9/vzuT5NUjl4s+1VL3KeQX6MGBWSmk2QERMAs4A2hbQZwC/Sikl4JGIGBoRO6WUFpYxV/dd/VN47RXeAvCnvMP0gsYmGDgYBg4qPg8YzNRX1jJ9cTPLol/2aOjH0ujHYW/dnU+cPj4rlku4anzdjPl86bqnWd2cTfeYv2Q1X7r2SYCy/uXx5eue5NePzHmz35LSm/2OirrrZsznS9c+2a2cnb1m2kuvc830+RWzvauvoad09/vdmS0ZB6nS+XMtVbdyFtCjgbaL6s5j06vLHR0zGqisArpSPpjW0ABNfbJHnz5Zsbqh39Sm36dPtiJFv37Qt3+2bFu//lm7/4Biu1+/wr4BG+/v2y+bg9zOuV+6hZZBm14Ff+C5xCc+NKLkL+OyyTPf/Etjg9XNLVw2eWZZ/+K46tGO13i+6tG5HRZ0W5Kzs9dsuApbKdvL/b2G7n+/O5PXz4tUTv5cS9UtUpmmBUTE+4CTUkofK/Q/CByWUvpkm2NuBr6eUnqg0L8L+HxKaXq7c10MXAwwcuTIQyZNmlSWzJ3ZfcYU+qxdzbr1LTT160uKhuzR0NFzFPsdHhNdvr61sZHWhiZS4bm1sZHWxkZSQxMp5+XYnpzf+coZB47ettfP012lvO+KFSsYPHhwycd35z0qTTm/11DZPy9tx1m1q5LHOa/fg7WoksdZPSPPMZ44ceL0lNL49tvLWUAfCXwlpXRSof8lgJTS19sc82NgSkrpqkJ/JjChqykc48ePT9OmTStL5s2ZMmUKEyZMyOW9K8GeX7plk6uZAI0RvPD1U0s+z9HfuJv5S1Zvsn300AE8+MXjtipjV0rJ33aMtyRnZ69pjOj0vfPYXu7vNVT2z0u9/1muF5U8znn9HqxFlTzO6hl5jnFEdFhAl/OS5lRgbETsHhF9gXOBG9odcwNwYWSOAJZW3Pxnvem8w3fp1vbOfO6kvRnQZ+NpMQP6NPK5k/be4myl6G7+LcnZ2WvOO3yXitpe7u81VP/Pi1RO/lxL1a1sc6BTSusj4lJgMtkydlemlJ6OiEsK+y8HbiFbwm4W2TJ2F5Urj7behnmrW7uqwob5fb396fPu5t+SnF29Zvxu21XU9nKr9p8XqZz8uZaqW9mmcJSLUzhUTo5xfXCc64PjXB8c59pXb1M4JEmSpJpjAS1JkiR1gwW0JEmS1A0W0JIkSVI3WEBLkiRJ3WABLUmSJHWDBbQkSZLUDRbQkiRJUjdYQEuSJEndYAEtSZIkdYMFtCRJktQNFtCSJElSN0RKKe8M3RIRi4GXcnr74cCrOb23eodjXB8c5/rgONcHx7n25TnGu6WURrTfWHUFdJ4iYlpKaXzeOVQ+jnF9cJzrg+NcHxzn2leJY+wUDkmSJKkbLKAlSZKkbrCA7p4r8g6gsnOM64PjXB8c5/rgONe+ihtj50BLkiRJ3eAVaEmSJKkbLKAlSZKkbrCALkFEnBwRMyNiVkR8Me886hkRsUtE3BMRz0bE0xHx6cL27SLijoh4vvA8LO+s2joR0RgRMyLipkLfMa4xETE0Iq6OiOcKf6aPdJxrT0T8XeH39VMRcVVE9Hecq19EXBkRr0TEU222dTquEfGlQk02MyJOyiOzBfRmREQj8APgFGA/4LyI2C/fVOoh64G/TyntCxwBfKIwtl8E7kopjQXuKvRV3T4NPNum7xjXnu8Ct6WU9gHeSjbejnMNiYjRwKeA8SmlA4BG4Fwc51rwC+Dkdts6HNfC39PnAvsXXvPDQq3WqyygN+8wYFZKaXZKaR0wCTgj50zqASmlhSmlxwrt5WR/4Y4mG99fFg77JXBmLgHVIyJiZ+BdwE/bbHaMa0hEbAO8A/gZQEppXUppCY5zLWoCBkREEzAQWIDjXPVSSvcBr7fb3Nm4ngFMSimtTSm9CMwiq9V6lQX05o0G5rbpzytsUw2JiDHAOOBRYGRKaSFkRTawQ47RtPX+B/g80Npmm2NcW/YAFgM/L0zV+WlEDMJxrikppfnAN4E5wEJgaUrpdhznWtXZuFZEXWYBvXnRwTbX/qshETEYuAb4TEppWd551HMi4jTglZTS9LyzqKyagIOBH6WUxgEr8b/xa05hDuwZwO7AKGBQRFyQbyrloCLqMgvozZsH7NKmvzPZfxmpBkREH7Li+TcppWsLmxdFxE6F/TsBr+SVT1vtaODdEfFXsulXx0XEr3GMa808YF5K6dFC/2qygtpxri0nAC+mlBanlJqBa4GjcJxrVWfjWhF1mQX05k0FxkbE7hHRl2zi+g05Z1IPiIggmzP5bErp22123QB8qND+EHB9b2dTz0gpfSmltHNKaQzZn927U0oX4BjXlJTSy8DciNi7sOl44Bkc51ozBzgiIgYWfn8fT/bZFce5NnU2rjcA50ZEv4jYHRgL/Km3w3knwhJExKlk8ygbgStTSl/LN5F6QkS8HbgfeJLi/Nh/JJsH/XtgV7Jf2O9LKbX/cIOqTERMAP4hpXRaRGyPY1xTIuJtZB8U7QvMBi4iu0jkONeQiPg34P1kqyjNAD4GDMZxrmoRcRUwARgOLAL+FbiOTsY1Iv4J+AjZz8FnUkq39npmC2hJkiSpdE7hkCRJkrrBAlqSJEnqBgtoSZIkqRssoCVJkqRusICWJEmSusECWpIqXESsyDuDJKnIAlqSJEnqBgtoSaoSETEhIqZExNUR8VxE/KZwRzYi4tCIeCginoiIP0XEkIjoHxE/j4gnI2JGREwsHPvhiLguIm6MiBcj4tKI+GzhmEciYrvCcXtGxG0RMT0i7o+IffL8+iWpUjTlHUCS1C3jgP2BBcCDwNER8Sfgd8D7U0pTI2IbYDXwaYCU0oGF4vf2iHhL4TwHFM7VH5gFfCGlNC4ivgNcSHb31SuAS1JKz0fE4cAPgeN66euUpIplAS1J1eVPKaV5ABHxODAGWAosTClNBUgpLSvsfzvwvcK25yLiJWBDAX1PSmk5sDwilgI3FrY/CRwUEYOBo4A/FC5yA/Qr75cmSdXBAlqSqsvaNu0Wst/jAaQOjo0OtnV0ntY2/dbCORuAJSmlt21xUkmqUc6BlqTq9xwwKiIOBSjMf24C7gPOL2x7C7ArMLOUExauYr8YEe8rvD4i4q3lCC9J1cYCWpKqXEppHfB+4HsR8QRwB9nc5h8CjRHxJNkc6Q+nlNZ2fqZNnA98tHDOp4Ezeja5JFWnSKmj//WTJEmS1BGvQEuSJEndYAEtSZIkdYMFtCRJktQNFtCSJElSN1hAS5IkSd1gAS1JkiR1gwW0JEmS1A3/H4GX0GKJRZ+kAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"X = np.linspace(0, 100, 500)\n",
"P_hat = 1 / (1 + np.exp(-log_reg.params[0] - log_reg.params[1] * X))\n",
"\n",
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.scatter(df[\"Annual_Family_Income\"], df[\"Own_House\"])\n",
"ax.plot(X, P_hat, color=\"tomato\", lw=3)\n",
"ax.grid()\n",
"ax.set_xlabel(\"Income\")\n",
"ax.set_ylabel(\"Owning a flat = 1, Not Owning = 0\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "1b4562a1-d3ec-42f5-b9e5-46005bdefea3",
"metadata": {},
"source": [
"To interpret estimated coefficients, we can calculate the marginal effect by taking derivative w.r.t certain $\\beta$, in this case $\\beta_1$.\n",
"$$\n",
"\\frac{dP_i}{d X_i} = \\frac{dP_i}{dZ_i}\\frac{dZ_i}{dX_i}=\\frac{e^{-Z_i}}{(1+e^{-Z_i})^2}\\beta_2\n",
"$$\n",
"Plot both accumulative and marginal effects"
]
},
{
"cell_type": "code",
"execution_count": 121,
"id": "3a56cd37-cf4a-4da6-8a67-f09058eeb09a",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAGrCAYAAACBqSI7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAB5QElEQVR4nO3dd5xcVfnH8c+zvSabCmkQOgRIQgmhE1pIIHSkS1FAVFRQUFDsDRVFEQUpivyAANKMUgIBYugQegklQCCVlM1me53z++Pczcxutu/M3NmZ7/v1mteec++de58pO/vsmeeea845REREREQkObLCDkBEREREJJMoARcRERERSSIl4CIiIiIiSaQEXEREREQkiZSAi4iIiIgkkRJwEREREZEkUgIu0gkzW2Jmh/XxvgeY2fvxjimezGy8mTkzywk7llhm9oiZnR2H/Uwzs2V9vK+Z2T/MbL2ZvRQs+6qZfW5m1WY2rL/x9VXYcbR/b/fn90SSL3jfbB12HCKZTgm4JJyZzQ8SmfywY0mUIJHdtrXvnHvaObdDmDH1V5BYNZrZ8HbLXw8e7/hEHNc5N9M5989E7DtW8BhqgoSk9fbdYPX+wOHAWOfcXmaWC/wBmO6cK3HOrevjMfv1T093ccTsP/YxvdGXY3UmXu9tM/uJmd0ej5jSgZmdY2ZvmVmtma0ys+vNrCzexwneNx/He7+90dnvgZndama/CCsukWRSAi4JFSRpBwAOOCbcaKQPPgFOa+2Y2a5AYV93lmqj7cCkICFpvf02WL4lsMQ5VxP0NwMKgHdCiTKqp3GUxTymSUmIS/rBzL4D/Aa4DBgM7I1/Dz5uZnlhxiYiiaEEXBLtLOAF4FagTVmBmY0zs/vNbI2ZrTOz62LWnW9mi8ysyszeNbPdg+VtRppjR0xaSw7M7LtmttrMVprZcWZ2pJl9YGblZvb9ju4be/+OHoSZ7WVmz5tZRbDf61r/MJrZgmCzN4IRx1Ni92Vml5vZve329yczuzZoDzazW4L9LjezX5hZdm/jiHl+LjSzD4NvHf5iZhasyzazq81srZl9DBzV8UvWxv/hX8NWZwO3tYvpKDN7zcwqzWypmf0kZl3rSNeXzewz4Mkgjt8HcXxiZhfFjoaZ/8bkvKB9jpk9E8S9Pth+Zsz+z415n3xsZl/pwWPqkpl9GbgZ2Cd4PWcDrSUXFWb2ZLDdjmb2ePC+et/MTo7ZR2HwGD81sw3BYygEFsTsp9rM9ung+Plm9kczWxHc/hgs276jOHr4mHryvvla8L6pMrOfm9k2wX0qzeyemPd7h78nZra5+dHbYTHL9jD/+53bgxg7fe8G6zv7TNgpeM9UmNk7ZnZMzH1uNbO/mi9rqjazZ4M4/xgc4z0z2y1m+9Fmdl8Q8ydm9s0u4h1sZrcF235qZleaWVawrsv3bbv9DAJ+CnzDOfeoc67JObcEOBmfhJ9pZgVmVmfBt1HBsZqD+2L+M+OPMY/5L2b2UPBcvWhm27R7nrft4bbTg/f2huB5/J8Fv5sdPI4u32N9YWbHBK9pRfAa79TR44h5LK1/C4ab2X+D+5Wb2dMxr02PX2ORhHLO6aZbwm7AYuBrwB5AE7BZsDwbeAO4BijGj+rtH6z7ArAcmAIYsC2wZbDOAdvG7P9W4BdBexrQDPwIyAXOB9YAdwKlwM5APbB1+/vG3H9ZTH8JcFjQ3gM/KpUDjAcWARfHbNs+ro37wv8RrQUGxTz2lcDeQf9B4G/B8zASeAn4SifPZ0/i+C9QBmwRPP4ZwboLgfeAccBQ4Klg+5xOjrUEOAyf9O0UxL00eDwOGB/zWHfF/0M/EfgcOC5YNz7Y9rbg8RUGcbwLjAWGAPNi4wDmA+cF7XPw75vzg+N/FVgBWLD+KGAb/PvkoOB53r2j17ODx9fmNWu37hzgmZh+6+NojbE4eC7ODV6L3YG1wM7B+r8Ej2NMEPe+QH77/XRy7J/h/2kdCYwAngN+3lEcHdy3w/X07H0zBxiE/z1pAJ4AtsaPyL4LnN2D35OHga/GrLsG+HMnsf4EuL2H790OPxPwv+eLge8DecAhQBWwQ8zv+Nrg8RcAT+K/1TkreF1+ATwVbJsFvIL//MgLHvvHwBGdxH8b8G/8Z8t44APgyz1537bbzwz859YmrynwT2B20F4AnBi0HwM+AmbGrDs+5jGXA3sFr/cdwF0dve+72hYYDlQCJwTrvhU8pvP68tnUw/fprUQ/z7cHavClYLnAd4PXOq8Hfwt+DdwQ3C8X/y2s9fY11k23RN40Ai4JY2b74/9I3uOcewX/B+P0YPVewGjgMudcjXOu3jn3TLDuPOC3zrmXnbfYOfdpDw/bBPzSOdcE3IX/I/In51yVc+4d/Ff3E3v7WJxzrzjnXnDONTs/OvU3fMLXk/t+CrwKHBcsOgSodc69YGabATPxf6hqnHOr8UnLqf2I4yrnXIVz7jN8kj05WH4y8Efn3FLnXDn+j1RPtI6CH45P4Je3i2m+c+4t51zEOfcmMLuDmH4SPL66II4/OeeWOefWA1d1c/xPnXM3Oeda8AnJKHwpBs65h5xzHwXvk//hE5MDevi4AF4NRslab0f08H6z8CUq/whei1eB+4CTgpG2LwHfcs4td861OOeec8419HDfZwA/c86tds6twY+OfrEXjwlgbcxjurSH75vfOOcqg9+Tt4HHnHMfO+c2AI8Au9G9fwJngv/GBV++9H+9iLuz925nnwl7AyXB/Rqdc0/ik/jTYvb5QPD464EHgHrn3G3B++numMc1BRjhnPtZsK+PgZvo4HcxeGynAFcEny1LgN/T9nXq9H3bznBgrXOuuYN1K4P1AP8DDjL/TdFE4NqgXxDE/nTM/e53zr0U7POOmOexI51teyTwjnPu/mDdtcCqznbSx8/I2PdpBdG/D+Cf34ecc48Hn+dX4/+B37ebfYL/OzAKP3DT5Px5C45evMYiiZZq9ZiSXs7G/xFfG/TvDJZdgx+F/bSTPzrj8Ml6X6wL/uAB1AU/P49ZX4f/g90r5r/+/wOwJ1CE/915pRe7uBOfFNyG/yNzZ7C8dRRvZcy37Vn40dW+xhH7R7KW6OMd3W6/Pf2n5v/wI2xb0a78JIhpKj6J3gU/qpQP/KvdZrHHbR9Hh481xsbH45yrDZ6nkuDYM4Ef40fLsvDPyVvdPaAYuzvnFvdi+1ZbAlODpKFVDv65Go4fbe3re3g0bV+bT4NlvTE89nerh++b9r8n7fub9+C4/wZuMD/LxvbABufcS72Iu7P3bmefCaOBpc65SMyyT/HfPLTq7nG1HmNLYHS71zSbtoltq+H493r71yn2uJ2+b9tZCww3s5wOPg9HBevBJ+B/wH/b8hbwOHAL/p+QxTGfs22OTdvnsSM9+rxwzjnrYlahPn5Gtn+f3hqzrs3vgXMuYmZLafscd+Z3+G9YHgue9xudc1fRu9dYJKE0Ai4JYb7e9WT8CM0qM1sFXAJMMrNJ+A/2Lazjk/KW4ssKOlKL/3Bv1ZOkoDM1vdjX9fjR3+2cc4PwX3lbF9u39y9gmpmNBY4nmoAvxX/dP9w5VxbcBjnndk5AHCvxiUyrLXpyp2Ck8RP8iNj9HWxyJ758YZxzbjD+q9/2Mbl2cYyN6Y+jD8zPqnMffmRsM+dcGb4EojevS18tBf4X85qVOX/C41fxCVM9Hb+HXQfL2luBTxRabREs64/+vn97JBhlvgc/iv9Fejf63ZXOPhNWAONa63sDW9DuW5peHOOTdq9pqXPuyA62XYsfZW3/OvXluM/jPwNOiF1oZsX4b8eeCBY9B+yA//z4n3Pu3eCYR+GT83hr83tqPpMd2/nmcX+Ptfk9CI4/juhz3OnfguBbie8457YGjga+bWaH0rvXWCShlIBLohwHtAAT8F9pTsbXET+NL2d4Cf8Bf5WZFQcnGe0X3Pdm4FLzJ3CZmW1rZq0fxK8Dp5s/kW8GPSwD6cTrwJFmNtTMNgcu7mLbUnw9ZLWZ7Yiv6Yz1Ob6esENBKcF84B/4PwCLguUr8WUTvzezQWaWZf7kt84eV3dxdOUe4JtmNtbMhgCX9+K+XwYOcdFZQdrHVO6cqzezvWj7NXJncXzLzMaYn2bte72II1braPsaoDkYDZ/ex3311n+B7c3si2aWG9ymmNlOwWjs34E/BCd8ZZvZPsE/DGuACF28V/AlPFea2QjzJ939COjvdH39ed/01m34Guhj6H/crTr7THgR/4/0d4PXYBo+4bqrD8d4Cag0s++ZP4k228x2MbMp7TcMvmW7B/ilmZUGsXybPjzeoMTnp8CfzWxG8DjG4/9pX0bwT4xzrhY/ovx1ogn3c8BXSEwC/hCwq/kT2XOC43Y1SBHv99g9wFFmdqj5k3i/g/9H5blg/et08rfAzGYF7xELYmoJbj1+jUUSTQm4JMrZwD+cc58551a13oDr8KNjhv9DuS3wGf4PzSkAzrl/Ab/Ej6xW4U9SHBrs91vB/SqC/TzYjxj/D38i6BJ8Enx3F9teik8sq/A1g+23/Qnwz6CW8WQ6dif+pMY72y0/C59MvgusB+7Ff/Xclzi6chMwF/+YX6Xj0ewOOV9nvbCT1V8DfmZmVfhk8Z4exPEY8CbwGn7Uuhn/B7LHnHNVwDeD463HPy9zerMPojPXtN7+2ItjT8fXjq7Af43/G/w/BOBfp7eAl/Enuf0GyAqSqF8Czwbvlb072P0vgIX45+ct/GvV37mR+/O+6RXn3LP4fzJeDWqB47HPDj8TnHON+ER/Jn5U+q/AWc659/pwjBb8Z8tk/Dc+a/GJ/+BO7vINfPL/MfBMENvfe3vc4Ni/xY8YX41PGF/Ej9Ye6tqeO/A/fMnaSzH9UqKz68RNUNLyBeC3wDr8YMpCfBLckbi+x5xz7+PPJ/gz/rU4Gjg6eM2h678F2+FP7q7Gf8PwV+fPVentayySMK0zCYiIhCIYub7BObdltxvLgGB+isQ7nXM3hx2LxEdQ5rMMOMM591TY8YgMdBoBF5GkCr76PdLMcsxsDP4kygfCjkviI/g6f3cSOMouyWFmR5hZWVA+1VrT/ULIYYmkBSXgIpJshq95XY8vQVmEL12RAc7M/on/6v/ioExHBrZ98LPPtJaAHOf8VKIi0k8qQRERERERSSKNgIuIiIiIJFFoF+LJyspyhYWFYR1eRERERDJEbW2tc86lzMBzaAl4YWEhNTUdTSksIiIiIhI/ZpZS5y+kzH8CIiIiIiKZQAm4iIiIiEgSKQEXEREREUkiJeAiIiIiIkmkBFxEREREJImUgIuIiIiIJJEScBERERGRJFICLiIiIiKSRErARURERESSSAm4iIiIiEgSKQEXEREREUkiJeAiIiIiIknUbQJuZn83s9Vm9nYn683MrjWzxWb2ppntHv8wRURERETSQ04PtrkVuA64rZP1M4HtgttU4Prgp4ikiQdfW87v5r7Pioo6RpcVctkRO3DcbmPCDqvf+vK44vVcdLafRC/vzJUPvsXsF5fS4hzZZpw2dRy/OG7XpDwXYe0/XvoSZ1+e70THNFCOG6/3fDpItcecavGkMnPOdb+R2Xjgv865XTpY9zdgvnNudtB/H5jmnFvZ1T6Li4tdTU1Nn4IWkeR58LXlXHH/W9Q1tWxcVpibza9P2DWcD1bnoHY+VN0H9S9C83IgG3LHQcE+MOg0KNgDzLrcTV8eV7yei872c+IeY7jvleUJW95ZnFc++Ba3v/DZJsvP3HuLTpPCRL8vUu5914m+xNmX5zvRMXXIOWhpgZbmmFtLzLKWNssXLFrFLf/7kEhzMzm0kOsiFGXDGVPGMmXcIL9dpAUiEX9zzvedi1kW2bQds/6j1ZW8sHgtkUiEbCJkOUduFmw9rJCl62pwkQjZLoIBuVmwy+hBjB6UDzi/H0fHbYJ+bJvguMQub3+/YF2H+4i5QcfbxD7XXb0OHfSrG1oor2kgNo/LMhhanEdxXk7b/bffzyaH60Ms7barb2qhqr4J52Djp69BaX4OBTkxBRfDRsKV13Z+jAQxs1rnXHHSD9yJeCTg/wWucs49E/SfAL7nnFvY1T6VgIsMDPtd9STLK+o2WT6mrJBnLz8kucHUzIPV34GGN7vervAAGHk1FO7V6SZ9eVzxei4620+2GS0dfCbHa3lncW5zxcOd7uejXx/Zq8cQr/dFSr3vutCXOPvyfLcRaYGGen+rr4OGumi/oY5f3v8ajTU1FLsmilwTBa6ZfFoYkguzdhoGTY3Q3ARNTb7d1AjNjW37TUHfRXr9nIh0acQo+PU/kn7YVEvAe1KC0p2Ohpk6zOrN7ALgAoC8vLw4HFpEEm1FB8lFV8sTIlIHqy+Bir/1bPu6p+HTfWDY92H4T8CyN9mkL48rXs9FZ9t3lJTFc3m8jtvVvuL1vkiJ910P9CXOlkiEIpoYGqljWKSOwa6BQa6BwZEGeLgKaqtjbjVQWxX8rPHJdlNjlzH9oLMV9cBLPXtcIpJY8UjAlwHjYvpjgRUdbeicuxG4EfwIeByOLSIJNrqssMMRvtFlhckJoHktLDsa6l+ILrNiGHwWlB4HeTv6UbrGt6HyHqicDTQDEVj3C6h/HcbMhqySTeLv7eOK13PR2X4SPQLeWZxd7acziX5fhP6+66HYOAtcE5tHahjVUs1OBY0w9z6oqoDKCv+zagNUVbBoXTmFNHe8w/vnJynyHrIsyM2FrGzIzobsnA5++vabq2qobYYmy6KZLJqDn7l5uRy2y+jofSwLsoJba9usk2XZbdb/6cmPKK9rpgXDYbSY0UIWZkaTMyIYEbON64cU5/OrE3YNStIsGDI032/TDt7rsdtZVtvl3e1j4/2s3fbt98Gm64LdtHvyY5rtVxon/+15Pq+qx7W742aDCvjXhft0vNPYbgf77HxdJ/eLaR973TOsqGzYZPPNBxXwn2/sH91/libgg/gk4HOAi8zsLvzJlxu6q/8WkYHjsiN26LCe9LIjdkj8wZtXw2fToHFRdFnpSbDZtZAzqu22eeOhZBYM/xGsOt/XiQPU/BeWHgXjHoas6LePfXlc8XouOttPomvAO4vztKnjOqxJPm3quA627voxxOt9Eer7riORCFSuhzUrYfVKKF8N69dyX+MyNlSsZGRLNUNcfXT7SuBfHe+q3/9CmEFeARQUQH4h5BcEt0IoKOSzmgjPfFZNpcuh1nKpJQeXk8esPbZkt21GQm6eT6pz8yCnk3brLXvTb48683EXtefEqW5/yxHLuaGT353/dHLeQ7yOnWpOP3rvDp/vb8/a1Zd5JNm5R+3ZYTxfPmpXGDw06fGkum4TcDObDUwDhpvZMuDHQC6Ac+4G4GHgSGAxUAucm6hgRST5Wk/aSvqZ7S2VsHR6TPJtMPKPMOQbXY/O5G0L456ANd+H8t/4ZXULYNnxMO4hsFygb48rXs9FV/vZc8uhCV3ekdYT/3ozK0ei3xehvO+cgw3lsPIzWPEZrF4Ba1bB2pX+ZwelH5sHt95qysphDQWss0I2ZBVQNnwIu2w3BgqLobgECkugqBiKS2OWFUNBkU+MuxhF3AIoem05f2n33O2W4N/ZZLxmffndSVehfTYPkHhSXY9OwkwEnYQpIp1yLbDsWKh5KFiQBaPvhEGn9G4/634Pay6N9ssugM17WEcu6W3Delj6kU+0WxPulZ/5uuu+ys6GsmEwZDiUDYeyoVBaBqWD/c9BrT/L/Gh1NzP1iEj8pONJmCIi8bXuVzHJN7D5zb1PvgGGfQciVbDup75fcSMUHgiDz4hPnJL6nIP1a+HTxfDZ4ujPinW931dxKYzYHEaMhuGb+UR7yIjg53CfaKu+VUR6QCPgIpJa6l6ET/cDgjrCod+DkVf1fX/OwcozgpMzgaxSGP865G3d30glFTU3w2cfwUfvwuJ34KNFvUu2C4pg9BYwagvYbAyMHOXraUeMgqKS7u8vIikp1UbAlYCLSOqIVMMnu0HTYt8v3A+2mA/Wzy/rWqpgye7R/RZMhS2f3lgPLgNYczN88h68+xq8/yYs+QAaN52JYRN5+TB2K38bvaVPuEeN8yPZKg0RSTuploCrBEVEUsfqS6NJclYpjPq//iffANmlvob8032BZn8FzXW/g+Hf7/++Jbmcg5VL4d1XYdFr8P5bUF/b9X0KimDLbWGLbWGLbXx787F+ijsRkRBoBFxEUkPdS/Dp1Gh/1G0w+IvxPca638Ka7/m2FcLW70HuFvE9hsRfczN8+Da8/ry/rVvd9fbDN4NtJsC2O8O2E2DMlkq2RTJcqo2AKwEXkfC5CHy6N9S/7PvFs2DsnPiXArhmWLInNLzh+6Unwph743sMiY+mRnj7FXjlaXjzpa5nJxk2EibsDjvtBtvv4mciERGJoQQ8oARcRDaq+Dus+rJvWz5s9Q7kbZOYY9U+C5/tH+2PmwvF0xNzLOmdSIsvKXnxKXj12c6T7sIi2HEyTNjNJ94jR6tuW0S6lGoJuGrARSRcLRtgzeXR/tBLE5d8AxTtB4POgsrbfP/zb8FWb8Wn1lz65vPl8PSj8PwT/kI4HRkyHCbvA7vtA9vvCjk6gVZEBi79xRGRcJVfAy1rfDtnHAy7IvHHHPkbqH7AzxHe+B5U3gGDz078cSWqqRFeew4WPALvvdHxNsM3g72mwR77+xMoNcotIgliZjOAPwHZwM3Ouavarbdg/ZH4K7+f45x7NWZ9NrAQWO6cm9Xd8ZSAi0h4Wsph/TXR/ohfQlYSviHM2dyPtK/9se+v/SkMOl3TEiZDxTp4co5PvKsrN11fWgZTDoSp02DrnZR0i0jCBcnzX4DDgWXAy2Y2xzn3bsxmM4HtgttU4PrgZ6tvAYuAQT05phJwEQlP+R8gEiRheTv4JDhZhlwM5X+CSDk0fQIb/uEvVS+JsfRjePx+eHE+tDS3XWdZMHEKHDATdp3iL+kuIpI8ewGLnXMfA5jZXcCxQGwCfixwm/MnT75gZmVmNso5t9LMxgJHAb8Evt2TA4aWgA8dOpT58+eHdXgRCVmubWDqiD+QE1y5+93VJ7P6s6eTGsO44pPYpvRGAOqXX8lLr48nQl5SY0h7DXVQUe7n6s4eDPseG12Xkwslg/wtOwcq6uHp5L4HRCRj5JjZwpj+jc65G4P2GGBpzLpltB3d7mybMcBK4I/Ad4HSHgfT0w3jrby8nGnTpoV1eBEJ2+oroLzOt/N2ZsLUnzDBspIbQ2QKfPQgtKymIHsNB076CIZ8NbkxpKsP34Z/3w7vvb7pum0mwPQT/AmVmp9bRJKj2Tm3ZyfrOqp1az9NYIfbmNksYLVz7hUzm9bTYFSCIiLJ11IJFX+N9of/xJchJFtWMQy7HFYH3xiW/8GXoZiSwj776F2feL/7atvllgV77AeHnwDb7BRObCIiHVsGjIvpjwVW9HCbk4BjzOxIoAAYZGa3O+fO7OqAmgdcRJKv/Jpo0pu3A2z1bjgJOEBLFXy0BUQqfH/M/VB6fDixDGRrVsK9t8Arz7RdnpUF+x4GR50GI0aFE5uIZLyu5gE3sxzgA+BQYDnwMnC6c+6dmG2OAi7Cz4IyFbjWObdXu/1MAy7VLCgiknpcsz/5sdWQS8JLvgGyS6HsQigPZpwq/50S8N6oq4GH7oJ5D0JzU3S5ZcE+h8Cs0/2FckREUpRzrtnMLgLm4qch/Ltz7h0zuzBYfwPwMD75XoyfhvDc/hxTI+AiklyVd8OKU307ezhs8xlkFYYbU/NKWLwlECSQWzwLRfuGGlLKi0TgmbnwwD+hqqLtur2mwTFnwuZjw4hMRGQTuhKmiGQu56D899F+2dfCT74BckbB4DP9VIQA5VdD0f3hxpTKln0Ct/0JPn6v7fKtd4RTL/Q/RUSkUxoBF5HkqXsePg1Gli3fj37njAw3plYN78InOwedLNhmCeSO6+oemaexAf57J8y9F1paosuHjoATv+RHvnXhHBFJQRoBF5HMtf6GaHvQ6amTfAPkT4CiQ6D2SSACFTfBiJ+FHVXq+GgR/P1q+Hx5dFl2Dsw82d/yC8KLTURkgNEIuIgkR0s5LB4NrsH3t3wZCjubkjUklffCii/4ds4o2OZTXZ6+uQnm3AGP3AMuEl2+3c7wxW/C6C3Di01EpIc0Ai4imWnDP6PJd8EeqZd8A5Qe6xPv5pX+VjUHBp0YdlThWb4Ebv6tv4x8q8Ii+ML5sP8RfopBERHpNX16ikjiOQcVMeUnZReGF0tXLBcGnxftV1wfXixhcs7PcPLLb7VNvnecBD+5Hg6cqeRbRKQfVIIiIolX8xQsPcS3swbBtiv8VShTUdNS+Gg8EJRbbP0+5G0fZkTJVV8Ht/8ZXngyuiw3z59kecgxSrxFZEBKtRIUfZKKSOJtuCnaHnRW6ibf4Gc+KTk62q/4R3ixJNuKz+AX32ibfI/eAq68Fg47Tsm3iEic6NNURBKrpQKqHoj2y84PLZQeG/ylaLvyNnAtnW+bLl5/AX51MaxaFl2233T4wbUwZnxYUYmIpCWdhCkiiVV1D7h6387fDQomhhtPT5TMhOyR0LIamldAzeNQMiPsqBLDOX8p+X/f5tsAeflw5jdg38PCjU1EJE1pBFxEEmvDrdH24HPCiqJ3LBcGnRntxz6GdNLUCH/7NTz4z2jyPWwkXHGNkm8RkQRSAi4iidPwvr/6JQC5/uI7A0XsPwvVD0LL+rAiSYzqSvj95bBwQXTZDhN9vfe4rcOLS0QkAygBF5HEqbwt2i6ZBTnDw4ultwp29fOVg5+/vPKucOOJpzWr4Kpvw+J3o8sOngWX/ApKy0ILS0QkUygBF5HEcC2wISYBH3x2eLH01eBzo+0N/wwvjnha8uGmJ1uefAGccRHk6LQgEZFkUAIuIolR+zQ0B0le9nB/YuNAM+hUNp6rXv8iNH7c5eYp78O34ervQVWF7+fkwoXfh+knhBqWiEimUQIuIolRNTvaLj0ZLC+8WPoqexgUHxHtD+QylHdehWt+APW1vl9UAt/5Nex5YLhxiYhkICXgIhJ/rhEq/xXtD6STL9uLjb3yjuhsIQPJa8/Dn38MjQ2+P2gIfPd3sN0u4cYlIpKhlICLSPzVPAaRYNaQnC2gcJ9w4+mP0mPACn278V1oeCvceHrrlWfg+p9Dc5PvDx3hk++xW4Ubl4hIBlMCLiLxVxlTfjLoVLAB/FGTVQIlx0b7sY8t1b3+Atz4a4hEfH/EKPje1bD52HDjEhHJcAP4r6KIpKRIDVQ9GO0POi20UOKmTRnKbHCR8GLpqbcXwg2/hJYW399sjE++h20WblwiIqIEXETirPo/4IIT/fJ2gvxJ4cYTDyVHQNYQ327+1M+IksoWvQ5/+Vm07GTEKLj0N1A2LNSwRETEUwIuIvHV5uTLU8EsvFjixfKg9Phov+q+8GLpzifvw3U/8ZeZBxg6Ei69CoYMoIsgiYikOSXgIhI/kRqoeSTaL/1CeLHEW+lJ0XbVfak5G8rny+FPP4KGet8fMtwn3yo7ERFJKUrARSR+qh8FV+fbeTtB/k7hxhNPxYdC1mDfbloCDa+GGs4mNpT7eb6rN/h+cam/tPzI0eHGJSIim1ACLiLxUx1TmlF6YnhxJILlQckx0X4qlaHU1cAffwhrV/l+Xj5882cweotw4xIRkQ4pAReR+Ig0QPV/o/10S8Ch7WOqvDc1ylBaWuD6X8LSj3w/K8tfXn6bNPr2QUQkzSgBF5H4qH0cIlW+nbt1esx+0l7xdD8vOEDTh9DwdrjxANx1A7wbUw5z1sUwcWpo4YiISPeUgItIfFS1Kz9Jh9lP2ssqhOKjov2wy1CenANP/SfaP/oM2H96ePGIiEiPKAEXkf5zTVD172g/HctPWrWZDeXe8OJ45xU/+t1qykFwzJnhxSMiIj2mBFxE+q92PkTW+3bOWCiYEmo4CVUyE6zQtxvfgYb3kx/DqqVww6+il5gfvz2c++30/NZBRCQNKQEXkf5rU35yAlgaf7RkFUPxzGg/2WUo9XXwl5/7mU/Az/V90Y/9zCciIjIgpPFfSRFJCtcCVQ9E++lcftIq9jEmMwF3Dm69BlZ+5vu5eT751iXmRUQGFCXgItI/dS9Ay2rfzh4JhfuFG08ylMwCcn274VVoWpac4z7+ACxcEO1/8Zuw5XbJObaIiMSNEnAR6Z/qOdF2ybFg2eHFkizZg6BoWrRf/Z9ON42bD96Ce2+O9qfNgn0PS/xxRUQk7pSAi0j/tEnAjw4vjmQrjbkqZuxzkAgb1rc96XLrHeGUCxJ7TBERSRgl4CLSd40fQuN7vm2FUHxouPEkU+w/G7VPQktVYo4TicDffweVwSwzpYPhwh/4+m8RERmQlICLSN/Fll4UHw5ZReHFkmy5W0av9ukaoeaxxBznsfvgnZgrXZ73PRg6IjHHEhGRpFACLiJ9V5Wh5SetShJchvLx+/DArdH+jJNh593jfxwREUkqJeAi0jct5VD3TLRfMiu8WMISm4DXPASuOX77rquBm66Clhbf33pHOO6s+O1fRERCowRcRPqm+hEgSA4LpkLO5qGGE4qC3SFntG+3rIO65+O37zv+AmtW+nZhEZx/OeTkxG//IiISGiXgItI3bWY/Oabz7dKZZbUtvYlXGcrCp+GFJ6P9sy6GERn4D46ISJpSAi4ivecaoeaRaD8T679bxbsOfEM53P7naH+fw2DKgf3fr4iIpAwl4CLSe7ULIBJMu5c7HvJ3CTWcUBUdAhbM/tL4ATS83/d9OQe3/QmqK31/6Ag47av9j1FERFKKEnAR6b325Sdm4cUStqwCKD4i2u/PKPizj8EbL0b7534bior7vj8REUlJPUrAzWyGmb1vZovN7PIO1g82s/+Y2Rtm9o6ZnRv/UEUkJTgH1f+N9jO5/KRV7HMQW5rTG+s+h9k3RPuHHAM77da/uEREJCV1m4CbWTbwF2AmMAE4zcwmtNvs68C7zrlJwDTg92amy7SJpKPG96HpE9/OKoEi1SdTMiParn0aWip7d//W0pOGOt/fbAyc+KX4xSciIimlJyPgewGLnXMfO+cagbuAY9tt44BSMzOgBCgH4jghroikjJqHo+2iw0H/a0POKMhvvUBOM9TO6939n38ierVLM/jSdyC/IK4hiohI6uhJAj4GWBrTXxYsi3UdsBOwAngL+JZzLtJ+R2Z2gZktNLOFzc3Kz0UGpOrY2U9mhhdHqol9Lqp7UYZSWQF3/y3aP+QY2Kb9l4wiIpJOepKAd3R2lWvXPwJ4HRgNTAauM7NBm9zJuRudc3s65/bM0QUlRAaeSDXULYj2i5WAb1R8ZLRd87AvK+mJ2ddDTTCjzLCRcPw5cQ9NRERSS08S8GXAuJj+WPxId6xzgfudtxj4BNgxPiGKSMqoedLPAQ6Qvyvkjg03nlRSOBWyhvp28wpoeLP7+7z+Arz8v2j/i9+EgsLExCciIimjJwn4y8B2ZrZVcGLlqUD7ebY+Aw4FMLPNgB2Aj+MZqIikgNj679gRXwHLhpKY6Qhjn6uO1NfCHddF+/scBrvsmZjYREQkpXSbgDvnmoGLgLnAIuAe59w7ZnahmV0YbPZzYF8zewt4Aviec25tooIWkRA4167+Wwn4JmJLcqq7ScDn3AHrg4/J0sFwygWJi0tERFKKuZ7WKcZZcXGxq6mpCeXYItIHDe/AJ8EVL7MGwXZrwXLDjSnVNK+GxZvjT5PJ8s9R9pBNt1v2Cfzs6xAJzlX/8mWwz6HJjFREJKOYWa1zLmWubKYrYYpIz8SO6BZPV/LdkZyRUDAl6ESg5vFNt3HOl560Jt/b7wp7H5K0EEVEJHxKwEWkZ2Kv8Kj6787FluZ0VIby3Dz48B3fzs6GMy/yc3+LiEjGUAIuIt1rqfRXeGwVe+VHaavNdISPQOwlEaqr4N6bo/3DT4TRWyYvNhERSQlKwEWke7VPsPHitvm7+Ss/SscK9oDsEb7dshrqX42u+/c/oWqDbw8dAUefnvz4REQkdErARaR7saUUuvpl1ywLimO+IWidjnD5Epgf8zyeeqEuNy8ikqGUgItI15xT/XdvtakDf8Q/h3f9LVqOMmF32G3fcGITEZHQKQEXka41vAXNy307a4i/4qN0rXg6Gz9e61+EN+bCotd837L8nN868VJEJGMpAReRrrUZ/Z4OlhNeLANF9lAo3CfoOFh4VXTdQUfCmPFhRCUiIilCCbiIdK3N1S9V/91jsXXgQz70PwuL4dgvhhOPiIikDCXgItK5liqoey7aLz4ivFgGmtjnasu1gIOjz/CXnRcRkYymBFxEOlf3P6DJt/MnQc7moYYzoBTsDs1Fvl3SCDsWwiFHhxuTiIikBCXgItK56rnRdvH08OIYiFYshcUxo90zxkFObnjxiIhIylACLiKdq30s2lb5Se/86xZYMizaH7QovFhERCSlaDoDEelY4xJo/MC3rRAK9ws1nAHl/Tfh7ZehJCYBr3saIrWQVRReXCIikhI0Ai4iHYsd/S46CLJ01cYecQ7u+7tvVxdA7WbB8kao/V94cYmISMpQAi4iHauJLT9R/XePvfYcfPyeb+fkwrDjoutq5nZ4FxERySxKwEVkU64Zap6I9lX/3TMtLXD/rdH+wbNg2PHRfuw/NSIikrFUAy4im6pfCJEK384ZA3k7hRrOgPHc47BqqW8XFsGRp0JRHlg+uAZoXARNSyF3XLhxiohIqDQCLiKbqmk3/aBZeLEMFI0NMOf2aP+Ik/xFd7IKofDA6HKVoYiIZDwl4CKyqRpNP9hrT86B9Wt9e9AQOPyE6LrY51BlKCIiGU8JuIi01VIBdS8GHYOiQ8OMZmCorYaH7472jz4d8mNmjSmJTcDngWtJXmwiItItM5thZu+b2WIzu7yD9WZm1wbr3zSz3YPlBWb2kpm9YWbvmNlPe3I8JeAi0lbtk0CQIBbsATnDQw1nQJh7r0/CAUaMggNmtl2ftzPkjPbtyHpfYy8iIinBzLKBvwAzgQnAaWY2od1mM4HtgtsFwPXB8gbgEOfcJGAyMMPM9u7umErARaQtlZ/0TnUlzPt3tH/cWZDT7vx2s7ZTOaoOXEQklewFLHbOfeycawTuAo5tt82xwG3OewEoM7NRQT8YgSE3uLnuDhjaLChDhw5l/vz5YR1eRDrkmDr83xQGnwyvfTCCDe/MDzWilFexDvY60rdz86DOoIPPtpEF45hQ5tsbVtzDa28fuMk2IiKSMDlmFvv1443OuRuD9hhgacy6ZcDUdvfvaJsxwMpgBP0VYFvgL865F+lGaAl4eXk506ZNC+vwItKRxg/h41W+nVXCbvt+FSwv3JhSWVUFXH4ONNT7/le+D1M6Saybd4HFvwAcg/PeY9oBkyG7LClhiogIzc65PTtZ19FUX+1HsTvdxjnXAkw2szLgATPbxTn3dlfBqARFRKJiy0+KDlby3Z2590WT7zHjYY/9O982Z7ivqQegJai1FxGRFLAMiL1Aw1hgRW+3cc5VAPOBGd0dUAm4iESp/rvnKiv81IOtjjkTsrr5SG0zHaHqwEVEUsTLwHZmtpWZ5QGnAnPabTMHOCuYDWVvYINzbqWZjQhGvjGzQuAw4L3uDqgEXEQ819R2VDb2pEHZ1Nx/+YvvAIzdCnbbt/v7xD6n1XPBdXuejoiIJJhzrhm4CJgLLALucc69Y2YXmtmFwWYPAx8Di4GbgK8Fy0cBT5nZm/hE/nHn3H+7O6a5kP4AFBcXu5qamlCOLSIdqF0Anx3k27njYeuPdQXMzmxYD1ecE03Av/ZD2H2/7u/nmuDDoRAJTpjf+n3I2z5hYYqIiGdmtc654rDjaKURcBHx2pefKPnu3KMxo9/jtunZ6DeA5ULRIdG+roopIpKRlICLiNcmAVf5Sac2lMP/Hor2jz2zd/+sxNaBV6sOXEQkEykBFxFoXhtzdcbstqO00tZj90VHv7fcDiZ1e8GztmIT8NqnwDXGLzYRERkQlICLCNQ+wcYpTwunan7qzlRXwvyY0e+jT+99qU7eNpC7tW+7Gqh9Ln7xiYjIgKAEXERUftJTT/y77bzfE9tfKK2HNB2hiEhGUwIukumca5sEKgHvWH2tT8BbHXlK9/N+d6ZNAq4TMUVEMo0ScJFM17gImpf7dlYZFEwJNZyU9dRDUBtMHzhydOeXnO+JooOBHN9ueBWaV/c7PBERGTiUgItkujblJ4eC5YQXS6pqbIDH74/2Z54MWdl931/2ICjcJ9qvmdf3fYmIyICjBFwk06n8pHvPzIXK9b49ZDjsc2j/9xn7XKsMRUQkoygBF8lkkXqo/V+0X6QEfBPNzfDovdH+ESdBTm7/9xubgNc+psvSi4hkECXgIpms7llwdb6dtz3kjQ81nJT04pNQHtRolw6GA2bEZ78Fe0DWUN9uXgkNb8dnvyIikvKUgItkstjSB41+byrSAg/fHe0ffjzkF8Rn35YNxYdF+ypDERHJGErARTKZ6r+79trz8HkwQ0xhMUw7Or77b1MHrvnARUQyhRJwkUzVvAoa3gg6uVB8cKjhpBznYG5M7ffBR0NRcXyPEZuA1y2ASF189y8iIilJCbhIpoqd+q5wX8gqCS+WVPThO/Dxe76dkwuHHhP/Y+SOg7ydfNs1QO2C+B9DRERSjhJwkUyl8pOuzf1XtL3PoTB4aGKOo6tiiohkHCXgIpnIRaDm8Wi/5IjOt81EKz6DN16M9qefmLhjtZ+OUERE0p4ScJFM1PAWtHzu29nDIH+3cONJNY/dF21P3htGjUvcsYoOBMvz7Ya3oWl54o4lIiIpQQm4SCaKLT8pOhxMHwUbVayDF56M9o84KbHHyyqGwgOi/dhvJkREJC3pr65IJoqtNVb9d1tP/geam3x76x1h250Tf0yVoYiIZBQl4CKZJlIDdU9H+0rAo+rrYP5/o/0jTgKzxB+3zXzgj/safRERSVtKwEUyTe3/wDX6dv4ukDsm3HhSyTNzobbat0eOht32Sc5x8ydC9ma+3bIWGl5LznFFRCQUSsBFMk2b+m+Nfm/U0gKP3x/tH34CZGUn59iWBcWHR/vVuiqmiEg6UwIukmli6781/WDUwgWwbrVvlwyG/Q7vevt403zgIiIZQwm4SCZp+gwag6s7WkHb2TcyWfvLzh9yNOTlJzeG4sOi7brnoKUquccXEZGkUQIukknalJ8cCFmF4cWSSt57Az77yLfz8uHgo5MfQ87mkD8p6DRB7fzkxyAiIkmhBFwkk7S5/LzKTzZ6/IFoe9/DoXRwOHGoDEVEJCMoARfJFK4Zap6I9pWAe58vhzdjLjt/2HGhhaL5wEVEMoMScJFMUf8yRCp8O2cM5E0INZyU8cS/o+2Je8HmY8OLpXA/sKAsqPEDaFwSXiwiIpIwSsBFMkXs1HbF05NzgZlUV1sNz8aMNIc5+g2QVQBF06J9jYKLiKQlJeAimUL135t6ei401Pv2mPGw026hhgO0uyqmEnARkXTUowTczGaY2ftmttjMLu9km2lm9rqZvWNm/4tvmCLSLy3rof6loGNtp7zLVC0tbctPDjsuNb4VaJOAz/O1+yIikla6TcDNLBv4CzATmACcZmYT2m1TBvwVOMY5tzPwhfiHKiJ9VvMEEPHtgimQPSzUcFLCa89BecyFd6YeHG48rfJ2gpygDj2ywdfui4hIWunJCPhewGLn3MfOuUbgLuDYdtucDtzvnPsMwDm3Or5hiki/1LSr/xaYFzP14LQjk3/hnc6YtX2NdFl6EZG005MEfAywNKa/LFgWa3tgiJnNN7NXzOysjnZkZheY2UIzW9jcrK9VRZLCOdV/t/fJ+7D4Xd/OzoFps8KNp73Y10gnYoqIpJ2cHmzTUVGk62A/ewCHAoXA82b2gnPugzZ3cu5G4EaA4uLi9vsQkURofA+ag/+hswZB4dRw40kF8x6Mtvc6CMpSrCSn+FD8R6+DuhehpQKyy8KNSURE4qYnI+DLgHEx/bHAig62edQ5V+OcWwssACYhIuGLnUmj6BCw3PBiSQXr18LCBdH+oceFFkqnsodBwZ5BJ9L2AkoiIjLg9SQBfxnYzsy2MrM84FRgTrtt/g0cYGY5ZlYETAUWxTdUEekTlZ+09dR//QwoANvtAuO3CzeezqgMRUQkbXWbgDvnmoGLgLn4pPoe59w7ZnahmV0YbLMIeBR4E3gJuNk593biwhaRHonUQ+38aD/TE/CGeljwcLR/+PHhxdKd9idiOlXtiYiki57UgOOcexh4uN2yG9r1fwf8Ln6hiUi/1T0Drs63c7eDvK3CjSdsLz4F1ZW+PXwzmLx3uPF0pXBvyCqFSBU0fwpNH0Le9mFHJSIicaArYYqks9j670yfftC5thfeOeQYyMoOL57uWK6v2W+lq2KKiKQNJeAi6Uz131Efvg3Ll/h2Xj7sPwCeD12WXkQkLSkBF0lXzSuh4c2gkwvFKXKlx7DEjn7vcygUlYQXS0+1ORHzKXCN4cUiIiJxowRcJF21mX5wP8gaAAlnopSv8Zeeb3XIMeHF0ht520Du1r4dqYa658ONR0RE4kIJuEi6Uv131P8egkjEt3eYCGPGhxpOr7QpQ9Fl6UVE0oEScJF05CJQ83i0n8n1302NsOCRaH+gjH63in3tqh8NLw4REYkbJeAi6aj+FWhZ49vZIyB/cqjhhGrhAqja4NtDhsPkfcKNp7eKDgGCq5c2vOZr+0VEZEBTAi6Sjmpipu0vngGWwb/qT/4n2p42C7JTeOrBjmQPgqL9o32NgouIDHgZ/FdZJI1Vx5RclBwZXhxh+/h9+OR9387JhQNnhBtPXxXHvIY1D3e+nYiIDAhKwEXSTfNaqH8p6GRl9gmYT82JtqccBKVloYXSLyUzo+2ax8E1hReLiIj0mxJwkXRTMxdwvl04FbKHhhpOaCor4OUF0f6hA+zky1h5EyBnC9+ObNB0hCIiA5wScJF0UxNTflKcweUnCx6B5mCkeOsdYfz24cbTH2ZtS4liS4xERGTAUQIukk5cpO1c0bGlC5mkpcXP/d1qoE092JHi2DIU1YGLiAxkSsBF0kn9QmhZ69vZm0H+buHGE5bXnoP1wfNQWgZ77N/l5gNC8SFgeb7d8CY0LQs3HhER6TMl4CLppFrTDwLwZMzJlwcdCbl54cUSL1klUHhgtF+j6QhFRAaqDP3rLJKmYuu/M7X8ZOnH8MFbvp2d7RPwdNGmDlxlKCIiA5UScJF00bwG6l8OOhk8/eBTMRfe2W0/f/XLdBFbB147D1xjeLGIiEifKQEXSRdtph/cB7KHhBpOKGqq4IUno/10OPkyVt4OkLuVb0eqoPbZcOMREZE+UQIuki7aTD+YoeUnzz4GjQ2+PW5r2G7ncOOJN7N2s6FoOkIRkYFICbhIOnAt7aYfTKO6556KtMBT/432Dz7aJ6zpRnXgIiIDXk7YAcRqampi2bJl1NfXhx2KpJiCggLGjh1Lbm5u2KGkpvqXoWWdb2dvDvmTQw0nFG8vhDUrfbuoBKYeHG48iVJ0MFg+uAZofAeaPoPcLcKOSiStKB8ZuAZKvpBSCfiyZcsoLS1l/PjxWDqOXEmfOOdYt24dy5YtY6uttgo7nNQUe2XEkhnpOfLbnSdiph484AjILwgvlkTKKoKiadFvPKofgSFfCTUkkXSjfGRgGkj5QkqVoNTX1zNs2DC92aUNM2PYsGEaiehKpl9+ftUyeOcV3zaDabPCjSfRYl9j1YGLxJ3ykYFpIOULKZWAA3qzS4f0vuhC86qY6QezofjwUMMJRezUgxP3ghGjwoslGWLneK+ZB5GG8GIRSVP6uzMwDZTXLeUS8FTwwAMPYGa89957YYeyifHjx7N27dout/nVr37Vpr/vvvvG5djZ2dlMnjx54+2qq64C4Omnn2bnnXdm8uTJ1NXVcdlll7Hzzjtz2WWX9foY7WOXHqh+KNouOgCyy0ILJRT1tfDc49F+uk092JG87SB3W992NVD3dLjxiEjcmRlf/OIXN/abm5sZMWIEs2bF7xu+I488koqKij7d9yc/+QlXX311h8vHjBnTJl9oPcZpp53GxIkTueaaa3jvvfeYPHkyu+22Gx999FGvjj1//nyee+65PsWdKpSAd2D27Nnsv//+3HXXXWGH0iftk9h4vUkLCwt5/fXXN94uv/xyAO644w4uvfRSXn/9dQoLC/nb3/7Gq6++yu9+97t+xy49UB0z80dxmpdedOT5J6Cu1rc3Hws77RZuPMkSOwperTIUkXRTXFzM22+/TV1dHQCPP/44Y8aM6dU+mpubu1z/8MMPU1ZW1tcQO3XJJZe0yRfKyspYtWoVzz33HG+++SaXXHIJDz74IMceeyyvvfYa22yzTa/2rwQ8DVVXV/Pss89yyy23tEnAW1pauPTSS9l1112ZOHEif/7znwF4+eWX2XfffZk0aRJ77bUXVVVV3HrrrVx00UUb7ztr1izmz58PQElJCd/73vfYY489OOyww3jppZeYNm0aW2+9NXPm+JPIurp/rOOOO4499tiDnXfemRtvvBGAyy+/nLq6OiZPnswZZ5yx8ZgAp5xyCg8/HJ227JxzzuG+++6jpaWFyy67jClTpjBx4kT+9re/9fj5uvnmm7nnnnv42c9+xhlnnMExxxxDTU0NU6dO5e6772bNmjWceOKJTJkyhSlTpvDss89ufJ7PPffcjc/nfffd12Hs0o1IPdQ8Fu2XHh1eLGFwDp6MKT+ZNguyMuRjrU0d+H87305EBqyZM2fy0EP+W87Zs2dz2mmnbVz30ksvse+++7Lbbrux77778v777wM+h/jCF77A0UcfzfTp06mtreXkk09m4sSJnHLKKUydOpWFCxcC0W/VlyxZwk477cT555/PzjvvzPTp0zcm/jfddBNTpkxh0qRJnHjiidTW1vbpsUyfPp3Vq1czefJkfvrTn/LHP/6Rm2++mYMP9jNW3X777ey1115MnjyZr3zlK7S0tADw6KOPsvvuuzNp0iQOPfRQlixZwg033MA111zD5MmTefrpgfkNYErNgtLGeTMSt++bH+101YMPPsiMGTPYfvvtGTp0KK+++iq77747N954I5988gmvvfYaOTk5lJeX09jYyCmnnMLdd9/NlClTqKyspLCwsMtD19TUMG3aNH7zm99w/PHHc+WVV/L444/z7rvvcvbZZ3PMMT3/+vzvf/87Q4cOpa6ujilTpnDiiSdy1VVXcd111/H6669vsv2pp57K3XffzZFHHkljYyNPPPEE119/PbfccguDBw/m5ZdfpqGhgf3224/p06dvcgZxa3Lc6oorruC8887jmWeeYdasWZx00kmAT/hbj3/66adzySWXsP/++/PZZ59xxBFHsGjRIn7+858zePBg3nrrLQDWr1/PiSee2Gns0ona+eCCD8Pc7SBv+1DDSbr33oCVn/l2fiHsl0H170XTwIr869/4gb9l2usvkgzvJbCmeEfX5epTTz2Vn/3sZ8yaNYs333yTL33pSxsTzh133JEFCxaQk5PDvHnz+P73v899990HwPPPP8+bb77J0KFDufrqqxkyZAhvvvkmb7/9dpu/47E+/PBDZs+ezU033cTJJ5/Mfffdx5lnnskJJ5zA+eefD8CVV17JLbfcwje+8Y0u477mmmu4/fbbARgyZAhPPfUUc+bMYdasWRv/xjvnKCkp4dJLL2XRokXcfffdPPvss+Tm5vK1r32NO+64g5kzZ3L++eezYMECttpqK8rLyxk6dCgXXnjhxvsOVKmbgIdk9uzZXHzxxYB/48+ePZvdd9+defPmceGFF5KT45+yoUOH8tZbbzFq1CimTJkCwKBBg7rdf15eHjNm+H8udt11V/Lz88nNzWXXXXdlyZIlvYr12muv5YEHHgBg6dKlfPjhhwwbNqzT7WfOnMk3v/lNGhoaePTRRznwwAMpLCzkscce48033+Tee+8FYMOGDXz44YebJOCtJSi9MW/ePN59992N/crKSqqqqpg3b16bbxiGDMnAy6bHQ3XM6G9Jho1+AzwZM/XgvodCYXF4sSRbVgEUT4fqB32/6j8w7DuhhiQi8TVx4kSWLFnC7NmzOfLItjNcbdiwgbPPPpsPP/wQM6OpqWnjusMPP5yhQ4cC8Mwzz/Ctb30LgF122YWJEyd2eKytttpqY3K+xx57bMxJ3n77ba688koqKiqorq7miCOO6DbuSy65pFfJ8RNPPMErr7yyMZ+qq6tj5MiRvPDCCxx44IEb85HWx5QOlIDHWLduHU8++SRvv/02ZkZLSwtmxm9/+1ucc5ucWdvRMoCcnBwikcjGfux0OLm5uRvvk5WVRX5+/sZ2a61WV/dvNX/+fObNm8fzzz9PUVER06ZN63banYKCAqZNm8bcuXO5++67N36V5Zzjz3/+c49+qXorEonw/PPPb/LNQGfPnfSCc23rv0syrP573Wp4/YVo/+AM/Aek5OhoAl49Rwm4SBo65phjuPTSS5k/fz7r1q3buPyHP/whBx98MA888ABLlixh2rRpG9cVF0cHI5zrepS9VWs+An7ShdYSlHPOOYcHH3yQSZMmceutt3ZYEttfzjnOPvtsfv3rX7dZPmfOnLTNFVI3Ae+iTCRR7r33Xs4666w2NdAHHXQQzzzzDNOnT+eGG25g2rRpG0tQdtxxR1asWMHLL7/MlClTqKqqorCwkPHjx/PXv/6VSCTC8uXLeemll3oVR0/uv2HDBoYMGUJRURHvvfceL7wQTURyc3Npamrq8CpQp556KjfffDMLFy7k1ltvBeCII47g+uuv55BDDiE3N5cPPviAMWPGtPkF7qvp06dz3XXXbZwR5fXXX2fy5Mkbl//xj38EfAnKkCFDuoxd2ml4C5qD8ouswVC0f7jxJNv/HgYX/KO642QYvWWo4YSi5CjAAAd1z0JLOWSnzwiRSEropkwk0b70pS8xePBgdt111zbJ74YNGzaelNn697wj+++/P/fccw8HH3ww77777sbSz56qqqpi1KhRNDU1cccdd/T6RNCeOPTQQzn22GO55JJLGDlyJOXl5VRVVbHPPvvw9a9/nU8++aRNCUppaSmVlZVxjyOZMuRspZ6ZPXs2xx9/fJtlJ554InfeeSfnnXceW2yxBRMnTmTSpEnceeed5OXlcffdd/ONb3yDSZMmcfjhh1NfX89+++3HVlttxa677sqll17K7rvv3qs4enL/GTNm0NzczMSJE/nhD3/I3nvvvXHdBRdcwMSJEzs8kXH69OksWLCAww47jLy8PADOO+88JkyYwO67784uu+zCV77ylQ7PnG6tAW+9tc6C0pVrr72WhQsXMnHiRCZMmMANN9wA+Dqy9evXs8suuzBp0iSeeuqpbmOXdtrMfjIDLIP+aWlqhAUxM38ckoGj3wA5m0HB1KDTotlQRNLQ2LFjN5aQxPrud7/LFVdcwX777bfxhMWOfO1rX2PNmjVMnDiR3/zmN0ycOJHBgwf3+Pg///nPmTp1Kocffjg77rhjj+7TeoJk6627EtsJEybwi1/8gunTpzNx4kQOP/xwVq5cyYgRI7jxxhs54YQTmDRpEqeccgoARx99NA888MCAPgnTevrVRLwVFxe7mpqaNssWLVrETjvtFEo8kvr0/mhnyT5QH3zzMep2GJxB/7Q8Pw9uCeafHToCfn0rZGeHGlJo1v0a1nzft0tPhjF3hxuPSBpIp783LS0tNDU1UVBQwEcffcShhx7KBx98sHEQLh119PqZWa1zLmVOFErdEhQR6Vzzaqh/MehkQUkCZw1KRe2nHszU5Bt8HXhrAl7zKLhGsPT9wyoivVNbW8vBBx9MU1MTzjmuv/76tE6+Bwol4CIDUfXDQPDtVeF+kN357Ddp5+P34RM/3y05uXBA/E8eHlDydobc8dC0BCKVULsAig8LOyoRSRGlpaUb5/2WzpnZDOBPQDZws3PuqnbrLVh/JFALnOOce9XMxgG3AZsDEeBG59yfujueasBFBqI20w9m2OwnT8U89r0OgtKy0EJJCWZQEnP9gNj3hoiIdMvMsoG/ADOBCcBpZjah3WYzge2C2wXA9cHyZuA7zrmdgL2Br3dw302kXAIeVk26pDa9L2JEGqA25uqXmTT/d1UFvPy/aP+Qnl+4Kq3Fvgeq/+OnqBSRftHfnYGpj6/bXsBi59zHzrlG4C7g2HbbHAvc5rwXgDIzG+WcW+mcezU4dhWwCOh2qpjQSlCGDh26yVySjY2NLF++fOPFbkRaNTc3E4lE+Pzzz8MOJXRD8l5m0tBqAOqaR/Pi86uADHleKtfDPkHSnV8AS1b4W4YzHPuNLCYnqwaaPuHlZ2+lpnmr7u8oIh0qKSlh2bJlDB48OG3noU5Hzjk2bNhATU1NR/OV55hZbC3Ojc65G4P2GGBpzLplwFTa6mibMcDK1gVmNh7YDXiRboSW6ZaXl7eZNB6gqamJZcuWdXtBGck8BQUFjB07VvODA3x+P6z3zcIRX2DaLgeHG0+ytLTAFedA+Rrf//JlsM+0MCNKLcuPgqp7AJiy0yoYdm7IAYkMXK35yPLly8MORXqpoKCASZMmdZQvNDvn9uzkbh39l9V+KL3LbcysBLgPuNg51+0k5Sk11Jybm7vJ5c9FJIZzUPXvaD+T6r/feCGafJcOhj0PCDeeVFNy9MYEnKo5MOyKcOMRGcCUj2ScZcC4mP5YoP3Xq51uY2a5+OT7Dufc/T05YMrVgItIFxpei7n6ZRkUHRRqOEkVO/XggTMhV9NotVFyJP7kffwUlc2rQw1HRGQAeRnYzsy2MrM84FRgTrtt5gBnmbc3sME5tzKYHeUWYJFz7g89PaAScJGBpOqBaLtkVuZc/XLFp/De676dlQUHHRVqOCkpe6ifkhIAB9UPhRqOiMhA4ZxrBi4C5uJPorzHOfeOmV1oZhcGmz0MfAwsBm4CvhYs3w/4InCImb0e3I7s7pgpVYIiIt2ITcBLjw8vjmSLnXpwt3391S9lUyVHQ90C365+EMpUBy4i0hPOuYfxSXbsshti2g74egf3e4aO68O7pBFwkYGi8UNofMe3rQCKM+QCNLU18Ny8aP/gDJp2sbdKj4u2ax6DSHVooYiISOeUgIsMFLGj38XTIas4vFiS6fl50BDMjDR6S9hhYrjxpLK8bSF/V9929VD9aLjxiIhIh5SAiwwUbeq/M6T8JBKBJ2POgznkGH/lR+lcyQnRdnWPTsYXEZEkUwIuMhA0rYD6F4JONpRmSBnGotfh82Ae3sIi2PuQUMMZEGLPDah+yF85VUREUooScJGBoDpm7u+iAyF7WHixJNNTMaPf+02HgsLwYhko8idC7ta+HamE2ifDjUdERDahBFxkIKh+MNrOlPKTtavgjZir+U7LoIsO9YcZlMaUoVSpDEVEJNUoARdJdS0VUBMzihk700U6m/+Qv/InwM57wOZjw41nIIlNwKsfBNcSWigiIrIpJeAiqa76IaDZtwv2hNxxXW6eFhob4OmYGTwOOSa8WAaigqmQM8q3W9ZC3TPhxiMiIm0oARdJddUZOPvJS/Ohpsq3h28Ou+4ZajgDjmW1fa+oDEVEJKUoARdJZZE6qH4k2s+E8hPnYN6D0f60oyArO7RwBqw2deAPRMt5REQkdErARVJZzePgan07b3vI2ynceJLh/Tdh2Se+nZcPB8wIN56BquhAyBrq281Lof6VcOMREZGNlICLpLLY0oGS4zPjIjSxo9/7Hg7FpaGFMqBZLpTEzBevi/KIiKQMJeAiqco1tp3/O7akIF2tWQlvvBDtH3pseLGkgzZlKPepDEVEJEUoARdJVTXzIFLh2zlbQsGUUMNJiifntJ16cFQGzPiSSMWHgxX7duMH0Lgo3HhERARQAi6SuirvibYHfSH9y0/qa+GZudH+YceFFkrayCqEkiOj/ap7w4tFREQ2UgIukooiDW2vfll6cmihJM1z86AuOOF087F+BFz6r/TEaLvqX+HFISIiGykBF0lFtfMgssG3c8f7C/Cks0gEnoipdz/0WMjSx1NclBwFVujbDW9Dw7vhxiMiIj1LwM1shpm9b2aLzezyLrabYmYtZnZS/EIUyUCx5SelGVB+8vbL8Ply3y4shn0OCzeedJJVAiWzov3Ku8OLRUREgB4k4GaWDfwFmAlMAE4zswmdbPcbYG77dSLSC5GGdrOffCG8WJJlXszjPWAGFBSGF0s6Kj0l2q66W7OhiIiErCcj4HsBi51zHzvnGoG7gI7mBvsGcB+wOo7xiWSe2sczq/xk+RJ491Xftiw45OguN5c+KJkZMxvK+9DwVrjxiIhkuJ4k4GOApTH9ZcGyjcxsDHA8cENXOzKzC8xsoZktbG5u7m2sIpmhMuZEuUwoP4mt/d5tHxi+eXixpKusIig9JtqvUhmKiEiYepKAd/TXv/33l38Evueca+lqR865G51zezrn9szJyelhiCIZJNNmP6mugheejPY19WDixJahVKoMRUQkTD3JgpcBsVfDGAusaLfNnsBd5kfqhgNHmlmzc+7BeAQpkjFqHoNIpW/njoeCNJ+K7+lHoLHBt8dtA9vtEm486az4CMga5N9fTR9Bw6vp//4SEUlRPRkBfxnYzsy2MrM84FRgTuwGzrmtnHPjnXPjgXuBryn5FumDyjuj7dJT0rv8pKUFnvpPtH/Ysen9eMOWVQAlMafvxM60IyIiSdVtAu6cawYuws9usgi4xzn3jpldaGYXJjpAkYwRqW47+8mg08OLJRleexbK1/h26WDYa1qo4WSEQbGzodyjMhQRkZD0qBDbOfcw8HC7ZR2ecOmcO6f/YYlkoKp/g6vz7fxdoGBiuPEk2uMPRNvTZkFuXnixZIriwyGrDCIV0LQE6l+Awn1CDkpEJPPoUnMiqSK2/CTdR78XvwsfLfLtnFw46Mhw48kUltf20vQb7ggvFhGRDKYEXCQVNK+BmphrWJWeGl4syTD33mh774OhbFh4sWSaQWdG21V3g2sKLxYRkQylBFwkFVT9Cwhm8SzcF/K2CjWchPp8Obz+fLR/+AnhxZKJig6EnLG+3bLWz7wjIiJJpQRcJBVkUvnJ4w9ET/7bdQqMGR9qOBnHstq+xypVhiIikmxKwEXC1rgE6p4NOtn+6pfpqqoCno0ZcT3ipNBCyWiDzoi2qx6ElqrQQhERyURKwEXCVnVXtF18OOSMDC+WRJv/EDQ1+vYW28IOaT7TS6oqmAj5u/q2q2t79VUREUk4JeAiYXKubQlA7MhkumlsgCdjruF1xIm68E6YYt9rlbeHF4eISAZSAi4SpobXoeFt37bCtlcqTDfPPwFVG3x76EjY44Bw48l0sXXgNfOgeVV4sYiIZBgl4CJh2vDPaLv0BMguDS+WRIpE4LH7ov3DjoOcHl0HTBIldxwUHhR0IlB5V5ebi4hI/CgBFwmLa2xbfjL4nNBCSbg3X/TTDwIUFsOBM8KNR7zBMXOCb/i/8OIQEckwSsBFwlL9sJ+HGfy8zEUHhxtPIs2NGf0+6EgoKAovFokqPQks37cbXoX6N8ONR0QkQygBFwlLbPnJ4LPAssOLJZE+fg8+DOrcs7Ph0DSucx9ossug5Phof8M/QgtFRCSTKAEXCUPzGqj+b7Q/6KzwYkm02MvO73UwDBkeXiyyqbJzo+3K231plIiIJJQScJEwVM4Gmn27cB/I3yHUcBJm1VJ49dlof7ouO59yig5te2n66ofCjUdEJAMoARcJw4Zbo+10Pvny0X9FLzu/yxQYt3W48cimLBsGnx3tqwxFRCThlICLJFv9m9Dwmm9bPpSeHG48iVK+xs/93eqoU8KLRboW+09g9cOaE1xEJMGUgIskW+zod8nx/kS4dPTYfdDS4tvb7Qzb7RJuPNK5vG2hsPXCSC2aklBEJMGUgIskU6QeKm+L9tO1/KSqAhY8Eu0feWpooUgPDY45GXPDP6KlQyIiEndKwEWSqfoBaFnn2zlbQvFh4caTKE/8GxobfHvcNrDLnuHGI90b9AWwYt9uXAT1L4Ubj4hIGlMCLpJMFTdF22VfTs+5v+tq4Ik50f6Rp4BZePFIz2SV+CS8VcXN4cUiIpLmlICLJEvjh1D7VNDJgsFfCjWchJn/sE/CATYbA3vsF2480nODz4u2K2dDy4bwYhERSWNKwEWSJXZEsWQW5I4JL5ZEaWyAx++P9md8AbLScJQ/XRXuC3k7+7argco7wo1HRCRNKQEXSQbX2HZ+5bLzw4slkZ59HCrX+/aQ4bDPoeHGI71jBkMujPYr/qaTMUVEEkAJuEgyVM2BljW+nTMWimeEG08iNDfBI/dE+0ecBDm54cUjfTPoTLBC3254E+pfDDceEZE0pARcJBkqboy2B38JLCe8WBLl2cehfLVvlw6GA9Lwn4xMkF0Gg06L9tffEFooIiLpSgm4SKI1fgy1jwcd87OfpJvmJnjormh/xhcgvyC8eKR/yr4SbVfdDS3rw4tFRCQNKQEXSbSK66Pt4hmQu0V4sSRK+9HvabPCjUf6p2AK5O/m265eV8YUEYkzJeAiiRSphYpbov0hXw8vlkRpboKHY0a/jzhJo98DnVnbUfCKG3QypohIHCkBF0mkyjshEnx9n7s1FM8MN55EeG4erAtGv0sGw8FHhxuPxMeg0/3FecBfGbN2fqjhiIikEyXgIoniHKy/Ltof8nWwNPuVa26Ch2ZH+xr9Th/ZpTDorGh//Z/Di0VEJM2kWTYgkkLqnoGGN3zbimDwueHGkwibjH6r9jutDLko2q7+NzQuCS0UEZF0ogRcJFFiR78HnwnZQ8KLJRE2qf0+EQoKw4tH4i9/JyieHnQiUPGXUMMREUkXSsBFEqFpOVTdF+2XXdT5tgPVc/Ng7ee+XTJItd/pasi3ou2KmyFSE14sIiJpQgm4SCJU/A1o8e3Cg6Bg11DDibumRvjPHdH+dI1+p63iGZC7rW9HKjQloYhIHCgBF4m3SH2QgAeGpOHo91P/gfVrfbu0DA45JtRwJIEsC4Z8I9pff62mJBQR6Scl4CLxVvl/0BKcmJgzFkqPDTeeeKurgYfvjvZnnabR73Q3+BzIKvXtxkVQOy/UcEREBjol4CLx5CJQ/vtof8jFYLmhhZMQj90P1ZW+PWwkHJiGc5tLW9mD2s7iU/7H0EIREUkHSsBF4qn6IWh837ezBkHZ+eHGE29VFT4Bb3XsWZCbF1o4kkRDLgLMt2sehoa3Qw1HRGQgUwIuEk/lV0fbZRf4kcN08vDd0FDn26O3hL0PDjceSZ687aDkuGh/3e9CC0VEZKBTAi4SL3UvQd2CoJMDQ74Zajhxt241PPXfaP/4syErO7x4JPmGfS/arrwTmpaGF4uIyACmBFwkXmJrvwedBrnjwoslEf5zu7/4DsDWO8LkfcKNR5KvcCoUHhh0mlULLiLSR0rAReKh8ROoujfaH/qd8GJJhJVL4dmYmS9OOBfMwotHwjPsu9H2hhuhZX14sYiIDFBKwEXiofz3QMS3iw6HgkmhhhN3993iZ3gBmLA77Jhmj096rvhIyN/FtyPVsP76cOMRERmAlICL9FfzSthwc7Q/7NLwYkmERa/D6y/4tpkf/ZbMZQZDL4v21/8JInXhxSMiMgApARfpr3VXg2vw7YIpfgQ8XURa4J6bov19DoXx24UXj6SGQaf6i0yBv+jUhn+EG4+IyACjBFykP5rXQMUN0f6wH6ZXbfRz82DpR76dlw/HnxNqOJIiLA+GfjvaX3cVRBrCi0dEZIBRAi7SH+uvAVfr2/mToGRWuPHEU30dPHBrtH/ESTBkeGjhSIopuwCyR/p281LYcGuo4YiIDCRKwEX6qqUc1l8X7Q+7Mr1Gvx/9F2wIZrgoGwYzvhBuPJJasorb1oKv+xW4xvDiERHpBzObYWbvm9liM7u8g/VmZtcG6980s91j1v3dzFabWY8vEawEXKSvyq+FSJVv502A0hPCjSeeytfAY/dF+8efA/kFoYUjKWrIVyF7hG83f6ZRcBEZkMwsG/gLMBOYAJxmZhPabTYT2C64XQDETgF1KzCjN8dUAi7SFy2VfvaHVsN+AJZGv073/wMag5reLbb1J1+KtJdVDENjZv1Zq1FwERmQ9gIWO+c+ds41AncBx7bb5ljgNue9AJSZ2SgA59wCoLw3B8yJQ9B9MnToUObPnx/W4UX6ZXzJPxhfUgFAbfNYXnptM2B+mCHFT2MD5A6BA070/c3GwoIF4cYkKSvbdmXqiMHkZW2A5k95/6Xvs7Iujc6FEJF0kWNmC2P6NzrnbgzaY4ClMeuWAVPb3b+jbcYAK/sUTF/uFA/l5eVMmzYtrMOL9F3zGvj4/uh1d8b9imm7pMkIcaQFfnUxLPnQ93fbF84+P9SQZABYdwWs8SWTOwy/jx22/pWfKUVEJHU0O+f27GRdRydwuT5s02Np9J25SJKs+7W/AiBA3s4w6PRw44mnp+dGk++cXDhZybf0wJCvQ/Yw325aAhU3dbm5iEiKWQaMi+mPBVb0YZseUwIu0htNS6Hir9H+iF+CZYcXTzxVV/ra71YzT4YRo8KLRwaOrBIYGjNpwNqfRf9JFRFJfS8D25nZVmaWB5wKzGm3zRzgrGA2lL2BDc65PpWfgBJwkd5Z+7OYq17uBSXHhBtPPN3/D6gJZnUZvrlPwEV6ashFba+OWf6HcOMREekh51wzcBEwF1gE3OOce8fMLjSzC4PNHgY+BhYDNwFfa72/mc0Gngd2MLNlZvbl7o5pzvW5fKVfiouLXU1NTSjHFumTxg/g4wlAi++PmwfFaVL7/cn7vva79fPgop/A5L3DjEgGooq/w6rg705WCWz9MeSMCDcmERHAzGqdc8Vhx9FKI+AiPbXmR2xMvosOSZ/kO9ICd1wXTb4n7qXkW/pm8Fl+TnzwJSjrfhluPCIiKUoJuEhP1L0AVXdH+yN+FV4s8fbkf9qeeHnqhV1vL9IZy2n7u7H+r9D4SXjxiIikKCXgIt1xDj6/JNovPREK208POkCtWw0P3BrtH3UqjBwdWjiSBkqOgcJ9g04TrL0y1HBERFKREnCR7lTdBfUv+LblwYjfhhtPvDjnS08a6n1/9BY68VL6zwxGXBXtV94Jdc+HF4+ISArqUQJuZjPM7H0zW2xml3ew/gwzezO4PWdmk+IfqkgIInWw+nvR/pCLIW/r0MKJq4VPw5svRftnXexLUET6q+gAKDkh2v/8W+Ai4cUjIpJiuk3AzSwb+AswE5gAnGZmE9pt9glwkHNuIvBz4EZE0kH5H6A5uPJs9ggY9oNw44mXmiqYfX20P20WbNv+11qkH0ZeDZbv2/UvQ+X/hRuPiEgK6ckI+F7AYufcx865RuAu4NjYDZxzzznn1gfdF/BXBxIZ2JpX+qtethr+c8geFF488XTPTVAZ/MqWDYMTzg03Hkk/eVvB0Euj/dWXQ0tVePGIiKSQniTgY4ClMf1lwbLOfBl4pKMVZnaBmS00s4XNzc09j1IkDKu/By6Yqz5/Fyjrdl79geHNF+HZx6L9078GRSkzNaqkk2GXQ05wUm/LKk1LKCIS6EkCbh0s6/DqPWZ2MD4B/15H651zNzrn9nTO7ZmTk9PzKEWSreaptl+Zj7zGT7E20FVXwT//FO1POQh23y+8eCS9ZZXAiN9E++uvgcbF4cUjIpIiepKALwPGxfTHAivab2RmE4GbgWOdc+viE55ICFwjfP61aL/0ZCg+LLx44umu62FDuW8PGgJnfD3ceCT9DTodCoILO7lGWPW16EWfREQyVE8S8JeB7cxsKzPLA04F5sRuYGZbAPcDX3TOfRD/MEWSaN3V0Pieb2eV+tHvdPDac/DCk9H+F78JJWlS0y6py7Jg8+vY+Oem9nGonB1qSCIiYes2AXfONQMXAXOBRcA9zrl3zOxCM2u9ZN6PgGHAX83sdTNbmLCIRRKp8RNY9/Nof/jPITcNLkxTVQG3XRvt73Mo7LZPaOFIhinYA4Z8I9pffTG0lIcWjohI2MyF9FVgcXGxq6mpCeXYIh1yDpYdDTUP+X7+ZBj/8sCv/XYOrvsJvPGi75cNg5/eAMWloYYlGaalCj6ZAM3LfH/wl2HUzeHGJCIZw8xqnXMpM+OAroQp0qrq7mjyjcHmNwz85Btg/n+jyTfA2Rcr+Zbkyy6Fza6L9jfcArULwotHRCRESsBFAJo/h88vivbLvgKFU8OLJ16WL/Fzfrc67DjYdUpY0UimKz0WSo6P9ld9BSL14cUjIhISJeAizvmZGVqCyXtytmg7ddpA1dQIN17lfwKM3QpO/FK4MYlsdq0/uRn8yc5rfxRuPCIiIVACLlJ1D1TfH+2Pujk9rnj5r5v9CDhAbh5ccLn/KRKm3LEw4rfRfvnVUPtsePGIiIRACbhktubV8HnMXNiDz4fiw8OLJ14WLoAnY2YLPeUrMHrL8OIRiVX2FShq/T1zsPJsiFSHGpKISDIpAZfM5RysujCm9GQcjLw63JjiYdVS+EfM3OW77QsHHRlePCLtmcGoWyBrsO83fQSrO7yAsohIWlICLplrw01Q/UC0nw6lJw31cP0voaHO90eMgnO/7RMekVSSO87Xg7eq+CvUPB5ePCIiSaQEXDJTw7vw+cXRftlXoXh6aOHEhXNw+5+jdd85ufDVK6GoJNSwRDo16ItQcmy0v/JsXxYmIpLmlIBL5onUw4rTwAWjxHk7w8jfhxtTPPzvYXj+iWj/jK/DFtuEF49Id8xg879B9gjfb17pk3AXCTcuEZEEUwIumWfN96DhTd+2fBg9G7IKw42pv957A2b/Ndrf93DY/4jw4hHpqZzNYNRt0X7No35mFBGRNKYEXDJL1X9gfUzd6cjfQ8Gu4cUTD2tWwQ2/hJYW399iGz/6rbpvGShKZsDQ70b7a34AdS+EF4+ISIIpAZfM0bgYVn4x2i85Bsq+Fl488VBfC9f9BKorfX/QELjoJ5BfEGZUIr034hdQsHfQaYblp0LL+lBDEhFJFCXgkhkiNbD8eIhs8P2cLWDzWwb2KHEkArdc3faky6//EIaOCDUskT6xXBhzF2SV+X7zp8G5Gi2hhiUikghKwCX9OQcrz4OGt33f8mHMfZAzPNy4+utfN8Frz0X7X/wGbDMhvHhE+it3Sxj1j2i/Zq4vRxERSTNKwCX9rf8TVN0V7W/2VyjcM7x44uGx++HxmDnMDz8B9hvg0yiKAJQeB8Niku7y30DlXZ1uLiIyECkBl/RW8xisvjTaL7sAyr4UXjzxsHCBH/1utft+8IUvhxePSLwN/xkUHxXtr/wS1L8eWjgiIvGmBFzSV8PbsPwkIKghLdgLRl7b5V1S3gdvw82/82U1ANtOgPO+C1nZ4cYlEk+WBaPvgLwdfN/VwbJj/TzhIiJpQAm4pKfmVbD0KIhU+X7OWBjzAGTlhxtXf3z6Ifz5R9Dc5Pubj/UznuQN4Mck0pnswTDmQcgq9f3mz2DpLIhUhxqWiEg8KAGX9BOphWVH+z/YAFklMPYhyB0dblz9sXwJXPMDqKv1/UFD4Fu/gJJBoYYlklD5O8Lou4HgG56GV2H5yeCaQw1LRKS/lIBLenFNfuqy+oXBgiwYfQ8UTAw1rH75fAX84fvRub6LSuCSX8KIzcONSyQZSmbC5tdH+zWPwKqvRsuwREQGICXgkj5cxJ+sVT0numyz6/wf8IFq3Wr4w+Wwodz38wvh4l/CuK3DjUskmcrObzszyoabYe1Pw4tHRKSflIBLenAOPv8WVN4eXTb0uzDkq+HF1F9rVsHvLvNJOEBuHnzzp7D1DuHGJRKG4T+HQTFXsl33U1h3dXjxiIj0gxJwSQ9rfwwV10X7ZRfAiKvCi6e/Pl/uk++1n/t+dg58/UewwwAupRHpDzMYdTMUx8x3v+YyWP/X8GISEekjJeAy8K39Faz7ebRfeoq/2M5Avcz8is/gt5dB+Rrfz8mFi34MuwzwiweJ9Jfl+dmMCg+MLvv861Dxj87vIyKSgpSAy8DlHKz5MayNqQ0tngmjbwMboPNif/oh/O670ZrvvHxfdrLrlHDjEkkVWUUw9r9QMDW6bNV5UHFraCGJiPSWEnAZmJzzife6n0WXFR0GY+71o2QD0buvwm+/C1UVvp9fAN/6OUzYPdSwRFJOdimMewTyJwcLIrDqXJWjiMiAoQRcBh7nfO3nul9HlxXPgLFz/OjYQPTiU/CnH0FDne8XlcAlv1LNt0hnsofAuMdiknB8Ocq634UWkohITykBl4HFNcHKs6H899FlxbOCK+YVhhZWnzkHc++Dm34DLcHFRYYMh+9d7S8zLyKdyxkBWzwJBXtHl635Lqz5oeYJF5GUZi6kD6ni4mJXU1MTyrFlgGqpghUnQc1j0WUlx8OYuwZm2UlTI9xxHTwT83hGb+Hn+R46Iry4RAaalipYfgzUzo8uG3wObP63gfnZICJxZ2a1zrnisONopRFwGRiaV8Fn09om34PPhzH3DMw/sFUV8Icr2ibf2+0M3/29km+R3souhbEPQ/GR0WUbboWlR0LLhtDCEhHpjBJwSX31r8KSvaDh1eiy4T8JRrdyQgurzz77CH75LfjwneiyfQ6Db/8aSkrDi0tkIMsqhLEPwuAvR5fVPgGf7gdNn4YWlohIR1SCIqmt8m5YeS644OREsmDzG/ylqQca5+DpR+HOv0Jzk19mBid9GaafOHDnLRdJJc7Bul/B2iujy7KHw+i7ofiQ8OISkVClWgmKEnBJTa4F1v6w7UwnWYNh9GwomRleXH3VUO/rvZ+bF11WUATnfw8mTe38fiLSNxtuh5VfAoJ/dsmCkb+DIZfon12RDKQEPKAEXDrVvApWnAG1T0aX5W0PY+ZA/g7hxdVXSz/2s5ysiPkafMx4+OoPYPNxoYUlkvZqn4XlJ0HLquiy0lNh8xt93biIZAwl4AEl4NKhmsdhxZnQsjq6rHgmjL4TsstCC6tPIhGY9wDcf2u05ARg38PgjIv8hXZEJLGaVvjZk+qejy7L3dZ/m1a4Z3hxiUhSKQEPKAGXNlwjrP0JrLsKaH1PGgy7Eob/eOBdWr58Dfz99/De69Fleflw6oVwwAx9BS6STJEGWP0tqPhbzMJcGPErGPptMM1HIJLulIAHlIDLRvVv+IvrNLwRXZa9GYy+A4oPDS+uvohE4H8Pw31/h/ra6PItt4PzvgujVHIiEpoNd8DnX4VIVXRZ0aGw+c2QNz60sEQk8ZSAB5SAC64Z1v3Wj3wTU6JRdBiMvh1yNgsrsr5ZtRT++Sf48O3oMsuCI0+Go8+AnNzwYhMRr/EjWHE61L8UXZZVAiN+C2Vf0Wi4SJpSAh5QAp7h6l6EVRdCw+vRZVYAI66CId8YWH8EG+rhkbvh0Xvb1npvNgbO+ba/wI6IpA7XBGt+BOW/BSLR5UXT/AmaeduFFZmIJIgS8IAS8AzVUgFrvg8VNxCt9QYKpsKofw6sWU6cg5f/B/+6GdavjS7PyoIZJ8PRp0PuALxKp0imqHvBX2eg8b3oMsuDoZfBsO9DVlF4sYlIXCkBDygBzzCuBTb8HdZc2XaGEyvwJ1kOvXRgXdXys8Uw+4a25SYA47eDsy6GLbYJJSwR6aVIPaz96aaj4TlbwmbXQMlxOmlaJA0oAQ8oAc8gNY/D6u9Aw1ttlxfPhM2ug7ytw4mrLz5fDv/+Pz/yHfu7U1oGJ54L+x7uR8BFZGCpfxVWfQ3qX2y7vHA/GPEbKNovnLhEJC6UgAeUgGeAuhdh7Y+hZm7b5TljYOSfoPSEgTOytO5z+M+d8NzjfqaTVtnZcOhxMOt0KEqZ32sR6QsXgQ3/gDXfg5Z1bdeVHOunLcyfEE5sItIvSsADSsDTWN3LQeL9SNvlVgTDLoeh3xk4tZVrVsFj98LTc9ueYAn+EvInnaepBUXSTUu5L0tZfz1tZmgiCwafBUMvH1jnq4iIEvBWSsDTUN0LsPYXUPNQuxUGg8+B4b+A3NFhRNZ7n34Ic++Fl5/2o2KxdpoMx50N2+wUSmgikiSNH8PaH0Llne1WGJSeBMOugILdQglNRHpHCXhACXiacM1Q9QCsv6btpZ4BMBh0Ggz7IeTvGEp4vRKJwDuvwGP3waLXN12/zU5w/Dmw46RkRyYiYap/DdZcDjWPbbqueIafNaXo4IFTUieSgZSAB5SAD3At632tZPm10Pxpu5UGpafA8B9B/gAYJd6wHp59DBY8AmtXbbp+p91gxkkwYXf9gRXJZLVPw7pfb1peB5A3AYZ8DQZ9EbIHJT82EemSEvCAEvAByEWg9n+w4Raoug9cfdv1lgeDTvejQal+olKkBd57w9d2v/ostDS3XW9ZMOUAOOIkfxl5EZFW9a/Buqug6l+0uZ4B+KtqDvoilF0A+ZP0T7tIilACHlACPoA0feprICtugaaPNl2fPRzKvupHf3I2T358PeUcLPsEXngCXpwPFes23aaoBPY7HA45Fkak8GMRkfA1fgDlf4LK2yBSven6/F1h0Fkw+AzIGZX8+ERkIyXgASXgKa5phR/dqbq7g9ruQP7uMORCGHQmZBUmN76ecg5WfAqvPQcvL4DlSzrebpudYNpRsMcBkJef1BBFZIBrqYTK/4P1f4HGRR1skAXF031pXunRkD0s6SGKZDol4AEl4CmocTFU/weq/g11C9jkq1WArME+4S77cuqe/R+JwCfv+6T7tef8xXM6UjoYphwEB8yAcQPoYkAikpqcg9qngjK9B8DVdbBRtj9hs/QEKD0+tb81FEkjSsADSsBTgGv0o9vVD/nEu/G9TjbMhuLDYdAZUHpiao52V6zzM5i88wq8+xpUV3a8XV4+7LYv7H2IP7kyJye5cYpIZmip9OfKVN4GtfM7365gDyg+wt8K9wHLTVqIIplECXhACXgIXIs/eaj2Sah9AmqfAVfbycYGRdNg0KlQcgLkDE9mpN2rroKP3oX33/RJd2elJQD5BbDrFJ94T5oKBQPkIkAikh4alwQlffdD/Qudb5dVCkWH+mS86CDI21EncYrEiRLwgBLwJIjUQv1Cf4Gcuuf8DCaRis63t0IoPgxKjva3VPlq1DkoXwMfvg0fvgOL3+k64QYoLfNJ9+77wYTdVNctIqmhaZkvT6m6D+qeAVo63zZ7OBTuD0UH+J8Fu2mEXKSPlIAHlIDHmWuCxveh/vUg4X4eGt6gyw93gNxtgqR7lh95Cbu8xDlfTvLph/DpYv9zyYdQub7r+2XnwHY7w857+NvYrSArKzkxi4j0RcsG/21kzVyontvBNRXasSIo2B0K9vSlKwV7QN72YNnJiVdkAFMCHlAC3g8t5VD/hk+wG97w7cZ3fE13d3JGQdEhwdech0DulomPtzN1NbDiMz9LSett2Sf+wjjdycqCLbb1SfdOk2H7iVCQgrXpIiI94ZwfRKmZG5QJPgOR8u7vl1UC+btB/kTI3zl600wrIm0oAQ8oAe+Ga4LGj/0HcuP7fr7Z1nbL6p7vJ28CFO7tT+4p3C/5NYXNTbD2c3+FydUrYM1KWLnUJ9vla3q+n/xCP1XgdjvDtjvD1jv62m4RkXTkIv7E+NqnfalK7dPdj5DHyt7MXxAtb2fI3xFyt4W8bfygi8pYJAMpAQ9kfALuGn0tYNOnwW2J/3Btirl1Vz7SXs4WUDDJfy1ZuA8U7AXZZQkIPkZTI6xf68tG1q/1SfWalbB6JaxZAeVr/R+S3sgvhC228VegHL+d/7nZGJWUiEhma1oB9a8Et4X+1vJ5L3eS7ZPw3G2ChHxryB0HOeMgdyzkjFaCLmlJCXggbRPwSD00r4KWVdC80rebg3ZLa3uFv3U0z3ZPWL4f1SiY5C91nD/Jt7OHxOkxRKC2BqoqoGqD/1lZAdUboon2+rWwfp1f1lfZOT6xHr0ljNnS/xy9JWw2GrJU0ygi0iXn/N+S+ld9GWJDcGtcBK6+jzs1fwJ+ztiYpHwUZI+EnJH+Z2s7SzNKycChBDyQ0gm4i/jLCkeqIFLpZw5pWedrr1vK27Yjsf11fvt4yRnnT7DJ2yG4Be3cLXp20k1TI9RW+2S6ttrXXMf2N95qoKYqSLgr/BzaLb0cfe+MGQwZDiNG+dvIUTBytE+0R47RPNwiIvHmWqDpk2hC3rQYGj/yP5tXxO84Vhwk5ZtBzgjIGuq/dc0eAllDgp9BP7ZthZpeUZJuQCbgZjYD+BOQDdzsnLuq3XoL1h8J1ALnOOde7WqfoSTg5X+A5tUxiXXsz9h2dRKCMcjaHLLHgY0GRoPbHNxm0LIZNA+HpixoqPe3xoagXRezrD7ajl1WXw/1tb7+OtGysqBsmE+yy4bBkBEwYvPgNhqGbwa5eYmPQ0REuhep9cl540fQ9JFvNy2H5qXQvMx/S9vXb2d7yvLASvwJpJ3dOlpvRZBVABbcYtub9PUtqrTVXQLen1y3u/t2eLzuEnAzywY+AA4HlgEvA6c5596N2eZI4BtBUFOBPznnpna131AS8HdHQNbaxB4jkgX1hVBXBHUFUJMP1flQnQtVOVCVDdV5UFXgt01VhUVQMhgGlfk5tUsH+5+DymDoCJ9wDxnOgx/V8rvHP2RFRR2jywq57IgdOG63MQkP78oH32L2i0tpcY5sM06bOo5fHLdrp9s/+Npyfjf3/V7F2dl9Um15MvT2+e5MmI9BJFHS6n3tmvwoedMyn5A3LfUn/jev9vXmzat9v2V1z2beCk3Opgm65UZv5Hbd78kysoJEP7tt27I6X2bB8th2V/uwrGB7i/40C/rtbp0uz+rl9rH9rF5sT8xyNv0Z8rceXSXg/cl1e3LfDo/ZgwR8H+Anzrkjgv4VAM65X8ds8zdgvnNudtB/H5jmnFvZ2X5DScBfGAxlPSwRacwObjnQkAv1uVAX/OyqXZ9L9E0XsuwcKCqBouLoz8ISXl7dwCtrmqi0fH/LymeD5bPXpK34+tF7+mS7B6PWD762nCvuf4u6pmi5SmFuNr8+YdeE/vG58sG3uP2FzzZZfubeW3SYFPYlzs7uc+IeY7jvleUpszzRzzX0/vnuTFjvF5FEytj3tXP+W+OW1dD8ObSsgZb1/hZZDy0Vwc/WZRXRdSmduEvC5W4N23yU9MN2k4D3OdcFxnd33470pAB3DLA0pr8Mn/l3t80YoNMEPBTv7QRN5dHEujEHGnJ8v6ld2yUwic7Kgpxcf8vN9cluaz8npp+b62cEyc+HvAI/7V5+gW8XFEbb+fnBusK26/PyO/yP89QrHqaleNN/vJ55z/H1s0f0+GH8bu77bf7oANQ1tfC7ue8n9A/P7BeXdrq8o4SwL3F2dp/WUeBUWZ7o5xp6/3x3Jqz3i0giZez72gyyB/tb3nY9v59z/gTRSHX05mra9tusq445J6vO37f1FqnvvJ/oMhpJN/3JdXty3030JAHvKBNt/87uyTaY2QXABQBjxoxh/vz5PTh8HJVd5k8uNIh+bUInPy1mu95uH3Nr34/XVzARoA6oawAagJ6f/HnxLp3XhvfmNTl1XBWM62hNVUJf297G35c4O79Pqknscw0D//0ikkh6X8dLDlAW3PrLYTSTZY3+RvDTmjGaseBnlrX4frAsixasdVn7bawZoyW6D5oxc0AEC25YTJsIZjHtoO9jawm2j2njgmPHbLNxe//Tt/1Pn2+09mm3PqZvrf1I19taV/uKua91fBwgZhnRZbbpP0J1dfW8GM7vRo6ZLYzp3+icuzFo9yfX7VEOvEkw3W2Az+RjP17GAu1Po+7JNgQP9EbwJSjTpk3rweEl3r58xcObjKYCZJvx0RnTeryfH1z1JMsr6jZZPqaskG/0Yj+91dv4+xJnZ/fJNuv02GEsT/RzDQP//SKSSHpfi/SCcxTimGahnAPX7Jzbs5N1/cl183pw30305Bl4GdjOzLYyszzgVGBOu23mAGeZtzewoav6bwnXaVM7HtrtbHlnLjtiBwpz255pXpibzWVH7NDn2Hqit/H3Jc7O7nPa1HEptTzRzzUM/PeLSCLpfS3SC9Z6UmjK6U+u25P7bqLbEXDnXLOZXQTMxZ+m+3fn3DtmdmGw/gbgYfxZoYvxU7Oc27PHK2Fordvt76wWrfWNyT77v7fx9yXOru6z55ZDU2p5og3094tIIul9LTLw9SfX7ey+3R1TF+IRERERkbSWahfiScnvAURERERE0pUScBERERGRJFICLiIiIiKSRErARURERESSSAm4iIiIiEgSKQEXEREREUkiJeAiIiIiIkmkBFxEREREJImUgIuIiIiIJJEScBERERGRJFICLiIiIiKSRErARURERESSyJxz4RzYLALUhXJwyAGaQzq2JIde48yg1zkz6HXODHqd01+Yr3Ghcy5lBp5DS8DDZGYLnXN7hh2HJI5e48yg1zkz6HXODHqd059e46iU+U9ARERERCQTKAEXEREREUmiTE3Abww7AEk4vcaZQa9zZtDrnBn0Oqc/vcaBjKwBFxEREREJS6aOgIuIiIiIhEIJuIiIiIhIEmVUAm5mM8zsfTNbbGaXhx2PxIeZjTOzp8xskZm9Y2bfCpYPNbPHzezD4OeQsGOV/jGzbDN7zcz+G/T1GqcZMyszs3vN7L3gd3ofvc7px8wuCT6v3zaz2WZWoNd54DOzv5vZajN7O2ZZp6+rmV0R5GTvm9kR4UQdjoxJwM0sG/gLMBOYAJxmZhPCjUripBn4jnNuJ2Bv4OvBa3s58IRzbjvgiaAvA9u3gEUxfb3G6edPwKPOuR2BSfjXW69zGjGzMcA3gT2dc7sA2cCp6HVOB7cCM9ot6/B1Df5OnwrsHNznr0GulhEyJgEH9gIWO+c+ds41AncBx4Yck8SBc26lc+7VoF2F/4M9Bv/6/jPY7J/AcaEEKHFhZmOBo4CbYxbrNU4jZjYIOBC4BcA51+icq0CvczrKAQrNLAcoAlag13nAc84tAMrbLe7sdT0WuMs51+Cc+wRYjM/VMkImJeBjgKUx/WXBMkkjZjYe2A14EdjMObcSfJIOjAwxNOm/PwLfBSIxy/Qap5etgTXAP4JSo5vNrBi9zmnFObccuBr4DFgJbHDOPYZe53TV2eua0XlZJiXg1sEyzcGYRsysBLgPuNg5Vxl2PBI/ZjYLWO2ceyXsWCShcoDdgeudc7sBNagMIe0ENcDHAlsBo4FiMzsz3KgkBBmdl2VSAr4MGBfTH4v/ykvSgJnl4pPvO5xz9weLPzezUcH6UcDqsOKTftsPOMbMluDLxw4xs9vRa5xulgHLnHMvBv178Qm5Xuf0chjwiXNujXOuCbgf2Be9zumqs9c1o/OyTErAXwa2M7OtzCwPX/g/J+SYJA7MzPA1o4ucc3+IWTUHODtonw38O9mxSXw4565wzo11zo3H/+4+6Zw7E73GacU5twpYamY7BIsOBd5Fr3O6+QzY28yKgs/vQ/Hn7uh1Tk+dva5zgFPNLN/MtgK2A14KIb5QZNSVMM3sSHwdaTbwd+fcL8ONSOLBzPYHngbeIlof/H18Hfg9wBb4D/wvOOfanxwiA4yZTQMudc7NMrNh6DVOK2Y2GX+ibR7wMXAufrBIr3MaMbOfAqfgZ7F6DTgPKEGv84BmZrOBacBw4HPgx8CDdPK6mtkPgC/h3wcXO+ceSX7U4cioBFxEREREJGyZVIIiIiIiIhI6JeAiIiIiIkmkBFxEREREJImUgIuIiIiIJJEScBERERGRJFICLiIiIiKSRErARURERESS6P8BQvnBQ0ovIgMAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"beta1_hat, beta2_hat = log_reg.params[0], log_reg.params[1]\n",
"X = np.linspace(0, 100, 500)\n",
"\n",
"Z = beta1_hat + beta2_hat * X\n",
"dPdX = np.exp(-Z) / (1 + np.exp(-Z)) ** 2 * beta2_hat\n",
"\n",
"fig, ax1 = plt.subplots(figsize=(12, 7))\n",
"ax1.scatter(df[\"Annual_Family_Income\"], df[\"Own_House\"])\n",
"ax1.plot(X, P_hat, color=\"tomato\", lw=3, label=\"Accumulative Effect\")\n",
"ax1.legend(loc=\"center left\")\n",
"\n",
"ax2 = ax1.twinx()\n",
"ax2.plot(X, dPdX, lw=3, color=\"Gold\", label=\"Marginal Effect\")\n",
"ax2.legend(loc=\"center right\")\n",
"ax2.set_title(\"Accumulative and Marginal Effect of Family Income on Owning a House\")\n",
"ax2.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "5d22484f-81f9-4c0b-a2e2-ec5f520dcab6",
"metadata": {},
"source": [
"As family income raises, the marginal probability will reach maximum around $40000$ pounds, to summarize the effect, we calculate the marginal effect at the mean value of the independent variables. "
]
},
{
"cell_type": "code",
"execution_count": 107,
"id": "8aee3923-548b-448d-ab30-193a6e3abe5e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.03297676614392626"
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X_mean = np.mean(df[\"Annual_Family_Income\"])\n",
"Z_mean_effect = beta1_hat + beta2_hat * X_mean\n",
"dPdX = np.exp(-Z_mean_effect) / (1 + np.exp(-Z_mean_effect)) ** 2 * beta2_hat\n",
"dPdX"
]
},
{
"cell_type": "markdown",
"id": "9ee4701f-f876-4246-a110-02280df28a81",
"metadata": {},
"source": [
"At the sample mean, $1000$ pounds increase in family income increases the probability of owning a house by $3.2 \\%$. "
]
},
{
"cell_type": "markdown",
"id": "3f46c537-b4cf-46d2-bec6-06a1f3676b80",
"metadata": {},
"source": [
"## More Than One Variable"
]
},
{
"cell_type": "code",
"execution_count": 138,
"id": "cf5d81fc-3eaf-473b-b26c-17cdd394b294",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.247296\n",
" Iterations 8\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: admitted No. Observations: 30\n",
"Model: Logit Df Residuals: 26\n",
"Method: MLE Df Model: 3\n",
"Date: Thu, 09 Sep 2021 Pseudo R-squ.: 0.6432\n",
"Time: 12:16:41 Log-Likelihood: -7.4189\n",
"converged: True LL-Null: -20.794\n",
"Covariance Type: nonrobust LLR p-value: 6.639e-06\n",
"===================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"const -16.3315 8.601 -1.899 0.058 -33.189 0.526\n",
"gmat 0.0025 0.018 0.141 0.888 -0.032 0.037\n",
"gpa 3.3208 2.397 1.385 0.166 -1.378 8.020\n",
"work_experience 0.9975 0.585 1.704 0.088 -0.150 2.145\n",
"===================================================================================\n"
]
}
],
"source": [
"df = pd.read_excel(\n",
" \"Basic_Econometrics_practice_data.xlsx\", sheet_name=\"GMAT_Work_Quali_Resp\"\n",
")\n",
"df.head()\n",
"\n",
"Y = df[[\"admitted\"]]\n",
"X = df[[\"gmat\", \"gpa\", \"work_experience\"]]\n",
"X = sm.add_constant(X)\n",
"log_reg_house = sm.Logit(Y, X).fit()\n",
"print(log_reg_house.summary())"
]
},
{
"cell_type": "markdown",
"id": "5b597b6b-ee22-4e90-9e26-6ba747c8104e",
"metadata": {},
"source": [
"The estimated model is\n",
"$$\n",
"\\hat{P_i} = \\frac{1}{1+e^{-(-16.3315+0.0025 X_{2i} + 3.3208 X_{3i}+.9975 X_{4i})}}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 131,
"id": "b3038355-3619-43f5-9be2-87a350195b4c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Marginal effect of GMAT, GPA and Working Exp are 0.0005751109753324387, 0.7627645112662916 and 0.22911309163038573 accordingly.\n"
]
}
],
"source": [
"beta1_hat, beta2_hat, beta3_hat, beta4_hat = (\n",
" log_reg.params[0],\n",
" log_reg.params[1],\n",
" log_reg.params[2],\n",
" log_reg.params[3],\n",
")\n",
"X2_mean = np.mean(df[\"gmat\"])\n",
"X3_mean = np.mean(df[\"gpa\"])\n",
"X4_mean = np.mean(df[\"work_experience\"])\n",
"Z_mean_effect = (\n",
" beta1_hat + beta2_hat * X2_mean + +beta3_hat * X3_mean + beta4_hat * X4_mean\n",
")\n",
"dPdX2 = np.exp(-Z_mean_effect) / (1 + np.exp(-Z_mean_effect)) ** 2 * beta2_hat\n",
"dPdX3 = np.exp(-Z_mean_effect) / (1 + np.exp(-Z_mean_effect)) ** 2 * beta3_hat\n",
"dPdX4 = np.exp(-Z_mean_effect) / (1 + np.exp(-Z_mean_effect)) ** 2 * beta4_hat\n",
"print(\n",
" \"Marginal effect of GMAT, GPA and Working Exp are {}, {} and {} accordingly.\".format(\n",
" dPdX2, dPdX3, dPdX4\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"id": "f137a00b-8f62-4cb6-92d6-f91e7a123cbf",
"metadata": {},
"source": [
"## Goodness of Fit"
]
},
{
"cell_type": "markdown",
"id": "7f128c30-2e12-4c1f-94ef-f31929e23579",
"metadata": {},
"source": [
"You might be wondering why log-likelihood is a negative number as in the report, recall that any likelihood is a joint distribution which must lie between $[0, 1]$, so the log likelihood must be negative, the plot below is a reminder."
]
},
{
"cell_type": "code",
"execution_count": 137,
"id": "6abb2172-52eb-43d7-8665-0fe96ca02ffa",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAGvCAYAAAC6pwjDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0FElEQVR4nO3deXhc1Z3m8fdo3/fVsmx5kxe8WxgIMcgGAiQEQodsMKGb7o5D95Okk0l6phM6k16m8/T0kp6ZhHSWgU5nw4EQICEkbInAJICxjXfLK7asfS8tJalUqjN/VFnYwotslepe1f1+nkePVFVXdX/yz1d+OZx7jrHWCgAAAPCSBKcLAAAAAGKNEAwAAADPIQQDAADAcwjBAAAA8BxCMAAAADyHEAwAAADPIQQDAADAcwjBABDHjDE/MMasjnz9OWPMVx0uCQBcIcnpAgAA0+qbkv7cGPOQpPdKusXhegDAFQw7xgFAfDPGvCwpTdId1toWp+sBADdgOgQATDNjTL4xpt8Yc/2E539gjPmZMcZM47mNpGFJvyYAA8DbCMEAMM2stT2SviPpc6efM8Z8WdIySf/FTu//kvtLSW9Kes90hm0AmGkIwQAQG1+TdLMxZoEx5kOSNkt6v7XWP9U3NsZ8whhzwBjjM8b8yhhTEnn+3ZI2SvqipFcl3TTVcwFAvGBOMADEiDHmYUmzJK2TdLO1dueE15+W9O7zfPsr1trbzvGeX5L0QUkfkdSg8I1wIUkPSHpG0vuste3GmGpJ/2St/UCUfhwAmNEIwQAQI8aY5ZL2SvqItfbRKLxfiaQTklZbaw9HnrtG0jettWum+v4AEM9YIg0AYidF0oikn0Xp/W6IvOe2M6b7GoXnAAMALoAQDACxs0rSPmtt8FwvGmN+JWnDeb53q7X21gnPFUh6wlr7oSjWCACewI1xABA7qyXtOt+L1tpbrbVZ5/mYGIAlaaekjcaYtZJkjMkxxtzBKhAAcHGEYACInVW6QAi+VNbaVyX9naTHjTEDkg5IumWal1wDgLjAjXEAAADwHEaCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmObJZRVFRkq6qqYna+wcFBZWZmxux8mBz64k70xZ3oizvRF/ehJ+7kZF927NjRaa0tnvi8IyG4qqpK27dvj9n56urqVFtbG7PzYXLoizvRF3eiL+5EX9yHnriTk30xxpw81/NMhwAAAIDnEIIBAADgOYRgAAAAeA4hGAAAAJ5DCAYAAIDnEIIBAADgOYRgAAAAeA4hGAAAAJ5DCAYAAIDnEIIBAADgOYRgAAAAeA4hGAAAAJ5DCAYAAIDnEIIBAADgOUlOFwAAAID4NBIcU0f/iJr6Q06X8g6EYAAAAFySQDCkjoERtfUNq71vRO39w+Nft/WPqL1vWO39I+oeDEiSSjKM7nm/w0VPQAgGAACAJGksZNU1MKK2vnDAbesfVptvOPy4P/y5vW9YXZFwe6bEBKOirBSV5aRpdn6G1s3NV0l2mkpzUtXZcNiBn+bCCMEAAABxzlqrvuFgONj2DavVFx6pbfUNjz/X1jeijoERjYXsWd+bYKSirFSV5qSpIi9NqyvzVJYTDrclOamRoJumgswUJSaYc56/zn88Fj/mJYlKCDbGPCzpNknt1trl0XhPAAAAXNzoWEgd/SNqjYTb08H29OPTAXdodOwd35ubnqzSnHDAXVSaPR5uS3PSVJKTprKcNBVlpSgpMf7WUojWSPD3JH1D0vej9H4AAACe5w8E1eILT0lo8b0dbM/83DkwInv24K1SEhNUkpOqspw0XVGRqxuWhgNtaW7aWUE3LTnRmR/MBaISgq21LxtjqqLxXgAAAPHu9PSEVt+wWnxDkc/hYNvSdzr0DqlvOPiO781NTx4PtMvKc8aDbVluONiWRaYmGHPuqQkIY04wAABAFFlr1TcUVEvfkFp6w+G2xTf0dsiNfO0PnD09wUTm3pbnpmluYYaunl9wRsBNU3luuspy0pSe4t3R22gyduL4+eW+UXgk+OnzzQk2xmyWtFmSSktL123ZsiUq552MgYEBZWVlxex8mBz64k70xZ3oizvRF/eJRU+Gglbdw1bdQ6Hw5/GPtx+PTJh+ayTlpRrlpxkVRD7y0xLO+NooL9Uo6Tw3ls10Tl4rGzdu3GGtrZn4fMxC8Jlqamrs9u3bo3Leyairq1NtbW3MzofJoS/uRF/cib64E31xn6n2ZCQ4plbfsJp7wyO2zb1DavYNq7k3PKrb7BtS/4QpCsZIJdmpKs9NV3lkxLY8N03leW9/XZKdGpc3l02Wk9eKMeacIZjpEAAAwBOsteocCISD7Rnh9szHHf0j7/i+wswUleelaU5kikJ5XjjYzop8Ls1JU7KHA+5MFa0l0h6RVCupyBjTKOkr1tqHovHeAAAAkzE8OqbWwZBeOdKp5t4hNY2H2yE19w6rqXdIgeDZ2/emJyeqIj8cZpeW56g8N12z8sIB93TI9fIKCvEsWqtDfCwa7wMAAHAup282a+z1q6lnaDzkNvUOqalnSE294aXCJElbX5cUnqZQmp2mWXlpumJWjt6zrHQ82Fbkp6siL1256cmsouBRTIcAAACOs9aqazCgxp5wqG3s8Z8RcIfU2DOkgZGz5+KmJiWoIjJie8OSHFXkp6u/9YQ2Xb1Ws/PTVZqTppQkping3AjBAABg2p2ej9vY41djz1DkI/x1OOT6NTx69lSF7LQkzc7P0Oz8DF09v1AVeenjI7gV+ekqPMdauHV1TbpmQWEsfzTMUIRgAAAwZdZadQ8GdOqMcHuq239W2B2ZMB83LyNZlfkZWlicpdrqYs3OT1dFfoYq8tI1uyBdOWnJDv008AJCMAAAmJS+4VGd6vbrVPfZQfdU5OuJmz/kZSRrdn66FpVka9OSksiobrpm52eoIj9dWanEEDiHv30AAECSFAiG1Nw7pIZIsG3o9qux++3Hvf7Rs47PSk3S7Px0zS3M1LsXFquyIP2MoJuubEZy4WKEYAAAPMJaqx7/qBq6wwH3VLdfDV3+8cctviGFzthDKyUxQRX56aosyNDK2bmqLMhQZX6GKgvSVZmfobwMVlbAzEUIBgAgjoyOhUdzT54Rbhu6/DoZCb0TV1gozk7VnIIMrZ9XEAm56ZpTkKHKggyV5qQpMU638QUIwQAAzDD+QFAnu/yRoDs4HnhPdA2quXdYY2cM56YkJYwH26siQXduQYbmFIZHddNT2AgC3kQIBgDAhXxDo2roCgfbk12DOtHlH/88cWvf3PRkzS3M0OrKfN2+Kl1zCzI1pzBDcwszVJqdpgRGc4F3IAQDAOCQXn9AJ7r8OtE5GAm7/vHP3YOBs44tzUnV3IJM1VYXq6ooU3MKwiF3bkGmcjO4AQ24VIRgAACmUd/wqE50DuqtyEc48IbD7pmrLRgjzcpNV1VRhm5ZXqaqwgzNKcjUvEjgZdoCEF2EYAAApmgoMKYTXW8H3TMDb9cZI7pnBt33rijXvMJMVRVlqqowfCNaWjJBF4gVQjAAAJMQHAupsWdIb3UO6ljHgN7qHNTOI0P64qsvqsU3fNaxJdmpmleUqZuWlaqqKDyae3pEl6ALuAMhGACAiNNb/x7vHNTxjgEd7xjUsY5BvdU5oIZuv0bH3l51ISctSUWp0jXzC8MhtzhTVYXhsJvJTmiA63GVAgA8Z3QspJNdfh0bD7oDOt4xoGMdg/INvT1PNyUxQXMLM7SgOEs3LSvT/OJMzS/K1PziLOVnJOull15Sbe1q534QAJeNEAwAiFu+oVEd6xjQsfZwwD3aHg67J7v9Z62lW5KdqgXFWbptZbnmF2eNh93Z+RlsFgHEKUIwAGBGs9aqtW9YR9sHxj+OdQzoaPugOgfeXk83JTFBVUUZWlyWrfeuKNeCkkzNLwoH3uw0lhgDvIYQDACYEcZCVo09fh1pG9CR8cDbr2Mdg2dtBZyTlqSFJVnauLhYC0uytLAkSwuKszQ7P11JiQkO/gQA3IQQDABwleBYSCe7w2H3aHu/jrQP6EhbeHR3JBgaP640J1ULS7L0wbUVWliarYXFWVpQkqnirFQZwxQGABdGCAYAOOLtsNuvw20DOtzWryNt4aXHAmNvh92KvHQtKs3SuxYUqro0WwtLwyO7uelMYQBw+QjBAIBpFQpZNfYM6XBbvw619Yc/t/breMfZYXd2frqqS7NVu6RY1SXZWhQJuyw3BmA68JsFABAV1lp19I/oUCTkHmoNB97DbQMaGh0bP+70yO511cVaVJIVHt0tIewCiC1+4wAALtnASFCH2/pV39KvQ619qo8E3h7/22vsFmWlanFZlj66vlKLS7NVXZatRSVZrMQAwBUIwQCA8xoLWZ3sGlR9a7/qW/p0sLVf9a19OtU9NH5MZkqiqsuydcvysvGwu7g0W4VZqQ5WDgAXRggGAEiSfP5RHWzt08GWPtW3hMPuobZ+DY+G5+0mGGleUaZWzs7TR2oqtbgsR0vKslWRl64ENpQAMMMQggHAY0Ihq4Zuvw629OlASzj0HmzpV1Pv26O7BZkpWlqerXuumqslZdlaWp6jhSVZSktOdLByAIgeQjAAxLHh0TEdbuvXgeZw4D3QHA69g4HwjWoJRppfnKW1c/N1z9VztLQ8R8vKc1SSzVq7AOIbIRgA4oRvaFQHmvu0v9kX+dynox0DGgtZSVJWapKWlmfrrnWzw2F3Vo6qS7MZ3QXgSYRgAJiB2vuHtb+5T/ubfNrf3Kd9zb6zblYrzUnVsvIc3bSsVFfMCgfeyvwM5u4CQAQhGABczFqrtr4Rvdke1JvPH9a+Jp/2NvnU3j8yfkxVYYZWVuTpY+vn6IpZuVpWnqPibFZmAIALIQQDgIu09Q1rT6NPext7tbfJp71NfeocCAfeBHNEC4qzdO3CIi2vyB0f4c1h3V0AuGSEYABwSOfAiPY09kZCr097mnzq6D8deKWFJVm6rrpIKypyFWw/rnved70yUvi1DQDRwG9TAIiBvuFR7Wv0aXejbzz4nl6SzBhpYXGWNiws0orZuVpRkatls3LOCrx1dScJwAAQRfxGBYAoGwmOqb6lX7sbe7X7lE+7G3t1rGNANrxIg+YUZGjNnDz90buqtHJ2rpZX5CozlV/HABBL/NYFgCmwNrzxxK5TvXqzoVe7TvXqQHOfAmPhXdaKslK1ujJXd6yapZWVeVpZkav8zBSHqwYAEIIB4BL0DY9qdyTwvtnQo92NPnUPBiRJ6cmJWjE7V/ddW6VVlXlaVZmnWblpbDoBAC5ECAaA8wiFrI52DGjnyZ5w6D3VoyPt4WkNp+fx3ri0RKsr87W6Mk/VpVlKSkxwumwAwCQQggEgon94VLtP+bTjZI92NPTozYYe9Q8HJUl5GclaU5mn21bO0to5+VpZmcvSZAAwgxGCAXiStVaNPUPacbJH2092a/uJHh1q6x8f5V1cmq33rwoH3rVz8jSvKJNpDQAQRwjBADwhOBbSwZZ+bTvRrR2R0Ht617Ws1CStmZOnW5aXad3cfK2qzGOUFwDiHCEYQFzyB4La1dCrN0706I0T3drZ0CN/YEySVJGXrmsWFKpmbr7WzS3Q4rJsJSYwygsAXkIIBhAXfP5RbT/ZrW1vdev1t7q1r8mnYMjKGGlJWY7uWjdbV1YVqKYqX+W56U6XCwBwGCEYwIzUOTASDrzHu7TtRI/qW/tkrZSSmKBVlbnafN18XTmvQOvm5jO1AQDwDoRgADNCe/+wXj/erdff6tLrx7t1pH1AUnht3nVz8/W5G6u1fl6BVlfmKS050eFqAQBuRwgG4EpdAyN67Xi3Xj3eqVePdelYx6AkKTMlUTVVBbpzbYWunl+oFRW5SmZtXgDAJSIEA3AFn39Ur73VpVePhT8OtfVLCofeK+cV6EM1lbp6fqGWz8phQwoAwJQRggE4Yigwpu0nu/W7o136/bFO7WvyKWSltOQEXVlVoNtXz9I1CxjpBQBMD0IwgJgYC1ntbfLplSMdeuVop3ae7FVgLKSkBKM1c/L06U2LdO3CIq2uzFNKEqEXADC9CMEAps3JrkFtPdKpV4506vfHOtUX2YJ4WXmO/ujaKr1rQaGurCpQZiq/igAAscW/PACipm94VK8e69LWIx16+XCnGrr9kqRZuWm6ZXmZ3r2oWNcuKFRhVqrDlQIAvI4QDOCyhUJW+5p9eulQh14+0qGdDb0aC1llpiTqmgWF+pN3z9OGRUWaV5QpY9iRDQDgHoRgAJeka2BEW4906qXDHXr5cIe6BgOSpBUVubr/+vnasKhYa+fkM68XAOBqUQnBxphbJP0fSYmS/p+19h+j8b4AnBeK3ND2m/p21R1q154mn6yVCjJTdN2iItUuLtGGRUVMcQAAzChTDsHGmERJD0q6SVKjpDeMMT+31h6Y6nsDcIZvaFRbj3Tot/UdeulwuzoHAjJGWl2Zp8/dWK3rq4u1oiJXCQlMcQAAzEzRGAleL+motfa4JBljtki6QxIhGJhB3uoc1IsH2/TTbUM68tzzGgtZ5aYn6/rqYm1cUqzrq0tUkJnidJkAAESFsdZO7Q2MuUvSLdbaP408/rikq6y1n5pw3GZJmyWptLR03ZYtW6Z03ksxMDCgrKysmJ0Pk0NfnDUWsjraG9Kb7WPa1RFU62D4d0F5htXa0hStKknUgtwEJTLa6wpcL+5EX9yHnriTk33ZuHHjDmttzcTnozESfK5/Id+RrK2135H0HUmqqamxtbW1UTj15NTV1SmW58Pk0JfY8weCevlwp5470Krf1rerxz+q5ESjq+cX6v4lJbphaamO7dlGX1yI68Wd6Iv70BN3cmNfohGCGyVVnvF4tqTmKLwvgCjoHBjRiwfb9Nz+Nr1ytFMjwZBy05O1aUmJblpWqg2LipSdljx+/DEHawUAIFaiEYLfkLTIGDNPUpOkj0q6OwrvC+AyNfb49ez+Nj27v1VvnOiWtdLs/HTdfdUc3bSsVFdWFSg5kSXMAADeNeUQbK0NGmM+JelZhZdIe9hau3/KlQG4JEfbB/Ts/lb9el+r9jb5JElLyrL1mU2LdPMVZVpans2GFQAARERlnWBr7TOSnonGewGYvCNt/frl3hY9s7dFh9sGJElr5uTpi7cu0c1XlKmqKNPhCgEAcCd2jANmmEOtbwffo+0DMka6cm6B/ub9y3Tz8jKV56Y7XSIAAK5HCAZmgOMdA3p6T4ue3tOsw23h4Lu+qkD33nGFbrmiTCU5aU6XCADAjEIIBlyqscevX+5p0S/2NGtfU5+kcPD9uzuu0C3Ly1SSTfAFAOByEYIBF+keDOiXe1v01JtN2n6yR1J4q+K/ft9SvW9lOVMdAACIEkIw4DB/IKjnD7TpqV3Nevlwh4Ihq+rSLP3lzYt1+6pZqizIcLpEAADiDiEYcMBYyOq14116fEejfr2/Vf7AmMpz0/QnG+bpA6srtKSM5cwAAJhOhGAgho629+vxnU168s0mtfiGlZ2WpNtXzdIH1lRofVWBEhIIvgAAxAIhGJhmPv+ontrdpJ/uaNSeRp8SE4yury7WA+9bqhuXliotOdHpEgEA8BxCMDANxkJWvzvaqcd2NOrZ/a0KBENaVp6jL9+2TLevmqXi7FSnSwQAwNMIwUAUner267Htp/TTHY1q9g0rNz1Zd6+fo7vWzdbyilynywMAABGEYGCKAsGQXjjYpke2NWjrkU4ZI21YVKwvMd0BAADXIgQDl+lE56C2vHFKP91xSp0DAc3KTdNnb1ykD9dUalYe6/kCAOBmhGDgEgTHwqO+P3jtpH53tEuJCUablpTo7vVzdF11sRJZ3QEAgBmBEAxMQnvfsB7ZdkqPbGtQa9+wZuWm6fM3VevDV1aqNIftiwEAmGkIwcB5WGv12vFu/fC1k3p2f6uCIavrqov19x9Yrk1LShj1BQBgBiMEAxMMj47pqV1N+o/fnVB9a79y05N137VVuuequaoqynS6PAAAEAWEYCCixTekH7x6Uo9sa1CPf1RLyrL1vz64QrevqlB6Cis8AAAQTwjB8Lxdp3r1/7Ye16/2tSpkrW5aWqr7rp2nq+cXyBimPAAAEI8IwfCkUMjqxfp2fffl49p2olvZqUn642urdO81VaosyHC6PAAAMM0IwfCU4dExPb6zUQ9tfUvHOwdVkZeuL9+2TB+5slJZqVwOAAB4Bf/qwxN8/lF9/9UT+t7vT6hrMKCVs3P19Y+t0a3Ly5SUmOB0eQAAIMYIwYhr7f3DeuiVt/Sj1xo0MBLUpiUl+uR187V+HvN9AQDwMkIw4lJDl1/ffvmYHtvRqOBYSLetnKU/q12gpeU5TpcGAABcgBCMuPJW56C+/psjempXsxKN0QfXzdYnr5vP+r4AAOAshGDEheMdA/rGb47qyV1NSklK0H3vqtInrpvPlsYAAOCcCMGY0Y5Fwu9TkfD7pxvm6xMb5qs4O9Xp0gAAgIsRgjEjNXT59b9fOKwndzUpNSlRn9gwX5+4br6Ksgi/AADg4gjBmFHa+ob19d8c0ZZtp5SYYPSnG+ZrM+EXAABcIkIwZoSewYC+9dIxfe/3JzQWsvro+kp9etMi5vwCAIDLQgiGq/kDQT209S19++XjGgwEdeeaCn32hmrNKWRrYwAAcPkIwXClsZDV4zsb9a/PHVJb34jes6xUX7h5sapLs50uDQAAxAFCMFxn65EO/cMvD6q+tV+rK/P04N1rVVNV4HRZAAAgjhCC4RqHWvv11WcO6qXDHaosSNc37l6j960oZ3tjAAAQdYRgOK7XH9C/PX9YP3jtpLLTkvXX71uqj18zV6lJiU6XBgAA4hQhGI4JWasfvX5S//LsIfmGRvXxq+fqczdVKy8jxenSAABAnCMEwxFvnOjW3/x+WA39+3T1/AJ95f1XaGl5jtNlAQAAjyAEI6ba+4f11V8e1JO7mlWQZvTg3Wv13hVlzPsFAAAxRQhGTIRCVj/e1qD/9et6jYyG9OlNC7U8oVk3ryx3ujQAAOBBhGBMu/rWPn3pZ3u1s6FX18wv1P+8c7kWFGeprq7F6dIAAIBHEYIxbfyBoP7Pi0f00Na3lJOerK99eJXuXFPB1AcAAOA4QjCmxe+Oduq/P75HjT1D+nDNbH3x1qXKz2TVBwAA4A6EYETVwEhQX33moH78eoPmF2XqJ5uv1lXzC50uCwAA4CyEYETNK0fCo7/NviF9YsM8ff49i5WWzIYXAADAfQjBmLL+4VF99Zl6PbItPPr70/uv0bq5BU6XBQAAcF6EYEzJq8e69IXHdqvFN6TN183Xf72pmtFfAADgeoRgXJZAMKR/e+GwvvXSMVUVZuqx+9+ldXPznS4LAABgUgjBuGTHOwb02Z/s0p5Gnz56ZaW+fNsyZabyVwkAAMwcJBdMmrVWj24/pb/5+QGlJCXo3+9Zq1tXsOMbAACYeQjBmBSff1R/9bM9+tW+Vr1rQaH+9cOrVJ6b7nRZAAAAl4UQjIva09irP//RTrX1Deuvbl2izRvmKyGBXd8AAMDMRQjGeVlr9eNtDfrbnx9QUVaKHv3kNVozh5vfAADAzJcwlW82xnzIGLPfGBMyxtREqyg4bygwps8/ulsPPLFPVy8o1NOf2UAABgAAcWOqI8H7JP2BpG9HoRa4xPGOAf3ZD3fqcHu/PnvjIn160yIlMv0BAADEkSmFYGvtQUkyhoAUL369r1VfeGy3khONvnffel1fXex0SQAAAFHHnGBICs//ffC3R/Uvzx3Wqso8ffOetarIY/UHAAAQn4y19sIHGPOCpLJzvPSAtfapyDF1kr5grd1+gffZLGmzJJWWlq7bsmXL5dZ8yQYGBpSVlRWz8800gTGr/9g3oldbxnR1eaL+eHmqUhKnf3SfvrgTfXEn+uJO9MV96Ik7OdmXjRs37rDWvuPetYuOBFtrb4xGAdba70j6jiTV1NTY2traaLztpNTV1SmW55tJ2vuHtfn7O7Srxa+/vHmx/rx2Qcymt9AXd6Iv7kRf3Im+uA89cSc39oXpEB62v9mnT/zndvX4R/Wt/7JWtyxn9zcAAOANU10i7U5jTKOkayT90hjzbHTKwnR7bn+r7vr3V2UlPXb/NQRgAADgKVNdHeIJSU9EqRbEyCPbGvTAE3u1oiJX3723RiU5aU6XBAAAEFNMh/CQM1eAqF1crG/es1YZKfwVAAAA3kMC8ohQyOrvnj6g7/3+hO5cU6F/umulkhOnNBsGAABgxiIEe0AgGNLnH9utX+xu1ic2zNMXb12qBHaAAwAAHkYIjnODI0Hd/8Md2nqkU1+8dYk+ef0Cp0sCAABwHCE4jvmGRvWHD2/T3iaf/vmulfpQTaXTJQEAALgCIThO9Q2P6t6Ht+lAs0//fs9aveeKc236BwAA4E2E4DjUNzyqex86HYDX6cZlpU6XBAAA4CosDxBn+ofDUyD2N/v0TQIwAADAORGC40h/ZArE3kafHrx7rW4iAAMAAJwTIThOnB4B3tvo04PMAQYAALggQnAc8AeCuu8/3tCeRp++cfda3UwABgAAuCBC8AwXHAvp0z9+UzsbevT1j63RLcsJwAAAABfD6hAzmLVWDzyxTy/Wt+sf7lyuW1eUO10SAADAjMBI8Az2by8c0U+2n9JnNi3UPVfNdbocAACAGYMQPEP96PWT+r8vHtFHair1uZuqnS4HAABgRiEEz0DP7W/Vl5/cp01LSvQPdy6XMcbpkgAAAGYUQvAMs+Nktz79yJtaMTtP37h7jZISaSEAAMClIkHNICe7BvUn/7lds/LS9fAf1igjhfsaAQAALgcheIbwB4L65A92SJL+8771KsxKdbgiAACAmYsQPANYa/XffrpHh9v69fWPrdGcwgynSwIAAJjRCMEzwHe3HtfTe1r0lzcv0YZFxU6XAwAAMOMRgl3ud0c79Y+/qtd7V5Tp/uvnO10OAABAXCAEu1hjj1+f+vFOLSjO0j/dtYql0AAAAKKEEOxSw6Njuv+HOxQcs/r2x9cpK5WVIAAAAKKFZOVC1lo98MQ+7Wvq00N/WKP5xVlOlwQAABBXGAl2oce2N+rxnY36ixsW6YalpU6XAwAAEHcIwS7T0OXX3/5iv66ZX6i/uGGR0+UAAADEJUKwi4yFrD7/2C4lGKN/+fAqJSRwIxwAAMB0YE6wi3x363G9caJH//qhVarIS3e6HAAAgLjFSLBLHGzp09eeO6xbrijTH6ytcLocAACAuEYIdoGR4Jg+95NdyklP1lf/YAXrAQMAAEwzpkO4wNeeP6z61n49/Ec1KshMcbocAACAuMdIsMPeONGt77x8XB9bX6lNS1gODQAAIBYIwQ4aGAnqvz66S5X5Gfrr9y1zuhwAAADPYDqEg/7l2UNq6hnSo5+8RplsiwwAABAzjAQ7pL61Tz947aTuuWquaqoKnC4HAADAUwjBDrDW6itP7VdOWpI+/55qp8sBAADwHEKwA57e06LX3+rWF25erLwMVoMAAACINUJwjPkDQX31mYO6YlaOPnrlHKfLAQAA8CTuxoqxB397VC2+YX39Y2uUmMCmGAAAAE5gJDiGTnQO6rsvv6U711RwMxwAAICDCMEx9PdPH1ByotEXb13idCkAAACeRgiOkd/Wt+vF+nZ95oZFKslJc7ocAAAATyMEx8BIcEx/+4v9ml+cqfuuned0OQAAAJ5HCI6B//jdCZ3o8usr779CKUn8kQMAADiNRDbNBkeC+vZLx3R9dbGury52uhwAAACIEDztfvjaSfX4R/UXNy5yuhQAAABEEIKn0VBgTN/delzvXliktXPynS4HAAAAEYTgafTItgZ1DgT06U0LnS4FAAAAZyAET5Ph0TF9++VjWj+vQFfNL3S6HAAAAJyBEDxNHtvRqLa+EX1mE3OBAQAA3IYQPA0CwZC+VXdMa+bk6dqFjAIDAAC4zZRCsDHmn40x9caYPcaYJ4wxeVGqa0Z78s0mNfUO6TObFskY43Q5AAAAmGCqI8HPS1purV0p6bCkL069pJktOBbSg3VHtaIiV7WLWRcYAADAjaYUgq21z1lrg5GHr0maPfWSZrZf7GnWyS6/PrVpIaPAAAAALhXNOcF/LOlXUXy/GWcsZPWN3xzVkrJs3bS01OlyAAAAcB7GWnvhA4x5QVLZOV56wFr7VOSYByTVSPoDe543NMZslrRZkkpLS9dt2bJlKnVfkoGBAWVlZU37eba1BPXN3SP681WpWl+eNO3nm+li1RdcGvriTvTFneiL+9ATd3KyLxs3btxhra2Z+PxFQ/DFGGP+UNL9km6w1von8z01NTV2+/btUzrvpairq1Ntbe20n+f9X39F/kBQz33ueiUmMBXiYmLVF1wa+uJO9MWd6Iv70BN3crIvxphzhuCprg5xi6T/Lun2yQbgeLW/2ae9TT7de00VARgAAMDlpjon+BuSsiU9b4zZZYz5VhRqmpEefeOUUpISdMfqWU6XAgAAgIuY0sRVa+3CaBUykw2PjunJXc26+Yoy5WWkOF0OAAAALoId46Lg2f2t8g2N6iM1lU6XAgAAgEkgBEfBY9sbNTs/Xe9awBbJAAAAMwEheIpOdfv1ytFOfWhdpRK4IQ4AAGBGIARP0WM7GmWMdFeN5zfLAwAAmDEIwVMwFrL66fZT2rCoWBV56U6XAwAAgEkiBE/BK0c71ewb5oY4AACAGYYQPAWPvnFK+RnJunFZidOlAAAA4BIQgi9T92BAzx1o1QfWVCg1KdHpcgAAAHAJCMGX6ck3mzQ6ZvWRK5kKAQAAMNMQgi+DtVaPbj+lVbNztaQsx+lyAAAAcIkIwZdhT6NP9a39+jCjwAAAADMSIfgy/GT7KaUlJ+j9q2Y5XQoAAAAuAyH4EgXHQnp6d7NuXV6unLRkp8sBAADAZSAEX6KdDb3qGw7qPctKnS4FAAAAl4kQfInqDrUrKcHo2kVFTpcCAACAy0QIvkR1hzq0dm4+UyEAAABmMELwJWjvG9aBlj5dX13sdCkAAACYAkLwJXjpcIckqXYxIRgAAGAmIwRfgrrDHSrJTtWycjbIAAAAmMkIwZMUHAvplSOdur66WMYYp8sBAADAFBCCJ2l3Y698Q6O6nqkQAAAAMx4heJLqDnUowUgbFhKCAQAAZjpC8CTVHerQ2jn5ys1gaTQAAICZjhA8CZ0DI9rb5GNpNAAAgDhBCJ6El8eXRitxuBIAAABEAyF4EuoOdagoK0VXzGJpNAAAgHhACL6IsZDV1iMdum5RsRISWBoNAAAgHhCCL2JPY696/CyNBgAAEE8IwRdxemm06xYRggEAAOIFIfgiXjrcoVWVecrPTHG6FAAAAEQJIfgCugcD2t3Yy9JoAAAAcYYQfAFbj3TIWpZGAwAAiDeE4AuoO9ShgswUrazIdboUAAAARBEh+DxCIauXD3dow6IilkYDAACIM4Tg8zjRNaiuwYCuXVDkdCkAAACIMkLwedS39kuSlrFLHAAAQNwhBJ9HfUufEoy0sCTL6VIAAAAQZYTg8zjQ0q/5xVlKS050uhQAAABEGSH4POpb+7SkLNvpMgAAADANCMHn0Dc8qsaeIS0tZz4wAABAPCIEn8OhyE1xS8sZCQYAAIhHhOBzqG/pkyQtKWMkGAAAIB4Rgs/hYGu/ctOTVZ6b5nQpAAAAmAaE4HOobwnfFGcMO8UBAADEI0LwBKGQVX1rPzfFAQAAxDFC8ASnevzyB8a4KQ4AACCOEYInONgSXhmCm+IAAADiFyF4goMtfTJGqi5lJBgAACBeEYInqG/t07zCTKWnsF0yAABAvCIET8BNcQAAAPGPEHyGgZGgTnb5taSMqRAAAADxjBB8htPbJS9hJBgAACCuTSkEG2P+3hizxxizyxjznDFmVrQKc0J9a3i7ZJZHAwAAiG9THQn+Z2vtSmvtaklPS/ofUy/JOQdb+pSdmqSKvHSnSwEAAMA0mlIIttb2nfEwU5KdWjnOqm/p15JytksGAACId8baqeVWY8w/SLpXkk/SRmttx3mO2yxpsySVlpau27Jly5TOeykGBgaUlZV1wWOstfqzF/y6tiJJH1+WGqPKvG0yfUHs0Rd3oi/uRF/ch564k5N92bhx4w5rbc3E5y8ago0xL0gqO8dLD1hrnzrjuC9KSrPWfuVixdTU1Njt27dfvOooqaurU21t7QWPOdXt14Z/+q2+eucK3X3VnNgU5nGT6Qtij764E31xJ/riPvTEnZzsizHmnCE46WLfaK29cZLn+LGkX0q6aAh2o4Mt4ZkdS7gpDgAAIO5NdXWIRWc8vF1S/dTKcU59a7+MkRazXTIAAEDcu+hI8EX8ozFmsaSQpJOS7p96Sc442NKnuQUZykyd6h8JAAAA3G5Kic9a+8FoFeK0+tZ+LSljkwwAAAAvYMc4Sf5AUCe6BpkPDAAA4BGEYEmH2wZkrbSU7ZIBAAA8gRCst1eGWMp0CAAAAE8gBEuqb+lTZkqiZuezXTIAAIAXEIIlHWzt15LyHCUksF0yAACAF3g+BFtrdbClT0vKuCkOAADAKzwfgpt9w+ofDmoJN8UBAAB4hudDcEOXX5I0rzDT4UoAAAAQK54Pwb3+gCQpPzPZ4UoAAAAQK4TgoVFJUn5GisOVAAAAIFY8H4J7IiPBeRmMBAMAAHiF50Nwr39UKUkJSk9OdLoUAAAAxAgh2B9QfkayjGGNYAAAAK/wfAju8Y8qL535wAAAAF7i+RDs848yHxgAAMBjPB+Ce/wBVoYAAADwGEIwI8EAAACe4+kQbK2VbyigPEaCAQAAPMXTIXgwMKbRMat8RoIBAAA8xdMhuGeQjTIAAAC8yNMh2BfZMpnpEAAAAN7i6RB8estkVocAAADwFo+H4NMjwUyHAAAA8BJPh2CfnznBAAAAXuTpEDw+Esy2yQAAAJ7i8RAcUGZKolKSPP3HAAAA4DmeTn8+/ygrQwAAAHiQp0Nwjz+g/EzmAwMAAHiNx0PwKPOBAQAAPMjTIdg3NMrKEAAAAB7k6RDc4w+wUQYAAIAHeTYEj4UsI8EAAAAe5dkQ3D88KmvF6hAAAAAe5NkQfHqjjHxGggEAADzHwyGYLZMBAAC8yrMh2Hd6y2SmQwAAAHiOZ0Pw6ZFgVocAAADwHg+H4MhIcDrTIQAAALzGsyHY5w/IGCmHEAwAAOA5ng3BPf5R5aYnKzHBOF0KAAAAYsyzIbh3aJSpEAAAAB7l3RDsD7AyBAAAgEd5NgT3+ANslAEAAOBRng3Bvf5RRoIBAAA8yuMhmJFgAAAAL/JkCA4EQxoYCbJRBgAAgEd5MgT7hk5vmcxIMAAAgBd5MgT3RrZMZk4wAACAN3kyBJ/eMpnVIQAAALzJkyF4fCQ4nZFgAAAAL/JoCGZOMAAAgJdFJQQbY75gjLHGmKJovN9064mMBOdnMhIMAADgRVMOwcaYSkk3SWqYejmx0Ts0qqQEo8yURKdLAQAAgAOiMRL8b5L+myQbhfeKiV5/QHkZKTLGOF0KAAAAHDClEGyMuV1Sk7V2d5TqiYmewVFWhgAAAPAwY+2FB3CNMS9IKjvHSw9I+pKk91hrfcaYE5JqrLWd53mfzZI2S1Jpaem6LVu2TKXuSzIwMKCsrKzxx/+4bUghK33pqvSY1YB3mtgXuAN9cSf64k70xX3oiTs52ZeNGzfusNbWTHz+oiH4fIwxKyS9KMkfeWq2pGZJ6621rRf63pqaGrt9+/bLOu/lqKurU21t7fjjW/73y6osyNB3733HnwdiaGJf4A70xZ3oizvRF/ehJ+7kZF+MMecMwUmX+4bW2r2SSs44wQldYCTYTXr8Aa2cnet0GQAAAHCIZ9cJZstkAAAA77rskeCJrLVV0Xqv6TQUGNNIMMRGGQAAAB7muZHg8Y0yGAkGAADwLM+F4PEtk9MZCQYAAPAqD4bg8Egwc4IBAAC8y3MhuCcyEpyfyUgwAACAV3kuBPcORUaC0xkJBgAA8CrvheDTc4JZHQIAAMCzPBeCewYDSk9OVFpyotOlAAAAwCGeC8G9Q6OMAgMAAHic90KwP8DKEAAAAB7nwRA8qnxGggEAADzNcyG4xx9gOgQAAIDHeS4E9/pHmQ4BAADgcZ4KwdZa9Q4xHQIAAMDrPBWC+0eCGgtZNsoAAADwOE+F4N5BNsoAAACA10JwZMvkfOYEAwAAeJqnQnAPWyYDAABAHgvBvf7wSDCrQwAAAHibx0JweCSY1SEAAAC8zVMhuCcyEpybTggGAADwMk+F4F7/qLLTkpSU6KkfGwAAABN4Kg32+gOsDAEAAABvheAe/ygrQwAAAMBbIbjXH2BlCAAAAHgsBA+NsjIEAAAAvBWCewYDymNlCAAAAM/zTAgOjoXUNxxkOgQAAAC8E4L7hoOS2CgDAAAAHgrBPWyZDAAAgAjPhODe8RDMSDAAAIDXeSgEj0piJBgAAAAeCsE9kRDMnGAAAAB4JgRfX12sH/3pVSrLTXO6FAAAADgsyekCYqU4O1XF2alOlwEAAAAX8MxIMAAAAHAaIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xhrbexPakyHpJMxPGWRpM4Yng+TQ1/cib64E31xJ/riPvTEnZzsy1xrbfHEJx0JwbFmjNlura1xug6cjb64E31xJ/riTvTFfeiJO7mxL0yHAAAAgOcQggEAAOA5XgnB33G6AJwTfXEn+uJO9MWd6Iv70BN3cl1fPDEnGAAAADiTV0aCAQAAgHFxFYKNMbcYYw4ZY44aY/7qHK8bY8z/jby+xxiz1ok6vWYSfak1xviMMbsiH//DiTq9xBjzsDGm3Riz7zyvc604YBJ94VqJMWNMpTHmt8aYg8aY/caYvzjHMVwvMTbJvnC9xJgxJs0Ys80YszvSl789xzGuuV6SnDpxtBljEiU9KOkmSY2S3jDG/Nxae+CMw26VtCjycZWkf498xjSZZF8kaau19raYF+hd35P0DUnfP8/rXCvO+J4u3BeJayXWgpI+b63daYzJlrTDGPM8/7Y4bjJ9kbheYm1E0iZr7YAxJlnSK8aYX1lrXzvjGNdcL/E0Erxe0lFr7XFrbUDSFkl3TDjmDknft2GvScozxpTHulCPmUxfEGPW2pcldV/gEK4VB0yiL4gxa22LtXZn5Ot+SQclVUw4jOslxibZF8RY5BoYiDxMjnxMvPnMNddLPIXgCkmnznjcqHdeEJM5BtE12T/zayL/++RXxpgrYlMaLoBrxb24VhxijKmStEbS6xNe4npx0AX6InG9xJwxJtEYs0tSu6TnrbWuvV7iZjqEJHOO5yb+18dkjkF0TebPfKfCWxoOGGPeK+lJhf83CZzDteJOXCsOMcZkSXpc0mettX0TXz7Ht3C9xMBF+sL14gBr7Zik1caYPElPGGOWW2vPvM/BNddLPI0EN0qqPOPxbEnNl3EMouuif+bW2r7T//vEWvuMpGRjTFHsSsQ5cK24ENeKMyJzGx+X9CNr7c/OcQjXiwMu1heuF2dZa3sl1Um6ZcJLrrle4ikEvyFpkTFmnjEmRdJHJf18wjE/l3Rv5M7EqyX5rLUtsS7UYy7aF2NMmTHGRL5er/Dfy66YV4ozca24ENdK7EX+vB+SdNBa+7XzHMb1EmOT6QvXS+wZY4ojI8AyxqRLulFS/YTDXHO9xM10CGtt0BjzKUnPSkqU9LC1dr8x5v7I69+S9Iyk90o6Kskv6T6n6vWKSfblLkl/ZowJShqS9FHLLi7TyhjziKRaSUXGmEZJX1H4BgauFQdNoi9cK7F3raSPS9obmecoSV+SNEfienHQZPrC9RJ75ZL+M7IyVIKkR621T7s1i7FjHAAAADwnnqZDAAAAAJNCCAYAAIDnEIIBAADgOYRgAAAAeA4hGAAAAJ5DCAYAAIDnEIIBAADgOYRgAAAAeM7/BxAjItFufjs+AAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"X = np.linspace(0.01, 3, 100)\n",
"Y = np.log(X)\n",
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.plot(X, Y)\n",
"ax.grid()\n",
"ax.set_title(\"$Y = e^X$\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "66f4222e-720d-49e2-a5cb-7ad2e0761b9e",
"metadata": {},
"source": [
"Because the estimation method is MLE, which means we are not seeking to minimize $RSS$, then it would be meaningless to report $R^2$. So the proper replacement is pseudo-$R^2$ as in report.\n",
"\n",
"Pseudo-$R^2$ is comparing the log-likelihood $\\ln{L}$ of the original model with log-likelihood $\\ln{L_0}$ that would have been obtained with only the intercept in the regression.\n",
"$$\n",
"\\text{pseudo-}R^2=1-\\frac{\\ln{L}}{\\ln{L_0}}\n",
"$$\n",
"Compare with $R^2$ in OLS, conceptually, $\\ln{L}$ is equivalent to $RSS$ and $\\ln{L_0}$ is equivalent to $TSS$ of the linear regression model.\n",
"$$\n",
"R^2 =1-\\frac{RSS}{TSS}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "c1f1b23b-9c3f-437a-b53d-da0b7d97fea6",
"metadata": {},
"source": [
"The results object has those data encapsulated, we just call the properties"
]
},
{
"cell_type": "code",
"execution_count": 147,
"id": "2f016611-443c-42a7-9fe3-81b8bcf4c269",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pseudo-R2: 0.6432266516492773\n"
]
}
],
"source": [
"Pseudo_R2 = 1 - log_reg_house.llf / log_reg_house.llnull\n",
"print(\"Pseudo-R2: {}\".format(Pseudo_R2))"
]
},
{
"cell_type": "markdown",
"id": "01916e25-baf5-4f55-a11d-15601bba1109",
"metadata": {},
"source": [
"The equivalent $F$-test in linear regression model is **likelihood ratio** (LR)\n",
"$$\n",
"2\\ln{(L/L_0)} = 2(\\ln{L}-\\ln{L_0})\n",
"$$\n",
"which follows $\\chi^2$ with $k-1$ degree of freedom, $k-1$ is the number of independent variables."
]
},
{
"cell_type": "code",
"execution_count": 146,
"id": "4df10d4b-61fd-45f4-bf7a-c4ec227099e0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Likelihood Ratio: 26.75104440310262\n"
]
}
],
"source": [
"LR = 2 * (log_reg_house.llf - log_reg_house.llnull)\n",
"print(\"Likelihood Ratio: {}\".format(LR))"
]
},
{
"cell_type": "markdown",
"id": "ba9beee0-3e0e-4710-8531-c31dfdf6891c",
"metadata": {},
"source": [
"Keep in mind, with qualitative response models, goodness of fit is not a main concern. "
]
},
{
"cell_type": "markdown",
"id": "8117ddc1-286f-4955-82e0-87662dcb2531",
"metadata": {},
"source": [
"# The Probit Model"
]
},
{
"cell_type": "markdown",
"id": "5faae836-eaf9-4b19-ae1f-23982f7dd975",
"metadata": {},
"source": [
"The **probit model** provides an alternative but similar approach to qualitative response model. We also start by defining variable $Z_i$\n",
"$$\n",
"Z_i = \\beta_1 + \\beta_2X_2+...+\\beta_k X_k\n",
"$$\n",
"But the sigmoid function is a accumulative standard normal distribution\n",
"$$\n",
"P_i = F(Z_i)\n",
"$$\n",
"The marginal effect of $X_i$ is \n",
"$$\n",
"\\frac{\\partial P_i}{\\partial X_i}=\\frac{dP_i}{dZ_i}\\frac{\\partial Z}{\\partial X_i}=f(Z)\\beta_i\n",
"$$\n",
"where $f(Z) = F'(Z)$, is the PDF function of standard normal distribution\n",
"$$\n",
"f(Z) = \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}Z^2}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 155,
"id": "c2607a1f-a36f-4523-b260-03013db47459",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"