{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Python: Static Panel Models with Fixed Effects" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we illustrate how the [DoubleML](https://docs.doubleml.org/stable/index.html) package can be used to estimate treatment effects for static panel models with fixed effects in a partially linear panel regression [DoubleMLPLPR](https://docs.doubleml.org/stable/guide/models.html#partially-linear-models-plm) model. The model is based on [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import optuna\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "from sklearn.base import clone\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.preprocessing import PolynomialFeatures\n", "from sklearn.compose import ColumnTransformer\n", "from sklearn.pipeline import make_pipeline\n", "from sklearn.base import BaseEstimator, TransformerMixin\n", "from sklearn.linear_model import LassoCV\n", "from lightgbm import LGBMRegressor\n", "\n", "from doubleml.data import DoubleMLPanelData\n", "from doubleml.plm.datasets import make_plpr_CP2025\n", "from doubleml import DoubleMLPLPR\n", "\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data\n", "\n", "We will use the implemented data generating process [make_plpr_CP2025](https://docs.doubleml.org/stable/api/datasets.html#dataset-generators) to generate data similar to the simulation in [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011). For exposition, we use the simple linear `dgp_type=\"dgp1\"`, with 150 units, 10 time periods per unit, and a true treatment effect of `theta=0.5`.\n", "\n", "We set `time_type=\"int\"` such that the time variable values will be integers. It's also possible to use `\"float\"` or `\"datetime\"` time variables with [DoubleMLPLPR](https://docs.doubleml.org/stable/api/dml_models.html#doubleml-plm)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idtimeydx1x2x3x4x5x6...x21x22x23x24x25x26x27x28x29x30
011-1.2904790.9083071.710715-1.853675-1.4739071.366514-0.3220242.944020...-1.828362-3.010547-0.840202-3.0851591.169952-0.954107-3.925198-0.779510-0.4307001.004298
112-2.850646-1.316777-0.3250434.178599-1.159857-0.139527-0.230115-0.631976...-0.724172-0.421045-2.012480-2.081784-2.734123-0.879470-2.1412184.598401-4.222797-2.523024
213-4.338502-1.756120-0.8975901.505972-0.9251891.511500-2.2065610.132579...1.766109-2.252858-2.919826-1.974066-0.7738810.244633-1.7275501.6654670.562291-1.553616
314-2.7132360.9348661.9878492.596228-0.220666-0.480717-3.966273-0.911226...0.8561240.727759-0.5015791.0775042.268052-3.8214221.629055-0.220834-1.185091-5.462884
415-5.782997-4.357881-3.0865593.796975-1.539641-2.425617-1.020599-1.666200...2.617215-1.231835-0.8913500.2469812.4896420.319735-2.8103660.5858263.6437490.147147
\n", "

5 rows × 34 columns

\n", "
" ], "text/plain": [ " id time y d x1 x2 x3 x4 \\\n", "0 1 1 -1.290479 0.908307 1.710715 -1.853675 -1.473907 1.366514 \n", "1 1 2 -2.850646 -1.316777 -0.325043 4.178599 -1.159857 -0.139527 \n", "2 1 3 -4.338502 -1.756120 -0.897590 1.505972 -0.925189 1.511500 \n", "3 1 4 -2.713236 0.934866 1.987849 2.596228 -0.220666 -0.480717 \n", "4 1 5 -5.782997 -4.357881 -3.086559 3.796975 -1.539641 -2.425617 \n", "\n", " x5 x6 ... x21 x22 x23 x24 x25 \\\n", "0 -0.322024 2.944020 ... -1.828362 -3.010547 -0.840202 -3.085159 1.169952 \n", "1 -0.230115 -0.631976 ... -0.724172 -0.421045 -2.012480 -2.081784 -2.734123 \n", "2 -2.206561 0.132579 ... 1.766109 -2.252858 -2.919826 -1.974066 -0.773881 \n", "3 -3.966273 -0.911226 ... 0.856124 0.727759 -0.501579 1.077504 2.268052 \n", "4 -1.020599 -1.666200 ... 2.617215 -1.231835 -0.891350 0.246981 2.489642 \n", "\n", " x26 x27 x28 x29 x30 \n", "0 -0.954107 -3.925198 -0.779510 -0.430700 1.004298 \n", "1 -0.879470 -2.141218 4.598401 -4.222797 -2.523024 \n", "2 0.244633 -1.727550 1.665467 0.562291 -1.553616 \n", "3 -3.821422 1.629055 -0.220834 -1.185091 -5.462884 \n", "4 0.319735 -2.810366 0.585826 3.643749 0.147147 \n", "\n", "[5 rows x 34 columns]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.random.seed(123)\n", "data = make_plpr_CP2025(num_id=150, num_t=10, dim_x=30, theta=0.5, dgp_type=\"dgp1\", time_type=\"int\")\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create a corresponding [DoubleMLPanelData](https://docs.doubleml.org/stable/api/generated/doubleml.data.DoubleMLPanelData.html) object, we need to set `static_panel=True` and specify `id_col` and `time_col` columns." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "data_obj = DoubleMLPanelData(data, y_col=\"y\", d_cols=\"d\", t_col=\"time\", id_col=\"id\", static_panel=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The partially linear panel regression (PLPR) model extends the partially linear model to panel data by introducing fixed effects $\\alpha_i^*$.\n", "\n", "The PLPR model takes the form\n", "\n", "$$\n", "\\begin{align*}\n", " Y_{it} &= \\theta_0 D_{it} + g_1(X_{it}) + \\alpha_i^* + U_{it}, \\\\\n", " D_{it} &= m_1(X_{it}) + \\gamma_i + V_{it},\n", "\\end{align*}\n", "$$\n", "\n", "where\n", "- $Y_{it}$ outcome, $D_{it}$ treatment, $X_{it}$ covariates, $\\theta_0$ causal treatment effect\n", "- $g_1$ and $m_1$ nuisance functions\n", "- $\\alpha_i^*$, $\\gamma_i$ unobserved individual heterogeneity, correlated with covariates\n", "- $U_{it}$, $V_{it}$ error terms\n", "\n", "Further note $\\mathbb{E}[U_{it} \\mid D_{it}, X_{it}, \\alpha_i^*] = 0$ and $\\mathbb{E}[V_{it} \\mid X_{it}, \\gamma_i]=0$, but $\\mathbb{E}[\\alpha_i^* \\mid D_{it}, X_{it}] \\neq 0$.\n", "\n", "Alternatively we can write the partialling-out PLPR as \n", "\n", "$$\n", "\\begin{align*}\n", " Y_{it} &= \\theta_0 V_{it} + \\ell_1(X_{it}) + \\alpha_i + U_{it}, \\\\\n", " V_{it} &= D_{it} - m_1(X_{it}) - \\gamma_i,\n", "\\end{align*}\n", "$$\n", "\n", "with nuisance function $\\ell_1$ and fixed effect $\\alpha_i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assumptions\n", "\n", "Define $\\xi_i$ as time-invariant heterogeneity terms influencing outcome and treatment and $L_{t-1}(W_i) = \\{ W_{i1}, \\dots, W_{it-1} \\}$ as lags of a random variable $W_{it}$ at wave $t$.\n", "\n", "- *No feedback to predictors*\n", "$$ X_{it} \\perp L_{t-1} (Y_i, D_i) \\mid L_{t-1} (X_i), \\xi_i $$\n", "- *Static panel*\n", "$$ Y_{it}, D_{it} \\perp L_{t-1} (Y_i, X_i, D_i) \\mid X_{it}, \\xi_i $$\n", "- *Selection on observables and omitted time-invariant variables*\n", "$$ Y_{it} (.) \\perp D_{it} \\mid X_{it}, \\xi_i $$\n", "- *Homogeneity and linearity of the treatment effect*\n", "$$ \\mathbb{E} [Y_{it}(d) - Y_{it}(0) \\mid X_{it}, \\xi_i] = d \\theta_0 $$\n", "- *Additive Separability*\n", "$$\n", "\\begin{align*}\n", "\\mathbb{E} [Y_{it}(0) \\mid X_{it}, \\xi_i] &= g_1(X_{it}) + \\alpha^*_i \\quad \\text{where } \\alpha^*_i = \\alpha^*(\\xi_i), \\\\\n", "\\mathbb{E} [D_{it} \\mid X_{it}, \\xi_i] &= m_1(X_{it}) + \\gamma_i \\quad \\text{where } \\gamma_i = \\gamma(\\xi_i)\n", "\\end{align*}\n", "$$\n", "\n", "For more information, see [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To estimate the causal effect, we can create a [DoubleMLPLPR](https://docs.doubleml.org/stable/api/dml_models.html#doubleml-plm) object. The model described in [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011) uses block-k-fold cross-fitting, where the entire time series of the sampled unit is allocated to one fold to allow for possible serial correlation within each unit which is common with panel data. Furthermore, cluster robust standard error are employed. [DoubleMLPLPR](https://docs.doubleml.org/stable/guide/models.html#partially-linear-models-plm) implements both aspects by using `id_col` as the cluster variable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimation Approaches\n", "\n", "[Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011) describes multiple estimation approaches, which can be set with the `approach` parameter. Depending on the type of `approach`, different data transformations are performed along the way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correlated Random Effects\n", "\n", "The correlated random effects (cre) approaches includes the a general approach (`cre_general`) and an approach relying on normality assumptions (`cre_normal`).\n", "\n", "#### `cre_general`\n", "\n", "The `cre_general` approach:\n", "\n", "- Learning $\\ell_1$ from $\\{ Y_{it}, X_{it}, \\bar{X}_i : t=1,\\dots, T \\}_{i=1}^N$.\n", "- First learning $\\tilde{m}_1$ from $\\{ D_{it}, X_{it}, \\bar{X}_i : t=1,\\dots, T \\}_{i=1}^N$, with predictions $\\hat{m}_{1,it} = \\tilde{m}_1 (X_{it}, \\bar{X}_i)$\n", " - Calculate $\\hat{\\bar{m}}_i = T^{-1} \\sum_{t=1}^T \\hat{m}_{1,it}$,\n", " - Calculate final nuisance part as $\\hat{m}^*_1 (X_{it}, \\bar{X}_i, \\bar{D}_i) = \\hat{m}_{1,it} + \\bar{D}_i - \\hat{\\bar{m}}_i$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "learner = LassoCV()\n", "ml_l = clone(learner)\n", "ml_m = clone(learner)\n", "\n", "dml_plpr_cre_general = DoubleMLPLPR(data_obj, ml_l=ml_l, ml_m=ml_m, approach=\"cre_general\", n_folds=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can look at the the transformed data using the `data_transform` property after the `DoubleMLPLPR` object was created." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idtimeydx1x2x3x4x5x6...x21_meanx22_meanx23_meanx24_meanx25_meanx26_meanx27_meanx28_meanx29_meanx30_mean
011-1.2904790.9083071.710715-1.853675-1.4739071.366514-0.3220242.944020...1.24018-0.52821-0.7341450.2274941.1647630.412979-1.2726080.459816-0.829863-1.145189
112-2.850646-1.316777-0.3250434.178599-1.159857-0.139527-0.230115-0.631976...1.24018-0.52821-0.7341450.2274941.1647630.412979-1.2726080.459816-0.829863-1.145189
213-4.338502-1.756120-0.8975901.505972-0.9251891.511500-2.2065610.132579...1.24018-0.52821-0.7341450.2274941.1647630.412979-1.2726080.459816-0.829863-1.145189
314-2.7132360.9348661.9878492.596228-0.220666-0.480717-3.966273-0.911226...1.24018-0.52821-0.7341450.2274941.1647630.412979-1.2726080.459816-0.829863-1.145189
415-5.782997-4.357881-3.0865593.796975-1.539641-2.425617-1.020599-1.666200...1.24018-0.52821-0.7341450.2274941.1647630.412979-1.2726080.459816-0.829863-1.145189
\n", "

5 rows × 64 columns

\n", "
" ], "text/plain": [ " id time y d x1 x2 x3 x4 \\\n", "0 1 1 -1.290479 0.908307 1.710715 -1.853675 -1.473907 1.366514 \n", "1 1 2 -2.850646 -1.316777 -0.325043 4.178599 -1.159857 -0.139527 \n", "2 1 3 -4.338502 -1.756120 -0.897590 1.505972 -0.925189 1.511500 \n", "3 1 4 -2.713236 0.934866 1.987849 2.596228 -0.220666 -0.480717 \n", "4 1 5 -5.782997 -4.357881 -3.086559 3.796975 -1.539641 -2.425617 \n", "\n", " x5 x6 ... x21_mean x22_mean x23_mean x24_mean x25_mean \\\n", "0 -0.322024 2.944020 ... 1.24018 -0.52821 -0.734145 0.227494 1.164763 \n", "1 -0.230115 -0.631976 ... 1.24018 -0.52821 -0.734145 0.227494 1.164763 \n", "2 -2.206561 0.132579 ... 1.24018 -0.52821 -0.734145 0.227494 1.164763 \n", "3 -3.966273 -0.911226 ... 1.24018 -0.52821 -0.734145 0.227494 1.164763 \n", "4 -1.020599 -1.666200 ... 1.24018 -0.52821 -0.734145 0.227494 1.164763 \n", "\n", " x26_mean x27_mean x28_mean x29_mean x30_mean \n", "0 0.412979 -1.272608 0.459816 -0.829863 -1.145189 \n", "1 0.412979 -1.272608 0.459816 -0.829863 -1.145189 \n", "2 0.412979 -1.272608 0.459816 -0.829863 -1.145189 \n", "3 0.412979 -1.272608 0.459816 -0.829863 -1.145189 \n", "4 0.412979 -1.272608 0.459816 -0.829863 -1.145189 \n", "\n", "[5 rows x 64 columns]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dml_plpr_cre_general.data_transform.data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the covariates inlcude the original $X_{it}$ and additionally the unit mean $\\bar{X}_i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After fitting the model, we can print the `DoubleMLPLPR` object. The data summary corresponds to the transformed data. Additional Information at the end also includes a pre-transformation data summary." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "================== DoubleMLPLPR Object ==================\n", "\n", "------------------ Data Summary ------------------\n", "Outcome variable: y\n", "Treatment variable(s): ['d']\n", "Covariates: ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'x13', 'x14', 'x15', 'x16', 'x17', 'x18', 'x19', 'x20', 'x21', 'x22', 'x23', 'x24', 'x25', 'x26', 'x27', 'x28', 'x29', 'x30', 'x1_mean', 'x2_mean', 'x3_mean', 'x4_mean', 'x5_mean', 'x6_mean', 'x7_mean', 'x8_mean', 'x9_mean', 'x10_mean', 'x11_mean', 'x12_mean', 'x13_mean', 'x14_mean', 'x15_mean', 'x16_mean', 'x17_mean', 'x18_mean', 'x19_mean', 'x20_mean', 'x21_mean', 'x22_mean', 'x23_mean', 'x24_mean', 'x25_mean', 'x26_mean', 'x27_mean', 'x28_mean', 'x29_mean', 'x30_mean']\n", "Instrument variable(s): None\n", "Time variable: time\n", "Id variable: id\n", "Static panel data: True\n", "No. Unique Ids: 150\n", "No. Observations: 1500\n", "\n", "------------------ Score & Algorithm ------------------\n", "Score function: partialling out\n", "Static panel model approach: cre_general\n", "\n", "------------------ Machine Learner ------------------\n", "Learner ml_l: LassoCV()\n", "Learner ml_m: LassoCV()\n", "Out-of-sample Performance:\n", "Regression:\n", "Learner ml_l RMSE: [[1.8336022]]\n", "Learner ml_m RMSE: [[0.99534683]]\n", "\n", "------------------ Resampling ------------------\n", "No. folds per cluster: 5\n", "No. folds: 5\n", "No. repeated sample splits: 1\n", "\n", "------------------ Fit Summary ------------------\n", " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.490173 0.026161 18.737077 2.468215e-78 0.438899 0.541447\n", "\n", "------------------ Additional Information -------------\n", "Cluster variable(s): ['id']\n", "\n", "Pre-Transformation Data Summary: \n", "Outcome variable: y\n", "Treatment variable(s): ['d']\n", "Covariates: ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'x13', 'x14', 'x15', 'x16', 'x17', 'x18', 'x19', 'x20', 'x21', 'x22', 'x23', 'x24', 'x25', 'x26', 'x27', 'x28', 'x29', 'x30']\n", "No. Observations: 1500\n", "\n" ] } ], "source": [ "dml_plpr_cre_general.fit()\n", "print(dml_plpr_cre_general)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### `cre_normal`\n", "\n", "The `cre_normal` approach:\n", "\n", "Under the assumption that the conditional distribution $D_{i1}, \\dots, D_{iT} \\mid X_{i1}, \\dots X_{iT}$ is multivariate normal (see [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011) for further details):\n", "- Learn $\\ell_1$ from $\\{ Y_{it}, X_{it}, \\bar{X}_i : t=1,\\dots, T \\}_{i=1}^N$,\n", "- Learn $m^*_{1}$ from $\\{ D_{it}, X_{it}, \\bar{X}_i, \\bar{D}_i: t=1,\\dots, T \\}_{i=1}^N$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.505807 0.027168 18.617927 2.299468e-77 0.45256 0.559055\n" ] } ], "source": [ "dml_plpr_cre_normal = DoubleMLPLPR(data_obj, ml_l=ml_l, ml_m=ml_m, approach=\"cre_normal\", n_folds=5)\n", "dml_plpr_cre_normal.fit()\n", "print(dml_plpr_cre_normal.summary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `cre_normal` approach uses additionally inlcudes $\\bar{D}_i$ in the treatment nuisance estimation. The corresponding data can be assesses by the `d_mean` property." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.94278854],\n", " [-0.94278854],\n", " [-0.94278854],\n", " ...,\n", " [ 0.15320478],\n", " [ 0.15320478],\n", " [ 0.15320478]], shape=(1500, 1))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dml_plpr_cre_normal.d_mean" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Transformation Approaches\n", "\n", "The transformation approaches include first differences (`fd_exact`) and within-group (`wg_approx`) transformations.\n", "\n", "#### `fd_exact`\n", "\n", "Consider first differences (FD) transformation (`fd_exact`) $Q(Y_{it})= Y_{it} - Y_{it-1}$, under the assumptions from above, [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011) show that $\\mathbb{E}[Y_{it}-Y_{it-1} | X_{it-1},X_{it}] =\\Delta \\ell_1 (X_{it-1}, X_{it})$ and $\\mathbb{E}[D_{it}-D_{it-1} | X_{it-1},X_{it}] =\\Delta m_1 (X_{it-1}, X_{it})$. Therefore, the transformed nuisance function can be learnt as\n", "\n", "- $\\Delta \\ell_1 (X_{it-1}, X_{it})$ from $\\{ Y_{it}-Y_{it-1}, X_{it-1}, X_{it} : t=2, \\dots , T \\}_{i=1}^N$,\n", "- $\\Delta m_1 (X_{it-1}, X_{it})$ from $\\{ D_{it}-D_{it-1}, X_{it-1}, X_{it} : t=2, \\dots , T \\}_{i=1}^N$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idtimey_diffd_diffx1x2x3x4x5x6...x21_lagx22_lagx23_lagx24_lagx25_lagx26_lagx27_lagx28_lagx29_lagx30_lag
012-1.560167-2.225084-0.3250434.178599-1.159857-0.139527-0.230115-0.631976...-1.828362-3.010547-0.840202-3.0851591.169952-0.954107-3.925198-0.779510-0.4307001.004298
113-1.487856-0.439343-0.8975901.505972-0.9251891.511500-2.2065610.132579...-0.724172-0.421045-2.012480-2.081784-2.734123-0.879470-2.1412184.598401-4.222797-2.523024
2141.6252662.6909861.9878492.596228-0.220666-0.480717-3.966273-0.911226...1.766109-2.252858-2.919826-1.974066-0.7738810.244633-1.7275501.6654670.562291-1.553616
315-3.069761-5.292747-3.0865593.796975-1.539641-2.425617-1.020599-1.666200...0.8561240.727759-0.5015791.0775042.268052-3.8214221.629055-0.220834-1.185091-5.462884
416-1.0947990.5510510.289315-2.823134-3.137179-1.425923-0.7301160.232687...2.617215-1.231835-0.8913500.2469812.4896420.319735-2.8103660.5858263.6437490.147147
\n", "

5 rows × 64 columns

\n", "
" ], "text/plain": [ " id time y_diff d_diff x1 x2 x3 x4 \\\n", "0 1 2 -1.560167 -2.225084 -0.325043 4.178599 -1.159857 -0.139527 \n", "1 1 3 -1.487856 -0.439343 -0.897590 1.505972 -0.925189 1.511500 \n", "2 1 4 1.625266 2.690986 1.987849 2.596228 -0.220666 -0.480717 \n", "3 1 5 -3.069761 -5.292747 -3.086559 3.796975 -1.539641 -2.425617 \n", "4 1 6 -1.094799 0.551051 0.289315 -2.823134 -3.137179 -1.425923 \n", "\n", " x5 x6 ... x21_lag x22_lag x23_lag x24_lag x25_lag \\\n", "0 -0.230115 -0.631976 ... -1.828362 -3.010547 -0.840202 -3.085159 1.169952 \n", "1 -2.206561 0.132579 ... -0.724172 -0.421045 -2.012480 -2.081784 -2.734123 \n", "2 -3.966273 -0.911226 ... 1.766109 -2.252858 -2.919826 -1.974066 -0.773881 \n", "3 -1.020599 -1.666200 ... 0.856124 0.727759 -0.501579 1.077504 2.268052 \n", "4 -0.730116 0.232687 ... 2.617215 -1.231835 -0.891350 0.246981 2.489642 \n", "\n", " x26_lag x27_lag x28_lag x29_lag x30_lag \n", "0 -0.954107 -3.925198 -0.779510 -0.430700 1.004298 \n", "1 -0.879470 -2.141218 4.598401 -4.222797 -2.523024 \n", "2 0.244633 -1.727550 1.665467 0.562291 -1.553616 \n", "3 -3.821422 1.629055 -0.220834 -1.185091 -5.462884 \n", "4 0.319735 -2.810366 0.585826 3.643749 0.147147 \n", "\n", "[5 rows x 64 columns]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dml_plpr_fd_exact = DoubleMLPLPR(data_obj, ml_l=ml_l, ml_m=ml_m, approach=\"fd_exact\", n_folds=5)\n", "dml_plpr_fd_exact.data_transform.data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the outcome and treatment variables are now labeled `y_diff` and `d_diff` to indicate the first-difference transformation. Moreover, lagged covariates $X_{it-1}$ are added and rows for the first time period are dropped." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_diff 0.511822 0.032746 15.630162 4.536510e-55 0.447641 0.576002\n" ] } ], "source": [ "dml_plpr_fd_exact.fit()\n", "print(dml_plpr_fd_exact.summary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### `wg_approx`\n", "\n", "For within-group (WG) transformation (`wg_approx`) $Q(X_{it})= X_{it} - \\bar{X}_{i}$, where $\\bar{X}_{i} = T^{-1} \\sum_{t=1}^T X_{it}$, approximate the model as\n", "\n", "$$\n", "\\begin{align*}\n", " Q(Y_{it}) &\\approx \\theta_0 Q(D_{it}) + g_1 (Q(X_{it})) + Q(U_{it}), \\\\\n", " Q(D_{it}) &\\approx m_1 (Q(X_{it})) + Q(V_{it}).\n", "\\end{align*}\n", "$$\n", "\n", "Similarly for the partialling-out PLPR\n", "\n", "$$\n", "Q(Y_{it}) \\approx \\theta_0 Q(V_{it}) + \\ell_1 (Q(X_{it})) + Q(U_{it}).\n", "$$\n", "\n", "- $\\ell_1$ can be learnt from transformed data $\\{ Q(Y_{it}), Q(X_{it}) : t=1,\\dots,T \\}_{i=1}^N$,\n", "- $m_1$ can be learnt from transformed data $\\{ Q(D_{it}), Q(X_{it}) : t=1,\\dots,T \\}_{i=1}^N$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idtimey_demeand_demeanx1_demeanx2_demeanx3_demeanx4_demeanx5_demeanx6_demean...x21_demeanx22_demeanx23_demeanx24_demeanx25_demeanx26_demeanx27_demeanx28_demeanx29_demeanx30_demean
0111.5435711.7606602.207607-2.039516-0.6428471.0142040.3841661.826013...-3.162933-2.4759420.000728-3.3442190.082829-1.351303-2.670511-1.2753440.4075962.187878
112-0.016596-0.4644240.1718493.992759-0.328797-0.4918370.476074-1.749982...-2.0587430.113560-1.171550-2.340844-3.821246-1.276666-0.8865314.102567-3.384501-1.339444
213-1.504452-0.903767-0.4006981.320131-0.0941291.159190-1.500371-0.985427...0.431538-1.718253-2.078895-2.233126-1.861004-0.152563-0.4728631.1696321.400587-0.370036
3140.1208141.7872192.4847412.4103870.610394-0.833027-3.260084-2.029232...-0.4784471.2623640.3393520.8184431.180930-4.2186182.883741-0.716668-0.346795-4.279304
415-2.948947-3.505528-2.5896673.611134-0.708582-2.777927-0.314410-2.784206...1.282645-0.697230-0.050420-0.0120801.402519-0.077461-1.5556790.0899914.4820451.330727
\n", "

5 rows × 34 columns

\n", "
" ], "text/plain": [ " id time y_demean d_demean x1_demean x2_demean x3_demean x4_demean \\\n", "0 1 1 1.543571 1.760660 2.207607 -2.039516 -0.642847 1.014204 \n", "1 1 2 -0.016596 -0.464424 0.171849 3.992759 -0.328797 -0.491837 \n", "2 1 3 -1.504452 -0.903767 -0.400698 1.320131 -0.094129 1.159190 \n", "3 1 4 0.120814 1.787219 2.484741 2.410387 0.610394 -0.833027 \n", "4 1 5 -2.948947 -3.505528 -2.589667 3.611134 -0.708582 -2.777927 \n", "\n", " x5_demean x6_demean ... x21_demean x22_demean x23_demean x24_demean \\\n", "0 0.384166 1.826013 ... -3.162933 -2.475942 0.000728 -3.344219 \n", "1 0.476074 -1.749982 ... -2.058743 0.113560 -1.171550 -2.340844 \n", "2 -1.500371 -0.985427 ... 0.431538 -1.718253 -2.078895 -2.233126 \n", "3 -3.260084 -2.029232 ... -0.478447 1.262364 0.339352 0.818443 \n", "4 -0.314410 -2.784206 ... 1.282645 -0.697230 -0.050420 -0.012080 \n", "\n", " x25_demean x26_demean x27_demean x28_demean x29_demean x30_demean \n", "0 0.082829 -1.351303 -2.670511 -1.275344 0.407596 2.187878 \n", "1 -3.821246 -1.276666 -0.886531 4.102567 -3.384501 -1.339444 \n", "2 -1.861004 -0.152563 -0.472863 1.169632 1.400587 -0.370036 \n", "3 1.180930 -4.218618 2.883741 -0.716668 -0.346795 -4.279304 \n", "4 1.402519 -0.077461 -1.555679 0.089991 4.482045 1.330727 \n", "\n", "[5 rows x 34 columns]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dml_plpr_wg_approx = DoubleMLPLPR(data_obj, ml_l=ml_l, ml_m=ml_m, approach=\"wg_approx\", n_folds=5)\n", "dml_plpr_wg_approx.data_transform.data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the outcome, treatment and covariate variables are now labeled `y_deman`, `d_deman`, `xi_deman` to indicate the within-group transformations." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_demean 0.495323 0.025841 19.167824 6.872435e-82 0.444675 0.545972\n" ] } ], "source": [ "dml_plpr_wg_approx.fit()\n", "print(dml_plpr_wg_approx.summary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the simple linear data generating process `dgp_type=\"dgp1\"`, we can see that all approaches lead to estimated close the true effect of `theta=0.5`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `data_original` property additionally includes the original data before any transformation was applied." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idtimeydx1x2x3x4x5x6...x21x22x23x24x25x26x27x28x29x30
011-1.2904790.9083071.710715-1.853675-1.4739071.366514-0.3220242.944020...-1.828362-3.010547-0.840202-3.0851591.169952-0.954107-3.925198-0.779510-0.4307001.004298
112-2.850646-1.316777-0.3250434.178599-1.159857-0.139527-0.230115-0.631976...-0.724172-0.421045-2.012480-2.081784-2.734123-0.879470-2.1412184.598401-4.222797-2.523024
213-4.338502-1.756120-0.8975901.505972-0.9251891.511500-2.2065610.132579...1.766109-2.252858-2.919826-1.974066-0.7738810.244633-1.7275501.6654670.562291-1.553616
314-2.7132360.9348661.9878492.596228-0.220666-0.480717-3.966273-0.911226...0.8561240.727759-0.5015791.0775042.268052-3.8214221.629055-0.220834-1.185091-5.462884
415-5.782997-4.357881-3.0865593.796975-1.539641-2.425617-1.020599-1.666200...2.617215-1.231835-0.8913500.2469812.4896420.319735-2.8103660.5858263.6437490.147147
\n", "

5 rows × 34 columns

\n", "
" ], "text/plain": [ " id time y d x1 x2 x3 x4 \\\n", "0 1 1 -1.290479 0.908307 1.710715 -1.853675 -1.473907 1.366514 \n", "1 1 2 -2.850646 -1.316777 -0.325043 4.178599 -1.159857 -0.139527 \n", "2 1 3 -4.338502 -1.756120 -0.897590 1.505972 -0.925189 1.511500 \n", "3 1 4 -2.713236 0.934866 1.987849 2.596228 -0.220666 -0.480717 \n", "4 1 5 -5.782997 -4.357881 -3.086559 3.796975 -1.539641 -2.425617 \n", "\n", " x5 x6 ... x21 x22 x23 x24 x25 \\\n", "0 -0.322024 2.944020 ... -1.828362 -3.010547 -0.840202 -3.085159 1.169952 \n", "1 -0.230115 -0.631976 ... -0.724172 -0.421045 -2.012480 -2.081784 -2.734123 \n", "2 -2.206561 0.132579 ... 1.766109 -2.252858 -2.919826 -1.974066 -0.773881 \n", "3 -3.966273 -0.911226 ... 0.856124 0.727759 -0.501579 1.077504 2.268052 \n", "4 -1.020599 -1.666200 ... 2.617215 -1.231835 -0.891350 0.246981 2.489642 \n", "\n", " x26 x27 x28 x29 x30 \n", "0 -0.954107 -3.925198 -0.779510 -0.430700 1.004298 \n", "1 -0.879470 -2.141218 4.598401 -4.222797 -2.523024 \n", "2 0.244633 -1.727550 1.665467 0.562291 -1.553616 \n", "3 -3.821422 1.629055 -0.220834 -1.185091 -5.462884 \n", "4 0.319735 -2.810366 0.585826 3.643749 0.147147 \n", "\n", "[5 rows x 34 columns]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dml_plpr_wg_approx.data_original.data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Feature preprocessing pipelines\n", "\n", "We can incorporate preprocessing pipelines. For example, when using Lasso, we may want to include polynomial and interaction terms. Here, we create a class that allows us to include, for example, polynomials of order 3 and interactions between all variables." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "class PolyPlus(BaseEstimator, TransformerMixin):\n", " \"\"\"PolynomialFeatures(degree=k) and additional terms x_i^(k+1).\"\"\"\n", "\n", " def __init__(self, degree=2, interaction_only=False, include_bias=False):\n", " self.degree = degree\n", " self.extra_degree = degree + 1\n", " self.interaction_only = interaction_only\n", " self.include_bias = include_bias\n", " self.poly = PolynomialFeatures(degree=degree, interaction_only=interaction_only, include_bias=include_bias)\n", "\n", " def fit(self, X, y=None):\n", " self.poly.fit(X)\n", " self.n_features_in_ = X.shape[1]\n", " return self\n", "\n", " def transform(self, X):\n", " X = np.asarray(X)\n", " X_poly = self.poly.transform(X)\n", " X_extra = X ** self.extra_degree\n", " return np.hstack([X_poly, X_extra])\n", "\n", " def get_feature_names_out(self, input_features=None):\n", " input_features = np.array(\n", " input_features\n", " if input_features is not None\n", " else [f\"x{i}\" for i in range(self.n_features_in_)]\n", " )\n", " poly_names = self.poly.get_feature_names_out(input_features)\n", " extra_names = [f\"{name}^{self.extra_degree}\" for name in input_features]\n", " return np.concatenate([poly_names, extra_names])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example we use the non-linear and discontinuous `dgp_type=\"dgp3\"`, with 30 covariates and a true treatment effect `theta=0.5`." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "dim_x = 30\n", "theta = 0.5\n", "\n", "np.random.seed(123)\n", "data_dgp3 = make_plpr_CP2025(num_id=500, num_t=10, dim_x=dim_x, theta=theta, dgp_type=\"dgp3\")\n", "dml_data_dgp3 = DoubleMLPanelData(data_dgp3, y_col=\"y\", d_cols=\"d\", t_col=\"time\", id_col=\"id\", static_panel=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can apply the polynomial and intercation transformation for specific sets of covariates. For example, for the `fd_exact` approach, we can apply it to the original $X_{it}$ and lags $X_{it-1}$ seperately using `ColumnTransformer`.\n", "\n", "To achieve this, we pass need to pass the corresponding indices for these two sets. `DoubleMLPLPR` stacks sets $X_{it}$ and $X_{it-1}$ column-wise. Given our example data has 30 covariates, this means that the first 30 features in the nuisance estimation correspond to the original $X_{it}$, and the last 30 correspond to lags $X_{it-1}$. Therefore we define the indices `indices_x` and `indices_x_tr` as below." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "indices_x = [x for x in range(dim_x)]\n", "indices_x_tr = [x + dim_x for x in indices_x]\n", "\n", "preprocessor = ColumnTransformer(\n", " [\n", " (\n", " \"poly_x\",\n", " PolyPlus(degree=2, include_bias=False, interaction_only=False),\n", " indices_x,\n", " ),\n", " (\n", " \"poly_x_tr\",\n", " PolyPlus(degree=2, include_bias=False, interaction_only=False),\n", " indices_x_tr,\n", " ),\n", " ],\n", " remainder=\"passthrough\",\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This preprocessor can be applied for approaches `cre_general` and `cre_normal` in the same fashion. In this case the two sets of covariates would be the original $X_{it}$ and the unit mean $\\bar{X}_i$.\n", "\n", "**Remark**: Note that we set `remainder=\"passthrough\"` such that all remaining features, not part of `indices_x` and `indices_x_tr`, would not be preprocessed but still included in the nuisance estimation. This is particularly important for the `cre_normal` approach, as $\\bar{D}_i$ is further added to $X_{it}$ and $\\bar{X}_i$ in the treatment nuisance model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can create the learner using a pipeline and fit the model." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Pipeline(steps=[('columntransformer',\n",
       "                 ColumnTransformer(remainder='passthrough',\n",
       "                                   transformers=[('poly_x', PolyPlus(),\n",
       "                                                  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9,\n",
       "                                                   10, 11, 12, 13, 14, 15, 16,\n",
       "                                                   17, 18, 19, 20, 21, 22, 23,\n",
       "                                                   24, 25, 26, 27, 28, 29]),\n",
       "                                                 ('poly_x_tr', PolyPlus(),\n",
       "                                                  [30, 31, 32, 33, 34, 35, 36,\n",
       "                                                   37, 38, 39, 40, 41, 42, 43,\n",
       "                                                   44, 45, 46, 47, 48, 49, 50,\n",
       "                                                   51, 52, 53, 54, 55, 56, 57,\n",
       "                                                   58, 59])])),\n",
       "                ('standardscaler', StandardScaler()),\n",
       "                ('lassocv', LassoCV(cv=2, n_alphas=20, n_jobs=5))])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
" ], "text/plain": [ "Pipeline(steps=[('columntransformer',\n", " ColumnTransformer(remainder='passthrough',\n", " transformers=[('poly_x', PolyPlus(),\n", " [0, 1, 2, 3, 4, 5, 6, 7, 8, 9,\n", " 10, 11, 12, 13, 14, 15, 16,\n", " 17, 18, 19, 20, 21, 22, 23,\n", " 24, 25, 26, 27, 28, 29]),\n", " ('poly_x_tr', PolyPlus(),\n", " [30, 31, 32, 33, 34, 35, 36,\n", " 37, 38, 39, 40, 41, 42, 43,\n", " 44, 45, 46, 47, 48, 49, 50,\n", " 51, 52, 53, 54, 55, 56, 57,\n", " 58, 59])])),\n", " ('standardscaler', StandardScaler()),\n", " ('lassocv', LassoCV(cv=2, n_alphas=20, n_jobs=5))])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ml_lasso = make_pipeline(\n", " preprocessor, StandardScaler(), LassoCV(n_alphas=20, cv=2, n_jobs=5)\n", ")\n", "\n", "ml_lasso" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_diff 0.516443 0.018066 28.585995 1.003247e-179 0.481034 0.551852\n" ] } ], "source": [ "plpr_lasso_fd = DoubleMLPLPR(dml_data_dgp3, clone(ml_lasso), clone(ml_lasso), approach=\"fd_exact\", n_folds=5)\n", "plpr_lasso_fd.fit(store_models=True)\n", "print(plpr_lasso_fd.summary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given that we apply the polynomial and interactions preprossing to two sets of 30 columns each, the number of features is 1050." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1050" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_lasso_fd.models[\"ml_m\"][\"d_diff\"][0][0].named_steps[\"lassocv\"].n_features_in_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As describes above, for the `cre_normal` approach adds $\\bar{X}_i$ to the features used in the treatment nuisance estimation." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.569761 0.028406 20.057831 1.724430e-89 0.514086 0.625436\n" ] }, { "data": { "text/plain": [ "1051" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_lasso_cre_normal = DoubleMLPLPR(dml_data_dgp3, clone(ml_lasso), clone(ml_lasso), approach=\"cre_normal\", n_folds=5)\n", "plpr_lasso_cre_normal.fit(store_models=True)\n", "print(plpr_lasso_cre_normal.summary)\n", "plpr_lasso_cre_normal.models[\"ml_m\"][\"d\"][0][0].named_steps[\"lassocv\"].n_features_in_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the `wg_approx` approach, there is only one set of features. We can create a similar learner for this setting." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "preprocessor_wg = ColumnTransformer(\n", " [\n", " (\n", " \"poly_x\",\n", " PolyPlus(degree=2, include_bias=False, interaction_only=False),\n", " indices_x,\n", " )\n", " ],\n", " remainder=\"passthrough\",\n", ")\n", "\n", "ml_lasso_wg = make_pipeline(\n", " preprocessor_wg, StandardScaler(), LassoCV(n_alphas=20, cv=2, n_jobs=5)\n", ")" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_demean 1.152402 0.014064 81.940739 0.0 1.124837 1.179966\n" ] }, { "data": { "text/plain": [ "525" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_lasso_wg = DoubleMLPLPR(dml_data_dgp3, clone(ml_lasso_wg), clone(ml_lasso_wg), approach=\"wg_approx\", n_folds=5)\n", "plpr_lasso_wg.fit(store_models=True)\n", "print(plpr_lasso_wg.summary)\n", "plpr_lasso_wg.models[\"ml_l\"][\"d_demean\"][0][0].named_steps[\"lassocv\"].n_features_in_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that for the more complicated data generating process `dgp3`, the approximation approach performs worse compared to the other approaches." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As another example, below we should how to select a specific covariate subset for preprocessing. This can be useful in case of the data includes dummy covariates, where adding polynomials might not be appropriate." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "39" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_cols = dml_data_dgp3.x_cols \n", "x_cols_to_pre = [\"x3\", \"x6\", \"x22\"]\n", "\n", "indices_x_pre = [i for i, c in enumerate(x_cols) if c in x_cols_to_pre]\n", "\n", "preprocessor_alt = ColumnTransformer(\n", " [\n", " (\n", " \"poly_x\",\n", " PolyPlus(degree=2, include_bias=False, interaction_only=False),\n", " indices_x_pre,\n", " )\n", " ],\n", " remainder=\"passthrough\",\n", ")\n", "ml_lasso_alt = make_pipeline(\n", " preprocessor_alt, StandardScaler(), LassoCV(n_alphas=20, cv=2, n_jobs=5)\n", ")\n", "\n", "plpr_lasso_wg.learner[\"ml_l\"] = ml_lasso_alt\n", "plpr_lasso_wg.learner[\"ml_m\"] = ml_lasso_alt\n", "\n", "plpr_lasso_wg.fit(store_models=True)\n", "plpr_lasso_wg.models[\"ml_l\"][\"d_demean\"][0][0].named_steps[\"lassocv\"].n_features_in_" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
ColumnTransformer(remainder='passthrough',\n",
       "                  transformers=[('poly_x', PolyPlus(), [2, 5, 21])])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
" ], "text/plain": [ "ColumnTransformer(remainder='passthrough',\n", " transformers=[('poly_x', PolyPlus(), [2, 5, 21])])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_lasso_wg.models[\"ml_l\"][\"d_demean\"][0][0].named_steps['columntransformer']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also look at the resulting features.\n", "\n", "**Remark**: Note, however, that the feature names here refer only to the corresponding `x_cols` indices, not the column names from the `pd.DataFrame` because [DoubleML](https://docs.doubleml.org/stable/index.html) uses `np.array`'s for fitting the model. Therefore the difference to the names from `x_cols_to_pre`." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['poly_x__x2', 'poly_x__x5', 'poly_x__x21', 'poly_x__x2^2',\n", " 'poly_x__x2 x5', 'poly_x__x2 x21', 'poly_x__x5^2',\n", " 'poly_x__x5 x21', 'poly_x__x21^2', 'poly_x__x2^3', 'poly_x__x5^3',\n", " 'poly_x__x21^3', 'remainder__x0', 'remainder__x1', 'remainder__x3',\n", " 'remainder__x4', 'remainder__x6', 'remainder__x7', 'remainder__x8',\n", " 'remainder__x9', 'remainder__x10', 'remainder__x11',\n", " 'remainder__x12', 'remainder__x13', 'remainder__x14',\n", " 'remainder__x15', 'remainder__x16', 'remainder__x17',\n", " 'remainder__x18', 'remainder__x19', 'remainder__x20',\n", " 'remainder__x22', 'remainder__x23', 'remainder__x24',\n", " 'remainder__x25', 'remainder__x26', 'remainder__x27',\n", " 'remainder__x28', 'remainder__x29'], dtype=object)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_lasso_wg.models[\"ml_l\"][\"d_demean\"][0][0].named_steps['columntransformer'].get_feature_names_out()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hyperparameter tuning\n", "\n", "In this section we will use the `tune_ml_models()` method to tune hyperparameters using the [Optuna](https://optuna.org/) package. More details can found in the [Python: Hyperparametertuning with Optuna](https://docs.doubleml.org/stable/examples/learners/py_optuna.html) example notebook.\n", "\n", "As an example, we use [LightGBM](https://lightgbm.readthedocs.io/en/stable/) regressors and compare the estimates for the different static panel model approaches, when applied to the non-linear and discontinuous `dgp3`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "dim_x = 30\n", "theta = 0.5\n", "\n", "np.random.seed(11)\n", "data_tune = make_plpr_CP2025(num_id=4000, num_t=10, dim_x=dim_x, theta=theta, dgp_type=\"dgp3\")\n", "dml_data_tune = DoubleMLPanelData(data_tune, y_col=\"y\", d_cols=\"d\", t_col=\"time\", id_col=\"id\", static_panel=True)\n", "ml_boost = LGBMRegressor(random_state=314, verbose=-1)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "# parameter space for both ml models\n", "def ml_params(trial):\n", " return {\n", " \"n_estimators\": 100,\n", " \"learning_rate\": trial.suggest_float(\"learning_rate\", 0.1, 0.4, log=True),\n", " \"max_depth\": trial.suggest_int(\"max_depth\", 2, 10),\n", " \"min_child_samples\": trial.suggest_int(\"min_child_samples\", 1, 5),\n", " \"reg_lambda\": trial.suggest_float(\"reg_lambda\", 1e-2, 5, log=True),\n", " }\n", "\n", "param_space = {\n", " \"ml_l\": ml_params,\n", " \"ml_m\": ml_params\n", "}\n", "\n", "optuna_settings = {\n", " \"n_trials\": 100,\n", " \"show_progress_bar\": True,\n", " \"verbosity\": optuna.logging.WARNING, # Suppress Optuna logs\n", "}" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "eb8bee93245c430fb5fce83c2c66eaf2", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/100 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errtP>|t|2.5 %97.5 %
d0.502560.00802562.6205170.00.486830.518289
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.50256 0.008025 62.620517 0.0 0.48683 0.518289" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_tune_cre_general = DoubleMLPLPR(dml_data_tune, clone(ml_boost), clone(ml_boost), approach=\"cre_general\", n_folds=5)\n", "\n", "plpr_tune_cre_general.tune_ml_models(\n", " ml_param_space=param_space,\n", " optuna_settings=optuna_settings,\n", ")\n", "\n", "plpr_tune_cre_general.fit()\n", "plpr_tune_cre_general.summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "0.509102" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "508dec2d1acd4dfe81fc4d2758acaaef", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/100 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errtP>|t|2.5 %97.5 %
d0.4617590.01030844.7968540.00.4415560.481962
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.461759 0.010308 44.796854 0.0 0.441556 0.481962" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_tune_cre_normal = DoubleMLPLPR(dml_data_tune, clone(ml_boost), clone(ml_boost), approach=\"cre_normal\", n_folds=5)\n", "\n", "plpr_tune_cre_normal.tune_ml_models(\n", " ml_param_space=param_space,\n", " optuna_settings=optuna_settings,\n", ")\n", "\n", "plpr_tune_cre_normal.fit()\n", "plpr_tune_cre_normal.summary" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1acbbcfc44ee4b9b8c995f6e6620f788", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/100 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errtP>|t|2.5 %97.5 %
d_diff0.5456980.00892861.123730.00.52820.563196
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_diff 0.545698 0.008928 61.12373 0.0 0.5282 0.563196" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_tune_fd = DoubleMLPLPR(dml_data_tune, clone(ml_boost), clone(ml_boost), approach=\"fd_exact\", n_folds=5)\n", "\n", "plpr_tune_fd.tune_ml_models(\n", " ml_param_space=param_space,\n", " optuna_settings=optuna_settings,\n", ")\n", "\n", "plpr_tune_fd.fit()\n", "plpr_tune_fd.summary" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "73577cb5b19b428dab75d207746df8db", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/100 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errtP>|t|2.5 %97.5 %
d_demean1.13650.005177219.5478460.01.1263541.146646
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_demean 1.1365 0.005177 219.547846 0.0 1.126354 1.146646" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plpr_tune_wg = DoubleMLPLPR(dml_data_tune, clone(ml_boost), clone(ml_boost), approach=\"wg_approx\", n_folds=5)\n", "\n", "plpr_tune_wg.tune_ml_models(\n", " ml_param_space=param_space,\n", " optuna_settings=optuna_settings,\n", ")\n", "\n", "plpr_tune_wg.fit()\n", "plpr_tune_wg.summary" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True treatment effect: 0.5\n", "\n", " Model theta se ci_lower ci_upper\n", "cre_general 0.502560 0.008025 0.486830 0.518289\n", " cre_normal 0.461759 0.010308 0.441556 0.481962\n", " fd_exact 0.545698 0.008928 0.528200 0.563196\n", " wg_approx 1.136500 0.005177 1.126354 1.146646\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAJOCAYAAABm7rQwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAovhJREFUeJzs3XmcjXX/x/H3ObOc2cyMGTPDTJgxJLINSraoaCy56U4kZSlKkRgqKoQ7WqxFkZBKIXXXHUUiFSlLlsiSvds2Y5ndrOf6/eE353bMMGe2M2a8no/HPB5zvuf7va7PdZZvec91fS+TYRiGAAAAAAAAACcyl3YBAAAAAAAAuPEQSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTEUoBAFBGmEwmvfLKK6VdBq4Dbdu2Vb169fLtd/ToUZlMJn3wwQclXxRuOB988IFMJpO2bt1a2qUAAMooQikAQJlx6NAhPfnkk6pRo4Y8PDzk6+urli1baubMmbp48WJpl4cbQE7Ik/Pj5uamSpUqqUWLFnrxxRd1/Pjx0i6x2Kxfv952nB9//HGefVq2bCmTyZQrIAsPD9d99913ze3369fP7rX09fVVw4YNNXXqVKWnp9v6vfLKKzKZTDp79qxDtea8LzVq1FCfPn10+PBhW78r3z+z2ayAgAB17NhRmzZtcuRlAQAAxci1tAsAAMARK1eu1IMPPiiLxaI+ffqoXr16ysjI0IYNG/Tcc89pz549eu+990q7zBJ18eJFubryn+7rQa9evdSpUydZrVZduHBBW7Zs0YwZMzRz5kzNnz9fDz30UGmXWGw8PDz0ySef6JFHHrFrP3r0qH755Rd5eHgUetsWi0Xvv/++JCk+Pl6ff/65Ro4cqS1btmjJkiUF3t7QoUN12223KTMzU7///rvee+89rVy5Un/88YdCQ0Nt/XLev+zsbB04cEDvvPOO7rrrLm3ZskX169cv9PEAAICC4f9sAQDXvSNHjuihhx5S9erVtW7dOlWpUsX23ODBg3Xw4EGtXLmyFCssOVarVRkZGfLw8CjSP/7LupSUFHl7e5d2GTaNGzfOFdIcO3ZM9957r/r27as6deqoYcOGpVRd8erUqZP+85//6OzZs6pUqZKt/ZNPPlFISIhq1aqlCxcuFGrbrq6udq/j008/rWbNmmnp0qWaNm2aXZDkiNatW6t79+6SpP79++vmm2/W0KFDtWjRIo0ePdrW78r3r3Xr1urYsaPeffddvfPOO4U6ltKSlpYmd3d3mc1cAAEAKHv4rxcA4Lr3xhtvKDk5WfPnz7cLpHLUrFlTzz77rO1xVlaWJk6cqMjISFksFoWHh+vFF1+0uyRI+t8lRuvXr1fTpk3l6emp+vXra/369ZKkL774QvXr15eHh4eaNGmi7du3243v16+ffHx8dPjwYUVHR8vb21uhoaGaMGGCDMOw6ztlyhS1aNFCgYGB8vT0VJMmTbR8+fJcx2IymTRkyBAtXrxYt956qywWi1atWmV77vI1pZKSkjRs2DCFh4fLYrEoODhY7du31++//263zc8++0xNmjSRp6enKlWqpEceeUQnTpzI81hOnDihbt26ycfHR0FBQRo5cqSys7Ov8s78z1dffaXOnTsrNDRUFotFkZGRmjhxYp5jf/vtN3Xq1EkVK1aUt7e3GjRooJkzZ+aq5dChQ+rUqZMqVKig3r17S7oUTo0YMUJVq1aVxWJR7dq1NWXKlFyv95o1a9SqVSv5+/vLx8dHtWvX1osvvmjX5+2339att94qLy8vVaxYUU2bNtUnn3yS77FeTfXq1fXBBx8oIyNDb7zxht1zhw8f1oMPPqiAgAB5eXnpjjvuyBWk5qzPc/ToUbv2nEvTcj6Xl9u2bZtatGghT09PRUREaM6cOQ7Vum/fPnXv3l0BAQHy8PBQ06ZN9Z///CfPvl27dpXFYtFnn31m1/7JJ5+oR48ecnFxcWifjjCbzWrbtq0k5XodCuPuu++WdCnYvpbWrVtLunSJcH7Onz+vkSNHqn79+vLx8ZGvr686duyonTt32vXLed+WLl2qF198UZUrV5a3t7f+8Y9/6O+//7brm7NGWH7vZ842lyxZopdffllhYWHy8vJSYmKiJMe+67t27VK/fv1sl0FXrlxZjz32mM6dO5frWE+cOKHHH3/c9r2OiIjQU089pYyMDLt+6enpiomJUVBQkLy9vXX//fcrLi4u1/a+/fZbtW7dWt7e3qpQoYI6d+6sPXv22PU5ffq0+vfvr5tuukkWi0VVqlRR165di+XzAAC4/nCmFADguvf111+rRo0aatGihUP9BwwYoEWLFql79+4aMWKEfvvtN02ePFl79+7Vv//9b7u+Bw8e1MMPP6wnn3xSjzzyiKZMmaIuXbpozpw5evHFF/X0009LkiZPnqwePXpo//79dmckZGdnq0OHDrrjjjv0xhtvaNWqVRo3bpyysrI0YcIEW7+ZM2fqH//4h3r37q2MjAwtWbJEDz74oFasWKHOnTvb1bRu3TotW7ZMQ4YMUaVKlRQeHp7ncQ4aNEjLly/XkCFDVLduXZ07d04bNmzQ3r171bhxY0mXgo7+/fvrtttu0+TJk3XmzBnNnDlTGzdu1Pbt2+Xv7293LNHR0WrWrJmmTJmi77//XlOnTlVkZKSeeuqpa77mH3zwgXx8fBQTEyMfHx+tW7dOY8eOVWJiot58801bvzVr1ui+++5TlSpV9Oyzz6py5crau3evVqxYkStYjI6OVqtWrTRlyhR5eXnJMAz94x//0A8//KDHH39cjRo10urVq/Xcc8/pxIkTmj59uiRpz549uu+++9SgQQNNmDBBFotFBw8e1MaNG23bnzdvnoYOHaru3bvr2WefVVpamnbt2qXffvtNDz/88DWP9VqaN2+uyMhIrVmzxtZ25swZtWjRQqmpqRo6dKgCAwO1aNEi/eMf/9Dy5ct1//33F2pfFy5cUKdOndSjRw/16tVLy5Yt01NPPSV3d3c99thjVx23Z88etWzZUmFhYRo1apS8vb21bNkydevWTZ9//nmuery8vNS1a1d9+umnts/Bzp07tWfPHr3//vvatWtXoeq/mpxgKDAw0Gnbygk8KlasmO82Dx8+rC+//FIPPvigIiIidObMGc2dO1dt2rTRn3/+mevsrldffVUmk0kvvPCCYmNjNWPGDLVr1047duyQp6enrV9B3s+JEyfK3d1dI0eOVHp6utzd3R3+rq9Zs0aHDx9W//79VblyZdulz3v27NGvv/4qk8kkSTp58qRuv/12xcfH64knntAtt9yiEydOaPny5UpNTZW7u7utnmeeeUYVK1bUuHHjdPToUc2YMUNDhgzR0qVLbX0++ugj9e3bV9HR0Xr99deVmpqqd999V61atdL27dtt89wDDzygPXv26JlnnlF4eLhiY2O1Zs0aHT9+/KpzIQCgDDMAALiOJSQkGJKMrl27OtR/x44dhiRjwIABdu0jR440JBnr1q2ztVWvXt2QZPzyyy+2ttWrVxuSDE9PT+PYsWO29rlz5xqSjB9++MHW1rdvX0OS8cwzz9jarFar0blzZ8Pd3d2Ii4uztaemptrVk5GRYdSrV8+4++677dolGWaz2dizZ0+uY5NkjBs3zvbYz8/PGDx48FVfi4yMDCM4ONioV6+ecfHiRVv7ihUrDEnG2LFjcx3LhAkT7LYRFRVlNGnS5Kr7uNrxGYZhPPnkk4aXl5eRlpZmGIZhZGVlGREREUb16tWNCxcu2PW1Wq25ahk1apRdny+//NKQZPzrX/+ya+/evbthMpmMgwcPGoZhGNOnTzck2b3+V+ratatx66235ntcVzpy5IghyXjzzTevuW1JRkJCgmEYhjFs2DBDkvHzzz/b+iQlJRkRERFGeHi4kZ2dbRiGYSxcuNCQZBw5csRuez/88EOuz16bNm0MScbUqVNtbenp6UajRo2M4OBgIyMjw67ehQsX2vrdc889Rv369W3vi2Fcev1btGhh1KpVK9d+P/vsM2PFihWGyWQyjh8/bhiGYTz33HNGjRo1bLVc+VpWr17d6Ny581VfI8O49D57e3sbcXFxRlxcnHHw4EFj0qRJhslkMho0aGDrN27cuHzfz5xaFyxYYMTFxRknT540Vq5caYSHhxsmk8nYsmWL3esxfvx4Iy4uzjh9+rTx888/G7fddpvtWPOTlpZme89yHDlyxLBYLHbfn5yawsLCjMTERFv7smXLDEnGzJkzbW2Ovp8526xRo4bdd64g3/W8vquffvqpIcn46aefbG19+vQxzGaz7bW7XM73Necz265dO7vv8PDhww0XFxcjPj7eMIxLn3d/f39j4MCBdts5ffq04efnZ2u/cOFCvt8vAED5wuV7AIDrWs5lKRUqVHCo/zfffCNJiomJsWsfMWKEJOW6ZKpu3bpq3ry57XGzZs0kXbrsp1q1arnaL7+TV44hQ4bYfs+5/C4jI0Pff/+9rf3KMyISEhLUunXrXJfaSVKbNm1Ut27dfI5U8vf312+//aaTJ0/m+fzWrVsVGxurp59+2m49qs6dO+uWW27Jcx2uQYMG2T1u3bp1nsd8pcuPLykpSWfPnlXr1q2Vmpqqffv2SZK2b9+uI0eOaNiwYXZnaEmynZ1xuSvPzvrmm2/k4uKioUOH2rWPGDFChmHo22+/lSTbtr/66itZrdY86/X399d///tfbdmyJd9jKygfHx9Jl16HnLpvv/12tWrVyq7PE088oaNHj+rPP/8s1H5cXV315JNP2h67u7vrySefVGxsrLZt25bnmPPnz2vdunXq0aOH7X06e/aszp07p+joaP3111+5LveSpHvvvVcBAQFasmSJDMPQkiVL1KtXr0LVfbmUlBQFBQUpKChINWvW1IsvvqjmzZvnOqPRUY899piCgoIUGhqqzp07KyUlRYsWLVLTpk3t+o0bN05BQUGqXLmyWrdurb1792rq1Km29aiuxWKx2M6WzM7O1rlz52yXiOb1fe7Tp4/d/NW9e3dVqVLFNlflKMj72bdvX7vvXEG+65ePS0tL09mzZ3XHHXdIkq1+q9WqL7/8Ul26dMn12km5v69PPPGEXVvr1q2VnZ2tY8eOSbp0dlZ8fLx69epl+8ydPXtWLi4uatasmX744Qdbbe7u7lq/fn2h1ykDAJQthFIAgOuar6+vpP/9Az8/x44dk9lsVs2aNe3aK1euLH9/f9s/knJcHjxJkp+fnySpatWqebZf+Q8ls9msGjVq2LXdfPPNkuzXxFmxYoXuuOMOeXh4KCAgQEFBQXr33XeVkJCQ6xgiIiLyO0xJl9ba2r17t6pWrarbb79dr7zyil2AlHOstWvXzjX2lltuyfVaeHh4KCgoyK6tYsWKDv3jcM+ePbr//vvl5+cnX19fBQUF2RaSzjnGnEup6tWrl+/2XF1dddNNN9m1HTt2TKGhobkCyjp16tiel6SePXuqZcuWGjBggEJCQvTQQw9p2bJldgHVCy+8IB8fH91+++2qVauWBg8ebHd5X1EkJydL+l+QeuzYsTzfgyvrLqjQ0NBci7/n9dm73MGDB2UYhsaMGWMLg3J+xo0bJ0mKjY3NNc7NzU0PPvigPvnkE/3000/6+++/i3SZYw4PDw+tWbNGa9assW1348aNub5Tjho7dqzWrFmjdevWadeuXTp58qQeffTRXP2eeOIJrVmzRl9//bWGDx+uixcvOrR2mnQpsJk+fbpq1aoli8WiSpUqKSgoSLt27crz+1yrVi27xyaTSTVr1sz1HhXk/bxyjijId/38+fN69tlnFRISIk9PTwUFBdm2l1N/XFycEhMTHfquSrnn0ZzLIHPmjr/++kvSpbD/ys/dd999Z/vMWSwWvf766/r2228VEhKiO++8U2+88YZOnz7tUB0AgLKHNaUAANc1X19fhYaGavfu3QUal9eZN3m52iLNV2s3rlhQ2xE///yz/vGPf+jOO+/UO++8oypVqsjNzU0LFy7Mc2Hty89kuJYePXqodevW+ve//63vvvtOb775pl5//XV98cUX6tixY4HrLOyC1fHx8WrTpo18fX01YcIERUZGysPDQ7///rteeOGFq56tdC2Xn41SUJ6envrpp5/0ww8/aOXKlVq1apWWLl2qu+++W999951cXFxUp04d7d+/XytWrNCqVav0+eef65133tHYsWM1fvz4Qu03x+7duxUcHGwLVB11tc+so2GJI3Lei5EjRyo6OjrPPlcGujkefvhhzZkzR6+88ooaNmzo0Nl8+XFxcVG7du2KvJ0c9evXd2h7tWrVsvW777775OLiolGjRumuu+7K88ygy02aNEljxozRY489pokTJyogIEBms1nDhg0r1Ge9MBydI/LSo0cP/fLLL3ruuefUqFEj+fj4yGq1qkOHDoWuP7/5Mme7H330kSpXrpyrn6vr//5JMmzYMHXp0kVffvmlVq9erTFjxmjy5Mlat26doqKiClUfAOD6RSgFALju3XfffXrvvfe0adMmu0vt8lK9enVZrVb99ddftjNRpEuLTcfHx6t69erFWpvVatXhw4dtZzRI0oEDByTJtijv559/Lg8PD61evVoWi8XWb+HChUXef5UqVfT000/r6aefVmxsrBo3bqxXX31VHTt2tB3r/v37bXchy7F///5iey3Wr1+vc+fO6YsvvtCdd95pa7/yjmeRkZGSLoU2hQkiqlevru+//15JSUl2Z0vlXB54+fGYzWbdc889uueeezRt2jRNmjRJL730kn744Qfbvr29vdWzZ0/17NlTGRkZ+uc//6lXX31Vo0ePtrsEqiA2bdqkQ4cO2c4Sy6lr//79ufpeWXfO2SXx8fF2/a52JtXJkyeVkpJid3bNlZ+9K+WcgeTm5lbg96BVq1aqVq2a1q9fr9dff71AY693L730kubNm6eXX37ZdrfLq1m+fLnuuusuzZ8/3649Pj5elSpVytU/5yyhHIZh6ODBg2rQoIFde2HezxyOftcvXLigtWvXavz48Ro7duxVawwKCpKvr2+B/xhwNTnf/eDgYIc+d5GRkRoxYoRGjBihv/76S40aNdLUqVP18ccfF0s9AIDrB5fvAQCue88//7y8vb01YMAAnTlzJtfzhw4d0syZMyVJnTp1kiTNmDHDrs+0adMkKded7orDrFmzbL8bhqFZs2bJzc1N99xzj6RLZxGYTCa7M16OHj2qL7/8stD7zM7OznWpUHBwsEJDQ5Weni5Jatq0qYKDgzVnzhxbm3Tptux79+4tttci5yyJy88iy8jI0DvvvGPXr3HjxoqIiNCMGTNyBS+OnIHWqVMnZWdn273ekjR9+nSZTCbb2WHnz5/PNbZRo0aSZHsdzp07Z/e8u7u76tatK8MwlJmZmW8teTl27Jj69esnd3d3Pffcc3Z1b968WZs2bbK1paSk6L333lN4eLjtjKOcf7j/9NNPtn7Z2dl677338txfVlaW5s6da3uckZGhuXPnKigoSE2aNMlzTHBwsNq2bau5c+fq1KlTuZ6Pi4u76vGZTCa99dZbGjduXJ6XxJVl/v7+evLJJ7V69Wrt2LHjmn1dXFxyfV4/++yzPNfikqQPP/zQ7vLj5cuX69SpU7nOZizM+5nD0e96Xt9VKfd8aTab1a1bN3399dfaunVrrv0V9IzR6Oho+fr6atKkSXl+v3I+d6mpqUpLS7N7LjIyUhUqVLA7LgBA+cGZUgCA615kZKQ++eQT9ezZU3Xq1FGfPn1Ur149ZWRk6JdfftFnn32mfv36SZIaNmyovn376r333rNdVrZ582YtWrRI3bp101133VWstXl4eGjVqlXq27evmjVrpm+//VYrV67Uiy++aFufqXPnzpo2bZo6dOighx9+WLGxsZo9e7Zq1qypXbt2FWq/SUlJuummm9S9e3c1bNhQPj4++v7777VlyxZNnTpV0qWzYV5//XX1799fbdq0Ua9evWy3iQ8PD9fw4cOL5TVo0aKFKlasqL59+2ro0KEymUz66KOPcv3D1Ww2691331WXLl3UqFEj9e/fX1WqVNG+ffu0Z88erV69+pr76dKli+666y699NJLOnr0qBo2bKjvvvtOX331lYYNG2YLdSZMmKCffvpJnTt3VvXq1RUbG6t33nlHN910k22x8XvvvVeVK1dWy5YtFRISor1792rWrFnq3LmzQ4vq//777/r4449ltVoVHx+vLVu26PPPP7cd++VnwYwaNUqffvqpOnbsqKFDhyogIECLFi3SkSNH9Pnnn9suU7z11lt1xx13aPTo0Tp//rxtYfGsrKw8awgNDdXrr7+uo0eP6uabb9bSpUu1Y8cOvffee3Jzc7tq7bNnz1arVq1Uv359DRw4UDVq1NCZM2e0adMm/fe//9XOnTuvOrZr167q2rVrvq+PdGn9qn/961+52qOiogoViE6bNk1eXl52bWazWS+++GKBt5WXZ599VjNmzNBrr72mJUuWXLXffffdpwkTJqh///5q0aKF/vjjDy1evPiq62AFBASoVatW6t+/v86cOaMZM2aoZs2aGjhwoF2/wr6fkuPfdV9fX9s6TZmZmQoLC9N3332X66xG6dJlit99953atGmjJ554QnXq1NGpU6f02WefacOGDbluVnAtvr6+evfdd/Xoo4+qcePGeuihhxQUFKTjx49r5cqVatmypWbNmqUDBw7onnvuUY8ePVS3bl25urrq3//+t86cOaOHHnrI4f0BAMqQ0rjlHwAAhXHgwAFj4MCBRnh4uOHu7m5UqFDBaNmypfH222/b3d4+MzPTGD9+vBEREWG4ubkZVatWNUaPHm3XxzCuftt6ScbgwYPt2nJuJX/5rcpzbml/6NAh49577zW8vLyMkJAQY9y4cbluGT9//nyjVq1ahsViMW655RZj4cKFtlvd57fvy58bN26cYRiXbhf/3HPPGQ0bNjQqVKhgeHt7Gw0bNjTeeeedXOOWLl1qREVFGRaLxQgICDB69+5t/Pe//7Xrk3MsV8qrxrxs3LjRuOOOOwxPT08jNDTUeP75543Vq1cbkowffvjBru+GDRuM9u3b2+pu0KCB8fbbb+dbi2FcurX88OHDjdDQUMPNzc2oVauW8eabb9rdjn7t2rVG165djdDQUMPd3d0IDQ01evXqZRw4cMDWZ+7cucadd95pBAYGGhaLxYiMjDSee+45IyEh4ZrHmfM5yPlxdXU1AgICjGbNmhmjR482jh07lue4Q4cOGd27dzf8/f0NDw8P4/bbbzdWrFiRZ7927doZFovFCAkJMV588UVjzZo1uV7HNm3aGLfeequxdetWo3nz5oaHh4dRvXp1Y9asWXnWu3Dhwlz76dOnj1G5cmXDzc3NCAsLM+677z5j+fLltj4//PCDIcn47LPPrvma5NRyuerVq9u9Tpf/PP7444ZhXPt9vlzOZzCvHxcXlwLVmtf3+HL9+vUzXFxcjIMHD151G2lpacaIESOMKlWqGJ6enkbLli2NTZs2GW3atDHatGlj65dT06effmqMHj3aCA4ONjw9PY3OnTvn+pw4+n7md5yOfNf/+9//Gvfff7/h7+9v+Pn5GQ8++KBx8uRJu/klx7Fjx4w+ffoYQUFBhsViMWrUqGEMHjzYSE9PNwzDMBYuXGhIMrZs2ZJnnVd+93/44QcjOjra8PPzMzw8PIzIyEijX79+xtatWw3DMIyzZ88agwcPNm655RbD29vb8PPzM5o1a2YsW7bsqu8HAKBsMxlGIVZsBQAA6tevn5YvX2672xoA5Fi/fr3uuusuffbZZ+revfs1+7Zt21Znz54ttjWcAAAoK1hTCgAAAAAAAE5HKAUAAAAAAACnI5QCAAAAAACA07GmFAAAAAAAAJyOM6UAAAAAAADgdIRSAAAAAAAAcDrX0i7A2axWq06ePKkKFSrIZDKVdjkAAAAAAADlimEYSkpKUmhoqMzmq58PdcOFUidPnlTVqlVLuwwAAAAAAIBy7e+//9ZNN9101edvuFCqQoUKki69ML6+vqVcTdFYrVbFxcUpKCjomskjAODqmEsBoOiYSwGg6MrTXJqYmKiqVavaMpirueFCqZxL9nx9fctFKJWWliZfX98y/4EFgNLCXAoARcdcCgBFVx7n0vyWTSofRwkAAAAAAIAyhVAKAAAAAAAATkcoBQAAAAAAAKe74daUclR2drYyMzNLu4xrslqtyszMVFpaWrm53vR64ebmJhcXl9IuAwAAAACAcotQ6gqGYej06dOKj48v7VLyZRiGrFarkpKS8l08DAXn7++vypUr89oCAAAAAFACCKWukBNIBQcHy8vL67oOJAzDUFZWllxdXa/rOssawzCUmpqq2NhYSVKVKlVKuSIAAAAAAMofQqnLZGdn2wKpwMDA0i4nX4RSJcfT01OSFBsbq+DgYC7lAwAAAACgmLEQ0WVy1pDy8vIq5UpwPcj5HFzva4sBAAAAAFAWEUrlgbOOIPE5AAAAAACgJBFKAQAAAAAAwOkIpeB0X375pWrWrCkXFxcNGzbsqm0AAAAAAKD8IpQqB8xms0wm01V/XnnlldIu0c6TTz6p7t276++//9bEiROv2lYU69evl8lkUnx8fJG3BQAAAAAAih933ysHTp48aVv/aOnSpRo7dqz2799ve97Hx8f2u2EYys7Olqtr6bz1ycnJio2NVXR0tEJDQ6/aBgAAAAAAyjfOlCoHKleubPvx8/OTyWSyPd63b58qVKigb7/9Vk2aNJHFYtGGDRvUr18/devWzW47w4YNU9u2bW2PrVarJk+erIiICHl6eqphw4Zavnz5NWtJT0/XyJEjFRYWJm9vbzVr1kzr16+XdOnspQoVKkiS7r77bplMpqu2SdKGDRvUunVreXp6qmrVqho6dKhSUlLs9vXCCy+oatWqslgsqlmzpubPn6+jR4/qrrvukiRVrFhRJpNJ/fr1K/wLDAAAAAAAih2h1A1i1KhReu2117R37141aNDAoTGTJ0/Whx9+qDlz5mjPnj0aPny4HnnkEf34449XHTNkyBBt2rRJS5Ys0a5du/Tggw+qQ4cO+uuvv9SiRQvbGVyff/65Tp06ddW2Q4cOqUOHDnrggQe0a9cuLV26VBs2bNCQIUNs++rTp48+/fRTvfXWW9q7d6/mzp0rHx8fVa1aVZ9//rkkaf/+/Tp16pRmzpxZ2JcOAAAAAACUAC7fc9CZn4/rzIbj+fbzCqugmn0a2rUd/HCnUk8k5Ts2pFU1hbSuVugar2XChAlq3769w/3T09M1adIkff/992revLkkqUaNGtqwYYPmzp2rNm3a5Bpz/PhxLVy4UMePH7ddhjdy5EitWrVKCxcu1KRJkxQcHCxJCggIUOXKlSUpz7bJkyerd+/etkXPa9Wqpbfeektt2rTRu+++q+PHj2vZsmVas2aN2rVrZ6svR0BAgG3b/v7+Dh83AAAAAABwDkIpB2WnZykzMT3ffll+ltxtyRkOjc1OzypUbY5o2rRpgfofPHhQqampuYKsjIwMRUVF5Tnmjz/+UHZ2tm6++Wa79vT0dAUGBhZo/zt37tSuXbu0ePFiW5thGLJarTpy5Ij++OMPubi45BmOAQAAAACA6x+hlINcLK5y880dOF3J1cc9zzZHxrpYSu7t8Pb2tntsNptlGIZdW2Zmpu335ORkSdLKlSsVFhZm189iyftYkpOT5eLiom3btsnFxcXuucsXW3dEcnKynnzySQ0dOjTXc9WqVdPBgwcLtD0AAAAAAHB9IZRyUEjrwl9ad+XlfNeDoKAg7d69265tx44dcnNzkyTVrVtXFotFx48fd/hspKioKGVnZys2NlatW7cuUn2NGzfWn3/+qZo1a+b5fP369WW1WvXjjz/aLt+7nLv7pXAwOzu7SHUAAAAAAICSQSh1g7r77rv15ptv6sMPP1Tz5s318ccfa/fu3bZL8ypUqKCRI0dq+PDhslqtatWqlRISErRx40b5+vqqb9++ubZ58803q3fv3urTp4+mTp2qqKgoxcXFae3atWrQoIE6d+7scH0vvPCC7rjjDg0ZMkQDBgyQt7e3/vzzT61Zs0azZs1SeHi4+vbtq8cee0xvvfWWGjZsqGPHjik2NlY9evRQ9erVZTKZtGLFCnXq1Emenp4FPlsLAAAAAIDidm7VNJ1bNc2uzTAMybAqOztbCS4uksksk8lk1yewQ4wCO8Q4s9QSx933blDR0dEaM2aMnn/+ed12221KSkpSnz597PpMnDhRY8aM0eTJk1WnTh116NBBK1euVERExFW3u3DhQvXp00cjRoxQ7dq11a1bN23ZskXVqhXsLLMGDRroxx9/1IEDB9S6dWtFRUVp7NixtgXUJendd99V9+7d9fTTT+uWW27RwIEDlZKSIkkKCwvT+PHjNWrUKIWEhNjdtQ8AAAAAgNKSfTFRWRdO2P1kx59UdsJpKTlO2QmnlR1/Mnefi4mlXXqxMxlXLixUziUmJsrPz08JCQny9fW1ey4tLU1HjhxRRESEPDw8SqlCxxmGoaysLLm6uuZKUFF0Ze3zAKBwrFarYmNjFRwcLLOZv9UAQGEwlwKA43LOlLJmpsmafF7StWIZk8w+ATK7eZSpM6Wulb1cjsv3AAAAAAAAnCSwQ4wsYbfq+LTOkskkXetcIZNJ1tR43RSzUj71o51XpJPwZwwAAAAAAAAnyU6J199vP3ApjDKs1+5sWCXD0N9vP6DslHin1OdMhFIAAAAAAABOEr9xkYz01PwDqRyGVUZGquI3fliyhZUCQikAAAAAAAAnMAxD59e8Xaix59e8pfK2LDihFAAAAAAAgBNkJ59TZuwhXXtx8zwYhjJjDyk75XyJ1FVaCKUAAAAAAACcwJqWXLTxF5OKqZLrA6EUAAAAAACAE5g9fIo23rNCMVVyfXAt7QLKg2k/HtL0nw7btRmGIashGTJkkklmk2Qymez6DL+zhmLaRDqzVAAAAAAAUEpcfALlFhypzLjDl+6+5yiTSW5BNeTiHVByxZUCQqlikJiWpRMJaYUaBwAAAAAAbgwmk0kB7Z/RmcXDCzw2oP3QXCe7lHVcvlcMfD1cFebnoUAvN+X38TBJCvRyU5ifh3w9yARL09GjR2UymbRjx47SLgUAAAAAcIPwb9lXJouXZHIwkjGZZXL3kn/LPiVbWCkglCoGMW0i9X6Phoq/mKX8QkuTSYq/mKX3ezTk0j0AAAAAAG4w8T8vkMnNQzKsjg0wrDK5eSj+5wUlW1gpIJQqBvEXM9V90VYZurSO1LXkrDPVfdFWxV/MdE6BecjMLL19l7SMjIzSLgEAAAAAgDxlX0yUNflcgcZYk88p+2JiCVVUegilisGirX8rNSM730Aqh9WQUjOy9eHWv4u1DqvVqjfeeEM1a9aUxWJRtWrV9Oqrr9ouU1u6dKnatGkjDw8PLV68WJL0/vvvq06dOvLw8NAtt9yid955x+H9/fLLL2rUqJE8PDzUtGlTffnll7kuh9u9e7c6duwoHx8fhYSE6NFHH9XZs2dtz7dt21ZDhw7V888/r4CAAFWuXFmvvPKK3X7i4+M1YMAABQUFydfXV3fffbd27txpe/6VV15Ro0aN9P777ysiIkIeHh6SpFWrVqlVq1by9/dXYGCg7rvvPh06dKgQrywAAAAAAMXDxdNXrhXD5FoxTC5+VWT29JPMLvadzC4ye/rJxb/K//p6+pZOwSWIRY2KyDAMzdpwpFBj395wRM+0iii2hcpGjx6tefPmafr06WrVqpVOnTqlffv22Z4fNWqUpk6dqqioKFswNXbsWM2aNUtRUVHavn27Bg4cKG9vb/Xt2/ea+0pMTFSXLl3UqVMnffLJJzp27JiGDRtm1yc+Pl533323BgwYoOnTp+vixYt64YUX1KNHD61bt87Wb9GiRYqJidFvv/2mTZs2qV+/fmrZsqXat28vSXrwwQfl6empb7/9Vn5+fpo7d67uueceHThwQAEBl+48cPDgQX3++ef64osv5OJy6cuckpKimJgYNWjQQMnJyRo7dqzuv/9+7dixQ2YzeSwAAAAAwPkCO8QosEOMXZthGMpKOqvY/x5V8E3hcq1Qqdwtap4XQqkiOpeaoUPnUgs8zpB06FyqzqdmKtDbvch1JCUlaebMmZo1a5YtUIqMjFSrVq109OhRSdKwYcP0z3/+0zZm3Lhxmjp1qq0tIiJCf/75p+bOnZtvKPXJJ5/IZDJp3rx58vDwUN26dXXixAkNHDjQ1icn7Jo0aZKtbcGCBapataoOHDigm2++WZLUoEEDjRs3TpJUq1YtzZo1S2vXrlX79u21YcMGbd68WbGxsbJYLJKkKVOm6Msvv9Ty5cv1xBNPSLp0yd6HH36ooKAg274eeOABu5oXLFigoKAg/fnnn6pXr57jLy4AAAAAACXIZDLJxSdQLgHZcvEJvCECKYlQqsiS07OLND4pPatYQqm9e/cqPT1d99xzz1X7NG3a1PZ7SkqKDh06pMcff9wuSMrKypKfn1+++9u/f78aNGhgu1ROkm6//Xa7Pjt37tQPP/wgHx+fXOMPHTpkF0pdrkqVKoqNjbVtIzk5WYGBgXZ9Ll68aHcpXvXq1e0CKUn666+/NHbsWP322286e/asrNZLi8gdP36cUAoAAAAAgFJGKFVEPhaX/DtdQwVL8bwFnp6e+fbx9va2/Z6cnCxJmjdvnpo1a2bXL+fyt6JKTk5Wly5d9Prrr+d6rkqVKrbf3dzc7J4zmUy2ACk5OVlVqlTR+vXrc23D39/f9vvlx5ajS5cuql69uubNm6fQ0FBZrVbVq1ePhdABAAAAALgOEEoVUaCXuyIDvXT4XKocXOdckmSSVCPQSwFebvn2dUStWrXk6emptWvXasCAAfn2DwkJUWhoqA4fPqzevXsXeH+1a9fWxx9/rPT0dNtldVu2bLHr07hxY33++ecKDw+Xq2vhPmqNGzfW6dOn5erqqvDwcIfHnTt3Tvv379e8efPUunVrSdKGDRsKVQMAAAAAACh+rPZcRCaTSUNaRRRqbHEucu7h4aEXXnhBzz//vD788EMdOnRIv/76q+bPn3/VMePHj9fkyZP11ltv6cCBA/rjjz+0cOFCTZs2Ld/9Pfzww7JarXriiSe0d+9erV69WlOmTJEk2zENHjxY58+fV69evbRlyxYdOnRIq1evVv/+/ZWd7dhlj+3atVPz5s3VrVs3fffddzp69Kh++eUXvfTSS9q6detVx1WsWFGBgYF67733dPDgQa1bt04xMTFX7Q8AAAAAAJyLUKoY9G1aVV7uLjI7mC+ZTZKXu4v6NK1arHWMGTNGI0aM0NixY1WnTh317NnTtjZTXgYMGKD3339fCxcuVP369dWmTRt98MEHiojIP2Tz9fXV119/rR07dqhRo0Z66aWXNHbsWEmyrTMVGhqqjRs3Kjs7W/fee6/q16+vYcOGyd/f3+G735lMJn3zzTe688471b9/f91888166KGHdOzYMYWEhFx1nNls1pIlS7Rt2zbVq1dPw4cP15tvvunQPgEAAAAAQMkzGYZRkKvOyrzExET5+fkpISFBvr6+ds+lpaXpyJEjioiIsFvA2xGr98fqvvc3y5Ah6zVeUbNJMsmklQNu1721gwtzCDaGYSgrK0uurq7Xxcr8ixcvVv/+/ZWQkODQGlfXu6J8HgCUHVarVbGxsQoODnY4MAcA2GMuBYCiK09z6bWyl8uV7aO8Tkz78ZAGLNspf09X5RfxGYbk7+mqx5ft1LQfD12783Xuww8/1IYNG3TkyBF9+eWXeuGFF9SjR49yEUgBAAAAAICSRShVDBLTsnQiIU3nUjPzXezckHQuNVMnEtKUmJbljPIKZdKkSfLx8cnzp2PHjpKk06dP65FHHlGdOnU0fPhwPfjgg3rvvfdKuXIAAAAAAFAWcPe9YuDr4aowP/vLuwzj0mV8hgyZZLp02d4Vl9j5ely/L/+gQYPUo0ePPJ/LORPq+eef1/PPP+/MsgAAAAAAQDlRqqnITz/9pDfffFPbtm3TqVOn9O9//1vdunW7av9Tp05pxIgR2rp1qw4ePKihQ4dqxowZTqv3amLaRCqmTWRpl1GsAgICFBAQUNplAAAAAACAcqpUL99LSUlRw4YNNXv2bIf6p6enKygoSC+//LIaNmxYwtUBAAAAAACgpJTqmVIdO3a0rU/kiPDwcM2cOVOStGDBgpIqCwAAAAAAACXs+l3UqJikp6crPT3d9jgxMVHSpVstWq1Wu75Wq1WGYdh+yoKcOstKvWVJzucgr88KgPIjZ+7new4AhcdcCgBFV57mUkePodyHUpMnT9b48eNztcfFxSktLc2uLTMzU1arVVlZWcrKcvzOeEk73lLyjrft2i6FRFbJMCSTSZI510LnPo2eUYVGQx3ez5UMw1B2drak3Iuoo+iysrJktVp17tw5ubm5lXY5AEqI1WpVQkKCDMOQ2cxNaQGgMJhLAaDoytNcmpSU5FC/ch9KjR49WjExMbbHiYmJqlq1qoKCguTr62vXNy0tTUlJSXJ1dZWrawFemqwUZaecLHhxWSkF289VEJiUDFdXV5nNZgUGBsrDwyP/AQDKJKvVKpPJpKCgoDL/H38AKC3MpQBQdOVpLnX039DlPpSyWCyyWCy52s1mc6432Wy+dDZTzo+jXCy+cvEJk5GVJmvaeUnXupTOJLNHgEyuHnKx+BbpDCfDMGzjOVMqf+vXr9ddd92lCxcuyN/fP9/+OZ+DvD4rAMoXvusAUHTMpQBQdOVlLnW0/nIfSjmDX+NhcguoqzNfdZVkUn6hlDU9XiEdvpJX+L1OqhAAAAAAAOD6UqrRW3Jysnbs2KEdO3ZIko4cOaIdO3bo+PHjki5detenTx+7MTn9k5OTFRcXpx07dujPP/90dul2stPiFbuy56X1o5TfYl6X1pmKXdlT2WnxTqgub5mZmaW272vJyMgo7RIAAAAAAIATlGootXXrVkVFRSkqKkqSFBMTo6ioKI0dO1aSdOrUKVtAlSOn/7Zt2/TJJ58oKipKnTp1cnrtl0ve+5GMzFTlH0jlsMrITFXy3o+LtQ6r1ao33nhDNWvWlMViUbVq1fTqq6/q6NGjMplMWrp0qdq0aSMPDw8tXrxYkvT++++rTp068vDw0C233KJ33nnHoX3lbPOLL77QXXfdJS8vLzVs2FCbNm2y6/f555/r1ltvlcViUXh4uKZOnWr3fHh4uCZOnKg+ffrI19dXTzzxhD744AP5+/trxYoVql27try8vNS9e3elpqZq0aJFCg8PV8WKFTV06FDbQu+S9NFHH6lp06aqUKGCKleurIcfflixsbFFfFUBAAAAAEBJKNXL99q2bfv/d6nL2wcffJCr7Vr9S4NhGErcMbtQYxN3zJJvo8HFth7U6NGjNW/ePE2fPl2tWrXSqVOntG/fPtvzo0aN0tSpUxUVFWULpsaOHatZs2YpKipK27dv18CBA+Xt7a2+ffs6tM+XXnpJU6ZMUa1atfTSSy+pV69eOnjwoFxdXbVt2zb16NFDr7zyinr27KlffvlFTz/9tAIDA9WvXz/bNqZMmaKxY8dq3LhxkqSff/5Zqampeuutt7RkyRIlJSXpn//8p+6//375+/vrm2++0eHDh/XAAw+oZcuW6tmzp6RLZ39NnDhRtWvXVmxsrGJiYtSvXz998803xfL6AgAAAACA4sOaUkVkTTunrITDhRhpKCvhsKxp5+XiGVjkOpKSkjRz5kzNmjXLFihFRkaqVatWOnr0qCRp2LBh+uc//2kbM27cOE2dOtXWFhERoT///FNz5851OJQaOXKkOnfuLEkaP368br31Vh08eFC33HKLpk2bpnvuuUdjxoyRJN188836888/9eabb9qFUnfffbdGjBhhe/zzzz8rMzNT7777riIjIyVJ3bt310cffaQzZ87Ix8dHdevW1V133aUffvjBFko99thjtm3UqFFDb731lm677TYlJyfLx8enIC8nAAAAAAAoYWV7OffrgDUjuYjjk4qljr179yo9PV333HPPVfs0bdrU9ntKSooOHTqkxx9/XD4+Praff/3rXzp06JDD+23QoIHt9ypVqkiS7ZK5vXv3qmXLlnb9W7Zsqb/++svusrvL68rh5eVlC6QkKSQkROHh4XbhUkhIiN3ledu2bVOXLl1UrVo1VahQQW3atJGkXJeAAgAAAACA0seZUkVkdi/aGThm9wrFUoenp2e+fby9vW2/JydfCtPmzZunZs2a2fVzcXFxeL9ubm6233MuQ7RaHV1bK3ddeW03Z9t5teXsKyUlRdHR0YqOjtbixYsVFBSk48ePKzo6msXTAQAAAAC4DhFKFZHZI1CufjWUlXBEUkHWuzLJ1S9CZo+AYqmjVq1a8vT01Nq1azVgwIB8+4eEhCg0NFSHDx9W7969i6WGK9WpU0cbN260a9u4caNuvvnmAgVfjti3b5/OnTun1157TVWrVpV0aSF9AAAAAABwfSKUKiKTySTfRoN1/seRBR7r22hIsS1y7uHhoRdeeEHPP/+83N3d1bJlS8XFxWnPnj1XvaRv/PjxGjp0qPz8/NShQwelp6dr69atunDhgmJiYopc04gRI3Tbbbdp4sSJ6tmzpzZt2qRZs2Y5fIe/gqhWrZrc3d319ttva9CgQdq9e7cmTpxY7PsBAAAAAADFgzWlioFPnUdlcvOS4y+nWSY3L/nUeaRY6xgzZoxGjBihsWPHqk6dOurZs6fdmktXGjBggN5//30tXLhQ9evXV5s2bfTBBx8oIiKiWOpp3Lixli1bpiVLlqhevXoaO3asJkyYYLfIeXEJCgrSBx98oM8++0x169bVa6+9pilTphT7fgAAAAAAQPEwGYZRkGvOyrzExET5+fkpISFBvr6+ds+lpaXpyJEjioiIkIeHR4G2m3r0O535qqtkGJKutaaSWTKZFNLtP/Kq3r7gB3AZwzCUlZUlV1fXYjvjCv9TlM8DgLLDarUqNjZWwcHBMpv5Ww0AFAZzKQAUXXmaS6+VvVyubB/ldSLh9xk6+/2TMlv8lf+6UobMFn+dXfOEEn6fUfLFAQAAAAAAXIcIpYqBNT1R2cknZE07J0dCKWvauUv90xOdUV6hTJo0ST4+Pnn+dOzYsbTLAwAAAAAAZRwLnRcDs8VXLj5hdm1GzmV8hiGZTJLMuS6xM1uufgpbaRs0aJB69OiR53Oenp5OrgYAAAAAAJQ3hFLFwK/xMPk1HlbaZRSrgIAABQQElHYZAAAAAACgnOLyPQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATseaUsVg+u4fNX3PT3ZthmHIKsO2zrlZplwLnQ+/9U4Nr9fGmaUCAAAAAABcFwilikFiZppOpCYUahwAAAAAAMCNiMv3ioGvm4fCvPwUaPGSKZ++JkmBFi+FefnJ182j2GowDENPPPGEAgICZDKZtGPHjlx92rZtq2HDhhXbPgEAAAAAAAqLUKoYDK/XRvNa9VB8xkWZ8omlTDIpPuOi5rXqUayX7q1atUoffPCBVqxYoVOnTqlevXrFtu3ScPTo0auGawAAAAAAoOwjlCoG8ekX9eC6RTIMySrjmn1z1pl6cN0ixadfLLYaDh06pCpVqqhFixaqXLmyXF25MhMAAAAAAFy/CKWKwYcHtyo1KyPfQCqHVYZSszL00aGtxbL//v3765lnntHx48dlMpkUHh6ulJQU9enTRz4+PqpSpYqmTp1aoG2mp6dr5MiRCgsLk7e3t5o1a6b169dLktLS0nTrrbfqiSeesPU/dOiQKlSooAULFkiSzp07p169eiksLExeXl6qX7++Pv30U7t9WK1WvfHGG6pZs6YsFouqVaumV199VZIUEREhSYqKipLJZFLbtm0L+eoAAAAAAIDrEaFUERmGoVl7NxRq7Nt/bpBhOBZkXcuMGTM0YcIE3XTTTTp16pS2bNmi5557Tj/++KO++uorfffdd1q/fr1+//13h7c5ZMgQbdq0SUuWLNGuXbv04IMPqkOHDvrrr7/k4eGhxYsXa9GiRfrqq6+UnZ2tRx55RO3bt9djjz0m6VJw1aRJE61cuVK7d+/WE088oUcffVSbN2+27WP06NF67bXXNGbMGP3555/65JNPFBISIkm2ft9//71OnTqlL774osivEwAAAAAAuH5wjVcRnUtP1aGkcwUeZ0g6lHRO59NTFejhXaQa/Pz8VKFCBbm4uKhy5cpKTk7W/Pnz9fHHH+uee+6RJC1atEg33XSTQ9s7fvy4Fi5cqOPHjys0NFSSNHLkSK1atUoLFy7UpEmT1KhRI/3rX//SgAED9NBDD+nYsWNasWKFbRthYWEaOXKk7fEzzzyj1atXa9myZbr99tuVlJSkmTNnatasWerbt68kKTIyUq1atZIkBQUFSZICAwNVuXLlIr0+AAAAAADg+kMoVUTJmelFGp+UmV7kUOpKhw4dUkZGhpo1a2ZrCwgIUO3atR0a/8cffyg7O1s333yzXXt6eroCAwNtj0eMGKEvv/xSs2bN0rfffmv3XHZ2tiZNmqRly5bpxIkTysjIUHp6ury8vCRJe/fuVXp6ui00AwAAAAAANxZCqSLycbMUaXyFIo4vCcnJyXJxcdG2bdvk4uJi95yPj4/t99jYWB04cEAuLi7666+/1KFDB9tzb775pmbOnKkZM2aofv368vb21rBhw5SRkSFJ8vT0dM7BAAAAAACA6xJrShVRoMVLkRUCZSrgOJOkyAqBCrB4FXtNkZGRcnNz02+//WZru3Dhgg4cOODQ+KioKGVnZys2NlY1a9a0+7n8UrrHHntM9evX16JFi/TCCy9o7969tuc2btyorl276pFHHlHDhg1Vo0YNu/3XqlVLnp6eWrt2bZ41uLu7S7p0xhUAAAAAACh/OFOqiEwmk4bUaaWYzV8VeOwzdVvJZCponJU/Hx8fPf7443ruuecUGBio4OBgvfTSSzKbHcsgb775ZvXu3Vt9+vTR1KlTFRUVpbi4OK1du1YNGjRQ586dNXv2bG3atEm7du1S1apVtXLlSvXu3Vu//vqr3N3dVatWLS1fvly//PKLKlasqGnTpunMmTOqW7euJMnDw0MvvPCCnn/+ebm7u6tly5aKi4vTnj179Pjjjys4OFienp5atWqVbrrpJnl4eMjPz6/YXysAAAAAAFA6OFOqGPSp2VReru4yO3i+lFkmebm669HIpiVW05tvvqnWrVurS5cuateunVq1aqUmTZo4PH7hwoXq06ePRowYodq1a6tbt27asmWLqlWrpn379um5557TO++8o6pVq0qS3nnnHZ09e1ZjxoyRJL388stq3LixoqOj1bZtW1WuXFndunWz28eYMWM0YsQIjR07VnXq1FHPnj0VGxsrSXJ1ddVbb72luXPnKjQ0VF27di2eFwYAAAAAAFwXTIZhGKVdhDMlJibKz89PCQkJ8vX1tXsuLS1NR44cUUREhDw8PAq03dUn9qvLmvdlGJJVV39JzTLJZJJWtB+ge8McW3j8agzDUFZWllxdXUvkjKsbXVE+DwDKDqvVqtjYWAUHBzt8RikAwB5zKQAUXXmaS6+VvVyubB/ldWL67h81cMMy+bt7yrhGICVJhgz5u3tqwIZlmr77RydVCAAAAAAAcH1hTalikJiZphOpCQ71NSSdS0+1jSsNP//8szp27HjV55OTk51YDQAAAAAAuBERShUDXzcPhXnZL8JtGIasMmQYksmUc9meKde40tC0aVPt2LGjVPYNAAAAAAAgEUoVi+H12mh4vTalXYbDPD09VbNmzdIuAwAAAAAA3MBYUwoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNOxplQxOLdqms6tmmbXZhiGZFhlW+ncZM610HlghxgFdohxZqkAAAAAAADXBUKpYpB9MVFZF04UahwAAAAAAMCNiFCqGLh4+sq1YpismWmyJp+XZFyjt0lmnwCZ3Tzk4unrrBIBAAAAAACuK4RSxSCwQ4wsYbfq+LTOly7VM64RSplMsqbG66aYlfKpH+28IgEAAAAAAK4jLHReDLJT4vX32w9cCqMM67U7//86U3+//YCyU+KdUl95kZmZWdolAAAAAACAYkIoVQziNy6SkZ6afyCVw7DKyEhV/MYPi2X/K1askL+/v7KzsyVJO3bskMlk0qhRo2x9BgwYoEceeUSSNG/ePFWtWlVeXl66//77NW3aNPn7+zu0r0OHDqlr164KCQmRj4+PbrvtNn3//fd2fcLDwzVx4kT16tVL3t7eCgsL0+zZs+36mEwmvfvuu+rYsaM8PT1Vo0YNLV++3Pb80aNHZTKZtHTpUrVp00YeHh5avHixrFarJkyYoJtuukkWi0WNGjXSqlWrJF1aXL5du3aKjo6+tNC8pPPnz+umm27S2LFjC/aiAgAAAACAEkUoVUSGYej8mrcLNfb8mrds4UlRtG7dWklJSdq+fbsk6ccff1SlSpW0fv16W58ff/xRbdu21caNGzVo0CA9++yz2rFjh9q3b69XX33V4X0lJyerU6dOWrt2rbZv364OHTqoS5cuOn78uF2/N998Uw0bNtT27ds1atQoPfvss1qzZo1dnzFjxuiBBx7Qzp071bt3bz300EPau3evXZ+csXv37lV0dLRmzpypqVOnasqUKdq1a5eio6P1j3/8Q3/99ZdMJpMWLVqkLVu26K233pIkDRo0SGFhYYRSAAAAAABcZ0xGcaQiZUhiYqL8/PyUkJAgX1/7hcbT0tJ05MgRRUREyMPDw6HtZSWd1YEhQYWu5+bZZ+XqE1iosYZhKCsrS66urmratKl69eqlkSNH6v7779dtt92m8ePH69y5c0pISNBNN92kAwcOaMyYMUpOTtaKFSts23nkkUe0YsUKxcfHF6qOevXqadCgQRoyZIikS2dK1alTR99++62tz0MPPaTExER98803ki6dKTVo0CC9++67tj533HGHGjdurHfeeUdHjx5VRESEZsyYoWeffdbWJywsTIMHD9aLL75oa7v99tt122232c7G+uyzz9SnTx8NGzZMb7/9trZv365atWoV+LgK83kAUPZYrVbFxsYqODhYZjN/qwGAwmAuBYCiK09z6bWyl8uV7aO8DljTkos2/mJSsdTRpk0brV+/XoZh6Oeff9Y///lP1alTRxs2bNCPP/6o0NBQ1apVS/v379ftt99uN/bKx9eSnJyskSNHqk6dOvL395ePj4/27t2b60yp5s2b53p85VlQjvRp2rSp7ffExESdPHlSLVu2tOvTsmVLu3EPPvig7r//fr322muaMmVKoQIpAAAAAABQsrj7XhGZPXyKNt6zQrHU0bZtWy1YsEA7d+6Um5ubbrnlFrVt21br16/XhQsX1KZNm2LZz8iRI7VmzRpNmTJFNWvWlKenp7p3766MjIxi2f6VvL29CzwmNTVV27Ztk4uLi/76668SqAoAAAAAABQVZ0oVkYtPoNyCIyWTqWADTSa5BUfKxTugWOrIWVdq+vTptgAqJ5Rav3692rZtK0mqXbu2tmzZYjf2ysfXsnHjRvXr10/333+/6tevr8qVK+vo0aO5+v3666+5HtepU6fAfS7n6+ur0NBQbdy4MVdNdevWtT0eMWKEzGazvv32W7311ltat26do4cHAAAAAACchDOlishkMimg/TM6s3h4gccGtB8qU0HDrKuoWLGiGjRooMWLF2vWrFmSpDvvvFM9evRQZmamLah65plndOedd2ratGnq0qWL1q1bp2+//dbhOmrVqqUvvvhCXbp0kclk0pgxY2S15r7r4MaNG/XGG2+oW7duWrNmjT777DOtXLnSrs9nn32mpk2bqlWrVlq8eLE2b96s+fPnX3P/zz33nMaNG6fIyEg1atRICxcu1I4dO7R48WJJ0sqVK7VgwQJt2rRJjRs31nPPPae+fftq165dqlixokPHCAAAAAAASh5nShUD/5Z9ZbJ4SSYHX06TWSZ3L/m37FOsdbRp00bZ2dm2s6ICAgJUt25dVa5cWbVr15Z0af2lOXPmaNq0aWrYsKFWrVql4cOHO7yQ97Rp01SxYkW1aNFCXbp0UXR0tBo3bpyr34gRI7R161ZFRUXpX//6l6ZNm6bo6Gi7PuPHj9eSJUvUoEEDffjhh/r000/tznjKy9ChQxUTE6MRI0aofv36WrVqlf7zn/+oVq1aiouL0+OPP65XXnnFVtP48eMVEhKiQYMGOXR8AAAAAADAObj73mWKcre15D9W6/i0zpJhSEbuM4dsTGbJZFK1mG/kU//ewhyCzeV33yvKGVcDBw7Uvn379PPPPxepnhzh4eEaNmyYhg0bdtU+JpNJ//73v9WtW7di2WdJ4O57wI2hPN3lBABKC3MpABRdeZpLufueE51bNU0n5z8us5f/pVDqWgxDZi9/nZz/mM6tmuaU+q40ZcoU7dy5UwcPHtTbb7+tRYsWqW/fvqVSCwAAAAAAuDGxplQxyL6YqKwLJxzsbciafE7W/x9XGjZv3qw33nhDSUlJqlGjht566y0NGDBAknTrrbfq2LFjeY6bO3euevfu7cxSAQAAAABAOUUoVQxcPH3lWjHMrs3IuYzPMC7dmc9kznWJnYvn1U9hK0nLli276nPffPONMjMz83wuJCTEoe3ndTe+K91gV40CAAAAAIArEEoVg8AOMQrsEFPaZRSL6tWrl3YJAAAAAADgBsCaUgAAAAAAAHC6Ug2lfvrpJ3Xp0kWhoaEymUz68ssv8x2zfv16NW7cWBaLRTVr1tQHH3xQ7HVxaRkkPgcAAAAAAJSkUg2lUlJS1LBhQ82ePduh/keOHFHnzp111113aceOHRo2bJgGDBig1atXF0s9bm5ukqTU1NRi2R7KtpzPQc7nAgAAAAAAFJ9SXVOqY8eO6tixo8P958yZo4iICE2dOlWSVKdOHW3YsEHTp09XdHR0ketxcXGRv7+/YmNjJUleXl65Fie/nhiGoaysLLm6ul7XdZY1hmEoNTVVsbGx8vf3l4uLS2mXBAAAAABAuVOmFjrftGmT2rVrZ9cWHR2tYcOGXXVMenq60tPTbY8TExMlSVarVVarNVf/4OBgGYahM2fOFE/RJcxqtcpsZmmwkuDv76/g4OA8PycAyg+r1SrDMPiuA0ARMJcCQNGVp7nU0WMoU6HU6dOnFRISYtcWEhKixMREXbx4UZ6enrnGTJ48WePHj8/VHhcXp7S0tDz3Yzab5e/vr+zs7OIpvIQYhqGkpCT5+PhwplQxc3FxkdlsVlxcXGmXAqCEWa1WJSQkyDAMQn4AKCTmUgAouvI0lyYlJTnUr0yFUoUxevRoxcTE2B4nJiaqatWqCgoKkq+vbylWVnRWq1VxcXEKCgoq8x9YACgtVqtVJpOJuRQAioC5FACKrjzNpR4eHg71K1OhVOXKlXNdVnfmzBn5+vrmeZaUJFksFlksllztZrO5zL/JkmQymcrNsQBAaWEuBYCiYy4FgKIrL3Opo/WXqaNs3ry51q5da9e2Zs0aNW/evJQqAgAAAAAAQGGUaiiVnJysHTt2aMeOHZKkI0eOaMeOHTp+/LikS5fe9enTx9Z/0KBBOnz4sJ5//nnt27dP77zzjpYtW6bhw4eXRvkAAAAAAAAopFINpbZu3aqoqChFRUVJkmJiYhQVFaWxY8dKkk6dOmULqCQpIiJCK1eu1Jo1a9SwYUNNnTpV77//vqKjo0ulfgAAAAAAABROqa4p1bZtWxmGcdXnP/jggzzHbN++vQSrAgAAAAAAQEkrU2tKAQAAAAAAoHwglAIAAAAAAIDTEUoBAAAAAADA6QilAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA4HaEUAAAAAAAAnI5QCgAAAAAAAE5HKAUAAAAAAACnI5QCAAAAAACA0xFKAQAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOB2hFAAAAAAAAJyOUAoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNMRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTEUoBAAAAAADA6QilAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA4HaEUAAAAAAAAnI5QCgAAAAAAAE5HKAUAAAAAAACnI5QCAAAAAACA0xFKAQAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOB2hFAAAAAAAAJyOUAoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNMRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTEUoBAAAAAADA6QilAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA43XURSs2ePVvh4eHy8PBQs2bNtHnz5qv2zczM1IQJExQZGSkPDw81bNhQq1atcmK1AAAAAAAAKKpSD6WWLl2qmJgYjRs3Tr///rsaNmyo6OhoxcbG5tn/5Zdf1ty5c/X222/rzz//1KBBg3T//fdr+/btTq4cAAAAAAAAhVXqodS0adM0cOBA9e/fX3Xr1tWcOXPk5eWlBQsW5Nn/o48+0osvvqhOnTqpRo0aeuqpp9SpUydNnTrVyZUDAAAAAACgsFxLc+cZGRnatm2bRo8ebWszm81q166dNm3alOeY9PR0eXh42LV5enpqw4YNV+2fnp5ue5yYmChJslqtslqtRT2EUmW1WmUYRpk/DgAoTcylAFB0zKUAUHTlaS519BhKNZQ6e/assrOzFRISYtceEhKiffv25TkmOjpa06ZN05133qnIyEitXbtWX3zxhbKzs/PsP3nyZI0fPz5Xe1xcnNLS0op+EKXIarUqISFBhmHIbC71k94AoExiLgWAomMuBYCiK09zaVJSkkP9SjWUKoyZM2dq4MCBuuWWW2QymRQZGan+/ftf9XK/0aNHKyYmxvY4MTFRVatWVVBQkHx9fZ1VdomwWq0ymUwKCgoq8x9YACgtzKUAUHTMpQBQdOVpLr3yCrerKdVQqlKlSnJxcdGZM2fs2s+cOaPKlSvnOSYoKEhffvml0tLSdO7cOYWGhmrUqFGqUaNGnv0tFossFkuudrPZXObfZEkymUzl5lgAoLQwlwJA0TGXAkDRlZe51NH6S/Uo3d3d1aRJE61du9bWZrVatXbtWjVv3vyaYz08PBQWFqasrCx9/vnn6tq1a0mXCwAAAAAAgGJS6pfvxcTEqG/fvmratKluv/12zZgxQykpKerfv78kqU+fPgoLC9PkyZMlSb/99ptOnDihRo0a6cSJE3rllVdktVr1/PPPl+ZhAAAAAAAAoABKPZTq2bOn4uLiNHbsWJ0+fVqNGjXSqlWrbIufHz9+3O60r7S0NL388ss6fPiwfHx81KlTJ3300Ufy9/cvpSMAAAAAAABAQZkMwzBKuwhnSkxMlJ+fnxISEsrFQuexsbEKDg4u89ebAkBpYS4FgKJjLgWAoitPc6mj2UvZPkoAAAAAAACUSYRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATleoUOrQoUN6+eWX1atXL8XGxkqSvv32W+3Zs6dYiwMAAAAAAED5VOBQ6scff1T9+vX122+/6YsvvlBycrIkaefOnRo3blyxFwgAAAAAAIDyp8Ch1KhRo/Svf/1La9askbu7u6397rvv1q+//lqsxQEAAAAAAKB8KnAo9ccff+j+++/P1R4cHKyzZ88WS1EAAAAAAAAo3wocSvn7++vUqVO52rdv366wsLBiKQoAAAAAAADlW4FDqYceekgvvPCCTp8+LZPJJKvVqo0bN2rkyJHq06dPSdQIAAAAAACAcqbAodSkSZN0yy23qGrVqkpOTlbdunV15513qkWLFnr55ZdLokYAAAAAAACUM64FHeDu7q558+ZpzJgx2r17t5KTkxUVFaVatWqVRH0AAAAAAAAohwocSuWoVq2aqlWrVpy1AAAAAAAA4AZR4FDqscceu+bzCxYsKHQxAAAAAAAAuDEUOJS6cOGC3ePMzEzt3r1b8fHxuvvuu4utMAAAAAAAAJRfBQ6l/v3vf+dqs1qteuqppxQZGVksRQEAAAAAAKB8K/Dd9/LciNmsmJgYTZ8+vTg2BwAAAAAAgHKuWEIpSTp06JCysrKKa3MAAAAAAAAoxwp8+V5MTIzdY8MwdOrUKa1cuVJ9+/YttsIAAAAAAABQfhU4lNq+fbvdY7PZrKCgIE2dOjXfO/MBAAAAAAAAUiFCqR9++KEk6gAAAAAAAMANpNjWlAIAAAAAAAAc5dCZUlFRUTKZTA5t8Pfffy9SQQAAAAAAACj/HAqlunXrVsJlAAAAAAAA4EbiUCg1bty4kq4DAAAAAAAANxDWlAIAAAAAAIDTFfjue9nZ2Zo+fbqWLVum48ePKyMjw+758+fPF1txAAAAAAAAKJ8KfKbU+PHjNW3aNPXs2VMJCQmKiYnRP//5T5nNZr3yyislUCIAAAAAAADKmwKHUosXL9a8efM0YsQIubq6qlevXnr//fc1duxY/frrryVRIwAAAAAAAMqZAodSp0+fVv369SVJPj4+SkhIkCTdd999WrlyZfFWBwAAAAAAgHKpwKHUTTfdpFOnTkmSIiMj9d1330mStmzZIovFUrzVAQAAAAAAoFwqcCh1//33a+3atZKkZ555RmPGjFGtWrXUp08fPfbYY8VeIAAAAAAAAMofh+++N2vWLD3yyCN67bXXbG09e/ZUtWrVtGnTJtWqVUtdunQpkSIBAAAAAABQvjh8ptRLL72k0NBQ9e7dW+vWrbO1N2/eXDExMQRSAAAAAAAAcJjDodTp06c1Z84cnTx5Uu3bt1dERIQmTpyov//+uyTrAwAAAAAAQDnkcCjl6empPn366IcfftBff/2lRx99VPPnz1dERIQ6dOigzz77TJmZmSVZKwAAAAAAAMqJAi90Lkk1atTQhAkTdOTIEX377bcKDAxUv379FBYWVtz1AQAAAAAAoBwqVCiVw2QyydXVVSaTSYZhcKYUAAAAAAAAHFKoUOrvv//WhAkTVKNGDbVv314nT57UvHnzdOrUqeKuDwAAAAAAAOWQq6MdMzIy9MUXX2jBggVat26dqlSpor59++qxxx5TjRo1SrJGAAAAAAAAlDMOh1KVK1dWamqq7rvvPn399deKjo6W2Vykq/8AAAAAAABwg3I4lHr55Zf16KOPKigoqCTrAQAAAAAAwA3A4VAqJiamJOsAAAAAAADADYTr7wAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOF2BQ6kJEyYoNTU1V/vFixc1YcKEYikKAAAAAAAA5VuBQ6nx48crOTk5V3tqaqrGjx9fLEUBAAAAAACgfCtwKGUYhkwmU672nTt3KiAgoFiKAgAAAAAAQPnm6mjHihUrymQyyWQy6eabb7YLprKzs5WcnKxBgwaVSJEAAAAAAAAoXxwOpWbMmCHDMPTYY49p/Pjx8vPzsz3n7u6u8PBwNW/evESKBAAAAAAAQPnicCjVt29fSVJERIRatGghNze3EisKAAAAAAAA5ZvDoVSONm3ayGq16sCBA4qNjZXVarV7/s477yy24gAAAAAAAFA+FXih819//VU1a9ZUnTp1dOedd6pt27a2n7vuuqtQRcyePVvh4eHy8PBQs2bNtHnz5mv2nzFjhmrXri1PT09VrVpVw4cPV1paWqH2DQAAAAAAAOcr8JlSgwYNUtOmTbVy5UpVqVIlzzvxFcTSpUsVExOjOXPmqFmzZpoxY4aio6O1f/9+BQcH5+r/ySefaNSoUVqwYIFatGihAwcOqF+/fjKZTJo2bVqRagEAAAAAAIBzFDiU+uuvv7R8+XLVrFmzWAqYNm2aBg4cqP79+0uS5syZo5UrV2rBggUaNWpUrv6//PKLWrZsqYcffliSFB4erl69eum3334rlnoAAAAAAABQ8gocSjVr1kwHDx4sllAqIyND27Zt0+jRo21tZrNZ7dq106ZNm/Ic06JFC3388cfavHmzbr/9dh0+fFjffPONHn300Tz7p6enKz093fY4MTFRkmS1WnOth1XWWK1WGYZR5o8DAEoTcykAFB1zKQAUXXmaSx09hgKHUs8884xGjBih06dPq379+rnuwtegQQOHt3X27FllZ2crJCTErj0kJET79u3Lc8zDDz+ss2fPqlWrVjIMQ1lZWRo0aJBefPHFPPtPnjxZ48ePz9UeFxdX5tehslqtSkhIkGEYMpsLvDwYAEDMpQBQHJhLAaDoytNcmpSU5FC/AodSDzzwgCTpscces7WZTCYZhiGTyaTs7OyCbrJA1q9fr0mTJumdd96xnbX17LPPauLEiRozZkyu/qNHj1ZMTIztcWJioqpWraqgoCD5+vqWaK0lzWq1ymQyKSgoqMx/YAGgtDCXAkDRMZcCQNGVp7nUw8PDoX4FDqWOHDlS4GKuplKlSnJxcdGZM2fs2s+cOaPKlSvnOWbMmDF69NFHNWDAAElS/fr1lZKSoieeeEIvvfRSrjfOYrHIYrHk2o7ZbC7zb7J0KRAsL8cCAKWFuRQAio65FACKrrzMpY7WX+BQqnr16gUu5mrc3d3VpEkTrV27Vt26dZN0KRlcu3athgwZkueY1NTUXAfn4uIiSTIMo9hqAwAAAAAAQMkpVPT20UcfqWXLlgoNDdWxY8ckSTNmzNBXX31V4G3FxMRo3rx5WrRokfbu3aunnnpKKSkptrvx9enTx24h9C5duujdd9/VkiVLdOTIEa1Zs0ZjxoxRly5dbOEUAAAAAAAArm8FPlPq3Xff1dixYzVs2DC9+uqrtjWk/P39NWPGDHXt2rVA2+vZs6fi4uI0duxYnT59Wo0aNdKqVatsi58fP37c7syol19+WSaTSS+//LJOnDihoKAgdenSRa+++mpBDwUAAAAAAAClxGQU8Jq3unXratKkSerWrZsqVKignTt3qkaNGtq9e7fatm2rs2fPllStxSIxMVF+fn5KSEgoFwudx8bGKjg4uMxfbwoApYW5FACKjrkUAIquPM2ljmYvBT7KI0eOKCoqKle7xWJRSkpKQTcHAAAAAACAG1CBQ6mIiAjt2LEjV/uqVatUp06d4qgJAAAAAAAA5VyB15SKiYnR4MGDlZaWJsMwtHnzZn366aeaPHmy3n///ZKoEQAAAAAAAOVMgUOpAQMGyNPTUy+//LJSU1P18MMPKzQ0VDNnztRDDz1UEjUCAAAAAACgnClwKCVJvXv3Vu/evZWamqrk5GQFBwcXd10AAAAAAAAoxwoVSuXw8vKSl5dXcdUCAAAAAACAG4RDoVTjxo21du1aVaxYUVFRUTKZTFft+/vvvxdbcQAAAAAAACifHAqlunbtKovFIknq1q1bSdYDAAAAAACAG4BDodS4cePy/B0AAAAAAAAoDHNBB2zZskW//fZbrvbffvtNW7duLZaiAAAAAAAAUL4VOJQaPHiw/v7771ztJ06c0ODBg4ulKAAAAAAAAJRvBQ6l/vzzTzVu3DhXe1RUlP78889iKQoAAAAAAADlW4FDKYvFojNnzuRqP3XqlFxdHVqiCgAAAAAAADe4AodS9957r0aPHq2EhARbW3x8vF588UW1b9++WIsDAAAAAABA+VTgU5umTJmiO++8U9WrV1dUVJQkaceOHQoJCdFHH31U7AUCAAAAAACg/ClwKBUWFqZdu3Zp8eLF2rlzpzw9PdW/f3/16tVLbm5uJVEjAAAAAAAAyplCLQLl7e2tJ554orhrAQAAAAAAwA3CoVDqP//5jzp27Cg3Nzf95z//uWbff/zjH8VSGAAAAAAAAMovh0Kpbt266fTp0woODla3bt2u2s9kMik7O7u4agMAAAAAAEA55VAoZbVa8/wdAAAAAAAAKAyzI50CAgJ09uxZSdJjjz2mpKSkEi0KAAAAAAAA5ZtDoVRGRoYSExMlSYsWLVJaWlqJFgUAAAAAAIDyzaHL95o3b65u3bqpSZMmMgxDQ4cOlaenZ559FyxYUKwFAgAAAAAAoPxxKJT6+OOPNX36dB06dEiSlJCQwNlSAAAAAAAAKDSHQqmQkBC99tprkqSIiAh99NFHCgwMLNHCAAAAAAAAUH4VeKHzu+66S+7u7iVaFAAAAAAAAMo3FjoHAAAAAACA07HQOQAAAAAAAJyuwAudm0wmFjoHAAAAAABAkbDQOQAAAAAAAJzOoVDqckeOHLH9npaWJg8Pj2ItCAAAAAAAAOWfQwudX85qtWrixIkKCwuTj4+PDh8+LEkaM2aM5s+fX+wFAgAAAAAAoPwpcCj1r3/9Sx988IHeeOMNubu729rr1aun999/v1iLAwAAAAAAQPlU4FDqww8/1HvvvafevXvLxcXF1t6wYUPt27evWIsDAAAAAABA+VTgUOrEiROqWbNmrnar1arMzMxiKQoAAAAAAADlW4FDqbp16+rnn3/O1b58+XJFRUUVS1EAAAAAAAAo3wp8972xY8eqb9++OnHihKxWq7744gvt379fH374oVasWFESNQIAAAAAAKCcKfCZUl27dtXXX3+t77//Xt7e3ho7dqz27t2rr7/+Wu3bty+JGgEAAAAAAFDOFPhMKUlq3bq11qxZU9y1AAAAAAAA4AZRqFBKkrZt26a9e/dKkm699VbWkwIAAAAAAIDDChxKxcbG6qGHHtL69evl7+8vSYqPj9ddd92lJUuWKCgoqLhrBAAAAAAAQDlT4DWlnnnmGSUlJWnPnj06f/68zp8/r927dysxMVFDhw4tiRoBAAAAAABQzhT4TKlVq1bp+++/V506dWxtdevW1ezZs3XvvfcWa3EAAAAAAAAonwp8ppTVapWbm1uudjc3N1mt1mIpCgAAAAAAAOVbgUOpu+++W88++6xOnjxpaztx4oSGDx+ue+65p1iLAwAAAAAAQPlU4FBq1qxZSkxMVHh4uCIjIxUZGamIiAglJibq7bffLokaAQAAAAAAUM4UeE2pqlWr6vfff9f333+vffv2SZLq1Kmjdu3aFXtxAAAAAAAAKJ8KHEpJkslkUvv27dW+ffvirgcAAAAAAAA3AIcv31u3bp3q1q2rxMTEXM8lJCTo1ltv1c8//1ysxQEAAAAAAKB8cjiUmjFjhgYOHChfX99cz/n5+enJJ5/UtGnTirU4AAAAAAAAlE8Oh1I7d+5Uhw4drvr8vffeq23bthVLUQAAAAAAACjfHA6lzpw5Izc3t6s+7+rqqri4uGIpCgAAAAAAAOWbw6FUWFiYdu/efdXnd+3apSpVqhRLUQAAAAAAACjfHA6lOnXqpDFjxigtLS3XcxcvXtS4ceN03333FaqI2bNnKzw8XB4eHmrWrJk2b9581b5t27aVyWTK9dO5c+dC7RsAAAAAAADO5+pox5dffllffPGFbr75Zg0ZMkS1a9eWJO3bt0+zZ89Wdna2XnrppQIXsHTpUsXExGjOnDlq1qyZZsyYoejoaO3fv1/BwcG5+n/xxRfKyMiwPT537pwaNmyoBx98sMD7BgAAAAAAQOlwOJQKCQnRL7/8oqeeekqjR4+WYRiSJJPJpOjoaM2ePVshISEFLmDatGkaOHCg+vfvL0maM2eOVq5cqQULFmjUqFG5+gcEBNg9XrJkiby8vAilAAAAAAAAyhCHQylJql69ur755htduHBBBw8elGEYqlWrlipWrFionWdkZGjbtm0aPXq0rc1sNqtdu3batGmTQ9uYP3++HnroIXl7exeqBgAAAAAAADhfgUKpHBUrVtRtt91W5J2fPXtW2dnZuc6wCgkJ0b59+/Idv3nzZu3evVvz58+/ap/09HSlp6fbHicmJkqSrFarrFZrISu/PlitVhmGUeaPAwBKE3MpABQdcykAFF15mksdPYZChVLXi/nz56t+/fq6/fbbr9pn8uTJGj9+fK72uLi4PBdtL0usVqsSEhJkGIbMZofXrAcAXIa5FACKjrkUAIquPM2lSUlJDvUr1VCqUqVKcnFx0ZkzZ+zaz5w5o8qVK19zbEpKipYsWaIJEyZcs9/o0aMVExNje5yYmKiqVasqKChIvr6+hS/+OmC1WmUymRQUFFTmP7AAUFqYSwGg6JhLAaDoytNc6uHh4VC/Ug2l3N3d1aRJE61du1bdunWTdOlNWLt2rYYMGXLNsZ999pnS09P1yCOPXLOfxWKRxWLJ1W42m8v8myxdWmi+vBwLAJQW5lIAKDrmUgAouvIylzpaf6lfvhcTE6O+ffuqadOmuv322zVjxgylpKTY7sbXp08fhYWFafLkyXbj5s+fr27duikwMLA0ygYAAAAAAEARlHoo1bNnT8XFxWns2LE6ffq0GjVqpFWrVtkWPz9+/HiuhG3//v3asGGDvvvuu9IoGQAAAAAAAEVkMgzDKO0inCkxMVF+fn5KSEgoF2tKxcbGKjg4uMyf2gcApYW5FACKjrkUAIquPM2ljmYvZfsoAQAAAAAAUCYRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTEUoBAAAAAADA6QilAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA4HaEUAAAAAAAAnI5QCgAAAAAAAE5HKAUAAAAAAACnI5QCAAAAAACA0xFKAQAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOB2hFAAAAAAAAJyOUAoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNMRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTEUoBAAAAAADA6QilAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA4HaEUAAAAAAAAnI5QCgAAAAAAAE5HKAUAAAAAAACnI5QCAAAAAACA0xFKAQAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOB2hFAAAAAAAAJyOUAoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNMRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACc7roIpWbPnq3w8HB5eHioWbNm2rx58zX7x8fHa/DgwapSpYosFotuvvlmffPNN06qFgAAAAAAAEXlWtoFLF26VDExMZozZ46aNWumGTNmKDo6Wvv371dwcHCu/hkZGWrfvr2Cg4O1fPlyhYWF6dixY/L393d+8QAAAAAAACiUUg+lpk2bpoEDB6p///6SpDlz5mjlypVasGCBRo0alav/ggULdP78ef3yyy9yc3OTJIWHhzuzZAAAAAAAABRRqV6+l5GRoW3btqldu3a2NrPZrHbt2mnTpk15jvnPf/6j5s2ba/DgwQoJCVG9evU0adIkZWdnO6tsAAAAAAAAFFGpnil19uxZZWdnKyQkxK49JCRE+/bty3PM4cOHtW7dOvXu3VvffPONDh48qKefflqZmZkaN25crv7p6elKT0+3PU5MTJQkWa1WWa3WYjwa57NarTIMo8wfBwCUJuZSACg65lIAKLryNJc6egylfvleQVmtVgUHB+u9996Ti4uLmjRpohMnTujNN9/MM5SaPHmyxo8fn6s9Li5OaWlpzii5xFitViUkJMgwDJnN18Wa9QBQ5jCXAkDRMZcCQNGVp7k0KSnJoX6lGkpVqlRJLi4uOnPmjF37mTNnVLly5TzHVKlSRW5ubnJxcbG11alTR6dPn1ZGRobc3d3t+o8ePVoxMTG2x4mJiapataqCgoLk6+tbjEfjfFarVSaTSUFBQWX+AwsApYW5FACKjrkUAIquPM2lHh4eDvUr1VDK3d1dTZo00dq1a9WtWzdJl96EtWvXasiQIXmOadmypT755BNZrVbbm3TgwAFVqVIlVyAlSRaLRRaLJVe72Wwu82+yJJlMpnJzLABQWphLAaDomEsBoOjKy1zqaP2lfpQxMTGaN2+eFi1apL179+qpp55SSkqK7W58ffr00ejRo239n3rqKZ0/f17PPvusDhw4oJUrV2rSpEkaPHhwaR0CAAAAAAAACqjU15Tq2bOn4uLiNHbsWJ0+fVqNGjXSqlWrbIufHz9+3C5hq1q1qlavXq3hw4erQYMGCgsL07PPPqsXXnihtA4BAAAAAAAABWQyDMMo7SKcKTExUX5+fkpISCgXa0rFxsYqODi4zJ/aBwClhbkUAIqOuRQAiq48zaWOZi9l+ygBAAAAAABQJhFKAQAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOkIpAAAAAAAAOB2hFAAAAAAAAJyOUAoAAAAAAABORygFAAAAAAAApyOUAgAAAAAAgNMRSgEAAAAAAMDpCKUAAAAAAADgdIRSAAAAAAAAcDpCKQAAAAAAADgdoRQAAAAAAACcjlAKAAAAAAAATkcoBQAAAAAAAKcjlAIAAAAAAIDTuZZ2AQAAAADKhum7f9T0PT/l+Vx2drZcXFzyfG74rXdqeL02JVkaAKAMIpQCAAAA4JDEzDSdSE0o1DgAAK5EKAUAAADAIb5uHgrz8rNrM2ToZGqiJCnU01cmkynPcQAAXIlQCgAAAIBDhtdrk+syvJTMdPl+/JIk6c/7n1cFCwEUAMAxLHQOAAAAoNAMw7D9fi4txe4xAADXwplSAAAAAAosPv2iPjy4VW/9+bOtLfKLyYqsEKghdVqpT82m8rd4lmKFAIDrHaEUAAAAgAJZfWK/Hly3SKlZGbmeO5x0TjGbv9LLv3+rz+7uq+iw2qVQIQCgLODyPQAAAAAOW31iv7qseV8XszJlSLryYr2ctotZmeqy5n2tPrHf+UUCAMoEQikAAAAADolPv6gH1y2SYUjWXHGUPasMGYb04LpFik+/6KQKAQBlCaEUAAAAAId8eHCrUrMy8g2kclhlKDUrQx8d2lrClQEAyiJCKQAAAAD5MgxDs/ZuKNTYt//cwF35AAC5EEoBAAAAyNe59FQdSjrn4DlS/2NIOpR0TufTU0uiLABAGUYoBQAAACBfyZnpRRqfVMTxAIDyh1AKAAAAQL583CxFGl+hiOMBAOUPoRQAAACAfAVavBRZIVCmAo4zSYqsEKgAi1dJlAUAKMMIpQAAAADky2QyaUidVoUa+0zdVjKZChpnAQDKO0IpAAAAAA7pU7OpvFzdZXbwfCmzTPJyddejkU1LuDIAQFlEKAUAAADAIf4WT312d1+ZTMo3mDLLJJNJWn53X/lbPJ1UIQCgLCGUAgAAAOCw6LDa+rr9AHm6uskk5Yqmcto8Xd20ov0A3RtW2/lFAgDKBNfSLgAAAABA2RIdVlvHe4zRR4e2auaen3Uk+bztuRoVAvVM3VbqU7Op/Nw5QwoAcHWEUgAAAAAKzN/iqWfqtlb/mrfJb/HLkqTDD7yoahUqsqg5AMAhXL4HAAAAoNAuD6ACLF4EUgAAh3GmFAAAAACHTN/9o6bv+cmuzZBh+73uv9/IM5QafuudGl6vTYnXBwAoWwilAAAAADgkMTNNJ1ITrvr8yYuJVx0HAMCVCKUAAAAAOMTXzUNhXn55PpednS0XF5erjgMA4EqEUgAAAAAcMrxemzwvw7NarYqNjVVwcLDMZpatBQA4hv9iAAAAAAAAwOkIpQAAAAAAAOB0hFIAAAAAAABwOtaUKgOm/XhI0386nOdz1uxsma+yoOTwO2sopk1kSZYGAAAAAABQKIRSZUBiWpZOJFzrNrqZVx0HAAAAAABwPSKUKgN8PVwV5md/G13DMHQyMV2SFOprkclkynMcAAAAAADA9YjUogyIaROZ6zK85LRM+b68SpK0YXALVQ/wzjOYAgAAAAAAuB4RSpUx8RcztWjr33rr5yO2thqTf1BkoJeGtIpQ36ZV5e/pVooVAgAAAAAA5I9QqgxZvT9W3RdtVWpGdq7nDp9LVcxXe/Tyt/u0vG9TRdcOLoUKAQAAAAAAHGMu7QLgmNX7Y3Xf+5t1MTNbhiTjiudz2i5mZuu+9zdr9f5Y5xcJAAAAAADgIEKpMiD+Yqa6L9oqQ4asV6ZRV7AakiFD3RdtVfzFvO/KBwAAAAAAUNqui1Bq9uzZCg8Pl4eHh5o1a6bNmzdfte8HH3wgk8lk9+Ph4XHV/uXBoq1/KzUjO99AKofVkFIzsvXh1r9LtjAAAAAAAIBCKvVQaunSpYqJidG4ceP0+++/q2HDhoqOjlZs7NUvP/P19dWpU6dsP8eOHXNixc5lGIZmbTiSf8c8vL3hiAzDwSQLAAAAAADAiUo9lJo2bZoGDhyo/v37q27dupozZ468vLy0YMGCq44xmUyqXLmy7SckJMSJFTvXudQMHTqXmmsNqfwYkg6dS9X5VC7hAwAAAAAA159SDaUyMjK0bds2tWvXztZmNpvVrl07bdq06arjkpOTVb16dVWtWlVdu3bVnj17nFFuqUhOz32nvYJISs8qpkoAAAAAAACKj2tp7vzs2bPKzs7OdaZTSEiI9u3bl+eY2rVra8GCBWrQoIESEhI0ZcoUtWjRQnv27NFNN92Uq396errS09NtjxMTEyVJf0z9RT4e3teszyusgiIfaWDXdujjXUo9kZTvsYW0rKbgVlVtj7PTs/TnjN/yHSdJNR6tL+9Q30s1uBUtN0z86Zh27j2bbz+/WwJVrestdm373tmizKSMfMfe1KGmKjb833uYdjZFf83f4VB9tzzVVG6+Ftvjs1tO6tS6/C9X9KjkpVqPR9m1HV22R0lH4vMdW6lpqKrcE2HX9sfrGx2qN7zHraoQ4W97nHQkXkeXORaK1n+hpd3jU2uP6OzWk/mOqxDhr/Aet9q1/TV/u9LOpuY7tsrdEap0W6jtcWZiuva9u9Whems93kgelf73Hbmw84z+u+pgvuPcKrjrlqdvs2s7/tU+Jew7l+/YgAYhCutY065tz4xfZXUgnK3W9Rb53RJoe5xyMlGHP/oj33GSVHdYM7lY/jcdxm74W2c2Hs933PUwR0hSwr5zOv5V3nPm5cwWF9067A67thPfHtT5XWfyHVue5whDhi76ZOp08gGZZGKOYI7IhTnixp4jrsQckfccceVcyhzBHHE55gjmiCsxR+Q9R1w+l3qH+ZbpOcJqtTq0v1INpQqjefPmat68ue1xixYtVKdOHc2dO1cTJ07M1X/y5MkaP358rvZkc6oMs+ma+8rMUq61rRKykpRmzv8L6pZ6QYr93wRozchWkgPjJOns+XNKcU2TdGlNqWrebvo7JbNAl/CZJFXzt8ianaRkB/ZrTXeTxxXHmmikKMuc/+V/55LPKzP2f69levxFh4819lyc3NLcbY/jky84NDbdas313sRnJirVgbHmtHi5XDHW4fcm/qwuxv7vP54p8YmOH+sV+zyfFu/Q2OxMc+7PoTVZGea0fMdaki/IGvu/r3lmSobD9cadPyuLNcX2ODH5vENjXY3MXPVeSE9UiiP7TU+Q25WfQ6XIMOc/oZ1LPKf02P/9B+Xi+RTHjzUuTmZ3F9vj86mOfQ6vhzlCkpITHfssmZT7s3Q+PcGhseV5jjBkKNNsVZY5UyaZmCOYI3Jhjrix54grMUfkvc8r51LmCOaIyzFHMEdciTki/7k0K8tUpueIpKT8wzOplEOpSpUqycXFRWfO2KfrZ86cUeXKlR3ahpubm6KionTwYN6p6ujRoxUTE2N7nJiYqKpVq8rH6iUfq9c1t+3lWkHBwcF2bUmup+XmQOAX4FXRbmx2epbi8tlfjkoBgfIO/l8yOaheqF76reCLuQ+7M1KBSdkyWdPz7etn8c11rOdNx5Rpzf+vF4E+Aap42dg0c4rirf91qMbgwCC7v16Yj2Up3Rqf7zgPs1euelPd4uTiQBob4OGfa+wZR98b/0qqEOxve5yU4q5ka/5/HZKUa5/ZHimyWvP/cldwy/3eJJhPKM2a/1l0AT4VVemysZmJ6TpvdeyzFBRQye6vF26nDF20ns93nJvJPVe9aZbzMlvz/5+OAItfrrFx8pbVmv9fLwJ9A+UXfNlfL7ISlWjN/69DkhQUFGT31wt5pSvTmv8ker3MEZbzLkqxxuU7ziyXXPVmWhIla/7/01Ge5whDhi5aM+VpdZNJJuYI5ohcmCNu7DniSswRec8RV86lzBHMEZdjjmCOuBJzRN5zxOVzqXcZnyM8PDwc2p/JKOXbszVr1ky333673n77bUmXTvGqVq2ahgwZolGjRuU7Pjs7W7feeqs6deqkadOm5ds/MTFRfn5+SkhIkK+vb779rwfxFzNVdeIaXczMltWBd8tskjzdXPT3mPby93Qr+QIBoAyz/v9fQ4ODg2U2l/r9PwCgTGIuBYCiK09zqaPZS6kfZUxMjObNm6dFixZp7969euqpp5SSkqL+/ftLkvr06aPRo0fb+k+YMEHfffedDh8+rN9//12PPPKIjh07pgEDBpTWIZQ4f083Le/b9NJlJde+4lBmk2SSSZ/3bUogBQAAAAAArlulvqZUz549FRcXp7Fjx+r06dNq1KiRVq1aZVv8/Pjx43YJ4YULFzRw4ECdPn1aFStWVJMmTfTLL7+obt26pXUIThFdO1grBtyu7ou2KjXj0umFl580lZNVebq56PO+TXVv7eBc2wAAAAAAALhelPrle85WFi/fu1z8xUx9uPVvzfz5iI6c/981wpGBXnqmVYT6Nq0qP86QAgCHlafTpAGgtDCXAkDRlae51NHspdTPlELB+Hu6aWjrGnrstqryfXmVJOnIi3epWkVvmUz5XNsHAAAAAABwnSCUKgOm/XhI0386bNd2+QluLWf9kmcgNfzOGoppE1ni9QEAAAAAABQUoVQZkJiWpRMJV7/V6snE9KuOAwAAAAAAuB4RSpUBvh6uCvPzyPM5a3a2zC4uVx0HAAAAAABwPSK1KANi2kTmeRleeVoEDQBKWsLvM5Tw+8zcTxhStjVb/zW7/O9Wppfxa/ys/BoPK/H6AAAAgBsNoRQA4IZgTU9UdvKJqz6ffY1xAAAAAIofoRQA4IZgtvjKxSfMvtEwlJ1yUpLk4l1FMuU+69RsufotbAEAAAAUHqEUAOCG4Nd4WK7L8KyZKTo2u6IkKbTPbrlaKpRCZQAAAMCNiYWIAAAAAAAA4HSEUgAAAAAAAHA6QikAAAAAAAA4HaEUAOCGZRiG7ffsi+fsHgMAAAAoWYRSAIAbTnZavBK2v62THzextZ384Gb994M6Stj+trLT4kuvOAAAAOAGQSgFALihpB79Tn/Pj9D5H0cqK/Go3XNZCUd0/seR+nt+hFKPflc6BQIAAAA3CEIpAMANI/XodzrzVVcZmRclGf//c7lLbUbmRZ35qivBFAAAAFCCCKWA/2vvzuOqKtf//782g4AgoIIDSuCMiiMamOBxHlIkhyzPMW0gLRXSyswsT6ViltnPqQKttDRLS0zzaKZkilYqjmnJB02cUEMUUAZh7/39g99eiUN1Tskg7+c/5dp77X0vHvtxrXVf931ft4hUCOa8S5xf9wBYrYDlD95tAauV8+se0FI+EREREZHbREkpERGpEC7/9BHWghz+OCFlY8FakMPln5bezmaJiIiIiFRYSkqJiMgdz2q1krVvwf90bta++dqVT0RERETkNlBSSkRE7niWvAsUZh7jxhpSf8RKYeYxLHkZt6NZIiIiIiIVmpJSIiJyx7NcvfwXz8/+m1oiIiIiIiI2SkqJiMgdz66S2188v8rf1BIREREREbFRUkpERO54ds7VcfCoD5j+yzNNOHjUx8652u1oloiIiIhIhaaklIiI3PFMJhPurcf8T+e6tx6LyfTfJrNEREREROSPKCklIiIVglvThzA5VubP3/rsMDlWxq3psNvZLBERERGRCktJKRERqRDsnT2p0fdTMJn449ufHZhM1Oi3AntnzxJonYiIiIhIxaOklIiIVBiV/XtSM+ILTI4uFNWXun5ZXtExk6MLNe9bQ2W/HiXfSBERERGRCsKhtBsgIiJSkir798T3sV+4/NNSsvbOozDrF+M1B496uLceS5VmD2Hn5FGKrRQRERERufMpKSUiIhWOvbMnHm3G4tb8YU68XbSzns8j/0cl97tU1FxEREREpIRo+Z6IiFRY1yag7J2rKSElIiIiIlKClJQSEREREREREZESp6SUiIiIiIiIiIiUONWUEhGRCiFzz/9H5p45xQ9arcb/nvkwEEw3jtV4tH0Kj7bjbnPrREREREQqHiWlRESkQrDkZ2G+fPqWr5uvpN3yPBERERER+fspKSUiIhWCnZM79m51bnzBCmaLGXs7e7hJnXM7J/fb3zgRERERkQpISSkREakQPNqOu+kyPIvFwvnz56lRowZ2diq1KCIiIiJSUvT0LSIiIiIiIiIiJU5JKRERERERERERKXFKSomIiIiIiIiISIlTUkpEREREREREREqcklIiIiIiIiIiIlLilJQSEREREREREZESp6SUiIiIiIiIiIiUOCWlRERERERERESkxCkpJSIiIiIiIiIiJU5JKRERERERERERKXFKSomIiIiIiIiISIlTUkpEREREREREREqcklIiIiIiIiIiIlLilJQSEREREREREZESp6SUiIiIiIiIiIiUOCWlRERERERERESkxCkpJSIiIiIiIiIiJU5JKRERERERERERKXFKSomIiIiIiIiISIlzKO0GlDSr1QpAVlZWKbfkr7NYLGRnZ+Ps7IydnfKLIiL/C8VSEZG/TrFUROSvu5NiqS3nYsvB3EqFS0plZ2cD4OvrW8otERERERERERG5c2VnZ+Ph4XHL103WP0pb3WEsFgtnzpyhSpUqmEym0m7OX5KVlYWvry8nT57E3d29tJsjIlIuKZaKiPx1iqUiIn/dnRRLrVYr2dnZ+Pj4/O6srwo3U8rOzo66deuWdjP+Vu7u7uX+BysiUtoUS0VE/jrFUhGRv+5OiaW/N0PKpnwvUhQRERERERERkXJJSSkRERERERERESlxSkqVY05OTvz73//GycmptJsiIlJuKZaKiPx1iqUiIn9dRYylFa7QuYiIiIiIiIiIlD7NlBIRERERERERkRKnpJSIiIiIiIiIiJQ4JaVEREREROR3WSyW0m6CiIjcgZSUKgdmzZpFeHh4aTdDRKRce//9941YqnKKIiL/HTu737oNiqEiIvJ3UVKqHPDy8uKbb77h119/Le2miIiUW/b29qxbt478/HxMJlNpN0dEpFxZs2YNTzzxBFarVTFURET+NkpKlTEWi+WG6dE9evTg6tWrHDx4sJRaJSJSvlit1htiaYcOHXBxcSExMbGUWiUiUn6YzWbMZrPxb1dXV+Li4khPTycxMZHY2NhSbJ2ISPl34sQJ5syZw5EjR4CKOwtVSakywtZ5srOzKzY9GqBOnTo0atSITZs2lUbTRETKDVssNZlMN8TS2rVr07x5c9asWQNU3Bu/iMifYW9vj729vfHvb7/9FpPJRN26dQkPD2ffvn1cvny5FFsoIlI+2J45jx07xqlTp4zj77//PnPnzqVevXpYLJYKOwtVSalScu3IE/y2Tv/bb79l1qxZbN26lcLCQuP1nj178vXXXxc7JiJSkd0sqWSLpYcOHeKDDz5gx44dxmuVK1emS5cubNy4scTaKCJSVt0shtqOZWZmsnjxYoYOHcq4ceOAoqS/u7s7Dz74IBcvXuSdd97Bzc2tJJssIlKuXLlyBSgaLD148CB9+/YlPDycU6dOUVhYyOeff85TTz1FpUqVbhhMrUgq7pWXMtvI09mzZwF47733aNiwIcOGDWP16tWMHDmSyMhI4/19+vTh4MGDnDlzplTaKyJSVtg6TdePJhUUFLBixQoCAwPp1KkTb7/9No899hjR0dFAUdzt0qULycnJXLhwocKORolIxWWxWIyB0ZvFQJPJxE8//UT//v1544038PLyonbt2uzevZtp06bx1FNPsW/fPvLz80u66SIi5UJKSgpPPPEEd911F127dmXatGmYzWZatGhBUlISNWrUICwsjDfffJP09HRtaIaSUreF1WqlsLDwplvnZmVlsXv3bjZu3IijoyNDhw4FirKoU6dO5eTJkyQmJrJq1So+/PBDvv/+ewCCg4Oxs7Nj165dJXotIiKlyWq1FptZaiuwm5eXR0JCAklJScXe/8svv/D0009z6tQpdu3aRWxsLLGxsezduxeApk2b4u7uzubNm0v0OkRESoPFYik2I8rOzs4YGN27dy+JiYnk5OQYr+fk5BATE0NeXh7btm1j3rx5PPPMM7Rr1w6Azp0789NPP3H69OmSvRARkTLqwIEDPPfcc6xbtw6r1cqUKVM4evQos2bNYvjw4bz66qtMmjSJX3/9lcqVK/P5558zYMAApk+fTkhICFWrVgUqdlkJJaVuA5PJhIODA3Z2dsU6Uzk5OTz33HP06dOHRYsWsXr1aj7++GMABg8ezNChQzl16hRz5sxh0qRJAKxfv56cnBw8PDxo06YNX331Valck4hISbg+mW8ymYrNLDWZTHz99dc0bdqUESNGMGjQID777DMKCwtxdHRk8ODBPProoxQUFPDpp5/y3nvvUVBQwKpVqwCoWbMm7dq1Y/Xq1UDFfgAQkTvP9THNzs6u2IyokydPEhUVRc2aNenRowfR0dEMHDiQ9PR0AC5dusT69et59tlnqVatGgAODg7G+UFBQTg7OxuDpoqhIlIR5eXl0b9/fz7++GNeeeUV9u/fj6enJ+vXr2fLli2MHTuWIUOGMGbMGOLi4ti0aZNRPsLNzY1hw4Zx+fJlDh48yPTp04Gbz16tKJSU+h9dvyPJtVJTU5k8eTJt2rShZ8+eLFq0iPz8fCpXrkyrVq0oKCigYcOG9O3bl9q1a2O1WvHx8WHJkiX07NmTFStWEBYWxoABA1i7di3Z2dkA9OrVi61bt5Kbm1uSlyoicltdG0uvXU9/4sQJLly4wGOPPYabmxv/+Mc/WLRoEUuWLOGdd94hNTWVtm3b8vbbb3PgwAEAGjRowOrVqwkJCWH69Ol4eHgwcOBAPv/8c65evYqzszPdu3fn22+/BSr2A4CI3Hmuj2m//voro0ePZvv27QBkZGSQnp7O8uXLOXfuHKtWrSIjI4OZM2dy9epVPDw8yMrKwtXVFbhxoKBKlSqEhITw5Zdf3vT7RETuNBcuXGDlypXs3LmTgoICAJydnTl16hTDhg2jffv2fPXVV3Ts2JGff/6ZatWqcc899xjnd+rUifr167Nt2zbjWFJSEvXq1SM2NpZFixbx8MMP33SVVUWhpNR/4dqO0/U7kkDRaJHVauWll15ix44djB49mgEDBvDSSy8ZGdC2bdvi5OSEj48PAIWFhZhMJn788UemTZvG/fffz9dff82zzz7Lfffdx9GjRzl58iQA9957L8nJyaSkpJTQFYuI/L1scfJatliam5vLsmXLeOGFF/jiiy9o164dDz/8MHXr1iUxMZFOnToRHR2Ns7MzvXv3xs7OjmeeeYaCggK2bNkCFK3jnzBhAgMGDGDr1q3MnTuXQYMGcezYMX755Reg6OEgLS2N5OTkEr12EZG/w7V1oa63fft2du7cafw7ISGB5cuX06RJEwD8/f2ZNWsWXbt25ezZs+zZs4esrCw2bdpEcnIyrq6uBAYGGiP6dnZ2Rsy27bT3z3/+k7Vr1/Lcc88xfPhwnn766dt5uSIiJW737t2sW7eOBQsW0KJFC8aPH88///lPJk2aZMTfgQMHUqVKlWI1oVq2bMnx48fJzc01Ymf9+vVxcnLi0qVLxrnz5s1j1KhRdOvWjQ0bNrB9+3aCg4M5fPhwyV9sGaCk1C1YrdYbspW2jlN+fj4rV65kyJAhTJw4kWPHjgFFo0XvvPMOKSkprFu3jscff5yxY8cycuRIFi5cyKZNm2jUqBHt2rUr9sAARTf9X375hcGDB1O5cmUyMjJYvHgx2dnZJCUlYbFYaNmyJQEBAcYUaxGR8uLawrrXj6x/9NFHTJo0ienTpzNv3jxcXV2pXr06HTt2ZNeuXQwaNIjWrVvz0ksv0aFDB3799Vfj3ICAALy9vY2aUQBHjx7l4YcfxtPTk5ycHD766COuXr1q7MRXr1496tevz5EjR0rgykVE/l7X1oU6fvw4eXl5AEYtk379+rFmzRoA5s6dy6hRo/Dy8sJqteLh4YGrqysjR46kXbt2xMTEEBAQgNlsNmLksGHDWLZsGV988YVRx++HH35g6dKlQFFSKiYmhl27duHs7MwDDzxQCn8FEZG/T1paWrFd7uPj4wkPD2ft2rXEx8dz9OhRnnnmGVauXMnrr78OQFhYWLGafAChoaEUFhaSkJBQ7Hk3KSmJJk2aYG9vzzfffEN6ejodO3YEimpHf/nllyxcuJBmzZpVyGXRDn/8lorFYrEY6++v7zh9+eWXJCQk4Obmxt69e6lduzYbNmwgMTGRxYsX06hRI44dO4aDgwMJCQnMmTOHAwcOcPXqVfr27UvNmjVxd3enRYsWbNq0CfhtnX6zZs3w8vJi7Nix9OnTh127djF06FDy8/PJzMzEYrFQqVKlCps9FZHyzdaBSkpK4uzZs4SEhFC9enUArl69ysyZM2nXrh0fffQRTZo0ITc3lwYNGrBr1y5atmwJgI+PDyEhIXz11VcUFhbi4OBA9erVjfedOXOGBg0a4OPjw+jRoxk8eDBJSUkEBweTk5NjJKFq1aqlGaciUmbs2bMHs9lM+/btjSTQ70lJSWHatGmsWbMGT09PWrduTVRUFF26dGHz5s08++yzTJgwgcWLF2OxWOjevTvw20YR7733Htu2bWPVqlV06NCB7OxsAgMDjWXQo0aN4vDhwzz22GN0796dlJQUzp49y5gxY7h69SqVKlUiKiqKqKio2/63ERH5O9niYHp6OhcvXuStt97igw8+oEGDBvTq1YuXX36ZKlWqMGLECGbMmEHt2rUJDg4G4Mknn+Tw4cP85z//ISoqipCQEJydndm/fz/NmzfHYrHg7OzMgw8+yKxZs7hy5QrDhw9n+fLlAEYS6siRIzz88MMEBQUZ7bLNZoWKuSy6Qs+UMpvNNy0ICUVT9hYuXFhs9D0/P5/FixcTHx/P1KlTiYuLY+HChWRnZxMfHw+Al5cX27dv5+WXXyYoKIgVK1aQmprK0qVLadGiBQ4ODjRr1ozc3Fz27NkDFHXIAFavXk3jxo355JNP8PX1JSIigsTERJ577rliRSZvNWVbRKS03Gx26bVWr15NnTp16NGjB5MnT6Zfv37GqHyfPn2oXr06zZs3N27KLi4uBAcHc/HiRc6dOwcUJfGbNm3K5cuXi+2617p1a65evcr333+PyWRiyZIl1KhRgxkzZpCbm8vQoUNJSEjgtddeK9YmxVIRKW1Xrlxh2LBhLFy4EPj9zogtxi5YsIBz587x5ZdfsmXLFjw9PZkxYwZbt24FICYmhhkzZrB582bOnTtHixYtgKJn3JycHLZv387dd99Nhw4dANiyZQtZWVn8+OOPXLhwAVdXV+Li4liyZAl16tRh+PDh7Nixg0mTJlGpUiWjPb+327SISFlkMpk4dOgQjRo14oUXXqCwsJCNGzfy7LPPsnjxYuNY48aNqVKlCj4+PsVmUN19991cuXKF/fv34+TkRLt27YyNyGzxOyYmhr59+/L222/j7+/PpEmTGDlypDFA8MQTTxATE4Ozs3OxtlXEGVI2FWqmlG1rcTs7u2JTn6+1b98+Hn74YU6fPk3jxo2ZP38+nTp1Yt68eXTq1Ak/Pz/q1KlD69atgaIfZtOmTdm9ezdXr17F19cXb29vPvzwQ5o1a2Z87okTJ0hNTSUsLIwmTZpQUFDA2rVradu2rdGOkJAQgoODb3ggMZvNxdp6s3aLiJQm2+zSzMxMTpw4QfPmzY0k/5kzZ3j++eeJjIxkwoQJHD9+nAkTJvDkk0+SmJiIj48Pfn5+ODo6kpWVhbu7OwCNGjXCy8uLDRs2MGLECKBoJKlGjRps2bLFGLkKCAggLy/PqBHVrVs3wsLCinWe4LeZsDaKpSJSErKysli6dClLly5l3rx5BAUFGRvmuLq6cs8993DhwgVOnTpF3bp1bzlbys7Ojm3btpGYmMj69evx8vIiLy+Pli1bsmLFCuLj4+nUqROVKlWif//+2Nvbk5eXx+OPP857772Ht7c3lStXpkaNGmzevJmlS5fi5OTE2rVr6dq1KxcvXiQ1NZXq1atjZ2dH37596du37y2vy7bbtIhIWZSQkEB8fDy5ubkMGzaMzp07A9C8eXPq1q3Lf/7zHzZt2kSHDh0ICwvDarUybdo0IiIi6N69O8HBwezfv5/MzExjdn+tWrUoKCgwJpWEh4czb948Ll++jJubG1A0s3/mzJkMHToUZ2dnAgICbmibLSdxbayviDOkbCrUTCnbzdPOzo78/HxOnTrFs88+a+xud+XKFZ5++mm6d+9OWloa27dv58MPPyQ2NpYtW7bg7e1N48aNcXJy4vz588bntmjRgjNnzpCcnEx4eDh+fn5ER0eTkJBAbm4uBw4cYPr06axevRqLxUL9+vUZO3YsXbp0AYp3jEwmExaLhcLCQiNbqo6TiJRlO3fuJC0tjYiICGrXrk14eDhjxowxZiIdOHCAjIwMhg8fjpubG4GBgbz22mucPXvWmGXavn17UlJSitWL8vHxoVWrVmzYsME45ufnh7u7u7HzExTNlFq3bh3PP/+8caxSpUo3jOJfm5ASESkpU6dOZcGCBYSGhuLn5wcUPdvZEudt2rTh119/5aeffgJ+f7T86tWrJCUlERMTQ8OGDfH29mbhwoWMHTvWSN4DfPbZZ9SqVYuPPvqI9PR0IiIijJ2fJk2aRO/evXn++ed5+umnCQ4OZsGCBWzevJm2bdsW+z7bgG5FHsEXkfLl/PnzPPjgg4waNYqzZ89iNpvp1asXixcvJj8/H4CgoCD8/PyoWbOmcV5oaCi+vr7Gc+egQYPYvHlzsV3z9uzZw+nTpwkJCQGgc+fOpKamkpaWVqwNVquV1q1bGwmp6+Oovb19hU5CXa9CPaHn5eUZN/FevXoxf/58Zs+ebdQZOX/+POnp6Tz//PNcvnyZd999l1mzZlFYWGhsHx4UFMT58+c5evSo8bnt27fHbDazc+dO3N3defvtt3FwcGDcuHEEBAQQEhLC2bNn6du3LyaTCW9vb6Kjo+nUqdNN22lnZ4eDg4N+qCJS5p07d46QkBAef/xxI7H0yiuvsHjxYhYtWgQU3cBtM0Rt6tevT7t27YzkUs+ePTl+/LixQx5A1apVCQoKYvXq1caxGjVqMGHCBObPn28cs7e3p3bt2je07dqBCBGR0pCZmUliYiIRERG8/vrreHl5AfDVV1/RvXt3IiIiOHbsGNnZ2UZS6vdilouLC1WrViU5OZnXXnuNgwcPcujQIV577TVat25tDAZ8+umntG3blm7duvHpp5/SpEkTBg0axJdffomfnx9vvPEGP//8MydPnmTMmDHUqlULk8l0Q/LJZDKp8yQi5UpBQQFdu3Zl8+bNrFy5kg8++IBJkyaxePFio3xOaGgoZrOZU6dOGefVqlWLWrVqcfr0aQAiIiLIz88nKiqKt956ixdeeIF3332XyZMn4+LiAkDbtm25fPkyjRo1KtYGW8y8dpKJ4uitVagn9TVr1vD+++8zduxY3nzzTY4fPw5g1DX54YcfyMnJITAwED8/P2JjY6lRowYbNmxg/PjxAHTo0IErV64UKzjepk0b7O3t2b59O1CUuNqwYQPz58/nk08+IScnhy+++IKuXbsW+zFqDb6IlCVms7lYXPozI+M1a9akR48e/PDDD/Tt2xcfHx9GjBjBQw89xJIlS4CiadJZWVmkpqYa57m6ulK3bl1ji/GuXbtitVqLxVZHR0fCwsJ4/PHHuXLlinH8H//4B61atfrL1ysicrt5eHhQv359EhISeOihh3jwwQfJyMggKioKf39/HnvsMU6dOsWhQ4c4fPhwseT9zdSsWZOqVavSvn17Bg8ejL+/PwCnTp1i2bJlHD9+nL1795KYmMjQoUMB8PX1JTY2lqSkJPr16weAs7Mzbm5uN8zOV6dJRMqqvLw8srKygN9/RvX29iYyMpIqVaowZ84c+vTpw8yZMzly5Aj79u0Diko9ZGZmkpCQYJxntVrZsWMHHTt2xGKxUKtWLe666y569erFiRMn2LdvHxMnTuTJJ58s9n2VK1e+ZXsUU/+ccrkQ3Gw2YzKZjJGkP9qpxFZHZP78+QQFBTFu3DgA4uLiyM7OZu3atURHR+Pr64ujoyP9+/fnrbfewsXFxfiO3NxcoCgbWqVKFX788Ufj82vUqEGHDh3w9/enoKAAR0dHrFZrsZlQtpGra5fiafReREqLrTD5zerV5eTkkJ6ezl133fWnPqtdu3akpqbSuHFj49iAAQNYunQpJ0+eJCwsDEdHR1atWkWvXr0AOHv2LKtXr+bFF1/EbDbj4eFB1apVSUlJ4cqVK7i6ugJFyaquXbvetP260YtIWRcdHc3y5cuxt7encePGjBo1ihkzZuDi4sKrr76Kj48P/fv3Jzw8nMOHD5OamkrDhg1vGeMaNGhAVFQU06dP5/z580RGRpKSksJnn32G1WqlY8eO2NvbM3r0aLp162ac5+joiK+v7w2fZ6uzKiJSFm3ZsoWdO3fSqFEjXn75Ze69915mzJiB2Wy+ZU27SpUqkZWVxaOPPsqFCxfo3LkzkydPJioqisOHD1NYWEi9evXw9/fnzTffpGrVqvTu3ZsVK1bg7OxMaGioERc7duxIamoqq1atokqVKrdsp55J/5oyfxeyrWW/lr29vbGDyIkTJ/7wR2BnZ8f58+eL3ZAtFgvu7u4MGTKEXbt2kZmZSWBgIPXr1+f48eO4uroaP8Zjx44RGRlJWloarq6ueHt7c/78eS5cuGB8xxtvvMGYMWNwdHQEbj5lT7WhRKS02eKpbUmGjdVq5eOPP6Z58+bUqVOHf/3rX0yZMsWIc783ItWzZ09++eWXYrX2OnXqRGFhId9//z3VqlVj4sSJfPzxx/zrX/8iNjaWJ598kmbNmjFgwACjHV999RVz5841ElI2FovlhpmluvmLSHkQHR3N4sWLqVOnDqNHj6ZLly7s27ePoKAgfHx8jJh8//33k5WV9afqSo0ZM4aYmBiOHTvGgAEDGD9+PF5eXkycOBF/f39atmzJK6+8UmxnJ8VMESnLbDHv9OnTRimHxx57jIiICJKSkvj88885ePCgMaP+j5LpcXFxHDhwgNmzZ/Pqq68SGhqKvb29kfyHoiV8bm5u7Nu3jwEDBhAXF8cLL7xAixYtjOfOBx54gD179hjL/P5oNqv8b8rsTCnbjnM36zgtX76c6dOnc+bMGQIDA+nSpQtPPfUU1atXv+XIUqVKlfDx8TGKkNl+yIGBgWRnZ7N79266devGlClTGDBgAMHBwdx7770kJyezd+9eGjVqZCwfWbhwIVWrVr1lm6+lhwARKS22ZM6169htMSorK4sFCxZw9OhR/v3vf5ORkcHs2bOJjIzkvvvuY+/evcyfP5+LFy8yb9683/2etm3bUqlSJb7//nvq1auH1WrF1dWVtm3bsmXLFgYOHMj999+Pi4sLn332GXPmzCE0NJRp06ZRt25d43O8vb1v+vkaxReR8qphw4aYzWZ8fHyIj48nJCSEgIAAdu7cCWDE6LCwMMaPH8/PP/9MeHj478Y9e3t7IiMjjaS+p6fnDe+xWCzGrqgiImVVRkYGZrMZb29v/u///o/hw4czePBgwsLC+Oijj9i0aROdOnXi5MmT5Obm8s0331BYWHjLWVK2FVInTpygUaNGeHh4APDJJ5+QlpZGZmYmu3btokGDBoSEhLBhwwbuv/9+4uLijMkl1+rRowcXL17k0KFDNG3a9Kbvkb+uTDzpX7+eHYp3nGbMmEFkZCQnT540Mp6RkZHs2bOH8ePHk5iYyMsvv/y73+Hp6UlAQAD79u0rNpr/3XffYbVa2bVrFwDBwcFs3LiRoUOHsn37dipXrsycOXNYvXq1MZ3alpC6fhRLM6FEpCy5dtOES5cukZGRwYYNG2jfvj3Lli1j06ZNtGrVCqvVyieffMLdd9/N+PHjqVevHg0aNMDJyYkFCxZw6dKl3+3YVKlSheDgYDZt2lTseIcOHVizZo1RN6pfv34sWrSIw4cPExcXR4sWLW7r9YuIlAW1atWiZcuWRu3RXr16kZSURGpqqtHByc3NJT8/n0OHDhm7Ql/LYrFgNpuNEhZms5nq1avj6el5w06jwA1bjYuIlAXXroD6+uuv8fLy4rXXXgPgrrvuIi0tjcDAQBISEmjbtq2xy52vry+RkZGYTCa2bNkC3Lw+s61/fu+993Ly5EmGDx9OeHg4cXFxTJkyhcGDB9OkSRMAunfvzpUrVzhx4oRRfqewsND4LIvFgouLC1u3buW+++67HX8O+f+ViaTU7e442X6c9913H46OjgwdOpTvv/+e3bt3c+jQIYKCgowdoCwWC4GBgYwbN46NGzeycOFCevToYexIcu3n62YvImXZ8ePHiY6Opm7duoSEhPDBBx9Qp04dkpKSeOutt3jnnXeIiorirrvuIj4+nosXL3LvvfdSvXp1evToQaVKlZg/f/4tR6Ou1blzZ5YvX87ly5eN2Dh58mQSExONUSoABweHm3agRETuVB4eHrRu3ZrU1FTS09Pp168fDRo04IknnmDHjh1cunSJBQsW4O3tzY8//khycjJQtEzE1oGzs7MzSkEUFBRgb29vxFDtNCoi5cH69evx9vY2Vh9lZmbi7+/PsmXLSE5OxsnJiZycHEwmEzk5OXh6ehbbHa9hw4b4+vqyfv164OZJKdskkV69evHpp5/SsmVLmjZtyhtvvGHU9GvTpg1WqxVPT088PDxISEggIyPDiKU2dnZ2WK1WQkND/9SzsPzvysTd63Z3nGwdpObNm7NgwQKcnJzo3bs3Xbp0wdfXl3HjxnHw4EGg+DKR6ztOSkKJSHlhNpuZMmUKR44cYfbs2axYsYJWrVrh5+eHv78/Xbp0KVaYvGXLlsTHxxMUFMTatWtJSUnhiy++YPTo0bi5uf3h9z344IO8+OKLxeKwl5cXfn5+N7xXHSgRqWiaNGmCq6sr33zzDQDvvvsujo6ODBkyhOrVq1NQUMB7773H22+/TZs2bYCi4uS2DtZ3333HxIkTad68ubHsWTFURMqT7t274+bmxqxZswA4d+4cDzzwAPXq1SMuLo60tDRatGjBiRMnCA0N5cyZM6SkpBjnV61alcLCQmPW6a1WKdmS+YGBgcTGxvL6668TFBQE/Fbawvae6OhoBg0adMsi5ur/l4xST/nZOk7nzp1j9uzZBAQEcP78+T/sOE2YMIEXX3yRwMBA3N3d//T3dezYkbVr13L06FHjcx955BHq169PZmZmsRH967OlIiLlxbJly9iwYQNLliyhT58+xWZ6tmnThpMnT2I2m40lHiEhIezYsYPHH3+82K573333HcnJyYwYMeKG7zCbzVitVhwcHGjYsCEvvPBCiV2fiEh54u/vT7Vq1Vi5ciX3338/Xbp0ISgoiF27dtGkSZNi9fWgaDnf/PnzWbt2LQcPHsTFxYW7776b6OhowsLCSukqRET+d46OjowbN45Vq1bRr18/vL292bhxI7NmzeKNN95gxYoV1KhRg5ycHNq3b4+bmxsffPABoaGhVK5cme3bt3Pp0iWOHz9ORkYG1apVw2q1FttR+vp61ACFhYXGTqO2ZL7tv48++miJ/x3kRqWecSnpjhNAamoqzs7OHD9+nC+++ILvvvuOl19+GQ8PD20zLiJ3hL1791K/fn369OkDFCXZbcUf+/Tpw+TJk8nKyjJq5D3yyCMsW7aMfv368corr+Dj48PGjRvZunUrERERFBQU4ODgcMtNKC5fvoybm9vvFp8UEamoatWqxYgRI6hcubJxzN3dnW7dugG/lZqwWq3GrtH79++nc+fOvPnmmwQGBuLi4lIqbRcR+bsMGTKEvXv3MmPGDKZOncqPP/5Ix44d2bZtGxs3buTYsWP06dOHqlWrEhMTw8CBA4mIiMDR0ZFLly4xYcIE4uLiSExMpH///saGDrYk09mzZ4mPj+enn35i5MiRBAYG6rm0HCj1eb+36jgB9OnTh6SkJLKysoxE0SOPPELNmjXp168f8fHx/PDDD0ydOpUpU6aQnZ1NQUGBsezOdoO3t7cv9mO8cOECTzzxBC1atCA2NpaoqCjCw8ON7xcRKe/8/f2NdfjX1iQB6Nq1K+np6cbW41A0Jfrjjz8mJCSEWbNm0atXL2NHkmHDhuHo6GjMHjWZTKSlpbFo0SL69u1LnTp1mDt3LoBu/CIiN+Hi4sKoUaN46KGHih23Pate37Hy8/Nj6dKlvPrqq7Rv314JKRG5I9StW5eRI0eydu1aPvvsM5o1a0ZGRgaRkZFkZ2dz5MgR472dO3dm8+bN3H333fj7+/P+++8TFBREYWEh+fn5AOTk5LBmzRoeeeQR6tWrR0BAAO+++y5VqlShZs2apXWZ8l8q9d6Dv78/K1euBDBG4G/WcbrnnnuA3zpOb775JrNmzeLQoUM0b96cESNGMGTIEGMXE1vHKC0tjXXr1hEfH8/evXsZP348EyZMICYmhgYNGtxy/aiISHnWpUsXxo8fz7Zt24ot9Th27BjVq1cnICCA9evXG7HVarXSpEkTYmNjOXPmDHXq1LnhM3fv3s3s2bP54YcfSE9Pp2HDhnTt2pWnn36aDh06lNi1iYiUV7ZnXRsNhopIRWK1WgkLC6N///7ExMTQr18/8vLy8PHxoX///iQmJpKRkWG8v02bNrRu3dqIlRMnTsRsNhMaGkpOTg6tWrUCIDQ0lJkzZxIWFkbt2rVL5drkf2ey2oZoSsmBAwdo3bo133777U07Tvfccw8DBw5k6tSpAMbyOqvV+l93nHr37k2HDh2KTZ0WEblT9ezZk4yMDJ588kkiIiLYv38/q1atIjIykri4ONLT041BgZu5dumzxWJh7ty5HDhwgN69e+umLyIiIiL/FVtffsuWLXTv3p3GjRtz8OBB7O3tuXz5Mnl5eXh5eRnvKygoYNKkSWRnZ/Ptt99SWFjIzJkzGTRoEACnT5+mdu3a2vihnCv1pBSo4yQicjukpKTw1ltvsXXrVtLS0igsLOShhx7i+eefp0aNGsbMUhERERGRkpKfn8/WrVvx9/enUaNGxV67dhKKyWTi9ddfJzk5mbCwMCIiIvD09CydRsttUyaSUuo4iYjcPnv37sXZ2ZmmTZsWO66NHUREREREpDSViaSUjTpOIiJ/r+vj57W7mYqIiIiIlBb18wXKUFJKHScRkdtHN30RERERESlrykxSykYdJxERERERERGRO1+ZK1OvhJSIiIiIiIiIyJ2vzCWlRERERERERETkzqeklIiIiIiIiIiIlDglpUREREREREREpMQpKSUiIiIiIiIiIiVOSSkRERERERERESlxSkqJiIiIiIiIiEiJU1JKRERERERERERKnJJSIiIiIiIiIiJS4pSUEhERERERERGREqeklIiIiIiIiIiIlLj/B34oC0LyDUxEAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "palette = sns.color_palette(\"colorblind\")\n", "\n", "ci_cre_general = plpr_tune_cre_general.confint()\n", "ci_cre_normal = plpr_tune_cre_normal.confint()\n", "ci_fd = plpr_tune_fd.confint()\n", "ci_wg = plpr_tune_wg.confint()\n", "\n", "comparison_data = {\n", " \"Model\": [\"cre_general\", \"cre_normal\", \"fd_exact\", \"wg_approx\"],\n", " \"theta\": [plpr_tune_cre_general.coef[0], plpr_tune_cre_normal.coef[0], plpr_tune_fd.coef[0], plpr_tune_wg.coef[0]],\n", " \"se\": [plpr_tune_cre_general.se[0], plpr_tune_cre_normal.se[0], plpr_tune_fd.se[0], plpr_tune_wg.se[0]],\n", " \"ci_lower\": [ci_cre_general.iloc[0, 0], ci_cre_normal.iloc[0, 0], ci_fd.iloc[0, 0], ci_wg.iloc[0, 0]],\n", " \"ci_upper\": [ci_cre_general.iloc[0, 1], ci_cre_normal.iloc[0, 1], ci_fd.iloc[0, 1], ci_wg.iloc[0, 1]]\n", "}\n", "df_comparison = pd.DataFrame(comparison_data)\n", "\n", "print(f\"True treatment effect: {theta}\\n\")\n", "print(df_comparison.to_string(index=False))\n", "\n", "# Create comparison plot \n", "plt.figure(figsize=(12, 6))\n", "\n", "for i in range(len(df_comparison)):\n", " plt.errorbar(i, df_comparison.loc[i, \"theta\"],\n", " yerr=[[df_comparison.loc[i, \"theta\"] - df_comparison.loc[i, \"ci_lower\"]],\n", " [df_comparison.loc[i, \"ci_upper\"] - df_comparison.loc[i, \"theta\"]]],\n", " fmt='o', capsize=5, capthick=2, ecolor=palette[i], color=palette[i],\n", " label=df_comparison.loc[i, \"Model\"], markersize=10, zorder=2)\n", "plt.axhline(y=theta, color=palette[4], linestyle='--',\n", " linewidth=2, label=\"True effect\", zorder=1)\n", "\n", "plt.title(\"Comparison across DoubleMLPLPR approaches\")\n", "plt.ylabel(\"Coefficient Value\")\n", "plt.xticks(range(4), df_comparison[\"Model\"], rotation=15, ha=\"right\")\n", "plt.legend()\n", "plt.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We again see that the `wg_approx` leads to a biased estimate in the non-linear and discontinuous `dgp3` setting. The approaches `cre_general`, `cre_normal`, `fd_exact`, in combination with [LightGBM](https://lightgbm.readthedocs.io/en/stable/) regressors, tuned using the [Optuna](https://optuna.org/) package, lead to estimate close to the true treatment effect.\n", "\n", "This is line with the simulation results in [Clarke and Polselli (2025)](https://doi.org/10.1093/ectj/utaf011), albeit only for only one dataset in this example." ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.12.3" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }