{
"cells": [
{
"cell_type": "markdown",
"id": "e401eb2c",
"metadata": {},
"source": [
"# R: DoubleML for Difference-in-Differences"
]
},
{
"cell_type": "markdown",
"id": "8ff5d0ae",
"metadata": {},
"source": [
"In this example, we demonstrate, how `DoubleML` can be used in combination with the [did package for R](https://bcallaway11.github.io/did/index.html) in order to estimate group-time average treatment effects in difference-in-difference (DiD) models with multiple periods."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e8e29a2f",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:43:43.070559Z",
"iopub.status.busy": "2024-08-31T09:43:43.069050Z",
"iopub.status.idle": "2024-08-31T09:43:43.573195Z",
"shell.execute_reply": "2024-08-31T09:43:43.571993Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"library(DoubleML)\n",
"library(did)\n",
"library(mlr3)\n",
"library(mlr3learners)\n",
"\n",
"# suppress messages during fitting\n",
"lgr::get_logger(\"mlr3\")$set_threshold(\"warn\")\n",
"\n",
"set.seed(1234)\n"
]
},
{
"cell_type": "markdown",
"id": "19aaa906",
"metadata": {},
"source": [
"# Demo Example from `did`\n",
"\n",
"We will demonstrate the use of `DoubleML` for DiD in the [introductory example](https://bcallaway11.github.io/did/articles/did-basics.html) of the `did` package. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c5d9b7ac",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:43:43.597534Z",
"iopub.status.busy": "2024-08-31T09:43:43.575619Z",
"iopub.status.idle": "2024-08-31T09:43:43.850466Z",
"shell.execute_reply": "2024-08-31T09:43:43.849368Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"text/html": [
"15916"
],
"text/latex": [
"15916"
],
"text/markdown": [
"15916"
],
"text/plain": [
"[1] 15916"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"A tibble: 6 × 7\n",
"\n",
"\t| G | X | id | cluster | period | Y | treat |
|---|
\n",
"\t| <dbl> | <dbl> | <int> | <int> | <dbl> | <dbl> | <dbl> |
|---|
\n",
"\n",
"\n",
"\t| 3 | -0.8762330 | 1 | 5 | 1 | 5.562556 | 1 |
\n",
"\t| 3 | -0.8762330 | 1 | 5 | 2 | 4.349213 | 1 |
\n",
"\t| 3 | -0.8762330 | 1 | 5 | 3 | 7.134037 | 1 |
\n",
"\t| 3 | -0.8762330 | 1 | 5 | 4 | 6.243056 | 1 |
\n",
"\t| 2 | -0.8738481 | 2 | 36 | 1 | -3.659387 | 1 |
\n",
"\t| 2 | -0.8738481 | 2 | 36 | 2 | -1.274099 | 1 |
\n",
"\n",
"
\n"
],
"text/latex": [
"A tibble: 6 × 7\n",
"\\begin{tabular}{lllllll}\n",
" G & X & id & cluster & period & Y & treat\\\\\n",
" & & & & & & \\\\\n",
"\\hline\n",
"\t 3 & -0.8762330 & 1 & 5 & 1 & 5.562556 & 1\\\\\n",
"\t 3 & -0.8762330 & 1 & 5 & 2 & 4.349213 & 1\\\\\n",
"\t 3 & -0.8762330 & 1 & 5 & 3 & 7.134037 & 1\\\\\n",
"\t 3 & -0.8762330 & 1 & 5 & 4 & 6.243056 & 1\\\\\n",
"\t 2 & -0.8738481 & 2 & 36 & 1 & -3.659387 & 1\\\\\n",
"\t 2 & -0.8738481 & 2 & 36 & 2 & -1.274099 & 1\\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"A tibble: 6 × 7\n",
"\n",
"| G <dbl> | X <dbl> | id <int> | cluster <int> | period <dbl> | Y <dbl> | treat <dbl> |\n",
"|---|---|---|---|---|---|---|\n",
"| 3 | -0.8762330 | 1 | 5 | 1 | 5.562556 | 1 |\n",
"| 3 | -0.8762330 | 1 | 5 | 2 | 4.349213 | 1 |\n",
"| 3 | -0.8762330 | 1 | 5 | 3 | 7.134037 | 1 |\n",
"| 3 | -0.8762330 | 1 | 5 | 4 | 6.243056 | 1 |\n",
"| 2 | -0.8738481 | 2 | 36 | 1 | -3.659387 | 1 |\n",
"| 2 | -0.8738481 | 2 | 36 | 2 | -1.274099 | 1 |\n",
"\n"
],
"text/plain": [
" G X id cluster period Y treat\n",
"1 3 -0.8762330 1 5 1 5.562556 1 \n",
"2 3 -0.8762330 1 5 2 4.349213 1 \n",
"3 3 -0.8762330 1 5 3 7.134037 1 \n",
"4 3 -0.8762330 1 5 4 6.243056 1 \n",
"5 2 -0.8738481 2 36 1 -3.659387 1 \n",
"6 2 -0.8738481 2 36 2 -1.274099 1 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Generate data, original code available at https://github.com/bcallaway11/did/blob/master/vignettes/did-basics.Rmd\n",
"time.periods <- 4\n",
"sp <- reset.sim()\n",
"sp$te <- 0\n",
"\n",
"set.seed(1814)\n",
"\n",
"# generate dataset with 4 time periods\n",
"time.periods <- 4\n",
"\n",
"# add dynamic effects\n",
"sp$te.e <- 1:time.periods\n",
"\n",
"# generate data set with these parameters\n",
"# here, we dropped all units who are treated in time period 1 as they do not help us recover ATT(g,t)'s.\n",
"dta <- build_sim_dataset(sp)\n",
"\n",
"# How many observations remained after dropping the ``always-treated'' units\n",
"nrow(dta)\n",
"#This is what the data looks like\n",
"head(dta)"
]
},
{
"cell_type": "markdown",
"id": "f508a08b",
"metadata": {},
"source": [
"### Comparison to `did` package\n",
"\n",
"By default, estimation in `did` is based on (unpenalized) linear and logistic regression. Let's start with this default model first."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "945902b5",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:43:43.854280Z",
"iopub.status.busy": "2024-08-31T09:43:43.853129Z",
"iopub.status.idle": "2024-08-31T09:43:45.084442Z",
"shell.execute_reply": "2024-08-31T09:43:45.083076Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"att_gt(yname = \"Y\", tname = \"period\", idname = \"id\", gname = \"G\", \n",
" xformla = ~X, data = dta)\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"Group-Time Average Treatment Effects:\n",
" Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band] \n",
" 2 2 0.9209 0.0640 0.7432 1.0986 *\n",
" 2 3 1.9875 0.0638 1.8102 2.1648 *\n",
" 2 4 2.9552 0.0636 2.7786 3.1318 *\n",
" 3 2 -0.0433 0.0659 -0.2264 0.1399 \n",
" 3 3 1.1080 0.0632 0.9325 1.2836 *\n",
" 3 4 2.0590 0.0628 1.8845 2.2335 *\n",
" 4 2 0.0023 0.0647 -0.1774 0.1820 \n",
" 4 3 0.0615 0.0645 -0.1176 0.2407 \n",
" 4 4 0.9523 0.0671 0.7660 1.1387 *\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"P-value for pre-test of parallel trends assumption: 0.60857\n",
"Control Group: Never Treated, Anticipation Periods: 0\n",
"Estimation Method: Doubly Robust\n"
]
}
],
"source": [
"# estimate group-time average treatment effects using att_gt method\n",
"example_attgt <- att_gt(yname = \"Y\",\n",
" tname = \"period\",\n",
" idname = \"id\",\n",
" gname = \"G\",\n",
" xformla = ~X,\n",
" data = dta\n",
" )\n",
"\n",
"# summarize the results\n",
"summary(example_attgt)"
]
},
{
"cell_type": "markdown",
"id": "ed9b4f0e",
"metadata": {},
"source": [
"### Using ML for DiD: Integrating `DoubleML` in `did`\n",
"\n",
"As described in our [Section on DiD models in the user guide](https://docs.doubleml.org/stable/guide/models.html#difference-in-differences-models-did), [Sant'Anna and Zhao (2020)](https://linkinghub.elsevier.com/retrieve/pii/S0304407620301901) have developed a doubly robust DiD model which is compatible with ML-based estimation. As this doubly robust model is internally used in `did`, it is possible to use `DoubleML` here to obtain valid point estimates and confidence intervals. For this, we need to write a wrapper around a `DoubleMLIRM` model and pass it to `did` as a custom estimation approach. Once this is implemented, we can use all the nice features and advantages of the `did` package.\n",
"\n",
"For now, let's abstract from using fancy ML algorithms to keep the comparison to the classic `did` implementation simple. Hence, we will use linear and logistic regression for the nuisance compontents in the DiD model."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a0e92c95",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:43:45.088233Z",
"iopub.status.busy": "2024-08-31T09:43:45.087364Z",
"iopub.status.idle": "2024-08-31T09:43:49.281909Z",
"shell.execute_reply": "2024-08-31T09:43:49.280780Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"att_gt(yname = \"Y\", tname = \"period\", idname = \"id\", gname = \"G\", \n",
" xformla = ~X, data = dta, est_method = doubleml_did_linear)\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"Group-Time Average Treatment Effects:\n",
" Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band] \n",
" 2 2 0.9145 0.0654 0.7418 1.0872 *\n",
" 2 3 1.9951 0.0655 1.8221 2.1681 *\n",
" 2 4 2.9561 0.0653 2.7838 3.1285 *\n",
" 3 2 -0.0418 0.0704 -0.2276 0.1441 \n",
" 3 3 1.1041 0.0649 0.9327 1.2754 *\n",
" 3 4 2.0533 0.0669 1.8768 2.2298 *\n",
" 4 2 -0.0028 0.0716 -0.1918 0.1862 \n",
" 4 3 0.0635 0.0675 -0.1147 0.2416 \n",
" 4 4 0.9609 0.0673 0.7833 1.1386 *\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"P-value for pre-test of parallel trends assumption: 0.64579\n",
"Control Group: Never Treated, Anticipation Periods: 0\n"
]
}
],
"source": [
"# DoubleML wrapper for did\n",
"set.seed(1234)\n",
"doubleml_did_linear <- function(y1, y0, D, covariates,\n",
" ml_g = lrn(\"regr.lm\"),\n",
" ml_m = lrn(\"classif.log_reg\"),\n",
" n_folds = 10, n_rep = 1, ...) {\n",
" \n",
" # warning if n_rep > 1 to handle mapping from psi to inf.func\n",
" if (n_rep > 1) {\n",
" warning(\"n_rep > 1 is not supported.\")\n",
" }\n",
" # Compute difference in outcomes\n",
" delta_y <- y1 - y0\n",
" # Prepare data backend\n",
" dml_data = DoubleML::double_ml_data_from_matrix(X = covariates, y = delta_y, d = D)\n",
" # Compute the ATT\n",
" dml_obj = DoubleML::DoubleMLIRM$new(dml_data, ml_g = ml_g, ml_m = ml_m, score = \"ATTE\", n_folds = n_folds)\n",
" dml_obj$fit()\n",
" att = dml_obj$coef[1]\n",
" # Return results\n",
" inf.func <- dml_obj$psi[, 1, 1]\n",
" output <- list(ATT = att, att.inf.func = inf.func)\n",
" return(output)\n",
"}\n",
"\n",
"example_attgt_dml_linear <- att_gt(yname = \"Y\",\n",
" tname = \"period\",\n",
" idname = \"id\",\n",
" gname = \"G\",\n",
" xformla = ~X,\n",
" data = dta,\n",
" est_method = doubleml_did_linear)\n",
"\n",
"\n",
"summary(example_attgt_dml_linear)"
]
},
{
"cell_type": "markdown",
"id": "344bfbf4",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"source": [
"Any differences from the default `did` implementation arise due to sampling randomness, because `DoubleML` uses cross-fitting internally, which is not necessary if classical parametric estimation methods are used.\n",
"\n",
"Next, let's demonstrate how we can use more complex ML learners. For this, we just have to pass another `mlr3` learner through the wrapper, for example a random forest. Please note that the original data generating process is linear, such that we don't expect random forest to lead to better results than the linear learners. We provide a variant of the wrapper that includes an evaluation of the nuisance predictions at the end of this notebook."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "23e26476",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:43:49.285302Z",
"iopub.status.busy": "2024-08-31T09:43:49.284393Z",
"iopub.status.idle": "2024-08-31T09:44:05.252074Z",
"shell.execute_reply": "2024-08-31T09:44:05.250911Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"att_gt(yname = \"Y\", tname = \"period\", idname = \"id\", gname = \"G\", \n",
" xformla = ~X, data = dta, est_method = doubleml_did_rf)\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"Group-Time Average Treatment Effects:\n",
" Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band] \n",
" 2 2 0.9647 0.0736 0.7603 1.1691 *\n",
" 2 3 2.1055 0.0894 1.8571 2.3538 *\n",
" 2 4 3.1295 0.1096 2.8250 3.4339 *\n",
" 3 2 0.0820 0.0719 -0.1177 0.2818 \n",
" 3 3 1.2075 0.0707 1.0112 1.4039 *\n",
" 3 4 2.2865 0.0916 2.0323 2.5408 *\n",
" 4 2 0.1807 0.0743 -0.0257 0.3871 \n",
" 4 3 0.2451 0.0662 0.0611 0.4290 *\n",
" 4 4 1.1401 0.0711 0.9425 1.3376 *\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"P-value for pre-test of parallel trends assumption: 2e-05\n",
"Control Group: Never Treated, Anticipation Periods: 0\n"
]
}
],
"source": [
"# DoubleML wrapper for did with random forest learner\n",
"set.seed(1234)\n",
"\n",
"doubleml_did_rf <- function(y1, y0, D, covariates,\n",
" ml_g = lrn(\"regr.ranger\"),\n",
" ml_m = lrn(\"classif.ranger\"),\n",
" n_folds = 10, n_rep = 1, ...) {\n",
"\n",
"# warning if n_rep > 1 to handle mapping from psi to inf.func\n",
"if (n_rep > 1) {\n",
"warning(\"n_rep > 1 is not supported.\")\n",
"}\n",
"# Compute difference in outcomes\n",
"delta_y <- y1 - y0\n",
"# Prepare data backend\n",
"dml_data = DoubleML::double_ml_data_from_matrix(X = covariates, y = delta_y, d = D)\n",
"# Compute the ATT\n",
"dml_obj = DoubleML::DoubleMLIRM$new(dml_data, ml_g = ml_g, ml_m = ml_m, score = \"ATTE\", n_folds = n_folds)\n",
"dml_obj$fit()\n",
"att = dml_obj$coef[1]\n",
"# Return results\n",
"inf.func <- dml_obj$psi[, 1, 1]\n",
"output <- list(ATT = att, att.inf.func = inf.func)\n",
"return(output)\n",
"}\n",
"\n",
"example_attgt_dml_rf <- att_gt(yname = \"Y\",\n",
" tname = \"period\",\n",
" idname = \"id\",\n",
" gname = \"G\",\n",
" xformla = ~X,\n",
" data = dta,\n",
" est_method = doubleml_did_rf)\n",
"\n",
"summary(example_attgt_dml_rf)"
]
},
{
"cell_type": "markdown",
"id": "25d43c1d",
"metadata": {},
"source": [
"We can see that the results are not dramatically different from the results before. We can observe from the larger standard errors that the default random forest learners seems to be a less precise prediction rule."
]
},
{
"cell_type": "markdown",
"id": "676c1a49",
"metadata": {},
"source": [
"### Exploiting the Functionalities of `did`\n",
"\n",
"The `did` package offers various tools for multi-period DiD models, for example plotting the group-time average treatment effects, which can be exploited just as in the native `did` usage."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d8b9596e",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:44:05.256160Z",
"iopub.status.busy": "2024-08-31T09:44:05.255311Z",
"iopub.status.idle": "2024-08-31T09:44:05.704679Z",
"shell.execute_reply": "2024-08-31T09:44:05.703497Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAIAAAByhViMAAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nOzde0BUdf7/8c/AMMDAcFcUBQFRSCHlspZpZZlrWplmqZWVGZutXaxNTdvd\n1i5mbvbtav2qrU3tm2ttpW1a6xZJeTfAO164yEXkKjBchBlm5vfH2Z3le84ogzIMnJ6Pfzrz\n5n3mvDHFl59zGY3NZhMAAADo/TzcPQAAAAC6BsEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ\n7AAAAFSCYAcAAKASBDsAAACVINgBAACohNbdAzhmtVqLiopKS0vr6upaW1s1Go2Pj09oaGh0\ndHS/fv3cPR0AAEBPpOmBHynW2Ni4fft2o9Ho8KtRUVFXXHGFRqPp5qkAAAB6uB63Ymc2mzMz\nM5uamoQQer0+JibGYDBYLJazZ88WFRW1tbUVFxcHBAQMGzbM3ZMCAAD0LD0u2J08eVJKdSEh\nIePGjdNq/z1hTEzM0KFDMzIy/Pz8rFarvX/37t3FxcUeHh7Tpk3bs2fPmTNnkpKShg4dKoSw\nWq35+flFRUVGo9Fms/n6+vbv3z8hIcHX11faNyMjo7q6WqPR3HHHHfY33Lhxo8lkCgwMnDhx\nohBi165dJSUlQoipU6fm5uaWlJS0tLT4+fklJCTExMR0168KAABAx3pcsCsqKpI2Ro4caU91\nEoPBMGXKFNlJWE9PTyGE1WrNzc0tLS0VQrS1tQkhLBbLjz/+WFVVJbVpNJrGxsaTJ08WFxeP\nGzcuMDDQyXmk9xdC7Nixw2g0+vv7t7S0NDQ07Nu3TwhBtgMAAD1Hzwp2ZrO5oaFBCOHr6xsW\nFqZsuMCldfn5+XFxcSEhIQaDQQhx9OhRKdVFRUWlpKRotdqCgoLs7OzW1ta9e/dOmDDByZHs\nRzSbzZMmTdLpdLW1td9//73Vaj106FB0dDRX+wEAgB6iZwW71tZWacPPz699fefOndJqnN2Y\nMWMGDBjQvjJw4MCUlBT7y4KCAiGETqcbNWqUh4eHECIuLq6srKy8vLy2tra+vt75RTtJQkKC\nTqcTQgQHB/fr16+srKylpcVoNHb2fQAAAFykZz3Hzr761f4qOicNHDjQvt3Y2ChlxJCQECnV\nSUJCQqSNurq6zr5/UFCQfVtaFBRCNDc3d/Z9AAAAXKRnrdh5e3tLG42Nje3rI0eOHD58uBCi\nsrIyJyfH4b72WyKEEGazWdrw8vJq32N/aW9wXvsL/qSlOyGEyWTq7PsAAAC4SM9asdNqtQEB\nAUIIk8lUXV1tr+v1+sDAwMDAQL1ef75926/M2QOcLHjZ85w9mQkhbDabfYHQarWeL/O1fyv7\ntj2JAgAAuF3PCnZCiOjoaGkjKyvLfsmdxGazVVZWOvMm/v7+UuSqqamRbpKVVFRUSBvSOVl7\nvLM/DLm8vPx8T2xuHzTt2xcImgAAAN2sxwW7IUOGSFew1dfXf/vtt4cOHSouLj516tTBgwe3\nbNly8uRJIYSHh0f7E68OxcbGCiHa2tqys7Pb2tqsVuuJEydqamqEEOHh4f7+/kIIaXVQCJGT\nk1NXV1dRUZGTkyN7xopdbm5uXV2dzWY7deqU9D56vd7+DgAAAG7Xs66xE0J4enpee+2127dv\nlz4lNjc3V9YQGBiYmppqvw3ifIYNG1ZdXV1VVXXq1Cnp2XjSUpxer09LS5N6YmNjT5w4YbVa\nq6qqtm7dKoQYPHhwTU2Nw1srQkNDt27d6uHhYT9vO2LEiEv7XgEAALpSjwt2Qgi9Xj9hwoRT\np06VlpbW1taaTCatVuvj4xMaGhoREREREeHMo+OkgJiXl1dcXCx98oSfn9+AAQPi4+PtZ2D9\n/f3HjRuXnZ1tNBq9vb2jo6MTExN/+OEH8Z+nHLeXlJQUEBBQWFjY2toqfaZZ+/twAQAA3E5z\nvkvKINm3b19hYaEQYvLkydIJXAAAgJ6px11jBwAAgItDsAMAAFAJgh0AAIBKcI0dAACASrBi\nBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAA\noBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIE\nOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqoXX3AJ1jtVqLiopKS0vr6upaW1s1Go2Pj09oaGh0\ndHS/fv3cPZ1T6uvrjx49WllZaTKZdDpdaGjo0KFD+/bt6+65AABAr6ex2WzunsFZjY2N27dv\nNxqNDr8aFRV1xRVXaDSabp6qU4xG47/+9S+LxSKrX3XVVQMHDnTLSAAAQDV6zYqd2WzOzMxs\namoSQuj1+piYGIPBYLFYzp49W1RU1NbWVlxcHBAQMGzYMHdPeiFZWVlSqrvsssv69u1bU1Nz\n+PBhIcSBAwcIdgAA4BL1mmB38uRJKdWFhISMGzdOq/335DExMUOHDs3IyPDz87Narfb+3bt3\nFxcXe3h4TJs2bc+ePWfOnElKSho6dKgQwmq15ufnFxUVGY1Gm83m6+vbv3//hIQEX19fad+M\njIzq6mqNRnPHHXfY33Djxo0mkykwMHDixIlCiF27dpWUlAghpk6dmpubW1JS0tLS4ufnl5CQ\nEBMT4/BbMJvN1dXVQoi+ffsmJSUJIcLDwysqKqqqqpqamqQzsy74lQMAAL8UvSbYFRUVSRsj\nR460pzqJwWCYMmWK7CSsp6enEMJqtebm5paWlgoh2trahBAWi+XHH3+sqqqS2jQaTWNj48mT\nJ4uLi8eNGxcYGOjkPNL7CyF27NhhNBr9/f1bWloaGhr27dsnhHCY7bRa7bRp04QQHh7/vWfF\nfirc/oYAAAAXp3cEO7PZ3NDQIITw9fUNCwtTNlzg0rr8/Py4uLiQkBCDwSCEOHr0qJTqoqKi\nUlJStFptQUFBdnZ2a2vr3r17J0yY4ORI9iOazeZJkybpdLra2trvv//earUeOnQoOjpaOZJG\no5FF0rKyMmkNLzw8nGAHAAAuUe8Idq2trdKGn59f+/rOnTul1Ti7MWPGDBgwoH1l4MCBKSkp\n9pcFBQVCCJ1ON2rUKGnlLC4urqysrLy8vLa2tr6+3vlFO0lCQoJ0CjU4OLhfv35lZWUtLS1G\no7HD9ykvL9+9e7cQwtPTs/2EAAAAF6d3PMfOvvrV/io6J7W/KaGxsVHKiCEhIe3Ph4aEhEgb\ndXV1nX3/oKAg+7a0KCiEaG5uvvBehYWFP/30U1tbm6en5+jRo+07AgAAXLTesWLn7e0tbTQ2\nNravjxw5cvjw4UKIysrKnJwch/vab4kQQpjNZmnDy8urfY/9pb3Bee3PrtrvfjCZTBfY5fDh\nw0ePHhVC+Pj4XHXVVQ5PLgMAAHRW7wh2Wq02ICDAaDSaTKbq6mp7EtLr9dKGLPC1135lzh7g\nZMHLnufa35dqs9msVqu0u9VqPV/mM5lM9jHsb2tPokpHjhyRUl1QUNDVV1/dPncCAABcit5x\nKlYIER0dLW1kZWXZL7mT2Gy2yspKZ97E399filw1NTXSTbKSiooKaUM6J2uPd/aHIZeXl5/v\nSc7S3Q+ybXvUkykpKTly5IgQok+fPtdffz2pDgAAdKHesWInhBgyZEhhYWFDQ0N9ff23334b\nGxsbGBhotVqNRmNJSYn0iDsPD48Oo1JsbGxubm5bW1t2dnZKSoqHh0deXl5NTY0QIjw83N/f\nXwgREBBQVlYmhMjJyUlOTm5tbc3JydFqte2zoF1ubm5YWFhgYGBRUZH0Pnq9PiAgQNlpNpvt\n54ujoqLsaVISEhJCzgMAAJei1wQ7T0/Pa6+9dvv27dKnxObm5soaAgMDU1NT7bdBnM+wYcOq\nq6urqqpOnTolPRtPWorT6/VpaWlST2xs7IkTJ6xWa1VV1datW4UQgwcPrqmpcXhrRWho6Nat\nWz08POw3dowYMcLhoSsrK1taWqTtrKws2Vf5VDEAAHCJek2wE0Lo9foJEyacOnWqtLS0trbW\nZDJptVofH5/Q0NCIiIiIiAhnPihWCoh5eXnFxcXSJ0/4+fkNGDAgPj7efgbW399/3Lhx2dnZ\nRqPR29s7Ojo6MTHxhx9+EP95ynF7SUlJAQEBhYWFra2t0meakc8AAIBbaM536RgubN++fYWF\nhUKIyZMnSydwAQAA3KvX3DwBAACACyPYAQAAqATBDgAAQCW4xg4AAEAlWLEDAABQCYIdAACA\nShDsAAAAVIJgBwAAoBIEOwAAAJXoQcGuqqqqoKDAbDa7exAAAIBeqQcFu0WLFg0ePLigoMDd\ngwAAAPRKPSjYAQAA4FIQ7AAAAFSCYAcAAKASBDsAAACVINgBAACohNbdAwAAALjKmeZzRY1N\nrj5KlL9fhN7X1UdxhquCXUFBwUsvvfT999+XlpYGBASMGTPm6aefHjVqlIsOBwAAoLS/pvaj\nE/muPsrdcTFTBg109VGc4ZJgd/z48TFjxjQ0NMyYMWPw4MF5eXmffvrpli1bMjMzR48e7Yoj\nAgAAKMUHBdwVF+N8f0FD4+6KqlF9Q+MCApzfa1hwYOdHcwmXBLtHH3307Nmz27Ztu+aaa6TK\n9OnTb7vtttdee41gBwAAuk2swT/W4O98/7YzFbsrqkaEBN8woL/rpnIdlwS7K664Ii0tzZ7q\nhBBTpkzx8vI6fvy4Kw4HAAAA4aJg9/zzz8sq5eXlZrM5JqYTa6EAAADoFJc/7qS5uXnbtm2T\nJ082GAy///3vlV81/ofFYtFoNK6eBwAAwKET9cbvTpcLIX4oqzheZ3T3OBfDtY87CQoKqq+v\nF0LMnj37yy+/jI2NlTUsWbJk586d9pf+/p04Cw4AANBVMs9UvH30hLSdZ2x4JuvAvMuGXB/R\nz71TdZZrg91vf/vbs2fPHj58+JNPPjl16tSaNWtk2S40NHTAgAHSdlVVlc1mc+k8AAAASk1t\nbR8elz8V5aMT+WlhoQE6L7eMdHE03ZOltm3bdvPNNw8ePDgnJ8fDw/H53zlz5qxZs+bYsWPx\n8fHdMBIAAIBkf03tiv2HlfVFI4alhYV2/zwXrZs+UmzcuHG33nrrwYMHuTEWAAD0NDbheJ2r\n151K7Ppgd/r06REjRtx7772yektLixCiqcnlH+sBAADQKYMNBi/FGUUvD48hgQa3zHPRuj7Y\nDRgw4OzZsxs2bNizZ4+9eOLEia1bt/r7+ycmJnb5EQEAAC5FgM5r9hD5Q9nuiosO0uncMs9F\nc8nNE2+++ebtt99+9dVXT58+ffDgwadPn/7ss8+ampreeustHx8fVxwRAADgUtw4MKKfr+8n\n+YVFDU2Rfvq7h8Qkh4a4e6hOc0mwmzp16o4dO1566aXvv//+s88+MxgMo0ePfuyxx2655RZX\nHA4AAODSjQwNrjOZ3jl64sbIiN6Y6oTrHndyxRVXfPnlly56cwAAACh1012xAAAAcDXXPqAY\nAADAjUoam47XNzjff6LeKITIrTMK0YmPOR0SaBjk79fp4VyAYAcAAFTrcG39RyfkHynRoe3l\nldvLK53vvzsuhmAHAADgWonBgb9JGOLqo/Scx90R7AAAgGpF+vtF9oy1tO7BzRMAAAAqQbAD\nAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQ\nCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYId\nAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCa27BwCA\nbmKx2VosFlcfRavx8Pbk38wA3INgB+CX4kBN7coDR1x9lAkD+qcnxLn6KADgEMEOwC9FsLfu\nyr5hzvcbzeajtfURet8ofz/n94ox+Hd+NADoGq4KdrW1tcuXL//ss8/Ky8sjIiJGjBixZMmS\nK6+80kWHA4AOxRj8n0i6zPn+3Lr6ZVkH0/qE3h0X47qpAKALuSTYnT17NjU19dSpUzfddNN9\n991XUFCwYcOGf/7zn3v37k1KSnLFEQEAAOCSYPfMM8+cOnXqzTfffOSRR6TKbbfdNn369CVL\nlmzevNkVRwQAAIBLbt3y8vIaP378vHnz7JVp06b5+voeOeLyy5YBoEtYbLYjtXVCiKLGpnqT\nyd3jAIBTXLJi9+qrr8oqJpOpra1t4MCBrjgcAHStmpbW5fsPn25qFkIcqKl9fFfWI8PjU8NC\n3D0XAHSgmx629O6775rN5lmzZnXP4QDgUrx99ISU6iTNbW2rjxyvY90OQI/XHY87yczMXLRo\n0dixYx966CHZl5566ql9+/ZJ201NTf7+PCYAgJudbTUdrq2TFZva2rKqzo4f0M8tIwGAk1we\n7NavX3///fcnJiZu2rRJq5Uf7ty5c0aj0f5So9G4eh4AuLAms9lhvbGtrZsnAYDOcmGws9ls\ny5Yte+6552688cZPP/3UYDAoe1566aW2//ysfPjhh7Ozs103DwA4o6+vj5eHh9lqldUH+und\nMg8AOM9V19jZbLb09PTnnnvu0Ucf/frrrx2mOiGEXq8P+A9PT0+bzeaieQDASd6enrdFR8qK\nCUGByaHBbpkHAJznqhW7J5544sMPP3zxxReXLl3qokMAgItMjY7UaDQbT5W0WCwaIcb063vv\nkFgPrhUB0OO5ZMXuiy++eP311xcsWECqA9AbeWg006Ijnxo5XAhxU9TAR4fHB+q83D0UAHTM\nJSt2ixcvFkJYrdYlS5bIvvTUU08FB3M6A0AvIC3QebBOB6D3cEmwy8/PF0K8+eabyi899NBD\nBDsAAABXcEmw4x4IAD2QxWZrsVic75eazVZrU2cedKLVeHh7dtOz3wFApjseUAwAPcGBmtqV\nBzr9idXflJR9U1LmfP+EAf3TE+I6exQA6BIEOwC/FAE6r6SQIFcfJcLP19WHAIDzIdgB+KWI\nCzD8ITnJ3VMAgAtxIQgAAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAA\nACpBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpB\nsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMA\nAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACrhwmBnNpuX\nLl3q6emZlpbmuqMAAABAonXR++bm5s6ePfvkyZMuen8AAADIuCTYGY3G1NTU4cOHZ2dnJyYm\nuuIQl+hM87nf79vfiR1sQmj++x8nDQ8OfPLyYZ2dDQAA4OK4JNi1tbXNnz9/xYoVXl5ernj/\nS+eh0fh7deJ7t9hEdUuLt6dnkK4Te/loPTs/GgAAwEVySbALCQlZtWqVK965q4T7+rxx1a+c\n7z/bavrt9j0jQ4N/l3SZ66YCAAC4FNwVCwAAoBKuunnCSatXrz5y5Ii0XVhYqNfr3TsPAABA\n7+XmYHf8+PG9e/faX3p69tCL0irPtQgh6k3mFovFp6cOCQAAfuHcHOwWLlzY1NQkbS9btmz/\n/s7cqdotLDbbe7knt52pEEIcq6tfsPPnh4YNSQ4NcfdcAAAAcm4OdlFRUfZtX19fi8XixmEc\n+nthsZTqJHUm0+uHj/15VEpfXx83TgUAAKDEzRMXYhPinyVlsuK5Nktmu6gHAADQQxDsLsRk\nsTa1tSnrZ1tN3T8MAADAhbnkVGxmZuY333wjbbe1tZ0+fXrJkiXSy0WLFoWGhrrioK7g7ekR\nqPOqN5lldc7DAgCAHsglwW7Xrl0rV660vywvL7e/TE9P70XBTghxc9TA/80rbF8xeHmN6x/u\nrnkAAADOxyWnYpcsWWI7j7i4OFcc0XVuiRo4ZdBAT82/PyK2v9534eXDgr117p0KAABAiWvs\nOqDRiLvjYl66IkUIkRQS9MqVqQlBAe4eCgAAwAGCnVP8tVohhF6rtS/dAQAA9DQEOwAAAJUg\n2AEAAKiEmz95wl3ONJ/7/b5OfHyZTQghRFZVzdzMXc7vNTw48MnLh3VyNAAAgIv0Cw12HhqN\nv1fnvndDJ/uFED5az87uAgAAcNF+ocEu3Nfnjat+5e4pAAAAuhLX2AEAAKgEwQ4AAEAlCHYA\nAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAq\nQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbAD\nAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAq4apgV1dX\n9/jjj0dHR+t0uoiIiPT09DNnzrjoWAAAABBCaF3xpiaTafz48dnZ2dOnT09JScnPz1+7dm1G\nRkZWVlZwcLArjggAAACXBLvVq1dnZ2evXLly8eLFUmXixIkzZ85cvnz5qlWrXHFEAAAAuORU\n7Nq1aw0Gw4IFC+yVGTNmxMXFrVu3zmazueKIAAAA6Ppg19LScujQoVGjRnl7e7evjx07trKy\nsrCwsMuPCAAAAOGKYFdSUmKxWCIjI2X1QYMGCSEKCgq6/IgAAAAQrrjGrqGhQQjh5+cnq/v7\n+9u/avfVV18VFRVJ2+Xl5d7e3u+//36fPn3a90yePDkpKUn2bj/99NPOnTtlxREjRtx4442y\nYkFBwWeffSYrhoWFPfDAA7Jia2vra6+9pvyOHn/8cdnqoxDigw8+qK6ulhVnzJgRExMjK377\n7bcHDhyQFceMGTN27FhZ8dChQ1u2bJEVo6OjZ86cKSuePXv2/ffflxW9vLx+97vfKed/4403\nzp07Jyvef//9ffv2lRU///zzvLw8WfHXv/51cnKyrLh79+7MzExZcfjw4TfffLOsWFxcvH79\nelkxKCho3rx5sqLFYnF4CeYjjzyi/O20Zs2a8vJyWXHatGlDhw6VFb/77rusrCxZcdSoUddd\nd52smJub+9VXX8mKAwYMmD17tqxoNBrfeecd5ahPPfWUsvj222/LftsLIe65556IiAhZcdOm\nTceOHZMVr7/++l/96ley4s8///z999/LivHx8VOnTpUVz5w5s3btWlnR39//4YcfVo768ssv\nW61WWfGhhx4KDAyUFT/55JOSkhJZ8ZZbbhk2bJismJmZuXv3blkxJSVlwoQJsmJeXt7nn38u\nK4aHh8+ZM0dWPHfu3BtvvKGc/8knn9Rq5T/T3n///bNnz8qKs2bNkv6p2d6WLVsOHTokK159\n9dVXXXWVrHjgwIFvv/1WVoyNjb3jjjtkxerq6g8++EBW9Pb2fvzxx5Xzv/baa62trbLiAw88\nEBYWJit++umnyrMfN95444gRI2TFHTt2bN++XVa8/PLLJ02aJCueOnVqw4YNsmJoaGh6erqs\naDab/+d//kc5/4IFC3x8fGTFv/71r5WVlbLi7bffPnjwYFlx69atOTk5suLo0aOvueYaWfHI\nkSNff/21rBgVFXXnnXfKinV1de+++66s6OnpuXDhQuX8b731VlNTk6x433339evXT1b88ssv\nT5w4ISvecMMNqampsuLevXt/+OEHWfGyyy6bMmWKrHj69OmPP/5YVgwICPjtb3+rHHXlypXK\n4vz58w0Gg6y4bt26srIyWfHWW29NSEiQFTMyMvbt2ycrpqWljR8/XlY8fvz4xo0bZcX+/fvf\ne++9smJjY+Pq1auVoy5atMjDQ76u9P/+3/+rr6+XFe+66y7lOtE//vGPo0ePyorXXnvtlVde\nKStmZ2f/61//khXj4uKmT58uK1ZUVHz00Ueyoq+v72OPPaac/5VXXmlra5MVf/Ob34SEhMiK\nf/vb3+whx+6mm25KTEyUFR2mmpEjR06cOFFWLCgoiIqKUv6s+y9bV5P+En3kkUdk9T//+c9C\niC+//LJ98dFHH01tR/mbUgjx17/+VXmUP/7xj8rOBx98UNmpjEpCiMTERGVnXV2dw1+i+vp6\nZfPw4cOVnd98842y8ze/+Y2y85lnnlF2fvjhh8rOG264Qdmp/OtfCKHX65WdNpstNDRU2ZyT\nk6PsvOWWW5Sdq1evVna++OKLys7Zs2crO5U/1IQQgwcPVnYq/0qTlJeXK5vT0tKUnV988YWy\n0+Efy0WLFik7lQFUCDFmzBhlp/IPqhBCo9EoO20228CBA5XNO3fuVHbOmDFD2blq1Spl56uv\nvqrsnD59urJzz549ys7+/fs7HNXhT4qCggJlp/LvWiHExx9/rOxcsmSJsvPhhx9Wdm7atEnZ\nmZycrOysqqpSdgohmpublc3KrC+E+O6775SdygQphHjuueeUncqsIISYNGmSslOZFIUQgYGB\nyk6bzRYQEKBsPnz4sLJT+S9YIcR7772n7Fy2bJmyc+7cucrOrVu3KjsTEhKUnY2NjcpOIURN\nTY2yWZk1hRD/+Mc/lJ0OE8zTTz+t7FT+W0UIMW7cOGWn8l+qQgidTqfstNls4eHhyuZ9+/Yp\nO6dNm6bsfP3115Wd0l98MrNmzVJ2KvO3ECIqKkrZqfzXl6SkpETZPHr0aGXnhg0blJ1PPvmk\nsvPxxx9Xdv79739Xdo4aNUrZqcyUErPZrGxWLosIITIzM5Wdd999t7JzxYoVys633npL2Tll\nyhRlZ3Z2trIzLCxM2Wmz2Xx9fZXNx48fV3Zef/31ys6PPvpI2fmHP/xB2Tlv3jxl5+bNmx3+\nrLPr+hU76WeTconCaDQKIWTRbdasWePGjbN/qwcPHnzttddkKxnK5QohxMyZM5XLeLGxscrO\n5OTkTz/9VFZUrkAIIfR6vbJTCOHwf+Err7wifUftjRw5Utk5b9485eKEcmFDCHHdddcpB3D4\nsyYiIkLZ6enpqewUQnzwwQcmk0lWjI6OVnY+9dRT99xzj6yYkpKi7Jw2bVpcXCSROCQAACAA\nSURBVJysqFwCEUIMHz5cOaq0fCuj1Wod/vo7/J+1YsWK2tpaWXHUqFHKzvvvv1+5OBofH6/s\nHDNmjHIA5WKJVHQ4qkPvvPOOcsXUYdp44oknbr/9dlnR4d+LN91004ABA2RFhwkyLi5OOapy\nWUWyfv16m+L2JuXKrhDi2WefVaarK664Qtl51113KX8LKX/zCCHS0tKUowYFBSk7DQaDw19/\nnU6nLL7++uvKH0fKnx5CiPnz50+ePFlWdPhPuBtuuEE5gHJdRwgRGRmp7PTy8lJ2CiHWrFlj\nNpuV76DsfPrpp+fOnSsrKpeLhBB33HGH8qeNw79BL7/8cuWoDv+x7e3t7fDX3+Gf65dffln5\nb2aHo6anpyvX0S+77DJl5zXXXKMcQHaqRxIeHq7sVK4VSd5///2WlhZZUbmyKIRYuHChcnXQ\n4c//KVOmKH/YOvx/Gh8frxxVr9crOzUajcNff4f/hn/hhRdqampkReXKlhDinnvuUf4RHjJk\niLLziiuuUA6gXKwSQgQHBzsc1eHfVg5XTB3+Bnj00UdvvfVWWdHhH+qJEycqB1CeLRFCxMTE\nKDuVJ+skH3/8scVikRX79++v7HzmmWceeughWfF8qebyyy+XFc+Xahz+rLPTKH+OXyKTyeTn\n5zd27FjZUs1dd921fv36oqKiqKgohzvOmTNnzZo1x44dc/iXLgAAAC6s62+e0Ol0qampe/fu\nbW5uthetVmtmZmZkZOT5Uh0AAAAukUueY/fAAw80Nze//PLL9sp7771XVlamvAgXAAAAXcUl\nnzwxd+7cdevWLVu2LCcnJyUlJTc3d8OGDUlJSQ5vRAIAAECXcMmKnaen55YtWxYuXLh///4X\nXnjhp59+mj9//rZt2xxeBwoAAIAu0fU3T1w06eaJq666ivwHAADg0Ndff32+O3ZFjwp2VVVV\nU6dOVT6gryfo06dPVFSUxWLZv3+/u2cB0E0GDRoUFhZmNpsPHjzo7lkAdJP4+Hh/f/+mpiaH\nj4ztCZqami6wBOaSa+wuTp8+fXbs2OHuKRxbu3btG2+8IT3W0t2zAOgmTz75ZGZmZlBQEH/w\ngV+O22+//dSpUwkJCbm5ue6e5WK45Bo7AAAAdD+CHQAAgEoQ7AAAAFSiB9080ZPV19fv2bPH\nYDA4/EBlAKpUVlZ2+PDh8PBwh5/YC0CV8vLyCgoKYmNjHX6qdc9HsAMAAFAJTsUCAACoBMEO\nAABAJQh2Hdi/f3+/fv00Go1Go/Hw8PDx8Xn66afdPRQA11q/fn1wcLCHh4f0Z9/Hx+eZZ55x\n91AAuk94eLhGowkKCnL3IJ1GsLuQrKyslJSUiooKnU4XFRWl1+tbW1tXrFjx6quvuns0AK7y\nl7/85a677qqrq/Pz8xs0aJCfn19ra+vzzz+/dOlSd48GoDs899xzlZWV7p7iInHzxIWEhYXV\n1NQkJydnZ2dLlalTp27atEmn07W2trp3NgAuotPpzGbzggULXnvtNaly2223ffnll3q9vqmp\nyb2zAXC1xsbGwMBAIYTVag0MDKyrq3P3RJ1DsLuQ8PDw2traqqoq6f+xEMJisWi1WiEEv26A\nWkVFRVmt1tLSUnvl3Llzer3ew8PDYrG4cTAA3SAhIeH48eP33nvv2rVre2Ow41TshVRUVJhM\nJnuqE0KcPXtWCOHp6em+oQC4VnFxcftUJ4TYtm2bEEKn07lnIADd5X//93+PHz8eEhKycOFC\nd89ykQh2nTNhwgQhxOWXX+7uQQB0h5KSkkceeeSmm24SQnCNHaB6999/vxDip59+cvcgF0/r\n7gF6k8cee+zAgQOenp4//viju2cB4HIajUba8Pf3f++99+688073zgPApa699lqz2TxlypRh\nw4YdOnTI3eNcJIKdsyZOnLh161YPD489e/b4+/u7exwALhcREdHc3NzQ0NDY2HjPPfcIIch2\ngFpt3779xx9/9Pb23rRpk7tnuSSciu2YxWKJjo7eunWrTqc7duxYamqquycC0B1Onz5dW1vb\n1tb28MMPWyyWu+++22QyuXsoAC4xceJEIcSGDRvcPcilIth1wGKxhIaGFhUVhYWF1dXVDRky\nxN0TAehub731lp+fn81m+/DDD909C4Cud++99zY3N0dFRWm12s2bN2/evFm6zK6trW3z5s3H\njh1z94CdwKnYDvTv37++vn7IkCEnTpxw9ywAXO6bb7656aab/Pz8Ghoa2telB5303meWAriA\nrVu3CiGKi4tvvvnm9vWmpqabb745IiLi9OnTbhqt0wh2FzJ16tSqqqo+ffqQ6oBfiEmTJtls\ntsbGxj/84Q8vvPCCVPzwww9bWlqEEPPmzXPrdABcYtmyZZmZme0r9fX133zzjZeX1/Tp06+7\n7jp3DXYReEDxhXh4eNhstrCwMOXzqzZv3jxy5Ei3TAXApaZPn/7FF18IIfR6fUhIiNFoNBqN\nQoiUlJSsrCx3TwegOxw6dOjyyy/nAcVqI6Xe6urqMoW8vDx3TwfAJT7//PPf//73Pj4+zc3N\npaWlRqPRy8tr5syZpDoAPR8rdgAAACrBih0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAA\nqATBDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATB\nDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAA\nQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwAAABUQuvuATrHarUWFRWVlpbW1dW1trZqNBofH5/Q\n0NDo6Oh+/fq5ezqn1NfXHzt2rLq6+ty5c15eXmFhYUOHDu3Tp4+75wIAAL2exmazuXsGZzU2\nNm7fvt1oNDr8alRU1BVXXKHRaLp5qk4pKyvbsWOH8td89OjRkZGRbhkJAACoRq9ZsTObzZmZ\nmU1NTUIIvV4fExNjMBgsFsvZs2eLiora2tqKi4sDAgKGDRvm7knPy2q1/vzzzzabzdPTMykp\nKTAwsKys7OTJk0KIQ4cOEewAAMAl6jXB7uTJk1KqCwkJGTdunFb778ljYmKGDh2akZHh5+dn\ntVrt/bt37y4uLvbw8Jg2bdqePXvOnDmTlJQ0dOhQIYTVas3Pzy8qKjIajTabzdfXt3///gkJ\nCb6+vtK+GRkZ1dXVGo3mjjvusL/hxo0bTSZTYGDgxIkThRC7du0qKSkRQkydOjU3N7ekpKSl\npcXPzy8hISEmJsbht2A0Gr28vLy8vKKjo6VJwsPDKysr6+vrGxsbbTZbD19uBAAAPVyvCXZF\nRUXSxsiRI+2pTmIwGKZMmSJLRZ6enkIIq9Wam5tbWloqhGhraxNCWCyWH3/8saqqSmrTaDSN\njY0nT54sLi4eN25cYGCgk/NI7y+E2LFjh9Fo9Pf3b2lpaWho2LdvnxDCYbYLCgqaNGmSrChN\npdfrSXUAAOAS9Y5gZzabGxoahBC+vr5hYWHKhgukovz8/Li4uJCQEIPBIIQ4evSolOqioqJS\nUlK0Wm1BQUF2dnZra+vevXsnTJjg5Ej2I5rN5kmTJul0utra2u+//95qtR46dCg6OvrCQa2i\nosJsNhcVFTU1NWk0msTERCePCwAAcD69I9i1trZKG35+fu3rO3fulFbj7MaMGTNgwID2lYED\nB6akpNhfFhQUCCF0Ot2oUaM8PDyEEHFxcWVlZeXl5bW1tfX19c4v2kkSEhJ0Op0QIjg4uF+/\nfmVlZS0tLUaj8cLvk5mZKW2Eh4cPGzaMu2IBAMCl6x3PsbOvfrW/is5JAwcOtG83NjZKGTEk\nJERKdZKQkBBpo66urrPvHxQUZN+WFgWFEM3NzU7uXlFRcfDgwYs4LgAAgEzvWLHz9vaWNhob\nG9vXR44cOXz4cCFEZWVlTk6Ow33tt0QIIcxms7Th5eXVvsf+0t7gvPYX/ElLd0IIk8l04b1m\nzJjR1tZWWVm5b9++mpqaH3/8cdKkSbKpAAAAOqV3rNhptdqAgAAhhMlkqq6uttf1en1gYGBg\nYKBerz/fvu1X5uzJSRa87HnOnsyEEDabzb5AaLVaz5f52r+VfdueRGUsFktzc7N0w4RWq42I\niIiPjxdCtLS0VFZWnu9bAAAAcEbvCHZCiOjoaGkjKyvLfsmdxGazOZmK/P39pchVU1MjpStJ\nRUWFtCGdk7XHO/vDkMvLy8/3JOf2QdO+7TBo5ufnf/75519//bX07DqJ/XtpPw8AAMBF6B2n\nYoUQQ4YMKSwsbGhoqK+v//bbb2NjYwMDA61Wq9FoLCkpkR5x5+Hh0f7Eq0OxsbG5ubltbW3Z\n2dkpKSkeHh55eXk1NTVCiPDwcH9/fyFEQEBAWVmZECInJyc5Obm1tTUnJ0er1TrMXrm5uWFh\nYYGBgUVFRdL76PV6aX1Rpl+/fhqNxmazHT161GazBQYG1tfX5+XlCSE0Go3Du30BAACc12uC\nnaen57XXXrt9+3bpU2Jzc3NlDYGBgampqfbbIM5n2LBh1dXVVVVVp06dkp6NJy3F6fX6tLQ0\nqSc2NvbEiRNWq7Wqqmrr1q1CiMGDB9fU1Di8xSE0NHTr1q0eHh7287YjRoxweGg/P7/ExMRD\nhw5ZLJbDhw+3/1J8fLzshl8AAIDO6jXBTgih1+snTJhw6tSp0tLS2tpak8mk1Wp9fHxCQ0Mj\nIiIiIiKcecavFBDz8vKKi4ulT57w8/MbMGBAfHy8/Qysv7//uHHjsrOzjUajt7d3dHR0YmLi\nDz/8IBydME1KSgoICCgsLGxtbZU+06z9fbgyl112WXBw8MmTJ2tqasxms5eXV3BwcGxsLJ8n\nBgAALp3mfJeO4cL27dtXWFgohJg8ebJ0AhcAAMC9es3NEwAAALgwgh0AAIBKEOwAAABUgmvs\nAAAAVIIVOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACASvSgYDdnzhyNRnP8\n+HF3DwIAANAr9aBgBwAAgEtBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgElp3DwAAAOAq9SZz\nTWurq48S6q0L1OlcfRRnEOwAAIBq7ayo+uhEvquPcndczJRBA119FGcQ7AAAgGpF+utvGNDP\n+f6y5nNHa+svCwoY4Kd3fq9og1/nR3MJgh0AAFCtxOCgxOAg5/u3nak4Wls/tl/fGwb0d91U\nrsPNEwAAACpBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEk4Fu9ra2oULFw4aNMjb2zsmJmbq\n1Km7d+++8C51dXWPP/54dHS0TqeLiIhIT08/c+ZMVwwMAADQ9cxW62cFRX89kS+EWHey8NOC\nIrPV6u6hOq3jx52cPXs2NTX11KlTN91003333VdQULBhw4Z//vOfe/fuTUpKcriLyWQaP358\ndnb29OnTU1JS8vPz165dm5GRkZWVFRwc3NXfAgAAwKVae7Jga+m/F6FaLJbPC4vrTabfJAxx\n71Sd1XGwe+aZZ06dOvXmm28+8sgjUuW2226bPn36kiVLNm/e7HCX1atXZ2dnr1y5cvHixVJl\n4sSJM2fOXL58+apVq7pqdAAAgC5Rca7FnursvjtdPjlyQKeeVOx2HZ+K9fLyGj9+/Lx58+yV\nadOm+fr6Hjly5Hy7rF271mAwLFiwwF6ZMWNGXFzcunXrbDbbJU4MAADQtUqams5Tb+7mSS5R\nx8Hu1Vdf/e6777y8vOwVk8nU1tY2cKDjz0RraWk5dOjQqFGjvL2929fHjh1bWVlZWFh4iRMD\nAAB0Lb2n43OYeq1nN09yiS7mrth3333XbDbPmjXL4VdLSkosFktkZKSsPmjQICFEQUHBRRwR\nAADAdYYEGsJ8vGXFUG/v+MBAt8xz0Tr9WbGZmZmLFi0aO3bsQw895LChoaFBCOHnJ/80XH9/\nf/tX7Z599tns7Gxpu6qqSuoBAADoTl4eHo8NT/jzwSON5jap4qfVPpoY7+3Zyx4M17lgt379\n+vvvvz8xMXHTpk1a7YX21Wg0sop0dZ2sXlNTc/r06QvsBQAA0A3igwJeG5320cmC7WcqrwoP\nmxsfZ2h3HVpv4WwOtdlsf/rTn+66667rrrtu27ZtISEh5+sMCAgQipU5IYTRaBRCGAyG9sWX\nXnop4z8uu+yyxsbGzo0PAADQRQxeXknBQUKI4cFBvTHVCSdX7Gw2W3p6+ocffvjoo4+++uqr\nnp4XupAwKipKq9UWFRXJ6vn5+UKIIUP+z/Ng9Pr/3kLs6enJPbMAAAAXzakVuyeeeOLDDz98\n8cUX33jjjQunOiGETqdLTU3du3dvc/N/7xC2Wq2ZmZmRkZFRUVGXNC8AAADOo+Ng98UXX7z+\n+usLFixYunSpw4aWlpb9+/dLC3KSBx54oLm5+eWXX7ZX3nvvvbKysvT09EufGAAAAA51fCpW\n+vQIq9W6ZMkS2Zeeeuqp4ODgvLy85OTk8ePHf/fdd1J97ty569atW7ZsWU5OTkpKSm5u7oYN\nG5KSkhYuXNjl3wAAAAAkHQc7aSnuzTffVH7poYcecvjZr56enlu2bHn22Wc/++yzLVu29O3b\nd/78+c8991z7K+oAAABczWy1mqxW5/tbrVYhhMlqbWprc34vnYeHl0ePeDCKpufcrzBnzpw1\na9YcO3YsPj7e3bMAAAA1+Kak7KMT+R33XZq742KmDHL8iVzdrNMPKAYAAOgt+ut9r+wb5uqj\nRPj5uvoQTiLYAQAA1RoZGjwy1MFlY2rVI84HAwAA4NIR7AAAAFSCYAcAAKASBDsAAACVINgB\nAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACo\nBMEOAABAJZwNdmazeenSpZ6enmlpaR02f/TRRxpHXnjhhUubFgAAAOeldaYpNzd39uzZJ0+e\ndPJN6+rqhBB33nlnVFRU+/qYMWM6Ox8AAACc1HGwMxqNqampw4cPz87OTkxMdOZNpWD3u9/9\nzpnlPQAAAHSJjk/FtrW1zZ8/f+fOnXFxcU6+qRTsgoKCLmk0AAAAdEbHwS4kJGTVqlVeXl7O\nv6k92FksltLS0urq6osfEAAAAM5xyV2x9fX1QojXXnutT58+kZGRffr0iY+P/+STT1xxLAAA\nAEicunmis6QVu/Xr1y9evHjAgAG5ubmrV6++++67Gxoa5s2b177zq6++KioqkrbLy8u9vb1d\nMQ8AAMAvgcZmsznf7ePjk5iY+PPPP1+4LSMjo7a29sYbb/Tz85MqR48eTUlJ8fPzO3PmjE6n\ns3c+9thjO3futL88ceLEvn374uPjO/MtAAAAQAgXrdhdf/31ssqwYcMmT5785ZdfHjhw4Fe/\n+pW9PmvWrHHjxknbH3300cGDB10xDwAAwC+BS4KdQ3379hVCNDY2ti9eddVV9u2vvvrKbDZ3\n2zwAAAAq0/U3TzQ2Nr7zzjvr16+X1Y8cOSKEGDRoUJcfEQAAAKJLVuxaWlqOHTtmMBgGDx4s\nhNDr9cuXL6+vr09OTk5ISJB6Nm3atH379uTk5NjY2Es/IgAAAJQ6DnaZmZnffPONtN3W1nb6\n9OklS5ZILxctWhQaGpqXl5ecnDx+/PjvvvtOCOHh4fH2229PnTo1LS1t1qxZERERhw8f3rhx\nY0BAwF/+8hfXfScAAAC/cB0Hu127dq1cudL+sry83P4yPT09NDRUucuUKVN27Njx/PPPf/75\n542NjX379r3nnnv++Mc/Ov/ZFQAAAOiszj3uxKXmzJmzZs2aY8eO8bgTAACAi+CST54AAABA\n9yPYAQAAqATBDgAAQCUIdgAAACpBsAMAAFCJ7vtIMQBwr+qW1v01ta4+ykA/34SgQFcfBQAc\nItgB+KUobmx6/9hJVx9lwoD+BDsA7kKwA/BLEWPwfyLpMuf7Tzc1f1pQlNYn9Op+fZ3fK9zX\np/OjAUDXINgB+KUI9tZd2TfM+f7cunohRITet1N7AYAbcfMEAACAShDsAAAAVIJgBwAAoBIE\nOwAAAJUg2AGAA01tbTsrqoQQx+qMpU3N7h4HAJxCsAMAudKm5id2/by19IwQ4kS9ccnenB/K\nKtw9FAB0jGAHAHJvHTlebzLbX5qt1r+eyKs81+LGkQDAGQQ7APg/qlpaCxsaZcVWi7UbPo4M\nAC6Rs8HObDYvXbrU09MzLS3Nmf66urrHH388Ojpap9NFRESkp6efOXPmEuYEgG7SYrF0qg4A\nPYdTnzyRm5s7e/bskyed/YxFk8k0fvz47Ozs6dOnp6Sk5Ofnr127NiMjIysrKzg4+BKmBQCX\n6+fr4+vpeU4R42IM/m6ZBwCc1/GKndFoTE1N9fDwyM7O9vLycuZNV69enZ2dvXLlyr///e9P\nP/30Bx988PHHHxcWFi5fvvySBwYA1/Ly8Lh7SIysmNYnNDEkyC3zAIDzOg52bW1t8+fP37lz\nZ1xcnJNvunbtWoPBsGDBAntlxowZcXFx69ats9lsFzkpAHSXCQP6Pzw8vo+vtxDCx9NzanTk\nY8PjNe6eCgA61HGwCwkJWbVqlZNrdUKIlpaWQ4cOjRo1ytvbu3197NixlZWVhYWFFzMmAHSv\na/r1fXhYvBDi1wP73zk42tvT090TAUDHnLrGrlNKSkosFktkZKSsPmjQICFEQUFBbGysvXj8\n+PH6+nppu6mpSavt+nkAAAB+Ibo+SDU0NAgh/Pz8ZHV/f3/7V+1Wr169c+dO+0tfX98unwcA\nAOAXwlUrZBqN/HIU6eo6WX3kyJF6vV7a3rdvX1tbm4vmAQAAUL2uD3YBAQFCsTInhDAajUII\ng8HQvjh37lz79pw5czIyMrp8HgAAgF+Irv/kiaioKK1WW1RUJKvn5+cLIYYMGdLlRwQAAIBw\nxYqdTqdLTU3du3dvc3Oz/TSr1WrNzMyMjIyMiorq8iMCgDMOnq197dAx5/stNpsQ4puSsu9P\nlzu/13UR4fcMie24DwBcoAuCXUtLy7FjxwwGw+DBg6XKAw888OCDD7788st/+tOfpMp7771X\nVlb27LPPXvrhAODi6Dw8w319XH0Ug9MPhwKALtdxsMvMzPzmm2+k7ba2ttOnTy9ZskR6uWjR\notDQ0Ly8vOTk5PHjx3/33XdSfe7cuevWrVu2bFlOTk5KSkpubu6GDRuSkpIWLlzoom8DADqU\nEBSwYlSyu6cAABfqONjt2rVr5cqV9pfl5eX2l+np6aGhocpdPD09t2zZ8uyzz3722Wdbtmzp\n27fv/Pnzn3vuOfuZWQAAAHQ5Tc/5jK85c+asWbPm2LFj8fHx7p4FAACg9+n6u2IBAADgFgQ7\nAAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAA\nlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwAAABUwqlgV1dX9/jjj0dHR+t0\nuoiIiPT09DNnzlyg/6OPPtI48sILL3TR2AAAAJDTdthhMpnGjx+fnZ09ffr0lJSU/Pz8tWvX\nZmRkZGVlBQcHO9ylrq5OCHHnnXdGRUW1r48ZM6ZLhgYAAIBSx8Fu9erV2dnZK1euXLx4sVSZ\nOHHizJkzly9fvmrVKoe7SMHud7/7XVpaWhfOCgAAgAvo+FTs2rVrDQbDggUL7JUZM2bExcWt\nW7fOZrM53EUKdkFBQV01JQAAADrUQbBraWk5dOjQqFGjvL2929fHjh1bWVlZWFjocC97sLNY\nLKWlpdXV1V01LgAAAM6ng2BXUlJisVgiIyNl9UGDBgkhCgoKHO5VX18vhHjttdf69OkTGRnZ\np0+f+Pj4Tz75pCsGBgAAgGMdXGPX0NAghPDz85PV/f397V9Vklbs1q9fv3jx4gEDBuTm5q5e\nvfruu+9uaGiYN29e+86nnnpq37590nZTU5P0tgAAALgIHd88IYTQaDSyinR1nbIu+eMf//jI\nI4/ceOON9kQ4e/bslJSUp59++v7779fpdPbOc+fOGY3GCxwIAAAATuog2AUEBAhHK3NSGjMY\nDA73uv7662WVYcOGTZ48+csvvzxw4MCvfvUre/2NN96wb8+ZMycrK8vpyQEAAPB/dHCNXVRU\nlFarLSoqktXz8/OFEEOGDHH+SH379hVCNDY2dnJCAAAAOKWDYKfT6VJTU/fu3dvc3GwvWq3W\nzMzMyMhI2fOHJY2Nje+888769etl9SNHjoj/3HUBAACALtfxc+weeOCB5ubml19+2V557733\nysrK0tPTpZctLS379++X1vCEEHq9fvny5Q8++OCxY8fsu2zatGn79u3JycmxsbFdOj8AAAD+\nreObJ+bOnbtu3bply5bl5OSkpKTk5uZu2LAhKSlp4cKFUkNeXl5ycvL48eO/++47IYSHh8fb\nb789derUtLS0WbNmRUREHD58eOPGjQEBAX/5y19c+90AAAD8gnW8Yufp6blly5aFCxfu37//\nhRde+Omnn+bPn79t2za9Xn++XaZMmbJjx45rrrnm888/X7FixZ49e+65556srKyUlJQuHR4A\nAAD/pTnfx4J1vzlz5qxZs+bYsWPx8fHungUAAKD36XjFDgAAAL0CwQ4AAEAlCHYAAAAqQbAD\nAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQ\nCYIdAACASmjdPYB7WBsbGjO+dfVRtH3C9aOvcfVRAAAAJL/UYNdy7tzPu119FN3goQQ7AADQ\nbX6hwc4zKCTkoSec77c1Ndaue98rdojh1zc7v5eHt0/nRwMAALhITgW7urq6ZcuWbdy4says\nLCwsbPLkyc8//3z//v27dpfupNFqvSIGOt9vNdYLITx8fDu1FwAAQHfqONiZTKbx48dnZ2dP\nnz49JSUlPz9/7dq1GRkZWVlZwcHBXbULAAAALlHHwW716tXZ2dkrV65cvwQSwwAAIABJREFU\nvHixVJk4ceLMmTOXL1++atWqrtoFAAAAl0hjs9ku3JGcnJyfn19VVeXt7W0vDhkyxGg0lpeX\nazSaLtlFCDFnzpw1a9YcO3YsPj7+or4XF7Ia66tWPec97PKgWfe5exYAAADHOniOXUtLy6FD\nh0aNGtU+ogkhxo4dW1lZWVhY2CW7AAAA4NJ1EOxKSkosFktkZKSsPmjQICFEQUFBl+wCAACA\nS9fBNXYNDQ1CCD8/P1nd39/f/tVL2eXDDz88ceKEtF1SUuLr6/vKK6+EhIS077njjjtSU1Nl\n77Z169aMjAxZcdSoUbfddpusePz48b/+9a+yYnh4+BNPyB930tLSsmzZMuV39MdHH27b9i8h\nhKkwr+mnDP3oazRa7auvvlpRUSHrnDt37tChQ2XFL774Yu/evbLi+PHjJ0yYICv+/PPPf//7\n32XFuLi49PR0WbGqquqVV16RFb28vJ5//nnl/M8991xzc7Os+Nhjj0VERMiKa9asyc3NlRWn\nTp165ZVXyorbtm379lv5E56Tk5NnzpwpKxYUFLz33nuyYmho6KJFi2RFi8Xy+9//Xjn/H/7w\nB+k3T3tvvfVWaWmprHjPPfcMHz5cVvzqq6927twpK15zzTWTJ0+WFQ8cOLB+/XpZcdCgQb/9\n7W9lxbq6updeeklW1Gg0K1asUM6/YsWK+vp6WXH+/PlRUVGy4ieffHLw4EFZ8aabbrr66qtl\nxe3bt3/99deyYmJi4uzZs2XFkpKS1atXy4oGg8HhL/XTTz9ttVplxcWLF8v+SAoh3n33XeXq\n+6xZs0aOHCkrfvvtt9u2bZMVR48efeutt8qKubm5a9askRUjIiIee+wxWbGpqcnhb/Xnn3/e\ny8tLVly1alV1dbWsmJ6eHhcXJyt+9tlnWVlZsuKvf/3r66+/Xlbcu3fvF198ISvGx8fff//9\nsmJFRcWrr74qK/r4+Dj8UfOnP/2ptbVVVnziiSfCw8NlxfY/Oe1uu+22UaNGyYrff//9v/71\nL1kxLS3t9ttvlxVPnjz5wQcfyIp9+vR58sknZUWTyfTMM88o53/mmWf0er2s+Prrr585c0ZW\nnDNnTkJCgqy4cePG3bvljxe97rrrJk6cKCvm5ORs2LBBVoyNjX3wwQdlxZqampdffllW9PT0\nXL58uXL+F154obGxUVZ85JFHBg6UPwxh3bp1R44ckRWnTJly1VVXyYo//vjjli1bZMURI0bc\neeedsmJRUdE777wjKwYFBS1ZskRWtNlsS5cuVc6/dOnSwMBAWfHtt98uLi6WFe+6667LL79c\nVty8efNPP/0kK44ZM+aWW26RFQ8fPvzxxx/LipGRkQ8//LCs2NDQ4PCX+sUXX/TwkK8r/fnP\nfz579qysOG/evJiYGFnxb3/72/79+2XFG2+8cdy4cbLirl27Nm3aJCtedtll990nv6qqrOz/\nt3fn8U2Vaf/H76Rp2qb7xlJoKdBCgRZoi6iIUqyACCJYFRdUBlEUVHQG3EaQQVxxF8ZHdFAR\nH4bfKMqjgDIV2WS1C7K0QBe6AKWlW7qnTfL7IzOdzrlDF2ibcvi8X/N6XunFdXKu1KfJN/fJ\nOTnzwQcfKIru7u6LFi2S51+0aFF9fb2iuGDBgoCAAEXxk08+yczMVBTvuuuumJgYRdFuqrn6\n6qunTZumKB4/frxfv37yc91/WJtle5p7/PHHFfU333xTCPHtt99e4iZPPPFEbBOenp7yhJ99\n9pm8F7u/60ceeUTulP+ohBCRkZFyZ1lZmd1fUcYLTxUs+mPj/0q++NhqscgBQgixZcsW+W4f\nfvhhuXPx4sVy5+rVq+XOm266Se5MT0+XOw0Gg9xptVr9/f3l5pSUFLlT/gMWQqxcuVLufPXV\nV+XOGTNmyJ2//PKL3Nm/f3+5U35JsykoKJCbR4wYIXdu2LBB7pRjgRBi4cKFcqec6oQQ1113\nndyZk5Mjd2o0GrnTarXKrwpCiD179sidd911l9z51ltvyZ1yVhBCJCQkyJ379++XO3v27Gl3\nVJ3Ozju9rKwsufOGG+xceXvt2rVyp/yyJISYN2+e3Ck//wohoqOj5c6ioiK5UwhRXV0tN8vv\ntYQQiYmJcufMmTPlzqVLl8qdH3/8sdw5ceJEufPw4cNyp7e3t9xptVq9vLzk5iNHjsidN998\ns9y5atUqudNugpw1a5bcuXXrVrkzIiJC7pTTj01xcbHcPGzYMLnz+++/lzvld1BCiBdeeEHu\nXLNmjdwZFxcnd2ZkZMider1e7rRarXKAFkIcPHhQ7pRfa4UQ77//vtxpe+FTuPvuu+XO3bt3\ny50hISFyp/zuyyYvL09uvvbaa+XO9evXy51yghdCPPXUU3KnvAAhhBg5cqTceebMGbuj1tfX\ny81ygBNC7NixQ+6877775M7XXntN7lyxYoXcOWXKFLkzOTlZ7gwICJA7rVarm5ub3Hz8+HG5\nU35bKIT4/PPP5c4XX3xR7pwzZ47cuWnTJrvPdY1aOHkiIyMjPDz8wQcf/Pzzz5vWFy1atGzZ\nssTExPj4+EvZJCUlpbi42HZ75cqVP/zww/r1620HbRuFhobK0eTMmTPyu8CAgADFtkKI8vJy\n+W/bzc1t8ODBiqLZbJbfBBg3rBuotTj99wkf3nfcl6VzrampUTSHh4fLT805OTnygkHPnj3l\nBbPi4uJTp04pil5eXuHh4YpibW2t/H5Rq9VGR0cLyaFDhxoaGhTFQYMGye+tMzIy5LWlkJCQ\nwMBARbGgoOD06dOKop+fn/yXWVFRIS8tuLi4REZGKooX+tMaOnSo/NYkLS1NXobs37+/j4+P\nopiXl1dYWKgodu/eXc5bJSUl8iqUh4eHfDaPyWSy+4ItLy0LIQ4fPmwymRTFiIgIeVU7Kyur\ntLRUUezdu7f8elNYWJiXl6co+vr69uvXT1GsqqqS3wY4OzvL79eFEMnJyfITQmRkpOLzskKI\n48ePy6/uffv2ldf2Tp8+XVBQoCgGBgbKC5ZlZWXyW1uDwTBo0CBFsaGh4dChQ/L80dHR8jLA\n0aNHa2trFcUBAwbIbyNPnTrV+HTUKCgoSL4A5/nz5+Vw7+3tLa8C1tTUHDt2TFF0cnKSlzaF\nEKmpqWazWVEcPHiw/Cpy8uRJo9GoKPbp00deMDh79qz84mr3qdJoNJ48eVJRdHV1ld/EWiyW\nlJQUef5hw4bJ7w2OHTsmP1WGhYXJa0u5ublyZO/Ro0evXr0URbtPlZ6ennKIr6urO3LkiKKo\n0Wjk9RIhxO+//y4vw9h9qszMzJQXAoKDg7t166Yonjt3Tj62YPepsrKy8vjx44qiXq+PioqS\nR5WXloUQUVFRer1eUUxPT6+qqlIU+/XrJ196LD8/Xz4M1a1bN/mDVaWlpfKnqtzd3eVV2Pr6\nevkohBAiJiZGPo3yyJEj8tv7gQMHykdssrOz5bW9Xr169ejRQ1EsKiqSFyx9fHz69++vKFZX\nV8sHrHQ6nd13JikpKXK8HjJkiKur8lsJTpw4IR/bbH2qsftUWV5e7unpKT/XNWoh2JlMJnd3\n99GjRyvWXe69995169bl5OTIu7yITWy64FmxVrO5cOmzQvoVuY0c5TU5wSEjAQAAXEgLJ0/o\n9frY2NgDBw40XR2xWCw7duwIDg62G9EuYpMuS6PRCHsXZ9FcOCkDAAA4SssB5aGHHqqurm76\n+dNVq1adOXOm8RP9tbW1qampTY+htLjJZUOr1YcqF2yFEPr+XWVNEQAAoFHLFyg2m81jx47d\ntWvXbbfdFhMTk5aWtn79+sjIyH379tk+eXDkyJGoqKj4+PjExMRWbmJXFzwUK4Qwl5wv+fg9\nS5PPiLgOi/VOuNeBIwEAANjV8oqdk5PT5s2bFyxYkJqaumzZsl27ds2dO3f79u3NRLSL2KTL\ncvIL8H/iWcNVo4QQWi9v7zvu875deZo6AABAV9Dyil2n6ZordjZ8pRgAAOj6OAkAAABAJQh2\nAAAAKkGwAwAAUAmCHQAAgEoQ7AAAAFTCznd+XwmsdXV1Gcov0Gyuv6ZaCGE2ltcetfMNlRei\ndffUhyq/uxMAAKCDXKHBzlxRXr5+TVu3asjPadNW+v4D9KFz2roXAACAi3OFBjutwd1j3KSO\n3ouTr19H7wIAAKDRlRvs3K+/0dFTAAAAtCdOngAAAFAJgh0AAIBKEOwAAABUgmAHAACgEl3u\n5InFixf7+Pg4egoAAICu6MMPP9Tr9Rf6V43Vau3MaZpx9OjR+++/PyUlxdGD2OHn5xcaGtrQ\n0PD77787ehYAnaR3797dunWrq6s7evSoo2cB0EnCwsK8vLwqKipOnjzp6Fnsq6qqMhgMF/rX\nLrRiN2TIkL///e8VFRWOHsSOn376acOGDS4uLr/99pujZwHQSVauXPn777/7+Pjwhw9cORYv\nXnzu3LmwsLB169Y5ehb7XF1dm/nXLhTshBADBgxw9Aj22d6vazSa2NhYR88CoJP4+voKIZyc\nnPjDB64cbm5uQggXF5fL9A+fkycAAABUgmAHAACgEl3rUGyXNXbs2Ly8vO7duzt6EACdZ8aM\nGX5+fuHh4Y4eBEDnmTt37t69e0eOHOnoQS5SFzorFgAAAJeCQ7EAAAAqQbADAABQCYJdC1JT\nU3v06KHRaDQajVardXV1feGFFxw9FICOtW7dOl9fX61Wa/vbd3V1Xbx4saOHAtB5unfvrtFo\nLsevwiLYNScpKSkmJubcuXN6vT4kJMRgMNTV1b322mvvvvuuo0cD0FE+/fTTe++9t6yszN3d\nvU+fPu7u7nV1dS+//PLzzz/v6NEAdIalS5cWFhY6eoqLxMkTzQkICCguLo6Ojk5OTrZVpk6d\nunHjRr1eX1dX59jZAHQQvV5fX18/f/789957z1a5/fbbv/32W4PBUFVV5djZAHS0yspKb29v\nIYTFYvH29i4rK3P0RG1DsGtO9+7dS0tLi4qKbP+NhRBms1mn0wkh+L0BahUSEmKxWPLz8xsr\nNTU1BoNBq9WazWYHDgagE0RERBw/fvyBBx5Ys2bN5RjsOBTbnHPnzplMpsZUJ4QoKSkRQjg5\nOTluKAAdKzc3t2mqE0Js375dCKHX6x0zEIDO8tVXXx0/ftzPz2/BggWOnuUiEezaZty4cUKI\noUOHOnoQAJ0hLy/v8ccfnzRpkhCCz9gBqveHP/xBCLFr1y5HD3Lx+OaJNnjyyScPHTrk5OS0\nc+dOR88CoMNpNBrbDQ8Pj1WrVt1zzz2OnQdAhxozZkx9ff2UKVMGDx58+PBhR49zkQh2rTVh\nwoStW7dqtdr9+/d7eHg4ehwAHS4oKKi6urqioqKysvL+++8XQpDtALXavXv3zp07XVxcNm7c\n6OhZLgmHYltmNptDQ0O3bt2q1+vT09NjY2MdPRGAznD69OnS0tKGhoZ58+aZzeb77rvPZDI5\neigAHWLChAlCiPXr1zt6kEtFsGuB2Wz29/fPyckJCAgoKyvj68CBK9CKFSvc3d2tVuvq1asd\nPQuA9vfAAw9UV1eHhITodLpNmzZt2rTJ9jG7hoaGTZs2paenO3rANuBQbAt69uxZXl4eHh5+\n4sQJR88CoMNt2bJl0qRJ7u7uFRUVTeu2C51cvtcsBdCMrVu3CiFyc3MnT57ctF5VVTV58uSg\noKDTp087aLQ2I9g1Z+rUqUVFRYGBgaQ64AoxceJEq9VaWVn54osvLlu2zFZcvXp1bW2tEGLO\nnDkOnQ5Ah1iyZMmOHTuaVsrLy7ds2eLs7JyQkDB27FhHDXYRuEBxc7RardVqDQgIkK9ftWnT\npuHDhztkKgAdKiEhYcOGDUIIg8Hg5+dnNBqNRqMQIiYmJikpydHTAegMhw8fHjp0KBcoVhtb\n6j1//vwZSUZGhqOnA9Ahvvnmmz//+c+urq7V1dX5+flGo9HZ2Xn69OmkOgBdHyt2AAAAKsGK\nHQAAgEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAA\ngEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ\n7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ7AAA\nAFRC5+gB2sZiseTk5OTn55eVldXV1Wk0GldXV39//9DQ0B49ejh6urYpLi7etm2b1WoVQowf\nP97Hx8fREwEAgMvb5RTsKisrd+/ebTQamxarqqqqqqpyc3NDQkKuvvpqjUbjqPHaxGKxHDx4\n0JbqAAAA2sVlE+zq6+t37NhRVVUlhDAYDH379vX09DSbzSUlJTk5OQ0NDbm5uV5eXoMHD3b0\npK1y7Ngxo9Ho7OxcX1/v6FkAAIBKXDbB7uTJk7ZU5+fnFxcXp9P9a/K+ffsOGDBg27Zt7u7u\nFoulsX/fvn25ublarXbatGn79+8/e/ZsVFTUgAEDhBAWiyUzMzMnJ8doNFqtVjc3t549e0ZE\nRLi5udm23bZt2/nz5zUazZ133tl4h999953JZPL29p4wYYIQYu/evXl5eUKIqVOnpqWl5eXl\n1dbWuru7R0RE9O3bt/nHUl5enp6ertPpwsLC0tLS2vPXBAAArmCXTbDLycmx3Rg+fHhjqrPx\n9PScMmWK4iCsk5OTEMJisaSlpeXn5wshGhoahBBms3nnzp1FRUW2No1GU1lZefLkydzc3Li4\nOG9v71bOY7t/IcSvv/5qNBo9PDxqa2srKioOHjwohGgm21mt1oMHD1oslqFDh7ZyXwAAAK1x\neQS7+vr6iooKIYSbm1tAQIDc0MxH6zIzM8PCwvz8/Dw9PYUQx44ds6W6kJCQmJgYnU6XlZWV\nnJxcV1d34MCBcePGtXKkxj3W19dPnDhRr9eXlpb+/PPPFovl8OHDoaGhFxrpxIkTJSUlfn5+\n4eHhJ0+ebOXuAAAAWnR5BLu6ujrbDXd396b1PXv22FbjGl133XW9evVqWundu3dMTEzjj1lZ\nWUIIvV4/cuRIrVYrhAgLCztz5kxBQUFpaWl5eXnrF+1sIiIi9Hq9EMLX17dHjx5nzpypra01\nGo1276eqqurIkSNarfaqq666XM7zAAAAl4vL4zp2jRmo6afoWql3796NtysrK20Z0c/Pz5bq\nbPz8/Gw3ysrK2nr/TS9TYlsUFEJUV1fbbT548KDZbI6IiGhrfAQAAGjR5bFi5+LiYrtRWVnZ\ntD58+PAhQ4YIIQoLC1NSUuxu23hKhBCi8RRUZ2fnpj2NP17EOapNP/BnW7oTQphMJrkzOzu7\nsLDQ09Pzcjl1FwAAXF4uj2Cn0+m8vLyMRqPJZDp//nzjx+wMBoPthiLwNdV0Za4xwCmCV2Oe\na0xmQgir1WqxWGybWyyWC2U+k8nUOEbj3TYm0aZsR40rKiq+/vprxT9t3bpVCHHnnXdyfBYA\nAFy0y+NQrBAiNDTUdiMpKanxI3c2Vqu1sLCwNXfi4eFhi1zFxcW2k2Rtzp07Z7thOybbGO8a\nL4ZcUFBwoYsJnz9/Xr7dGPUAAAA6zeWxYieECA8Pz87OrqioKC8v//HHH/v16+ft7W2xWIxG\nY15enu0Sd1qttumBV7v69euXlpbW0NCQnJwcExOj1WozMjKKi4uFEN27d/fw8BBCeHl5nTlz\nRgiRkpISHR1dV1eXkpKi0+maZsFGaWlpAQEB3t7eOTk5tvsxGAxeXl5y5zXXXGM2m5tWsrKy\njhw5IoQYM2aMt7c3y3UAAOBSXDbBzsnJacyYMbt377Z9S6x8XV9vb+/Y2NjG0yAuZPDgwefP\nny8qKjp16pTt2ni2pTiDwTBixAhbT79+/U6cOGGxWIqKimwHSfv3719cXGz31Ap/f/+tW7dq\ntdrGEzuGDRtmd9fOzs6Kz/Y1fj7PxcXF1dW1+ckBAACad9kEOyGEwWAYN27cqVOn8vPzS0tL\nTSaTTqdzdXX19/cPCgoKCgpqzYqXLSBmZGTk5ubavnnC3d29V69eAwcObDwC6+HhERcXl5yc\nbDQaXVxcQkNDIyMjf/nlF/Hvqxw3FRUV5eXllZ2dXVdXZ/tOs6bn4QIAAHQaDd9Df3EOHjyY\nnZ0thLjllltsB3ABAAAc67I5eQIAAADNI9gBAACoBMEOAABAJfiMHQAAgEq0dsWuvr7++eef\nd3JyarwmSPNKS0sXLFjQp08fFxeXvn37Tp06dd++fZcwJwAAAFrQqsudpKWlzZgx4+TJk628\n05KSktjY2FOnTk2aNOnBBx/Myspav379Tz/9dODAgaioqEuYFgAAABfU8oqd0WiMjY3VarXJ\nycmK6+teyOLFi0+dOvXhhx/+8MMPS5cuXbt27fr162tra5977rlLHhgAAAD2tRzsGhoa5s6d\nu2fPnrCwsFbeqbOzc3x8/Jw5cxor06ZNc3NzO3r06EWOCQAAgJa07eQJV1fXyMjI3377ra27\nqaur8/T0HDly5O7du9u6LQAAAFqjk75S7OOPP66vr7/77rsV9ePHj5eXlzf+GBER4eXl1Tkj\nAQAAqExnBLsdO3YsXLhw9OjRjz76qOKfVq5cuWfPnsYfKyoqtmzZ0q9fv06YCgAAQGU6PNit\nW7fuD3/4Q2Rk5MaNG3U65e4GDhzY0NBgu3306NG0tLT6+vqOHgkAAECVOjDYWa3WJUuWLF26\n9Oabb/5//+//eXp6yj3z5s1rvD1z5sydO3d23DwAAADq1lHBzmq1zp49e/Xq1U888cS7777r\n5OTUQTsCAACATUd9V+zTTz+9evXqV1999YMPPiDVAQAAdIJ2CHa1tbWpqamZmZmNlQ0bNrz/\n/vvz589//vnnL/3+AQAA0BotH4rdsWPHli1bbLcbGhpOnz7d+AUSCxcu9Pf3z8jIiI6Ojo+P\nT0xMtNWfeeYZIYTFYpG/auLZZ5/19fVtt/EBAAAurN5iMVksHb0XvVbrrO2oo6Bt0nKw27t3\n7xtvvNH4Y0FBQeOPs2fP9vf3lzexrd59+OGH8j89+uijBDsAANA5Ek8XfH4is+W+S3NfWN8p\nfXp39F5ao+Vg99xzzzX/Ha+RkZGKr69o07dZAAAAdJAAV5coP5/W95eaTPmV1b3c3fxcXFq/\nVTc317aP1iE66ZsnAAAAOt9Vgf5XBdo5ungh28+e++jYiVuCe93Uq2fHTdVxusTxYAAAAFw6\ngh0AAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0A\nAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQAghzFZrbmWVECKvqqrBYnX0OBdD\n5+gBAAAAHO90VfXbh9NOV1ULIX7MO3OouGzB0MG93Q2OnqttWLEDAABXOrPV+v6RdFuqE0II\noTlbXfPu4TSz9TJbtyPYAQCAK12msSKnskpRzK+qPl5mdMg8F41gBwAArnRGU73dernJ1MmT\nXCKCHQAAuNJ1d3OValYhRE+DW+cPcykIdgAA4EoX7OE+MtD/v0pWTWyAXx9PDwdNdJEIdgAA\nAOLRwQOu79FNI4QQQiPE6J7d5g4eqHHwUG3G5U4AAACEu073+JCB4T5eq9Mz7h/Qb1JwL0dP\ndDFYsQMAAPgXF6228f9eji7XuQEAAKBAsAMAAFAJgh0AAIBKEOwAAABUgmAHAACgEgQ7AAAA\nlWhVsCsrK3vqqadCQ0P1en1QUNDs2bPPnj3b/Cbp6en3339/z549nZ2dAwMDp02bduDAgfYY\nGAAAAPa1fIFik8kUHx+fnJyckJAQExOTmZm5Zs2abdu2JSUl+fr62t3k6NGj1157rbOz8+OP\nPx4WFpaTk7Ny5crrrrvup59+uvHGG9v7IQAAAECI1gS7lStXJicnv/HGG88884ytMmHChOnT\np7/yyitvvfWW3U1effXVioqKbdu2jR071laZMmXKsGHDXn75ZYIdAABAB2n5UOyaNWs8PT3n\nz5/fWLnrrrvCwsK+/PJLq9Vqd5PMzEwhxOjRoxsrQ4cO9fLyOnXq1KXOCwAAgAtoIdjV1tYe\nPnx45MiRLi4uTeujR48uLCzMzs62u1VERIQQ4vjx442V8+fPV1ZWDho06JIHBgAAgH0tBLu8\nvDyz2RwcHKyo9+nTRwiRlZVld6tnn33W19d3xowZu3fvLigoSElJufvuu11dXV966aV2GRoA\nAACyFj5jV1FRIYRwd3dX1D08PBr/VTZo0KC9e/fefvvt119/va0SEhKSmJh49dVXKzqfffbZ\ngwcP2m5XVVXZ7hYAAAAXoVWXO9FoNIqK7dN1ct0mLS1t4sSJFRUVb7/99vfff/+3v/3N09Nz\n4sSJiYmJis6amhrjv5nN5gvdIQAAAFrUwoqdl5eXsLcyZzQahRCenp52t5o1a9a5c+dOnDjR\nq1cvW+Xuu+8eMGDAzJkzs7OznZ2dGzs/+OCDxtszZ85MSkpq+0MAAACwL/F0wf9m2D8lwK4G\ni1UIseZE1v9mnGr9Vnf163NzcFBbZ+sILQS7kJAQnU6Xk5OjqNvOew0PD5c3qays3L9/f1xc\nXGOqE0IYDIb4+Pg1a9acOHFiyJAhlzw2AABAyww6p+5urh2+F+eWrx/XOVqYQ6/Xx8bGHjhw\noLq62mAw2IoWi2XHjh3BwcEhISHyJjU1NVartba2VlG3VeQ6AABABxnVPXBU90BHT9F5Wv6M\n3UMPPVRdXb18+fLGyqpVq86cOTN79mzbj7W1tampqbY1PCFEYGBg3759f/vttxMnTjRuUlZW\nlpiY6OXlFRkZ2a7zAwAA4F9aXjmcNWvWl19+uWTJkpSUlJiYmLS0tPXr10dFRS1YsMDWkJGR\nER0dHR8f33huxNtvv33HHXeMGjXq0Ucf7d+//9mzZz/99NOSkpKVK1cqrocHAACA9tJysHNy\nctq8efNf/vKXf/zjH5s3b+7WrdvcuXOXLl3aeGRWNm3atN27d7+jzMKHAAAfWUlEQVT55pur\nVq0qLS319PSMjY1dsWLFLbfc0q7DAwAA4D80F/pasM43c+bML774Ij09feDAgY6eBQAA4PLT\nquvYAQAAoOsj2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIE\nOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAA\nAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwAAoBIEOwAAAJUg\n2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVKJVwa6srOypp54KDQ3V6/VB\nQUGzZ88+e/Zs6/fxxz/+UaPRzJ49+2KHBAAAQMt0LXaYTKb4+Pjk5OSEhISYmJjMzMw1a9Zs\n27YtKSnJ19e3xc1/++23Dz74oD1GBQAAQHNaXrFbuXJlcnLyG2+88fXXX7/wwgt/+9vf1q5d\nm52d/corr7S4bUNDw8MPPxwZGdkeowIAAKA5LQe7NWvWeHp6zp8/v7Fy1113hYWFffnll1ar\ntflt33777UOHDr3++uuXOiYAAABa0kKwq62tPXz48MiRI11cXJrWR48eXVhYmJ2d3cy2mZmZ\nf/nLXx599NFrrrmmHSYFAABAs1oIdnl5eWazOTg4WFHv06ePECIrK6uZbefMmePj4/Paa69d\n4ogAAABojRZOnqioqBBCuLu7K+oeHh6N/2rX559//vPPP3/99dfe3t5lZWUXalu9evWJEyds\nt/Py8tzc3Fo5NwAAABRaPitWCKHRaBQV26fr5LpNYWHhn/70p8mTJyckJDR/z6mpqXv27PnP\nNLpWzQMAAABZC0HKy8tL2FuZMxqNQghPT0+7W82fP99kMq1cubLF3c+bN2/GjBm228uXLz90\n6FCLmwAAAMCuFoJdSEiITqfLyclR1DMzM4UQ4eHh8iZbtmz5+9//vmjRIq1Wm5+fL/6dAqur\nq/Pz8728vGxh0WbgwIGNt93d3RsaGi72gQAAAFzpWjh5Qq/Xx8bGHjhwoLq6urFosVh27NgR\nHBwcEhIib/Lzzz8LIV5++eXgfxsyZIgQYt26dcHBwa+++mq7zg8AAIB/afkzbQ899NAjjzyy\nfPnyl156yVZZtWrVmTNn/vKXv9h+rK2tTU9P9/T07N+/v60/Li6u6T1UVVXdfffd48ePf+KJ\nJ8LCwtr5EQAAAEAI0ZpgN2vWrC+//HLJkiUpKSkxMTFpaWnr16+PiopasGCBrSEjIyM6Ojo+\nPj4xMVEIMWjQoEGDBjW9B9tZscHBwZMnT+6AhwAAAAAhWvPNE05OTps3b16wYEFqauqyZct2\n7do1d+7c7du3GwyGTpgPAAAAraRp8WvBOs3MmTO/+OKL9PT0pmdUAAAAoJVaXrEDAADAZYFg\nBwAAoBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVKLl74oF\nAHWwVBhNudkdvRcnP3/nnr07ei8AYBfBDsCVov7s6fL1azp6L25XjXK+lWAHwDEIdgCuFLqA\nbh7jJrW+31xaUvPbXn1oP334oNZv5dyzV9tHA4D2QbADcKVw8vN3v/7G1vebTmXV/LbXuXef\nNm0FAA7EyRMAAAAqQbADAABQCYIdAACASlyhn7EzlxSXfvFxR+/FuU9f79vv6ei9AAAA2LBi\n1zoWi7m02FxZ4eg5AAAALugKXbFz8vMPePqF1vdbjOVFby11CY/wufvBjpsKAADgUrBiBwAA\noBIEOwAAAJUg2AEAAKgEwQ4AAEAlCHYAAAAqQbADAABQCYIdAACAShDsAAAAVIJgBwB2WGtr\nao+kCiFMOdkNReccPQ4AtArBDgCUGorOnX//9ZoDvwoh6vNOlXz0Tk3KAUcPBQAtI9gBgJLx\n6/+1VFU2/mhtaKj44VtzWYkDRwKA1iDYAcB/MZeV1p/NVxSt9aa6E2kOmQcAWo9gBwD/xWqq\nu0Dd1MmTAEBbtSrYlZWVPfXUU6GhoXq9PigoaPbs2WfPnm33TQCgK3DyC9C4uMh15569On8Y\nAGgTXYsdJpMpPj4+OTk5ISEhJiYmMzNzzZo127ZtS0pK8vX1ba9NAKCL0Oh0nuMnG7//pmnR\nJWKIvl+4o0YCgFZqOditXLkyOTn5jTfeeOaZZ2yVCRMmTJ8+/ZVXXnnrrbfaaxMA6Drcrhql\ncdZX/vyTubxEo9cbrh7tPuYmodE4ei4AaIHGarU23xEdHZ2ZmVlUVOTS5NhEeHi40WgsKCjQ\n2Humu4hNhBAzZ8784osv0tPTBw4ceFGPpaPU52RX/vMHU+4podcbRlzrETde4+rq6KEAdDjT\nqazS1SvdR4/1GD/Z0bMAQKu08Bm72traw4cPjxw50uW/P3EyevTowsLC7OzsdtmkKzOdyir5\n2wpT7ikhhDCZqvfsKF37ibBYHDwWAACApIVDsXl5eWazOTg4WFHv06ePECIrK6tfv36Xssme\nPXsKCgpst0tKSpydnX/88cfff/+96YZXXXVVaGio4t6OHj167NgxRbFfv36xsbGKYkFBwa5d\nuxRFb2/v8ePHK4r19fXfffedoljx0/fje/g7a/+TgOtzT9Wk/razqNRoNCqar7/++h49eiiK\nSUlJWVlZiuLgwYOHDBmiKJ46dergwYOKYvfu3W+44QblVBUVP/74o6Lo5OR0++23C8nGjRtN\n0tl848aN8/HxURR//fXXM2fOKIoxMTH9+/dXFNPT0w8fPqwo9unTZ+TIkYpiUVHR9u3bFUUP\nD4+JEycqihaL5ZtvvhGSW2+91VVaIk1MTCwtLVUUR40a1auX8uPtqampJ0+eVBQHDhw4dOhQ\nRTEvL2/fvn2KYkBAwNixYxXF6urqTZs2yaPeeeedcvGHH36oqalRFG+88UZ/f39Fcd++fXl5\neYrisGHDBgwYoCiePHkyNTVVUezdu/e1116rKJaUlPz888+Koqur66233iqP+vXXX8tL+Lfc\ncou7u7uiuH379qKiIkXx6quvDgkJURQPHz6cnp6uKIaFhUVHRyuKZ86c+fXXXxVFHx+fcePG\nKYp1dXX/93//J89/++23Ozk5KYo//vhjRUWFojhmzJhu3bopigcPHjx16lTTSkPRud5FxcpB\nhcjKykpKSlIUe/Tocf311yuK5eXlW7duVRSdnZ2nTp0qz//dd9/V19crihMmTPDy8lIUd+3a\n1fjM2Sg2NlZ+Qj527NjRo0cVxb59+44YMUJRPHfu3M6dOxVFT0/Pm2++WVFsaGj49ttv5flv\nu+02vV6vKP7zn/8sKytTFEePHt2zZ09FMTk5OTMzU1EcNGhQZGSkopiTk3PggPJ60YGBgXFx\ncYpiZWXlli1bFEWtVpuQkCDP//3339fW1iqKN910k/zR8D179pw+fVpRHD58eHi48lOYx48f\nV7yiCSGCg4OvueYaRfH8+fO//PKLomgwGCZNmiSP+o9//EMuTp482c3NTVHctm1bcXGxonjN\nNdfIr9GHDh06ceKEohgeHj58+HBFMT8/f+/evYqin59ffHy8olhbW/v999/Lo95xxx3ysbvN\nmzdXVVUpinFxcYGBgYri/v37c3NzFcWoqKiIiAhFMSMjIyUlRVEMCgq67rrrFMWysrJ//vOf\niqKLi8uUKVPk+Tds2GA2mxXFm2++2dPTU1HcsWNHYWGhomg31Rw5ciQtTXlBJbup5uzZs926\ndZOf6/7D2izbk9e8efMU9eXLl9se2yVu8sQTT8Q2If9ShBCfffaZvJdFixbJnY888ojcuXnz\nZrkzMjJS7pSffWxOPvN4waI/Nv1f+f/9Q45lQogtW7bId/vwww/LnYsXL5Y7V69eLXfedNNN\ncqf8SimEMBgMcqfVapUDhBAiJSVF7rT7Yr9y5Uq589VXX5U7Z8yYIXfKT1VCiP79+8uddXX2\nrzFRUFAgN8svS+IC/w/55JNPyp0LFy6UO9etWyd3XnfddXJnTk6O3KnRaOROq9Xau3dvuXnP\nnj1y51133SV3vvXWW3Lnu+++K3cmJCTInfv375c7e/bsaXdUnc7OO72srCy5U36zIYRYu3at\n3Pncc8/JnfLzg9Vq3bhxo9wZHR0td8qZ0qa6ulpulmOxECIxMVHunDlzptz5TNyoip++V3R+\n/PHHcufEiRPl+5Tf/wghvL295U6r1SoHOCHEkSNH5E45bAkhVq1aJXcuWbJE7pw1a5bcKQdQ\nIURERITcWVlZKXcKIYqLi+XmYcOGyZ3ff6/8lVqt1scee0zufOGFF+TONWvWyJ1xcXFyZ0ZG\nhtyp1+vlTqvV2r17d7n54MGDcue0adPkzvfff1/ufPPNN+XOu+++W+7cvXu33BkSEiJ3Wi5w\nyCgvL09ult/sCSHWr18vd/7pT3+SO5966im58+uvv5Y7R44cKXfKKwU29fX1cnPfvn3lzh07\ndsid9913n9z52muvyZ0rVqyQO6dMmSJ3Jicny50BAQFyp9VqlQO0EOL48eNy54033ih3fv75\n53Lniy++KHfOmTNH7ty0aZPd57pGLZ88IYSQk7XVarVbb+smN910U+NbnE2bNh05cmTp0qWK\nv65Ro0bJ9z9lyhT59XLQoEFyZ2RkpPws7OfnJ3e6ubkpO61W46YNrjplLtbonF966SV5xchu\n2psxY4acQmJiYuTOUaNGyaPKS1BCiO7du8uddl+VhRDvvPOO/DZUfrsmhJg3b97kycrPEo0e\nPVrunDhxopwX5XerQogBAwbIo9p9AdPpdHZfL+02v/DCC/Kru7wIJISYPn26/N/F7ovNVVdd\nJQ8gL8EKIfz8/OyOatfrr78uvw2VV1aEEA8//LD8lld+Zy+EiI+Plwew+7QYGhoqdxoMBruj\nfvTRR/JrRkBAgNy5YMEC+bn16quvljunTZsmD2b3L2XYsGHyqHb37uHhYff37+zsLBeXLl1a\nXl6uKMrv7IUQDz74oOJVsKH4/IAsOxclvv766+UB7P5NBQUFyZ3yspbN+++/Ly+uBwUFyZ1P\nPvmknC3k9UIhxOTJk+W1MbufYx48eLA8qryuL4RwcXGx+/uXV3aFEIsWLZJXjOT1ciHEvffe\nKy8O2f2jvuaaa+QB5IcphAgMDJQ7tVr7n0Favny5vLhuO9ak8Oijj8rZ2u5L1fjx4729vRVF\n+RiIrSiP6uHhIXdqNBq7v3+7F5149tlnz51TftOx3TfGd9xxh/wuKCoqSu6MiYmRB5CXwIUQ\n3t7edke1+5/glVdekRfX7b6szJo1S35vKR8vEkLExcXJA9j9bxocHCx3ygeLbFasWNHQ0KAo\n2n1j8PTTT0+fPl1RtJu2p0yZIj+H2E01UVFRdp/rGrVw8kRGRkZ4ePiDDz74+eefN60vWrRo\n2bJliYmJ8uvQRWxi0zVPnij739V16coDGb4zH9P3C3PIPAA6DSdPALjstHDyREhIiE6nkw88\n2T4JYTdKX8QmXZnX5ASt13+93zKMuoFUBwAAuqAWgp1er4+NjT1w4EB1dXVj0WKx7NixIzg4\nWP6g9MVt0pVpvbwDnnjW/YabhBBaX3/fB+d43nybo4cCAACwo+WvFHvooYeqq6ttpz7YrFq1\n6syZM7Nnz7b9WFtbm5qa2vRsphY3ubxoXFwMI0cJIZx79tL3t/NBbAAAgK6g5ZMnZs2a9eWX\nXy5ZsiQlJSUmJiYtLW39+vVRUVELFiywNWRkZERHR8fHxycmJrZyEwAAALS7llfsnJycNm/e\nvGDBgtTU1GXLlu3atWvu3Lnbt2+/0Il1F7cJAAAALlGrLnfi4eGxfPnypodWm7JdE65NmwAA\nAKDdtbxiBwAAgMsCwQ4AAEAlCHYAAAAq0arP2AGAGpjNFuk7u5phqasTQlgbGizSN001Q6Nz\n0jjb/9IwAOhoBDsAV4q6zBNlaz9t61bV+3ZV79vV+n63q0Z53ZrQ1r0AQLsg2AG4Umg9PV2H\nDOvovTgH9eroXQDAhRDsAFwpnHv29p7+gKOnAIAOxMkTAAAAKkGwAwAAUAmCHQAAgEoQ7AAA\nAFSCYAcAAKASV+pZsWaz2VjW+nZLZYUQwlpvMpcWt34rjc5Z6+nV5tkAAAAuyhUa7BpKi4s/\neKOtW5lOpp9/99XW9+v7D/B9cE5b9wIAAHBxrtBgp3V1cxtxTUfvRRfYvaN3AQAA0OhKDXYe\nnl5T7nT0FAAAAO2JkycAAABUgmAHAACgEgQ7AAAAlSDYAQAAqESXO3nik08+CQwMdPQUAAAA\nXdEf//hHZ2fnC/6ztcv45ZdfIiIiOvE30wY+Pj7Dhw8fMmSIowcB0Hl69uw5fPjwgQMHOnoQ\nAJ2nT58+w4cP79u3r6MHuaCqqqpm0pTGarU6esL/2LFjR2FhoaOnsGP//v3bt2/XaDTPPPOM\no2cB0Em++eabjIwMvV7/9NNPO3oWAJ3kk08+KSkp8fX1feSRRxw9i3233367k5PThf61ax2K\nHTNmjKNHsK+mpmb79u1arfbOO7n6HXCl2LdvX0ZGhrOzM3/4wJVj/fr1JSUlHh4el+kfPidP\nAAAAqATBDgAAQCW61qHYLmvEiBFDhgwJCgpy9CAAOs/NN998/vz5qKgoRw8CoPNMmTLl559/\nHjt2rKMHuUhd6+QJAAAAXDQOxQIAAKgEwQ4AAEAlCHYtSE1N7dGjh0aj0Wg0Wq3W1dX1hRde\ncPRQADrWunXrfH19tVqt7W/f1dV18eLFjh4KQOfp3r27RqPx8fFx9CBtRrBrTlJSUkxMzLlz\n5/R6fUhIiMFgqKure+211959911Hjwago3z66af33ntvWVmZu7t7nz593N3d6+rqXn755eef\nf97RowHoDEuXLu2aX5fQGpw80ZyAgIDi4uLo6Ojk5GRbZerUqRs3btTr9XV1dY6dDUAH0ev1\n9fX18+fPf++992yV22+//dtvvzUYDFVVVY6dDUBHq6ys9Pb2FkJYLBZvb++ysjJHT9Q2BLvm\ndO/evbS0tKioyPbfWAhhNpt1Op0Qgt8boFYhISEWiyU/P7+xUlNTYzAYtFqt2Wx24GAAOkFE\nRMTx48cfeOCBNWvWXI7BjkOxzTl37pzJZGpMdUKIkpISIUQz39EG4HKXm5vbNNUJIbZv3y6E\n0Ov1jhkIQGf56quvjh8/7ufnt2DBAkfPcpEIdm0zbtw4IcTQoUMdPQiAzpCXl/f4449PmjRJ\nCMFn7ADV+8Mf/iCE2LVrl6MHuXh880QbPPnkk4cOHXJyctq5c6ejZwHQ4TQaje2Gh4fHqlWr\n7rnnHsfOA6BDjRkzpr6+fsqUKYMHDz58+LCjx7lIBLvWmjBhwtatW7Va7f79+z08PBw9DoAO\nFxQUVF1dXVFRUVlZef/99wshyHaAWu3evXvnzp0uLi4bN2509CyXhEOxLTObzaGhoVu3btXr\n9enp6bGxsY6eCEBnOH36dGlpaUNDw7x588xm83333WcymRw9FIAOMWHCBCHE+vXrHT3IpSLY\ntcBsNvv7++fk5AQEBJSVlYWHhzt6IgCdbcWKFe7u7lardfXq1Y6eBUD7e+CBB6qrq0NCQnQ6\n3aZNmzZt2mT7mF1DQ8OmTZvS09MdPWAbcCi2BT179iwvLw8PDz9x4oSjZwHQ4bZs2TJp0iR3\nd/eKioqmdduFTi7fa5YCaMbWrVuFELm5uZMnT25ar6qqmjx5clBQ0OnTpx00WpsR7JozderU\noqKiwMBAUh1whZg4caLVaq2srHzxxReXLVtmK65evbq2tlYIMWfOHIdOB6BDLFmyZMeOHU0r\n5eXlW7ZscXZ2TkhIGDt2rKMGuwhcoLg5Wq3WarUGBATI16/atGnT8OHDHTIVgA6VkJCwYcMG\nIYTBYPDz8zMajUajUQgRExOTlJTk6OkAdIbDhw8PHTqUCxSrjS31nj9//owkIyPD0dMB6BDf\nfPPNn//8Z1dX1+rq6vz8fKPR6OzsPH36dFIdgK6PFTsAAACVYMUOAABAJQh2AAAAKkGwAwAA\nUAmCHQAAgEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmC\nHQAAgEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAA\ngEoQ7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ\n7AAAAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ7AAA\nAFSCYAcAAKASBDsAAACVINgBAACoBMEOAABAJQh2AAAAKkGwAwAAUAmCHQAAgEoQ7ABc3qZO\nnappIjAwMD4+fufOnY6e61Kp9XEB6FAEOwCdy2o1l5dZa2vb8S779+//y7998MEHDQ0NY8eO\nPXjwYDvuojVK6kxVDQ3teIft/rhWrlw5c+bM9hsQQJejc/QAAK4gNUn7KhO3WKoqhRD6vv09\nJyfoArtf+t16eHjExcU1/jht2rSwsLD3339/7dq1l37nrbG7oHBtRnZpnUkIMcDb66GBYaGe\n7pd+t+3+uJKSki59KgBdGSt2ADpJ7eEU48Z/2FKdEMKUnVm25hNLTU2778jV1XXo0KEZGRm2\nHwMDA997771bbrnF1dW1vLy8oaFhyZIlERERbm5uAwYM+Oijjy5xdynFJR8ePW5LdUKIE+XG\nV1MPl5lMl3i3MsXjqqurW7hwYXBwsF6v79Onz5///OeGf68X7tq164YbbvDx8fH09Bw9erTt\nAG5cXNxnn332xRdfaDSa1NTUdh8PQFfAih2AdmYuOW+xd6S14qfvlZ3lpZXbfnSLvkpudvL1\n17q5XfQM2dnZQ4YMsd3W6/WrVq269dZbFy1a5O7uvnDhwlWrVn300UejRo1KTEycP3++Xq9/\n6KGHWrzPwpraSntHWteczFZUyk316zJPTegdJDcHurp4Oju3/QH9S9PHNXfu3O++++6vf/3r\niBEj9u3b99hjj9XU1LzzzjtVVVWTJ0++5557Pv74Y6vVumLFiokTJ+bn52/cuDE+Pj48PPzD\nDz/09fW96BkAdGUEOwDtrGLTt3Un01vZXLN/d83+3XLd5+4HXQYPbf1OGxerzp079+GHH6an\np7///vu2ikajMRgMb7zxhhDCaDT+9a9/ff755x944AEhRFhYWFJS0uuvv96aYLcu89Sec0Wt\nnGf7mXPbz5yT648OGjA2qA1Hny/0uIqLi9esWbN8+fLp06cLIfr375+Wlvbee++9/vrrubm5\nRqNxxowZgwYNEkJ88MEH06dPd3FxMRgMOp3OxcUlICCg9QMAuLwQ7AC0M5fIYboedharqvfu\ntEorXs69Q/R9w+Rmp4Burd/joUOHnJssg/n6+q5evXr8+PGNlWuvvdZ2IzU11WQyNf2nuLi4\nTz/9tLKy0sPDo/m9jAj0D3R1ketb88/WmM2KYh8P9+H+dlbF+ni04bN3zTyuQ4cONTQ0XHPN\nNf8Zb8SIqqqqkydPRkREDBw48L777nvsscfGjx8fHR09ZsyY1u8UwGWNYAegnblFj7Rbt9ab\nqvf91+KcRqfzuvUOXc9el7jHAQMGfPXVV7bbBoMhPDzc+b8Pd3p7e9tuGI1GIcTYsWM1Go2t\nYrFYhBAFBQVhYXbyZVPXdQ+8rnugXHfSajdk5yqKswb2j/DxbvMj+W/NPC7bA/Hy8mps9vT0\nFEJUVFQ4OTnt2rXrzTff/OSTT55//vmQkJBly5bdf//9lzgMgMsCwQ5AJ/EYN6mhsMCU9a/P\n/mt0Os+Jt116qhNCuLm5jRgxojWdtoS3du3aqKiopvXg4OCL3vsdfUNyKiqTzpcIIYRVODtp\np/fvc+mpTjT7uGwPxBbvbGy3bfXAwMDly5cvX7782LFjb7/99gMPPDB48ODY2NhLHwlAF0ew\nA9BJNM5635mPmbJO1p/O07oZ9GEDnHz8OnmGYcOGubi4FBYWRkRE2CpFRUVardbFxc4x1lZy\n0mieGTYkraz8ZHmFm84p0tenp+HiT/topWHDhul0ul9//bXxaOzevXu9vb3Dw8Ozs7N///33\n2267TQgxePDg//mf//niiy+OHDliC3ZWq7WjZwPgQAQ7AJ1K3y9c3y/cUXv38vJ65JFHXnrp\npYCAgJEjR+bk5Dz11FO9e/f+4YcfLvGeB/l4D2qPVbpW8vPzmzVr1muvvRYWFhYdHb19+/aV\nK1cuXLhQp9Pl5uYmJCS88cYbkydP1mg0X331lVartX3K0NfXNyUlJTU1NTg42N/fv9OmBdBp\nCHYArizvvPOOj4/PM888c/bs2R49ekyZMuXVV1919FAX48MPP/T09Jw7d25hYWFwcPCLL774\n3HPPCSHGjBnz2Wefvf3224sXL9bpdEOGDPn2228HDBgghHjyySfvv//+0aNHf/PNNxMmTHD0\nIwDQ/jQsywMAAKgD3zwBAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0A\nAIBKEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBK\nEOwAAABUgmAHAACgEgQ7AAAAlSDYAQAAqATBDgAAQCUIdgAAACpBsAMAAFAJgh0AAIBKEOwA\nAABUgmAHAACgEgQ7AAAAlSDYAQAAqMT/B35v0jtA7GExAAAAAElFTkSuQmCC",
"text/plain": [
"plot without title"
]
},
"metadata": {
"image/png": {
"height": 420,
"width": 420
}
},
"output_type": "display_data"
}
],
"source": [
"# Plot group-time average treatment effects\n",
"ggdid(example_attgt_dml_linear)"
]
},
{
"cell_type": "markdown",
"id": "d4e701ec",
"metadata": {},
"source": [
"It's also possible to calculate aggregated effect estimates. Please note that, the results are again very close to those in [the original notebook](https://bcallaway11.github.io/did/articles/did-basics.html#aggregating-group-time-average-treatment-effects)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "8f163147",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:44:05.707983Z",
"iopub.status.busy": "2024-08-31T09:44:05.707126Z",
"iopub.status.idle": "2024-08-31T09:44:05.777325Z",
"shell.execute_reply": "2024-08-31T09:44:05.776247Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"aggte(MP = example_attgt_dml_linear, type = \"simple\")\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"\n",
" ATT Std. Error [ 95% Conf. Int.] \n",
" 1.6586 0.0368 1.5864 1.7308 *\n",
"\n",
"\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"Control Group: Never Treated, Anticipation Periods: 0\n"
]
}
],
"source": [
"agg.simple <- aggte(example_attgt_dml_linear, type = \"simple\")\n",
"summary(agg.simple)"
]
},
{
"cell_type": "markdown",
"id": "fb050302",
"metadata": {},
"source": [
"### Details on Predictive Performance\n",
"\n",
"We can add an evaluation functionality to the wrapper to assess how the predictive performance from the linear and logistic regression differ from that of the random forest learner."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "a8bb468f",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:44:05.780722Z",
"iopub.status.busy": "2024-08-31T09:44:05.779819Z",
"iopub.status.idle": "2024-08-31T09:44:05.800366Z",
"shell.execute_reply": "2024-08-31T09:44:05.799324Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"In order to avoid name clashes, do not attach 'mlr3measures'. Instead, only load the namespace with `requireNamespace(\"mlrmeasures\")` and access the measures directly via `::`, e.g. `mlr3measures::auc()`.\n",
"\n"
]
}
],
"source": [
"library(mlr3measures)\n",
"\n",
"# Add a wrapper that computes the RMSE and accuracy for the nuisance components, can be customized by providing custom measures\n",
"eval_preds = function(y, d, predictions, params_names, custom_measures = NULL) {\n",
" measures_res = list()\n",
"\n",
" if (!is.null(custom_measures)) {\n",
" # Alternatively provide a named list with custom evaluation functions\n",
" measure_funcs = list()\n",
" } else {\n",
" measure_funcs = list()\n",
" measure_funcs[['ml_m']] = mlr3measures::acc\n",
" measure_funcs[['ml_g0']] = mlr3measures::rmse\n",
" }\n",
"\n",
" for (param_name in params_names) {\n",
" preds = predictions[[param_name]][, 1, 1]\n",
"\n",
" if (param_name == \"ml_m\") {\n",
" obs = d\n",
" # map probability predictions to binary\n",
" preds = as.factor(ifelse(preds > 0.5, 1, 0))\n",
" obs = as.factor(preds)\n",
" }\n",
"\n",
" else if (param_name == \"ml_g0\") {\n",
" obs = y[d == 0]\n",
" preds = preds[d == 0]\n",
" }\n",
"\n",
" if (param_name == \"ml_g1\") {\n",
" next\n",
" }\n",
"\n",
" else {\n",
" measure_func = measure_funcs[[param_name]]\n",
" measure_pred = measure_func(obs, preds)\n",
"\n",
" measures_res[[param_name]] = measure_pred\n",
" }\n",
"\n",
" }\n",
" return(measures_res)\n",
"}\n",
"\n",
"# evaluate learner performance: linear models\n",
"doubleml_did_eval_linear <- function(y1, y0, D, covariates,\n",
" ml_g = lrn(\"regr.lm\"),\n",
" ml_m = lrn(\"classif.log_reg\"),\n",
" n_folds = 10, n_rep = 1, ...) {\n",
"\n",
" # warning if n_rep > 1 to handle mapping from psi to inf.func\n",
" if (n_rep > 1) {\n",
" warning(\"n_rep > 1 is not supported.\")\n",
" }\n",
" # Compute difference in outcomes\n",
" delta_y <- y1 - y0\n",
" # Prepare data backend\n",
" dml_data = DoubleML::double_ml_data_from_matrix(X = covariates, y = delta_y, d = D)\n",
" # Compute the ATT\n",
" dml_obj = DoubleML::DoubleMLIRM$new(dml_data, ml_g = ml_g, ml_m = ml_m, score = \"ATTE\", n_folds = n_folds)\n",
" dml_obj$fit(store_predictions = TRUE)\n",
" att = dml_obj$coef[1]\n",
" # Return results\n",
" inf.func <- dml_obj$psi[, 1, 1]\n",
"\n",
" # Evaluate learner performance\n",
" predictions = dml_obj$predictions\n",
" params_names = dml_obj$params_names()\n",
" eval_predictions = eval_preds(delta_y, D, predictions, params_names)\n",
" print(eval_predictions)\n",
"\n",
" output <- list(ATT = att, att.inf.func = inf.func)\n",
" return(output)\n",
" }\n",
"\n",
" library(mlr3measures)\n",
"\n",
"# evaluate learner performance: random forest\n",
"doubleml_did_eval_rf <- function(y1, y0, D, covariates,\n",
" ml_g = lrn(\"regr.ranger\"),\n",
" ml_m = lrn(\"classif.ranger\"),\n",
" n_folds = 10, n_rep = 1, ...) {\n",
"\n",
" # warning if n_rep > 1 to handle mapping from psi to inf.func\n",
" if (n_rep > 1) {\n",
" warning(\"n_rep > 1 is not supported.\")\n",
" }\n",
" # Compute difference in outcomes\n",
" delta_y <- y1 - y0\n",
" # Prepare data backend\n",
" dml_data = DoubleML::double_ml_data_from_matrix(X = covariates, y = delta_y, d = D)\n",
" # Compute the ATT\n",
" dml_obj = DoubleML::DoubleMLIRM$new(dml_data, ml_g = ml_g, ml_m = ml_m, score = \"ATTE\", n_folds = n_folds)\n",
" dml_obj$fit(store_predictions = TRUE)\n",
" att = dml_obj$coef[1]\n",
" # Return results\n",
" inf.func <- dml_obj$psi[, 1, 1]\n",
"\n",
" # Evaluate learner performance\n",
" predictions = dml_obj$predictions\n",
" params_names = dml_obj$params_names()\n",
" eval_predictions = eval_preds(delta_y, D, predictions, params_names)\n",
" print(eval_predictions)\n",
"\n",
" output <- list(ATT = att, att.inf.func = inf.func)\n",
" return(output)\n",
" }\n"
]
},
{
"cell_type": "markdown",
"id": "b0bca913",
"metadata": {},
"source": [
"Running the evaluation wrappers helps to see that the random forest learner has a higher RMSE for predicting the outcome $E[\\Delta Y|D=1,X]$. Both models predict individuals' treatment (group) status with an accuracy of $1$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "de15e86b",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:44:05.803609Z",
"iopub.status.busy": "2024-08-31T09:44:05.802736Z",
"iopub.status.idle": "2024-08-31T09:44:09.793389Z",
"shell.execute_reply": "2024-08-31T09:44:09.792200Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"current period: 2 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.425103\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.409154\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.397313\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 2 \n",
"current group: 3 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.425493\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 3 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.40676\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 3 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.421083\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 2 \n",
"current group: 4 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.426055\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 4 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.40583\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 4 \n",
"set pretreatment period to be 3 \n",
"$ml_g0\n",
"[1] 1.423951\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"att_gt(yname = \"Y\", tname = \"period\", idname = \"id\", gname = \"G\", \n",
" xformla = ~X, data = dta, est_method = doubleml_did_eval_linear, \n",
" print_details = TRUE)\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"Group-Time Average Treatment Effects:\n",
" Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band] \n",
" 2 2 0.9145 0.0654 0.7418 1.0872 *\n",
" 2 3 1.9951 0.0655 1.8221 2.1681 *\n",
" 2 4 2.9561 0.0653 2.7838 3.1285 *\n",
" 3 2 -0.0418 0.0704 -0.2276 0.1441 \n",
" 3 3 1.1041 0.0649 0.9327 1.2754 *\n",
" 3 4 2.0533 0.0669 1.8768 2.2298 *\n",
" 4 2 -0.0028 0.0716 -0.1918 0.1862 \n",
" 4 3 0.0635 0.0675 -0.1147 0.2416 \n",
" 4 4 0.9609 0.0673 0.7833 1.1386 *\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"P-value for pre-test of parallel trends assumption: 0.64579\n",
"Control Group: Never Treated, Anticipation Periods: 0\n"
]
}
],
"source": [
"# Run estimation with evaluation: Linear model\n",
"set.seed(1234)\n",
"example_attgt_dml_eval_linear <- att_gt(yname = \"Y\",\n",
" tname = \"period\",\n",
" idname = \"id\",\n",
" gname = \"G\",\n",
" xformla = ~X,\n",
" data = dta,\n",
" est_method = doubleml_did_eval_linear,\n",
" print_details = TRUE)\n",
"\n",
"\n",
"summary(example_attgt_dml_eval_linear)\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "13fe0be0",
"metadata": {
"execution": {
"iopub.execute_input": "2024-08-31T09:44:09.797370Z",
"iopub.status.busy": "2024-08-31T09:44:09.796527Z",
"iopub.status.idle": "2024-08-31T09:44:25.338660Z",
"shell.execute_reply": "2024-08-31T09:44:25.337087Z"
},
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"current period: 2 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.571778\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.916236\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 2 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 2.404318\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 2 \n",
"current group: 3 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.570486\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 3 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.529782\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 3 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.88629\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 2 \n",
"current group: 4 \n",
"set pretreatment period to be 1 \n",
"$ml_g0\n",
"[1] 1.570562\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 3 \n",
"current group: 4 \n",
"set pretreatment period to be 2 \n",
"$ml_g0\n",
"[1] 1.529405\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n",
"current period: 4 \n",
"current group: 4 \n",
"set pretreatment period to be 3 \n",
"$ml_g0\n",
"[1] 1.560689\n",
"\n",
"$ml_m\n",
"[1] 1\n",
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Call:\n",
"att_gt(yname = \"Y\", tname = \"period\", idname = \"id\", gname = \"G\", \n",
" xformla = ~X, data = dta, est_method = doubleml_did_eval_rf, \n",
" print_details = TRUE)\n",
"\n",
"Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. \"Difference-in-Differences with Multiple Time Periods.\" Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. , \n",
"\n",
"Group-Time Average Treatment Effects:\n",
" Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band] \n",
" 2 2 0.9647 0.0736 0.7603 1.1691 *\n",
" 2 3 2.1055 0.0894 1.8571 2.3538 *\n",
" 2 4 3.1295 0.1096 2.8250 3.4339 *\n",
" 3 2 0.0820 0.0719 -0.1177 0.2818 \n",
" 3 3 1.2075 0.0707 1.0112 1.4039 *\n",
" 3 4 2.2865 0.0916 2.0323 2.5408 *\n",
" 4 2 0.1807 0.0743 -0.0257 0.3871 \n",
" 4 3 0.2451 0.0662 0.0611 0.4290 *\n",
" 4 4 1.1401 0.0711 0.9425 1.3376 *\n",
"---\n",
"Signif. codes: `*' confidence band does not cover 0\n",
"\n",
"P-value for pre-test of parallel trends assumption: 2e-05\n",
"Control Group: Never Treated, Anticipation Periods: 0\n"
]
}
],
"source": [
"# Run estimation with evaluation: Linear model\n",
"set.seed(1234)\n",
"example_attgt_dml_eval_rf <- att_gt(yname = \"Y\",\n",
" tname = \"period\",\n",
" idname = \"id\",\n",
" gname = \"G\",\n",
" xformla = ~X,\n",
" data = dta,\n",
" est_method = doubleml_did_eval_rf,\n",
" print_details = TRUE)\n",
"\n",
"\n",
"summary(example_attgt_dml_eval_rf)"
]
},
{
"cell_type": "markdown",
"id": "c5734056",
"metadata": {},
"source": [
"### Acknowledgements and Final Remarks\n",
"\n",
"We'd like to thank the authors of the `did` package for R for maintaining a flexible interface for multi-period DiD models.\n",
"\n",
"We'd like to note that the implementation presented is here is very similar to the one implemented in the Python package. For more details, we would like to reference to the [DiD](https://docs.doubleml.org/stable/examples/py_double_ml_did.html) and [pretesting](https://docs.doubleml.org/stable/examples/py_double_ml_did_pretest.html) examples in Python."
]
},
{
"cell_type": "markdown",
"id": "76593791",
"metadata": {},
"source": [
"### References\n",
"\n",
"Callaway, Brantly, and Pedro HC Sant’Anna. \"Difference-in-differences with multiple time periods.\" Journal of Econometrics 225.2 (2021): 200-230.\n",
"\n",
"Sant’Anna, Pedro HC, and Jun Zhao. \"Doubly robust difference-in-differences estimators.\" Journal of Econometrics 219.1 (2020): 101-122."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "R",
"language": "R",
"name": "ir"
},
"language_info": {
"codemirror_mode": "r",
"file_extension": ".r",
"mimetype": "text/x-r-source",
"name": "R",
"pygments_lexer": "r",
"version": "4.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 5
}