{ "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_demean`, `d_demean`, `xi_demean` 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": 17, "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_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_jobs=5))])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ml_lasso = make_pipeline(\n", " preprocessor, StandardScaler(), LassoCV(cv=2, n_jobs=5)\n", ")\n", "\n", "ml_lasso" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_diff 0.517696 0.018044 28.690074 5.074285e-181 0.482329 0.553062\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": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1050" ] }, "execution_count": 19, "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": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.560037 0.028292 19.794714 3.306228e-87 0.504585 0.615488\n" ] }, { "data": { "text/plain": [ "1051" ] }, "execution_count": 20, "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": 21, "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(cv=2, n_jobs=5)\n", ")" ] }, { "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_demean 1.151822 0.014172 81.277037 0.0 1.124047 1.179598\n" ] }, { "data": { "text/plain": [ "525" ] }, "execution_count": 22, "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": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "39" ] }, "execution_count": 23, "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(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": 24, "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": 24, "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": 25, "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": 25, "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": 26, "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": 27, "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": 28, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "0e715d42cedf4f24b082a3742b6ba7a1", "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.4950910.00749166.0884830.00.4804080.509774
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.495091 0.007491 66.088483 0.0 0.480408 0.509774" ] }, "execution_count": 28, "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": 29, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c36de95d4af446f29d9f84c53133e9ed", "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.4897590.01000248.9682710.00.4701560.509362
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d 0.489759 0.010002 48.968271 0.0 0.470156 0.509362" ] }, "execution_count": 29, "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": 30, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "39ffcc1fef014cdc8010286cf9fe469c", "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.5492950.00873662.8788430.00.5321740.566417
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_diff 0.549295 0.008736 62.878843 0.0 0.532174 0.566417" ] }, "execution_count": 30, "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": 31, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "072390e43ee2426caf19fa0ffbb35cc4", "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.1339620.005109221.9708920.01.1239491.143975
\n", "" ], "text/plain": [ " coef std err t P>|t| 2.5 % 97.5 %\n", "d_demean 1.133962 0.005109 221.970892 0.0 1.123949 1.143975" ] }, "execution_count": 31, "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": 32, "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.495091 0.007491 0.480408 0.509774\n", " cre_normal 0.489759 0.010002 0.470156 0.509362\n", " fd_exact 0.549295 0.008736 0.532174 0.566417\n", " wg_approx 1.133962 0.005109 1.123949 1.143975\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAJOCAYAAABm7rQwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAorhJREFUeJzs3Xd0VNXexvFnJmXSSEJiEkgEEkIRpAVQpAkqGIpcsAAiShMUBRECFlRA4AoqRVCsCIiKSrFdQVEEUUGuUkSkS/fSEko6qXPeP3gzMiSQSTKZkPD9rJW1MvvsPed3pmzlyTn7mAzDMAQAAAAAAAC4kLmsCwAAAAAAAMDVh1AKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgCAcsJkMun5558v6zJwBWjfvr0aNGhQaL9Dhw7JZDLpvffeK/2icNV57733ZDKZtGnTprIuBQBQThFKAQDKjf379+vhhx9WzZo15eXlJX9/f7Vu3VqzZ8/WuXPnyro8XAXyQp68Hw8PD11zzTVq1aqVnnnmGR05cqSsS3SatWvX2o7zww8/LLBP69atZTKZ8gVkkZGRuuOOOy77/AMGDLB7Lf39/dW4cWPNmDFDmZmZtn7PP/+8TCaTTp065VCtee9LzZo11a9fPx04cMDW7+L3z2w2KygoSJ07d9aGDRsceVkAAIATuZd1AQAAOGLFihXq2bOnLBaL+vXrpwYNGigrK0vr1q3TE088oR07duidd94p6zJL1blz5+Tuzn+6rwR9+vRRly5dZLVadfbsWW3cuFGzZs3S7NmzNW/ePN17771lXaLTeHl56aOPPtL9999v137o0CH98ssv8vLyKvZzWywWvfvuu5KkxMREffrppxozZow2btyoTz75pMjPN2LECN1www3Kzs7Wli1b9M4772jFihX6888/FR4ebuuX9/7l5uZq7969euONN3TLLbdo48aNatiwYbGPBwAAFA3/ZwsAuOIdPHhQ9957r2rUqKE1a9aoatWqtm3Dhg3Tvn37tGLFijKssPRYrVZlZWXJy8urRP/4L+/S0tLk6+tb1mXYNG3aNF9Ic/jwYd1+++3q37+/6tWrp8aNG5dRdc7VpUsX/ec//9GpU6d0zTXX2No/+ugjhYWFqXbt2jp79myxntvd3d3udXz00UfVokULLV68WDNnzrQLkhzRtm1b3XPPPZKkgQMHqk6dOhoxYoQWLlyosWPH2vpd/P61bdtWnTt31ptvvqk33nijWMdSVjIyMuTp6SmzmQsgAADlD//1AgBc8V5++WWlpqZq3rx5doFUnlq1aunxxx+3Pc7JydHkyZMVHR0ti8WiyMhIPfPMM3aXBEn/XGK0du1aNW/eXN7e3mrYsKHWrl0rSfrss8/UsGFDeXl5qVmzZvr999/txg8YMEB+fn46cOCAYmNj5evrq/DwcE2aNEmGYdj1nT59ulq1aqXg4GB5e3urWbNmWrZsWb5jMZlMGj58uBYtWqTrr79eFotFK1eutG27cE2plJQUjRw5UpGRkbJYLAoNDVXHjh21ZcsWu+dcunSpmjVrJm9vb11zzTW6//77dfTo0QKP5ejRo+rRo4f8/PwUEhKiMWPGKDc39xLvzD++/PJLde3aVeHh4bJYLIqOjtbkyZMLHPvrr7+qS5cuqly5snx9fdWoUSPNnj07Xy379+9Xly5dVKlSJfXt21fS+XBq9OjRqlatmiwWi+rWravp06fne71XrVqlNm3aKDAwUH5+fqpbt66eeeYZuz6vvfaarr/+evn4+Khy5cpq3ry5Pvroo0KP9VJq1Kih9957T1lZWXr55Zftth04cEA9e/ZUUFCQfHx8dNNNN+ULUvPW5zl06JBde96laXmfywtt3rxZrVq1kre3t6KiovTWW285VOvu3bt1zz33KCgoSF5eXmrevLn+85//FNi3e/fuslgsWrp0qV37Rx99pF69esnNzc2hfTrCbDarffv2kpTvdSiOW2+9VdL5YPty2rZtK+n8JcKFOXPmjMaMGaOGDRvKz89P/v7+6ty5s/744w+7fnnv2+LFi/XMM8+oSpUq8vX11b/+9S/9/fffdn3z1ggr7P3Me85PPvlEzz33nCIiIuTj46Pk5GRJjn3Xt23bpgEDBtgug65SpYoGDRqk06dP5zvWo0eP6sEHH7R9r6OiovTII48oKyvLrl9mZqbi4uIUEhIiX19f3XnnnUpISMj3fN98843atm0rX19fVapUSV27dtWOHTvs+pw4cUIDBw7UtddeK4vFoqpVq6p79+5O+TwAAK48nCkFALjiffXVV6pZs6ZatWrlUP/Bgwdr4cKFuueeezR69Gj9+uuvmjp1qnbt2qXPP//cru++fft033336eGHH9b999+v6dOnq1u3bnrrrbf0zDPP6NFHH5UkTZ06Vb169dKePXvszkjIzc1Vp06ddNNNN+nll1/WypUrNWHCBOXk5GjSpEm2frNnz9a//vUv9e3bV1lZWfrkk0/Us2dPLV++XF27drWrac2aNVqyZImGDx+ua665RpGRkQUe59ChQ7Vs2TINHz5c9evX1+nTp7Vu3Trt2rVLTZs2lXQ+6Bg4cKBuuOEGTZ06VSdPntTs2bO1fv16/f777woMDLQ7ltjYWLVo0ULTp0/X999/rxkzZig6OlqPPPLIZV/z9957T35+foqLi5Ofn5/WrFmj8ePHKzk5WdOmTbP1W7Vqle644w5VrVpVjz/+uKpUqaJdu3Zp+fLl+YLF2NhYtWnTRtOnT5ePj48Mw9C//vUv/fDDD3rwwQfVpEkTffvtt3riiSd09OhRvfLKK5KkHTt26I477lCjRo00adIkWSwW7du3T+vXr7c9/9y5czVixAjdc889evzxx5WRkaFt27bp119/1X333XfZY72cli1bKjo6WqtWrbK1nTx5Uq1atVJ6erpGjBih4OBgLVy4UP/617+0bNky3XnnncXa19mzZ9WlSxf16tVLffr00ZIlS/TII4/I09NTgwYNuuS4HTt2qHXr1oqIiNDTTz8tX19fLVmyRD169NCnn36arx4fHx91795dH3/8se1z8Mcff2jHjh169913tW3btmLVfyl5wVBwcLDLnisv8KhcuXKhz3ngwAF98cUX6tmzp6KionTy5Em9/fbbateunXbu3Jnv7K4XXnhBJpNJTz31lOLj4zVr1ix16NBBW7dulbe3t61fUd7PyZMny9PTU2PGjFFmZqY8PT0d/q6vWrVKBw4c0MCBA1WlShXbpc87duzQf//7X5lMJknSsWPHdOONNyoxMVEPPfSQrrvuOh09elTLli1Tenq6PD09bfU89thjqly5siZMmKBDhw5p1qxZGj58uBYvXmzr88EHH6h///6KjY3VSy+9pPT0dL355ptq06aNfv/9d9s8d/fdd2vHjh167LHHFBkZqfj4eK1atUpHjhy55FwIACjHDAAArmBJSUmGJKN79+4O9d+6dashyRg8eLBd+5gxYwxJxpo1a2xtNWrUMCQZv/zyi63t22+/NSQZ3t7exuHDh23tb7/9tiHJ+OGHH2xt/fv3NyQZjz32mK3NarUaXbt2NTw9PY2EhARbe3p6ul09WVlZRoMGDYxbb73Vrl2SYTabjR07duQ7NknGhAkTbI8DAgKMYcOGXfK1yMrKMkJDQ40GDRoY586ds7UvX77ckGSMHz8+37FMmjTJ7jliYmKMZs2aXXIflzo+wzCMhx9+2PDx8TEyMjIMwzCMnJwcIyoqyqhRo4Zx9uxZu75WqzVfLU8//bRdny+++MKQZPz73/+2a7/nnnsMk8lk7Nu3zzAMw3jllVcMSXav/8W6d+9uXH/99YUe18UOHjxoSDKmTZt22eeWZCQlJRmGYRgjR440JBk///yzrU9KSooRFRVlREZGGrm5uYZhGMaCBQsMScbBgwftnu+HH37I99lr166dIcmYMWOGrS0zM9No0qSJERoaamRlZdnVu2DBAlu/2267zWjYsKHtfTGM869/q1atjNq1a+fb79KlS43ly5cbJpPJOHLkiGEYhvHEE08YNWvWtNVy8WtZo0YNo2vXrpd8jQzj/Pvs6+trJCQkGAkJCca+ffuMKVOmGCaTyWjUqJGt34QJEwp9P/NqnT9/vpGQkGAcO3bMWLFihREZGWmYTCZj48aNdq/HxIkTjYSEBOPEiRPGzz//bNxwww22Yy1MRkaG7T3Lc/DgQcNisdh9f/JqioiIMJKTk23tS5YsMSQZs2fPtrU5+n7mPWfNmjXtvnNF+a4X9F39+OOPDUnGTz/9ZGvr16+fYTabba/dhfK+r3mf2Q4dOth9h0eNGmW4ubkZiYmJhmGc/7wHBgYaQ4YMsXueEydOGAEBAbb2s2fPFvr9AgBULFy+BwC4ouVdllKpUiWH+n/99deSpLi4OLv20aNHS1K+S6bq16+vli1b2h63aNFC0vnLfqpXr56v/cI7eeUZPny47fe8y++ysrL0/fff29ovPiMiKSlJbdu2zXepnSS1a9dO9evXL+RIpcDAQP366686duxYgds3bdqk+Ph4Pfroo3brUXXt2lXXXXddgetwDR061O5x27ZtCzzmi114fCkpKTp16pTatm2r9PR07d69W5L0+++/6+DBgxo5cqTdGVqSbGdnXOjis7O+/vprubm5acSIEXbto0ePlmEY+uabbyTJ9txffvmlrFZrgfUGBgbqf//7nzZu3FjosRWVn5+fpPOvQ17dN954o9q0aWPX56GHHtKhQ4e0c+fOYu3H3d1dDz/8sO2xp6enHn74YcXHx2vz5s0Fjjlz5ozWrFmjXr162d6nU6dO6fTp04qNjdVff/2V73IvSbr99tsVFBSkTz75RIZh6JNPPlGfPn2KVfeF0tLSFBISopCQENWqVUvPPPOMWrZsme+MRkcNGjRIISEhCg8PV9euXZWWlqaFCxeqefPmdv0mTJigkJAQValSRW3bttWuXbs0Y8YM23pUl2OxWGxnS+bm5ur06dO2S0QL+j7369fPbv665557VLVqVdtclaco72f//v3tvnNF+a5fOC4jI0OnTp3STTfdJEm2+q1Wq7744gt169Yt32sn5f++PvTQQ3Ztbdu2VW5urg4fPizp/NlZiYmJ6tOnj+0zd+rUKbm5ualFixb64YcfbLV5enpq7dq1xV6nDABQvhBKAQCuaP7+/pL++Qd+YQ4fPiyz2axatWrZtVepUkWBgYG2fyTluTB4kqSAgABJUrVq1Qpsv/gfSmazWTVr1rRrq1OnjiT7NXGWL1+um266SV5eXgoKClJISIjefPNNJSUl5TuGqKiowg5T0vm1trZv365q1arpxhtv1PPPP28XIOUda926dfONve666/K9Fl5eXgoJCbFrq1y5skP/ONyxY4fuvPNOBQQEyN/fXyEhIbaFpPOOMe9SqgYNGhT6fO7u7rr22mvt2g4fPqzw8PB8AWW9evVs2yWpd+/eat26tQYPHqywsDDde++9WrJkiV1A9dRTT8nPz0833nijateurWHDhtld3lcSqampkv4JUg8fPlzge3Bx3UUVHh6eb/H3gj57F9q3b58Mw9C4ceNsYVDez4QJEyRJ8fHx+cZ5eHioZ8+e+uijj/TTTz/p77//LtFljnm8vLy0atUqrVq1yva869evz/edctT48eO1atUqrVmzRtu2bdOxY8f0wAMP5Ov30EMPadWqVfrqq680atQonTt3zqG106Tzgc0rr7yi2rVry2Kx6JprrlFISIi2bdtW4Pe5du3ado9NJpNq1aqV7z0qyvt58RxRlO/6mTNn9PjjjyssLEze3t4KCQmxPV9e/QkJCUpOTnbouyrln0fzLoPMmzv++usvSefD/os/d999953tM2exWPTSSy/pm2++UVhYmG6++Wa9/PLLOnHihEN1AADKH9aUAgBc0fz9/RUeHq7t27cXaVxBZ94U5FKLNF+q3bhoQW1H/Pzzz/rXv/6lm2++WW+88YaqVq0qDw8PLViwoMCFtS88k+FyevXqpbZt2+rzzz/Xd999p2nTpumll17SZ599ps6dOxe5zuIuWJ2YmKh27drJ399fkyZNUnR0tLy8vLRlyxY99dRTlzxb6XIuPBulqLy9vfXTTz/phx9+0IoVK7Ry5UotXrxYt956q7777ju5ubmpXr162rNnj5YvX66VK1fq008/1RtvvKHx48dr4sSJxdpvnu3btys0NNQWqDrqUp9ZR8MSR+S9F2PGjFFsbGyBfS4OdPPcd999euutt/T888+rcePGDp3NVxg3Nzd16NChxM+Tp2HDhg49X+3atW397rjjDrm5uenpp5/WLbfcUuCZQReaMmWKxo0bp0GDBmny5MkKCgqS2WzWyJEji/VZLw5H54iC9OrVS7/88oueeOIJNWnSRH5+frJarerUqVOx6y9svsx73g8++EBVqlTJ18/d/Z9/kowcOVLdunXTF198oW+//Vbjxo3T1KlTtWbNGsXExBSrPgDAlYtQCgBwxbvjjjv0zjvvaMOGDXaX2hWkRo0aslqt+uuvv2xnokjnF5tOTExUjRo1nFqb1WrVgQMHbGc0SNLevXslybYo76effiovLy99++23slgstn4LFiwo8f6rVq2qRx99VI8++qji4+PVtGlTvfDCC+rcubPtWPfs2WO7C1mePXv2OO21WLt2rU6fPq3PPvtMN998s6394jueRUdHSzof2hQniKhRo4a+//57paSk2J0tlXd54IXHYzabddttt+m2227TzJkzNWXKFD377LP64YcfbPv29fVV79691bt3b2VlZemuu+7SCy+8oLFjx9pdAlUUGzZs0P79+21nieXVtWfPnnx9L6477+ySxMREu36XOpPq2LFjSktLszu75uLP3sXyzkDy8PAo8nvQpk0bVa9eXWvXrtVLL71UpLFXumeffVZz587Vc889Z7vb5aUsW7ZMt9xyi+bNm2fXnpiYqGuuuSZf/7yzhPIYhqF9+/apUaNGdu3FeT/zOPpdP3v2rFavXq2JEydq/Pjxl6wxJCRE/v7+Rf5jwKXkffdDQ0Md+txFR0dr9OjRGj16tP766y81adJEM2bM0IcffuiUegAAVw4u3wMAXPGefPJJ+fr6avDgwTp58mS+7fv379fs2bMlSV26dJEkzZo1y67PzJkzJSnfne6cYc6cObbfDcPQnDlz5OHhodtuu03S+bMITCaT3Rkvhw4d0hdffFHsfebm5ua7VCg0NFTh4eHKzMyUJDVv3lyhoaF66623bG3S+duy79q1y2mvRd5ZEheeRZaVlaU33njDrl/Tpk0VFRWlWbNm5QteHDkDrUuXLsrNzbV7vSXplVdekclksp0ddubMmXxjmzRpIkm21+H06dN22z09PVW/fn0ZhqHs7OxCaynI4cOHNWDAAHl6euqJJ56wq/u3337Thg0bbG1paWl65513FBkZaTvjKO8f7j/99JOtX25urt55550C95eTk6O3337b9jgrK0tvv/22QkJC1KxZswLHhIaGqn379nr77bd1/PjxfNsTEhIueXwmk0mvvvqqJkyYUOAlceVZYGCgHn74YX377bfaunXrZfu6ubnl+7wuXbq0wLW4JOn999+3u/x42bJlOn78eL6zGYvzfuZx9Lte0HdVyj9fms1m9ejRQ1999ZU2bdqUb39FPWM0NjZW/v7+mjJlSoHfr7zPXXp6ujIyMuy2RUdHq1KlSnbHBQCoODhTCgBwxYuOjtZHH32k3r17q169eurXr58aNGigrKws/fLLL1q6dKkGDBggSWrcuLH69++vd955x3ZZ2W+//aaFCxeqR48euuWWW5xam5eXl1auXKn+/furRYsW+uabb7RixQo988wztvWZunbtqpkzZ6pTp0667777FB8fr9dff121atXStm3birXflJQUXXvttbrnnnvUuHFj+fn56fvvv9fGjRs1Y8YMSefPhnnppZc0cOBAtWvXTn369LHdJj4yMlKjRo1yymvQqlUrVa5cWf3799eIESNkMpn0wQcf5PuHq9ls1ptvvqlu3bqpSZMmGjhwoKpWrardu3drx44d+vbbby+7n27duumWW27Rs88+q0OHDqlx48b67rvv9OWXX2rkyJG2UGfSpEn66aef1LVrV9WoUUPx8fF64403dO2119oWG7/99ttVpUoVtW7dWmFhYdq1a5fmzJmjrl27OrSo/pYtW/Thhx/KarUqMTFRGzdu1Keffmo79gvPgnn66af18ccfq3PnzhoxYoSCgoK0cOFCHTx4UJ9++qntMsXrr79eN910k8aOHaszZ87YFhbPyckpsIbw8HC99NJLOnTokOrUqaPFixdr69ateuedd+Th4XHJ2l9//XW1adNGDRs21JAhQ1SzZk2dPHlSGzZs0P/+9z/98ccflxzbvXt3de/evdDXRzq/ftW///3vfO0xMTHFCkRnzpwpHx8fuzaz2axnnnmmyM9VkMcff1yzZs3Siy++qE8++eSS/e644w5NmjRJAwcOVKtWrfTnn39q0aJFl1wHKygoSG3atNHAgQN18uRJzZo1S7Vq1dKQIUPs+hX3/ZQc/677+/vb1mnKzs5WRESEvvvuu3xnNUrnL1P87rvv1K5dOz300EOqV6+ejh8/rqVLl2rdunX5blZwOf7+/nrzzTf1wAMPqGnTprr33nsVEhKiI0eOaMWKFWrdurXmzJmjvXv36rbbblOvXr1Uv359ubu76/PPP9fJkyd17733Orw/AEA5Uha3/AMAoDj27t1rDBkyxIiMjDQ8PT2NSpUqGa1btzZee+01u9vbZ2dnGxMnTjSioqIMDw8Po1q1asbYsWPt+hjGpW9bL8kYNmyYXVvereQvvFV53i3t9+/fb9x+++2Gj4+PERYWZkyYMCHfLePnzZtn1K5d27BYLMZ1111nLFiwwHar+8L2feG2CRMmGIZx/nbxTzzxhNG4cWOjUqVKhq+vr9G4cWPjjTfeyDdu8eLFRkxMjGGxWIygoCCjb9++xv/+9z+7PnnHcrGCaizI+vXrjZtuusnw9vY2wsPDjSeffNL49ttvDUnGDz/8YNd33bp1RseOHW11N2rUyHjttdcKrcUwzt9aftSoUUZ4eLjh4eFh1K5d25g2bZrd7ehXr15tdO/e3QgPDzc8PT2N8PBwo0+fPsbevXttfd5++23j5ptvNoKDgw2LxWJER0cbTzzxhJGUlHTZ48z7HOT9uLu7G0FBQUaLFi2MsWPHGocPHy5w3P79+4177rnHCAwMNLy8vIwbb7zRWL58eYH9OnToYFgsFiMsLMx45plnjFWrVuV7Hdu1a2dcf/31xqZNm4yWLVsaXl5eRo0aNYw5c+YUWO+CBQvy7adfv35GlSpVDA8PDyMiIsK44447jGXLltn6/PDDD4YkY+nSpZd9TfJquVCNGjXsXqcLfx588EHDMC7/Pl8o7zNY0I+bm1uRai3oe3yhAQMGGG5ubsa+ffsu+RwZGRnG6NGjjapVqxre3t5G69atjQ0bNhjt2rUz2rVrZ+uXV9PHH39sjB071ggNDTW8vb2Nrl275vucOPp+FnacjnzX//e//xl33nmnERgYaAQEBBg9e/Y0jh07Zje/5Dl8+LDRr18/IyQkxLBYLEbNmjWNYcOGGZmZmYZhGMaCBQsMScbGjRsLrPPi7/4PP/xgxMbGGgEBAYaXl5cRHR1tDBgwwNi0aZNhGIZx6tQpY9iwYcZ1111n+Pr6GgEBAUaLFi2MJUuWXPL9AACUbybDKMaKrQAAQAMGDNCyZctsd1sDgDxr167VLbfcoqVLl+qee+65bN/27dvr1KlTTlvDCQCA8oI1pQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuBxrSgEAAAAAAMDlOFMKAAAAAAAALkcoBQAAAAAAAJdzL+sCXM1qterYsWOqVKmSTCZTWZcDAAAAAABQoRiGoZSUFIWHh8tsvvT5UFddKHXs2DFVq1atrMsAAAAAAACo0P7++29de+21l9x+1YVSlSpVknT+hfH39y/jakrGarUqISFBISEhl00eAQCXxlwKACXHXAoAJVeR5tLk5GRVq1bNlsFcylUXSuVdsufv718hQqmMjAz5+/uX+w8sAJQV5lIAKDnmUgAouYo4lxa2bFLFOEoAAAAAAACUK4RSAAAAAAAAcDlCKQAAAAAAALjcVbemlKNyc3OVnZ1d1mVcltVqVXZ2tjIyMirM9aZXCg8PD7m5uZV1GQAAAAAAVFiEUhcxDEMnTpxQYmJiWZdSKMMwZLValZKSUujiYSi6wMBAValShdcWAAAAAIBSQCh1kbxAKjQ0VD4+Pld0IGEYhnJycuTu7n5F11neGIah9PR0xcfHS5KqVq1axhUBAAAAAFDxEEpdIDc31xZIBQcHl3U5hSKUKj3e3t6SpPj4eIWGhnIpHwAAAAAATsZCRBfIW0PKx8enjCvBlSDvc3Clry0GAAAAAEB5RChVAM46gsTnAAAAAACA0kQoBQAAAAAAAJcjlILLffHFF6pVq5bc3Nw0cuTIS7YBAAAAAICKi1CqAjCbzTKZTJf8ef7558u6RDsPP/yw7rnnHv3999+aPHnyJdtKYu3atTKZTEpMTCzxcwEAAAAAAOfj7nsVwLFjx2zrHy1evFjjx4/Xnj17bNv9/PxsvxuGodzcXLm7l81bn5qaqvj4eMXGxio8PPySbQAAAAAAoGLjTKkKoEqVKrafgIAAmUwm2+Pdu3erUqVK+uabb9SsWTNZLBatW7dOAwYMUI8ePeyeZ+TIkWrfvr3tsdVq1dSpUxUVFSVvb281btxYy5Ytu2wtmZmZGjNmjCIiIuTr66sWLVpo7dq1ks6fvVSpUiVJ0q233iqTyXTJNklat26d2rZtK29vb1WrVk0jRoxQWlqa3b6eeuopVatWTRaLRbVq1dK8efN06NAh3XLLLZKkypUry2QyacCAAcV/gQEAAAAAgNMRSl0lnn76ab344ovatWuXGjVq5NCYqVOn6v3339dbb72lHTt2aNSoUbr//vv1448/XnLM8OHDtWHDBn3yySfatm2bevbsqU6dOumvv/5Sq1atbGdwffrppzp+/Pgl2/bv369OnTrp7rvv1rZt27R48WKtW7dOw4cPt+2rX79++vjjj/Xqq69q165devvtt+Xn56dq1arp008/lSTt2bNHx48f1+zZs4v70gEAAAAAgFLA5XsOOvnzEZ1cd6TQfj4RlVSrX2O7tn3v/6H0oymFjg1rU11hbasXu8bLmTRpkjp27Ohw/8zMTE2ZMkXff/+9WrZsKUmqWbOm1q1bp7ffflvt2rXLN+bIkSNasGCBjhw5YrsMb8yYMVq5cqUWLFigKVOmKDQ0VJIUFBSkKlWqSFKBbVOnTlXfvn1ti57Xrl1br776qtq1a6c333xTR44c0ZIlS7Rq1Sp16NDBVl+eoKAg23MHBgY6fNwAAAAAAMA1CKUclJuZo+zkzEL75QRY8relZjk0Njczp1i1OaJ58+ZF6r9v3z6lp6fnC7KysrIUExNT4Jg///xTubm5qlOnjl17ZmamgoODi7T/P/74Q9u2bdOiRYtsbYZhyGq16uDBg/rzzz/l5uZWYDgGAAAAAACufIRSDnKzuMvDP3/gdDF3P88C2xwZ62YpvbfD19fX7rHZbJZhGHZt2dnZtt9TU1MlSStWrFBERIRdP4ul4GNJTU2Vm5ubNm/eLDc3N7ttFy627ojU1FQ9/PDDGjFiRL5t1atX1759+4r0fAAAAAAA4MpCKOWgsLbFv7Tu4sv5rgQhISHavn27XdvWrVvl4eEhSapfv74sFouOHDni8NlIMTExys3NVXx8vNq2bVui+po2baqdO3eqVq1aBW5v2LChrFarfvzxR9vlexfy9DwfDubm5paoDgAAAAAAUDoIpa5St956q6ZNm6b3339fLVu21Icffqjt27fbLs2rVKmSxowZo1GjRslqtapNmzZKSkrS+vXr5e/vr/79++d7zjp16qhv377q16+fZsyYoZiYGCUkJGj16tVq1KiRunbt6nB9Tz31lG666SYNHz5cgwcPlq+vr3bu3KlVq1Zpzpw5ioyMVP/+/TVo0CC9+uqraty4sQ4fPqz4+Hj16tVLNWrUkMlk0vLly9WlSxd5e3sX+WwtAAAAAACc7fTKmTq9cqZdm2EYkmFVbm6uktzcJJNZJpPJrk9wpzgFd4pzZamljrvvXaViY2M1btw4Pfnkk7rhhhuUkpKifv362fWZPHmyxo0bp6lTp6pevXrq1KmTVqxYoaioqEs+74IFC9SvXz+NHj1adevWVY8ePbRx40ZVr160s8waNWqkH3/8UXv37lXbtm0VExOj8ePH2xZQl6Q333xT99xzjx599FFdd911GjJkiNLS0iRJERERmjhxop5++mmFhYXZ3bUPAAAAAICyknsuWTlnj9r95CYeU27SCSk1QblJJ5SbeCx/n3PJZV2605mMixcWquCSk5MVEBCgpKQk+fv7223LyMjQwYMHFRUVJS8vrzKq0HGGYSgnJ0fu7u75ElSUXHn7PAAoHqvVqvj4eIWGhsps5m81AFAczKUA4Li8M6Ws2Rmypp6RdLlYxiSzX5DMHl7l6kypy2UvF+LyPQAAAAAAABcJ7hQnS8T1OjKzq2QySZc7V8hkkjU9UdfGrZBfw1jXFeki/BkDAAAAAADARXLTEvX3a3efD6MM6+U7G1bJMPT3a3crNy3RJfW5EqEUAAAAAACAiySuXygjM73wQCqPYZWRla7E9e+XbmFlgFAKAAAAAADABQzD0JlVrxVr7JlVr6qiLQtOKAUAAAAAAOACuamnlR2/X5df3LwAhqHs+P3KTTtTKnWVFUIpAAAAAAAAF7BmpJZs/LkUJ1VyZSCUAgAAAAAAcAGzl1/JxntXclIlVwb3si6gIpj543698tMBuzbDMGQ1JEOGTDLJbJJMJpNdn1E311Rcu2hXlgoAAAAAAMqIm1+wPEKjlZ1w4Pzd9xxlMskjpKbcfINKr7gyQCjlBMkZOTqalFGscQAAAAAA4OpgMpkU1PExnVw0qshjgzqOyHeyS3nH5XtO4O/lrogALwX7eKiwj4dJUrCPhyICvOTvRSZYlg4dOiSTyaStW7eWdSkAAAAAgKtEYOv+Mll8JJODkYzJLJOnjwJb9yvdwsoAoZQTxLWL1ru9GivxXI4KCy1NJinxXI7e7dWYS/cAAAAAALjKuPkGqtpjn54PCAoLpkxmyWRStcc+k5tvoEvqcyVCKSdIPJetexZukqHz60hdTt46U/cs3KTEc9muKbAA2dllt+/SlpWVVdYlAAAAAABQoNMrZ+rYvAdl9gksfF0pw5DZJ1DH5g3S6ZUzXVKfKxFKOcHCTX8rPSu30EAqj9WQ0rNy9f6mv51ah9Vq1csvv6xatWrJYrGoevXqeuGFF2yXqS1evFjt2rWTl5eXFi1aJEl69913Va9ePXl5eem6667TG2+84fD+fvnlFzVp0kReXl5q3ry5vvjii3yXw23fvl2dO3eWn5+fwsLC9MADD+jUqVO27e3bt9eIESP05JNPKigoSFWqVNHzzz9vt5/ExEQNHjxYISEh8vf316233qo//vjDtv35559XkyZN9O677yoqKkpeXl6SpJUrV6pNmzYKDAxUcHCw7rjjDu3fv78YrywAAAAAAM6Rey5ZOWePypp6WlJhQYIha+pp5Zw9qtxzya4oz6UIpUrIMAzNWXewWGNfW3dQRlFW2y/E2LFj9eKLL2rcuHHauXOnPvroI4WFhdm2P/3003r88ce1a9cuxcbGatGiRRo/frxeeOEF7dq1S1OmTNG4ceO0cOHCQveVnJysbt26qWHDhtqyZYsmT56sp556yq5PYmKibr31VsXExGjTpk1auXKlTp48qV69etn1W7hwoXx9ffXrr7/q5Zdf1qRJk7Rq1Srb9p49eyo+Pl7ffPONNm/erKZNm+q2227TmTNnbH327dunTz/9VJ999pktFEtLS1NcXJw2bdqk1atXy2w2684775TVai3OywsAAAAAQIm5efvLvXKE3Y9bYLjcAqpIfiFyC6git8Dw/H28/cu6dKdjpe0SOp2epf2n04s8zpC0/3S6zqRnK9jXs8R1pKSkaPbs2ZozZ4769+8vSYqOjlabNm106NAhSdLIkSN111132cZMmDBBM2bMsLVFRUVp586devvtt23PcSkfffSRTCaT5s6dKy8vL9WvX19Hjx7VkCFDbH3mzJmjmJgYTZkyxdY2f/58VatWTXv37lWdOnUkSY0aNdKECRMkSbVr19acOXO0evVqdezYUevWrdNvv/2m+Ph4WSwWSdL06dP1xRdfaNmyZXrooYcknb9k7/3331dISIhtX3fffbddzfPnz1dISIh27typBg0aOP7iAgAAAADgJMGd4hTcKS5fu9VqVXx8vEJDQ2U2Xx3nEBFKlVBqZm6Jxqdk5jgllNq1a5cyMzN12223XbJP8+bNbb+npaVp//79evDBB+2CpJycHAUEBBS6vz179qhRo0a2S+Uk6cYbb7Tr88cff+iHH36Qn59fvvH79++3C6UuVLVqVcXHx9ueIzU1VcHBwXZ9zp07Z3cpXo0aNewCKUn666+/NH78eP366686deqU7QypI0eOEEoBAAAAAFDGCKVKyM/iVqLxlSzOeQu8vb0L7ePr62v7PTU1VZI0d+5ctWjRwq6fm1vJjunCfXTr1k0vvfRSvm1Vq1a1/e7h4WG3zWQy2QKk1NRUVa1aVWvXrs33HIGBgbbfLzy2PN26dVONGjU0d+5chYeHy2q1qkGDBiyEDgAAAADAFYBQqoSCfTwVHeyjA6fTC12e7EImSTWDfRTk41FoX0fUrl1b3t7eWr16tQYPHlxo/7CwMIWHh+vAgQPq27dvkfdXt25dffjhh8rMzLRdVrdx40a7Pk2bNtWnn36qyMhIubsX76PWtGlTnThxQu7u7oqMjHR43OnTp7Vnzx7NnTtXbdu2lSStW7euWDUAAAAAAADnuzouUixFJpNJw9tEFWvsY22iZDKZnFKHl5eXnnrqKT355JN6//33tX//fv33v//VvHnzLjlm4sSJmjp1ql599VXt3btXf/75pxYsWKCZMwu/zeR9990nq9Wqhx56SLt27dK3336r6dOnS5LtmIYNG6YzZ86oT58+2rhxo/bv369vv/1WAwcOVG6uY5c9dujQQS1btlSPHj303Xff6dChQ/rll1/07LPPatOmTZccV7lyZQUHB+udd97Rvn37tGbNGsXF5b9mFwAAAAAAlA1CKSfo37yafDzdZHYwXzKbJB9PN/VrXs2pdYwbN06jR4/W+PHjVa9ePfXu3du2NlNBBg8erHfffVcLFixQw4YN1a5dO7333nuKiio8ZPP399dXX32lrVu3qkmTJnr22Wc1fvx4SbKtMxUeHq7169crNzdXt99+uxo2bKiRI0cqMDDQ4UXbTCaTvv76a918880aOHCg6tSpo3vvvVeHDx+2u7Pgxcxmsz755BNt3rxZDRo00KhRozRt2jSH9gkAAAAAAEqfyTCMolx1Vu4lJycrICBASUlJ8ve3v51iRkaGDh48qKioKLsFvB3x7Z543fHubzJkyHqZV9RskkwyacXgG3V73dDiHIKNYRjKycmRu7u70864KolFixZp4MCBSkpKcmiNqytdST4PAMqPq/EuJwDgbMylAFByFWkuvVz2cqHyfZRXiJk/7tfgJX8o0NtdhUV8hiEFervrwSV/aOaP+y/f+Qr3/vvva926dTp48KC++OILPfXUU+rVq1eFCKQAAAAAAEDpIpRyguSMHB1NytDp9OxCFzs3JJ1Oz9bRpAwlZ+S4orximTJlivz8/Ar86dy5syTpxIkTuv/++1WvXj2NGjVKPXv21DvvvFPGlQMAAAAAgPKAu+85gb+XuyIC7C/vMozzl/EZMmSS6fxlexddYufvdeW+/EOHDlWvXr0K3JZ3JtSTTz6pJ5980pVlAQAAAACACqJMU5GffvpJ06ZN0+bNm3X8+HF9/vnn6tGjxyX7Hz9+XKNHj9amTZu0b98+jRgxQrNmzXJZvZcS1y5ace2iy7oMpwoKClJQUFBZlwEAAAAAACqoMr18Ly0tTY0bN9brr7/uUP/MzEyFhIToueeeU+PGjUu5OgAAAAAAAJSWMj1TqnPnzrb1iRwRGRmp2bNnS5Lmz59fWmUBAAAAAACglF25ixo5SWZmpjIzM22Pk5OTJZ2/1aLVarXra7VaZRiG7ac8yKuzvNRbnuR9Dgr6rACoOPLmfr7nAFB8zKUAUHIVaS519BgqfCg1depUTZw4MV97QkKCMjIy7Nqys7NltVqVk5OjnBzH74yXsvVVpW59za7tfEhklQxDMpkkmfMtdO7X5DFVajLC4f1czDAM5ebmSsq/iDpKLicnR1arVadPn5aHh0dZlwOglFitViUlJckwDJnN3JQWAIqDuRQASq4izaUpKSkO9avwodTYsWMVFxdne5ycnKxq1aopJCRE/v7+dn0zMjKUkpIid3d3ubsX4aXJSVNu2rGiF5eTVrT9XAKBSelwd3eX2WxWcHCwvLy8Ch8AoFyyWq0ymUwKCQkp9//xB4CywlwKACVXkeZSR/8NXeFDKYvFIovFkq/dbDbne5PN5vNnM+X9OMrN4i83vwgZORmyZpyRdLlL6UwyewXJ5O4lN4t/ic5wMgzDNp4zpQq3du1a3XLLLTp79qwCAwML7Z/3OSjoswKgYuG7DgAlx1wKACVXUeZSR+uv8KGUKwQ0HSmPoPo6+WV3SSYVFkpZMxMV1ulL+UTe7qIKAQAAAAAArixlGr2lpqZq69at2rp1qyTp4MGD2rp1q44cOSLp/KV3/fr1sxuT1z81NVUJCQnaunWrdu7c6erS7eRmJCp+Re/z60epsMW8zq8zFb+it3IzEl1QXcGys7PLbN+Xk5WVVdYlAAAAAAAAFyjTUGrTpk2KiYlRTEyMJCkuLk4xMTEaP368JOn48eO2gCpPXv/Nmzfro48+UkxMjLp06eLy2i+UuusDGdnpKjyQymOVkZ2u1F0fOrUOq9Wql19+WbVq1ZLFYlH16tX1wgsv6NChQzKZTFq8eLHatWsnLy8vLVq0SJL07rvvql69evLy8tJ1112nN954w6F95T3nZ599pltuuUU+Pj5q3LixNmzYYNfv008/1fXXXy+LxaLIyEjNmDHDbntkZKQmT56sfv36yd/fXw899JDee+89BQYGavny5apbt658fHx0zz33KD09XQsXLlRkZKQqV66sESNG2BZ6l6QPPvhAzZs3V6VKlVSlShXdd999io+PL+GrCgAAAAAASkOZXr7Xvn37/79LXcHee++9fG2X618WDMNQ8tbXizU2eesc+TcZ5rT1oMaOHau5c+fqlVdeUZs2bXT8+HHt3r3btv3pp5/WjBkzFBMTYwumxo8frzlz5igmJka///67hgwZIl9fX/Xv39+hfT777LOaPn26ateurWeffVZ9+vTRvn375O7urs2bN6tXr156/vnn1bt3b/3yyy969NFHFRwcrAEDBtieY/r06Ro/frwmTJggSfr555+Vnp6uV199VZ988olSUlJ011136c4771RgYKC+/vprHThwQHfffbdat26t3r17Szp/9tfkyZNVt25dxcfHKy4uTgMGDNDXX3/tlNcXAAAAAAA4D2tKlZA147Rykg4UY6ShnKQDsmackZt3cInrSElJ0ezZszVnzhxboBQdHa02bdro0KFDkqSRI0fqrrvuso2ZMGGCZsyYYWuLiorSzp079fbbbzscSo0ZM0Zdu3aVJE2cOFHXX3+99u3bp+uuu04zZ87UbbfdpnHjxkmS6tSpo507d2ratGl2odStt96q0aNH2x7//PPPys7O1ptvvqno6GhJ0j333KMPPvhAJ0+elJ+fn+rXr69bbrlFP/zwgy2UGjRokO05atasqVdffVU33HCDUlNT5efnV5SXEwAAAAAAlLLyvZz7FcCalVrC8SlOqWPXrl3KzMzUbbfddsk+zZs3t/2elpam/fv368EHH5Sfn5/t59///rf279/v8H4bNWpk+71q1aqSZLtkbteuXWrdurVd/9atW+uvv/6yu+zuwrry+Pj42AIpSQoLC1NkZKRduBQWFmZ3ed7mzZvVrVs3Va9eXZUqVVK7du0kKd8loAAAAAAAoOxxplQJmT1LdgaO2bOSU+rw9vYutI+vr6/t99TU82Ha3Llz1aJFC7t+bm5uDu/Xw8PD9nveZYhWq6Nra+Wvq6DnzXvugtry9pWWlqbY2FjFxsZq0aJFCgkJ0ZEjRxQbG8vi6QAAAAAAXIEIpUrI7BUs94Caykk6KKko612Z5B4QJbNXkFPqqF27try9vbV69WoNHjy40P5hYWEKDw/XgQMH1LdvX6fUcLF69epp/fr1dm3r169XnTp1ihR8OWL37t06ffq0XnzxRVWrVk3S+YX0AQAAAADAlYlQqoRMJpP8mwzTmR/HFHmsf5PhTlvk3MvLS0899ZSefPJJeXp6qnXr1kpISNCOHTsueUnfxIkTNWLECAUEBKhTp07KzMzUpk2bdPbsWcXFxZW4ptGjR+uGG27Q5MmT1bt3b23YsEFz5sxx+A5/RVG9enV5enrqtdde09ChQ7V9+3ZNnjzZ6fsBAAAAAADOwZpSTuBX7wGZPHzk+MtplsnDR3717ndqHePGjdPo0aM1fvx41atXT71797Zbc+ligwcP1rvvvqsFCxaoYcOGateund577z1FRUU5pZ6mTZtqyZIl+uSTT9SgQQONHz9ekyZNslvk3FlCQkL03nvvaenSpapfv75efPFFTZ8+3en7AQAAAAAAzmEyDKMo15yVe8nJyQoICFBSUpL8/f3ttmVkZOjgwYOKioqSl5dXkZ43/dB3Ovlld8kwJF1uTSWzZDIprMd/5FOjY9EP4AKGYSgnJ0fu7u5OO+MK/yjJ5wFA+WG1WhUfH6/Q0FCZzfytBgCKg7kUAEquIs2ll8teLlS+j/IKkbRllk59/7DMlkAVvq6UIbMlUKdWPaSkLbNKvzgAAAAAAIArEKGUE1gzk5WbelTWjNNyJJSyZpw+3z8z2RXlFcuUKVPk5+dX4E/nzp3LujwAAAAAAFDOsdC5E5gt/nLzi7BrM/Iu4zMMyWSSZM53iZ3ZculT2Mra0KFD1atXrwK3eXt7u7gaAAAAAABQ0RBKOUFA05EKaDqyrMtwqqCgIAUFBZV1GQAAAAAAoILi8j0AAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC7HmlJO8Mr2H/XKjp/s2gzDkFWGbZ1zs0z5Fjofdf3NGtWgnStLBQAAAAAAuCIQSjlBcnaGjqYnFWscAAAAAADA1YjL95zA38NLET4BCrb4yFRIX5OkYIuPInwC5O/h5bQaDMPQQw89pKCgIJlMJm3dujVfn/bt22vkyJFO2ycAAAAAAEBxEUo5wagG7TS3TS8lZp2TqZBYyiSTErPOaW6bXk69dG/lypV67733tHz5ch0/flwNGjRw2nOXhUOHDl0yXAMAAAAAAOUfoZQTJGaeU881C2UYklXGZfvmrTPVc81CJWaec1oN+/fvV9WqVdWqVStVqVJF7u5cmQkAAAAAAK5chFJO8P6+TUrPySo0kMpjlaH0nCx9sH+TU/Y/cOBAPfbYYzpy5IhMJpMiIyOVlpamfv36yc/PT1WrVtWMGTOK9JyZmZkaM2aMIiIi5OvrqxYtWmjt2rWSpIyMDF1//fV66KGHbP3379+vSpUqaf78+ZKk06dPq0+fPoqIiJCPj48aNmyojz/+2G4fVqtVL7/8smrVqiWLxaLq1avrhRdekCRFRUVJkmJiYmQymdS+fftivjoAAAAAAOBKRChVQoZhaM6udcUa+9rOdTIMx4Ksy5k1a5YmTZqka6+9VsePH9fGjRv1xBNP6Mcff9SXX36p7777TmvXrtWWLVscfs7hw4drw4YN+uSTT7Rt2zb17NlTnTp10l9//SUvLy8tWrRICxcu1Jdffqnc3Fzdf//96tixowYNGiTpfHDVrFkzrVixQtu3b9dDDz2kBx54QL/99pttH2PHjtWLL76ocePGaefOnfroo48UFhYmSbZ+33//vY4fP67PPvusxK8TAAAAAAC4cnCNVwmdzkzX/pTTRR5nSNqfclpnMtMV7OVbohoCAgJUqVIlubm5qUqVKkpNTdW8efP04Ycf6rbbbpMkLVy4UNdee61Dz3fkyBEtWLBAR44cUXh4uCRpzJgxWrlypRYsWKApU6aoSZMm+ve//63Bgwfr3nvv1eHDh7V8+XLbc0RERGjMmDG2x4899pi+/fZbLVmyRDfeeKNSUlI0e/ZszZkzR/3795ckRUdHq02bNpKkkJAQSVJwcLCqVKlSotcHAAAAAABceQilSig1O7NE41OyM0scSl1s//79ysrKUosWLWxtQUFBqlu3rkPj//zzT+Xm5qpOnTp27ZmZmQoODrY9Hj16tL744gvNmTNH33zzjd223NxcTZkyRUuWLNHRo0eVlZWlzMxM+fj4SJJ27dqlzMxMW2gGAAAAAACuLoRSJeTnYSnR+EolHF8aUlNT5ebmps2bN8vNzc1um5+fn+33+Ph47d27V25ubvrrr7/UqVMn27Zp06Zp9uzZmjVrlho2bChfX1+NHDlSWVlZkiRvb2/XHAwAAAAAALgisaZUCQVbfBRdKVimIo4zSYquFKwgi4/Ta4qOjpaHh4d+/fVXW9vZs2e1d+9eh8bHxMQoNzdX8fHxqlWrlt3PhZfSDRo0SA0bNtTChQv11FNPadeuXbZt69evV/fu3XX//fercePGqlmzpt3+a9euLW9vb61evbrAGjw9PSWdP+MKAAAAAABUPJwpVUImk0nD67VR3G9fFnnsY/XbyGQqapxVOD8/Pz344IN64oknFBwcrNDQUD377LMymx3LIOvUqaO+ffuqX79+mjFjhmJiYpSQkKDVq1erUaNG6tq1q15//XVt2LBB27ZtU7Vq1bRixQr17dtX//3vf+Xp6anatWtr2bJl+uWXX1S5cmXNnDlTJ0+eVP369SVJXl5eeuqpp/Tkk0/K09NTrVu3VkJCgnbs2KEHH3xQoaGh8vb21sqVK3XttdfKy8tLAQEBTn+tAAAAAABA2eBMKSfoV6u5fNw9ZXbwfCmzTPJx99QD0c1LraZp06apbdu26tatmzp06KA2bdqoWbNmDo9fsGCB+vXrp9GjR6tu3brq0aOHNm7cqOrVq2v37t164okn9MYbb6hatWqSpDfeeEOnTp3SuHHjJEnPPfecmjZtqtjYWLVv315VqlRRjx497PYxbtw4jR49WuPHj1e9evXUu3dvxcfHS5Lc3d316quv6u2331Z4eLi6d+/unBcGAAAAAABcEUyGYRhlXYQrJScnKyAgQElJSfL397fblpGRoYMHDyoqKkpeXl5Fet5vj+5Rt1XvyjAkqy79kpplkskkLe84WLdHOLbw+KUYhqGcnBy5u7uXyhlXV7uSfB4AlB9Wq1Xx8fEKDQ11+IxSAIA95lIAKLmKNJdeLnu5UPk+yivEK9t/1JB1SxTo6S3jMoGUJBkyFOjprcHrluiV7T+6qEIAAAAAAIArC2tKOUFydoaOpic51NeQdDoz3TauLPz888/q3LnzJbenpqa6sBoAAAAAAHA1IpRyAn8PL0X42C/CbRiGrDJkGJLJlHfZninfuLLQvHlzbd26tUz2DQAAAAAAIBFKOcWoBu00qkG7si7DYd7e3qpVq1ZZlwEAAAAAAK5irCkFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlWFPKCU6vnKnTK2fatRmGIRlW2VY6N5nzLXQe3ClOwZ3iXFkqAAAAAADAFYFQyglyzyUr5+zRYo0DAAAAAAC4GhFKOYGbt7/cK0fImp0ha+oZScZleptk9guS2cNLbt7+rioRAAAAAADgikIo5QTBneJkibheR2Z2PX+pnnGZUMpkkjU9UdfGrZBfw1jXFQkAAAAAAHAFYaFzJ8hNS9Tfr919PowyrJfv/P/rTP392t3KTUt0SX0VRXZ2dlmXAAAAAAAAnIRQygkS1y+UkZleeCCVx7DKyEpX4vr3nbL/5cuXKzAwULm5uZKkrVu3ymQy6emnn7b1GTx4sO6//35J0ty5c1WtWjX5+Pjozjvv1MyZMxUYGOjQvvbv36/u3bsrLCxMfn5+uuGGG/T999/b9YmMjNTkyZPVp08f+fr6KiIiQq+//rpdH5PJpDfffFOdO3eWt7e3atasqWXLltm2Hzp0SCaTSYsXL1a7du3k5eWlRYsWyWq1atKkSbr22mtlsVjUpEkTrVy5UtL5xeU7dOig2NjY8wvNSzpz5oyuvfZajR8/vmgvKgAAAAAAKFWEUiVkGIbOrHqtWGPPrHrVFp6URNu2bZWSkqLff/9dkvTjjz/qmmuu0dq1a219fvzxR7Vv317r16/X0KFD9fjjj2vr1q3q2LGjXnjhBYf3lZqaqi5dumj16tX6/fff1alTJ3Xr1k1Hjhyx6zdt2jQ1btxYv//+u55++mk9/vjjWrVqlV2fcePG6e6779Yff/yhvn376t5779WuXbvs+uSN3bVrl2JjYzV79mzNmDFD06dP17Zt2xQbG6t//etf+uuvv2QymbRw4UJt3LhRr776qiRp6NChioiIIJQCAAAAAOAKYzKckYqUI8nJyQoICFBSUpL8/e0XGs/IyNDBgwcVFRUlLy8vh54vJ+WU9g4PKXY9dV4/JXe/4GKNNQxDOTk5cnd3V/PmzdWnTx+NGTNGd955p2644QZNnDhRp0+fVlJSkq699lrt3btX48aNU2pqqpYvX257nvvvv1/Lly9XYmJisepo0KCBhg4dquHDh0s6f6ZUvXr19M0339j63HvvvUpOTtbXX38t6fyZUkOHDtWbb75p63PTTTepadOmeuONN3To0CFFRUVp1qxZevzxx219IiIiNGzYMD3zzDO2thtvvFE33HCD7WyspUuXql+/fho5cqRee+01/f7776pdu3aRj6s4nwcA5Y/ValV8fLxCQ0NlNvO3GgAoDuZSACi5ijSXXi57uVD5PsorgDUjtWTjz6U4pY527dpp7dq1MgxDP//8s+666y7Vq1dP69at048//qjw8HDVrl1be/bs0Y033mg39uLHl5OamqoxY8aoXr16CgwMlJ+fn3bt2pXvTKmWLVvme3zxWVCO9GnevLnt9+TkZB07dkytW7e269O6dWu7cT179tSdd96pF198UdOnTy9WIAUAAAAAAEoXd98rIbOXX8nGe1dySh3t27fX/Pnz9ccff8jDw0PXXXed2rdvr7Vr1+rs2bNq166dU/YzZswYrVq1StOnT1etWrXk7e2te+65R1lZWU55/ov5+voWeUx6ero2b94sNzc3/fXXX6VQFQAAAAAAKCnOlCohN79geYRGSyZT0QaaTPIIjZabb5BT6shbV+qVV16xBVB5odTatWvVvn17SVLdunW1ceNGu7EXP76c9evXa8CAAbrzzjvVsGFDValSRYcOHcrX77///W++x/Xq1Stynwv5+/srPDxc69evz1dT/fr1bY9Hjx4ts9msb775Rq+++qrWrFnj6OEBAAAAAAAX4UypEjKZTArq+JhOLhpV5LFBHUfIVNQw6xIqV66sRo0aadGiRZozZ44k6eabb1avXr2UnZ1tC6oee+wx3XzzzZo5c6a6deumNWvW6JtvvnG4jtq1a+uzzz5Tt27dZDKZNG7cOFmt+e86uH79er388svq0aOHVq1apaVLl2rFihV2fZYuXarmzZurTZs2WrRokX777TfNmzfvsvt/4oknNGHCBEVHR6tJkyZasGCBtm7dqkWLFkmSVqxYofnz52vDhg1q2rSpnnjiCfXv31/btm1T5cqVHTpGAAAAAABQ+jhTygkCW/eXyeIjmRx8OU1mmTx9FNi6n1PraNeunXJzc21nRQUFBal+/fqqUqWK6tatK+n8+ktvvfWWZs6cqcaNG2vlypUaNWqUwwt5z5w5U5UrV1arVq3UrVs3xcbGqmnTpvn6jR49Wps2bVJMTIz+/e9/a+bMmYqNjbXrM3HiRH3yySdq1KiR3n//fX388cd2ZzwVZMSIEYqLi9Po0aPVsGFDrVy5Uv/5z39Uu3ZtJSQk6MEHH9Tzzz9vq2nixIkKCwvT0KFDHTo+AAAAAADgGtx97wIludta6p/f6sjMrpJhSEb+M4dsTGbJZFL1uK/l1/D24hyCzYV33yvJGVdDhgzR7t279fPPP5eonjyRkZEaOXKkRo4ceck+JpNJn3/+uXr06OGUfZYG7r4HXB0q0l1OAKCsMJcCQMlVpLmUu++50OmVM3Vs3oMy+wSeD6UuxzBk9gnUsXmDdHrlTJfUd7Hp06frjz/+0L59+/Taa69p4cKF6t+/f5nUAgAAAAAArk6sKeUEueeSlXP2qIO9DVlTT8v6/+PKwm+//aaXX35ZKSkpqlmzpl599VUNHjxYknT99dfr8OHDBY57++231bdvX1eWCgAAAAAAKihCKSdw8/aXe+UIuzYj7zI+wzh/Zz6TOd8ldm7elz6FrTQtWbLkktu+/vprZWdnF7gtLCzMoecv6G58F7vKrhoFAAAAAAAXIZRyguBOcQruFFfWZThFjRo1yroEAAAAAABwFWBNKQAAAAAAALgcoVQBuLQMEp8DAAAAAABKE6HUBTw8PCRJ6enpZVwJrgR5n4O8zwUAAAAAAHCeMl1T6qefftK0adO0efNmHT9+XJ9//rl69Ohx2TFr165VXFycduzYoWrVqum5557TgAEDnFKPm5ubAgMDFR8fL0ny8fHJtzj5lcQwDOXk5Mjd3f2KrrO8MQxD6enpio+PV2BgoNzc3Mq6JAAAAAAAKpwyDaXS0tLUuHFjDRo0SHfddVeh/Q8ePKiuXbtq6NChWrRokVavXq3BgweratWqio2NdUpNVapUkSRbMHUlMwxDVqtVZnP+O/uh5AIDA22fBwAAAAAA4FxlGkp17txZnTt3drj/W2+9paioKM2YMUOSVK9ePa1bt06vvPKK00Ipk8mkqlWrKjQ0VNnZ2U55ztJitVp1+vRpBQcHy2zmSkxn8vDw4AwpAAAAAABKUZmGUkW1YcMGdejQwa4tNjZWI0eOdPq+3NzcrvhQwmq1ysPDQ15eXoRSAAAAAACgXClXodSJEycUFhZm1xYWFqbk5GSdO3dO3t7e+cZkZmYqMzPT9jg5OVnS+UDHarWWbsGlzGq12i7hAwAUD3MpAJQccykAlFxFmksdPYZyFUoVx9SpUzVx4sR87QkJCcrIyCiDipzHarUqKSlJhmFwphQAFBNzKQCUHHMpAJRcRZpLU1JSHOpXrkKpKlWq6OTJk3ZtJ0+elL+/f4FnSUnS2LFjFRcXZ3ucnJysatWqKSQkRP7+/qVab2mzWq0ymUwKCQkp9x9YACgrzKUAUHLMpQBQchVpLvXy8nKoX7kKpVq2bKmvv/7arm3VqlVq2bLlJcdYLBZZLJZ87Wazudy/ydL5hdkryrEAQFlhLgWAkmMuBYCSqyhzqaP1l+lRpqamauvWrdq6dask6eDBg9q6dauOHDki6fxZTv369bP1Hzp0qA4cOKAnn3xSu3fv1htvvKElS5Zo1KhRZVE+AAAAAAAAiqlMQ6lNmzYpJiZGMTExkqS4uDjFxMRo/PjxkqTjx4/bAipJioqK0ooVK7Rq1So1btxYM2bM0LvvvqvY2NgyqR8AAAAAAADFU6aX77Vv316GYVxy+3vvvVfgmN9//70UqwIAAAAAAEBpK98XKQIAAAAAAKBcIpQCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFyOUAoAAAAAAAAuRygFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlCKUAAAAAAADgcoRSAAAAAAAAcDlCKQAAAAAAALgcoRQAAAAAAABcjlAKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFyOUAoAAAAAAAAuRygFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlCKUAAAAAAADgcoRSAAAAAAAAcDlCKQAAAAAAALgcoRQAAAAAAABcjlAKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFyOUAoAAAAAAAAuRygFAAAAAAAAl7siQqnXX39dkZGR8vLyUosWLfTbb79dsm92drYmTZqk6OhoeXl5qXHjxlq5cqULqwUAAAAAAEBJlXkotXjxYsXFxWnChAnasmWLGjdurNjYWMXHxxfY/7nnntPbb7+t1157TTt37tTQoUN155136vfff3dx5QAAAAAAACiuMg+lZs6cqSFDhmjgwIGqX7++3nrrLfn4+Gj+/PkF9v/ggw/0zDPPqEuXLqpZs6YeeeQRdenSRTNmzHBx5QAAAAAAACiuMg2lsrKytHnzZnXo0MHWZjab1aFDB23YsKHAMZmZmfLy8rJr8/b21rp160q1VgAAAAAAADiPe1nu/NSpU8rNzVVYWJhde1hYmHbv3l3gmNjYWM2cOVM333yzoqOjtXr1an322WfKzc0tsH9mZqYyMzNtj5OTkyVJVqtVVqvVSUdSNqxWqwzDKPfHAQBlibkUAEqOuRQASq4izaWOHkOZhlLFMXv2bA0ZMkTXXXedTCaToqOjNXDgwEte7jd16lRNnDgxX3tCQoIyMjJKu9xSZbValZSUJMMwZDaX+ZWYAFAuMZcCQMkxlwJAyVWkuTQlJcWhfmUaSl1zzTVyc3PTyZMn7dpPnjypKlWqFDgmJCREX3zxhTIyMnT69GmFh4fr6aefVs2aNQvsP3bsWMXFxdkeJycnq1q1agoJCZG/v7/zDqYMWK1WmUwmhYSElPsPLACUFeZSACg55lIAKLmKNJdevOzSpZRpKOXp6almzZpp9erV6tGjh6Tzb8Lq1as1fPjwy4718vJSRESEsrOz9emnn6pXr14F9rNYLLJYLPnazWZzuX+TJclkMlWYYwGAssJcCgAlx1wKACVXUeZSR+sv88v34uLi1L9/fzVv3lw33nijZs2apbS0NA0cOFCS1K9fP0VERGjq1KmSpF9//VVHjx5VkyZNdPToUT3//POyWq168skny/IwAAAAAAAAUARlHkr17t1bCQkJGj9+vE6cOKEmTZpo5cqVtsXPjxw5YpewZWRk6LnnntOBAwfk5+enLl266IMPPlBgYGAZHQEAAAAAAACKymQYhlHWRbhScnKyAgIClJSUVCHWlIqPj1doaGi5P7UPAMoKcykAlBxzKQCUXEWaSx3NXop1lPv379dzzz2nPn36KD4+XpL0zTffaMeOHcWrFgAAAAAAAFeVIodSP/74oxo2bKhff/1Vn332mVJTUyVJf/zxhyZMmOD0AgEAAAAAAFDxFDmUevrpp/Xvf/9bq1atkqenp6391ltv1X//+1+nFgcAAAAAAICKqcih1J9//qk777wzX3toaKhOnTrllKIAAAAAAABQsRU5lAoMDNTx48fztf/++++KiIhwSlEAAAAAAACo2IocSt1777166qmndOLECZlMJlmtVq1fv15jxoxRv379SqNGAAAAAAAAVDBFDqWmTJmi6667TtWqVVNqaqrq16+vm2++Wa1atdJzzz1XGjUCAAAAAACggnEv6gBPT0/NnTtX48aN0/bt25WamqqYmBjVrl27NOoDAAAAAABABVTkUCpP9erVVb16dWfWAgAAAAAAgKtEkUOpQYMGXXb7/Pnzi10MAAAAAAAArg5FDqXOnj1r9zg7O1vbt29XYmKibr31VqcVBgAAAAAAgIqryKHU559/nq/NarXqkUceUXR0tFOKAgAAAAAAQMVW5LvvFfgkZrPi4uL0yiuvOOPpAAAAAAAAUME5JZSSpP379ysnJ8dZTwcAAAAAAIAKrMiX78XFxdk9NgxDx48f14oVK9S/f3+nFQYAAAAAAICKq8ih1O+//2732Gw2KyQkRDNmzCj0znwAAAAAAACAVIxQ6ocffiiNOgAAAAAAAHAVcdqaUgAAAAAAAICjHDpTKiYmRiaTyaEn3LJlS4kKAgAAAAAAQMXnUCjVo0ePUi4DAAAAAAAAVxOHQqkJEyaUdh0AAAAAAAC4irCmFAAAAAAAAFyuyHffy83N1SuvvKIlS5boyJEjysrKstt+5swZpxUHAAAAAACAiqnIZ0pNnDhRM2fOVO/evZWUlKS4uDjdddddMpvNev7550uhRAAAAAAAAFQ0RQ6lFi1apLlz52r06NFyd3dXnz599O6772r8+PH673//Wxo1AgAAAAAAoIIpcih14sQJNWzYUJLk5+enpKQkSdIdd9yhFStWOLc6AAAAAAAAVEhFDqWuvfZaHT9+XJIUHR2t7777TpK0ceNGWSwW51YHAAAAAACACqnIodSdd96p1atXS5Iee+wxjRs3TrVr11a/fv00aNAgpxcIAAAAAACAisfhu+/NmTNH999/v1588UVbW+/evVW9enVt2LBBtWvXVrdu3UqlSAAAAAAAAFQsDp8p9eyzzyo8PFx9+/bVmjVrbO0tW7ZUXFwcgRQAAAAAAAAc5nAodeLECb311ls6duyYOnbsqKioKE2ePFl///13adYHAAAAAACACsjhUMrb21v9+vXTDz/8oL/++ksPPPCA5s2bp6ioKHXq1ElLly5VdnZ2adYKAAAAAACACqLIC51LUs2aNTVp0iQdPHhQ33zzjYKDgzVgwABFREQ4uz4AAAAAAABUQMUKpfKYTCa5u7vLZDLJMAzOlAIAAAAAAIBDihVK/f3335o0aZJq1qypjh076tixY5o7d66OHz/u7PoAAAAAAABQAbk72jErK0ufffaZ5s+frzVr1qhq1arq37+/Bg0apJo1a5ZmjQAAAAAAAKhgHA6lqlSpovT0dN1xxx366quvFBsbK7O5RFf/AQAAAAAA4CrlcCj13HPP6YEHHlBISEhp1gMAAAAAAICrgMOhVFxcXGnWAQAAAAAAgKsI198BAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHC5IodSkyZNUnp6er72c+fOadKkSU4pCgAAAAAAABVbkUOpiRMnKjU1NV97enq6Jk6c6JSiAAAAAAAAULEVOZQyDEMmkylf+x9//KGgoCCnFAUAAAAAAICKzd3RjpUrV5bJZJLJZFKdOnXsgqnc3FylpqZq6NChpVIkAAAAAAAAKhaHQ6lZs2bJMAwNGjRIEydOVEBAgG2bp6enIiMj1bJly1IpEgAAAAAAABWLw6FU//79JUlRUVFq1aqVPDw8Sq0oAAAAAAAAVGwOh1J52rVrJ6vVqr179yo+Pl5Wq9Vu+8033+y04gAAAAAAAFAxFTmU+u9//6v77rtPhw8flmEYdttMJpNyc3OdVhwAAAAAAAAqpiKHUkOHDlXz5s21YsUKVa1atcA78QEAAAAAAACXU+RQ6q+//tKyZctUq1at0qgHAAAAAAAAVwFzUQe0aNFC+/btK41aAAAAAAAAcJUo8plSjz32mEaPHq0TJ06oYcOG+e7C16hRI6cVBwAAAAAAgIqpyKHU3XffLUkaNGiQrc1kMskwDBY6BwAAAAAAgEOKHEodPHiwNOoAAAAAAADAVaTIoVSNGjVKow4AAAAAAABcRYq80LkkffDBB2rdurXCw8N1+PBhSdKsWbP05ZdfFquI119/XZGRkfLy8lKLFi3022+/Xbb/rFmzVLduXXl7e6tatWoaNWqUMjIyirVvAAAAAAAAuF6RQ6k333xTcXFx6tKlixITE21rSAUGBmrWrFlFLmDx4sWKi4vThAkTtGXLFjVu3FixsbGKj48vsP9HH32kp59+WhMmTNCuXbs0b948LV68WM8880yR9w0AAAAAAICyUeRQ6rXXXtPcuXP17LPPys3NzdbevHlz/fnnn0UuYObMmRoyZIgGDhyo+vXr66233pKPj4/mz59fYP9ffvlFrVu31n333afIyEjdfvvt6tOnT6FnVwEAAAAAAODKUeRQ6uDBg4qJicnXbrFYlJaWVqTnysrK0ubNm9WhQ4d/CjKb1aFDB23YsKHAMa1atdLmzZttIdSBAwf09ddfq0uXLkXaNwAAAAAAAMpOkRc6j4qK0tatW/MteL5y5UrVq1evSM916tQp5ebmKiwszK49LCxMu3fvLnDMfffdp1OnTqlNmzYyDEM5OTkaOnToJS/fy8zMVGZmpu1xcnKyJMlqtcpqtRap3iuN1WqVYRjl/jgAoCwxlwJAyTGXAkDJVaS51NFjKHIoFRcXp2HDhikjI0OGYei3337Txx9/rKlTp+rdd98tcqFFtXbtWk2ZMkVvvPGGWrRooX379unxxx/X5MmTNW7cuHz9p06dqokTJ+ZrT0hIKPeLo1utViUlJckwDJnNxVqzHgCuesylAFByzKUAUHIVaS5NSUlxqJ/JMAyjqE++aNEiPf/889q/f78kKTw8XBMnTtSDDz5YpOfJysqSj4+Pli1bph49etja+/fvr8TExALv5te2bVvddNNNmjZtmq3tww8/1EMPPaTU1NR8b1xBZ0pVq1ZNZ8+elb+/f5HqvdJYrVYlJCQoJCSk3H9gAaCsMJcCQMkxlwJAyVWkuTQ5OVmVK1dWUlLSZbOXIp8pJUl9+/ZV3759lZ6ertTUVIWGhharSE9PTzVr1kyrV6+2hVJWq1WrV6/W8OHDCxyTnp6e783JW3C9oHzNYrHIYrHkazebzeX+TZYkk8lUYY4FAMoKcykAlBxzKQCUXEWZSx2tv1ihVB4fHx/5+PiU5CkUFxen/v37q3nz5rrxxhs1a9YspaWlaeDAgZKkfv36KSIiQlOnTpUkdevWTTNnzlRMTIzt8r1x48apW7dudncDBAAAAAAAwJXLoVCqadOmWr16tSpXrqyYmBiZTKZL9t2yZUuRCujdu7cSEhI0fvx4nThxQk2aNNHKlStti58fOXLELmF77rnnZDKZ9Nxzz+no0aMKCQlRt27d9MILLxRpvwAAAAAAACg7DoVS3bt3t10Cd+HaT84yfPjwS16ut3btWrvH7u7umjBhgiZMmOD0OgAAAAAAAOAaDoVSFwZAhEEAAAAAAAAoqSKvnLVx40b9+uuv+dp//fVXbdq0ySlFAQAAAAAAoGIrcig1bNgw/f333/najx49qmHDhjmlKAAAAAAAAFRsRQ6ldu7cqaZNm+Zrj4mJ0c6dO51SFAAAAAAAACq2IodSFotFJ0+ezNd+/Phxubs7tEQVAAAAAAAArnJFDqVuv/12jR07VklJSba2xMREPfPMM+rYsaNTiwMAAAAAAEDFVORTm6ZPn66bb75ZNWrUUExMjCRp69atCgsL0wcffOD0AgEAAAAAAFDxFDmUioiI0LZt27Ro0SL98ccf8vb21sCBA9WnTx95eHiURo0AAAAAAACoYIq1CJSvr68eeughZ9cCAAAAAACAq4RDodR//vMfde7cWR4eHvrPf/5z2b7/+te/nFIYAAAAAAAAKi6HQqkePXroxIkTCg0NVY8ePS7Zz2QyKTc311m1AQAAAAAAoIJyKJSyWq0F/g4AAAAAAAAUh9mRTkFBQTp16pQkadCgQUpJSSnVogAAAAAAAFCxORRKZWVlKTk5WZK0cOFCZWRklGpRAAAAAAAAqNgcunyvZcuW6tGjh5o1aybDMDRixAh5e3sX2Hf+/PlOLRAAAAAAAAAVj0Oh1IcffqhXXnlF+/fvlyQlJSVxthQAAAAAAACKzaFQKiwsTC+++KIkKSoqSh988IGCg4NLtTAAAAAAAABUXEVe6PyWW26Rp6dnqRYFAAAAAACAio2FzgEAAAAAAOByLHQOAAAAAAAAlyvyQucmk4mFzgEAAAAAAFAiLHQOAAAAAAAAl3MolLrQwYMHbb9nZGTIy8vLqQUBAAAAAACg4nNoofMLWa1WTZ48WREREfLz89OBAwckSePGjdO8efOcXiAAAAAAAAAqniKHUv/+97/13nvv6eWXX5anp6etvUGDBnr33XedWhwAAAAAAAAqpiKHUu+//77eeecd9e3bV25ubrb2xo0ba/fu3U4tDgAAAAAAABVTkUOpo0ePqlatWvnarVarsrOznVIUAAAAAAAAKrYih1L169fXzz//nK992bJliomJcUpRAAAAAAAAqNiKfPe98ePHq3///jp69KisVqs+++wz7dmzR++//76WL19eGjUCAAAAAACgginymVLdu3fXV199pe+//16+vr4aP368du3apa+++kodO3YsjRoBAAAAAABQwRT5TClJatu2rVatWuXsWgAAAAAAAHCVKFYoJUmbN2/Wrl27JEnXX38960kBAAAAAADAYUUOpeLj43Xvvfdq7dq1CgwMlCQlJibqlltu0SeffKKQkBBn1wgAAAAAAIAKpshrSj322GNKSUnRjh07dObMGZ05c0bbt29XcnKyRowYURo1AgAAAAAAoIIp8plSK1eu1Pfff6969erZ2urXr6/XX39dt99+u1OLAwAAAAAAQMVU5DOlrFarPDw88rV7eHjIarU6pSgAAAAAAABUbEUOpW699VY9/vjjOnbsmK3t6NGjGjVqlG677TanFgcAAAAAAICKqcih1Jw5c5ScnKzIyEhFR0crOjpaUVFRSk5O1muvvVYaNQIAAAAAAKCCKfKaUtWqVdOWLVv0/fffa/fu3ZKkevXqqUOHDk4vDgAAAAAAABVTkUMpSTKZTOrYsaM6duzo7HoAAAAAAABwFXD48r01a9aofv36Sk5OzrctKSlJ119/vX7++WenFgcAAAAAAICKyeFQatasWRoyZIj8/f3zbQsICNDDDz+smTNnOrU4AAAAAAAAVEwOh1J//PGHOnXqdMntt99+uzZv3uyUogAAAAAAAFCxORxKnTx5Uh4eHpfc7u7uroSEBKcUBQAAAAAAgIrN4VAqIiJC27dvv+T2bdu2qWrVqk4pCgAAAAAAABWbw6FUly5dNG7cOGVkZOTbdu7cOU2YMEF33HGHU4sDAAAAAABAxeTuaMfnnntOn332merUqaPhw4erbt26kqTdu3fr9ddfV25urp599tlSKxQAAAAAAAAVh8OhVFhYmH755Rc98sgjGjt2rAzDkCSZTCbFxsbq9ddfV1hYWKkVCgAAAAAAgIrD4VBKkmrUqKGvv/5aZ8+e1b59+2QYhmrXrq3KlSuXVn0AAAAAAACogIoUSuWpXLmybrjhBmfXAgAAAAAAgKuEwwudAwAAAAAAAM5CKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuNwVEUq9/vrrioyMlJeXl1q0aKHffvvtkn3bt28vk8mU76dr164urBgAAAAAAAAlUeah1OLFixUXF6cJEyZoy5Ytaty4sWJjYxUfH19g/88++0zHjx+3/Wzfvl1ubm7q2bOniysHAAAAAABAcZV5KDVz5kwNGTJEAwcOVP369fXWW2/Jx8dH8+fPL7B/UFCQqlSpYvtZtWqVfHx8CKUAAAAAAADKkTINpbKysrR582Z16NDB1mY2m9WhQwdt2LDBoeeYN2+e7r33Xvn6+pZWmQAAAAAAAHAy97Lc+alTp5Sbm6uwsDC79rCwMO3evbvQ8b/99pu2b9+uefPmXbJPZmamMjMzbY+Tk5MlSVarVVartZiVXxmsVqsMwyj3xwEAZYm5FABKjrkUAEquIs2ljh5DmYZSJTVv3jw1bNhQN9544yX7TJ06VRMnTszXnpCQoIyMjNIsr9RZrVYlJSXJMAyZzWV+JSYAlEvMpQBQcsylAFByFWkuTUlJcahfmYZS11xzjdzc3HTy5Em79pMnT6pKlSqXHZuWlqZPPvlEkyZNumy/sWPHKi4uzvY4OTlZ1apVU0hIiPz9/Ytf/BXAarXKZDIpJCSk3H9gAaCsMJcCQMkxlwJAyVWkudTLy8uhfmUaSnl6eqpZs2ZavXq1evToIen8m7B69WoNHz78smOXLl2qzMxM3X///ZftZ7FYZLFY8rWbzeZy/yZLkslkqjDHAgBlhbkUAEqOuRQASq6izKWO1l/ml+/FxcWpf//+at68uW688UbNmjVLaWlpGjhwoCSpX79+ioiI0NSpU+3GzZs3Tz169FBwcHBZlA0AAAAAAIASKPNQqnfv3kpISND48eN14sQJNWnSRCtXrrQtfn7kyJF8CduePXu0bt06fffdd2VRMgAAAAAAAErIZBiGUdZFuFJycrICAgKUlJRUIdaUio+PV2hoaLk/tQ8AygpzKQCUHHMpAJRcRZpLHc1eyvdRAgAAAAAAoFwilAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFyOUAoAAAAAAAAuRygFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlCKUAAAAAAADgcoRSAAAAAAAAcDlCKQAAAAAAALgcoRQAAAAAAABcjlAKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFyOUAoAAAAAAAAuRygFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlCKUAAAAAAADgcoRSAAAAAAAAcDlCKQAAAAAAALgcoRQAAAAAAABcjlAKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXuyJCqddff12RkZHy8vJSixYt9Ntvv122f2JiooYNG6aqVavKYrGoTp06+vrrr11ULQAAAAAAAErKvawLWLx4seLi4vTWW2+pRYsWmjVrlmJjY7Vnzx6Fhobm65+VlaWOHTsqNDRUy5YtU0REhA4fPqzAwEDXFw8AAAAAAIBiKfNQaubMmRoyZIgGDhwoSXrrrbe0YsUKzZ8/X08//XS+/vPnz9eZM2f0yy+/yMPDQ5IUGRnpypIBAAAAAABQQmV6+V5WVpY2b96sDh062NrMZrM6dOigDRs2FDjmP//5j1q2bKlhw4YpLCxMDRo00JQpU5Sbm+uqsgEAAAAAAFBCZXqm1KlTp5Sbm6uwsDC79rCwMO3evbvAMQcOHNCaNWvUt29fff3119q3b58effRRZWdna8KECfn6Z2ZmKjMz0/Y4OTlZkmS1WmW1Wp14NK5ntVplGEa5Pw4AKEvMpQBQcsylAFByFWkudfQYyvzyvaKyWq0KDQ3VO++8Izc3NzVr1kxHjx7VtGnTCgylpk6dqokTJ+ZrT0hIUEZGhitKLjVWq1VJSUkyDENm8xWxZj0AlDvMpQBQcsylAFByFWkuTUlJcahfmYZS11xzjdzc3HTy5Em79pMnT6pKlSoFjqlatao8PDzk5uZma6tXr55OnDihrKwseXp62vUfO3as4uLibI+Tk5NVrVo1hYSEyN/f34lH43pWq1Umk0khISHl/gMLAGWFuRQASo65FABKriLNpV5eXg71K9NQytPTU82aNdPq1avVo0cPSeffhNWrV2v48OEFjmndurU++ugjWa1W25u0d+9eVa1aNV8gJUkWi0UWiyVfu9lsLvdvsiSZTKYKcywAUFaYSwGg5JhLAaDkKspc6mj9ZX6UcXFxmjt3rhYuXKhdu3bpkUceUVpamu1ufP369dPYsWNt/R955BGdOXNGjz/+uPbu3asVK1ZoypQpGjZsWFkdAgAAAAAAAIqozNeU6t27txISEjR+/HidOHFCTZo00cqVK22Lnx85csQuYatWrZq+/fZbjRo1So0aNVJERIQef/xxPfXUU2V1CAAAAAAAACgik2EYRlkX4UrJyckKCAhQUlJShVhTKj4+XqGhoeX+1D4AKCvMpQBQcsylAFByFWkudTR7Kd9HCQAAAAAAgHKJUAoAAAAAAAAuRygFAAAAAAAAlyOUAgAAAAAAgMsRSgEAAAAAAMDlCKUAAAAAAADgcoRSAAAAAAAAcDlCKQAAAAAAALgcoRQAAAAAAABcjlAKAAAAAAAALkcoBQAAAAAAAJcjlAIAAAAAAIDLEUoBAAAAAADA5QilAAAAAAAA4HKEUgAAAAAAAHA597IuAAAAAED58Mr2H/XKjp8K3Jabmys3N7cCt426/maNatCuNEsDAJRDhFIAAAAAHJKcnaGj6UnFGgcAwMUIpQAAAAA4xN/DSxE+AXZthgwdS0+WJIV7+8tkMhU4DgCAixFKAQAAAHDIqAbt8l2Gl5adKf8Pn5Uk7bzzSVWyEEABABzDQucAAAAAAABwOUIpAAAAAMVmGIbt99MZaXaPAQC4HC7fAwAAAFBkiZnn9P6+TXp158+2tujPpiq6UrCG12ujfrWaK9DiXYYVAgCudIRSAAAAAIrk26N71HPNQqXnZOXbdiDltOJ++1LPbflGS2/tr9iIumVQIQCgPODyPQAAAAAO+/boHnVb9a7O5WTLkHTxxXp5bedystVt1bv69uge1xcJACgXCKUAAAAAOCQx85x6rlkow5Cs+eIoe1YZMgyp55qFSsw856IKAQDlCaEUAAAAAIe8v2+T0nOyCg2k8lhlKD0nSx/s31TKlQEAyiNCKQAAAACFMgxDc3atK9bY13au4658AIB8CKUAAAAAFOp0Zrr2p5x28BypfxiS9qec1pnM9NIoCwBQjhFKAQAAAChUanZmicanlHA8AKDiIZQCAAAAUCg/D0uJxlcq4XgAQMVDKAUAAACgUMEWH0VXCpapiONMkqIrBSvI4lMaZQEAyjFCKQAAAACFMplMGl6vTbHGPla/jUymosZZAICKjlAKAAAAgEP61WouH3dPmR08X8osk3zcPfVAdPNSrgwAUB4RSgEAAABwSKDFW0tv7S+TSYUGU2aZZDJJy27tr0CLt4sqBACUJ4RSAAAAABwWG1FXX3UcLG93D5mkfNFUXpu3u4eWdxys2yPqur5IAEC54F7WBQAAAAAoX2Ij6upIr3H6YP8mzd7xsw6mnrFtq1kpWI/Vb6N+tZorwJMzpAAAl0YoBQAAAKDIAi3eeqx+Ww2sdYMCFj0nSTpw9zOqXqkyi5oDABzC5XsAAAAAiu3CACrI4kMgBQBwGKEUAAAAAAAAXI7L9wAAAAA45JXtP+qVHT/ZtRkybL/X//zlAs+UGnX9zRrVoF2p1wcAKF8IpQAAAAA4JDk7Q0fTky65/di55EuOAwDgYoRSAAAAABzi7+GlCJ+AArfl5ubKzc3tkuMAALgYoRQAAAAAh4xq0K7Ay/CsVqvi4+MVGhoqs5llawEAjuG/GAAAAAAAAHA5QikAAAAAAAC4HKEUAAAAAAAAXI5QCgAAAAAAAC5HKAUAAAAAAACXI5QCAAAAAACAyxFKAQAAAAAAwOUIpQAAAAAAAOByhFIAAAAAAABwOUIpAAAAAAAAuByhFAAAAAAAAFzOvawLQOFm/rhfr/x0oMBt1txcmd3cCtw26uaaimsXXZqlAUC5kbRllpK2zM6/wZByrbn6n9lNMuXfHND0cQU0HVnq9QEAAABXG0KpciA5I0dHkzIu0yP7kuMAAOdZM5OVm3r0kttzLzMOAAAAgPMRSpUD/l7uigjwsmszDEPHkjMlSeH+FplM+f+87+/F2wsAecwWf7n5Rdg3GoZy045Jktx8q0qm/Fe1my3+rigPAAAAuOqYDMMwyrqI119/XdOmTdOJEyfUuHFjvfbaa7rxxhsL7Pvee+9p4MCBdm0Wi0UZGZc7k+gfycnJCggIUFJSkvz9y+8/NFIzsuX/3EpJ0oGxt6hGkG+BwRQA4NKs2Wk6/HplSVK1R07L3VKpjCsCgPLJarUqPj5eoaGhMptZthYAiqMizaWOZi9lfpSLFy9WXFycJkyYoC1btqhx48aKjY1VfHz8Jcf4+/vr+PHjtp/Dhw+7sOKylXguW7N/PqDGM3+ytdWc+oPqvLhGs38+oMRzBV/KBwAAAAAAcCUp81Bq5syZGjJkiAYOHKj69evrrbfeko+Pj+bPn3/JMSaTSVWqVLH9hIWFubDisvPtnnhVm7xKcV/u0KEz6XbbDpxOV9yXO1Rt8ip9u+fSgR4AAAAAAMCVoEwXHcrKytLmzZs1duxYW5vZbFaHDh20YcOGS45LTU1VjRo1ZLVa1bRpU02ZMkXXX399gX0zMzOVmZlpe5ycfH7B2j9n/CI/L9/L1ucTUUnR9zeya9v/4TalH00p9NjCWldXaJtqtse5mTnaOevXQsdJUs0HGso3/J/T25J2n9YnS//UYylpMiQVdL1lXtu57Fzd8e6v+s+gGxVbN0RHv9mnM9tOFrrPgOuCVb37dXZtu9/YqOyUrELHXtuplio3/icYzDiVpr/mbS10nCRd90hzefhbbI9PbTym42sOFjrO6xof1X4wxq7t0JIdSjmYWOjYa5qHq+ptUXZtf7603qF6I3tdr0pRgbbHKQcTdWjJDofGNnyqtd3j46sP6tSmY4WOqxQVqMhe9p/vv+b9roxT6ZcY8Y+qt0bpmhvCbY+zkzO1+81NDtVb+8Em8rrmn+/I2T9O6n8r9xU6zqOSp6579Aa7tiNf7lbS7tOFjg1qFKaIzrXs2nbM+q+smZdagvof1btfp4Drgm2P044l68AHfxY6TpLqj2whN8s/02H8ur91cv2RQsddSXPEkS93FzrObHHT9SNvsmtjjkiXv/f533bM/EZSFV3TPII5wgHMEcwRF6q4c0R+/H9EwXOEIUPn/LJ1InWvTDIxRzBH2GGOYI64GHNEwXPEhXOpb4R/uZ4jrFarQ/sr01Dq1KlTys3NzXemU1hYmHbvLvig69atq/nz56tRo0ZKSkrS9OnT1apVK+3YsUPXXnttvv5Tp07VxIkT87WnmtNlmC+/BlN2jvJdRpiUk6IMc+FfUI/0s1L8PxOgNStXKQ6Mk6RTZ04rzf2fNbKOnTqt0ZcJpC5kNc6f/tbz/U3a/HBjZWQmObRfa6aHvC461mQjTTnmwi8HPJ16Rtnx/7yWmYnnHD7W+NMJ8sjwtD1OTD3r0NjM/7/W9kKJ2clKd2CsOSNRbheNdfi9STylc/H//MczLTHZ8WO9aJ9nMhIdGpubbc7/ObSmKstc+DpqltSzssb/8zXPTstyuN6EM6dksabZHiennnForLuRna/es5nJSnNkv5lJ8rj4c6g0GebCJ7TTyaeVGf/Pf1DOnUlz/FgTEmT2dLM9PpPu2OfwSpkjUpMd+yyZlP+zdOYqnSNMSpWP2xr5mf9j2+7vfb9yrFWUcuJemf/3oEyeAbZtzBH5MUcwR1yoos0Rl8P/RxS8T0OGss1W5ZizZZKJOYI5wg5zBHPExZgjCp9Lc3JM5XqOSEkpPDyTyuHd91q2bKmWLVvaHrdq1Ur16tXT22+/rcmTJ+frP3bsWMXFxdkeJycnq1q1avKz+sjP6nPZffm4V1JoaKhdW4r7CXk4EPgF+VS2G5ubmaOEQvaX55qgYPmG/pNMvvdzvDJVeCCVxyrpXLZVK49k6i5LgGQtfEIJsPjnO9YzpsPKthb+14tgvyBVvmBshjlNidb/OVRraHCI3V8vzIdzlGlNLHScl9knX73pHglycyCNDfIKzDf2pKPvTeA1qhQaaHuckuapVGvhfx2SlG+fuV5psloL/3JX8sj/3iSZjyrDWvjVt0F+lXXNBWOzkzN1xurYGmwhQdfY/fXC47ihc9YzhY7zMHnmqzfDckZma+H/0xFkCcg3NkG+sloL/+tFsH+wAkIv+OtFTrKSrYX/dUiSQkJC7P56IZ9MZVsLn0SvlDnCcsZNadaEQseZ5Zav3mxL8lU3R7iZN8rH83lJmfnGuZlOyv30bGV+NVchXT6Rd42OkpgjCsIcwRxxoYo0RxSG/48oeI4wZOicNVveVg+ZZGKOYI6wwxzBHHEx5oiC54gL51Lfcj5HeHl5ObS/Mr37XlZWlnx8fLRs2TL16NHD1t6/f38lJibqyy+/dOh5evbsKXd3d3388ceF9i1vd98zDEN1XlyjA6fTHQ6lJMkkqWawj/Y+fSt35QOA/5d+6Dud/LK7ZBg6H+FfilkymRTW/Uv5RN7uqvIAoNyqSHeMAoCyUpHm0nJx9z1PT081a9ZMq1evtrVZrVatXr3a7myoy8nNzdWff/6pqlWrllaZZep0epb2FzGQks6fVbX/dLrOpHM3PgCQpNyMRMWv6O1AIKXz2w1D8St6Kzcj0QXVAQAAAFefMo/e4uLiNHfuXC1cuFC7du3SI488orS0NA0cOFCS1K9fP7uF0CdNmqTvvvtOBw4c0JYtW3T//ffr8OHDGjx4cFkdQqlKdWDRtctJycxxUiUAUL6l7vpARna6Cg+k8lhlZKcrddeHpVkWAAAAcNUq8zWlevfurYSEBI0fP14nTpxQkyZNtHLlStvi50eOHLE7be3s2bMaMmSITpw4ocqVK6tZs2b65ZdfVL9+/bI6hFLlZ3ErvNNlVLKU+VsMAGXOMAwlb329WGOTt86Rf5NhXAoNAAAAOFmZrilVFlhTCgCuPrnnTunI2+GFd7yE6g8fl5t3cOEdAeAqVZHWQQGAslKR5tJysaYUCmcymTS8TVSxxj7WJopACgAkWbNSSzjesVvaAgAAAHAcoVQ50L95Nfl4usnsYL5kNkk+nm7q17xa6RYGAOWE2dOvhOMrOakSAAAAAHkIpcqBQG8PLevfXCaZCg2mzCbJJJM+7d9cgd4erikQAK5wZq9guQfU1PmLm4vCJPf/a+/O46qq9v+Pvw6DjAIqoKIIzhOOqGCC13nIgRyyvNe0gbRUSC0zs7yVU5bpz6kSrbQ0S0sc8mqmZIpWKo5pyUUTJ9QQBRRBOOf8/uB79hWnurdkkPfzn3Kfs89Zm8d5fPZen/VZa3nWwM65/L1oloiIiIhIqaakVAnRta4vX0W2wsXRHhO3dqtsx1wc7Vkf2YoudX0Lv5EiIsWUyWTCo+mI/+lcj6YjNRVaREREROQeUFKqBOla15dTr3ZmVkRDAsu7FnitRgVXZkU05PSrnZWQEhG5Dff6j2FydOWP3/rsMDm64l5/0L1sloiIiIhIqeVQ1A2Q/46XiyPR4TV4sqU/Hq9sBODXl9tTrZybRvJFRO7C3tkL3x6fc35NBFjtAMtd3m0HJhO+PVdg7+xVSC0UERERESldlJQqAWZ+d4xZ244XOGa1Wo3/bzNv520TUqPb1mDM32re8/aJiJQUroFdqBixhgvrH8Gam/V/R603vCM/lpocXfDtuQLXgM6F3kYRERERkdJCSakSICM7jzPp2Xd8/WxGzh3PExGRglwDu+D/1K9c+XkpGfvmkpfxq/Gag2d1PJqOpGyDx7Bz8izCVoqIiIiI3P+UlCoBPJwdqOLpfNvXLGYzdvb2dzxPRERuZe/shWezkbg3fJyT7+bvrOf3xL8p41FNU6FFRERERAqJshYlwJi/1bztNDyLxcKFCxfw9fXFzk5r1ouI/LduTEDZO5dXQkpEREREpBApkyEiIiIiIiIiIoVOSSkRERERERERESl0SkqJiIiIiIiIiEih05pSIiJSKqTv/X+k751d8KDVavzv2Y+DwHTrWI1n8+fwbD7qHrdORERERKT0UVJKRERKBUtOBuYrZ+74uvlqyh3PExERERGRv56SUiIiUirYOXlg717l1hesYLaYsbezh9tsvmfn5HHvGyciIiIiUgopKSUiIqWCZ/NRt52GZ7FYuHDhAr6+vtjZaalFEREREZHCoqdvEREREREREREpdEpKiYiIiIiIiIhIoVNSSkRERERERERECp2SUiIiIiIiIiIiUuiUlBIRERERERERkUKnpJSIiIiIiIiIiBQ6JaVERERERERERKTQKSklIiIiIiIiIiKFTkkpEREREREREREpdEpKiYiIiIiIiIhIoVNSSkRERERERERECp2SUiIiIiIiIiIiUuiUlBIRERERERERkUKnpJSIiIiIiIiIiBQ6JaVERERERERERKTQKSklIiIiIiIiIiKFTkkpEREREREREREpdEpKiYiIiIiIiIhIoVNSSkRERERERERECp1DUTegsFmtVgAyMjKKuCV/nsViITMzE2dnZ+zslF8UEflfKJaKiPx5iqUiIn/e/RRLbTkXWw7mTkpdUiozMxMAf3//Im6JiIiIiIiIiMj9KzMzE09Pzzu+brL+XtrqPmOxWDh79ixly5bFZDIVdXP+lIyMDPz9/Tl16hQeHh5F3RwRkRJJsVRE5M9TLBUR+fPup1hqtVrJzMzEz8/vrlVfpa5Sys7OjqpVqxZ1M/5SHh4eJf4HKyJS1BRLRUT+PMVSEZE/736JpXerkLIp2ZMURURERERERESkRFJSSkRERERERERECp2SUiWYk5MT//znP3FycirqpoiIlFiKpSIif55iqYjIn1caY2mpW+hcRERERERERESKniqlRERERERERESk0CkpJSIiIiIiIiIihU5JKRERERERuSuLxVLUTRARkfuQklIlwIwZM+jVq1dRN0NEpET78MMPjViq5RRFRP47dnb/6TYohoqIyF9FSakSwNvbm2+//ZbffvutqJsiIlJi2dvbs379enJycjCZTEXdHBGREmXt2rU888wzWK1WxVAREfnLKClVzFgsllvKozt37sz169c5dOhQEbVKRKRksVqtt8TS1q1b4+LiQnx8fBG1SkSk5DCbzZjNZuPfbm5uxMTEkJqaSnx8PAsWLCjC1omIlHwnT55k9uzZHD16FCi9VahKShUTts6TnZ1dgfJogCpVqlC7dm02b95cFE0TESkxbLHUZDLdEksrV65Mw4YNWbt2LVB6b/wiIn+Evb099vb2xr+/++47TCYTVatWpVevXuzfv58rV64UYQtFREoG2zPn8ePHOX36tHH8ww8/ZM6cOVSvXh2LxVJqq1CVlCoiN448wX/m6X/33XfMmDGDbdu2kZeXZ7zepUsXvvnmmwLHRERKs9sllWyx9PDhw3z00Ufs3LnTeM3V1ZX27duzadOmQmujiEhxdbsYajuWnp7O4sWLGThwIKNGjQLyk/4eHh48+uijXLp0iffeew93d/fCbLKISIly9epVIH+w9NChQ/To0YNevXpx+vRp8vLy+PLLL3nuuecoU6bMLYOppUnpvfIiZht5OnfuHAAffPABtWrVYtCgQaxevZqhQ4cSGRlpvL979+4cOnSIs2fPFkl7RUSKC1un6ebRpNzcXFasWEFQUBBt27bl3Xff5amnniI6OhrIj7vt27cnMTGRixcvltrRKBEpvSwWizEwersYaDKZ+Pnnn+nduzdvv/023t7eVK5cmT179jB58mSee+459u/fT05OTmE3XUSkREhKSuKZZ56hWrVqdOjQgcmTJ2M2m2nUqBEJCQn4+voSHh7OO++8Q2pqqjY0Q0mpe8JqtZKXl3fbrXMzMjLYs2cPmzZtwtHRkYEDBwL5WdRJkyZx6tQp4uPjWbVqFR9//DE//PADACEhIdjZ2bF79+5CvRYRkaJktVoLVJbaFtjNzs4mLi6OhISEAu//9ddfGTNmDKdPn2b37t0sWLCABQsWsG/fPgDq16+Ph4cHW7ZsKdTrEBEpChaLpUBFlJ2dnTEwum/fPuLj48nKyjJez8rKYurUqWRnZ7N9+3bmzp3L888/T4sWLQBo164dP//8M2fOnCncCxERKaYOHjzIiy++yPr167FarUycOJFjx44xY8YMBg8ezBtvvMH48eP57bffcHV15csvv6RPnz5MmTKF0NBQypUrB5TuZSWUlLoHTCYTDg4O2NnZFehMZWVl8eKLL9K9e3cWLVrE6tWr+fTTTwHo378/AwcO5PTp08yePZvx48cDsGHDBrKysvD09KRZs2Z8/fXXRXJNIiKF4eZkvslkKlBZajKZ+Oabb6hfvz5DhgyhX79+fPHFF+Tl5eHo6Ej//v158sknyc3N5fPPP+eDDz4gNzeXVatWAVCxYkVatGjB6tWrgdL9ACAi95+bY5qdnV2BiqhTp04RFRVFxYoV6dy5M9HR0fTt25fU1FQALl++zIYNG3jhhRcoX748AA4ODsb5wcHBODs7G4OmiqEiUhplZ2fTu3dvPv30U15//XUOHDiAl5cXGzZsYOvWrYwcOZIBAwYwYsQIYmJi2Lx5s7F8hLu7O4MGDeLKlSscOnSIKVOmALevXi0tlJT6H928I8mNkpOTmTBhAs2aNaNLly4sWrSInJwcXF1dadKkCbm5udSqVYsePXpQuXJlrFYrfn5+LFmyhC5durBixQrCw8Pp06cP69atIzMzE4CuXbuybds2rl27VpiXKiJyT90YS2+cT3/y5EkuXrzIU089hbu7O3/7299YtGgRS5Ys4b333iM5OZnmzZvz7rvvcvDgQQBq1qzJ6tWrCQ0NZcqUKXh6etK3b1++/PJLrl+/jrOzM506deK7774DSvcDgIjcf26Oab/99hvDhw9nx44dAKSlpZGamsry5cs5f/48q1atIi0tjenTp3P9+nU8PT3JyMjAzc0NuHWgoGzZsoSGhvLVV1/d9vtERO43Fy9eZOXKlezatYvc3FwAnJ2dOX36NIMGDaJly5Z8/fXXtGnThl9++YXy5cvzwAMPGOe3bduWGjVqsH37duNYQkIC1atXZ8GCBSxatIjHH3/8trOsSgslpf4LN3acbt6RBPJHi6xWK6+++io7d+5k+PDh9OnTh1dffdXIgDZv3hwnJyf8/PwAyMvLw2Qy8dNPPzF58mQefvhhvvnmG1544QUeeughjh07xqlTpwB48MEHSUxMJCkpqZCuWETkr2WLkzeyxdJr166xbNkyXn75ZdasWUOLFi14/PHHqVq1KvHx8bRt25bo6GicnZ3p1q0bdnZ2PP/88+Tm5rJ161Ygfx7/2LFj6dOnD9u2bWPOnDn069eP48eP8+uvvwL5DwcpKSkkJiYW6rWLiPwVblwX6mY7duxg165dxr/j4uJYvnw5devWBSAwMJAZM2bQoUMHzp07x969e8nIyGDz5s0kJibi5uZGUFCQMaJvZ2dnxGzbTnt///vfWbduHS+++CKDBw9mzJgx9/JyRUQK3Z49e1i/fj3z58+nUaNGjB49mr///e+MHz/eiL99+/albNmyBdaEaty4MSdOnODatWtG7KxRowZOTk5cvnzZOHfu3LkMGzaMjh07snHjRnbs2EFISAhHjhwp/IstBpSUugOr1XpLttLWccrJyWHlypUMGDCAcePGcfz4cSB/tOi9994jKSmJ9evX8/TTTzNy5EiGDh3KwoUL2bx5M7Vr16ZFixYFHhgg/6b/66+/0r9/f1xdXUlLS2Px4sVkZmaSkJCAxWKhcePG1KtXzyixFhEpKW5cWPfmkfVPPvmE8ePHM2XKFObOnYubmxsVKlSgTZs27N69m379+tG0aVNeffVVWrduzW+//WacW69ePXx8fIw1owCOHTvG448/jpeXF1lZWXzyySdcv37d2ImvevXq1KhRg6NHjxbClYuI/LVuXBfqxIkTZGdnAxhrmfTs2ZO1a9cCMGfOHIYNG4a3tzdWqxVPT0/c3NwYOnQoLVq0YOrUqdSrVw+z2WzEyEGDBrFs2TLWrFljrOP3448/snTpUiA/KTV16lR2796Ns7MzjzzySBH8FURE/jopKSkFdrmPjY2lV69erFu3jtjYWI4dO8bzzz/PypUreeuttwAIDw8vsCYfQFhYGHl5ecTFxRV43k1ISKBu3brY29vz7bffkpqaSps2bYD8taO/+uorFi5cSIMGDUrltGiH339L6WKxWIz59zd3nL766ivi4uJwd3dn3759VK5cmY0bNxIfH8/ixYupXbs2x48fx8HBgbi4OGbPns3Bgwe5fv06PXr0oGLFinh4eNCoUSM2b94M/GeefoMGDfD29mbkyJF0796d3bt3M3DgQHJyckhPT8disVCmTJlSmz0VkZLN1oFKSEjg3LlzhIaGUqFCBQCuX7/O9OnTadGiBZ988gl169bl2rVr1KxZk927d9O4cWMA/Pz8CA0N5euvvyYvLw8HBwcqVKhgvO/s2bPUrFkTPz8/hg8fTv/+/UlISCAkJISsrCwjCVWpUiVVnIpIsbF3717MZjMtW7Y0kkB3k5SUxOTJk1m7di1eXl40bdqUqKgo2rdvz5YtW3jhhRcYO3YsixcvxmKx0KlTJ+A/G0V88MEHbN++nVWrVtG6dWsyMzMJCgoypkEPGzaMI0eO8NRTT9GpUyeSkpI4d+4cI0aM4Pr165QpU4aoqCiioqLu+d9GROSvZIuDqampXLp0iVmzZvHRRx9Rs2ZNunbtymuvvUbZsmUZMmQI06ZNo3LlyoSEhADw7LPPcuTIEf71r38RFRVFaGgozs7OHDhwgIYNG2KxWHB2dubRRx9lxowZXL16lcGDB7N8+XIAIwl19OhRHn/8cYKDg4122apZoXROiy7VlVJms/m2C0JCfsnewoULC4y+5+TksHjxYmJjY5k0aRIxMTEsXLiQzMxMYmNjAfD29mbHjh289tprBAcHs2LFCpKTk1m6dCmNGjXCwcGBBg0acO3aNfbu3Qvkd8gAVq9eTZ06dfjss8/w9/cnIiKC+Ph4XnzxxQKLTN6pZFtEpKjcrrr0RqtXr6ZKlSp07tyZCRMm0LNnT2NUvnv37lSoUIGGDRsaN2UXFxdCQkK4dOkS58+fB/KT+PXr1+fKlSsFdt1r2rQp169f54cffsBkMrFkyRJ8fX2ZNm0a165dY+DAgcTFxfHmm28WaJNiqYgUtatXrzJo0CAWLlwI3L0zYoux8+fP5/z583z11Vds3boVLy8vpk2bxrZt2wCYOnUq06ZNY8uWLZw/f55GjRoB+c+4WVlZ7Nixg1atWtG6dWsAtm7dSkZGBj/99BMXL17Ezc2NmJgYlixZQpUqVRg8eDA7d+5k/PjxlClTxmjP3XabFhEpjkwmE4cPH6Z27dq8/PLL5OXlsWnTJl544QUWL15sHKtTpw5ly5bFz8+vQAVVq1atuHr1KgcOHMDJyYkWLVoYG5HZ4vfUqVPp0aMH7777LoGBgYwfP56hQ4caAwTPPPMMU6dOxdnZuUDbSmOFlE2pqpSybS1uZ2dXoPT5Rvv37+fxxx/nzJkz1KlTh3nz5tG2bVvmzp1L27ZtCQgIoEqVKjRt2hTI/2HWr1+fPXv2cP36dfz9/fHx8eHjjz+mQYMGxueePHmS5ORkwsPDqVu3Lrm5uaxbt47mzZsb7QgNDSUkJOSWBxKz2Vygrbdrt4hIUbJVl6anp3Py5EkaNmxoJPnPnj3LSy+9RGRkJGPHjuXEiROMHTuWZ599lvj4ePz8/AgICMDR0ZGMjAw8PDwAqF27Nt7e3mzcuJEhQ4YA+SNJvr6+bN261Ri5qlevHtnZ2cYaUR07diQ8PLxA5wn+Uwlro1gqIoUhIyODpUuXsnTpUubOnUtwcLCxYY6bmxsPPPAAFy9e5PTp01StWvWO1VJ2dnZs376d+Ph4NmzYgLe3N9nZ2TRu3JgVK1YQGxtL27ZtKVOmDL1798be3p7s7GyefvppPvjgA3x8fHB1dcXX15ctW7awdOlSnJycWLduHR06dODSpUskJydToUIF7Ozs6NGjBz169Ljjddl2mxYRKY7i4uKIjY3l2rVrDBo0iHbt2gHQsGFDqlatyr/+9S82b95M69atCQ8Px2q1MnnyZCIiIujUqRMhISEcOHCA9PR0o7q/UqVK5ObmGkUlvXr1Yu7cuVy5cgV3d3cgv7J/+vTpDBw4EGdnZ+rVq3dL22w5iRtjfWmskLIpVZVStpunnZ0dOTk5nD59mhdeeMHY3e7q1auMGTOGTp06kZKSwo4dO/j4449ZsGABW7duxcfHhzp16uDk5MSFCxeMz23UqBFnz54lMTGRXr16ERAQQHR0NHFxcVy7do2DBw8yZcoUVq9ejcVioUaNGowcOZL27dsDBTtGJpMJi8VCXl6ekS1Vx0lEirNdu3aRkpJCREQElStXplevXowYMcKoRDp48CBpaWkMHjwYd3d3goKCePPNNzl37pxRZdqyZUuSkpIKrBfl5+dHkyZN2Lhxo3EsICAADw8PY+cnyK+UWr9+PS+99JJxrEyZMreM4t+YkBIRKSyTJk1i/vz5hIWFERAQAOQ/29kS582aNeO3337j559/Bu4+Wn79+nUSEhKYOnUqtWrVwsfHh4ULFzJy5EgjeQ/wxRdfUKlSJT755BNSU1OJiIgwdn4aP3483bp146WXXmLMmDGEhIQwf/58tmzZQvPmzQt8n21AtzSP4ItIyXLhwgUeffRRhg0bxrlz5zCbzXTt2pXFixeTk5MDQHBwMAEBAVSsWNE4LywsDH9/f+O5s1+/fmzZsqXArnl79+7lzJkzhIaGAtCuXTuSk5NJSUkp0Aar1UrTpk2NhNTNcdTe3r5UJ6FuVqqe0LOzs42beNeuXZk3bx4zZ8401hm5cOECqampvPTSS1y5coX333+fGTNmkJeXZ2wfHhwczIULFzh27JjxuS1btsRsNrNr1y48PDx49913cXBwYNSoUdSrV4/Q0FDOnTtHjx49MJlM+Pj4EB0dTdu2bW/bTjs7OxwcHPRDFZFi7/z584SGhvL0008biaXXX3+dxYsXs2jRIiD/Bm6rELWpUaMGLVq0MJJLXbp04cSJE8YOeQDlypUjODiY1atXG8d8fX0ZO3Ys8+bNM47Z29tTuXLlW9p240CEiEhRSE9PJz4+noiICN566y28vb0B+Prrr+nUqRMREREcP36czMxMIyl1t5jl4uJCuXLlSExM5M033+TQoUMcPnyYN998k6ZNmxqDAZ9//jnNmzenY8eOfP7559StW5d+/frx1VdfERAQwNtvv80vv/zCqVOnGDFiBJUqVcJkMt2SfDKZTOo8iUiJkpubS4cOHdiyZQsrV67ko48+Yvz48SxevNhYPicsLAyz2czp06eN8ypVqkSlSpU4c+YMABEREeTk5BAVFcWsWbN4+eWXef/995kwYQIuLi4ANG/enCtXrlC7du0CbbDFzBuLTBRH76xUPamvXbuWDz/8kJEjR/LOO+9w4sQJAGNdkx9//JGsrCyCgoIICAhgwYIF+Pr6snHjRkaPHg1A69atuXr1aoEFx5s1a4a9vT07duwA8hNXGzduZN68eXz22WdkZWWxZs0aOnToUODHqDn4IlKcmM3mAnHpj4yMV6xYkc6dO/Pjjz/So0cP/Pz8GDJkCI899hhLliwB8sukMzIySE5ONs5zc3OjatWqxhbjHTp0wGq1Foitjo6OhIeH8/TTT3P16lXj+N/+9jeaNGnyp69XRORe8/T0pEaNGsTFxfHYY4/x6KOPkpaWRlRUFIGBgTz11FOcPn2aw4cPc+TIkQLJ+9upWLEi5cqVo2XLlvTv35/AwEAATp8+zbJlyzhx4gT79u0jPj6egQMHAuDv78+CBQtISEigZ8+eADg7O+Pu7n5Ldb46TSJSXGVnZ5ORkQHc/RnVx8eHyMhIypYty+zZs+nevTvTp0/n6NGj7N+/H8hf6iE9PZ24uDjjPKvVys6dO2nTpg0Wi4VKlSpRrVo1unbtysmTJ9m/fz/jxo3j2WefLfB9rq6ud2yPYuofUyIngpvNZkwmkzGS9Hs7ldjWEZk3bx7BwcGMGjUKgJiYGDIzM1m3bh3R0dH4+/vj6OhI7969mTVrFi4uLsZ3XLt2DcjPhpYtW5affvrJ+HxfX19at25NYGAgubm5ODo6YrVaC1RC2UaubpyKp9F7ESkqtoXJb7deXVZWFqmpqVSrVu0PfVaLFi1ITk6mTp06xrE+ffqwdOlSTp06RXh4OI6OjqxatYquXbsCcO7cOVavXs0rr7yC2WzG09OTcuXKkZSUxNWrV3FzcwPyk1UdOnS4bft1oxeR4i46Oprly5djb29PnTp1GDZsGNOmTcPFxYU33ngDPz8/evfuTa9evThy5AjJycnUqlXrjjGuZs2aREVFMWXKFC5cuEBkZCRJSUl88cUXWK1W2rRpg729PcOHD6djx47GeY6Ojvj7+9/yebZ1VkVEiqOtW7eya9cuateuzWuvvcaDDz7ItGnTMJvNd1zTrkyZMmRkZPDkk09y8eJF2rVrx4QJE4iKiuLIkSPk5eVRvXp1AgMDeeeddyhXrhzdunVjxYoVODs7ExYWZsTFNm3akJyczKpVqyhbtuwd26ln0j+n2N+FbHPZb2Rvb2/sIHLy5Mnf/RHY2dlx4cKFAjdki8WCh4cHAwYMYPfu3aSnpxMUFESNGjU4ceIEbm5uxo/x+PHjREZGkpKSgpubGz4+Ply4cIGLFy8a3/H2228zYsQIHB0dgduX7GltKBEparZ4apuSYWO1Wvn0009p2LAhVapU4R//+AcTJ0404tzdRqS6dOnCr7/+WmCtvbZt25KXl8cPP/xA+fLlGTduHJ9++in/+Mc/WLBgAc8++ywNGjSgT58+Rju+/vpr5syZYySkbCwWyy2Vpbr5i0hJEB0dzeLFi6lSpQrDhw+nffv27N+/n+DgYPz8/IyY/PDDD5ORkfGH1pUaMWIEU6dO5fjx4/Tp04fRo0fj7e3NuHHjCAwMpHHjxrz++usFdnZSzBSR4swW886cOWMs5fDUU08RERFBQkICX375JYcOHTIq6n8vmR4TE8PBgweZOXMmb7zxBmFhYdjb2xvJf8ifwufu7s7+/fvp06cPMTExvPzyyzRq1Mh47nzkkUfYu3evMc3v96pZ5X9TbCulbDvO3a7jtHz5cqZMmcLZs2cJCgqiffv2PPfcc1SoUOGOI0tlypTBz8/PWITM9kMOCgoiMzOTPXv20LFjRyZOnEifPn0ICQnhwQcfJDExkX379lG7dm1j+sjChQspV67cHdt8Iz0EiEhRsSVzbpzHbotRGRkZzJ8/n2PHjvHPf/6TtLQ0Zs6cSWRkJA899BD79u1j3rx5XLp0iblz5971e5o3b06ZMmX44YcfqF69OlarFTc3N5o3b87WrVvp27cvDz/8MC4uLnzxxRfMnj2bsLAwJk+eTNWqVY3P8fHxue3naxRfREqqWrVqYTab8fPzIzY2ltDQUOrVq8euXbsAjBgdHh7O6NGj+eWXX+jVq9dd4569vT2RkZFGUt/Ly+uW91gsFmNXVBGR4iotLQ2z2YyPjw///ve/GTx4MP379yc8PJxPPvmEzZs307ZtW06dOsW1a9f49ttvycvLu2OVlG2G1MmTJ6lduzaenp4AfPbZZ6SkpJCens7u3bupWbMmoaGhbNy4kYcffpiYmBijuORGnTt35tKlSxw+fJj69evf9j3y5xWLJ/2b57NDwY7TtGnTiIyM5NSpU0bGMzIykr179zJ69Gji4+N57bXX7vodXl5e1KtXj/379xcYzf/++++xWq3s3r0bgJCQEDZt2sTAgQPZsWMHrq6uzJ49m9WrVxvl1LaE1M2jWKqEEpHi5MZNEy5fvkxaWhobN26kZcuWLFu2jM2bN9OkSROsViufffYZrVq1YvTo0VSvXp2aNWvi5OTE/PnzuXz58l07NmXLliUkJITNmzcXON66dWvWrl1rrBvVs2dPFi1axJEjR4iJiaFRo0b39PpFRIqDSpUq0bhxY2Pt0a5du5KQkEBycrLRwbl27Ro5OTkcPnzY2BX6RhaLBbPZbCxhYTabqVChAl5eXrfsNArcstW4iEhxcOMMqG+++QZvb2/efPNNAKpVq0ZKSgpBQUHExcXRvHlzY5c7f39/IiMjMZlMbN26Fbj9+sy2/vmDDz7IqVOnGDx4ML169SImJoaJEyfSv39/6tatC0CnTp24evUqJ0+eNJbfycvLMz7LYrHg4uLCtm3beOihh+7Fn0P+T7FISt3rjpPtx/nQQw/h6OjIwIED+eGHH9izZw+HDx8mODjY2AHKYrEQFBTEqFGj2LRpEwsXLqRz587GjiQ3fr5u9iJSnJ04cYLo6GiqVq1KaGgoH330EVWqVCEhIYFZs2bx3nvvERUVRbVq1YiNjeXSpUs8+OCDVKhQgc6dO1OmTBnmzZt3x9GoG7Vr147ly5dz5coVIzZOmDCB+Ph4Y5QKwMHB4bYdKBGR+5WnpydNmzYlOTmZ1NRUevbsSc2aNXnmmWfYuXMnly9fZv78+fj4+PDTTz+RmJgI5E8TsXXg7OzsjKUgcnNzsbe3N2KodhoVkZJgw4YN+Pj4GLOP0tPTCQwMZNmyZSQmJuLk5ERWVhYmk4msrCy8vLwK7I5Xq1Yt/P392bBhA3D7pJStSKRr1658/vnnNG7cmPr16/P2228ba/o1a9YMq9WKl5cXnp6exMXFkZaWZsRSGzs7O6xWK2FhYX/oWVj+d8Xi7nWvO062DlLDhg2ZP38+Tk5OdOvWjfbt2+Pv78+oUaM4dOgQUHCayM0dJyWhRKSkMJvNTJw4kaNHjzJz5kxWrFhBkyZNCAgIIDAwkPbt2xdYmLxx48bExsYSHBzMunXrSEpKYs2aNQwfPhx3d/ff/b5HH32UV155pUAc9vb2JiAg4Jb3qgMlIqVN3bp1cXNz49tvvwXg/fffx9HRkQEDBlChQgVyc3P54IMPePfdd2nWrBmQvzi5rYP1/fffM27cOBo2bGhMe1YMFZGSpFOnTri7uzNjxgwAzp8/zyOPPEL16tWJiYkhJSWFRo0acfLkScLCwjh79ixJSUnG+eXKlSMvL8+oOr3TLCVbMj8oKIgFCxbw1ltvERwcDPxnaQvbe6Kjo+nXr98dFzFX/79wFHnKz9ZxOn/+PDNnzqRevXpcuHDhdztOY8eO5ZVXXiEoKAgPD48//H1t2rRh3bp1HDt2zPjcJ554gho1apCenl5gRP/mbKmISEmxbNkyNm7cyJIlS+jevXuBSs9mzZpx6tQpzGazMcUjNDSUnTt38vTTTxfYde/7778nMTGRIUOG3PIdZrMZq9WKg4MDtWrV4uWXXy606xMRKUkCAwMpX748K1eu5OGHH6Z9+/YEBweze/du6tatW2B9Pcifzjdv3jzWrVvHoUOHcHFxoVWrVkRHRxMeHl5EVyEi8r9zdHRk1KhRrFq1ip49e+Lj48OmTZuYMWMGb7/9NitWrMDX15esrCxatmyJu7s7H330EWFhYbi6urJjxw4uX77MiRMnSEtLo3z58lit1gI7St+8HjVAXl6esdOoLZlv+++TTz5Z6H8HuVWRZ1wKu+MEkJycjLOzMydOnGDNmjV8//33vPbaa3h6emqbcRG5L+zbt48aNWrQvXt3ID/Jblv8sXv37kyYMIGMjAxjjbwnnniCZcuW0bNnT15//XX8/PzYtGkT27ZtIyIigtzcXBwcHO64CcWVK1dwd3e/6+KTIiKlVaVKlRgyZAiurq7GMQ8PDzp27Aj8Z6kJq9Vq7Bp94MAB2rVrxzvvvENQUBAuLi5F0nYRkb/KgAED2LdvH9OmTWPSpEn89NNPtGnThu3bt7Np0yaOHz9O9+7dKVeuHFOnTqVv375ERETg6OjI5cuXGTt2LDExMcTHx9O7d29jQwdbkuncuXPExsby888/M3ToUIKCgvRcWgIUed3vnTpOAN27dychIYGMjAwjUfTEE09QsWJFevbsSWxsLD/++COTJk1i4sSJZGZmkpuba0y7s93g7e3tC/wYL168yDPPPEOjRo1YsGABUVFR9OrVy/h+EZGSLjAw0JiHf+OaJAAdOnQgNTXV2Hoc8kuiP/30U0JDQ5kxYwZdu3Y1diQZNGgQjo6ORvWoyWQiJSWFRYsW0aNHD6pUqcKcOXMAdOMXEbkNFxcXhg0bxmOPPVbguO1Z9eaOVUBAAEuXLuWNN96gZcuWSkiJyH2hatWqDB06lHXr1vHFF1/QoEED0tLSiIyMJDMzk6NHjxrvbdeuHVu2bKFVq1YEBgby4YcfEhwcTF5eHjk5OQBkZWWxdu1annjiCapXr069evV4//33KVu2LBUrViyqy5T/UpH3HgIDA1m5ciWAMQJ/u47TAw88APyn4/TOO+8wY8YMDh8+TMOGDRkyZAgDBgwwdjGxdYxSUlJYv349sbGx7Nu3j9GjRzN27FimTp1KzZo17zh/VESkJGvfvj2jR49m+/btBaZ6HD9+nAoVKlCvXj02bNhgxFar1UrdunVZsGABZ8+epUqVKrd85p49e5g5cyY//vgjqamp1KpViw4dOjBmzBhat25daNcmIlJS2Z51bTQYKiKlidVqJTw8nN69ezN16lR69uxJdnY2fn5+9O7dm/j4eNLS0oz3N2vWjKZNmxqxcty4cZjNZsLCwsjKyqJJkyYAhIWFMX36dMLDw6lcuXKRXJv870xW2xBNETl48CBNmzblu+++u23H6YEHHqBv375MmjQJwJheZ7Va/+uOU7du3WjdunWB0mkRkftVly5dSEtL49lnnyUiIoIDBw6watUqIiMjiYmJITU11RgUuJ0bpz5bLBbmzJnDwYMH6datm276IiIiIvJfsfXlt27dSqdOnahTpw6HDh3C3t6eK1eukJ2djbe3t/G+3Nxcxo8fT2ZmJt999x15eXlMnz6dfv36AXDmzBkqV66sjR9KuCJPSoE6TiIi90JSUhKzZs1i27ZtpKSkkJeXx2OPPcZLL72Er6+vUVkqIiIiIlJYcnJy2LZtG4GBgdSuXbvAazcWoZhMJt566y0SExMJDw8nIiICLy+vomm03DPFIimljpOIyL2zb98+nJ2dqV+/foHj2thBRERERESKUrFIStmo4yQi8te6OX7euJupiIiIiEhRUT9foBglpdRxEhG5d3TTFxERERGR4qbYJKVs1HESEREREREREbn/Fbtl6pWQEhERERERERG5/xW7pJSIiIiIiIiIiNz/lJQSEREREREREZFCp6SUiIiIiIiIiIgUOiWlRERERERERESk0CkpJSIiIiIiIiIihU5JKRERERERERERKXRKSomIiIiIiIiISKFTUkpERERERERERAqdklIiIiIiIiIiIlLolJQSEREREREREZFC9/8Bmfj7vlbpKVMAAAAASUVORK5CYII=", "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 }