{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparing Effect Estimators on Fast-food Employment Data\n",
"\n",
"We look at the famous Card and Krueger study[1] of the impact of a minimum wage change on employment levels in fast food restaurants near the border of New Jersey and Pennsylvania. The treatment in this case was whether the restaurant was in a state which enacted a minimum wage increase (in this case NJ), the outcome is the total number of employees per branch after the change of the law and many covariates were collected for controlling, as detailed below. \n",
"\n",
"The naive estimation would show that raising the minimum wage causes a decrease in number of employees hired per restaurant. However, once we correct for the covariates we find that in fact there was no such effect and that there may have been an increase in employment following the minimum wage increase.\n",
"\n",
"This notebook demonstrates how `causallib` enables causal effect estimation using different methods easily, by utilizing a consistent API.\n",
"\n",
"[1] : Card, D. & Krueger, A. Minimum Wages and Employment: A Case Study of the Fast Food Industry in New Jersey and Pennsylvania. w4509 http://www.nber.org/papers/w4509.pdf (1993) doi:10.3386/w4509."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data \n",
"\n",
"We begin by obtaining the data from Card's website."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import warnings\n",
"\n",
"\n",
"def load_card_krueger(select_covariates=[\"EMPTOT\", \"WAGE_ST\", \"INCTIME\", \"BK\", \"KFC\", \"ROYS\", \"WENDYS\"], outcome=\"EMPTOT2\", dropna=True, include_managers=False):\n",
" \"\"\"\n",
" Load data on the effect of a change in minimum wage laws on fast food employment\n",
"\n",
" Loads data from the study by Card and Krueger [1]. Data downloads from Card's website.\n",
"\n",
" Args:\n",
" select_covariates (list) : List of covariates to return. The default list is based on the analysis by Imbens and Rubin [2].\n",
" To return all covariates pass an empty list: []\n",
" outcome (str) : Name of column to use as outcome. Defaults to total employees after minimum wage change.\n",
" dropna (bool) : Whether to remove columns with NaN in the select_covariates or outcome columns.\n",
" include_managers (bool) : Whether to include managers in employee count. Card & Krueger include them, Imbens & Rubin\n",
" do not include them. (default: False)\n",
"\n",
" Returns:\n",
" databunch (Bunch) : Bunch with member `x`,`a`, `y` and descriptors corresponding to covariates, \n",
" treatment (New Jersey vs Pennsylvania), outcome, and descriptors as parsed \n",
" from the codebook file distributed by Card.\n",
"\n",
" [1] : Card, D. & Krueger, A. Minimum Wages and Employment: A Case Study of the Fast Food Industry in New Jersey and Pennsylvania. \n",
" w4509 http://www.nber.org/papers/w4509.pdf (1993) doi:10.3386/w4509.\n",
" [2] : Imbens, G.W. & Rubin D.B., Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction, Chapter 18, \n",
" (Cambridge University Press, 2015)\n",
"\n",
"\n",
" \"\"\"\n",
" import requests\n",
" import zipfile\n",
" from io import BytesIO, TextIOWrapper\n",
" from sklearn.utils import Bunch\n",
"\n",
" r = requests.get(\"https://davidcard.berkeley.edu/data_sets/njmin.zip\")\n",
" z = zipfile.ZipFile(BytesIO(r.content))\n",
" df = pd.read_csv(z.open(\"public.dat\"), engine=\"python\", sep=\"\\s+\",\n",
" header=None).applymap(lambda x: pd.to_numeric(x, errors=\"coerce\"))\n",
"\n",
" # Load column names and descriptions from `codebook`\n",
" codebook = [repr(line)\n",
" for line in TextIOWrapper(z.open(\"codebook\"), \"cp437\")]\n",
" # part of the codebook is not relevant\n",
" codes = codebook[7:11] + codebook[13:19] + \\\n",
" codebook[21:38] + codebook[40:59]\n",
" cols = [i.strip(\"'\\\" \").split()[0] for i in codes]\n",
" descriptions = [\" \".join(i.strip(\"'\\\" \").split()[\n",
" 4:]).rstrip('\\\\n') for i in codes]\n",
" column_descriptions = dict(zip(cols, descriptions))\n",
" df.columns = cols\n",
"\n",
" # Calculate total employment following Card and Krueger's method\n",
" if include_managers:\n",
" df = df.assign(EMPTOT=df.EMPFT + 0.5 * df.EMPPT + df.NMGRS)\n",
" column_descriptions[\"EMPTOT\"] = \"Initial total employment (full_time_employees + managers + 0.5*part_time_employees)\"\n",
" df = df.assign(EMPTOT2=df.EMPFT2 + 0.5 * df.EMPPT2 + df.NMGRS2)\n",
" column_descriptions[\"EMPTOT2\"] = \"Final total employment (full_time_employees + managers + 0.5*part_time_employees)\"\n",
" else:\n",
" # Calculate total employment following Imbens and Rubin's method\n",
" df = df.assign(EMPTOT=df.EMPFT + 0.5 * df.EMPPT)\n",
" column_descriptions[\"EMPTOT\"] = \"Initial total employment (full_time_employees + 0.5*part_time_employees)\"\n",
" df = df.assign(EMPTOT2=df.EMPFT2 + 0.5 * df.EMPPT2)\n",
" column_descriptions[\"EMPTOT2\"] = \"Final total employment (full_time_employees + 0.5*part_time_employees)\"\n",
"\n",
" #df = pd.get_dummies(df,prefix=\"name\",drop_first=True,columns=[\"CHAIN\"]).rename(columns=dict(name_1=\"KFC\",name_2=\"ROYS\",name_3=\"WENDYS\"))\n",
" df = pd.get_dummies(df, prefix=\"name\", drop_first=False, columns=[\"CHAIN\"]).rename(\n",
" columns=dict(name_1=\"BK\", name_2=\"KFC\", name_3=\"ROYS\", name_4=\"WENDYS\"))\n",
" column_descriptions.update(dict(BK=\"is branch of Burger King\", KFC=\"is branch of KFC\",\n",
" ROYS=\"is branch of Roy's\", WENDYS=\"is branch of Wendy's\"))\n",
" select_covariates = [i for i in select_covariates if i in df.columns]\n",
" if dropna:\n",
" if select_covariates:\n",
" df = df.dropna(subset=select_covariates + [\"STATE\", outcome])\n",
" else:\n",
" df = df.dropna()\n",
"\n",
" # Divide data into pre-intervention and post-intervention variables\n",
" x = df[[c for c in df.columns if \"2\" not in c]]\n",
"\n",
" y = df[[c for c in df.columns if \"2\" in c]]\n",
" y = df.pop(outcome)\n",
"\n",
" a = x.pop(\"STATE\")\n",
"\n",
" if select_covariates:\n",
" x = x[select_covariates]\n",
" return Bunch(X=x, a=a, y=y, descriptors=pd.Series(column_descriptions))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Card and Krueger studied the impact of an increase in the minimum wage in NJ on levels of employment in NJ and the adjacent area in PA. They included many covariates in the survey. Here we will use the subset used by Imbens and Rubin in their analysis of this data in \"Causal Inference for Statistics, Social, and Biomedical Sciences\"."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
EMPTOT
\n",
"
WAGE_ST
\n",
"
INCTIME
\n",
"
BK
\n",
"
KFC
\n",
"
ROYS
\n",
"
WENDYS
\n",
"
STATE
\n",
"
EMPTOT2
\n",
"
\n",
" \n",
" \n",
"
\n",
"
3
\n",
"
30.0
\n",
"
5.00
\n",
"
26.0
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
18.0
\n",
"
\n",
"
\n",
"
4
\n",
"
19.0
\n",
"
5.50
\n",
"
52.0
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
29.5
\n",
"
\n",
"
\n",
"
6
\n",
"
67.5
\n",
"
5.00
\n",
"
26.0
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
24.0
\n",
"
\n",
"
\n",
"
7
\n",
"
18.5
\n",
"
5.00
\n",
"
52.0
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
30.5
\n",
"
\n",
"
\n",
"
8
\n",
"
6.0
\n",
"
5.25
\n",
"
13.0
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
0
\n",
"
9.0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" EMPTOT WAGE_ST INCTIME BK KFC ROYS WENDYS STATE EMPTOT2\n",
"3 30.0 5.00 26.0 0 0 0 1 0 18.0\n",
"4 19.0 5.50 52.0 0 0 0 1 0 29.5\n",
"6 67.5 5.00 26.0 1 0 0 0 0 24.0\n",
"7 18.5 5.00 52.0 1 0 0 0 0 30.5\n",
"8 6.0 5.25 13.0 0 1 0 0 0 9.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = load_card_krueger()\n",
"X, a, y = data.X, data.a, data.y\n",
"data.X.join(data.a).join(data.y).head()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"EMPTOT Initial total employment (full_time_employees ...\n",
"WAGE_ST starting wage ($/hr)\n",
"INCTIME months to usual first raise\n",
"BK is branch of Burger King\n",
"KFC is branch of KFC\n",
"ROYS is branch of Roy's\n",
"WENDYS is branch of Wendy's\n",
"STATE 1 if NJ; 0 if Pa \n",
"EMPTOT2 Final total employment (full_time_employees + ...\n",
"dtype: object"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data.descriptors.loc[data.X.join(data.a).join(data.y).columns]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After loading the data, we begin by checking the association between treatment and outcome, without any causal inference."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"NJ fast food restaurants had a change of -0.21974 employees compared to PA fast food restaurants.\n"
]
}
],
"source": [
"naive_ate = data.y[data.a == 1].mean() - data.y[data.a == 0].mean()\n",
"print(\n",
" f\"NJ fast food restaurants had a change of {naive_ate:.5f} employees compared to PA fast food restaurants.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Based on this, we would conclude that the increase in minimum wage in fact decreased employment levels in fast food restaurants, because the level of employment in the treated group (NJ) was lower than the untreated group (PA). However, this inference would only be justified if the restaurants in NJ and PA were otherwise equivalent. Let us utilize the provided `calculate_covariate_balance` function to check that now:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"covariate\n",
"ROYS 0.001462\n",
"BK 0.015325\n",
"WAGE_ST 0.038154\n",
"INCTIME 0.080553\n",
"WENDYS 0.127842\n",
"KFC 0.142869\n",
"EMPTOT 0.195605\n",
"Name: unweighted, dtype: float64"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from causallib.evaluation.weight_evaluator import calculate_covariate_balance\n",
"\n",
"calculate_covariate_balance(X, a, w=a*0+1)[\"unweighted\"].sort_values()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
EMPTOT
\n",
"
WAGE_ST
\n",
"
INCTIME
\n",
"
BK
\n",
"
KFC
\n",
"
ROYS
\n",
"
WENDYS
\n",
"
\n",
"
\n",
"
STATE
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
20.172794
\n",
"
4.623382
\n",
"
19.051471
\n",
"
0.426471
\n",
"
0.132353
\n",
"
0.250000
\n",
"
0.191176
\n",
"
\n",
"
\n",
"
1
\n",
"
17.273297
\n",
"
4.604803
\n",
"
17.691756
\n",
"
0.415771
\n",
"
0.207885
\n",
"
0.250896
\n",
"
0.125448
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" EMPTOT WAGE_ST INCTIME BK KFC ROYS WENDYS\n",
"STATE \n",
"0 20.172794 4.623382 19.051471 0.426471 0.132353 0.250000 0.191176\n",
"1 17.273297 4.604803 17.691756 0.415771 0.207885 0.250896 0.125448"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X.groupby(a).mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we see that the untreated state (PA) had quite a bit higher employment levels `EMPTOT` than the treated state (NJ) prior to the intervention. Could it be that the measured decline in employment in NJ is in fact explained by the prior lower number of employees per branch? In order to answer this question we need to estimate the difference in counterfactual outcomes:\n",
"$$ATE = E[Y^{a=1}] - E[Y^{a=0}]$$\n",
"The fundamental problem of causal inference is that for each sample we need _both_ treated and untreated outcomes, but in reality no sample was both treated and untreated. Causal inference has many ways to estimate this quantity and `causallib` implements several of them, using a uniform API. To assess the effects, we will load many estimators and compare their results. Note how the same API calls can be used, even though the underlying estimation may be very different. Here we do not focus on the details of the implementation and how each method may be tweaked. To understand the details, you can reference the documentation for each method. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.linear_model import LogisticRegression, LinearRegression\n",
"from matplotlib import pyplot as plt\n",
"from causallib.estimation import (\n",
" DoublyRobustIpFeature,\n",
" DoublyRobustJoffe,\n",
" DoublyRobustVanilla,\n",
" IPW,\n",
" MarginalOutcomeEstimator,\n",
" Standardization,\n",
" StratifiedStandardization,\n",
" Matching,\n",
" PropensityMatching,\n",
")\n",
"from causallib.preprocessing.transformers import MatchingTransformer\n",
"\n",
"\n",
"def learner(): return LogisticRegression(\n",
" solver=\"liblinear\", class_weight=\"balanced\", max_iter=5000)\n",
"\n",
"\n",
"def makeipw(): return IPW(learner=learner())\n",
"def makestd(): return StratifiedStandardization(learner=LinearRegression())\n",
"\n",
"\n",
"all_estimators = [\n",
" MarginalOutcomeEstimator(learner=LinearRegression()),\n",
" Matching(with_replacement=False),\n",
" PropensityMatching(learner=learner(), caliper=0.1),\n",
" IPW(learner=learner()),\n",
" Standardization(learner=LinearRegression()),\n",
" StratifiedStandardization(learner=LinearRegression()),\n",
" DoublyRobustIpFeature(makestd(), makeipw()),\n",
" DoublyRobustJoffe(makestd(), makeipw()),\n",
" DoublyRobustVanilla(makestd(), makeipw()),\n",
"]\n",
"\n",
"\n",
"def get_all_estimates(X, a, y):\n",
" estimates = {}\n",
" for estimator in all_estimators:\n",
" estimator.fit(X, a, y)\n",
" estimator_name = estimator.__class__.__name__\n",
" estimates[estimator_name] = estimator.estimate_population_outcome(\n",
" X, a, y)\n",
" estimates_df = pd.concat(estimates, axis=1).T\n",
" return estimates_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The estimators in `causallib` share a similar API, all of them implementing `fit` and `estimate_population_outcome`. This allows us to easily try them all to see whether they are in agreement as to the causal effect. There is also a dummy estimator `MarginalOutcomeEstimator` which simply implements the naive estimate shown above. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f, ax = plt.subplots(figsize=(8, 6))\n",
"all_estimates.T.plot(ax=ax)\n",
"ax.set_ylabel(y.name)\n",
"ax.set_xlabel(a.name)\n",
"plt.xticks([0, 1], [\"PA\", \"NJ\"]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The most striking thing that we see in this figure is that all of the estimators except one show an increase in total employment under the intervention. The only one which does not show an increase is `MarginalOutcomeEstimator`. As its name suggests, this is an estimator that assumes marginal exchangeability, ie it does not correct for any confounders.\n",
"\n",
"All of the causal inference methods work by estimating what would have happened if the same sample with the same covariates received a treatment other than the one it actually received. Each of the causal methods shifted the estimate of `EMPTOT2` down for PA and, to a lesser extent, up for NJ. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Reweighting with IPW\n",
"Above we saw that the original sample was imbalanced in the sense that the distribution of covariates (especially `EMPTOT`) was highly unequal between the treated and control groups. Balancing the two groups is one of the central challenges of causal inference and it can be accomplished by either reweighting the samples or by restricting the analysis to a subset of the data which is balanced. The first method, reweighting is accomplished by the IPW estimator:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
abs_smd
\n",
"
weighted
\n",
"
unweighted
\n",
"
\n",
"
\n",
"
covariate
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
EMPTOT
\n",
"
0.000816
\n",
"
0.195605
\n",
"
\n",
"
\n",
"
WAGE_ST
\n",
"
0.031575
\n",
"
0.038154
\n",
"
\n",
"
\n",
"
INCTIME
\n",
"
0.008160
\n",
"
0.080553
\n",
"
\n",
"
\n",
"
BK
\n",
"
0.006652
\n",
"
0.015325
\n",
"
\n",
"
\n",
"
KFC
\n",
"
0.006516
\n",
"
0.142869
\n",
"
\n",
"
\n",
"
ROYS
\n",
"
0.009877
\n",
"
0.001462
\n",
"
\n",
"
\n",
"
WENDYS
\n",
"
0.009438
\n",
"
0.127842
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"abs_smd weighted unweighted\n",
"covariate \n",
"EMPTOT 0.000816 0.195605\n",
"WAGE_ST 0.031575 0.038154\n",
"INCTIME 0.008160 0.080553\n",
"BK 0.006652 0.015325\n",
"KFC 0.006516 0.142869\n",
"ROYS 0.009877 0.001462\n",
"WENDYS 0.009438 0.127842"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ipw = IPW(learner=LogisticRegression(solver=\"liblinear\", class_weight=\"balanced\"))\n",
"ipw.fit(X, a)\n",
"w = ipw.compute_weights(X, a)\n",
"calculate_covariate_balance(X, a, w)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that most of the variables are now balanced, especially `EMPTOT`. This gives us confidence that the estimated outcome is indeed an \"apples to apples\" comparison, after weighting."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Restricting the sample with Matching\n",
"Another way to deal with imbalanced covariates is to only include comparisons between treated and control samples which are very similar. There are a variety of ways in which this can be accomplished. One of the most common is matching by propensity score. Here we will consider the case without replacement, where the data set is restricted to a sequence of pairs where each pair has exactly one treatment and one control, and each sample can appear at most once.\n",
"\n",
"We do not want to match on the covariates, but rather on the propensity score. To do this we use the `PropensityTransformer` object as a `propensity_transform` within the `MatchingTransformer`."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"from causallib.preprocessing.transformers import PropensityTransformer\n",
"propensity_transform = PropensityTransformer(learner=learner(), include_covariates=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The propensity model will be autmatically fit when the `MatchingTransformer` is fit."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Transforming data from 347 samples to 136 samples by matching without replacement.\n"
]
}
],
"source": [
"mt = MatchingTransformer(with_replacement=False, propensity_transform=propensity_transform)\n",
"Xm, am, ym = mt.fit_transform(X, a, y)\n",
"print(f\"Transforming data from {X.shape[0]} samples to {Xm.shape[0]} samples by matching without replacement.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now check the covariate balance of the resulting data directly and also with propensity weights:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
abs_smd
\n",
"
weighted
\n",
"
unweighted
\n",
"
\n",
"
\n",
"
covariate
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
EMPTOT
\n",
"
0.005456
\n",
"
0.039639
\n",
"
\n",
"
\n",
"
WAGE_ST
\n",
"
0.028619
\n",
"
0.030309
\n",
"
\n",
"
\n",
"
INCTIME
\n",
"
0.003309
\n",
"
0.047924
\n",
"
\n",
"
\n",
"
BK
\n",
"
0.001717
\n",
"
0.021077
\n",
"
\n",
"
\n",
"
KFC
\n",
"
0.001242
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
ROYS
\n",
"
0.001137
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
WENDYS
\n",
"
0.004414
\n",
"
0.026073
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"abs_smd weighted unweighted\n",
"covariate \n",
"EMPTOT 0.005456 0.039639\n",
"WAGE_ST 0.028619 0.030309\n",
"INCTIME 0.003309 0.047924\n",
"BK 0.001717 0.021077\n",
"KFC 0.001242 0.000000\n",
"ROYS 0.001137 0.000000\n",
"WENDYS 0.004414 0.026073"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ipw_m = IPW(learner=learner())\n",
"ipw_m.fit(Xm, am)\n",
"wm = ipw_m.compute_weights(Xm, am)\n",
"calculate_covariate_balance(Xm, am, wm)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point we have substantially smaller covariate imbalance, but we are using only ~40% of the data. We can now check the estimates directly on the matched subset."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
0
\n",
"
1
\n",
"
ATE
\n",
"
\n",
" \n",
" \n",
"
\n",
"
MarginalOutcomeEstimator
\n",
"
17.544118
\n",
"
19.084559
\n",
"
1.540441
\n",
"
\n",
"
\n",
"
Matching
\n",
"
17.544118
\n",
"
19.084559
\n",
"
1.540441
\n",
"
\n",
"
\n",
"
PropensityMatching
\n",
"
17.645221
\n",
"
17.856618
\n",
"
0.211397
\n",
"
\n",
"
\n",
"
IPW
\n",
"
17.438962
\n",
"
19.157859
\n",
"
1.718897
\n",
"
\n",
"
\n",
"
Standardization
\n",
"
17.431120
\n",
"
19.197557
\n",
"
1.766437
\n",
"
\n",
"
\n",
"
StratifiedStandardization
\n",
"
17.431684
\n",
"
19.198245
\n",
"
1.766561
\n",
"
\n",
"
\n",
"
DoublyRobustIpFeature
\n",
"
17.500168
\n",
"
18.502095
\n",
"
1.001927
\n",
"
\n",
"
\n",
"
DoublyRobustJoffe
\n",
"
17.429642
\n",
"
19.192642
\n",
"
1.763000
\n",
"
\n",
"
\n",
"
DoublyRobustVanilla
\n",
"
17.429092
\n",
"
19.193340
\n",
"
1.764248
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1 ATE\n",
"MarginalOutcomeEstimator 17.544118 19.084559 1.540441\n",
"Matching 17.544118 19.084559 1.540441\n",
"PropensityMatching 17.645221 17.856618 0.211397\n",
"IPW 17.438962 19.157859 1.718897\n",
"Standardization 17.431120 19.197557 1.766437\n",
"StratifiedStandardization 17.431684 19.198245 1.766561\n",
"DoublyRobustIpFeature 17.500168 18.502095 1.001927\n",
"DoublyRobustJoffe 17.429642 19.192642 1.763000\n",
"DoublyRobustVanilla 17.429092 19.193340 1.764248"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"all_estimates_matched = get_all_estimates(Xm, am, ym)\n",
"all_estimates_matched.assign(ATE=all_estimates_matched[1] - all_estimates_matched[0])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFzCAYAAADfQWsjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACD1ElEQVR4nOzdd1yV5fvA8c/Nkr0EFVQU3IqACqi4R2ZqZpmZLbVpln1bjqYzc2allrlzZI7ShtowNTUn7r1xgROQPc/9+wM4P1QQUOCAXO/Xy5fnPPN+jsh1nnFdl9JaI4QQQoiywczUAxBCCCFE8ZHAL4QQQpQhEviFEEKIMkQCvxBCCFGGSOAXQgghyhAJ/EIIIUQZYmHqARQHNzc3Xb16dVMPQwghhCgWu3fvvq61ds9pXpkI/NWrVyc0NNTUwxBCCCGKhVLqXG7ziuxSv1JqrlLqqlLqULZp/kqpbUqpg0qp35RSjjmsV1UptUEpdUQpdVgp9b9s80YopS4ppfZl/ulSVOMXQgghHkRFeY9/PtD5tmmzgWFa64bASmBwDuulAe9presDzYA3lFL1s82forUOyPyzpgjGLYQQQjywiizwa603AZG3Ta4NbMp8/TfQM4f1IrTWezJfxwJHgcpFNU4hhBCiLCnue/yHgceAVUAvoOrdFlZKVQcaATuyTX5TKfUCEErGlYGoexlIamoqFy9eJCkp6V5WF6JYWVtbU6VKFSwtLU09FCFEKVfcgf9F4Gul1CfAr0BKbgsqpeyBn4C3tdYxmZO/BUYDOvPvyZnbzGn9V4FXAby8vO6Yf/HiRRwcHKhevTpKqXs+ICGKmtaaGzducPHiRby9vU09HCFEKVesefxa62Na605a6ybAEuB0TssppSzJCPqLtdY/Z1v/itY6XWttAGYBwXfZ10ytdaDWOtDd/c6MhqSkJMqXLy9BX5R4SinKly8vV6eEEIWiWAO/UqpC5t9mwMfAjByWUcAc4KjW+ovb5nlke/s4cIj7IEFflBbysyqEKCxFmc63BNgG1FFKXVRKvQT0UUqdAI4B4cC8zGU9lVJZT+i3AJ4H2ueQtjchMxXwANAOeKeoxl8clFI899xzxvdpaWm4u7vTrVu3QttHly5diI6Ovqd1R4wYwaRJk4CMy81jxoyhVq1a1K5dm3bt2nH48OE8t7Fq1SqOHDlyT/u/X/Pnz8fd3Z2AgADjn7uNZezYsbe8DwkJKZRxbNy4ka1btxbKtoQQ4n4V2T1+rXWfXGZ9lcOy4UCXzNdbgBxPb7TWzxfaAEsAOzs7Dh06RGJiIjY2Nvz9999UrlywBIa0tDQsLHL/Z1yzpnAyHqdPn87WrVvZv38/tra2/PXXX3Tv3p3Dhw9jbW2d63qrVq2iW7du1K9fP9dlilLv3r2ZNm1avpYdO3YsH374ofF9YQXrjRs3Ym9vX6AvEnn9uwohxL2SWv0m1qVLF1avXg3AkiVL6NPn/78v7dy5k+bNm9OoUSNCQkI4fvw4kHEm2717d9q3b0+HDh1ISEjgqaeeon79+jz++OM0bdrUWKmwevXqXL9+nbCwMOrVq8crr7xCgwYN6NSpE4mJiQDMmjWLoKAg/P396dmzJwkJCXeMc/z48UybNg1bW1sAOnXqREhICIsXLwbA3t7euOyKFSvo168fW7du5ddff2Xw4MEEBARw+vRpTp06RceOHfH396dx48acPn0arTWDBw/G19eXhg0bsnTpUiAjYLZp04bHHnsMHx8fhg0bxuLFiwkODqZhw4acPp3xiMi1a9fo2bMnQUFBBAUF8d9//931M4+IiKB169YEBATg6+vL5s2bGTZsGImJiQQEBPDss8/eckz5Hcdvv/1G06ZNadSoER07duTKlSuEhYUxY8YMpkyZQkBAAJs3byYsLIz27dvj5+dHhw4dOH/+PAD9+vVjwIABNG3alCFDhuTvB0gIIQpITimAkb8d5kh4TN4LFkB9T0eGP9ogz+WefvppRo0aRbdu3Thw4AAvvvgimzdvBqBu3bps3rwZCwsL1q1bx4cffshPP/0EwJ49ezhw4ACurq5MmjQJFxcXjhw5wqFDhwgICMhxXydPnmTJkiXMmjWLp556ip9++onnnnuOJ554gldeeQWAjz/+mDlz5jBo0CDjejExMcTHx+Pj43PL9gIDA+96uT8kJITu3bvTrVs3nnzySQCaNm3KsGHDePzxx0lKSsJgMPDzzz+zb98+9u/fz/Xr1wkKCqJ169YA7N+/n6NHj+Lq6oqPjw8vv/wyO3fu5KuvvmLq1Kl8+eWX/O9//+Odd96hZcuWnD9/nocffpijR48CsHTpUrZs2WIc07Zt2/jhhx94+OGH+eijj0hPTychIYFWrVoxbdo09u3bl+Ox5GccLVu2ZPv27SilmD17NhMmTGDy5MkMGDAAe3t73n//fQAeffRR+vbtS9++fZk7dy5vvfUWq1atAjKyTbZu3Yq5uXmun6sQQtwPCfwm5ufnR1hYGEuWLKFLl1srEN+8eZO+ffty8uRJlFKkpqYa5z300EO4uroCsGXLFv73v4zKxr6+vvj5+eW4L29vb+OXgiZNmhAWFgbAoUOH+Pjjj4mOjiYuLo6HH364kI8yQ2xsLJcuXeLxxx8HMN4i2LJlC3369MHc3JyKFSvSpk0bdu3ahaOjI0FBQXh4ZDzTWaNGDTp16gRAw4YN2bBhAwDr1q275d59TEwMcXFxQM6X+oOCgnjxxRdJTU2lR48euX5Run2dvMZx8eJFevfuTUREBCkpKbmm3m3bto2ff85IVnn++edvObvv1auXBH0hypir5y5x/Xw49VsFFcv+JPBDvs7Mi1L37t15//332bhxIzdu3DBO/+STT2jXrh0rV64kLCyMtm3bGufZ2dkVeD/lypUzvjY3Nzde6u/Xrx+rVq3C39+f+fPns3HjxlvWc3R0xM7OjjNnztxy1r97927atGkD3PrUeWGmnWUfs5mZmfG9mZkZaWlpABgMBrZv337XZw2ya926NZs2bWL16tX069ePd999lxdeeOG+xzFo0CDeffddunfvzsaNGxkxYkS+jzPLvfy7CiFKr+0r/2TvtqlYOStqN/sZi2Io0iX3+EuAF198keHDh9OwYcNbpt+8edP4sN/8+fNzXb9FixYsW7YMgCNHjnDw4MEC7T82NhYPDw9SU1ON9+xvN3jwYN566y3jl4V169axZcsWnnnmGQAqVqzI0aNHMRgMrFy50rieg4MDsbGxxtdVqlQxXtZOTk42XmZfunQp6enpXLt2jU2bNhEcnGuJhjt06tSJqVOnGt/ndrk+y7lz56hYsSKvvPIKL7/8Mnv27AHA0tLylqsqBZX93+v77783Ts/+GUDGLZAff/wRgMWLF9OqVat73qcQonRKTkxk0SejOXHxc2o9fJSKfldQxRSRJfCXAFWqVOGtt966Y/qQIUP44IMPaNSokfGsMicDBw7k2rVr1K9fn48//pgGDRrg5OSU7/2PHj2apk2b0qJFC+rWrZvjMoMGDSIoKIiGDRtSp04dRo8ezS+//IKNjQ0A48aNo1u3boSEhBgviUPGMwwTJ06kUaNGnD59moULF/L111/j5+dHSEgIly9f5vHHH8fPzw9/f3/at2/PhAkTqFSpUr7H//XXXxMaGoqfnx/169dnxoz/Lw+xdOnSW9L5tm7dysaNG/H396dRo0YsXbrUeJvk1Vdfxc/Pz/hwX0GNGDGCXr160aRJE9zc3IzTH330UVauXGl8uG/q1KnMmzcPPz8/Fi5cyFdf3ZHoIoR4gJ07dILZQ17Czv8XKje4QHycPy1b/oa5efGU5FZa62LZkSkFBgbqrKfcsxw9epR69eqZaESFKz09ndTUVKytrTl9+jQdO3bk+PHjWFlZmXpoohA9SD+zQpRFWmvWz1vG6bB5VGsRBgrKlXuR1q3ex8yscM/DlVK7tdaBOc2Te/wPgISEBNq1a0dqaipaa7755hsJ+kIIUYLERUWz9LPPsKm7nRptw4mPr4Bvgy/x8Wla7GORwP8AcHBw4PYrGkIIIUqGo/+FsmHZ51Rtfxpbh1gSE1rTscNX2No6mmQ8EviFEEKIIpCensbqr+dwPXk5NbudJz3dEgf7j2nfrt8tmVCxGzaQdvkyLn1yK3hbuCTwCyGEEIXsxqXLrBg/HNfgA1TzvEpcbDWCg7/Fw6OOcRlDYiJXJ04k6oclWDdsiHOvXqhiKNUtgV8IIYQoRKG/r2fnv1Oo3vksllbJpKU9TpcuY7G0/P9nr5KOHuXS+4NJOX0a1/79cX/n7WIJ+iCBXwghhCgUqUnJrBg3hRSXP6nd+QLJyQ5Uqvg5DRs+alxGGwxEzpvP1S+/xMLZGa+5c7ArpE6g+SV5/CZ0r2159+3bl6+ue9kb52Q3Y8YMFixYULDBCiGEyNXFY6eZ9f5LWDf4iSoNzxMf14CQ5mtuCfqpV65w/qWXuDpxIg5t2+D96y/FHvRBzvhN6l7b8u7bt4/Q0NA7avvn14ABA+5pPSGEELfSWrNp8UqOHZ+Fz6NhmJkZMDd7hUcfHXJLbn7Mn38R8emn6JQUPMaMxqlnz1se8CtOcsZvYgVty5uSksKnn35qrEi3dOlS4uLi6N+/Pw0bNsTPz8/YwQ/go48+wt/fn2bNmnHlyhUgo8LcpEmTAGjbti1Dhw4lODiY2rVrGzsD3q3VrxBCCEiIieX7oR9yKeFLarY7QWqKIzV85tO27TBj0DfExxP+0Udc+t//sPLywmflzzg/+aTJgj7IGX+GtcPgcsHq2+epUkN4ZFyei91LW95Ro0YRGhpq7Do3dOhQnJycjDX6o6KiAIiPj6dZs2Z89tlnDBkyhFmzZvHxxx/fMYa0tDR27tzJmjVrGDlyJOvWreObb77JV6tfIYQoi07s3M/fi8dQtd0p7BxjSExoSfv2X2Nn9//l0hMPHODS4MGknr9A+QGv4f7GG6hiaMKTFwn8JnavbXmzW7dunbHpC4CLiwsAVlZWxucFmjRpwt9//53j+k888YRxmaxWvflt9SuEEGWJNhhYM30el+N+pFa3c6SnW2Jv9yHt271oPIvX6encmDmTa9OmY1GxAtUWLsA2MMfquSYhgR/ydWZelO6lLW9+WFpaGn8Qzc3Nc230k9Vi9m7LCCFEWRd1+SrLPx+Fc5M9VG9whbg4L4ICv8HT8/97aKRcvET40KEk7t6NY9euVBr+KeaOpqnQlxu5x18CFLQt7+1tXh966CGmT59ufJ91qf9+3G+rXyGEeJDs+2szSya9QuVO/+LieZXU1Ed5pPMftwT9m7/9xtkePUg+fhzPiROoPHlSiQv6IIG/RChoW9527dpx5MgR48N9H3/8MVFRUfj6+uLv78+GDRvue0z32+pXCCEeBGmpKSwbPYn9R4ZT+5HDoMypWHEynR/+EiurjKul6bGxXHp/MOGDh1Cudm28V63E6dFH89iy6UhbXpEjafVb8sjPrBDFK+JUGCu/Gk6lVkdxcrtBXFx9WrWcgavr/6ddJ4SGEj5kKKlXruD2xkDcXn21wBX4rl++QGT4GWo3blNoY5e2vKLApNWvEKIs27L0Nw4f+BafrmcxM0/HTPWnW9cPMDc3B0CnpnJt+nRuzJyFZZUqVP9hMTb+/gXahzYYCF01ldoHxmOrHEn3O4S51OoXpiKtfoUQZVFiXBzLPhuPRfVN1OxwkYSE8tSr+QU1a7Y0LpMSFsalwUNIOngQp55PUPGDDzG3tyvQfi6cOkjMsjcIStnPEUtf7J6cXixBHyTwCyGEEACc3nOIP78fTeW2J7F3uklCfAjt20/Fzs4ZyKjSd/Onn7g89nOUpSWVv/wSx84PF2gfqSnJhC4ZSeMzM3FSluxo8ClBT7yNWeaVhOIggV8IIUSZprXmr+8WceH6Qmp2C8NgsMDOdijt271iTIlOi4ri8qfDif37b2ybNsVz/DgsK1Uq0H6Oh67Has07NDeEsce+NV7PTqOpZ7WiOKS7ksAvhBCizLp57QbLx47CIWAX3iFXiI+rQpMm06lc2de4TPy2bYQPHUZaVBQVBr+Pa//+KLP8J8XFxURxaOFggq+u4LpyYW/IdBp3ei7vFYuIBH4hhBBl0sEN29j8+3iqdThNOZtEUlK60rnzBKysrAEwpKRwbcqXRM6bh5WPD94zvsW6fv0C7WP/+h+puOkjgvUNdrk/Tv3nJ9PIybUoDiffJI/fhMzNzQkICMDX15devXqRkJBg6iEZhYaGGmsLbNy4ka1bt+a5Tr9+/bC1tb2luNDbb7+NUorr16/fdd2xY8fma/srVqy4Y3p4eDhPPvlknusLIQRAWmoqP4+byu59H1HrkUMoMzPc3SbwSOevjUE/+dQpwp7qTeS8eTj3eRrvn1YUKOhfv3yB3ZN74L/pNZLMbDjRdTlN35yHg4mDPkjgNykbGxv27dvHoUOHsLKyYsaMGbfMN2X53MDAQL7++msg/4EfoGbNmvzyyy8AGAwG1q9fn69Ww/kJ/Lnx9PTM8QuBEELc7mrYRWa/NwDtvZCq/udISKhN06ar8ffP6FmitSbyhx842/NJ0q5epco33+AxfDhmNjb52r42GNj181dYzmhKw5jNbKs2AM+hodQNfqgoD6tAJPCXEK1ateLUqVNs3LiRVq1a0b17d+rXr09SUpKx5W6jRo2MVfnmz5/PY489Rtu2balVqxYjR440bmvRokUEBwcTEBDAa6+9Rnp6OgD29vY5tuldvny5sepf69atgYxg361bN8LCwpgxYwZTpkwhICCAzZs34+3tbWwYFBMTc8v7p59+mqVLlxq30aJFCyyypaj06NGDJk2a0KBBA2bOnAnAsGHDSExMJCAggGeffRaABQsW4Ofnh7+/P88//7xx/U2bNhESEoKPj48x2IeFheHr62v8XJ544gk6d+5MrVq1GDJkiHHdOXPmULt2bYKDg3nllVd48803C+XfTghROmz/+Q9WznyNal224eAajeIFunX9jfLlqwCQduMGFwe8zpVRo7ENDsbnl1U4tG+X7+1fOHWQI+PaEHTgUy5ZehPxzD807z8eq3LWRXVI90Tu8QPjd47nWOSxQt1mXde6DA0emq9l09LSWLt2LZ07dwZgz549HDp0CG9vbyZPnoxSioMHD3Ls2DE6derEiRMnANi5cyeHDh3C1taWoKAgunbtip2dHUuXLuW///7D0tKSgQMHsnjxYl544YVc2/SOGjWKP//8k8qVKxMdHX3L2KpXr86AAQOwt7fn/fffB6Bt27asXr2aHj168OOPP/LEE09gmdlqsnbt2vz6669ERUWxZMkSnnvuOdauXWvc3ty5c3F1dSUxMZGgoCB69uzJuHHjmDZtGvv27QPg8OHDjBkzhq1bt+Lm5kZkZKRx/YiICLZs2cKxY8fo3r17jpf49+3bx969eylXrhx16tRh0KBBmJubM3r0aPbs2YODgwPt27fHv4DFNoQQpVNyQiLLx05AVVlPzY4XSUhwpWaNSdSq9f+V8uL+/ZfwDz/CEBtLxY8/xuXZZ4xP9OclNSWZ3T+MpNHZjBS9nb7DCXz8f8WaolcQEvhNKOssFzLO+F966SW2bt1KcHAw3t7eQEZ73EGDBgFQt25dqlWrZgz8Dz30EOXLlwcyWutu2bIFCwsLdu/eTVBQkHEfFSpUAHJv09uiRQv69evHU089ZWzRezcvv/wyEyZMoEePHsybN49Zs2bdMv+JJ57gxx9/ZMeOHXz33Xe3zPv6669ZuXIlABcuXODkyZPGY8iyfv16evXqhZubGwCurv9/T6xHjx6YmZlRv3594xWL23Xo0MHYV6B+/fqcO3eO69ev06ZNG+O2evXqZfwchRAPrnMHj7F61ggqtzuRmZsfTNt203Gwz/hdYEhK4uqEiUT98APl6tTBc95crGvXzvf2s1L0mhnC2OPQGq9nphFsghS9gijSwK+Umgt0A65qrX0zp/kDMwB7IAx4Vmsdk8O6nYGvAHNgttZ6XOZ0b+BHoDywG3hea51yP+PM75l5Ycu6x387O7v8VYC6/duoUgqtNX379uXzzz+/Y/nc2vTOmDGDHTt2sHr1apo0acLu3bvvut8WLVoQFhbGxo0bSU9PN15mz9K7d2+aNGlC3759McuW8rJx40bWrVvHtm3bsLW1pW3btiQlJeXrWLNktRCGjHtxeS0jrYaFKJu01vwzdylh4fMycvO1GTbW79G+3evG34NJx45x6f33STl1Gte+fXF/9x3Msv3+uJvbU/T2tfiWxg89U5SHVGiK+h7/fKDzbdNmA8O01g2BlcDg21dSSpkD04FHgPpAH6VU1uOU44EpWuuaQBTwUtEMvWRo1aoVixcvBuDEiROcP3+eOnXqAPD3338TGRlJYmIiq1atokWLFnTo0IEVK1Zw9epVACIjIzl37txd93H69GmaNm3KqFGjcHd358KFC7fMv70NMMALL7zAM888Q//+/e/YXrVq1fjss88YOHDgLdNv3ryJi4sLtra2HDt2jO3btxvnWVpaGp8TaN++PcuXL+fGjRvGY7hfQUFB/Pvvv0RFRZGWlsZPP/1039sUQpRMcVHRzH3/fa6bT8O7xSmSkirQoMFSQkIGZpwgGQzcmDefsF5PYbgZQ9XZs6n4wbB8B/19//xI3BeBBF9dwS73x7F5ZzcBpSToQxEHfq31JuD239q1gU2Zr/8GeuawajBwSmt9JvNs/kfgMZXxNa09kPUI9/dAj8Ied0kycOBADAYDDRs2pHfv3syfP994RhscHEzPnj3x8/OjZ8+eBAYGUr9+fcaMGUOnTp3w8/PjoYceIiIi4q77GDx4MA0bNsTX15eQkJA77n0/+uijrFy50vhwH8Czzz5LVFQUffr0yXGbr732GjVq1LhlWufOnUlLS6NevXoMGzaMZs2aGee9+uqr+Pn58eyzz9KgQQM++ugj2rRpg7+/P++++26BP7fbVa5cmQ8//JDg4GBatGhB9erVpc2wEA+gI1t2sWDMi1Rs9zduVS+TkvwwnTv/TdUqfgCkXrnChZdf5ur48di1aY33r79g37JFvrZ9/fJ5dk96jIDNr5FkZsuJbitKTIpegWiti/QPUB04lO39VqBH5ut3gdgc1nmSjMv7We+fB6YBbmR8IciaXjX7tnP706RJE327I0eO3DGtNJk3b55+4403TLb/5cuX6+eee85k+78XsbGxWmutU1NTdbdu3fTPP/9s4hEVTGn/mRWiKKWlpupVk77Rcye113/9XVOvXtNQ79277JZlbv71lz4e3FQfDWikI5cu1QaDIV/bNqSn6x0rpuibwz108qeueuvcITo5KbEoDqPQAKE6l5hoiof7XgS+Vkp9AvwK3Nf9+dwopV4FXgXw8vIqil2UWYMGDWLt2rWsWbPG1EMpkBEjRrBu3TqSkpLo1KkTPXr0MPWQhBCF4PrFCFZMHI578/14VbhOXGxtWrSYgZtbxkN2hvh4rowbR/TyFVg3aIDnpImUy3yAOi8ZXfQGEpxygCNWDbF/chrNawcU4dEUPaVzeUCq0HagVHXgd535cN9t82oDi7TWwbdNbw6M0Fo/nPn+g8xZ44BrQCWtddrty+UmMDBQ395i9ujRo9SrV+8ej0qI4ic/s0Lcaddv/7B76xSqtT6LhWUKZmZ9aNP6U2P9kMSDBwl/fzAp589T/pVXcH/zDZSVVZ7bTUlOYs+SUTQ6O5NkZckx38ElOkXvdkqp3VrrwJzmFfsZv1Kqgtb6qlLKDPiYjCf8b7cLqJX5BP8l4GngGa21VkptIONWwI9AX+CXYhq6EEKIEiIlKYkV4yaTXuFvaj10gcQEF3xqTKVO7fYA6PR0bsyazbVp07Bwd8fr+/nYBQfnsdUMGSl6b9PMcK7UpOgVRFGn8y0B2gJuSqmLwHDAXin1RuYiPwPzMpf1JOO+fpfMs/k3gT/JSOebq7U+nLnOUOBHpdQYYC8wpyiPQQghRMly4egpfp8xHI82x3BwjiY+PpB2bb/BwSGjJkjqpUtcGjqUxNDdOHZ5hErDh2Oej4d542KiOLzwfYKu/lTqUvQKokgDv9Y650e+M/Lzb182HOiS7f0a4I6byFrrM2Q89S+EEKIM0VqzceHPnAqbiU/XMNBmWFu/Tft2bxpz82/+vprLI0eCwYDn+HE4du+erwp8+/75kUqbPyJI32CX+xPUf34SAaXtaf18ksp9QgghSrz4mzEs+2wMNvW24tMygvg4Dxo1mkbVqgEApMfGcnn0aGJ+/Q2bRo3wnDgBqypV8tzu9cvnObdoEE3iNhJm5sWJR76jaVDHIj4a05ImPSZkb28PZDSZsbGxISAggPr16zNgwAAMBgOPP/44q1atMi5fp04dxowZY3zfs2dPfv755+IethBCFKvj2/fx/cgXcW/1J25eEaQkP8TDD/9lDPoJu3dz9rEexKxeg9ugN6m2cEGeQV8bDOz86UusZjSjYeyWzC56u6j7gAd9kDP+EqNGjRrs27ePtLQ02rdvb6zEt3XrVnr06MGNGzews7Nj27ZtxnW2bdvG9OnTTThqIYQoOtpgYM20uVxN+ZFaXc6RlmqNq+tnNG70dMb81FSuf/st12d8h2XlylRfvAibzP4nd3Ph5H5ilr9BcMrBByZFryDkjL+EsbCwICQkhFOnThESEsLWrVsB2Lp1K48++ijXrl1Da83Zs2exsbGhUqVKJh6xEEIUvsiIq8x67w0S3Wbj1fgsiQneBAX9agz6KefOEfbcc1z/5lucHnsM75Ur8wz6KclJbJs/jAqLOlA15TQ7fYdTd+i/eJWhoA9yxg/A5bFjST5auG15y9WrS6UPPyzwegkJCfzzzz+MGjWKJk2acOjQIVJSUti6dStt2rThzJkzHD16lL179xISElKoYxZCiJJg7x//sv3fyVTrdBpLyxS0oTdduozCwsICrTU3f17J5c8+Q1lYUHnKFzg+8kie28xK0WtuOMcehzZ4PTP1gUrRKwgJ/CXE6dOnCQgIQCnFY489xiOZP8gNGjRgz549bN++nSFDhnDmzBm2bt3K3r17adEif/WlhRCiNEhNTuHnCV+S4rqGWh0vkJTkTNUqX1G3bsZ99/ToaCKGjyD2zz+xDQ7Gc/w4LD087rrN7Cl615Qr+1rOoHHH3BLOygYJ/HBPZ+aFLese/+1atGjBpk2biI2NxcXFhWbNmjFt2jT27t3La6+9VvwDFUKIIhB+Ioxfpg/Ho/Vh3FyiiI9rRNu23+Lo6A5A/PbthA8dRlpkJBXefw/X/v1ReVTRK0spegUhgb+ECwkJ4b333qNt27YA+Pn5sX37dq5cuYKv7x1VkIUQotTZvORXjh7/Fp+uZwEoV24Q7dv9L6OFbkoKV7/6isi587CqXp3q06dj49vgrtu7fvk85xe9SeO4f8tMil5BSOAv4UJCQjhz5gwffJDRrsDCwoIKFSpQtWpVzMzk2UwhROmVGBfH8s8+x7LmJmq0Dic+viIB/tPw8moMQPLp01waPJjkI0dxfro3FYcMwczWNtftaYOBXSu/ou7BifjqZLZVH0CTZ0ZiVc66uA6pVJDAb0JxcXEAVK9enUOHDuW4TIUKFbi9kdLGjRuLemhCCFGkTu8+xJ+LR1K17Qls7OJITupAp4emYG1th9aa6KVLuTJuPGY2NlT5ZjoO7dvfdXvZU/QOWzXEoYyl6BWEBH4hhBDFRhsM/PHdQiJiF1HrkTDS0svh4jKKJo2fBSDtxg0iPv6EuA0bsGvZEs/Px2Lh7p7r9lKSk9i9ZASNz87GSVmxs+EIAnu8VWq66JmCBH4hhBDF4ua1Gyz9fAQuTUKpVucqcXHeNG82kwoVfACI27SJ8A8/whATQ8UPP8TluWdRd7mleSz0H8qteef/U/Sem0ZwJa/iOpxSSwK/EEKIIndg/Ta2/jUOr46nsbRKxmDoRdcuY7CwsMCQlMTVSZOJWrSIcrVq4TlnDtZ1aue6rbiYKA4veI+gaz9Lit49kMAvhBCiyKSlprJq0nQS7H+h5kPnSU5ypErlydSr9zAAScePE/7++ySfPIXLC89T4b33MCtXLtft7Vu3BI8tHxGkI9lVoScNnp9EgKNLcR3OA0ECvxBCiCJx9exFfv7qEyq2PEhl1yji4/xp2/Y7HB3d0QYDkQsWcG3yF5g5O1F11izsW7XMdVvZU/TOmlXjRJdZNA3sUIxH8+CQwC+EEKLQbftpDQcPTMP7kTOgwMryddo/+h5KKVKvXCXigw+I37oV+/bt8RgzGgvXnAvrZE/Ra6BT2Ob9Ok36jJAUvfsgieAm9tlnn9GgQQP8/PwICAhgx44dfPnllyQkJBTaPqpXr87169fvef358+fz5ptvAjBjxgwWLFhQ4G1s3LjR2HDofrYjhCjZkuMTWPjxcM5GjaNGm+OkpLhQt84iWrV6H6UUsevWcfaxx0jYs4dKI0ZQZfq0XIP++RP7ODKuNcEHR3DByoerz/1D837jJOjfJznjN6Ft27bx+++/s2fPHsqVK8f169dJSUmhd+/ePPfcc9jepVBFUUpPT8c8l1SYAQMG3NM2N27ciL29vbGx0L1uRwhRcp3df5S134+kSptj2DnEkpTYho4dv8LGxgFDQgJXxo0netkyrOvXx3PSRMr5+OS4nZTkJHb/MJzGYXNwlhS9QieB34QiIiJwc3OjXOaDLG5ubnz99deEh4fTrl073Nzc2LBhA6+//jq7du0iMTGRJ598kpEjRwIZZ/J9+/blt99+IzU1leXLl1O3bl1u3LhBnz59uHTpEs2bN7+lAFCPHj24cOECSUlJ/O9//+PVV18FwN7entdee41169Yxffp0Tp48yeeff46zszP+/v7GMY4YMQJ7e3ueeeYZunTpYtzuwYMHOXPmDAcOHGDMmDGkpKRQvnx5Fi9eTGJiIjNmzMDc3JxFixYxdepU/vnnH+zt7Xn//ffZt28fAwYMICEhgRo1ajB37lxcXFxo27YtTZs2ZcOGDURHRzNnzhxatWpVXP88Qoh80lqzbs4Szl+fT81HwkhPt8TJaTgd2r8AQOLBQ4QPHkzKuXOUf+Vl3AcNQllZ5bit7Cl6ux3aUu25qZKiV8gk8AObl53g+oW4Qt2mW1V7Wj2VezoKQKdOnRg1ahS1a9emY8eO9O7dm7feeosvvviCDRs24ObmBmTcDnB1dSU9PZ0OHTpw4MAB/Pz8Mvbj5saePXv45ptvmDRpErNnz2bkyJG0bNmSTz/9lNWrVzNnzhzjPufOnYurqyuJiYkEBQXRs2dPypcvT3x8PE2bNmXy5MlERETwzDPPsHv3bpycnGjXrh2NGjW6Zeyenp7GpkLTp0/n33//pVq1ajg6OrJ9+3aUUsyePZsJEyYwefJkBgwYYAz0AP/8849xWy+88AJTp06lTZs2fPrpp4wcOZIvv/wSgLS0NHbu3MmaNWsYOXIk69atu69/FyFE4Yq9EcWyz0fhELCD6sFXiI+rRtOm31GxYi10ejo3Zs/h2tSpWLi54TV/PnZNg3PcTk4pek0kRa9ISOA3IXt7e3bv3s3mzZvZsGEDvXv3Zty4cXcst2zZMmbOnElaWhoREREcOXLEGPifeOIJAJo0acLPP/8MwKZNm4yvu3btiovL/6e6fP3116xcuRKACxcucPLkScqXL4+5uTk9e/YEYMeOHbRt2xb3zGpZvXv35sSJEzkew3///cesWbPYsmULABcvXqR3795ERESQkpKCt7f3XT+DmzdvEh0dTZs2bQDo27cvvXr1Ms7PfnxhYWF33ZYQongd2bSLf1ePpVq7k1iWS8KQ/gRduozFwsKS1PBwwocMJSE0FIfOnfEYOQJzJ6cctyMpesVLAj/keWZelMzNzWnbti1t27alYcOGfP/997fMP3v2LJMmTWLXrl24uLjQr18/kpKSjPOzLsGbm5uTlpZ2131t3LiRdevWsW3bNmxtbWnbtq1xW9bW1rne189NREQEL730Er/++iv29vYADBo0iHfffZfu3buzceNGRowYUaBt3q4gxyeEKB7paWn8OmUGMeV+plan8yQnO1DZcwL162fc/otZs4aI4SMgPR2Pzz/HqcdjKKXu2M718HOc/2GQpOgVM3mq34SOHz/OyZMnje/37dtHtWrVcHBwIDY2FoCYmBjs7OxwcnLiypUrrF27Ns/ttm7dmh9++AGAtWvXEhUVBWScXbu4uGBra8uxY8fYvn17jus3bdqUf//9lxs3bhifHbhdamoqvXr1Yvz48dSu/f9fnG7evEnlypUBbvkSk/2YsnNycsLFxYXNmzcDsHDhQuPZvxCi5Ll+IZzZ7w+A6t9TxfccCfENaNFiLfXrdyE9Lo7wocO49O57lPPxwXvVSpwf73FH0Dekp7NzxRdYzWxGg9itbK/+BpWH7qSuBP1iIWf8JhQXF8egQYOIjo7GwsKCmjVrMnPmTJYsWULnzp3x9PRkw4YNNGrUiLp161K1alVatGiR53aHDx9Onz59aNCgASEhIXh5ZTwY07lzZ2bMmEG9evWoU6cOzZo1y3F9Dw8PRowYQfPmzXF2diYgIOCOZbZu3UpoaCjDhw9n+PDhAKxZs4YRI0bQq1cvXFxcaN++PWfPZvTXfvTRR3nyySf55ZdfmDp16i3b+v77740P9/n4+DBv3ryCfIxCiGKy67d17Nk1hWqPnMHMzICl5Wt07foe5ubmJOzZS/iQIaSGh+P2xhu4vT4AZXFniDl/Yh9xK97M7KLnh2OvaTSr5W+Coym71O0tXx9EgYGBOjQ09JZpR48epV69eiYakRAFJz+zwlRSkpJYMW4SuvLfVKpxkYQEN3x9v8K7ejN0WhrXv53B9W+/xdLTE88JE7Bt3OjObWRL0UtWVhxrOISgx9+6axMece+UUru11oE5zZMzfiGEELm6cOQkv88aTuW2R7FziCEpsSUdO0zDxsaBlAsXCH9/MIn79+P02GNU/ORjzDOf98nu2K51WK99h+aG85KiVwJI4BdCCHEHrTXrF6zgbMRsanQ5iyHdEkeHj+jQ/kW01kSvXMWV0aPB3JzKX0zGMVtdjyyxNyM5svA9gq6tzEjRa/UdTTo8bYKjEdlJ4BdCCHGL+OgYlo0djV3D//AOvkJ8nFdmbn5t0m/eJGL4CGL/+APboCA8x4/D0tPzjm3s+/sHPP/7iCAdxc4KT+L7/ERJ0SshJPALIYQwOr5tL+tXjsarzUnKWSeSnvYYjzwyDktLK+J37CR86FDSrl/H/d13Kf/Si6jb0oAzUvTepHHcJs6aVSeqyxyaBbY30dGInEjgF0IIQXp6GqunziGSZdR6+BwpyXZ4VJpCgwaPolNSuDp5Mjdmz8GqWjWqL1mCTUPfW9Y3pKcTuvIr6h6aSAOdynbvN2jyzHAsrcqZ6IhEbiTwCyFEGRcZfpXlkz7FvdkeqrrdID6uPq1bz8TZ2YPkM2cJf/99ko4cwfmpp6g4bChmtzUQkxS90kXyKEysMNvyzp8/n/DwcOP7l19+mSNHjgCwfPly6tWrR7t27QgNDeWtt94q0Lbbtm1LVkrk3LlzadiwIX5+fvj6+vLLL7/kuP/7lX2f92Ljxo1069YNgF9//TXHcsh52bdvH2vWrDG+v9ftCFFS7V67keXfvkS1Tv9i7xKNhflLdO26CienSkT9uJSzTzxBang4VaZPw2PUyFuCfkpyEtvmDaXS4g5USTnDzoYjqT/sX6pK0C/R5IzfhO6lLe/dWubOnz8fX19fPDMftJk9e7Zx3pw5c5g1axYtW7YEIDAwx/TOPF28eJHPPvuMPXv24OTkRFxcHNeuXctx/8Xtbp9N9+7d6d69e4G3uW/fPkJDQ42dCO91O0KUNKnJKfw0fgrpldZQs91FEhLK06D2l/j4hJAWGcnFjz8hbv167Fq0wGPsWCwrVrhl/VtS9BzbUe3ZryVFr5QosjN+pdRcpdRVpdShbNMClFLblVL7lFKhSqk72jQppdplzs/6k6SU6pE5b75S6my2eQFFNf7ikFNb3hUrVhjb8rZr1w7IaObz3nvv4e/vz7Zt2xg1ahRBQUH4+vry6quvorVmxYoVhIaG8uyzzxIQEEBiYqLxjHnUqFFs2bKFl156icGDB99yJhwfH8+LL75IcHAwjRo1Mp69JyYm8vTTT1OvXj0ef/xxEhMTAbh69SoODg7G2vz29vZ4e3vnuP+cxgkZZ/JDhw4lODiY2rVrG8v15rZPgNdff53AwEAaNGhgrBQIGa2Jhw4dSuPGjVm+fDl//PEHdevWpXHjxsZGRZDxpeTNN98EICAgwPjHxsaGf//9l507d9K8eXMaNWpESEgIx48fJyUlhU8//ZSlS5cSEBDA0qVLb9lOWFgY7du3x8/Pjw4dOnD+/HkA+vXrx1tvvUVISAg+Pj6sWLGiMH9shLhv4SfCmP3BS9j4LaVSzYskJjSnY4d1+PiEELd5C2cee4z4zZup+MEwqs6aeUvQj70ZyY5p/an9+5NYGxLY3/o7mry3CjcJ+qVGUZ7xzwemAQuyTZsAjNRar1VKdcl83zb7SlrrDUAAgFLKFTgF/JVtkcFa60L9Tbph/kyunjtTmJukQjUf2vV79a7L5Lctb/aWuQD169fn008/BeD555/n999/58knn2TatGlMmjTpjrP5Tz/9lPXr1xvnbdy40Tjvs88+o3379sydO5fo6GiCg4Pp2LEj3333Hba2thw9epQDBw7QuHFjAPz9/alYsSLe3t506NCBJ554wliO9/b9v/nmm3eM89FHHwVybrf77bff5rjPrHHm1pq4fPny7Nmzh6SkJGrVqsX69eupWbMmvXv3zvFzz2on/NtvvzFhwgRCQkJITExk8+bNWFhYsG7dOj788EN++uknRo0aRWhoKNOmTQMyvkBkGTRoEH379qVv377MnTuXt956i1WrVgEZX+q2bNnCsWPH6N69O08++eRdfxaEKA5aazYt+YWTYTOo8cgZDAYLHOyH0aH9KxiSk7k8dixRCxZSrlZNPGfPxrpOnVvWzylFr5Kk6JU6RXbGr7XeBETePhlwzHztBOR1Q/hJYK3WuuA3vEuBrLa8M2fOxN3dnd69e98SWLJkb5kLsGHDBpo2bUrDhg1Zv349hw8fvucx/PXXX4wbN46AgABjt77z58+zadMmnnvuOQD8/PyMQdbc3Jw//viDFStWULt2bd55551cO/DdbZw5tdvNbZ+Q0Zq4cePGNGrUiMOHDxufXQCMAf7YsWN4e3tTq1YtlFLGbeXk5MmTDB48mGXLlmFpacnNmzfp1asXvr6+vPPOO/n6TLdt28YzzzwDZHyxyWpNDNCjRw/MzMyoX78+V65cyXNbQhS1hNg4vv/wQy6nTcS72UmSEj3w9/+J4OBXSDp+grAnexG1YCEuzz9P9eXLbwn618PPsWfiowT89zrxZo6cfPRnmr0xG3sJ+qVScd/jfxv4Uyk1iYwvHSF5LP808MVt0z5TSn0K/AMM01on3++g8jozL0p5teWFW1vmJiUlMXDgQEJDQ6latSojRoy4pU1vQWmt+emnn6hz2zf7u1FKERwcTHBwMA899BD9+/e/I/jnNc6CtNvNqzWxnZ1dvscOGc2RnnrqKWbNmoWHhwcAn3zyCe3atWPlypWEhYXRtm3bAm3zdlnHB1AW+mGIku1U6EH+XjoSrzYnKGeTQFpqNx55ZAIW5pZELljA1UmTMXN0pOrM77Bv3dq4niE9nV0/f0m9w5MkRe8BUtxP9b8OvKO1rgq8A8zJbUGllAfQEPgz2+QPgLpAEOAKDL3L+q9mPkcQmvXwWUmTn7a8t8sKeG5ubsTFxd1y//hu6+Xm4YcfZurUqcbgtHfvXuDW1r6HDh3iwIEDAISHh7Nnz547xnz7/u82ztzkts/8tiauW7cuYWFhnD59GoAlS5bkuNyLL75I//79adWqlXFa9nbC2a+63O0zDQkJ4ccffwRg8eLFt2xPiJJAGwysmT6PzZvfoebD+1DmUKniJB5++CuIusmFV1/jytjPsWvRAp9fVt0S9M8d38exca1pengUF6xqcfX5DTTrN1aC/gOguM/4+wL/y3y9HJh9l2WfAlZqrVOzJmitIzJfJiul5gHv57ay1nomMBMyuvPdz6CLSn7b8mbn7OzMK6+8gq+vL5UqVSIoKMg4r1+/fgwYMAAbGxu2bduWrzF88sknvP322/j5+WEwGPD29ub333/n9ddfp3///tSrV4969erRpEkTAFJTU3n//fcJDw/H2toad3d3ZsyYkeP+cxtnbnLbp7+/f75aE1tbWzNz5ky6du2Kra0trVq1uiNonzt3jhUrVnDixAnmzp0LZGQ/DBkyhL59+zJmzBi6du1qXL5du3bGWyEffPDBLduaOnUq/fv3Z+LEibi7u0s7YVGi3Lx6g2UTPsE1aDde7teJj6tDq1azcHGpTOz69UR89DGGxEQqjRiOc+/eKKWA/++i1yRsNomqHLv8RhHYY5B00XuAFGlbXqVUdeB3rbVv5vujwOta641KqQ7ABK11k1zW3Q58kPmwX9Y0D611hMr4CZ0CJGmth+U1DmnLKx4E8jMr8mv/uv/YvmE8Xq1OYW6ehpXV87Rq+SEqOZkr4ycQvXQp5erXo/LEiZSrUcO4XlaKXnXDeXY7tKPas1Nxq1TVhEci7pVJ2vIqpZaQ8cS+m1LqIjAceAX4SillASQBr2YuGwgM0Fq/nPm+OlAV+Pe2zS5WSrkDCtgHDCiq8QshRGmTlprCyonTSHb9lZodLpCY4ErtWl9Qo0YrEg8dJnzwYFLCwij/8ku4v/UWysoKuLWL3lVVPiNFr7100XtQFVng11r3yWXWHWf4WutQ4OVs78OAyjksJ50ehBAiB5fPnmfVtE/waHkIV6doEuKb0r79N9haO3B91iyuffU1FuXL4zVvLnbNmhnX2/vXIipv/URS9MoQqdwnhBCl3Nblazh8fCreD59BazPs7QbTvt1rpF2+zPnX/0fCzp04PPwwHiNHYO7sDGR20Vv8Jo3jM7roRXebR7PGbU16HKJ4SOAXQohSKik+nmWfj8Oq5gZ8QiKIj/MkMPBbPD19ifnjDyI+HQ5paXiMHYvT4z1QSt2RorfN5w0C+0iKXlkigV8IIUqhs/uO8OfikVRpfRRr2wRSUx7hkUcmYZacTviwD7i5ahXW/n5UnjABq8yU23PH9xG/4g2aph7isJU/jk9Np3nNhiY+ElHcJPALIUQporXmr1kLCY9bSI2Hw0hNtaFChfH4NexJwt69hA8ZSuqlS7gNfB23119HWVpmpOgt/pQm5+aQqMqx0380QY+9KSl6ZZT8q5uQubk5AQEBNGjQAH9/fyZPnozBYLjn7WU1zrldv3798iygU716dWOr3TZt2nDu3Lm7Lp+fbebH7a1889OKt1+/fnh7exsb7Xz99dcF3m9YWJixWJAQpUXsjSjmDvsfsY7f4BVwhsSEGjRrtpqG9R7j2vTpnHvueTAYqLZoYcZT+5aWHNv5N+HjA2l+/jsOOrYmdcAOgh9/S4J+GSZn/CZkY2NjbBhz9epVnnnmGWJiYhg5cqRJxpPVGGj48OGMGTOGWbNmFfk+77WV78SJE++r8U1W4M+qtZ9fd2v9K0RROrRpB1v+/JxqbU9iYZmCmXqerl0/IT08nHMDXyBx714cuz9KpU8+wdzB4bYUPTdJ0RNG8pWvhKhQoQIzZ85k2rRpaK1JSkqif//+NGzYkEaNGhkr+GVvCwvQrVu3W7rtvfPOOzRo0IAOHTpwe6ni9evX06NHD+P7v//+m8cff/yOsTRv3pxLly4BubeeBVi3bh2BgYHUrl2b33///a7jS09Pp1+/fvj6+tKwYUOmTJmSYyvf7Ozt7e96PNnl1l44LCyMVq1a0bhxYxo3bszWrVsBGDZsGJs3byYgIIApU6bc9XO9vS3yokWLCA4OJiAggNdee4309PRcxyXE/UpPS+PniV+z9/AwanU4SLrBGu/qs2nbdjixv/3G2R6Pk3zqFJ6TJlF5wgTMHRzY+9ciEqc0IejaSnZW7IXDe6H4S9AXmeSMH4j+7TQp4fGFuk0rTzucH62R94LZ+Pj4kJ6eztWrV1m0aBFKKQ4ePMixY8fo1KkTJ06cuOv68fHxBAYGMmXKFEaNGsXIkSON7WQho/zswIEDuXbtmrHE7IsvvnjHdv744w/jF4S7tZ4NCwtj586dnD59mnbt2nHq1Klcx7Zv3z4uXbrEoUOHAIiOjsbZ2TnXVsJ5Hc/gwYMZM2YMAAsXLmTJkiU5theuUKECf//9N9bW1pw8eZI+ffoQGhrKuHHjmDRp0i1fWO72uWa1RT569Cjjx4/nv//+w9LSkoEDB7J48WJeeOGF3P9hhLhHV89dYuXXn1KxxT48naNJiA+kffsZWKeZEf7ee8SsWYtNYBMqjx+PZeXKXAsP48LiQTSO38QZSdETuZDAX0Jt2bKFQYMGARnNZ6pVq5Zn4DczMzO2qH3uueeMrW+zKKV4/vnnWbRoEf3792fbtm0sWLDAOL9du3ZERkZib2/P6NGjgYzWsz///DOQ0Xp2yJAhxuWfeuopzMzMqFWrFj4+Phw7dizXsfn4+HDmzBkGDRpE165d6dSpU56fwd2O5/ZL/f379+fXX39l0qRJAMb2wp6enrz55pvs27cPc3PzPD/DnGRvi/zPP/+we/duY++BxMREKlSoUOBtCpGXHav+4sDBr6j+8ClAYWf7Lu3bDSRh1y7ODB1G2rVruL/9NuVfeRkN7Fg+iXqHJ1Nfp7LN500C+3wqKXoiRxL4ocBn5kXlzJkzmJub3zWQWFhY3PIA4N1a8mY13ciuf//+PProo1hbW9OrVy8sLP7/R2DDhg04Ozvz7LPPMnz4cL744vaOyHffvlIq1/G5uLiwf/9+/vzzT2bMmMGyZcuMTXLyK6fjyZJbe+ERI0ZQsWJF9u/fj8FgwNraOsf17/a5Zm+LrLWmb9++fP755wUauxD5lZKUxPJx4zGrtg6fFuHEx3nQpMk3eLrX5dqUL7kxaxZWXl5UX/IDNg0bSoqeKDC5x19CXLt2jQEDBvDmm2+ilKJVq1YsXrwYgBMnTnD+/Hnq1KlD9erV2bdvHwaDgQsXLrBz507jNgwGg/FJ+x9++IGWLVvesR9PT088PT0ZM2YM/fv3v2O+hYUFX375JQsWLCAyMvKurWeXL1+OwWDg9OnTnDlz5q7ju379OgaDgZ49ezJmzBhja9+7tb3Nz/Fkya298M2bN/Hw8MDMzIyFCxca78ffvt+7fa7ZdejQgRUrVnD16lUAIiMj88yAECK/zh86yZxP+uPU5GfcvcJJTelE585/45bsQFifZ7gxcybOT/bE++efMK9di21zh+DxQwcqp55lp/9o6g/bSFUJ+iIPcsZvQomJiQQEBJCamoqFhQXPP/887777LgADBw7k9ddfp2HDhlhYWDB//nzKlStHixYt8Pb2pn79+tSrV4/GjRsbt2dnZ8fOnTsZM2YMFSpUYOnSpTnu99lnn+XatWu5dnrz8PCgT58+TJ8+/a6tZ728vAgODiYmJoYZM2ZgbW2d6/guXbpE//79jWfVWWfMd2slnN/jgdzbCw8cOJCePXuyYMECOnfujJ2dHQB+fn6Ym5vj7+9Pv379ePvtt3P9XLOrX78+Y8aMoVOnThgMBiwtLZk+fTrVMgukCHEvtNb8M28pFyLnUvPhs6SlWuPu/jl+DXsRvXw5Vz4fh5mVFZWnfo3jQw9xbOffWP/xDs0NF9jt2J5qz35NsHTRE/lUpG15Swppy3urN998k0aNGvHSSy+Zeih3ZW9vT1xcnKmHUWKU5Z/ZB1lc9E2WfT4Sx4CtuFS8RnxcDVq0mIWTmSMRn3xC3Lp/sAtpjsfn40gsZ87RRe8bU/SutB6Lf/unTH0IogQySVteUTI1adIEOzs7Jk+ebOqhCFHmHf1vN/+uHoNX2xNYWqageIauXUeQuG07Zz/4gPToaCoMHYpr3xfYt+4HKm/9hEAdxc6KT9Hw+QlUcnA29SGIUkgCfxmze/duUw8h3+RsXzyo0tPTWP31LKKtllOzwzmSkpyoXm0qNbxacG3CRCK//x6rmjWoOmsmsY7W7J3cncbxmzNT9ObTrHEbUx+CKMUk8AshRDGKDL/Cii8+xr35Xiq7RJEQ35h27WZgGRFJ2FO9ST5+HJdnn8Xt3XfYvfob6h+eTH3SJEVPFBoJ/EIIUUx2r97A7tDJVHvoFCiwsXmLdm0HEb34By5OnIiZgwNVv5vBjUpOnPzyIZqmHuaQdQDOvabTvKavqYcvHhAS+IUQooilJiWzYsJkdJU/qNHqEgnxFWnU+BsqWVXm4msDiN+8Gfs2bSg//GP2//EVTdbPlS56oshI4BdCiCJ06fgZfpv9KZVbHcLGLpaUpA506vQlKVt3cOajNzDEx1Px00+44lMew+zONDdcINSxA9Wf/UpS9ESRkK+RJiRtee+tLe/y5cupV68e7dq1A6BPnz74+fkxZcqU+x6PEIVFa82GRT/x5y+vU+PhHVhYpeFWfgwPt/uKyHETufj6QCwqVsR9/mxO3viTen88RTlDEvtbzyLwvZ9xk6Avioic8ZuQtOW9t7a8c+bMYdasWbRs2ZLLly+za9euuzYIEqK4JcTEsmzcGOx8N1Mt8Arxcd6EhMzC/moiZ3s+ScqZM7i++CKXGrhh/0svAnUU2yv2lhQ9USzkjL+EkLa8d7blXbJkCQ0bNsTX15ehQ4cCMGrUKLZs2cJLL73E4MGD6dSpE5cuXSIgIIDNmzdz+vRpOnfuTJMmTWjVqtVdGwcJURRO7jrAgnH9qNByNc4VroF+ii6PrEH9sp6zvZ/GEBeH48SxhKuNNN71NnFmjpx+7BeaDZyJnQR9UQzkjB9Yu3Ytly9fLtRtVqpUiUceeaRA60hb3v8XHh7O0KFD2b17Ny4uLnTq1IlVq1bx6aefsn79euM6b7zxBt26dTNeOenQoQMzZsygVq1a7Nixg4EDB7J+/fq8Pnoh7ps2GFg9fS43WELNDudITnKgmtcUvB0bcumV10jYvh37hx4iook7lUJfp7yk6AkTkTP+EmrLli0899xzwL235d2yZcst87O35Y2Ojmbbtm23fDlp164dlStXZu3atfTp0wfIaMv7zDPPABltebNv817b8v7xxx84Ojre9Vh27dpF27ZtcXd3x8LCgmeffZZNmzbddZ24uDi2bt1Kr169CAgI4LXXXiMiIuKu6whRGKIuX2P20IGkenxHlQZhJMb70abN31Q6Y+DMYz1IPHAAqzdfIa7iLoLDJnHOug7Xn99I876fSdAXxU7O+KHAZ+ZFRdry3h+DwYCzs7Px7F+I4rDv7y3s/G881TqeBKWxLvcGbYJe5erYz7n588+U823AtWbuNLr2WWaK3hiCHntDUvSEychPXgkhbXlvFRwczL///sv169dJT09nyZIltGlz9zKljo6OeHt7s3z5ciDjqer9+/ffdR0h7lVaagrLxk7kyLmPqNH6CCkpLtSv9wON7TsQ1rMnN1etgu6dMPM9StO4Hzjg2Ia013cS/PggCfrCpOSM34SkLe+dbXnT0tIoV64cHh4ejBs3jnbt2qG1pmvXrjz22GN5fqaLFy/m9ddfZ8yYMaSmpvL000/j7++fj38NIfLv8unzrPruIzxbHMLOIYbkpDY81P5LYr9fRNj0bzB3dyf+sToElpvPZdzZ32Y2ge16mXrYQgDSltdEIzKtktqWNzk5mZo1a3Lo0CGcnJxMPZwSpyz/zJYkm5f9xomz06na+Azp6ZZUqvgBDdzbEj5kKIl79pDWpAHu1XbjbhnFrswuevK0vihu0pZXGJXUtryhoaE8//zzDBw4UIK+KJGS4uNZ9vlnWNfdSPWgK8THVSMkZBaWWw9z9uUeaDQpbdzx9/ib0+bexHSVLnqiZJLAX8aU1La8gYGBHD161NTDECJHZ/Ye4a9lw6na6ihWVkloQ086txrGtTFjubZ6NWleFanidwIbuxS2+bxF4NMfy9P6osSSwC+EELnQWvPndwu4krKQmh3DSE62x6vqNKrGunGu55OkXb4CAdb41t7LYRvpoidKBwn8QgiRg5vXI1k+8RNcA3dRpfwNEuJ8adPqW5LmL+PcrFmkO1rj1f4aBjdLdkmKnihFJPALIcRtDm3cwX8bP6Na+5OYmRkoZzWAFr49iXj1bZIOHsS8hqJOwGn2uLbH+7mvCa5YxdRDFiLfijTwK6XmAt2Aq1pr38xpAcAMwBpIAwZqrXfmsG46cDDz7XmtdffM6d7Aj0B5YDfwvNY6pSiPQwhRNqSnpbHqi2kkuaykZuuLJMS74e8/Fced5wl7sycGnUqVFpHEV3XgYJtZkqInSqWivi41H+h827QJwEitdQDwaeb7nCRqrQMy/3TPNn08MEVrXROIAkpWTloBSFveW9vyjhw5kg8++OCW+fv27bunFLbw8HCefPJJADZu3Ei3bt2M+8veREiILFfPXWLWB69gWWcBlWpcJCmxFe2DfsJs4kIiPv4Ec8d4fDqFcySwB47vheIvQV+UUkUa+LXWm4DI2ycDWYXanYBw8kll1IhtD2RFnO+BHvc3StPJast7+PBh/v77b9auXWuylryQUbL3wIEDtG3bljFjxhTLPrMH/j59+txRdOjHH3809g0oCE9Pz0L5YiLKhm0r/+S3Ja/g0+k/rKyTcHb+hLYur3Kh59PE/LOOCgE3UR1tON9nJc1e/07y8kWhMSSkErv5EjHr7n6yVZhM8STK28BEpdQFYBLwQS7LWSulQpVS25VSPTKnlQeitdZpme8vApVzWlkp9Wrm+qG3t6ctiaQtbwBVq1bFxcWFHTt2GNdftmwZffr0YdasWQQFBeHv70/Pnj1JSEgAMq48vPXWW4SEhODj42MM9mFhYfj63v3p6t9++42mTZvSqFEjOnbsyJUrV+66vHjwJCcmsmjEcM7HjaJ68HGSEj1p1HAFlf+6zPkXX8I89SqVO0ZysstLVBu2k9qSly8KScrFWCKXnyB87E5urj5D8rkYiqugnike7nsdeEdr/ZNS6ilgDtAxh+Wqaa0vKaV8gPVKqYPAzfzuRGs9E5gJGZX77rbsiROjiY0r3BxyB/t61K79SYHWkba8GWf9P/74I02bNmX79u24urpSq1YtXF1deeWVVwD4+OOPmTNnDoMGDQIgIiKCLVu2cOzYMbp37268xJ+Xli1bsn37dpRSzJ49mwkTJpS4wkai6Jw7dJy1iz+lassjlLNOxJDWg461X+L8i++gw87jUjOe602rE/fsSpr5NDD1cMUDQKemk7D/OnHbw0m9GIeyMsOuSQXsmnpg5ZnzrdqiYIrA3xf4X+br5cDsnBbSWl/K/PuMUmoj0Aj4CXBWSllknvVXAS4V+YhNYMuWLcbAdq9teZ944olb5mdvy9u/f3+2bdvGggULjPPbtWtHZGQk9vb2jB49Gshoy/vzzz8DGW15hwwZYlz+Xtvydu3alU6dOuW4XO/evQkJCWHy5Mm3XOY/dOgQH3/8MdHR0cTFxfHwww8b1+nRowdmZmbUr1+/QGftFy9epHfv3kRERJCSkoK3t3e+1xWll9aaf+Yu4WLcPGp2PEtKsh2eHl/itvM6Ye/2wtw8DedWSVzoOoTA7gMlRU/ct7TricTtiCA+9Ao6MQ2LCjY4d6+BbeMKmFkXfxg2ReAPB9oAG8m4X3/y9gWUUi5AgtY6WSnlBrQAJmittVJqA/AkGU/29wV+ud8BFfTMvKhIW16oWrUq3t7e/Pvvv/z0009s27YNyLikv2rVKvz9/Zk/f/4ttzfKlfv/CmkFuVQ2aNAg3n33Xbp3787GjRsZMWJEvtcVpVNcVDRLx4/AufE2qnpfJyGuPi38xnNl6Cfc2H8Ee48krnbww+2V6QRJip64D9qgSToWSdz2CJJPRIGZwqZBeeyaeVDOxynH38/FpUi/yiqllgDbgDpKqYtKqZeAV4DJSqn9wFjg1cxlA5VSWWf/9YDQzGU2AOO01kcy5w0F3lVKnSLjnv+cojyG4iJtef9fnz59eOedd/Dx8aFKlYxfvrGxsXh4eJCammr8XO7XzZs3qVw54xGR77//vlC2KUquo/+FsvjLflRu+xcOLlFYWbxCkNnLRDzeG3XoMDaN04n6ZDyBH/9CeQn64h6lx6YQs+E8lyfs4saCI6Rejsexoxcew4Io/2w9rGs4mzToQxGf8Wutc3scu0kOy4YCL2e+3go0zGWbZ4DgwhqjKUlb3jvb8trY2NCrVy/eeustpk6datzX6NGjadq0Ke7u7jRt2vSOLwv3YsSIEfTq1QsXFxfat2/P2bNn73ubouRJT0/jt6++I85uOTXbXCAxoTy+NScTO34GUZu/x9oplRs9m+P19jSq20uDKFFwWmtSwmKI2x5B4qHrkK4pV9MZ524+WNcrjzI3baC/nbTlLYNKaltecXdl+Wf2Xt24dJmfvv6Iis33Yu94k6SEEOrb9ibhww8xi07FrK456sOvqB3cwdRDFaWQITmNhL1XidsWQdqVBJS1OXZNKmLXzANLd1uTjk3a8gqjktqWV4jCtmv1P+w/MJnqHU5j0OY42g/D9ZfNJP39PpZWBmL6tMX/g6nSRU8UWOrleOK2R5Cw5yo6JR1LTztcetbCxt8dMytzUw8vTxL4y5iS2pZXiMKSmpTMiokTMKv2J95NI4iPq4yXzSuUGzIBLqehvayxGf8ttRs1N/VQRSmi0wwkHr5B3PZwUs7GgIXC1s8du2YeWFV1MPl9+4KQwC+EeGBcPHaa3+d/QpWWh7C2SSAtpQsVt9zA5o8x6HRF4hPt8R89FTPzkn9WJkqGtOhk4ndGEL/zMoa4VMxdrXF6xBvbwIqY21maenj3pEwHfq11qfqWJsqusvAszv3QWrNh4QrORc6iRsczpKbYYh7/Il7fLST9NKS621LhyxlUaNLU1EMVpYA2aJJPRxO3LYKkozcAsK7jil1zD6xruaDMSnfcKLOB39ramhs3blC+fHkJ/qJE01pz48YNrK2tTT2UEikhJpal40bi6L8ZL7/rxMXUpvw+ayr9Po/UOHPSO7ehwcRpKMvSeXYmio8hIZX43VeJ3xFB2vVEzOwscGhTBbtgDyxcH5z/f3cN/EqphsAsMurhrwWGaq2jMuft1FqX2rS6KlWqcPHixTvq2QtREllbWxtrGoj/d2Lnftb/PpJqrY9hbp5G3LV2+K3YQNIhS5Lt7fCY+SWurduaepiihEu5GEvctggS9l+DNANW1Rxx7eCFTUM3lMWDV7kxrzP+b4ERwHYycuy3KKW6a61PA6X667OlpaWUaBWilNIGA79Pn0205Y/UanuOxAQXrPe5UXv1ehKvlUM1D6LuV9Mwd3TMe2OiTCopdfNNIa/A76C1/iPz9SSl1G7gD6XU82S01xVCiGIVdfkqy6d8TMVmu/F0iib6ah2C1h0lYUcsCcqOSmNH4vL443ILT+SopNXNN4U8j1Ip5aS1vgmgtd6glOpJRrMc16IenBBCZLf3r82E7hpPtQ4nQZuReqgazVYfJuacLeb16uAzdRpWcktE3KYk1803hbwC/3gy6uZvz5qgtT6glOoAlIzONkKIB15aago/T/wCg+dqfJqHE3fTHd+t10j/5xI3k+xwH/QGbq+9hrIoG2dsIn/SY1OI33WZ+J2XSY9OxszRCseOXtgFV8LcsewWbsrrf4mv1vqH2ydqrc+T0WxHCCGKVMSpc/wy50OqtDiEtU08CWGVCVp9lugj9lh4euI9eTI2AQGmHqYoIUpb3XxTyCvwdwY+LI6BCCHE7Tb9uIrT4d9So8MZUlPLUWFnOax/v8rNG/Y4P/EEFT/8EHN7O1MPU5QAOdXNt2/mUSLq5pc0eQV+c6WUC5DjVyStdWThD0kIUdYlxsWxbPxo7Br8i1fANeKuVSDgn8vEb7NF2zhR+cuROHbubOphihKgtNfNN4W8An9dYDc5B34N+BT6iIQQZdrpPYdYt3IEXiFHsbBMwexYeWqtvEnCBWvsmgbhOX4clpUqmXqYwoQepLr5ppBX4D+itW5ULCMRQpRpWmv+mDGf6yykRttzJCU4UGUzmP2dhko0p8Lgt3Ht3x9l9uAVVBH58yDWzTcFeQRWCGFyN6/dYPkXn+AWtBNP5ygSL1WkyppULPYmY+VdCc9JE7Fp0MDUwxQmoA2a5FPRxG1/MOvmm0Jegf+rrBdKKXsArXVckY5ICFGmHFy/je1bx+LV7iQA5nvcqL7aHPOIKJz7PE3FIUMws7Ex8ShFccuom3+F+O0RpN1IemDr5pvCXQO/1nq+UmogMAywA5RSKhYYr7X+pjgGKIR4MKWlprLqi6mkVliFT8glEm664LzJBbd/rmJuZ47HN9/g0L6dqYcpilmOdfM7Vntg6+abQl5Nej4CWgBttdZnMqf5AF8ppVy11mOKYYxCiAfM1XOXWPntB3i2OICzbSwJZytT7S8XLPefwK5VKzzHfoaFu7uphymKSVmum28KeV3qfwHw11onZU3QWp9RSj0F7Ack8AshCuS/FWs4fu5rfDqeJi21HKm7G1J71VV0XBgVPvoIl+eelaeyywipm28aeX2yOnvQzzYxUSllKKIxCSEeQMkJifwwZjj2DTdTrdFVYq9XxCu0ATa/bsGydm0858/HunZtUw9TFDGpm296eQX+S0qpDlrrf7JPVEq1ByKKblhCiAfJ2f1H+HPpJ3i1OoalZQoxZ/zxXZlE2uktuPbti/u772BWruzWTi8LpG5+yZFX4H8L+EUptYWMQj4AgWTc93+sKAcmhCj9tNasnPwVcTa/UbNDGEmJDjiceoSqM/9COztRdfZs7Fu2MPUwRRGRuvklU15P9R9WSvkCzwBZSbSbgNdyugUghBBZrodHsGLCMCq2OkJll0huXq2F7/qKpG9ajV2HDniMGY2Fi4uphymKgNTNL9nyeqr/L611J2BuMY1HCPEA+HveHM5eWIX3IydBgbrcDd+vdmBIuUqlUSNx7tVL7uU+gKRufumQ16V+yacRQuRb5LXL/DL6I8oFnKJGy4vEx7rhcygE84V/YNmgAZ4TJ1LOx9vUwxSFSOrmlz55BX4npdQTuc3UWv9cyOMRQpRC2mBgy9IZHN61gcodT2JrH0tSZDPqzL2B4eyflH/1VdzffANlZWXqoYpCInXzS688Az/Qjdy780ngF6KMu3LxNNu+/ITISpHU6HKatDQr7M72wHPKX5i5u1Pl+/nYBQebepiiEEjd/AdDXoH/nNb6xWIZiRCiVDGkp7P9x/Gc3roH6xbhVKt0hfiYatRYWwmzf9fg2OURKg0fjrmTk6mHKu6T1M1/sOQV+OXrmxDiDmeP7CJi0Ucc1HZU7XoSS6tk0q+2p9YX+zBLj6HS+HE4du8u93dLOamb/2DKK/D3VUr1AGoCB7XWfxb9kIQQJVVSYjx7Fn1Ewv7DXKqrqVHvIMlJDrjsDMFuwWZsAgLwnDgBq6pVTT1UcY+kbv6DL6/AP4CM/P2twGilVLDWenTRD0sIUdIc2bYW87XDOBZVG+e2N6nieoP4qPr4zEvC/OxW3N58E7cBr6EspMZ6aSR188uOvP41W5PRpCddKWULbAbyFfiVUnPJeDDwqtbaN3NaADADsAbSgIFa6523rRcAfAs4AunAZ1rrpZnz5gNtgJuZi/fTWu/Lz3iEEPfmZtR1ji98G6dzh/jXsRpe3Q6hlAGzM+2p9cVWLD098Vy0ENtGjUw9VFFAUje/bMor8KdordMBtNYJqmA/BfOBacCCbNMmACO11muVUl0y37e9bb0E4AWt9UmllCewWyn1p9Y6OnP+YK31igKMQwhxD7TBwN4/v6fq9hFEXPPjZGNnavgcJiHOjcp/VMFm/RacevSg4scfYW4vl4BLE6mbX7blFfjrKqUOZL5WQI3M94qMzn1+ua2otd6klKp++2QyzuQhI1UwPIf1TmR7Ha6UukpGIaHoPMYqhCgkly+cImLJIKpHH2RZYgieD5+kkkMMiVcb4zPlLBb6Ih5fTMaxSxdTD1Xkk9TNF1nyCvz1Cnl/bwN/KqUmAWZAyN0WVkoFA1bA6WyTP1NKfQr8AwzTWifnsu6rwKsAXl5e9z9yIcqA9LQ0QldMxPfol4RF12JlZT982u3DkG6J5c7meM7fjW1QEJ4TxmPp4WHq4Yp8kLr54nZKa110G8844/892z3+r4F/tdY/KaWeAl7VWnfMZV0PYCPQV2u9Pdu0y2R8GZgJnNZaj8prHIGBgTo0NLQQjkiIB9fZI7tI+flNfJJPsCiqK3bNz1Le4zLxN6vg9b055U5dw/1/b1H+xRdR5lJ3vaTLqW6+fXNPqZtfRiildmutA3Oal1eTnlgyLs8bJ2W+z7rU75jjirnrC/wv8/VyYHYu+3UEVgMfZQV9MnYYkfkyWSk1D3i/gPsXQtwmKTGevYs/psmF7zmbVIU5Zp2o+sh+rMolkhoWRI1J+7H2qo7njz9i49sg7w0Kk5G6+SI/8rrU/w9QiYzSvD9qrc/f5/7CyXgqfyPQHjh5+wJKKStgJbDg9of4lFIeWuuIzIcMewCH7nM8QpRph7euwfHv92iuw1ka15PE6heoWX83ycl22P3eAJc/9uPcuzcVhw7BzFYuC5dUUjdfFMRdA7/WuodSygl4ApillLIGlpLxJSDybusqpZaQ8cS+m1LqIjAceAX4SillASSReQ9eKRUIDNBavww8RUYaYXmlVL/MzWWl7S1WSrmTccVhHxl1BoQQBXQz8hrHF75NcNTvnE33ZFbys5Rvthu38teJu16Tal9fxzb9Gh7Tp+HQoYOphytyIHXzxb3K9z1+pZQZ8DTwNTBWa/1FUQ6sMMk9fiEyaIOBPX98T7WdI3DWMfxl9SwX08LxanYUM7N01AE/PGcdwqFFSzw+H4tlhQqmHrK4TU518+2CKkndfHGLe77Hn7lyCNAHaAVsAR7XWm8u3CEKIYra5QunuLzkTZokbOOkWQ3WmQ8gveJf1PA5T0JceSr8YIfTkZNU+PBDXJ57DmUmtdhLEqmbLwpLXg/3hZGRP/8jGZfl0zKnNwbQWu8p2uEJIe5X9hQ9JwxsrPwWJ8OuUiF4OfaON4m/WA/vSWewq1YJz+XfYV2ntqmHLDJJ3XxRFPI64w8j4yn+hzP/ZKfJeEBPCFFCZaXoNU07xgGbJkRWf4PT5+ZRrd0JDOkWmP1Vl1qrTuPywvNUeO89zMpJ1baSQOrmi6KU18N9bYtpHEKIQpSUGM/eRR8ReHEBccqOHf5jOXc8Dp06hepBEcRFe1BlWhIOqTF4zpqFfauWph5ymSd180VxyetS/xCt9YTM17201suzzRurtf6wqAcohCiY7Cl6u5wfxqXtB+z76QuqNN9POesEEo/Wo8b00zi1bY/HmNFYuLqaeshlmtTNF8Xtrk/1K6X2aK0b3/46p/clmTzVL8qC7Cl64aoiN9qNJ+aKJSfDp1LF9zQpybY4LHfBfW8MFYcNw7n3U3IWaSK51c23b+YhdfNFobifp/pVLq9zei+EMIHsKXqNdQzbPZ6lbs8RrJs+Bbv666na8DpxV6rhNeUaTlXKU/nneZTz8TH1sMskqZsvSoK8Ar/O5XVO74UQxezyhVNc/uENmiRu55R5DaK7/0AlVZ5l0wfi1ewIZuZppP5Xi5o/nMP9pZdwf+stlJWVqYdd5uRUN9+lZy2pmy9MIq/A76+UiiHj7N4m8zWZ76VShBAmkp6Wxq7lE2h47Csc0Wyv9Q5NnvqA9d//yDX9CTVanyMhzgW32Wa4xSTjOW8+ds2amnrYZYrUzRclVV5P9ctXUSFKmLOHd5CychDN0o5zwCYQt97TaeBYke9HvYt74HY8HKOJO+ODz9eXcGn/MB4jR2Du5GTqYZcZUjdflHSSECpEKXF7il5o4/E06fYqhzbtZM22wXi1OYHBYIb+xYu6myOpOHIsTj0ekzPLYiB180VpIoFfiFLg9hS9ms99RSMXd36e8iWp7ivxbhpOXFQFPKcm4l6xIp6rvsfKy8vUw37g3Vk33xKHNlWxC64kdfNFiSWBX4gSLCNF738ER63mkqrIwfbzCWr9ONfOhzPnq5fxbLYHJ5sEEvb7UGPOJSq+NhC31wegLOS/dlHKsW7+Q9Ww8ZW6+aLkk98OQpRAGSl686i2c5QxRc//+fFUtnNg+y9/cjRsCt5tT5GaYoPFPDfqRWg8FyzCtnGpKK1RKkndfPGgkMAvRAmTPUXvpHlNorv/QDP/FqQkJbF4zKeUq/UXXg2vEXu5MtW+iqJi2zZUnPEJ5vYSfIqC1M0XDxr5qRWihLgjRa/2uwQ+9QEWllZcOHyKtcs/pGqzQ5ibp5G80Ys6a+PwHDEOp65dTT30B47UzRcPMgn8QpQAOaXoNfOui9aadfN+ICJlDjVah5EY54TjLEUVpyp4/jIeS09PUw/9gSJ180VZIIFfCBPKSNH7kMCLC29J0VNmZsRHx7D0i49xbfQfnk7RxJ6uSo0ZkXi8NojyL72IMpcyG4Uht7r5zt18pG6+eCBJ4BfCRA7/txrHde9npuh1ptbzXxHoVgmAY9v2sGXjSKq2PAbajPSfKtHwlDWe83/ApqGviUf+YJC6+aKsksAvRDG7M0VvAUGtHwMgPT2N36Z9S6LzcrybXiI+0g2P6clUadmBiiuHYWYrAel+Sd18UdZJ4BeimNyeorfN8zkCnhtHZTsHACLDr/LTN8PwaLobd9s44vZUpfbPSVQZNQmHjh1NPPrSTermC/H/JPALUQwunz/J5SVv3pKi19y/hXH+rjXrOXh8ItXbniI1xRrzeW40sqqBx8rPsaxYwYQjL92kbr4Qd5LAL0QRykjRG0/DY1/fkaIHkJqcworJ47D0WUs1/6vERlSi+rcJVH35dVxfeAFlJlXgCkrq5gtxdxL4hSgit6boBeH29HSaVa9jnB9+IozffxhG1aYHsbBMIXGDJ777nKgyew7WdeuacOSlk9TNFyJ/JPALUciyp+jFKjtCm0ygSddXjGfvWms2/vATF2Jm4NPqLEnxjjjOsKZhUGcqLH8PM2sJUgUhdfOFKBgJ/EIUorul6AEkxMSxbPJwnP3/pbJHFLGnKlNzqab68M+wb93ahCMvXaRuvhD3TgK/EIXg5o0rHF/0To4pellOhR5gw1+f4tXiGGhIW1mJJqkN8Fg2Bovy5U008tJF6uYLcf/kf4oQ90EbDOxZO5fqu0bRWMeyzfMFAp4ba0zRy1pm9TeziLX7AZ9mF4mPcqXSLEWtvm/i/PTTkkqWB6mbL0ThksAvxD26fP4kV5a8QZPEHZy0qEV09x9p7hdyyzLRV66zfNowKgWH4m4bS9zeKtTf6orX9EmUq1HDRCMvHaRuvhBFQwK/EAV0Z4reewQ+NcyYopdl71+b2Xvoc6q3OUlaajnMFrjTzLcb7kv+h5mVVS5bL9ukbr4QRU8CvxAFcPbwDlJXvkmztBM5pugBpKWm8NMXkzHz+o1qAVeIjaiAz4/W1Pjoc+yaNTPRyEs2qZsvRPGRwC9EPiQlxGWk6F1alGOKXpbLZ8/z6/fDqNL0AJaWyST8W5nGUQFUXjgKc2dn0wy+BJO6+UIUPwn8QuTh0H+/4bxuMM11BLtcHqHWc1/ekqKXZdOyXwi7Pg2fVmdISnDAdr4jQU+9jdMTj8sDaNlI3XwhTKtIA79Sai7QDbiqtfbNnBYAzACsgTRgoNZ6Zw7r9gU+znw7Rmv9feb0JsB8wAZYA/xPa62L8jhE2XTzxhWOL3yb4Og1XFSVONRhAUGtHrtjueT4BH6cNAJH3w1UrhtJzGlP6m+shPeESVhVq2aCkZdMUjdfiJKhqM/45wPTgAXZpk0ARmqt1yqlumS+b5t9JaWUKzAcCAQ0sFsp9avWOgr4FngF2EFG4O8MrC3awxBlSW4pelWypehlObvvCOt+/4SqzQ6jlCb1Vw9a+vTCbf7rKEsJZlI3X4iSp0gDv9Z6k1Kq+u2TAcfM105AeA6rPgz8rbWOBFBK/Q10VkptBBy11tszpy8AeiCBXxSS/KToQcbT53/MnEuU5UJ8Qi4QH+1C5eXO1Ht/PLZNmphg5CWL1M0XouQyxT3+t4E/lVKTADPgzt+qUBm4kO39xcxplTNf3z79DkqpV4FXAby8vO570OLBlp6Wxq5l4/A7/jWOwPba7xHU+0PMLe78LxJzI4plX31AhcAdVLSPIXZ/ZRpHBFN55nDMHe68KlCWSN18IUo+UwT+14F3tNY/KaWeAuYAHQt7J1rrmcBMgMDAQHkGQOTqzKEdpK26e4pelgMbtrFrz2d4tTpBepolZkur0K7bYJze6VbMoy45pG6+EKWLKQJ/X+B/ma+XA7NzWOYSt973rwJszJxe5bbplwp9hKJMuCNFL3AiTbq8fEeKHmRcEVj51ZcYPFZSvdFlYi+7U2tjdWoPn4xl5RwvOj3wpG6+EKWTKf53hgNtyAjk7YGTOSzzJzBWKeWS+b4T8IHWOlIpFaOUakbGw30vAFOLfsjiQZPfFD2Aq+cusWruMKo03YulVTIJW6oQ4v407t++ijIvW7nmUjdfiNKvqNP5lpBx5u6mlLpIxpP6rwBfKaUsgCQy78MrpQKBAVrrlzMD/GhgV+amRmU96AcM5P/T+dYiD/aJAshvil6WrStXczL8K3xanyE50Q775T60fGMiNn5+xThq05O6+UI8OFRZSIEPDAzUoaGhph6GMCFtMLB77Ry8d43GSceyy/M5Gj3/Oda2Od+DTklKYtmkUdjWXYeT6w1iznoQcL4F1YZ8ipmdXTGP3jRyq5tv38xD6uYLUcIppXZrrQNzmic34sQDLyNFbyCBiTs5YVGb6MeW0bxh7jXzzx8+yR8rP8Ar6DBmykDaX9Xo0H4Yji91KsZRm47UzRfiwSaBXzyw7kjRq/M+QU99kGOKHmSc4f49bxHX9FxqhJwnPtqZqpvr4DtsCpYVKxbv4E1A6uYLUTZI4BcPpOwpevttgqjQZzrNquWcogcQF32TZVM+wq3xVio53CT2UBWa2T1Lxck5P+X/oJC6+UKUPRL4xQMlI0XvAwIvLc4zRS/LkS272LZ9FFVbHseQboH52jp0fnEy1vXqFePIi5fUzRei7JLALx4Yh7b8ivM/Q/KVogeQnp7GL1Onkeq+guqNI4i96k6Ds+2oMXo4ZtYPXllZqZsvhAAJ/OIBcEeKXseFBLXsftd1boRf4aeZQ6kctBuHcokk7KpOu6APcXq6QzGNuvhI3XwhRHYS+EWplT1Fr7GOZVvlF2j03OdUySVFL8uO39dx9OwkfFqdIjnRFseNTWn39ldYuLkV08iLh9TNF0LkRAK/KJUizh3n6o9v5jtFDyA1KZllX4zFptYfVG1wnZhzHjS17IvH6JcfmIfYpG6+ECIvEvhFqZKRovc5fsen4gRsrzOYoKeG5Zqil+XS8TOsXj4Mr8CDmJmnk7a1Lo/0/hLrWrWKZ+BFTOrmCyHyS34jiFLj9MHtpP8yKCNFzzaYCk9Pu2uKHmTk5q9ftJSI5JnUCDlHwk0nfM52ot6QUZhZWRXTyIuG1M0XQtwLCfyixMtK0Qu6tIibyoHQoEk0eeSlPPPrE2LiWDrlQ1wDtuDhcJPYY160bfgJzo+3L6aRFw2pmy+EuB8S+EWJlj1Fb6dLF+o8/yWB5fOuond8+z42bxqOV8hRDOkWWO1oRrcB07Bwcclz3ZIot7r5zt18pG6+EKJAJPCLEin6+mVOLHonM0XPg0MdFxKcR4oeZDzp/+s335DktBTvwHDirrkRYPYSXsNeKZWXvqVuvhCisEngFyWKNhjYvWY2PqGjaaTj2Va5L42eG5tnih5A1OVrrJgxFM+gUMpbJ5B4sDYPPzoVa5+axTDywiV184UQRUUCvygxIs4d5+qSNwhM2pWRotdjGs19m+Zr3d1/beTAsfFUb3mS1GQbXI73oPHAz1GWpaf8rNTNF0IUBwn8wuTuNUUPIC01heWTx2HlswYv32vEXvSkZa0RuHUtPRX4pG6+EKI4SeAXJnUvKXpZLp8+z6+Lh1A1+ADm5mkYDgXSre9MLJycinjU9y/Huvl1XbFv5kE5qZsvhChCEviFSdyeorc7aDKNH3kx3y1wN/74ExdjvqFGyzASYpyoYXiJ2m+9UcSjvn9SN18IkZiWSGRSJJGJkdxIukFkUiRW5lZ08+lWLPuXwC+K3aHNv+C8fgjN9WVjil6TfKToASTGx7N04se4BGzCo0I0cWe86dDhW+y9S3YFPqmbL8SDy6AN3Ey+mRHMkyK5kXjDGNBvJGb+nXSDyMSM+QlpCXdso6ZzTQn84sETff0yJxe9TVD02owUvYcWEdzi0Xyvf3rPYTb89RFeLY5iMJhhfaYL7ftNQeXjWQBTkLr5QpReyenJxkB9I+nGrQH8trP1qKQo0nX6HdswU2a4lHPB1cYVV2tX/Nz9cLV2pbxNecpbl8fV2tX43sW6+GqMlMzfmOKBog0Gdq+ehc/uMQQUMEUPMorXrP5uJnG2i/EOvkTcDTeCqgzH8+UuRTzyeyN184UoebTWxKTEGM+8jQE821l59tdxqXE5bsfGwiYjWFuXx9POk4ZuDW8J4FnzXG1ccbJywtys5KXfym8hUaTuJ0UP4Ob1SJZ9PRiP4FDcbOJJOt2QLr2/x9KxZD3AJ3XzhSh+qempOQbwWwJ5tkvsaTrtjm0oFC7WLsbg3aB8A+MZuvGs3Ob/X9talv7CWRL4RZFIT0tj19Kx+J2YVuAUvSz71//HnoOfUb31CVKTbXC/+Rr+rwwuukHfg9vr5ps7WuH4UDXsgiph7li6mwAJUdy01sSmxt5yif2W17fdM49Nic1xO+XMyxkDdQXbCtRzrXfnWXnm387lnLEwK1uhsGwdrSgWpw9ux/DLmzRLO1ngFD3I+NKwfNI4LLxX49XwKrGXK9O2+Tc4e/sW4ajzT+rmC5F/qYZUopKi7rgvnv2+efaAnmpIzXE7zuWcjcG7jmudHC+vZ11it7WwlatsdyGBXxSapIQ49i78gKDwe0vRA7h67hKr5r9P1ab7sTBPhfNt6f78TMzMTX+fTOrmC5HxxTchLeGWM+9cL7EnRXIz+WaO27E0szQGbVdrV2o517rlknpWEC9vXR5na2cszaSYVWGRwC8Kxf2k6GXZvGwVYVHTqNHqLIlxjtRxH4lPv15FNOL8k7r54kGXZkgjOjk6x9SznJ5iT05PznE7jlaOxkBe07nm/5+J25S/4wzd3tJezspNRAK/uC/R1y9zcuH/CLr5xz2l6AEkJyby47iPcArYhGetKOIv1qbTYwuwdnIvolHnTermi9IuITUhx6Cd01Ps0cnRaPQd27Aws7jlITcfZ587Anj2+ZbmclZeGkjgF/fkflP0spw9cIx/Vn9A1RZHQJthG/0MHV4YXUSjzltadBLxOy4Tv0vq5ouSJd2Qzs2Um3k+uZ4V4BPTEnPcjoOlg/Gp9eqO1WlcofEtl9izB3VHK0f5kvsAksAvCiw87DjXf7z3FD3IuE/4x7ffEW23BO+mF4mPcqep31dUrFWw7RQGqZsvTCUpLSnXJ9dvP0OPTo7GoA13bMNcmRsDtqu1K16OXjk/+JZZJKaceTkTHKkoSSTwi3xLS00hdNk4/E5MwxnYXncIQb2GFihFDyAuKpofJ79PpWa7cLeJIzk8iC5PL8DConjT36RuvihsBm0gJjkmxyfXc3qCPafSrQC2FrbGoF3Fvgr+7v53PPCWdYbuWM4RMyVln0X+SeAX+XL6wFYMvwyiWfop9ts2pWKf6TTzKnh9/IPr/2PX/rFUa3uCtJRyVLL4AN/nXi6CEedO6uaLgsitdGtOl9gLq3SrjYWNCY5UlBUS+MVdJcbHsm/RMILCf7jnFD2A9PQ0fho3FrMaf1DN/wpx16rQruMCHMtXK6KR30rq5osst5duze3J9YKUbvWw88DXzTfnh95KcOlWUTZJ4Be5OrjpF1w23F+KHsD18+GsnPMelZsdxNIyGaK68mivLzEr4JeHeyF188uG1PTUHKu75fg6OZI0g5RuFWVXkf3mU0rNBboBV7XWvpnTlgJZJdycgWitdcBt69UBlmab5AN8qrX+Uik1AngFuJY570Ot9ZqiOoayKnuK3gXleU8peln+++EnTsfMwLv1GZISHKnpNRnvjg8X8ohvpdOz6uaHk3wyWurml0Jaa+JS43IM4PdautXd1p26rnXvPCsvw6VbRdlUlD/l84FpwIKsCVrr3lmvlVKTgTtKOmmtjwMBmcuYA5eAldkWmaK1nlQkIy7j7kjRq9KPRs9+RtUCpugBpCYls2T0EByabKVypUjir9ajU4/FWFsXXXMdY938HZdJvyl180uaVEMq0UnRd3RHu+O9lG4VokgVWeDXWm9SSlXPaZ7K+J/4FNA+j810AE5rrc8V8vDEbQojRS/L+b2H+HPNJ3i1OQqAXfqrdHh6aGEO1yjXuvmPSt38ona30q055ZlL6VYhSgZTXddqBVzRWp/MY7mngSW3TXtTKfUCEAq8p7WOymlFpdSrwKsAXl5e9zncB1dGit7n+J2Yfl8pepARCP78aiqRzj/j0/wC8TfdaNZ0FhU8/Qp93FI3v2jcXro1r0Ix+SndWsO5BkHWQVK6VYgSQml9Z5nGQtt4xhn/71n3+LNN/xY4pbWefJd1rYBwoIHW+krmtIrAdUADowEPrfWLeY0jMDBQh4aG3vNxPKiyUvRqpZ9iv01Gil6le0jRA4iPjGbJxHepGLIHG9tYUqNCeKjHbCwsCrdYSE518+2be0rd/LtISE3Is8VpQUu35hTApXSrECWHUmq31jowp3nFfsavlLIAngCa5LHoI8CerKAPkP21UmoW8HuRDPIBlz1FL1o5sjv4Cxp37l/gFL0sh9b8w85jU6jW4ThpqVZUdB5Nww7PFNp4pW7+rW4v3ZrbU+xSulUIkRNTXOrvCBzTWl/MY7k+3HaZXynlobWOyHz7OHCoCMb3QDu46RdcNwymub7CTteu1Hn+K5q43lszHG0wsGLEp1BvA9UCLhMX6UW7TotwdKxcKGMtS3Xzk9KScrysfnsQj0yMJCo5Skq3CiHuWVGm8y0B2gJuSqmLwHCt9RxyuG+vlPIEZmutu2S+twMeAl67bbMTlFIBZFzqD8thvshFRoreWwTd/JMLypPDD/1AcIuu97y9G6fD+GneMKqEHMLSKhmV+DiPPjHhvnPzH5S6+VmlW40B/G4PviXeyLV0q52lnTF4V7WvKqVbhRD3rUjv8ZcUZfkevzYY2P37TGrs+Qx7HU9olRdo9NxnWNvY3fM2t876nlNpi/GsfZrkRAca1P+S6j5t72ucOdXNtwuqVKLq5qekp9zRQCW3YB6VFEWavrNIzO2lW7PfM5fSrUKIwlKi7vGL4pM9Re+4RR2iH59K8wb33v0uNT6eJSOG4NBsF5VdbpAQ2YBOjy6kXLl7z803Zd38rNKt+WlxeiPxRv5Kt9pL6VYhRMkmgf8BlD1FzwnF9rpDCeo15J5S9LKc27yNvzZOpGrHIyilsbN8kw5PvnNP2yrKuvm3l27N68G3/JRure9aX0q3CiEeGBL4HzCn9v+H/vUtmqWfYp9tMyo9M51mVWve8/a01vw5ehyRVf/Cp8V5EmLdaNp8LhUqNCjwtu6lbn5W6dacAriUbhVCiIKT324PiMT4WPYvHEpgxJJCSdEDiLsUzo9TPqRC6/1UtI8hJaYVj3SbgYVF/u+551Q3v1x9FwyN7YiumMSZpDBuXNqdY+nWrKfYUwwpOW5bSrcKIUTBSeB/AGSl6DXTV9jp2i2ji949puhlOfTDCnaEz8fr4ROkp1tSsfxYfNv3znHZrNKtWUH7RtINYqOisT+sqXrSCbtEa25axbOpyl7WOG0mLP0C7LpzO/kp3Wo8K5fSrUIIcU8k8JdiUdciOLXof4WWogeQnpzMyg8/wNBkJ9UbRxAbXQUXv/c5Z2bG3iMLc32KPSk9CTQ0SKxBt6jWtIhphCUWHLA/ydZah7jkEYmLrStNrUPoIqVbhRDCZCSdrxS61xS9u5VujUyMRB29hPfZFKq2OYaVVSKhlyqzRN9AZwvGOZVurWDuRv0IL6qfLI9NlDm6nMI8wAmX5tWwqeRY1B+HEEKI20g63wMk/Owxri99g8CkUI5Y1Cal2+e4V6vO/uhDRF7O4cn1xHyUbrWw59FtlfCsm0DlTqdJSrTjjFkvvJv6M+kupVtzrJvfXurmCyFESSaBv4TJKt16ewC/Hn+Nkyc2o+NPEelqzrVytYlTKRj2vwP7b91GQUq32t5M4ZfPR2HX9iBOrtdJuNmAjp2/x8bGJcfx6TQDCYevS918IYQopSTwF7GcSrfe8uT6bU+xx6fG57gdGwO4padSzsKBSp6NaORa9b5Lt5796Rd+PbQAr25HMFMGbK3eoH2Pd3IM3mWpbr4QQjzIJPDfA2Pp1jzqrxekdGtDt4Z3lG61w5ora76jdcQqknHgXPAIGnfue18pegCGpCT+/HA4kQ12UaPVORLiyhPcdDYVK/rdstyDUjdfCCHE/5PAX0BT905l5oGZOc4rzNKtBzetpPyGoQTqK+x0KZwUPYCb+w+wYu4E3NsfoZLDTZLjW/JIlxlYZKsJn1PdfIc2VUtU3XwhhBD3RgJ/AQVXCr6lClxhl269I0Wv0xKCQ7rc93a1wcD+L6YSaliHV7fjGNItcC8/Gr/2zxiXMWXdfCGEEMVDAn8BNfVoSlOPe290k5uMFL3vqLFnLAE6nm1VX6TRs2Ooeh9d9LKkXL7Cr59+Qnqbo3h7hBMfU5U27Rbg5ORVpHXzhRBClDwS+EuA7Cl6xy3qEv3ENJrXDyqUbV/+9Xd+3bKQyo8doZx1IqQ9Trfu4zFEphC9+kyB6uYLIYQo/eQ3vAmlpaYQunQsfie/wQnFjnrDCHxy8H110ctiSEjgvw9HcrrGfnwePkVKkh3eVb/FM7UxN+YdMdbNt2lQHrtmHpTzcZJUPCGEKAMk8JvIqf1b4Ne3aJZ+mn12zanUZxpN76OLXnZx+/azcvpkbDudoGr5ayTcrEdT+89J/ymRGzePYO5oheND1bALqoS5o1Wh7FMIIUTpIIG/mN3ZRe/LQknRA9Dp6Zz6cjobYv7Fq+dRzMzTsAx/loCjD5GSHk25ms44P+qDdb3yKHM5uxdCiLJIAn8xykrRa6avsLP8o9R5bkqhpOgBpIaHs27YcCJbn6VGo7MkxZWn8r5BuBpqYNesInbNPLB0v/+sAyGEEKWbBP5ikJGi9xZBN/8q1BS9LNd//ZU/1/2GY899VHKIJPVCMAE338L5EW+pmy+EEOIWEviLUFGm6AGk3Yzh2MhZnK52CvdHN6EN5jhcfgv/dv2lbr4QQogcSeAvIhkpegMJTNpd6Cl6adFJRP+yj8vHznM66B/sKx4n4WYVQprOprxnrULZhxBCiAeTBP5ClpGi9xn+J7/BCbNCS9Ez1s3fdonEo5Gcc9rDzZbzsLOOx5DyKF0fm4SZmfxzCiGEuDuJFIUoe4reXrsQPPpMve8Uvdvr5ienx7Lb+3tcau9Gpdji5TWV2rUeKaQjEEII8aCTwF8Ibk/R29P0Sxo9fH8perfXzTezT2XvzbWktNtCeberJMTUoX2n77GzLZysACGEEGWDBP77dPDfnym/cdgtKXqN7zFFL6e6+bYNnbn63zI2WR6m6mNHKGeeRjnzl2n/2DB5eE8IIUSBSeC/R1HXIji9cBCBMX/fd4pe2vVE4nZE3FE3n/TzbJo8ghtdo6lR4wyJ8S74NZxB5cqBhXw0QgghygoJ/AWkDQZCf5tBzb1j8dcJ95yip9M1ScciidsefkfdfKsqNlz+6mv+OLUb135n8XSMIimuGZ06z8TKqnBSAYUQQpRNEvgLaPusQTSPWMQxi3pEPzG1wCl66bEpxO+6TPyOy6TfTL6jbn7ymTPseeYj9gUmUvXJo2hthovjhzRu/1IRHZEQQoiyRAJ/AXm2fpEdB6oS9ORgzMzzVxFPa01KWAxx2yNIPHQd0vUddfO11kT9uJR1Py4lrecNvCtfJD7Wk5at5uHqWjjNe4QQQggJ/AVUrV4TqtVrkq9lDclpJOy9Sty2CNKuJKCszbFv5nFH3fy0yEjCPvyYDTbX8Xj5NNY28aQnd6Frty8wN7csqkMRQghRBkngLwKpl+OJ2x5Bwp6r6JR0LD3tcOlZK8e6+XGbN7NrzHjCHtFU9z1Baoo1VSt/RZ26XU00eiGEEA8yCfyFRKcZSDx8g7jt4aScjQELha2fe8bDejnUzTckJ3Nl0iT+ObQX6xcj8XK/QnxMLTo89D12dhVNdBRCCCEedEUW+JVSc4FuwFWttW/mtKVAncxFnIForXVADuuGAbFAOpCmtQ7MnO4KLAWqA2HAU1rrqKI6hvxIi04ifsdl4nddxhCXirmrNU6PeGMbWBFzu5wv0ycdP8HJwcPY1tCcqi8cxcIyBUvVl0cf+0Ry84UQQhSpojzjnw9MAxZkTdBa9856rZSaDNy8y/rttNbXb5s2DPhHaz1OKTUs8/3QQhtxPhnr5m+PIOnoDQCs67pi38yDcrVcUGY5B29tMBC1aBE7FvzI9V5p1Kh1msQEZxr6zqVK1abFeQhCCCHKqCIL/FrrTUqp6jnNUxmntU8B7Qu42ceAtpmvvwc2UsyBP2H/VWL+OkfajSTM7CxxaFMVu+BKWLha33W91KtXufjBh2xUsbi8cRVPpygSYwPp9MgcrKzsi2n0QgghyjpT3eNvBVzRWp/MZb4G/lJKaeA7rfXMzOkVtdYRma8vA7neDFdKvQq8CuDl5VU4owZ0sgEzeytcH6qGja8byiLvevyx//zDkVHjOPCwHV6BRwGFk/0QOrR/rdDGJYQQQuSHqQJ/H2DJXea31FpfUkpVAP5WSh3TWm/KvoDWWmd+MchR5peFmQCBgYG5LldQtkEVsQuulK9lDQkJXBk/ga2795DyYio+VQ8SH1uJlq3m4+paq7CGJIQQQuTbvbePu0dKKQvgCTIe0suR1vpS5t9XgZVAcOasK0opj8zteABXi3a0d8rvw3eJhw5zvGcvfo8+itVr53GvcoG0xIfp2m2jBH0hhBAmU+yBH+gIHNNaX8xpplLKTinlkPUa6AQcypz9K9A383Vf4JciHmuB6fR0rs+axfZX32BbJ0X1R/dgYZFG5UqTebjrN1KQRwghhEkVZTrfEjIexHNTSl0Ehmut5wBPc9tlfqWUJzBba92FjPv2KzPPrC2AH7TWf2QuOg5YppR6CThHxgOCJUZqRASXhgxla0o8lm+lUq3CUeJjatDhoYWSmy+EEKJEUFoX2u3vEiswMFCHhoYW6T5i1q7lzKjP2NnOkyqtjmBpmYyV6kPr9iMlN18IIUSxUkrtzqqBczup3Hef0uPiuDLmMw5sC+VqXzt8au0mKdGJOvW/w6taC1MPTwghhLiFBP77kLB3LxeGDGW7lwuO/4umsnMkiTGN6NRlHlZWDqYenhBCCHEHCfz3QKelcX3Gd5yZt4ADj1fGK3g/AA6279Gh/UATj04IIYTInQT+Akq5cIHwwUPYFxNHwiA7fLz2Ex9bgZCQebi51zX18IQQQoi7ksBfQNeWLOGfClDh2QtUsI0lNaEDXbtNw9zcytRDE0IIIfIkgb+A1tsewqt7KGmp5fCsOIH6vj1NPSQhhBAi3yTwF1At/xCu3ginXYcFODh4mno4QgghRIFI4C+gZq0GoPVrkpsvhBCiVDJFyd5ST4K+EEKI0koCvxBCCFGGSOAXQgghyhAJ/EIIIUQZIoFfCCGEKEMk8AshhBBliAR+IYQQogyRwC+EEEKUIRL4hRBCiDJEAr8QQghRhkjgF0IIIcoQCfxCCCFEGSKBXwghhChDlNba1GMockqpa8C5QtykG3C9ELcnhBCibCvsuFJNa+2e04wyEfgLm1IqVGsdaOpxCCGEeDAUZ1yRS/1CCCFEGSKBXwghhChDJPDfm5mmHoAQQogHSrHFFbnHL4QQQpQhcsYvhBBClCES+POglEpXSu1TSh1SSi1XStlmTrdQSl1TSo0z9RiFEEKULkoprZSanO39+0qpEZmvRyil3i+qfUvgz1ui1jpAa+0LpAADMqc/BJwAeimllMlGJ4QQojRKBp5QSrkV944l8BfMZqBm5us+wFfAeaC5yUYkhBCiNEoj44G+d4p7xxL480kpZQE8AhxUSlkDHYHfgCVkfAkQQgghCmI68KxSyqk4dyqBP282Sql9QCgZZ/dzgG7ABq11IvAT0EMpZW66IQohhChttNYxwALgreLcr0Vx7qyUStRaB2SfoJTqA7RUSoVlTioPtAf+Lt6hCSGEKOW+BPYA84prh3LGX0BKKUegFeClta6uta4OvIFc7hdCCFFAWutIYBnwUnHtUwJ/wT3+f+3dPWvWUBiH8esWUYQiWnDo4OhbQajYRVDsB1A6CNZ20sWhFFzcHPwAoovg4uhiJ3ERKSji0FkQpIriC7QOPuCiSMF6OySFUNs+AfvUyrl+kCHJTXIy/cnJyTnA08xcbBx7CJyNiJ3/qE2SpP/XTarV+ZZtpxr13xPO3CdJ0hYSEQ+Au5n5qBfX941fkqQtIiJeAr+AmZ7dwzd+SZLK4Ru/JEkFMfglSSqIwS9JUkGcwEdSKxFxDZgAlqgGH30F9gJ9wD7gfV06mZmz9YyXc5l5ISIuAVfq84PA6/o6j4E54AYw37jdRGa+6u0TSWVycJ+kriLiBHALGMnMxXpFsR2ZuRARI8DVzDzTqD9CNSlJP3AwM783zn0AhjOzU+9frPenNulxpKLZ1S+pjQGgszxxVWZ2MnNhnfpx4B7VL0mjm9A+SS0Z/JLamAH2R8SbiLgTEae71I8B92m/euVYRLxobLv+tsGSVmfwS+oqM78Bx4HLwBdguu6i/0NEDFP1DnwCngDHIqK/yy2mM3Oosf3YwOZLajD4JbWSmUuZ+SwzrwNTwLk1SseBw/W3/HfA7nVqJW0yg19SVxFxKCIONA4NAR9XqdsGnAeONlavHMXVK6Utw9/5JLXRB9yOiD3AT+AtVbf/SqeA+RUD/54DgxExkJmf17j+WEScbOxPZubsBrRb0gr+zidJUkHs6pckqSAGvyRJBTH4JUkqiMEvSVJBDH5Jkgpi8EuSVBCDX5Kkghj8kiQV5DfrF6UwaTlI8gAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f, ax = plt.subplots(figsize=(8, 6))\n",
"all_estimates_matched.T.plot(ax=ax)\n",
"plt.ylabel(y.name)\n",
"plt.xlabel(a.name)\n",
"plt.xticks([0, 1], [\"PA\", \"NJ\"]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case `Matching` and `MarginalOutcomeEstimator` are in fact identical because the solution to the matching problem is the same on the matched data. Once again, we see that all of the methods find a slight positive effect."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}