{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "36722953",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import statsmodels.formula.api as smf\n",
    "from toolz import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f4484b",
   "metadata": {},
   "source": [
    "# Why Prediction Metrics are Dangerous For Causal Models\n",
    " \n",
    "A common misconception I often hear is that, if we have random data, to evaluate a causal model we could just evaluate the predictive performance of such model on the random dataset, using a metric like $R^2$. Unfortunately things are not that simple and I'll try to explain why. \n",
    " \n",
    "Generally speaking, we can say any outcome is a function of the treatment and covariates\n",
    " \n",
    "$$\n",
    "Y = F(x, t)\n",
    "$$\n",
    " \n",
    "Let's say we can decompose this function into two additive pieces. One piece that doesn't depend on the treatment and another that depends only on the treatment and possible interactions.\n",
    " \n",
    "$$\n",
    "Y = g(x) + f(t,x)\n",
    "$$\n",
    " \n",
    "This additive structure places some restriction on the functional form but not much, so we can argue it is a pretty general way of describing a Data Generating Process (DGP).\n",
    " \n",
    "The point is that if the treatment effect is weaker than the covariates effect, then, **even if we have random data**, we can have a model which has higher predictive power but is bad for causal inference. All that model has to do is approximate $g(x)$ while disregarding $f(t,x)$.\n",
    " \n",
    "In other words, **predictive performance on a random dataset does not translate our preference for how good a model is for causal inference**.\n",
    " \n",
    "To show that, let's use some simulated data.\n",
    " \n",
    "## Simulating Data\n",
    " \n",
    "In the following DGP, we have covariates $X$ that have a high predictive power but don't interact with the treatment. In other words, they don't dictate the treatment effect heterogeneity. We also have features $W$ that only impact the outcome through the treatment (are not confounders). Since the treatment has low predictive power, $W$ also doesn't have much predictive power.\n",
    " \n",
    "$$\n",
    "Y_i = g(X_i) + f(T_i,W_i) + e_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe6e5192",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100000\n",
    "n_features = 20\n",
    "n_heter = 10\n",
    "\n",
    "np.random.seed(12321)\n",
    "\n",
    "X = np.random.normal(1, 10, (n, n_features))\n",
    "nuisance = np.random.uniform(-1,1, (n_features, 1))\n",
    "\n",
    "W = np.random.normal(1, 10, (n, n_heter))\n",
    "heter_y = np.random.uniform(-1,1, (n_heter, 1))\n",
    "\n",
    "T = np.random.normal(10, 2, (n, 1)) # T is random!\n",
    "Y = np.random.normal(T + T*W.dot(heter_y) + 20*X.dot(nuisance), 0.1)\n",
    "\n",
    "df = pd.concat([\n",
    "    pd.DataFrame(X, columns=[f\"f{f}\" for f in range(n_features)]),\n",
    "    pd.DataFrame(W, columns=[f\"w{f}\" for f in range(n_heter)])\n",
    "], axis=1).assign(T=T, Y=Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e943c9d8",
   "metadata": {},
   "source": [
    "Now, let's break that dataset into a training and a test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "48674d02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((50000, 32), (50000, 32))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "train, test = train_test_split(df, test_size=0.5)\n",
    "train.shape, test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85db5a8",
   "metadata": {},
   "source": [
    "and train two models for treatment effect heterogeneity, $M1$ and $M2$. $M1$ will include the highly predictive features that don't affect the treatment heterogeneity and $M2$ will include the low predictive features that do affect treatment heterogeneity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "12490bcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "m1 = smf.ols(\"Y~T*(\" + \"+\".join([f\"f{f}\" for f in range(n_features)])+\")\", data=df).fit()\n",
    "m2 = smf.ols(\"Y~T*(\" + \"+\".join([f\"w{f}\" for f in range(n_heter)])+\")\", data=df).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "642418c0",
   "metadata": {},
   "source": [
    "If we look at the predictive power of both models using the $R^2$, indeed, $M1$ is much better than $M2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5b4facc6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M1: 0.9160516511287358\n",
      "M2: 0.08378351037639298\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import r2_score\n",
    "\n",
    "print(\"M1:\", r2_score(test[\"Y\"], m1.predict(test)))\n",
    "print(\"M2:\", r2_score(test[\"Y\"], m2.predict(test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e2c5510",
   "metadata": {},
   "source": [
    "Now, let's calculate the cumulative elasticity curve for both models. For that, we will need Conditional Average Treatment Effect predictions. Since all the models here are linear, we can use the following formula to get CATE predictions\n",
    " \n",
    "$$\n",
    "\\hat{CATE_i} = M(X, W, t) - M(X, W, t-1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1afeeee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "@curry\n",
    "def elast(data, y, t):\n",
    "        return (np.sum((data[t] - data[t].mean())*(data[y] - data[y].mean())) /\n",
    "                np.sum((data[t] - data[t].mean())**2))\n",
    "    \n",
    "\n",
    "def cumulative_gain(dataset, prediction, y, t, min_periods=30, steps=100):\n",
    "    size = dataset.shape[0]\n",
    "    ordered_df = dataset.sort_values(prediction, ascending=False).reset_index(drop=True)\n",
    "    n_rows = list(range(min_periods, size, size // steps)) + [size]\n",
    "    return np.array([elast(ordered_df.head(rows), y, t) * (rows/size) for rows in n_rows])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4261421b",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_pred = test.assign(\n",
    "    cate_1 = m1.predict(test) - m1.predict(test.assign(T=test[\"T\"]-1)),\n",
    "    cate_2 = m2.predict(test) - m2.predict(test.assign(T=test[\"T\"]-1))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "945b269f",
   "metadata": {},
   "source": [
    "Once we have those CATE predictions, we can evaluate them using the cumulative elasticity curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f2605203",
   "metadata": {},
   "outputs": [],
   "source": [
    "cumelast_1 = cumulative_gain(test_pred, \"cate_1\", \"Y\", \"T\", steps=100)\n",
    "cumelast_2 = cumulative_gain(test_pred, \"cate_2\", \"Y\", \"T\", steps=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2d2407dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABQc0lEQVR4nO3dd3hURdvA4d+kN5JAeiX0EiCU0IuA0gRREOlIEfC1gKiIfnZf9bX3giICoghIURBBpfcaCL0lkJAAqUAgCSmbzPfHSUICSVhCtmQz93XtlezZ3XOew4ZnZ+fMPCOklCiKoiiWx8rUASiKoiiGoRK8oiiKhVIJXlEUxUKpBK8oimKhVIJXFEWxUDamDqA4T09PGRISYuowFEVRqoyIiIgUKaVXaY+ZVYIPCQlh3759pg5DURSlyhBCxJb1mOqiURRFsVAqwSuKolgoleAVRVEslErwiqIoFkoleEVRFAulEryiKIqFUgleURTFQqkEr1QvUsLBRRCv5lsols9gCV4I0UgIEVnsdlUIMc1Qx1OU28rLhRVPw++Pw+x74Y+nID3Z1FEpisEYbCarlPIk0BJACGENnAd+N9TxFAWA7HSwcwYhSm7PugpLxkL0Bug6HfJzYee3cPxPeOBzaDa47H1Kqd2s1BdepWoxVqmCe4FoKWWZU2oV5a6lRsN3XcDFB1qOhOZD4FoCRK2Ho8vhciwM/ApaP6o9v+VoWPk0LHsMrKyh6YO37jP3Osx/CFy8YNgvRj0dRblbxkrww4GFpT0ghJgMTAYIDg42UjiKRVr7OiDALRA2vqvdAIQ1BLaF/p9AvZ43nu/VEEYvh18Gw9LHYLgjNOx943EpYdWzELdLux+1Hurfa7TTUZS7JQy9JqsQwg64AIRKKRPLe254eLhUxcaUConZDvPuh56vQrcX4Mo5OPEXuAZAnW7g6F72a7PS4KeBkHQcHvwaQgeDtQ3s+QFWT4cuz8GRZWDvCo9v1lr7imImhBARUsrw0h4zRqdiP2D/7ZK7olRYfj7887KWzDs8pW1zD4YOT0DTgeUndwAHNxjzu9aiXz4JvgiDf1+Fv1+Chn2h52tw7+uQeBgO/XbjddevQEbq3cefeQkOLdG6gxSlEhkjwY+gjO4ZRakUh5fAxUi49w2wc6rYPpxqweTNMGwB1KoDO77SPiQGfa9dXA0dDP6tYMM7Wot/+xfweXP4rvPdJfn8fFgyDpZPhC9bw7652mgfRakEBu2iEUI4A+eAulLKtNs9X3XRKHrLvAQX9sOFA1pXiqs/TNxQeSNdUqK0lr1LsXUUzm6Bnx4AuxqQcw3qdofYHVC3B4xcfOvIHdC6dk7+Da3HQEjXW5+z/UtY+xp0mgLndkP8HvBqDI/9qx1fUW6jvC4ag15klVJmAB6GPIZSzSQcge2fw5HlIPO0bZ4NYcBnlTuM0bP+rdvqdIPmj0DyCbjvLe2C6+7vYc0M2DUTOj5Z8vnHVsKyiYCAw7+Bd1NoNwmaDwV7F7gQCev/C40HQK+3tdccWaaN6on8VetiUpS7YPCLrHdCteCVMl1LhJVT4PQ/YOcCbcZBwz7g2+L2feyGJCUsGgmn18LEtVo3Dmgjbn4dBgGtYcQiOLla+zBIOKRdrA0bDtEbIScDntiudREVmt0LMlPh6X1q7L1yW+W14FWCV8zflTiY/yBcu6iNaGk3ERxrmjqqGzIvwczOkJ6ofZvwaQonVmvfAsauuvEBJCXE74W9s+Ho71pf+6MroO49Jfd3aInWJz96GdS/z+ino1QtKsErVVdqtJbcs67CqCUQ3N7UEZUuJQoO/qp1ISUegRq+MGJxyT784jJS4OoF8Gtx62O6HPgsVGv9j1xs2LiVKs9kffCKUiZdttaX7RdW9nPS4mFuP8jXwbg/y3+uqXnW14ZS6svZU7uVxsZO64La8hFcOquN6lGUClAdfIrx5efD0gnwfTc49W/Zz9v0Ply/DOP+Mu/kbgjh47UJVXtna99ejiyDdW9B4jFTR6ZUIaoFr1S+yzHa0MCA1uDXUmuRFrfpf3BilXaxcc0MbXSKrUPJ56RGayNJ2k0C7ybGitx8uPpDkwe0IaC7v9eKowFs+1SrmXPPi+ATatoYFbOnWvBKxV06o40Nv9nfL8PfL8KPveD9IG3s+P752gShw0u1rodWY2DofLh8FnZ8ees+Nr0P1nbaRdXqqstz4Nsc2j8O4/+GF6K1MgxRG7SLuv+8ArlZ2nPz87XW/qdNteGZFXHxkPZtQbEY6iKrUjFSwg89IPEoPHNQa3GCNuLlixYQ/pjWMj+3SxvamBoFNg4g87XCX2P+0Fr2v42FU3/DU7uhZoi2j6QT8G0H6DwVev3XVGdovjIvwYa3Yd8c8A6F+97UPiRjtmqTsJAwaQN4NdJ/nzHbtA/i0EEwZI6hIlcMwNS1aBRLFLVOm0Wal6PNxiwUMVf72XmqVgem7/+08dwT12ut9pAuMPTnG902ff6nVXtcXVAgLD9P68Kxc4HO04x+WlWCUy1tYtfIJZCRDL8+ok2aeuBLeGqX9kG6eIxWG/9mGSlw4BetC6xQerJWTVPma/XxK6O+jmIWVAteuXNSwuz7ID0JgjtoSWHaIW1q/WehWgt9xB2UH9r+RUGpX8DaHvKytT7mHi8bJn5LkpECB36GZkPAPUjbdmYz/PyQ1lff623IvqrNITi4CI6t0D6UbZ2g73vah+4vD8O5nVqt/OWTtA/djk+Z9LQU/alx8Erlilqv1VAf8LlWX+WbttBpqtZfvOyxO5+gI6XWlZNyUuvXz0rTEpODq8FOweJt+wzWvVlyW+EM2qYPwuYP4exmbWJWyil44AttaOYP90JOOjy5q/TaOorZUePglcojJWz+AFwDoeUorasldLB2ga9WXe1Wt+ft91OcEFC7o3ZTKkfnaVCzjtZ6t3fVPiwD22k1cACCO8Gub7RaOC2GQ+ux2vY2Y7WSEHF7zHdSmaI3leAV/Ump1VSJ262tjlTYj971eTiyVKuz0vtdVT/FHAgBoQ+V/biVlVbBsuUocHC/0VoPHQx//5826kkl+CpP/U9Ubi/ltHYR9MuWWmEtt2Ct77aQT1NtzLatk7YWqlJ1ONUq+YFs76ItQH50uRoyaQFUglfKl58HC0fA/p/Bqwn0/xQmrQcb+5LPG/g1TN5UsiqiUjW1Hge5mbDza218vVJlqS4apXzH/4TU0/DIPG2MdFkc3U1btlepPAGtod692rWWqHXaqJrgDqaOSqkA1YJXyialNjXeoz40GWjqaBRjEQJGLYWHZmoVL+f0gWWTtPHyhXIyYOsnBWUU8kwXq1Iu1YJXyha9Hi4e1LpfrKxNHY1iTFZW2vWUpg9qQy63fQ5Ra7WL6DJfW5s2PUF77tE/YNDMGzORFbOhWvBK2bZ+Bq4B0GKYqSNRTMXOGXq+Cv/Zpo2ZX/EkrHxam1Q14V946Dut/v3MztoyiopZMWgLXgjhDswGmgESmCCl3GnIYyqV5NwuiN0Gfd+/tRqkUv14N9YKnh1drpVCaNxf68oJbg8hnWHJeFjxtFaKwsX7xus2fwS7vi0Yhlkw32HYLyY7jerG0C34L4C/pZSNgTDguIGPp1QGXbZWOsCxFrR+1NTRKObCygqaD4EmA0rOcnUPhsGztBITm96/sf38ftj4rlbWOHSwdi3n+J9q+KURGSzBCyHcgG7AjwBSyhwp5RVDHU+pJFIWzGTcDfd/pH1FV5Tb8agHbcZDxDxt3kR+Hqx6VmvND18A/T+GrgWlnxOPmDTU6sSQLfg6QDIwVwhxQAgxWwhxS7YQQkwWQuwTQuxLTk6+dS+KcW16Dw4t1vpdmw8xdTRKVXLPi2DrqNXA2fsjXIzUCpo5uGmPF67KdfGQqSKsdgyZ4G2A1sBMKWUrIAN46eYnSSlnSSnDpZThXl5lLFCsGMe+udrY51ajoet0U0ejVDUuXtD5GW21rrWvQ90eWtdMoRq+4OytjcxSjMKQCT4eiJdS7i64vxQt4SvmJjsdVjwFq6ZpE1wGfK4qCSoV0/EpcPHRhlL2/+TWvyO/MK1mkWIUBhtFI6VMEELECSEaSSlPAvcCasVgc3PxECwdry0A0XU6dH8JrG1NHZVSVdk5w6glcP2K1i9/M78WEL1BW2rw5nV4lUpn6IlOU4AFQgg74Aww3sDHU+5EehLMf1Ab9jZ2pbbEnqLcrcK+9rIek3mQdEwriaAYlEETvJQyEii1EL1iYlJqoxxyMmDC33e2fqeiVJRvC+1nwiGV4I1AzWStro4s0y6G9XhZJXfFeGqGgL2butBqJCrBV0fXEmH1dAgI1xZ9UBRjEULrh1dDJY1CJfjq6O8XISdTqxaoiogpxubbQpvslKczdSQWTyX46ibrqjZdvN0k8Gpo6miU6sgvDHRZ2joDikGpBF/dxGyDfB007GvqSJTqyq/gQqvqhzc4leCrm+j1YOsMQWpBZcVEPBqAjaPqhzcCleCrm+gNUKerKgGsmI61jVZhUrXgDU4l+Ork0lm4dAbq9TR1JEp1598KzkdAoprcbkgqwVcnZzZqP1WCV0yty7Nalclfh2rDdhWDUAm+OoneAG7B2sILimJKbgEwchFkpsLCYdqwXaXSqQRfXeTp4MwWqNdDVYpUzIN/K3j4R7gQCYtGwKl/VaKvZLdN8EIID2MEohjY+QjITlPdM4p5aXy/VlY4bg/8+gh8WAeWTVSToCqJPi34XUKIJUKI+4VQTb8qK3oDCCuoe4+pI1GUkto+BjPOwujl0GIoHF4CBxeaOiqLoE+CbwjMAsYAp4UQ/xNCqCmQVUl6MhxfCQFtwLGmqaNRlFvZOkD9e+GBL8G/tbaymC7b1FFVebdN8FKzVko5ApgEjAX2CCE2CyE6GjxCpeJyMmDzR/BlS0g+CW0nmjoiRSmfENp6wGlxsH++qaOp8m5bD76gD340Wgs+EW0Rj5VAS2AJ2uLairk5HwGLH4Wr8dB4ANz3Jng2MHVUinJ79XpC7c6w5SNoOQrsnEwdUZWlTxfNTsAVeEhK2V9KuVxKqZNS7gO+M2x4SoUc+AXm9NP63Mf/DcMXqOSuVB2Frfj0RNg729TRVGn6JPhXpZRvSynjCzcIIR4BkFJ+YLDIlDsnJayeoS2gHdwBJm+C2qoXTamCanfSFoDf9hlkpZk6mipLnwT/Uinb/q+yA1Eqwcb/wZ7vocOT2ogEZzXCVanC7n0drl+CbZ+bOpIqq8w+eCFEP+B+IEAI8WWxh1wBvQapCiFigGtAHqCTUqr1WQ0lciFs+RBajYY+/1OTmZSqz78ltBgGu76F8AngHqRtP7Ea1v8Xxq0CZ0+ThmjuymvBXwD2AVlARLHbSqDPHRyjh5SypUruBhSzHVZOgZCu0P8zldwVy9HzNa3rccM72v3z+2HpBEg+Dmc2mTS0qqDMFryU8iBwUAixQEqpppWZq6sXYfFobTHjYT+rMsCKZXEPgo5Pan3xTR6Av54DZy+4fhlit0PzIaaO0KyV2YIXQvxW8OsBIcShm2967l8C/wohIoQQk8s4zmQhxD4hxL7k5OQ7DL+ak1K7oJp7HUYsVJOYFMvU5Vlw8oDFo7S/9VFLIKgdxO40dWRmr7xx8M8U/BxwF/vvIqU8L4TwBtYKIU5IKbcUf4KUchbaTFnCw8PlXRyr+tk7W1uh6f6P1TBIxXI5uMG9b8CaF2HoT+DdWBsdtuEdyLwETrVMHaHZKrMFL6W8WOw5iVLKWCllLJAE6NXJK6U8X/AzCfgdaHd34SpFUqLg39e0oWRqhqpi6dqMhRdjbhTLq91Z+3lul8lCqgr0GSa5BMgvdj+vYFu5hBDOQogahb8DvYEjFQlSuYkuG36fDDb28OA36qKqUj3YOtz43b81WNtp/fBKmW5bqgCwkVLmFN6RUuYIIfS5kucD/F5QgNIG+FVK+XfFwlSK5Om0UQTnI+CRn8DVz9QRKYrx2TpoxfPOqX748uiT4JOFEAOllCsBhBAPAim3e5GU8gwQdpfxKcXl58OKJ+HEKuj3IYQ+ZOqIFMV0gjvCji+1onp2zqaOxizp00XzH+BlIcQ5IUQc8CLwuGHDUkrIz9eqQf45FQ4t1up0tFdvgVLN1e4M+TqI32vqSMzWbVvwUspooIMQwqXgfrrBo1I0GanaMMjY7ZB9VdvWeRp0nW7SsBTFLAS10wrqxe6Aut1NHY1ZKq9UwWgp5S9CiOdu2g6AlPJTA8dWvUkJK5/WhkG2HAWB4RDYFrwamToyRTEPDq7g00xL8EqpymvBF3Zq1SjlMTVe3dD2zYGTq7W6Mh2fMnU0imKeaneCiJ9Al6NmcZeivFIF3xf8uk5KWWIskhCis0Gjqu6ST8I/r2hjfts/YepoFMV81e4Mu7+Dc6qbpjT6XGT9Ss9tSmXITodlj2mr2Dw0E6z0eYsUpZpq0Eub6Rrxk6kjMUvl9cF3BDoBXjf1w7sC1oYOrFpKOgG/PQqpp2H4Qqjha+qIFMW82TpC2AjY+6O2uLyLl6kjMivlNQ/tABe0D4EaxW5XAVXCrbIdWgI/9NAWOBjzBzTqa+qIFKVqaDMe8nMhcoGpIzE75fXBbwY2CyHmFdSgQQhhBbhIKa8aK8Bq4dxuWD5Rm7gxZK6anaood8K7MQR3goh50Gmq6tYsRp9/ifeEEK4F9WSOAMeEEC8YOK7qJWqdNp535GKV3BWlIsInwOWzcHaTqSMxK/ok+KYFLfaHgDVAHWCMIYOqdmK3g28L7WKRoih3rulAcKwF++aaOhKzok+CtxVC2KIl+JVSylzUOPjKk5sF8fsgpIupI1GUqsvGHlqN0uaOXL1g6mjMhj4J/nsgBm3i0xYhRG20C61KZTi/D/Kyb9S3VhSlYsIf07o6//4/U0diNm6b4KWUX0opA6SU90tNLNDDCLFVDzHbAaGtUKMoSsXVqgP3vAjH/oDjf5o6GrOgT7lghBD9gVCgWMV9/muQiKqb2G3g20ytp6oolaHzM1qC/+t5rduzmv+/um0LXgjxHTAMmIK2VN8jQG0Dx1U96HIgbi/UVv3vilIprG1h4NeQkQL/vmrqaExOnz74TlLKR4HLUsq3gI5AQ8OGVU1c2A+66xCi+t8VpdL4t4TOU+HALxC3x9TRmJQ+Cf56wc9MIYQ/kAuowdqVIWab9jO4k2njUBRL0+0FsHXWknw1pk+CXyWEcAc+AvajjahZaMCYqo/Y7eDdFJw9TB2JolgWO2do3F/rj9dlmzoak9FnFM3bUsorUsplaH3vjaWUr+l7ACGEtRDigBBi1d0EanHycrUSBWp4pKIYRouhkJUGp9eaOhKTKa+a5OByHkNKuVzPYzwDHEerQqkUit0BuRmq/11RDKVuD3DyhMO/QZMBpo7GJMobJvlAOY9J4LYJXggRCPQH3gWeu83Tq4+MVFjxNLgGaIt6KIpS+axtoNlgrVZ8Vlq1LAVSXjXJ8ZWw/8+BGZS+7B8AQojJwGSA4ODgSjikmcvTwbIJkJ4IE9ZUyz86RTGa5kNhzyxt4lOr0aaOxujK7IMXQnxe7Pdnbnps3u12LIQYACRJKSPKe56UcpaUMlxKGe7lVQ2K9W/4L5zZBP0/gYA2po5GUSxbYDjUrAOHl5g6EpMo7yJrt2K/j73psRZ67LszMFAIEQMsAnoKIar3mKWYbbD9C620aWtVkFNRDE4IaP4InN1SLYuQlZfgRRm/60VK+X9SykApZQgwHNggpax+35GKO7YSbByhz/9MHYmiVB9hw7UiZL8OhWsJJR9LPqU1vKLWwal/ISfTNDEaSHkXWa2EEDXRPgQKfy9M9GpN1oqI3qCNmrF1NHUkilJ9eNSDEYu19Y5n94LRSyHzEmz+AM5sLPnckK4wejnY2Jkm1kpWXoJ3AyK4kdT3F3vsjurBSyk3AZvu5DUW50qctph2eGVcu1YU5Y40uA/G/wULhsJ3XSAvB5y94L63IKA12DjAxYOwejqsmgYPfqN171Rx5Y2iCTFiHJavsKVQV1VaVhST8G8FE9dq9eJDumiLdds53Xg8qB1kpsKm97TSw92q/sqkepULVipB9EZw8QXvJqaORFGqr5ohMKKcSiv3vAiXzsCGd8CnGTTqZ7TQDEEtP24M+XlaC75eT4v42qcoFksIGPgV1KoH2780dTR3TSV4Y7h4EK5fhnqqe0ZRzJ6NvTYp6twOSI02dTR3Ra8EL4ToIoQYX/C7lxCijmHDsjBF/e/dTRqGoih6KhxaebBqF87VZ0WnN4AXgcKVbG2B6j1h6U5FbwSf5uDibepIFEXRh6u/1qUauRDy800dTYXp04IfBAwEMgCklBcop7aMcpOcDDi3S3XPKEpV03IkXI2HmC2mjqTC9EnwOVJKScHYdyGEs2FDsjDHVkB+rkrwilLVNOoP9m4Q+aupI6kwfRL8b0KI7wF3IcQkYB3wg2HDshD752tlgf3C1MIeilLV2DpA84e1EiNZaaaOpkL0WdHpY2ApsAxoBLwupfzK0IFVaVLCxv/ByinahdVxf2lX5hVFqVpajgbddTii7/pG5uW2E52EEM8Bi6WU1Xfdqzu1+QPt1mo0DPgcrG1NHZGiKBUR0Fr7Br71E21kTRWrI6VPF00N4F8hxFYhxNNCCB9DB1WlndkMm96HsBEw8GuV3BWlKhMCer8DaXGw8xtTR3PH9OmieUtKGQo8BfgBm4UQ6wweWVWUngTLJ4FnQ21BDzVrVVGqvjrdoPEA2PYZXEs0dTR35E5msiYBCUAqoAZ03yw/T0vuWVfhkXlgpwYbKYrF6PVf0GXDhrdNHckd0Wei05NCiE3AesADmCSl1GdFp+pDSvj3VW0pvvs/BJ+mpo5IUZTK5FEP2j8OB36Bi4dMHY3e9GnBBwHTpJShUso3pZTHDB1UlSKlVn5017fQ/glopZbiUxSL1O0FcKwJ/76i/b+vAspbdNu14NePgHNCiFrFb8YJz8xJCatfgN0zocNT0Pc91e+uKJbK0R26/5+2vuvpf00djV7Ka8EXTt+KAPYV/Iwodl/Z+TXs/QE6TYE+76rkrgCwfH880xYdIFuXZ+pQlMoWPl4rJfzva5CnM3U0t1Xeik4DCn5WqHKkEMIB2ALYFxxnqZTyjYrsyyzl5cLOb6HOPdDrbZXcFQB0efl88PcJEq9mY2NtxUdDWiDU34blsLbVLrguHgX7f4K2j5k6onLpc5F1vT7bSpEN9JRShgEtgb5CiA53HKG5OrEKrl2ADk+o5K4U2XgymcSr2XSs68HSiHjmbI8xdUhKZWvcXys9svF/2qg5M1ZeH7xDQV+7pxCiZrH+9xAg4HY7lpr0gru2BbeqcWVCH3t+APfa0KC3qSNRzMivu2PxrmHPTxPa0SfUh3f/OsaWU8mmDkupTIWTnzJTYNunpo6mXOW14B9H629vTMn+9xXA1/rsXAhhLYSIRBtDv1ZKufuuojUXCUcgdju0nQhW1qaOxuD2n7vM9qgUoxxry6lkHvpmO0cvVL3iTuevXGfTqWSGtQ3CzsaKT4e2pKFPDaYtjiQrV/XHW5SA1tps9R1fQ8ppU0dTpvL64L8AvhBCTKlocTEpZR7QUgjhDvwuhGgmpTxS/DlCiMnAZIDg4OCKHMb49swCG0et1kwVkJ8vGTNnN1FJ6dRwsMXVwQY/d0caeLtQ39uFy5m57Iu5RETsZdrUrsnHj4Rha6199kclXWP07N1YC8GeV+7D0c5wH2g7o1OZNH8f2bp8Jszby+9Pdsbf/Ubtj9y8/KK4zNHiPecAGNY2CABnexteH9CUkbN3s/rwRQa3DjRleEpl6/VfOLEa/noeHl1hll21+pQq+EoI0UwIMVQI8Wjh7U4OIqW8AmwE+pby2CwpZbiUMtzLy+tOdmsamZfg0G/QYig4VY3RoquPXGR7VCrNA9xp6OOCo501h+PT+GL9aZ7+9QCv/XGEHdGp1PF0ZkXkBV5efhgpJRnZOv7zy36khGvZOtYcuVipcaVn68jP13rt9sVc4rGf9hJcy4lfJ7UnMzuPcXP3kHY9lwtXrvPcb5E0ee1vNp5MqtQYKosuL5/F++K4p6EXgTWdirZ3qOtBiIcTCwuSv2JBXLzh3tfg7GY4sszU0ZRKn2qSbwDdgabAaqAfsA2Yf5vXeQG5UsorQghHoBfwwd0GbHIHF2rlQ9tNNnUkesnPl3y5/jT1vV34fkwbrK1utDKu5+QRnZxODQcbgms5IYTgs7Wn+GL9aXzdHIi7lEl0cjo/T2jPq38cZtHeuEprhc7aEs17a05ga2WFn7sDydey8XV1YMGk9njXcOD7MW0YO3cPD8/cQdylTCTg6mjL+6tP0K2BV4nzMAcbTiSReDWbtx8s+S3Uykowol0w7605wenEazTwUYuhWZTwCdrs1n9e0a7HObje/jVGpM/33SHAvUCClHI8EAa46fE6P2CjEOIQsBetD35VhSM1FydWg29z8G1m6kj08vfRBE4lpjOlZ/1bkqKjnTXNAtyo7eFcNJRv2n0NGN42iK82RPFH5AWeu68hXRp4MrRtEHvOXuJMcnrR6/PzJYlXs9Dl3dmalbvPpPL+mhN0beDF+C4hNA9wo0cj76LkDtCpvicfPNyCmJQM7m/ux4bn7+G/D4ZyMvEaKyLP3+W/SuU6m5LBW38ew9fVgZ6Nby3TNKRNILbWgl9VK97yWFnDgE8hPRG2fGTqaG5x2xY8cF1KmS+E0BXMbk1CK19QLinlIaDV3QZoVrKvQdxu6PiUqSPRS2Hrva6XMwNa+Ov1GiEE7zzUjGxdPvlS8lSP+gAMaR3IJ/+e4rd98bzUrzHZujzGz93LjuhUrAR41bCnbUgtPn4kDAfbG/30P++MYcHuczzbqyG9m/qQmpHDlIUHqO3hzDcjW1HDoexyyoNbBzKghT92Nlo7xN/NkWYB0Xy69lSJ7cYUm5pBeraOJr6uWFkJjpxPY+ycPUhg3vi22JRyjcDDxZ4+ob4s33+eF/s2LvHvo1iAgDYQOkhbwa3nq2a1uI8+CX5fwUXSH9BG0aQDOw0ZlNmK2V6wvmpPU0eil3+PJXAi4RqfDQu7oy4NG2srPhvWssQ274LW6dKIeJ7t1YDnfzvIjuhUnuxeDxsrQdzl6/wReZ4cXT4zR2tdQX8fSeD1lUdxsrXm8Z8j6NnYm+s5eaRdz2Xe+HblJvdCxZO4lZXghT6NGTtnDwv3nGNspxC9z+luSCnZcjqFudvPsumkNuTRw9mOjvU82HQyGVcHG36e2J56Xi5l7mNk+2BWHbpY6sXWv49c5P+WH6ZDXQ/ua+JDz8be1HS2M+g5KZWs1Wg4uhxOroHQh0wdTZHbJngp5ZMFv34nhPgbcC1onVc/ZzZqo2eCzX++Vn6+5Iv1UdTxdOYBPVvvtzMsPIi1xxIZ9cNu9sVe5qV+jfnPPfWKHg8LdOPNP4/x2oojPNImkGmLDxAW6M7Pj7Vj8d44Plt7ioycPN4f3Jym/hXrq+zWwJMOdWvx1YbTDGkTiLO9Pm2UitsZncrbq45x7OJVPF3smXZfAwJrOrHtdDLbolII8XTih0fD8XMrf6WfjnU9qOPpzK+7z92S4OcWTIbaf+4ya44k4GBrxbzx7ehQ18NQp6VUtrrdoYa/do2uKiR4IUTr8h6TUu43TEhmLHoDhHQ2q69gZVm2P57jF6/yxfCWpXYbVET3Rl5417BnX+xlxnaszePd6pZ4fFznOiRey2bmpmiWRcTj7WrP7LHh1HCwZWLXugxo4c/h82nc16TiywkIobXiH565gyX74hjXuUKVNG7rbEoG7685zj9HEwlwd+TjR8J4IMwPexute2VIm0CklHqXIRBCMLJdMO+uPs7RC2mE+muXsc5fuc7us5d4rldDpvSsz+HzaTz320Em/bSPxY93rPAHoWJkVtYQNgy2f6kt/ONiHktmlNf8+aScxyRQNfopKktaPKScgjbjTB3JbWVk6/jon5OEBblXWusdtK6bVwc05eiFNGb0aVxqcpvRpxGX0nNYfyKRuePa4ely48PQ180BXzeHu46jTe2ahAW6sXBPHGM7hVRarZfMHB1rDiewJCKOXWcu4WhrzfO9GjKpW91S+83v9LhDw4P4bN0p5m6P4eNHwgD444B2wXhQqwCEELQIdGf+hHY8PHMHY+fuYfkTnQiq5VTebhVzETZSW/Xp0G/Q6WlTRwOUP9GphzEDMXvRG7Wfdc3vn+VaVi7OdjZYFfSzf785mqRr2cwc3aZoW2UZGObPwLCyPzSEEHwwpIXBJyWNbB/Mi8sOs//cZdrUvvv5CJk5Oh74ahvRyRnU9nDi+V4NGdY2CG/Xu/9AKuTmZMuQNoEs2hPHi30b4+lix+8HztM2pGaJJO7v7sj8Ce0Y8t1Oxvy4m18ndSgx4as8qenZeLiY/zdMi+TVEALCIXKBNhDDDCY+6VNs7NHSbsYIzqxEbwAXX/BuYupISjideI12765n4Dfb2H0mlfNXrvP9ljM8EOZPm9o1TRaXoWecDmjhj4u9Db/ujquU/X38zymikzOYOao1m6Z3Z8q9DSo1uRca37kOOXn5/LIrlqMXrhKVlM6gVrfOLWjgU4M549qSmp7D4G93cCKh/KJWWbl5vL7iCG3eWVf0rUAxgZYjIOkYXDxo6kgA/cbBty126wq8CQw0YEzmJz9fW46vXg+z+FQulKPLZ9riSBxsrbiUnsOwWbsY/O12AF7s28jE0RmWs70ND7b0Z9WhC6Rl5t7VviJiLzF3x1nGdKhNv+Z+Bi3vW8fTmXsbe7NgdyyL9p7DztqK/s39Sn1um9o1+e0/HZFIHpm5kx3RpdcDikq6xkPfbGf+zlhcHWz4bnM0soqsOGRxmj0M1nba5CczoM8ominF7xcMmVxkqIDMUsJBuH7J7IZHfrn+NEcvXOX7MW24p6EXs7eeYeamaJ7qUb/EdHlLNaJdMAt2n+P3A/EVvtialZvHC0sP4e/myIv9GldyhKWb0KUOo2bv5pdd5+gb6oubU9nDRZv4ubL8yc6Mm7OH0bN3U8/LhVB/V0I8nYm7dJ2TiVc5mXANVwdb5o5rS9K1LF5cdpid0al0qu9plPNRinGsqZUxiZgLrR8FP9MuX12RMWYZgGGGLpirqILy93W7mzSM4iJiL/PtpigeaRNIn1BfAJ7u2YD/3FPP7KbxG0qzALcKXWyNu5TJ8YtXuZqlY8upZM4kZzB/QjtcDDzkslCneh408qnBycRrDGp928rbBLg7svQ/nZi74yyH49PYdeYSf0RewKuGPY19azC+cx0e61IHH1cHsnLz+ODvk8zZflYleFPp9Tac+hf+eBImb9QWCTERfWrR/MmNOu5WaDVpfjNkUGYlT6fNUAtsZzZDn3adSeWFpQfxd3fk9QealnissoZEVhUj2gXz0vLDfLH+NCPbBxeVOihL/OVMen+2hevFyveO6xRCt4bGK3QnhOD53g35cdtZejTS72/KzcmWafc1LLqflZtX6sgeB1trRrcP5quNUZxNyaCOp/MdxZZ2PZctp5LZfCoZKwH1CyqOdqrnqWbg6supFjzwOSwaCVs/he4vmiwUcbu+OiHEPcXu6oBYKWW8IYIJDw+X+/aZ2XKvh36D5ZNg+K/aSi4mFBF7mU/XnmR7VCreNeyZObqNSS+kmoPMHB2T5u9je1Qq1laCHo28ePuhZmVOPHpyQQQbTiQxb3w7/N0ccXW0wd3JsmaNJl3LovP7GxjZLpi3HtSvZlJmjo7nFh9k7fFE8vIl7k622FgJUtJzALi3sTc/jmtryLAtz7KJcPR3mLxJq19lIEKICClleGmP6dMHv7lgJ66FzxdC1JJSXqrUKM1Rfr72CezVBBr2M1kYUkq+3RTNR/+cxNPFjlf7N2F0h9qqRQU42dmwYGIHopLSWRoRz9ztZ/nk31NF48yL2xGVwurDCTzXq6FFzxL1ruHAA2H+LImI57nejXBzLL+LIFuXx+M/R7A9KoWJXevSJ9SHlkE1sbYSXMnM4YetZ/hmYzQRsZcqZUhqtdHvQ21wxpoXYfxqk4SgzzDJyUKIBOAQsA+tHo2ZNbMN5NQaSD4OXZ8DK9N0feTo8pm+5BAf/XOSgWH+bJnRg4ldS594U53V93bhpX6NebhNIH8evMDljJwSj+vy8nnrz2ME1nRk8k0zcC3RY13qkJmTx9O/7ic9W1fiseRr2WQUbNPl5fPMwki2nk7h/cEtePn+JrSpXavoOo67kx1P9aiPh7Mdn62985WLcnT5XM26u1FOVZZTLej4tLb6m4lWfdIna70ANJNShkgp60op60gpLf9/iJSw9RNt3dXQwSYJIfFqFqN/3M2y/fFMu68BXwxviZOdcS4EVlWPdqxNti6fJRElx8cv2H2Ok4nXeLV/02rx4Rjq78aHQ1qwIzqVEbN2kXwtm/jLmTz/20Ha/W8dzd78h/s+3cyQ73by99EEXhvQlKFtSy8S62RnwxPd67EtKoXdZ1L1jiE9W8fQ73cS/s46Ziw9eNux/BYpbDgIa23ykwnok+CjgUxDB2J2zm6G8xHQZRpYGzepSilZEXme3p9t4VD8FT4f1pJp9zU06PhsS9HY15V2IbX4Zde5otWiYlMz+Pjfk3Sp70mfUB8TR2g8Q8OD+OHRNkQlpdP/y630/Hgzfx66wITOdXjm3gaEeDhxOTOHl/o15rEu5Q+MG9W+Nl417Pl07albxthLKXlm0QEGf7udU4nXAO0i8KSf9nH4fBp9Qn1ZefACfT/fytg5e6pXoq/hC/Xvg4OLIN/46/Lqk7n+D9ghhNgNZBdulFJONVhUppZ7XVuhxcVXqy9hRNdz8pi+9CB/HbpIq2B3PnkkjLrllKFVbjWmY22mLDzA5lPJtKtTi8nzI7ASgv8Nal7tPiR7NvZh4eQOPLPoAPc09OLZXg31LntQnKOdNU91r8ebfx67ZYz9isgLrIi8gL2NFQO+2sbzvRqyN+Yyu86m8tnQljzUKoArmTks2H2OWVvOcP8XWxnWNpjnezcsUavIYrUaBb/9o82Gb9DLqIfWZxTNHrQl+g4DRUv3SCl/quxgzGYUzZ/TtIkKI5dAw95GO6yUkulLDrH8QDzTezfi8W51q92wx8qQo8un8wcbCPV3xcHGmn+PJTB/Qnu6NFDjwu9GVm4ePT7ehLWV4NeJHQj2cCIlPZten26mjqcz341uw6t/HOHfY4kAvP1gKGM6hpTYx5XMHD5fd5pfdsViYy0Y0MKfke2DaRHgxrGLV9l1JpXUjByGhgeVW1+/StHlwKeNIaQrDK30tFnuKBp9EvwBKaVRVmYyiwR/eCkseww6T4Nebxn10Iv3nuPFZYeZ2rM+z/W27FIDhvbp2lN8uV67sPVq/yZM7Gr5l42M4XB8GmPm7MbBxppfJrbns3WnWHs0kdXPdKG+dw2klPx1+CI5uvxy1++NSkrnx21nWRl5noycPOxsrMjRae1HaytBXr7kviY+PNG9rmWM3FnzEuz7EZ4/qV18rUR3m+D/B8QAf1Kyi6bSh0maPMGnRMGse8CnGYxbZdQZaMcuXGXQt9sJD6nJ/Antq81sVENJSMuix8eb6NvMl0+HhlW7rhlDOpFwldGz95Cdm8e1bB3Tezfk6Z4NKrSv9GwdKyLPczoxnda1a9KhTi2srATzd8by884YLmfmMqFzHV7q19gkSzRWmoTD8F0Xbehk+8crddd3m+DPlrJZ3m4kjRAiCJgP+KDNhJ0lpfyivNeYNMFLCT/2htQo+M9WcCu79VHZrmXl8sBX27iem8dfU7tWj35JI7iUkUNNJ1uV3A3gTHI6o2bvxtPFnuVPdjJI9dDrOXl88PcJ5u2IoUWgG1+NaEVtjzubmWtWvuuqLRY0cV2l7vZuJzpVtO6MDnheSrlfCFEDiBBCrJVSHqvg/gzr2B8QvwcGfm3U5A7w7aZoYlIz+e3xjiq5V6Jaal1Tg6nr5cLG6d0Bw5WGdrSz5s2BoXSo68GMpQfp/dkWujbwpHsjb3o29i71YvGdrLJldHXvgd3fQ16u0XoH9KlFU2rtdynl/PJeJ6W8CFws+P2aEOI4EACYX4LXZcPaN7SumZbGHTVz4cp15mw7y0Mt/WlXxwL6GpVqw1jzCfo286VZgCuztpxhw4kk1h1PwsZK8M2o1kWF9gBWHbrA6yuO8vXIVnSqZ4YX1H3DIC9HWxnOJ9QohzRKPXghRAjQCthdymOThRD7hBD7kpOT72S3lWfPLLgSC73f1tZWNKJP/j2FlDC9j7qoqihlCazpxH8fbMbWGT1Y91w3QgPcmLrwAPtitEuBW04l8+ziSC5l5PDiskNFM3XNSmE9mouHjHbI2yZ4KeWUYrdJQGtA7/FLQggXYBkwTUp5ywwHKeUsKWW4lDLcy8t4Ff2KZF6CLR9B/V5Gr/d+7MJVlh+IZ1znkGpRv11R7pYQgvreNZg7ri0B7o5MmLeXpRHx/OeXCOp5ufDj2HDiLl3no39OmjrUW3k2ABtH7YKrkVSk80zvevBCCFu05L5ASrm8AscyvK2fQPY1rfVuZO+tOY6rgy1Pda9v9GMrSlVWy9mOnya0w97WmulLDuLpYs/8Ce24t4kPj3aszU87Y4pa92bDylrrmkkwoxa8EOJPIcTKgtsq4CTwux6vE8CPwHEp5ad3H6oB5Odr5YCbDDT6WqsRsZfYejqFKT3rl7uij6IopQuq5cT8Ce24v7kv8ye0K1pDd0bfxvi7OTJj2SGyco1fHqBcvs21BG+kJRX1acF/DHxScHsP6CalfEmP13UGxgA9hRCRBbf7Kx6qASQchIwkaGT8UsA7orSiTY+El17gSVGU22vi58q3o9oQUmxhExd7G94d1IwzyRks2Vc5i7JXGr8WkJUGV84Z5XBljqIRQtQHfArrwRfb3lkIYS+ljC5vx1LKbYCZjlcqcLpgPGq9e41+6Mi4K9T3drltrW5FUe7cPQ29aBHoxk87Yxndobb5DJ30LVinIOEQ1Kxt8MOV14L/HCit7NvVgseqvqi14N8KXIx7cVdKSWTcFVoGuRv1uIpSXQghGNsxhKikdLZH6V/i2OC8m4CwMtpImvISvI+U8pbLvQXbQgwWkbFkXoL4vdroGSOLv3yd1IwcleAVxYAGhPnh4WzHvB0xpg7lBjsn8GxotJE05SV493Ieu/N6o+bmzEaQ+UYv3wmw/9xlAFoFuxv92IpSXdjbWDOiXTDrTyQSd8mMlrTwbWG0kTTlJfh9QohJN28UQkxEW7avaju9DhxrQkAbox86Mu4KjrbWNPKpYfRjK0p1MqpDMFZC8POuWIMfK1uXx8qDFxj5wy5e/ePwLQujFPFtDlfPQ4bhu47KK1UwDfhdCDGKGwk9HLADBhk4LsPKz4eoddrFVSPPXAUtwTcPcFO13hXFwPzcHOkb6suiPeeYdl8Dgy15uXDPOT7+5ySpGTl4ONuxIzqV4FpOTO5Wr5SgWmg/Ew5BvR4GiadQmRlGSpkopewEvIVWLjgGeEtK2VFKmWDQqAwt4ZA2PNIE3TPZujyOnr9KS9U9oyhGMb5zCFezdLy+4mjRMo6V6UpmDm+sPEqwhxM/TWjHnlfu4/7mvry/5gRbT5dSfsW3WII3MH1KFWyUUn5VcNtg8IiM4fRa7acJhkcev3iNnLx8WqkLrIpiFOEhtZh2XwOWRsTz8u+HKz3J/37gPDm6fN55qBn3NPTC2krw0ZAwGnjXYMrCA7f2/zvVAtdAo4ykMe5q0uYier1JhkcCRBZcYLXEFnxubi7x8fFkZWWZOhRFTw4ODgQGBmJra9nzMZ65twG6PMnXG6MQQtC5vgeH4tM4k5zBhM4hJdaYvRNSSn7dfY6wIHdC/d2Ktjvb2zDr0TYM/Ho7zy6OZMl/OpYci+/XAs7v02a0GnCMfvVL8LnXIX4fdHzSJIc/EHcFH1d7/Nyq/kCkm8XHx1OjRg1CQkLMZ2KJUiYpJampqcTHx1OnTkWXfagahBA837shufn5fL/5DAv3nMPOxgoXext2RKfw66QOFRq2HBF7mdNJ6Xz4cItbHqvt4cz0Po147Y8j7DpziY71PG482Oh+OLka4vZAcPu7OLPyVb+rfOcjID8Xanc2yeEj467QKqimSY5taFlZWXh4eKjkXkUIIfDw8Kg237iEELzUtzELJ3Xgz6e7cOTNPvw9rSseLnaMn7uH6OT0O97nr7vPUcPehgFhfqU+/kibQDxd7Jm5+aaJ/6GDwM4FDpS7rMZdq34JPnYHICDIcJ+aZbmUkUNsaqZFds8UUsm9aqlu75cQgo71PGge6IadjRXeNRz4eUJ7rITg0R/3sDfmEnl69tFfycxh1eGLPNQqoMzROQ621jzWpQ5bTiVz5HzajQfsXbQkf+R3rZqtgVTDBL9dW7nJ0d3ohy68oq5msCqK+QjxdOanCe24mpXLI9/tpM07a5m68ABRSeUn3uX7tYurI9oFl/u8UR2CqWFvw8xNWis+KzePN1ce5bNLHSA3A47etjhvhVWvBJ+Xq/V51e5k9EMnX8vm7VXHaOjjQutgy+yiMQdCCEaPHl10X6fT4eXlxYABAwA4ceIEHTt2xN7eno8//thUYSpmplmAG9te7Mk3I1tzXxMfNp1MYvy8vaRl5pb6/PRsHXN3nKVlkDtN/V3L3bergy1jOtZm9ZGL/HM0gYFfb2Pejhi+OOlOkn0Icv/PhjgloLol+IuHIDdTrwS/MzqV11ccYeycPXT/aCPj5+6pcG3p/HzJC0sPcjVLx1cjWmNnU73+2Y3J2dmZI0eOcP36dQDWrl1LQEBA0eO1atXiyy+/ZPr06aYKUTFTbo629G/hx8ePhDFvQjsS0rJ4fknkLcMqpZS8tOwQ5y9f56V+jfXa9/jOdbCztuLxnyO4lJHD/Ant+M899ZmV3hkRvweSDbMCVfUaRRO7XftZToLPys3j/TUnmLcjBhd7G0I8najvXYN1xxN568+jvDf41qvltzNvRwybTibz3wdDaeRbPcoTvPXnUY5dKK0YacU19XfljQduv1jx/fffz19//cWQIUNYuHAhI0aMYOvWrQB4e3vj7e3NX3/9VamxKZaldXBNXrm/CW/+eYzvt5zhie43ZqT+siuWVYcu8kKfRnSo61HOXm7wqmHPs70acvh8Gm880BTvGg50ru/Js3EDiD72K9ZrZxIy8vNKP49qluB3gEd9cPEu9eFTidd4asF+TielM65TCC/1a1y0cvwHf59g5qZo2tWpxaBWgXof8tiFq7y/5gT3NfFmTAfD139WYPjw4fz3v/9lwIABHDp0iAkTJhQleEXR19hOIeyNvcxH/5wgR5dPi0A3hIC3Vx2nRyMvnrinlDIExeTl5ZGQkEBcXBxxcXG0CAjgPyM7kZGRQYcOHYiLiyMhIYGv8vPZ9dQivB76H85Olbs2c/VJ8Pn5cG4nNB1Y5lOe/+0glzNz+GlCO+5pWHIS1PO9GhIRc5mXlx+hgXcNrmblsufsJaKTM8iXEiklAe6OPN+7UdGHwvWcPKYuOoC7ky0fDgmrViMW9GlpG0qLFi2IiYlh4cKF3H+/eS0iplQdQgg+eLgFF69c57N1p4q2B7g78skjYSQnJxUl77i4OHx9fRk6dCgATZo0ISoqCp1OV/S6sWPH0qlTJ5ycnKhVqxZNmjQhKCgIV1sdWW3CcHas/Lkx1SfBJx+HrCsQXHr3zMmEaxw+n8brA5rektwBbKyt+HJEK/p/uZUBX20DwEpo60LaWAmEEKw+nEBsaibfjmqNjbUV7/x1jKikdH55rD21nO0MeXbKTQYOHMj06dPZtGkTqalmtOCDUiVcvnyZ2NhY4uLiuI84QnQxZFvZ0/nhSbSrU4tuHdpw9OjREq/p06dPUYLv378/dnZ2BAUFFd1q19a+wQshWL16tVHOo/ok+Ngd2s8y+t+X7Y/HxkrwYEv/Mnfh6+bAD2PDWX88kTa1axIeUgtXhxtTvH/aEcMbK48yY+kheof6sGD3OR7vVpcuDSo2DVqpuAkTJuDu7k7z5s3ZtGmTqcNRzMzZs2c5efJkiRa4lJJ58+YBMGTIEDZsuFF6y8bGhq5du/LNh28DMHXqVLKzs0skcC+vGw1DcxmhZbAEL4SYAwwAkqSUzQx1HL3F7tAK/LjfOmZVl5fP7wfO06OxNx4u9uXupnVwzTKHOY7tFMK1rFw+/vcUf0Sep1mAK8/3blQp4St3JjAwkKlTp96yPSEhgfDwcK5evYqVlRWff/45x44dw9W1/KFuStVy5MgR9uzZQ3x8fFECT05OZt++fQgheOutt/jpp58AsLKyws/Pj4YNGxa9fsaMGTz55JMEBgYSFBSEj48P1tY3SotPnjzZ6OdUEYZswc8DvgYMOxdXH1Jq/e8hXUot7LP1dArJ17J5uLX+F0/L8lSP+mTm5LF4bxxfDG+lhkQaWXr6rdPNu3fvTvfu3QHw9fUlPj7eyFEplSE7OxsbGxusra2JjIxk9erVJVrgcXFxnD59Gg8PDxYtWsS7774LgI+PD4GBgQQHB5OVlYWjoyPPP/88kyZNIigoCD8/v1uKrfXp08cUp1jpDJbgpZRbhBAhhtr/HbkcA9cuQnDHUh9euj+emk629Gxc+uiaOyGEYEbfxkzv3Qgrq+pzUVVR7kZubi4XLlzA09MTZ2dnDh48yOzZs4mLiytqhSclJXHgwAFatmzJ7t27eeWVV/Dw8CAoKIjg4GA6d+5MXp42V+Wpp55iwoQJBAQEYG9/67fy5s2bG/sUTcLkffBCiMnAZIDg4PKn/FbYuV3az1ISfFpmLmuPJTKyXXCltrZVclcUTfHhgsHBwfj7+3P8+HFee+21opZ3QkICUkr+/PNPBgwYwPnz5/n555+L+rfbtGlDYGAgHh7auPMxY8YwZswYnMoYVujnV3rxr+rG5AleSjkLmAUQHh5e+cutgNY94+AGXrfOOlt1+AI5uvxK6Z5RlOomPz+f5OTkokRdv359mjdvTlxcHMOHDycuLo4LFy4Utay//fZbnnjiCfLz8zly5AhBQUH07du3KJG3aKFNJOzXrx9Xrlwp87hlJXalJJMneKM4twuCOoDVrS303/efp6GPC80C1EU2RSlOSsmlS5dK9HHXr1+f3r17k5GRQYsWLYiPjycnJ6foNS+//DLNmzenRo0a2NnZ0b179xIjTVq1agVAaGgoJ06cKPPY1WnOiCFZfoLPSIWUk9ByxC0PXUy7zr7Yyzzfq6H6g1KqHZ1OR25uLtevX2fWrFnExcURFBRUNELE19eXpKSkEq8ZM2YMvXv3xsnJia5du+Lj41MigdetWxcAd3d3Nm7caPRzUkoy5DDJhUB3wFMIEQ+8IaX80VDHK1Nc2f3vaw5ra4ff30L11ymWJzs7m+zsbHJycopuNjY2BAZq3ZEnTpwgKyuLlJQUHn/8caysrBg8eHBRgn/22WdxcHAokcB9fHwArYVdOGZcMV+GHEVza5PZFM7tBGt7bQ3Wm6w+fJHGvjWo5+VigsAUQxBCMGrUKH755RdAa6X6+fnRvn17Vq1axYIFC/jggw+QUlKjRg1mzpxJWFiYiaOumMzMTDIzM8nNzS1K4AANGjQAIDY2lqtXbxR8s7W1xdnZuei+v782qc/GxoZz587h5+eHjc2NlPDSSy8Z4zQUA7L8LppzuyCgNdiUHCqVkJZV1D2jWI7i5YIdHR1vKRdcp04dNm/eTM2aNVmzZg2TJ09m9+7dJoy4pPz8fITQSl+kp6eTnp5eogWu0+lo3rw5QgiSkpJISUkBtCRtZ2dXYkigv78/fn5+2NnZYWtri9VN16Bq1aoFQGJiIkFBQcY7ScVoLDvB52TChQPQacotD605chFQ3TMGs+YlSDhcufv0bQ793r/t08orF9yp041SFR06dDDqpCcpJbm5udjY2GBlZUV6ejqXLl0q0QLPzc0lLCwMW1tb0tLSuHjxItbW1tjZ2WFnZ4ezszNSSoQQ+Pn54evri52d3S3JG8DFRX0zre4sO8Gfj4B8XakFxv46pLpnLJW+5YJ//PFH+vXrVynHlFKi0+nIycnB3t4eGxsb0tPTSUxMLErchV0oTZo0wdnZmevXr5OSklKUvN3c3LCzsyu64O/j44Ovr2+JKfLFlTaBR1GKs+wEf24X2gLbbUtsVt0zRqBHS9tQ9CkXvHHjRn788Ue2bdt22/0VJu/CJO3o6Ii9vT2ZmZmcO3euaLuU2jSOevXqUbNmTfLz88nMzMTOzq5o2KCtrS12dlplUU9PTzw9PcscwVW8P1xRKsKy/4LO7QTvpuBYsjiY6p6xfOWVCz506BATJ05kzZo11KpVq6jlXTyBu7i4kJ2dzalTp8jNzSU/P7/o9cHBwXh7exd1izg7O1OzZs0S3SgArq6u5U6JV0NzFUOz3ASfna5VkGz9aInNWbl5LNxzTnXPWLjCcsFNmzZl3bp16HQ6kpOTiYmJYdiwYfzyyy/Uq1ePAwcOlEjeoI3/dnFxwcbGBicnp6LEXXgr7BpxcHCgcWP91uRUFFOw3AR/cg3orkPooBKb3151jFOJ6fw4NtxEgSmVqfhY79zcXKSUJCQkFJULPnDgADExMWRkZBAbG8u7775LamoqTz75JKCNWvn3339LJPDCyoLW1tbUq1f+smyKYs4sN8EfWQquARDUvmjTisjz2iIc99Tl3iY+JgxO0VdGRgZZWVklulBsbGwICQkB4PTp02RlZRU9f/v27SXuF5aJHTlyJHZ2dixfvrzUESeKYoksM8FnXoKo9dD+8aL6M1FJ6fzf8sOE167JdLUIh0kVDvMDuHbtGhkZGSWGCYI20gTg/PnzRZN1CocLFr/4WFiBtLD1fXPyLr7KjqJUN5aZ4E+sgvxcaD4EgPXHE3ntjyM42Frz1chW2FqrFpyhFI71trW1RQjB1atXSUtLK5HAdTodrVq1QghBamoqKSkpWFlZ3dLHDRRNwLGzsyt1uKBaiUlRymaZCf7IMqhVl/OOjXhz/j7WHkukgbcL344Ow8+t8lcury6klKSmphIXF0fjxo1xdHRk/fr1zJkzh7i4OF5//XUyMzORUtKiRQvs7OxIT08nKSmpKHkXDhcsbMUHBAQQGBiItbV1qaNKHA2w0ryiVBeWl+CvJcLZLeR0fJaHZ+4k7XouL/VrzITOddTyeeWQUpKWllZUFrZNmzb4+PiwdetW3njjjaKVdQr7t/fs2UPbtm25ePEiO3fuJCgoCHt7e3x8fEp0lfj6+uLn51fmkMCbl0pTFKXyWF6CP7YCZD7zr4WTcDWLZU90ok3t0hfJrk6uXbtWYvmzjh070rRpU/bv38+oUaOIi4sjIyOj6PnLli1j8ODBCCHIysqidevWPPjgg0VVBQtHl4wePZrRo0cDcPz48aJKhYXUBU1FMR3LS/BHlpLj0YQP9wuGtAmoFsn9+vXrJVaPj4uLo1OnTvTs2ZOzZ8/SqlUr0tLSSrzm888/p2nTptSqVYumTZvSp0+fEmVhmzZtCkCXLl3YsWOHKU5LUZS7ZFkJPvkUxO1mZc2J2NlYMaOvZYyWOXv2bInkHR8fT3h4OOPGjSMjI6PUolKvvfYaPXv2xMfHh9GjR5dI3oGBgUUVFkNCQli2bJmxT8lgrK2tad68OTqdjjp16vDzzz/j7u5+1/udN28e+/bt4+uvv777IIvp3r07Z86cITY2tqgb66GHHmLdunWkp6frvZ9x48YxYMAAhgwZclfPUSyLZSX4iHnkCxvev9iGKf3q413DwdQR6WXv3r2cOXOmRAJv0qQJ77zzDgBt2rTh8uXLRc93d3cvuvjo7OzMhx9+WGJlncDAQBwctHN3cnKq9KRkzhwdHYmMjARg7NixfPPNN7zyyiumDeo23N3d2b59O126dOHKlStcvHjR1CEpFsJyEnxuFjLyVzZbtcPV05/xneuYNJziY73Xr1/P4cOHS7TC/f39i1rOEydO5NChQ4BW4jUoKIjatWsX7Wv27NlF24OCgm5psb/wwgtGOqs7071791u2DR06lCeffJLMzMxSC4GNGzeOcePGkZKScktLc9OmTXd0/I4dOxb9u+7Zs4dnnnmGrKwsHB0dmTt3Lo0aNWLevHmsXLmSzMxMoqOjGTRoEB9++CEAc+fO5b333sPd3Z2wsLCi4ZsxMTFMmDCBlJQUvLy8mDt3LsHBwYwbNw5HR0cOHDhAUlISc+bMYf78+ezcuZP27duXuQLS8OHDWbRoEV26dGH58uUMHjyYo0ePAtrf0YwZM1izZg1CCF599VWGDRuGlJIpU6awdu1agoKCigqYAURERPDcc8+Rnp6Op6cn8+bNw89P1V2qjiwmwctjfyCyLjMntwfvjWpu0BEzUkouX75ctGDCX3/9xZYtW0okcBsbG6KiogD49NNPWb16dYnlzwrXrgStbK29vT1BQUG4ubndMuJk8ODBBjsXS5WXl8f69et57LHHAGjcuDFbt27FxsaGdevW8fLLLxd9wEZGRnLgwAHs7e1p1KgRU6ZMwcbGhjfeeIOIiAjc3Nzo0aNH0YLRU6ZMYezYsYwdO5Y5c+YwdepU/vjjDwAuX77Mzp07WblyJQMHDmT79u3Mnj2btm3bEhkZScuWLW+J9d5772XSpEnk5eWxaNEiZs2axdtvvw3A8uXLiYyM5ODBg6SkpNC2bVu6devGzp07OXnyJMeOHSMxMZGmTZsyYcIEcnNzmTJlCitWrMDLy4vFixfzyiuvMGfOHMP/oytmx2ISfMKGmWTn+3BPn8G0r+tR4f1IKbly5QpxcXGEhoZibW3NypUrWbZsWYkulNzcXLKzs7G2tubPP/9kzpw5BAYGEhQURJcuXahT58Y3iO+//x5HR0dq1apV6nDB8HDLrItTXovbycmp3Mc9PT3vuMUO2gXnli1bcv78eZo0aUKvXr0ASEtLY+zYsZw+fRohRNGMWdASrJubGwBNmzYlNjaWlJQUunfvXjQTdtiwYZw6dQqAnTt3snz5ckBbhHrGjBlF+3rggQcQQtC8eXN8fHyKqkmGhoYSExNTaoK3tramS5cuLFq0iOvXrxeVYQDYtm0bI0aMwNraGh8fH+655x727t3Lli1birb7+/vTs2dPAE6ePMmRI0eKzjsvL0+13qsxgyZ4IURf4AvAGpgtpTRIkfD9+3bSOi2S370n81jX8otDFQ4XLBxpUqNGDVatWsUXX3xRNBKlcLhgfHw8AQEBHD9+nI0bNxIUFESbNm146KGHCAoKQqfTYW1tzWeffca3335b5pDAm4cOKoZT2AefmZlJnz59+Oabb5g6dSqvvfYaPXr04PfffycmJqZE91HxmbPW1tbodLoKH79wX1ZWViX2a2VlVe5+hw8fzqBBg3jzzTcrfGzQGiihoaHs3LnzrvajWAaD9WMIIayBb4B+QFNghBCiaWUf50pmDif++opcbOg66AlOnz7N+vXrmTdvHnFxcQCsW7eO5s2b4+7ujqurK6GhofTt27eonzM7O5urV68SGhrKpEmT+Pjjj1m8eHHRNPgXX3yRc+fOsX37dhYtWsRHH33E1KlTi/4DOzo6qvHeZsbJyYkvv/ySTz75BJ1OR1paWtHIobL6wotr3749mzdvJjU1ldzcXJYsWVL0WKdOnVi0aBEACxYsoGvXrncdb9euXfm///s/RowYccv2xYsXk5eXR3JyMlu2bKFdu3Z069ataPvFixfZuHEjAI0aNSI5Obkowefm5hb9nSvVjyFb8O2AKCnlGQAhxCLgQeBYZR7E3TaP2hfX4bsgg0tvlWy9L1myhKCgINzd3alXrx7du3cvGmUSFBREs2bNAHj44Yd5+OGHKzMsxQy0atWKFi1asHDhQmbMmMHYsWN555136N+//21f6+fnx5tvvknHjh1xd3cv0bXy1VdfMX78eD766KOii6x3SwjB9OnTb9k+aNAgdu7cSVhYGEIIPvzwQ3x9fRk0aBAbNmygadOmBAcH07FjR0Cr2bN06VKmTp1KWloaOp2OadOmERoaetcxKlWPKFxmrNJ3LMQQoK+UcmLB/TFAeynl0zc9bzIwGSA4OLhNbGzsnR0oT8f5LfN5Z+4aghq3KjFUsHD6vGIcx48fL6oCqVQd6n2r2oQQEVLKUi/kmfwiq5RyFjALIDw8/M4/baxtCOgxgZk9JlR2aIqiKFWaITuOzwNBxe4HFmxTFEVRjMCQCX4v0EAIUUcIYQcMB1Ya8HiKGTBUl59iGOr9smwGS/BSSh3wNPAPcBz4TUqpLudbMAcHB1JTU1XSqCIK6/sXlrVQLI9B++CllKuB1YY8hmI+AgMDiY+PJzk52dShKHpycHBQ8zQsmMkvsiqWw9bWtsQMXkVRTEvNzlEURbFQKsEriqJYKJXgFUVRLJTBZrJWhBAiGbjDqaxFPIGUSgynKlDnbPmq2/mCOuc7VVtK6VXaA2aV4O+GEGJfWdN1LZU6Z8tX3c4X1DlXJtVFoyiKYqFUglcURbFQlpTgZ5k6ABNQ52z5qtv5gjrnSmMxffCKoihKSZbUglcURVGKUQleURTFQlX5BC+E6CuEOCmEiBJCvGTqeAxBCBEkhNgohDgmhDgqhHimYHstIcRaIcTpgp81TR1rZRNCWAshDgghVhXcryOE2F3wfi8uKEVtMYQQ7kKIpUKIE0KI40KIjpb+Pgshni34uz4ihFgohHCwtPdZCDFHCJEkhDhSbFup76vQfFlw7oeEEK0retwqneCNtbC3GdABz0spmwIdgKcKzvMlYL2UsgGwvuC+pXkGrdx0oQ+Az6SU9YHLwGMmicpwvgD+llI2BsLQzt1i32chRAAwFQiXUjYDrNHWjrC093ke0PembWW9r/2ABgW3ycDMih60Sid4ii3sLaXMAQoX9rYoUsqLUsr9Bb9fQ/tPH4B2rj8VPO0n4CGTBGggQohAoD8wu+C+AHoCSwueYlHnLIRwA7oBPwJIKXOklFew8PcZraqtoxDCBnACLmJh77OUcgtw6abNZb2vDwLzpWYX4C6E8KvIcat6gg8A4ordjy/YZrGEECFAK2A34COlvFjwUALgY6q4DORzYAaQX3DfA7hSsJgMWN77XQdIBuYWdEvNFkI4Y8Hvs5TyPPAxcA4tsacBEVj2+1yorPe10vJaVU/w1YoQwgVYBkyTUl4t/pjUxrtazJhXIcQAIElKGWHqWIzIBmgNzJRStgIyuKk7xgLf55poLdY6gD/gzK1dGRbPUO9rVU/w1WZhbyGELVpyXyClXF6wObHwq1vBzyRTxWcAnYGBQogYtK63nmj90+4FX+XB8t7veCBeSrm74P5StIRvye/zfcBZKWWylDIXWI723lvy+1yorPe10vJaVU/w1WJh74K+5x+B41LKT4s9tBIYW/D7WGCFsWMzFCnl/0kpA6WUIWjv6wYp5ShgIzCk4GmWds4JQJwQolHBpnuBY1jw+4zWNdNBCOFU8HdeeM4W+z4XU9b7uhJ4tGA0TQcgrVhXzp2RUlbpG3A/cAqIBl4xdTwGOscuaF/fDgGRBbf70fqk1wOngXVALVPHaqDz7w6sKvi9LrAHiAKWAPamjq+Sz7UlsK/gvf4DqGnp7zPwFnACOAL8DNhb2vsMLES7xpCL9k3tsbLeV0CgjQ6MBg6jjTCq0HFVqQJFURQLVdW7aBRFUZQyqASvKIpioVSCVxRFsVAqwSuKolgoleAVRVEslErwiqIoFkoleEVRFAv1/8+kMZhjFJ+HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(cumelast_1)), cumelast_1, label=\"M1\")\n",
    "plt.plot(range(len(cumelast_2)), cumelast_2, label=\"M2\")\n",
    "plt.plot([0, 100], [0, elast(test_pred, \"Y\", \"T\")], linestyle=\"--\", label=\"Random Model\", color=\"black\")\n",
    "plt.legend()\n",
    "plt.ylabel(\"Cumulative Elasticity\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10ac5842",
   "metadata": {},
   "source": [
    "As we can see, now $M1$ is much worse than $M2$, even though it has a higher $R^2$ on the test set. This shows that predictive power doesn't translate directly to a good causal model, even if we use random data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e6aab85",
   "metadata": {},
   "source": [
    "## Predictive Metric For Causal Inference\n",
    " \n",
    "This doesn't mean we can't come up with a way to correctly evaluate a causal model using predictive metrics. In order to do so, let's go back to our additive assumption about de DGP.\n",
    " \n",
    "$$\n",
    "Y = g(x) + f(t,x)\n",
    "$$\n",
    " \n",
    "To use a predictive metric, we need to somehow transform the outcome in order to remove the $g(x)$ component from it. That way, all the remaining predictive power will necessarily be used to learn the causal relationship.\n",
    " \n",
    "$$\n",
    "\\tilde{Y} = Y - g(x) = f(t,x)\n",
    "$$\n",
    " \n",
    "One way of doing that is by orthogonalizing $Y$. We can use any ML model to estimate $g$ and get out of fold residuals\n",
    " \n",
    "$$\n",
    "\\tilde{Y} = Y - \\hat{g}(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "aa28a557",
   "metadata": {},
   "outputs": [],
   "source": [
    "denoise_m = smf.ols(\"Y~\"+\n",
    "                    \"+\".join([f\"w{f}\" for f in range(n_heter)])+\n",
    "                    \"+\"+\n",
    "                    \"+\".join([f\"f{f}\" for f in range(n_features)]), data=train).fit()\n",
    "\n",
    "test_res = test.assign(Y_res = test[\"Y\"] - denoise_m.predict(test) + test[\"Y\"].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6302a6a",
   "metadata": {},
   "source": [
    "Once we do that, we can evaluate the power of each model in predicting the residualized outcome $\\tilde{Y}$. The predictive performance on this new outcome will directly translate to a better causal model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0c6ed700",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M1: -281.7950185092476\n",
      "M2: -24.04912481654577\n"
     ]
    }
   ],
   "source": [
    "print(\"M1:\", r2_score(test_res[\"Y_res\"], m1.predict(test_res)))\n",
    "print(\"M2:\", r2_score(test_res[\"Y_res\"], m2.predict(test_res)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f12479",
   "metadata": {},
   "source": [
    "The downside of this approach is that it depends on how well you can estimate $g(x)$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "causal-glory",
   "language": "python",
   "name": "causal-glory"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}