Chapter 19 - Causal Inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 19. Causal Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we discuss causation.  Roughly speaking \"$X$ causes $Y$\" means that changing the value of $X$ will change the distribution of $Y$.  When $X$ causes $Y$, $X$ and $Y$ will be associated but the reverse is not, in general, true.\n",
    "\n",
    "We will consider two frameworks for discussing causation.  The first uses notation of **counterfactual** random variables.  The second, used in the next chapter, uses **directed acyclic graphs**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19.1 The Counterfactual Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $X$ is a binary treatment variable where $X = 1$ means \"treated\" and $X = 0$ means \"not treated\".\n",
    "\n",
    "Let $Y$ be some outcome variable such as the presence or absence of disease.  To distinguish the statement \"$X$ is associated with $Y$\" from the statement \"$X$ causes $Y$\" we need to enrich our probabilistic vocabulary.  We will decompose the response $Y$ into a more fine-grained object.\n",
    "\n",
    "We introduce two new random variables $(C_0, C_1)$ called **potential outcomes** with the following interpretation: $C_0$ is the outcome if the subject is not treated ($X = 0$) and $C_1$ is the outcome if the subject is treated ($X = 1$).  Then,\n",
    "\n",
    "$$ Y = \\begin{cases}\n",
    "C_0 & \\text{if } X = 0 \\\\\n",
    "C_1 & \\text{if } X = 1\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can express the relationship between $Y$ and $(C_0, C_1)$ more succintly by\n",
    "\n",
    "$$ Y = C_X $$\n",
    "\n",
    "This equation is called the **consistency relationship**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a toy dataset to make the relationship clear:\n",
    "\n",
    "$$\n",
    "\\begin{array}{cccc}\n",
    "X & Y & C_0 & C_1 \\\\\n",
    "\\hline\n",
    "0 & 4 & 4 & * \\\\\n",
    "0 & 7 & 7 & * \\\\\n",
    "0 & 2 & 2 & * \\\\\n",
    "0 & 8 & 8 & * \\\\\n",
    "\\hline\n",
    "1 & 3 & * & 3 \\\\\n",
    "1 & 5 & * & 5 \\\\\n",
    "1 & 8 & * & 8 \\\\\n",
    "1 & 9 & * & 9\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The asterisks denote unobserved values.  When $X = 0$ we don't observe $C_1$ in which case we say that $C_1$ is a **counterfactual** since it is the outcome you would have had if, counter to the fact, you had been treated ($X = 1$).  Similarly, when $X = 1$ we don't observe $C_0$ and we say that $C_0$ is counterfactual."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that there are four types of subjects:\n",
    "\n",
    "$$\n",
    "\\begin{array}{lcc}\n",
    "\\text{Type} & C_0 & C_1 \\\\\n",
    "\\hline\n",
    "\\text{Survivors}       & 1 & 1 \\\\\n",
    "\\text{Responders}      & 0 & 1 \\\\\n",
    "\\text{Anti-responders} & 1 & 0 \\\\\n",
    "\\text{Doomed}          & 0 & 0\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Think of all of the potential outcomes $(C_0, C_1)$ as hidden variables that contain all the relevant information about the subject."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the **average causal effect** or **average treatment effect** to be\n",
    "\n",
    "$$ \\theta = \\mathbb{E}(C_1) - \\mathbb{E}(C_0) $$\n",
    "\n",
    "The parameter $\\theta$ has the following interpretation:  $\\theta$ is the mean if everyone were treated ($X = 1$) minus the mean if everyone were not treated ($X = 0$).  There are other ways of measuring the causal effect.  For example, if $C_0$ and $C_1$ are binary, we define the **causal odds ratio**\n",
    "\n",
    "$$ \\frac{\\mathbb{P}(C_1 = 1)}{\\mathbb{P}(C_1 = 0)} \\div \\frac{\\mathbb{P}(C_0 = 1)}{\\mathbb{P}(C_0 = 0)}$$\n",
    "\n",
    "and the **causal relative risk**\n",
    "\n",
    "$$ \\frac{\\mathbb{P}(C_1 = 1)}{\\mathbb{P}(C_0 = 1)} $$\n",
    "\n",
    "The main ideas will be the same whatever causal effect we use.  For simplicity, we shall work with the average causal effect $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the **association** to be\n",
    "\n",
    "$$ \\alpha = \\mathbb{E}(Y | X = 1) - \\mathbb{E}(Y | X = 0)$$\n",
    "\n",
    "Again we could use the odds ratio or other summaries if we wish."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 19.1 (Association is not equal to Causation)**.  In general, $\\theta \\neq \\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 19.3**.  Suppose we randomly assign subjects to treatment and that $\\mathbb{P}(X = 0) > 0$ and $\\mathbb{P}(X = 1) > 0$.  Then $\\alpha = \\theta$.  Hence, any consistent estimator of $\\alpha$ is a consistent estimator of $\\theta$.  In particular, a consistent estimator is\n",
    "\n",
    "$$ \\hat{\\theta} = \\hat{\\mathbb{E}}(Y | X = 1) - \\hat{\\mathbb{E}}(Y | X = 0) = \\overline{Y}_1 - \\overline{Y}_0 $$\n",
    "\n",
    "is a consistent estimator of $\\theta$, where\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\hat{Y}_1 = \\frac{1}{n_1} \\sum_{i: X_i = 1} Y_i\n",
    "&\n",
    "\\hat{Y}_0 = \\frac{1}{n_0} \\sum_{i: X_i = 0} Y_i \\\\\n",
    "n_1 = \\sum_{i=1}^n X_i\n",
    "&\n",
    "n_0 = \\sum_{i=1}^n (1 - X_i)\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Since $X$ is randomly assigned, $X$ is independent of $(C_0, C_1)$.  Hence,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\theta &= \\mathbb{E}(C_1) - \\mathbb{E}(C_0) \\\\\n",
    "&= \\mathbb{E}(C_1 | X = 1) - \\mathbb{E}(C_0 | X = 0) \\\\\n",
    "&= \\mathbb{E}(Y | X = 1) - \\mathbb{E}(Y | X = 0) \\\\\n",
    "&= \\alpha\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The consistency follows from the law of large numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $Z$ is a covariate, we define the **conditional causal effect** by\n",
    "\n",
    "$$ \\theta_z = \\mathbb{E}(C_1 | Z = z) - \\mathbb{E}(C_0 | Z = z) $$\n",
    "\n",
    "In a randomized experiment, $\\theta_z = \\mathbb{E}(Y | X = 1, Z = z) - \\mathbb{E}(Y | X = 0, Z = z)$ and we can estimate the conditional causal effect using appropriate sample averages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary**\n",
    "\n",
    "- Random variables: $(C_0, C_1, X, Y)$\n",
    "- Consistency relationship: $Y = C_X$\n",
    "- Causal Effect: $\\theta = \\mathbb{E}(C_1) - \\mathbb{E}(C_0)$\n",
    "- Association: $\\alpha = \\mathbb{E}(Y | X = 1) - \\mathbb{E}(Y | X = 0)$\n",
    "- Random Assignment:  $(C_0, C_1) \\text{ ⫫ } X \\Longrightarrow \\theta = \\alpha$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19.2 Beyond Binary Treatments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $X \\in \\mathcal{X}$. The counterfactual vector $(C_0, C_1)$ now becomes the **counterfactual process**\n",
    "\n",
    "$$ \\{ C(x) : x \\in \\mathcal{X} \\} $$\n",
    "\n",
    "where $C(x)$ is the outcome a subject would have if subjected to treatment $x$.  The observed response is given by the consistency relation\n",
    "\n",
    "$$ Y \\equiv C(X) $$\n",
    "\n",
    "The **causal regression function** is \n",
    "\n",
    "$$ \\theta(x) = \\mathbb{E}(C(x)) $$\n",
    "\n",
    "The regression function, which measures association, is $r(x) = \\mathbb{E}(Y | X = x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 19.4**.  In general, $\\theta(x) \\neq r(x)$. However, when $X$ is randomly assigned, $\\theta(x) = r(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19.3 Observational Studies and Confounding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A study in which treatment (or exposure) is not randomly assigned is called an **observational study**. In these studies, subjects select their own value of the exposure $X$.  Association and causation could be quite different.  This discrepancy occurs in non-randomized studies because the potential outcome $C$ is not independent of treatment $X$.\n",
    "\n",
    "However, suppose we could find groupings of subjects such that, within each group, $X$ and $\\{ C(x) : x \\in \\mathcal{X} \\}$ are independent.  This would happen if the subjects are very similar within groups.  For example, suppose we find people who are very similar in age, gender, educational background and ethnic background.  Among those people we might feel it is reasonable to assume that the choice of $X$ is essentially random.  These other variables are called **confounding variables**.  If we denote these variables collectively as $Z$, then we can express this idea by saying that\n",
    "\n",
    "$$ \\{ C(x) : x \\in \\mathcal{X} \\} \\text{ ⫫ } X | Z$$\n",
    "\n",
    "If this holds and we observe $Z$ then we say there is **no unmeasured confounding**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 19.17**.  Suppose that $ \\{ C(x) : x \\in \\mathcal{X} \\} \\text{ ⫫ } X | Z$.  Then,\n",
    "\n",
    "$$ \\theta(x) = \\int \\mathbb{E}(Y | X = x, Z = z) d F_Z(z) dz $$\n",
    "\n",
    "If $\\hat{r}(x, z)$ is a consistent estimate of the regression function $\\mathbb{E}(Y | X = x, Z = z)$, then a consistent estimate of $\\theta(x)$ is\n",
    "\n",
    "$$ \\hat{\\theta}(x) = \\frac{1}{n} \\sum_{i=1}^n \\hat{r}(x, Z_i) $$\n",
    "\n",
    "In particular, if $r(x, z) = \\beta_0 + \\beta_1 x + \\beta_2 z$ is linear, then a consistent estimate of $\\theta(x)$ is\n",
    "\n",
    "$$ \\hat{\\theta}(x) = \\hat{\\beta}_0 + \\hat{\\beta}_2 x + \\hat{\\beta}_2 \\overline{Z}_n $$\n",
    "\n",
    "where $(\\hat{\\beta}_0, \\hat{\\beta}_1, \\hat{\\beta}_2)$ are the least squares estimators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Epidemiologists call this definition of $\\theta(x)$ the **adjusted treatment effect**.  The process of computing adjusted treatment effects is called **adjusting (or controlling) for confounding**.  The selection of what confounders $Z$ to measure and control for requires scientic insight.  Even after adjusting for confounders, we cannot be sure that there are not other confounding variables that we missed.  This is why observational studies must be treated with healthy skepticism.  Results from observational studies start to become believable when: (i) the results are replicated in many studies, (ii) each of the studies controlled for plausible confounding variables, (iii) there is a plausible scientic explanation for the existence of a causal relationship.\n",
    "\n",
    "A good example is smoking and cancer.  Numerous studies have shown a relationship between smoking and cancer even after adjusting for many confouding variables.  Moreover, in laboratory studies, smoking has been shown to damage lung cells.  Finally, a causal link between smoking and cancer has been found in randomized animal studies.  It is this collection of evidence over many years that makes this a convincing case.  One single observational study is not, by itself, strong evidence.  Remember that when you read the newspaper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19.4 Simpson's Paradox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X$ be a binary treatment variable, $Y$ a binary outcome and $Z$ a third binary variable such as gender.  Suppose the joint distribution of $X$, $Y$, $Z$ is\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccccc}\n",
    "       & Y = 1 & Y = 0 & Y = 1 & Y = 0 \\\\\n",
    "\\hline\n",
    "X = 1 & .1500 & .2250 & .1000 & .0250 \\\\\n",
    "X = 0 & .0375 & .0875 & .2625 & .1125 \\\\\n",
    "\\hline\n",
    "      & Z = 1 & & Z = 0 &\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The marginal distribution for $(X, Y)$ is\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc|c}\n",
    "       & Y = 1 & Y = 0 & \\\\\n",
    "\\hline\n",
    "X = 1 & .25 & .25 & .50 \\\\\n",
    "X = 0 & .30 & .20 & .50 \\\\\n",
    "\\hline\n",
    "      & .55 & .45 & 1\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now,\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(Y = 1 | X = 1) - \\mathbb{P}(Y = 1 | X = 0) \n",
    "&= \\frac{\\mathbb{P}(Y = 1, X = 1)}{\\mathbb{P}(X = 1)} - \\frac{\\mathbb{P}(Y = 1, X = 0)}{\\mathbb{P}(X = 0)} \\\\\n",
    "&= \\frac{.25}{.50} - \\frac{.30}{.50} \\\\\n",
    "&= -0.1\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We might naively interpret this to mean that the treatment is bad for you since $\\mathbb{P}(Y = 1 | X = 1) < \\mathbb{P}(Y = 1 | X = 0) $.  Furthermore, among men ($Z = 1$),\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(Y = 1 | X = 1, Z = 1) &- \\mathbb{P}(Y = 1 | X = 0, Z = 1) \\\\\n",
    "&= \\frac{\\mathbb{P}(Y = 1, X = 1, Z = 1)}{\\mathbb{P}(X = 1, Z = 1)} - \\frac{\\mathbb{P}(Y = 1, X = 0, Z = 1)}{\\mathbb{P}(X = 0, Z = 1)} \\\\\n",
    "&= \\frac{.15}{.3570} - \\frac{.0375}{.1250} \\\\\n",
    "&= 0.1\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, among women ($Z = 0$),\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{P}(Y = 1 | X = 1, Z = 0) &- \\mathbb{P}(Y = 1 | X = 0, Z = 0) \\\\\n",
    "&= \\frac{\\mathbb{P}(Y = 1, X = 1, Z = 0)}{\\mathbb{P}(X = 1, Z = 0)} - \\frac{\\mathbb{P}(Y = 1, X = 0, Z = 0)}{\\mathbb{P}(X = 0, Z = 0)} \\\\\n",
    "&= \\frac{.1}{.1250} - \\frac{.2625}{.3750} \\\\\n",
    "&= 0.1\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To summarize, it seems we have:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\text{Mathematical Statement} & \\text{English Statement?} \\\\\n",
    "\\hline\n",
    "\\mathbb{P}(Y = 1 | X = 1) < \\mathbb{P}(Y = 1 | X = 0) & \\text{treatment is harmful} \\\\\n",
    "\\mathbb{P}(Y = 1 | X = 1, Z = 1) > \\mathbb{P}(Y = 1 | X = 0, Z = 1) & \\text{treatment is beneficial to men} \\\\\n",
    "\\mathbb{P}(Y = 1 | X = 1, Z = 0) > \\mathbb{P}(Y = 1 | X = 0, Z = 0) & \\text{treatment is beneficial to women}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Clearly something is amiss.  We can't have a treatment that is good for men, good for women, but bad overall.  The problem is with the set of English statements.  **The inequality $\\mathbb{P}(Y = 1 | X = 1) < \\mathbb{P}(Y = 1 | X = 0)$ does not mean that the treatment is harmful**.\n",
    "\n",
    "The phrase \"treatment is harmful\" should be written as $\\mathbb{P}(C_1 = 1) < \\mathbb{P}(C_0 = 1)$.  The phrase \"treatment is harmful for men\" should be written $\\mathbb{P}(C_1 = 1 | Z = 1) < \\mathbb{P}(C_0 = 1 | Z = 0)$.  The three mathematical statements in the table are not contradictory, but the english statements are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now show that a real Simpson's paradox cannot happen, that is, there cannot be a treatment that is beneficial for men and women but harmful overall.  Suppose treatment is beneficial for each subgroup.  Then\n",
    "\n",
    "$$ \\mathbb{P}(C_1 = 1 | Z = z) > \\mathbb{P}(C_0 = 1 | Z = z) $$\n",
    "\n",
    "for all $z$.  It then follows that\n",
    "\n",
    "$$\\mathbb{P}(C_1 = 1) = \\sum_z \\mathbb{P}(C_1 = 1 | Z = z) \\mathbb{P}(Z = z) > \\sum_z \\mathbb{P}(C_1 = 0 | Z = z) \\mathbb{P}(Z = z) = \\mathbb{P}(C_1 = 0) $$\n",
    "\n",
    "so $\\mathbb{P}(C_1 = 1) > \\mathbb{P}(C_1 = 0)$ and the treatment is overall beneficial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19.6 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 19.6.1**.  Create an example like Example 19.2 in which $\\alpha > 0$ and $\\theta < 0$.\n",
    "\n",
    "Example:\n",
    "\n",
    "$$\n",
    "\\begin{array}{cccc}\n",
    "X & Y & C_0 & C_1 \\\\\n",
    "\\hline\n",
    "0 & 0 & 0 & 0* \\\\\n",
    "0 & 0 & 0 & 0* \\\\\n",
    "0 & 0 & 0 & 0* \\\\\n",
    "0 & 0 & 0 & 0* \\\\\n",
    "\\hline\n",
    "1 & 1 & 1* & 1 \\\\\n",
    "1 & 1 & 1* & 1 \\\\\n",
    "1 & 1 & 1* & 1 \\\\\n",
    "1 & 1 & 1* & 1\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\theta &= \\mathbb{E}(C_1) - \\mathbb{E}(C_0) = 1/2 - 1/2 = 0 \\\\\n",
    "\\alpha &= \\mathbb{E}(Y | X = 1) - \\mathbb{E}(Y | X = 0) = 1 - 0 = 1\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**. Consider the following data for the whole population:\n",
    "\n",
    "$$\n",
    "\\begin{array}{cccc}\n",
    "X & Y & C_0 & C_1 \\\\\n",
    "\\hline\n",
    "0 & 0 & 0 & 0* \\\\\n",
    "0 & 1 & 1 & 0* \\\\\n",
    "\\hline\n",
    "1 & 1 & 1* & 1\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\theta &= \\mathbb{E}(C_1) - \\mathbb{E}(C_0) = 1/3 - 2/3 = -1/3 \\\\\n",
    "\\alpha &= \\mathbb{E}(Y | X = 1) - \\mathbb{E}(Y | X = 0) = 1 - 1/2 = 1/2\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 19.6.2**.  Prove Theorem 19.4.\n",
    "\n",
    "In general, $\\theta(x) \\neq r(x)$. However, when $X$ is randomly assigned, $\\theta(x) = r(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\theta(x) &= \\mathbb{E}(C(x)) \\\\\n",
    "r(x) &= \\mathbb{E}(Y | X = x)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We can show the inequality by providing an example where it holds:  \n",
    "\n",
    "Let $C_i(x) \\equiv i$ for each population sample $i$ be constant.  Then $\\theta(x) = \\frac{1}{n} \\sum_i C_i(x) = n (n+1) / 2$ is a constant that does not depend on $x$.  On the other hand, if we have samples $(X_i, Y_i) = (i, i)$, then $r(x) = x$, which does depend on $x$, and the inequality holds.\n",
    "\n",
    "Now, when $X$ is randomly assigned, we have:\n",
    "\n",
    "$$\n",
    "\\theta(x) = \\mathbb{E}(C(x)) = \\mathbb{E}(C(x) | X = x) = \\mathbb{E}(Y | X = x) = r(x)\n",
    "$$\n",
    "\n",
    "where $\\mathbb{E}(C(x)) = \\mathbb{E}(C(x) | X = x)$ since $X \\text{ ⫫  } \\{ C(x) : x \\in \\mathcal{X} \\}$ (i.e. $X$ is randomly assigned) and $\\mathbb{E}(C(x) | X = x) = \\mathbb{E}(Y | X = x)$ since $Y \\equiv C(X)$ by definition.\n",
    "\n",
    "Note that this demonstration is analogous to the demonstration for the binary case in Theorem 19.3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 19.6.3**.  Suppose you are given data $(X_1, Y_1), \\dots, (X_n, Y_n)$ from an observational study, where $X_i \\in \\{0, 1\\}$ and $Y_i \\in \\{0, 1\\}$.  Although it is not possible to estimate the causal effect $\\theta$, it is possible to put bounds on $\\theta$.  Find upper and lower bounds on $\\theta$ that can be consistently estimated from the data.  Show that the bounds have width 1.  \n",
    "\n",
    "Hint:  note that $\\mathbb{E}(C_1) = \\mathbb{E}(C_1 | X = 1) \\mathbb{P}(X = 1) + \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\theta &= \\mathbb{E}(C_1) - \\mathbb{E}(C_0) \\\\\n",
    "&= \\mathbb{E}(C_1 | X = 1) \\mathbb{P}(X = 1) + \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0)\n",
    " - \\Big( \\mathbb{E}(C_0 | X = 1) \\mathbb{P}(X = 1) + \\mathbb{E}(C_0 | X = 0) \\mathbb{P}(X = 0) \\Big) \\\\\n",
    "&= \\mathbb{E}(C_1 | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(C_0 | X = 0) \\mathbb{P}(X = 0)\n",
    " + \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0) - \\mathbb{E}(C_0 | X = 1) \\mathbb{P}(X = 1) \\\\\n",
    "&= \\Big( \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) \\Big)\n",
    " + \\Big( \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0) - \\mathbb{E}(C_0 | X = 1) \\mathbb{P}(X = 1) \\Big)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first terms in this expression can be estimated: $Y \\equiv C_X$, so $\\mathbb{E}(C_x | X = x) = \\mathbb{E}(Y | X = x)$.  On the other hand, the unobserved terms $\\mathbb{E}(C_1 | X = 0)$ and $\\mathbb{E}(C_0 | X = 1)$ cannot be estimated from data.  Since $C_j \\in \\{0, 1\\}$, however, we can bound these expected values to the interval $[0, 1]$.  This means that:\n",
    "\n",
    "$$ 0 \\cdot \\mathbb{P}(X = 0) - 1 \\cdot \\mathbb{P}(X = 1) \\leq \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0) - \\mathbb{E}(C_0 | X = 1) \\mathbb{P}(X = 1) \\leq 1 \\cdot \\mathbb{P}(X = 0) - 0 \\cdot \\mathbb{P}(X = 1)$$\n",
    "\n",
    "or, simplifying,\n",
    "\n",
    "$$ - \\mathbb{P}(X = 1) \\leq \\mathbb{E}(C_1 | X = 0) \\mathbb{P}(X = 0) - \\mathbb{E}(C_0 | X = 1) \\mathbb{P}(X = 1) \\leq \\mathbb{P}(X = 0) $$\n",
    "\n",
    "which is a bound of width 1 (since $\\mathbb{P}(X = 0) + \\mathbb{P}(X = 1) = 1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we have\n",
    "\n",
    "$$ \\Big( \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) \\Big) - \\mathbb{P}(X = 1)\n",
    " \\leq \\theta \\leq\n",
    " \\Big( \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) \\Big) + \\mathbb{P}(X = 0)\n",
    "$$\n",
    "\n",
    "Given $Y$ is binary, we have:\n",
    "\n",
    "$$ \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) =  \\mathbb{P}(Y = 1 | X = 1) \\mathbb{P}(X = 1) =  \\mathbb{P}(X = 1, Y = 1) \\\\\n",
    "\\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) =  \\mathbb{P}(Y = 0 | X = 0) \\mathbb{P}(X = 0) =  \\mathbb{P}(X = 0, Y = 0) $$\n",
    "\n",
    "and so the estimates for the bounds are:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Big( \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) \\Big) - \\mathbb{P}(X = 1)\n",
    "&\\approx \\frac{1}{n} \\sum_{i = 1}^n \\left( I(X_i = 1, Y_i = 1) - I(X_i = 0, Y_i = 0) - I(X_i = 1)\\right) \\\\\n",
    "\\Big( \\mathbb{E}(Y | X = 1) \\mathbb{P}(X = 1) - \\mathbb{E}(Y | X = 0) \\mathbb{P}(X = 0) \\Big) + \\mathbb{P}(X = 0)\n",
    "&\\approx \\frac{1}{n} \\sum_{i = 1}^n \\left( I(X_i = 1, Y_i = 1) - I(X_i = 0, Y_i = 0) + I(X_i = 0)\\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 19.6.4**.  Suppose that $X \\in \\mathbb{R}$ and that, for each subject $i$, $C_i(x) = \\beta_{1i}x$.  Each subject has their own slope $\\beta_{1i}$.  Construct a joint distribution on $(\\beta_1, X)$ such that $\\mathbb{P}(\\beta_1 > 0) = 1$ but $\\mathbb{E}(Y | X = x)$ is a decreasing function of $x$, where $Y = C(X)$.  Interpret.\n",
    "\n",
    "Hint: Write $f(\\beta_1, x) = f(\\beta_1)f(x | \\beta_1)$.  Choose $f(x | \\beta_1)$ so that when $\\beta_1$ is large, $x$ is small and when $\\beta_1$ is small $x$ is large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "Let's construct a very simple joint distribution to develop this intuition.  Let $\\beta_1$ have a probability mass function with mass 1/2 on values 1/2 and 2.  Now let's pick points $(X, \\beta_1)$ with the following mass distribution:\n",
    "\n",
    "$$\\begin{array}{c|cc}\n",
    "\\mathbb{P} & X = 1 & X = 2 \\\\\n",
    "\\hline\n",
    "\\beta_1 = 1/2 & 0 & 0.5 \\\\\n",
    "\\beta_1 = 2 & 0.5 & 0 \\\\\n",
    "\\end{array}\n",
    "\\quad \\text{or} \\quad\n",
    "\\begin{array}{c|cc}\n",
    "\\mathbb{P} & X = 1 & X = 2 \\\\\n",
    "\\hline\n",
    "Y = 1 & 0 & 0.5 \\\\\n",
    "Y = 2 & 0.5 & 0 \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "For this toy example, $\\mathbb{E}(Y | X = x) = 3 - x$ is decreasing in $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAHgCAYAAABJt8A9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1hW5ePH8fdhCShuTX+5slyARUpDU9PcmYmjTM2RIC4cOXJVaoYjzay03ImKE7WSHOXISSammWmucqWWWxSQdX5/4NcyUUF4OPDweV0Xl3Ge59znY34vv59u7ue+DdM0ERERERERcLA6gIiIiIhIVqFyLCIiIiJyk8qxiIiIiMhNKsciIiIiIjepHIuIiIiI3KRyLCIiIiJyk5PVAf6tcOHCZpkyZayOISIiIiJ2bNeuXedN0yyS0mtZqhyXKVOGyMhIq2OIiIiIiB0zDOP43V7TsgoRERERkZtUjkVEREREblI5FhERERG5KUutOU5JfHw8p06dIjY21uoo8h+urq6UKFECZ2dnq6OIiIiIZIgsX45PnTqFh4cHZcqUwTAMq+PITaZpcuHCBU6dOsUjjzxidRwRERGRDJHll1XExsZSqFAhFeMsxjAMChUqpBl9ERERsStZvhwDKsZZlP5cRERExN5ki3IsIiIiIpIZVI7v49ixY3h7e6fpnjlz5nD69OkHet7SpUvx8vLCwcEhww5Eud+YVatW5erVqzRp0oSKFSvi5eXF4MGDM+TZIiIiItmJyrENpKcce3t7s3z5cmrVqpVhee415rFjx3j44YdxcnJiwIAB/Pbbb+zevZtt27axevXqDMsgIiIikh1k+d0q/q1v3w3s2fN3ho7p41OUSZNeuOd7EhIS6NixI7t376Z8+fLMnTsXd3d3du3aRb9+/bh27RqFCxdmzpw5bNu2jcjISNq1a4ebmxsRERGMHz+elStXEhMTQ/Xq1Zk2bdpd1+tWqlQpQ39/9xtz9erVNGrUCHd3d+rUqQOAi4sLVapU4dSpUxmeRURERCQr08xxKhw8eJDAwED27t1L3rx5+eyzz4iPj6dXr16EhYWxa9cuOnfuzLBhw2jVqhW+vr6EhoayZ88e3NzcCAoKYufOnezbt4+YmBjCw8PTlScqKgofH58Uv/bv35+msdasWUOjRo1uu3b58mVWrlxJ3bp105VTREREJLvJVjPH95vhtZWSJUvy3HPPAfD666/zySef0KhRI/bt20f9+vUBSExMpHjx4inev3HjRj744AOio6O5ePEiXl5eNG3a9IHzeHh4sGfPnge+/3/i4uI4deoUZcuWvXUtISGBNm3a0Lt379uui4iIiOQE2aocW+W/SyAMw8A0Tby8vIiIiLjnvbGxsfTo0YPIyEhKlizJiBEj0r03cFRUFDVr1kzxtQULFuDp6ZmqcbZs2UKNGjVuuxYYGEi5cuXo27dvujKKiIiIZEc2LceGYbwJBAAm8Avwhmma2e7UiBMnThAREUG1atVYuHAhNWrUoEKFCpw7d+7W9fj4eA4dOoSXlxceHh5ERUUB3CrChQsX5tq1a4SFhdGqVat05cmomeM1a9bQuHHjW9+//fbbXLlyhZkzZ6Z7bBEREZHsyGZrjg3DeBjoDfiapukNOAKv2ep5tlSpUiVCQkJ4/PHHuXjxIt27d8fFxYWwsDAGDRrEE088gY+PD9u3bwegU6dOdOvWDR8fH3LlykWXLl2oXLkyfn5+PPXUU/d81ooVKyhRogQRERE0adKEhg0bpjv/3cb8/vvvef7554HkY7qDg4PZv38/VapUwcfHRyVZREREchzDNE3bDJxcjn8AngCuAl8Cn5im+e3d7vH19TX/uw/vgQMHbLKDQ0536tQpunTpku7t2vTnIyIiIg/i1MVTFPYojKuza6Y/2zCMXaZp+qb0ms1mjk3T/BOYAJwAzgBX7lWMJXOVKFFC+xiLiIhIpktKSmLq91PxHO7JmFVjrI5zB1suqygANAMeAf4PyG0YxuspvC/QMIxIwzAiz507Z6s4WU7Pnj3v2Ibtiy++sDqWiIiIiM0c/uswL3z4At1Du/PMI8/QqXonqyPdwZYfyKsH/GGa5jkAwzCWA9WB+f9+k2ma04HpkLyswoZ5spQpU6ZYHUFEREQkUyQkJvDRuo9496t3yeWUi1kdZ/HGc2/c9VA0K9myHJ8AnjUMwx2IAeoCkfe+RURERETsyc8nf8Y/xJ9dx3fh5+PHlHZT+L/8/2d1rLuyWTk2TXOHYRhhwE9AArCbmzPEIiIiImLfbsTf4P1v3mfsmrEUdC/Ikq5LaFW1VZacLf43m+5zbJrmcGC4LZ8hIiIiIllLxNEI/EP8OXDmAB2qdWDiqxMplKeQ1bFSRSfkiYiIiEiGuH7jOsNWDOOTDZ9QokAJVvVeRePKje9/YxZis90qrBIaCmXKgIND8q+hoekfMyYmhueff57ExEQiIyPx9vYmLi4OgKNHj1K2bFmuXr2a4r2maVKjRo3btk1bsmQJjRo1Ii4ujlq1apGQkJD+kMCkSZOYO3cuAAMHDqRixYo8/vjjNG/enMuXL9/3/s6dO1O0aFG8vb1vuz5gwAA2bNiQIRlFRETEPq3bvw7v4d58vP5jetTuwa8jf812xRjsrByHhkJgIBw/DqaZ/GtgYPoL8uzZs2nRogWOjo74+vpSq1YtJkyYACRvyRYcHEzevHlTvNcwDKZOnUq/fv2IjY3l+vXrDBs2jClTpuDi4kLdunVZvHhx+gICCQkJzJ49m7Zt2wJQv3599u3bx969eylfvjxjxtx/H8FOnTqxZs2aO6736tWLsWPHpjujiIiI2J9L1y/hP8ef+h/Vx8XJhc0DNzO57WQ8XD2sjvZA7GpZxbBhEB19+7Xo6OTr7do9+LihoaEsWLDg1vejR4+mSpUqODk5ER8fT5s2be55v7e3N02bNmXcuHFcv36dDh068OijjwLg5+fHkCFDaJeegMCGDRtuZQJo0KDBrdeeffZZwsLC7jtGrVq1OHbs2B3XS5cuzYULFzh79izFihVLV04RERGxHyt+WkGPBT04F3WOwY0HM7zpcEtOvMtIdlWOT5xI2/XUiIuL4/fff6dMmTK3ruXPn59BgwbRo0cP9u/fn6pxhg8fTpUqVXBxceHfR2R7e3uzc+fOFO+pWbMmUVFRd1yfMGEC9erVu+3atm3bqFq1aorjzJ49m9atW6cq591UqVKFbdu20bJly3SNIyIiItnfX1f/oteCXizdtRSfkj580+sbqpSuYnWsDGFX5bhUqeSlFCldf1Dnz58nf/78d1xfvXo1Dz30EPv376dChQr3HSd37ty0bt2aPHnykCtXrlvXHR0dcXFxISoqCg+P23/8sGXLllTnPHPmDJUqVbrjenBwME5OTumemS5atCinT59O1xgiIiKSvZmmybyIefRd3JfrcdcJ9gtmYMOBODs5Wx0tw9hVOQ4OTl5j/O+lFe7uydcflJubG7GxsbddCw8P58qVK6xdu5bmzZvTsGFD3N3d7zuWg4MDDg53LvO+ceMGrq53/ggiLTPHKeUMCQkhPDyc9evXp3tPwdjYWNzc3NI1hoiIiGRfJy6coOv8rqzZt4bqj1ZnVsdZVCxe0epYGc6uyvH/JkeHDUteSlGqVHIxTs+kaYECBUhMTCQ2NhZXV1diYmLo378/K1aswNPTk2bNmhEcHExwcDB//vknHTp0YP369ake/8KFCxQpUgRn5zv/iystM8eVKlXiyJEjt75fs2YN48aNY9OmTbcV9wfJCHDo0CFeeeWVNN0jIiIi2V9SUhKfb/qcwcsGY2LyaZtP6VG7R4oTfvbA7n5X7drBsWOQlJT8azpXEwDJH27bunUrAKNGjcLPzw9PT08ARowYwaJFizh8+DBnzpy59YG41Nq4cSMvvvhiujM2btyYzZs33/o+KCiIqKgo6tevj4+PD926dQO4Z8Y2bdpQrVo1Dh48SIkSJZg1axYA8fHxHDlyBF9f33TnFBERkezj4NmDPD/+eYIWBFH90ersG7GPoBeC7LYYg53NHNtKUFAQEydOpF69eowePfq21zw8PDh69CgAkydPpmfPnncdZ8SIEXdcW7BgQaq2Wbuf0qVLU6hQIQ4fPky5cuVum0X+tx9++OGuGRcuXJji9fDwcFq1apXm4i8iIiLZU3xCPB9+9yEjvh6Bu4s7c96YQ4dqHbL80c8ZQW0nFZ588knq1KlDYmIijo6Od31fUFBQmsaNi4vDz88vVR/oS42xY8dy5swZypUrd9f3pDUjJO+h3L9///REExERkWxi94nd+If4s/vEblpWacnktpMpli/nbOVqmKZpdYZbfH19zX9vcwZw4MCBFHdhkKxBfz4iIiL2ITY+llHhoxi3ZhyF8xRmStsptKxqn1u4GoaxyzTNFNeLauZYREREJIfbdmQb/iH+HDx7kDeee4MJr0ygYO6CVseyhMqxiIiISA51LfYaQ1cMZfLGyZQqWIq1fdfSwKvB/W+0YyrHIiIiIjnQ2n1rCZwXyMlLJ+n1Qi+C/YLJ45rH6liWUzkWERERyUEuXr9Iv8X9CIkIoWKximx5awvPPfac1bGyDPvdpC6DHDt2DG9v7zTdM2fOnAc+annp0qV4eXnh4ODAfz+c+KAGDhxIxYoVefzxx2nevDmXL1/OkHFFREQke1m2axme73oyf8d8hr04jN3v7lYx/g+VYxtITzn29vZm+fLl1KpVK8Py1K9fn3379rF3717Kly+fIfsqi4iISPZx5vIZWn7eklZTW/F/+f+PyGGRvN/8fVydXa2OluVkq2UVfRf1Zc/JPRk6pk9JHya9Nume70lISKBjx47s3r2b8uXLM3fuXNzd3dm1axf9+vXj2rVrFC5cmDlz5rBt2zYiIyNp164dbm5uREREMH78eFauXElMTAzVq1dn2rRpd91E2xbbojVo8M/C+meffZawsLAMf4aIiIhkPaZpErI9hDeXvElMXAxjW4ylf4P+ODlmqwqYqTRznAoHDx4kMDCQvXv3kjdvXj777DPi4+Pp1asXYWFh7Nq1i86dOzNs2DBatWqFr68voaGh7NmzBzc3N4KCgti5cyf79u0jJiaG8PDwdOWJiorCx8cnxa/9+/ff897Zs2fTuHHjdD1fREREsr5j54/RcFJD3pjzBpUfrsze4XsZ1HiQivF9ZKt/O/eb4bWVkiVL8txzyetxXn/9dT755BMaNWrEvn37qF+/PgCJiYkUL148xfs3btzIBx98QHR0NBcvXsTLy4umTZs+cB4PDw/27En7DHpwcDBOTk60a9fugZ8tIiIiWVtiUiJTNk5h6IqhGBhMaTuFbs93w8FBc6Kpka3KsVX+uwTCMAxM08TLy4uIiIh73hsbG0uPHj2IjIykZMmSjBgxgtjY2HTliYqKombNmim+tmDBAjw9Pe+4HhISQnh4OOvXr88R56KLiIjkRAfOHCAgJIDtR7fTyLsR016fRqlCpayOla2oHKfCiRMniIiIoFq1aixcuJAaNWpQoUIFzp07d+t6fHw8hw4dwsvLCw8PD6KiogBuFeHChQtz7do1wsLCaNWqVbrypHXmeM2aNYwbN45Nmzbh7u6ermeLiIhI1hOfEM8Haz/gvfD3yJMrD3M7z+X1Z1/XhNgD0Px6KlSqVImQkBAef/xxLl68SPfu3XFxcSEsLIxBgwbxxBNP4OPjw/bt2wHo1KkT3bp1w8fHh1y5ctGlSxcqV66Mn58fTz311D2ftWLFCkqUKEFERARNmjShYcOG6c4fFBREVFQU9evXx8fHh27duqV7TBEREckafjr+E08FP8XbX76Nn48f+9/bT/tq7VWMH5BhmqbVGW7x9fU1/7u374EDB2yyg4NkDP35iIiIWCMmLoaRK0cy4dsJFPUoymftPsPvST+rY2ULhmHsMk3TN6XXtKxCREREJJvZcmgLAXMDOPTXIfxr+DPhlQnkd89vdSy7oHJskZ49e7Jt27bbrvXp04c33njDokQiIiKS1V2NucqQ5UP47PvPeKTwI6zrt466lepaHcuuqBxbZMqUKVZHEBERkWxk9S+r6Tq/K6cunaJvvb687/c+uXPltjqW3VE5FhEREcnCLly7wJuL32TeD/PwLO7J9kHbefbRZ62OZbdUjkVERESyINM0WRq5lKCFQVyKvsQ7L73DsBeHkcs5l9XR7JrKsYiIiEgWc/ryaXqG9uTLPV/iW9qXdf3W8XiJx62OlSOoHIuIiIhkEaZpMnvrbPov7c+NhBuMbzWevvX64uSoypZZ7O4QkNBfQikzqQwOIx0oM6kMob+EpnvMmJgYnn/+eRITE4mMjMTb25u4uDgAjh49StmyZbl69WqK95qmSY0aNVi9evWta0uWLKFRo0bExcVRq1YtEhIS0p0RYNKkScydOxeApUuX4uXlhYODA//dOzolJ0+epE6dOlSqVAkvLy8+/vjjW68NGDCADRs2ZEhGERERSdnv536n3sR6BMwNwKekD7+M+IUBDQeoGGcyuyrHob+EErgykONXjmNicvzKcQJXBqa7IM+ePZsWLVrg6OiIr68vtWrVYsKECUDylmzBwcHkzZs3xXsNw2Dq1Kn069eP2NhYrl+/zrBhw5gyZQouLi7UrVuXxYsXpysfQEJCArNnz6Zt27YAeHt7s3z5cmrVqpWq+52cnPjwww85cOAAP/zwA1OmTGH//v0A9OrVi7Fjx6Y7o4iIiNwpMSmRSesmUXlEZXYe28nU16eyof8GHiv6mNXRciS7+k+RYeuHER0ffdu16Phohq0fRrvK7R543NDQUBYsWHDr+9GjR1OlShWcnJyIj4+nTZs297zf29ubpk2bMm7cOK5fv06HDh149NFHAfDz82PIkCG0a/fg+QA2bNhwKxOQ5lPrihcvTvHixQHw8PCgUqVK/Pnnn3h6elK6dGkuXLjA2bNnKVasWLpyioiIyD9+/fNX/EP82fHHDppUbsLU16dSomAJq2PlaHZVjk9cOZGm66kRFxfH77//TpkyZW5dy58/P4MGDaJHjx63ZlfvZ/jw4VSpUgUXF5fbljl4e3uzc+fOFO+pWbMmUVFRd1yfMGEC9erVu+3atm3bqFq1aqqy3M+xY8fYvXs3zzzzzK1rVapUYdu2bbRs2TJDniEiIpKTxSXEMW7NOEaFjyKvW15CA0Jp83QbDMOwOlqOZ1fluFS+Uhy/cjzF6w/q/Pnz5M9/53GMq1ev5qGHHmL//v1UqFDhvuPkzp2b1q1bkydPHnLl+mcLFkdHR1xcXIiKisLDw+O2e7Zs2ZLqnGfOnEnzbHFKrl27RsuWLZk0adJtS0WKFi3K6dOn0z2+iIhITrfzj534h/jzy5+/0ObpNnz82scU8ShidSy5ya7WHAfXDcbd2f22a+7O7gTXDX7gMd3c3IiNjb3tWnh4OFeuXGHt2rUMHDiQ6Ojou9x9OwcHBxwc7vxXfuPGDVxdXe+4XrNmTXx8fO74WrduXapyplV8fDwtW7akXbt2tGjR4rbXYmNjcXNzS9f4IiIiOVn0jWgGLh3Is2Oe5eL1i3wd9DULuixQMc5i7Grm+H/rioetH8aJKycola8UwXWD07XeuECBAiQmJhIbG4urqysxMTH079+fFStW4OnpSbNmzQgODiY4OJg///yTDh06sH79+lSPf+HCBYoUKYKzs/Mdr6Vl5rhSpUocOXLkvu+7W0bTNPH396dSpUr069fvjvsOHTrEK6+8kuo8IiIi8o/vD35Pl7ldOPL3EQJrBfJByw/I557P6liSAruaOYbkgnys7zGShidxrO+xdBXj/2nQoAFbt24FYNSoUfj5+eHp6QnAiBEjWLRoEYcPH+bMmTO3PhCXWhs3buTFF19Md8bGjRuzefPmW9+vWLGCEiVKEBERQZMmTWjYsCHAXTNu27aNefPmsWHDhlsz1KtWrQKSZ5SPHDmCr69vunOKiIjkJFeir9BtXjfqTKiDaZps6L+Bae2nqRhnYTabOTYMowLw7z3KygLvmqY5yVbPtJWgoCAmTpxIvXr1GD169G2veXh4cPToUQAmT55Mz5497zrOiBEj7ri2YMECxowZk+6MpUuXplChQhw+fJhy5crRvHlzmjdvfsf7fvjhhxQz1qhRA9M0Uxw7PDycVq1apbn4i4iI5GTf7P2GrvO6cubKGfo36M97L7+Hey73+98olrJZ2zFN8yDgA2AYhiPwJ7DCVs+zpSeffJI6deqQmJiIo6PjXd8XFBSUpnHj4uLw8/NL1Qf6UmPs2LGcOXOGcuXK3fU9ac0IyXso9+/fPz3RREREcoxzUefou6gvC35cgPfD3izvsZynH3na6liSSsbdZgsz9CGG0QAYbprmc/d6n6+vr/nf09wOHDiQIbswiG3oz0dERCSZaZos+nERvRf15krMFd5u8jaDGw/GxcnF6mjyH4Zh7DJNM8X1opn1c/LXgIWZ9CwRERGRTHXq4im6h3YnfG84Tz/yNLM6zsL7YW+rY8kDsHk5NgzDBXgZGHKX1wOBQIBSpR58P2IRERGRzJaUlMTMrTMZGDaQ+MR4Jr46kd51e+PocPdlmJK1ZcbMcWPgJ9M0/0rpRdM0pwPTIXlZRSbkEREREUm3I38focvcLnx/8HvqVKjDjA4zeLToo1bHknTKjHLcBi2pEBERETuRmJTIpHWTeOerd3B2dGZGhxn41/DX0c92wqb7HBuG4Q7UB5bb8jm25ujoeNsJdWPHjgWgTJkyVK5cmcjISIYOHcqgQYNu3XP8+HHKli3L5cuXbxtrzpw5KW7pFhkZibe3N3FxcQAcPXqUsmXLcvXqVZv8ntq1a0fBggUJCwuzyfgiIiL2aN+f+6g2phoDlg6gfqX67B+5n4CaASrGdsSmM8emaUYDhWz5jMzg5ubGnj17Unxt48aNFC5cGC8vL5588kk6depEpUqV6NOnD6NGjSJ//vypeoavry+1atViwoQJDB06lJ49exIcHEzevHkz8rdyS2hoKJ06dbLJ2CIiIvbmRvwNxqwew+hVo8nvnp9FgYt41fdVlWI7lL1OdejbF+5SUh+Yjw9MSv+5JG5ubkycOJEePXrw1ltvERUVRbt2aTudb/To0VSpUgUnJyfi4+Np06ZNujIlJCRQrVo1xo8fT+3atRkyZAgODg4EBwena1wREZGcZMfvO/AP8efX07/y+rOv89GrH1HYo7DVscRGslc5tkhMTAw+Pj63vh8yZAitW7e+430vvvgis2bNokOHDreOm06L/PnzM2jQIHr06MH+/ftTfE9UVBQ1a9ZM8bUFCxbcOtYawMnJiTlz5tCqVSs++eQT1qxZw44dO9KcS0REJCe6fuM673z5DpPWT+Lh/A8T3iucJo83sTqW2Fj2KscZMMP7IO61rOK/evbsSUxMzAOferd69Woeeugh9u/fn+IYHh4eqc4C4OXlRfv27WnatCkRERG4uGgjchERkfvZcGADXeZ14fdzv9O9dnfGthhLXjfbLHWUrCV7leNswMHBAQeHB/ucY3h4OFeuXGHt2rU0b96chg0b4u5++xnsaZk5/p9ffvmF/Pnz89dfKe6mJyIiIjddjr7MwLCBzNwyk8eKPsb3A77n+QrPWx1LMpHKcRYRExND//79WbFiBZ6enjRr1ozg4OA71gendeZ4+fLlXLhwgc2bN/PSSy/x448/pvpDgiIiIjnJ13u+pntod85eOctbDd9ixMsjcHNxszqWZDKbbuVmL/635vh/X4MHD87wZ4waNQo/P79bM78jRoxg0aJFHD58+IHHPH/+PIMHD2bWrFmUL1+eoKAg+vTpk1GRRURE7MLfV//mtemv0WxKMwrnKcyOoTsY12qcinEOpZnjVEhMTEz1e2vXrk3t2rXT/IzRo0ff9r2HhwdHjx5N8zj/VrhwYQ4dOnTr+969e6drPBEREXtimiahO0Lps6gP125cY1SzUQxqNAhnJ2ero4mFNHOcDkWKFKFu3bpERkZaHSXN2rVrx6ZNm3B1dbU6ioiISKY7efEkL336Eu1ntaf8Q+XZ/c5u3n7pbRVj0cxxeuzcuTPN9/j4+FCmTJmMD5NGoaGhVkcQERHJdElJSUzbPI1BywYlHwPdehJBLwTh6OBodTTJIrJFOTZN025OoPn3fsnZnWmaVkcQERFJtcN/HSZgbgCbD22mXqV6TG8/nUeKPGJ1LMlisnw5dnV15cKFCxQqVMhuCrI9ME2TCxcuaFmGiIhkeQmJCUz8biLDvx6Oq7MrszvNplP1TuoVkqIsX45LlCjBqVOnOHfunNVR5D9cXV0pUaKE1TFERETu6ueTP+Mf4s+u47to/mRzprSdQvH8xa2OJVlYli/Hzs7OPPKIfuQhIiIiqXcj/gbvf/M+Y9eMpaB7QZZ2W0rLKi01Wyz3leXLsYiIiEhaRByNwD/EnwNnDtChWgcmvjqRQnkKWR1LsgmVYxEREbEL12Kv8faXb/PJhk8oWaAkq/usppF3I6tjSTajciwiIiLZ3nf7vyNwbiDHLhwjqE4Qo1uMxsPVw+pYkg2pHIuIiEi2den6JQYsHcDsbbOpUKwCW97aQo1yNayOJdmYyrGIiIhkSyt+WkGPBT04F3WOIY2H8G7Td3F11hajkj4qxyIiIpKtnL1yll4LexG2Kwyfkj580+sbqpSuYnUssRMqxyIiIpItmKbJvIh59F3cl+i4aEY3H82ABgNwdnK2OprYEZVjERERyfKOXzhO13ldWfvrWqo/Wp1ZHWdRsXhFq2OJHVI5FhERkSwrKSmJzzd9zuBlgzEx+bTNp/So3QMHBwero4mdUjkWERGRLOng2YMEhASw9chWGno1ZFr7aZQuVNrqWGLnVI5FREQkS4lPiGfCtxMYuXIk7i7uzHljDh2qddDRz5IpVI5FREQky9h9Yjf+If7sPrGbVlVb8WmbTymWr5jVsSQHUTkWERERy8XGx/Leyvf4YO0HFM5TmGXdl9GiSgurY0kOpHIsIiIiltp2ZBv+If4cPHuQN557gw9f+ZACuQtYHUtyKJVjERERsURUbBRDlw9lyvdTKF2wNGv7rqWBVwOrY0kOp3IsIiIimW7tvrUEzgvk5KWT9HqhF8F+weRxzYDPeI8AACAASURBVGN1LBGVYxEREck8F69f5M3FbzI3Yi4Vi1Vk61tbqf5YdatjidyiciwiIiKZImxXGD1De3Ix+iLDXhzG2y+9jauzq9WxRG6jciwiIiI2debyGYIWBrH8p+VUKVWFtX3X4lPKx+pYIilSORYRERGbME2TOdvn0G9JP2LiYhjbYiz9G/THyVH1Q7Iu/a9TREREMtyx88cInBfId/u/o2a5mszsMJPyxcpbHUvkvlSORUREJMMkJiUyZeMUhq4YioHBZ+0+o2utrjg4OFgdTSRVVI5FREQkQxw4cwD/EH8ijkbQ2LsxU1+fSqlCpayOJZImKsciIiKSLvEJ8Xyw9gPeC3+PPLnyMM9/Hu2eaYdhGFZHE0kzlWMRERF5YLuO76LznM7sPbWXV31f5dM2n1I0b1GrY4k8MJVjERERSbOYuBhGrhzJhG8nUNSjKCt6rMDvST+rY4mkm8qxiIiIpMnmQ5sJCAng8N+HCagZwPhW48nvnt/qWCIZQuVYREREUuVqzFUGLx/M599/ziOFH2Fdv3XUrVTX6lgiGcqm5dgwjPzATMAbMIHOpmlG2PKZIiIikvFW/bKKbvO7cerSKd6s9yaj/EaRO1duq2OJZDhbzxx/DKwxTbOVYRgugLuNnyciIiIZ6HzUed5c8ibzf5iPZ3FPtg/azrOPPmt1LBGbsVk5NgwjL1AL6ARgmmYcEGer54mIiEjGMU2TpZFLCVoYxKXoS7z70rsMfXEouZxzWR1NxKZsOXNcFjgHfGEYxhPALqCPaZrXbfhMERERSafTl0/TI7QHX+35Ct/Svqzrt47HSzxudSyRTGHLsxydgCrA56ZpPglcBwb/902GYQQahhFpGEbkuXPnbBhHRERE7sU0TWZumYnnu56s/XUtE16ZQMSQCBVjyVFsOXN8CjhlmuaOm9+HkUI5Nk1zOjAdwNfX17RhHhEREbmL38/9Tpe5Xdjw2waeL/88MzvO5LGij1kdSyTT2awcm6Z51jCMk4ZhVDBN8yBQF9hvq+eJiIhI2iUmJfLJ+k8Y9uUwnBycmNZ+GgE1AnBwsOUPl0WyLlvvVtELCL25U8XvwBs2fp6IiIik0q9//op/iD87/thBk8pNmPr6VEoULGF1LBFL2bQcm6a5B/C15TNEREQkbeIS4hi7eizvf/M++dzysSBgAa89/RqGYVgdTcRyOiFPREQkB9n5x078Q/z55c9faPt0Wya9NokiHkWsjiWSZagci4iI5ADRN6IZ/vVwJn43keL5ivN10Nc0faKp1bFEshyVYxERETv3/cHvCQgJ4Oi5o3St1ZVxLceRzz2f1bFEsiSVYxERETt1JfoKby17i+mbp/NokUfZ0H8DdSrWsTqWSJamciwiImKHwn8Op9v8bpy5coYBDQYw8uWRuOdytzqWSJanciwiImJHzkWdo8+iPiz8cSHeD3uzvMdynn7kaatjiWQbKsciIiJ2wDRNFv24iN6LenMl5gojXx7J4MaDcXFysTqaSLaiciwiIpLNnbp4iu6h3QnfG84zjzzDrI6z8HrYy+pYItmSyrGIiEg2lZSUxIwtMxgYNpCEpAQmvjqR3nV74+jgaHU0kWxL5VhERCQbOvL3EbrM7cL3B7/nhYovMKPDDMoWKWt1LJFsT+VYREQkG0lITGDSukm889U7uDi5MKPDDPxr+OvoZ5EMonIsIiKSTfxy6hf8Q/zZeWwnLz/xMp+1+4yHCzxsdSwRu6JyLCIiksXdiL/B6FWjGb16NAXcC7A4cDGv+L6i2WIRG1A5FhERycJ2/L4D/xB/fj39K68/+zqTWk+iUJ5CVscSsVsqxyIiIlnQ9RvXeefLd5i0fhIP53+Yb3p/w4uVX7Q6lojdUzkWERHJYtYfWE+XuV344/wfdK/dnbEtxpLXLa/VsURyBJVjERGRLOJy9GUGhg1k5paZlCtajk0DN1GrfC2rY4nkKCrHIiIiWcBXe76i+/zu/B31N4MaDWJ40+G4ubhZHUskx1E5FhERsdDfV/+m96LeLN65mCdKPMHKXiupWrqq1bFEciyVYxEREQuYpknojlD6LOrDtRvXeN/vfd5q+BbOTs5WRxPJ0VSORUREMtmJCyfoNr8bq/etptqj1ZjVcRaVileyOpaIoHIsIiKSaZKSkpi2eRpvhb1FkpnEx699TM86PXF0cLQ6mojcpHIsIiKSCQ6dPUTA3AC2HN5CvUr1mN5+Oo8UecTqWCLyHyrHIiIiNpSQmMDE7yYy/OvhuDq7MrvTbDpV76Sjn0WyKJVjERERG/n55M90ntOZn078RPMnmzOl7RSK5y9udSwRuQeVYxERkQwWGx/L++HvM27tOArlLkRYtzBaVm1pdSwRSQWVYxERkQy0/ch2/EP8+e3sb3Ss1pGJrSdSMHdBq2OJSCqpHIuIiGSAa7HXGPblMD7d8CklC5RkTZ81NPRuaHUsEUkjlWMREZF0+m7/dwTODeTYhWME1QlidIvReLh6WB1LRB6AyrGIiMgDunT9Ev2X9ueLbV9QoVgFtry1hRrlalgdS0TSQeVYRETkASz/aTk9F/TkXNQ5hjQewrtN38XV2dXqWCKSTirHIiIiaXD2ylmCFgSx7Kdl+JT0YVXvVTxZ6kmrY4lIBlE5FhERSQXTNJkbMZc3F79JdFw0o5uPZkCDATg7OVsdTUQykMqxiIjIfRy/cJyu87qy9te1PPfYc8zsMJOKxStaHUtEbEDlWERE5C6SkpL47PvPGLx8MACftvmUHrV74ODgYHEyEbEVlWMREZEUHDx7EP8Qf7Yd2UZDr4ZMaz+N0oVKWx1LRGxM5VhERORf4hPimfDtBEauHIm7izshb4TQvlp7DMOwOpqIZAKVYxERkZt2n9hN5zmd2XNyD62qtmJy28k8lPchq2OJSCZSORYRkRwvNj6WkStHMn7teIp4FGFZ92W0qNLC6lgiYgGVYxERydG2Ht6Kf4g/h/46xBvPvcGHr3xIgdwFrI4lIhaxaTk2DOMYEAUkAgmmafra8nkiIiKpFRUbxZDlQ5iycQplCpXh2ze/pb5nfatjiYjFMmPmuI5pmucz4TkiIiKpsnbfWgLnBXLy0kn61O3D+37vk8c1j9WxRCQL0LIKERHJMS5cu0C/Jf2YGzGXSsUrsW3QNqo9Ws3qWCKShdi6HJvAt4ZhmMA00zSn2/h5IiIidzBNk2U/LaNnaE8uRl/k7SZv83aTt8nlnMvqaCKSxdi6HD9nmuZpwzCKAt8ZhvGbaZqb//0GwzACgUCAUqVK2TiOiIjkNGcun6Hngp6s2L2CqqWr8u2b3/JEySesjiUiWZRNz780TfP0zV//BlYAT6fwnummafqapulbpEgRW8YREZEcxDRNvtj2BZ7DPVm9bzXjWo7jhyE/qBiLyD3ZbObYMIzcgINpmlE3/7kB8J6tniciIvI/f5z7g8B5gaw7sI5a5Wsxo/0Myhcrb3UsEckGbLms4iFgxc3jNp2ABaZprrHh80REJIdLTEpkysYpDFk+BEcHRz5v9zmBtQJxcLDpD0pFxI7YrBybpvk7oJ9diYhIpth/ej8BcwOIOBpBY+/GTGs/jZIFS1odS0SyGW3lJiIi2Vp8Qjzj1oxj1Dej8HD1YL7/fNo+05abP7kUEUkTlWMREcm2dh3fRec5ndl7ai+tn2rNJ699QtG8Ra2OJSLZmMqxiIhkOzFxMYz4egQTvp3AQ3kf4sueX9LMp5nVsUTEDqgci4hItrL50GYCQgI4/PdhutTswgetPiC/e36rY4mInVA5FhGRbOFqzFUGLx/M599/TtkiZVnfbz0vVHrB6lgiYmdUjkVEJMtb9csqus7ryunLp+lXvx/vNXuP3LlyWx1LROyQyrGIiGRZ56PO03dxX0J3hOJZ3JOwwWE8U/YZq2OJiB1TORYRkSzHNE2WRC6h18JeXIq+xPCmwxnSeAi5nHNZHU1E7JzKsYiIZCmnL5+m+/zufP3z1/iW9mV9v/VULlHZ6lgikkOoHIuISJZgmiazts5iwNIBxCXGMeGVCfSp2wcnR/1flYhkHv2NIyIiljv691EC5wWy4bcN1K5QmxkdZvBY0cesjiUiOZDKsYiIWCYxKZGP133M21+9jbOjM9PaTyOgRgAODg5WRxORHErlWERELLHvz334h/jz4x8/8tLjL/F5u88pUbCE1bFEJIdTORYRkUwVlxDHmFVjCF4VTD63fCwIWMBrT7+GYRhWRxMRUTkWEZHMs/OPnXQO6cy+P/fR9um2THptEkU8ilgdS0TkFpVjERGxuegb0bz79bt89N1HFM9XnJVBK3npiZesjiUicgeVYxERsamNv22ky9wuHD13lK61ujKu5TjyueezOpaISIpUjkVExCauRF/hrWVvMX3zdB4t8igbB2ykdoXaVscSEbknlWMREclwK39eSbf53Th75SwDGgxg5Msjcc/lbnUsEZH7UjkWEZEMcy7qHH0W9WHhjwup/HBlvuzxJU898pTVsUREUk3lWERE0s00TRb+uJDei3pzNeYqI18eyeDGg3FxcrE6mohImqgci4hIupy6eIruod0J3xvOM488w6yOs/B62MvqWCIiD0TlWEREHkhSUhIztsxgYNhAEpMS+aj1R/R6oReODo5WRxMReWAqxyIikmaH/zpMl7ld2HRoE3Ur1WV6++mULVLW6lgiIummciwiIqmWkJjApHWTeOerd8jllIuZHWbSuUZnHf0sInZD5VhERFJl76m9+M/xJ/J4JM18mvFZu8/4v/z/Z3UsEZEMpXIsIiL3dCP+BqNXjWb06tEUcC/A4sDFvOL7imaLRcQuqRyLiMhd/XD0B/xD/Nl/Zj/tn23PR60/olCeQlbHEhGxGZVjERG5w/Ub13nny3eYtH4SJQqUYFXvVTSu3NjqWCIiNqdyLCIit1l/YD1d5nbhj/N/0KN2D8a0GENet7xWxxIRyRQqxyIiAsDl6MsMWDqAWVtnUa5oOTYN3ESt8rWsjiUikqlUjkVEhK/2fEX3+d35O+pvBjUaxPCmw3FzcbM6lohIplM5FhHJwf66+he9F/ZmSeQSnijxBCt7raRq6apWxxIRsYzKsYhIDmSaJvN/mE/fxX25duMawX7BDGw4EGcnZ6ujiYhYSuVYRCSHOXHhBN3md2P1vtVUe7QaszrOolLxSlbHEhHJElSORURyiKSkJKZumsqgZYMwMfnktU/oUacHjg6OVkcTEckyVI5FRHKAQ2cPETA3gC2Ht1Dfsz7T20+nTOEyVscSEclyVI5FROxYQmICH377IcO/Tt594otOX9Cxekcd/SwichcqxyIidurnkz/TeU5nfjrxEy2qtGBym8kUz1/c6lgiIlmayrGIiJ2JjY/l/fD3Gbd2HIVyFyKsWxgtq7a0OpaISLagciwiYke2H9mOf4g/v539jY7VOjKx9UQK5i5odSwRkWzD5uXYMAxHIBL40zTNl2z9PBGRnOha7DWGrhjK5I2TKVmgJGv6rKGhd0OrY4mIZDt3LceGYawCepimeSydz+gDHADypnMcERFJwbe/fkvgvEBOXDxBz9o9Gd1iNB6uHlbHEhHJlhzu8doc4FvDMIYZhvFARyYZhlECaALMfJD7RUTk7i5dv8QbX7xBw0kNcXV2ZfPAzXza9lMVYxGRdLjrzLFpmksMw/gGeBeINAxjHpD0r9cnpmL8ScBbwF3/pjYMIxAIBChVqlQqY4uI5GzLf1pOzwU9ORd1jqEvDuWdl97B1dnV6lgiItne/dYcxwPXgVwkF9yke7/9H4ZhvAT8bZrmLsMwat/tfaZpTgemA/j6+pqpHV9EJCc6e+UsQQuCWPbTMp4s9SSre6/Gp5SP1bFEROzGvdYcNwImAl8DVUzTjE7j2M8BLxuG8SLgCuQ1DGO+aZqvP3BaEZEcyjRNQraH0G9JP6LjohnTYgz96/fH2emBVr2JiMhd3GvmeBjwimmavz7IwKZpDgGGANycOR6gYiwiknbHzh+j67yufLv/W2o8VoOZHWdSoVgFq2OJiNile605rpmZQURE5HZJSUlM2TiFISuGYGAwue1kuj/fHQeHe32WWkRE0iNTDgExTfN74PvMeJaIiD347cxvBMwNYNuRbTT0asi09tMoXai01bFEROyeTsgTEclC4hPiGb92PCPDR5InVx5C3gihfbX2GIZhdTQRkRxB5VhEJIv46fhP+If4s+fkHl6p+gqftv2Uh/I+ZHUsEZEcReVYRMRiMXExvBf+HuPXjqeIRxGWd19O8yrNrY4lIpIjqRyLiFho6+Gt+If4c+ivQ3R+rjMTXplAgdwFrI4lIpJjqRyLiFggKjaKIcuHMGXjFMoUKsN3b35HPc96VscSEcnxVI5FRDLZmn1r6DqvKycvnaRvvb6MajaKPK55rI4lIiKoHIuIZJoL1y7Qb0k/5kbMpVLxSmwbtI1qj1azOpaIiPyLyrGIiI2ZpknYrjCCFgRxMfoi77z0DsNeHEYu51xWRxMRkf9QORYRsaEzl8/QI7QHX+75kqqlq/Ltm9/yRMknrI4lIiJ3oXIsImIDpmnyxbYv6LekHzcSbvBBqw94s96bODnqr10RkaxMf0uLiGSwP879QeC8QNYdWEet8rWY0X4G5YuVtzqWiIikgsqxiEgGSUxKZPKGyQxdMRRHB0c+b/c5gbUCcXBwsDqaiIikksqxiEgG2H96P/4h/vzw+w+8WPlFpr4+lZIFS1odS0RE0kjlWEQkHeIS4hi3Zhzvf/M+Hq4ezPefT9tn2mIYhtXRRETkAagci4g8oMhjkfiH+LP31F5ee+o1Pn7tY4rmLWp1LBERSQeVYxGRNIqJi2H418P58NsPKZavGF/1/IqXfV62OpaIiGQAlWMRkTTYdHATAXMDOPL3EbrU7MIHrT4gv3t+q2OJiEgGUTkWEUmFqzFXGbRsEFM3TaVskbKs77eeFyq9YHUsERHJYCrHIiL38c3eb+g2vxunL5+mX/1+jGo2Cvdc7lbHEhERG1A5FhG5i/NR5+m7uC+hO0Lx+j8vwrqF8UzZZ6yOJSIiNqRyLCLyH6ZpsnjnYnot7MWVmCsMbzqcoS8OxcXJxepoIiJiYyrHIiL/8uelP+kR2oOvf/6ap8o8xayOs6hcorLVsUREJJOoHIuIkDxbPHPLTAaEDSA+MZ4Jr0ygb72+ODo4Wh1NREQykcqxiOR4R/8+Spe5Xdh4cCO1K9RmRocZPFb0MatjiYiIBVSORSTHSkxK5ON1H/P2V2/j7OjM9PbTCagZoKOfRURyMJVjEcmR9v25D/8Qf37840eaPtGUz9t9zsMFHrY6loiIWEzlWERylLiEOMasGkPwqmDyueVjYZeFtH6qtWaLRUQEUDkWkRzkxz9+xD/En31/7qPt0235+LWPKexR2OpYIiKShagci4jdi74Rzbtfv8tH331E8XzFWRm0kpeeeMnqWCIikgWpHIuIXdv420YC5gbw+7nf6fZ8N8a2GEs+93xWxxIRkSxK5VhE7NKV6Cu8tewtpm+ezmNFH+P7Ad/zfIXnrY4lIiJZnMqxiNidlT+vpNv8bpy9cpaBDQcyoukI3HO5Wx1LRESyAZVjEbEb56LO0XthbxbtXETlhyvzVc+v8C3ja3UsERHJRlSORSTbM02ThT8upPei3lyNucp7zd5jUKNBuDi5WB1NRESyGZVjEcnWTl48Sff53fnml294tuyzzOwwE6+HvayOJSIi2ZTKsYhkS0lJSczYMoOBYQNJTEpkUutJBL0QhKODo9XRREQkG1M5FpFs5/Bfh+kytwubDm2ibqW6TG8/nbJFylodS0RE7IDKsYhkGwmJCXy07iPe/epdcjnlYlbHWbzx3Bs6+llERDKMyrGIZAt7T+3Ff44/kccjaebTjM/afcb/5f8/q2OJiIidUTkWkSztRvwNglcFM2b1GAq6F2RJ1yW0qtpKs8UiImITNivHhmG4ApuBXDefE2aa5nBbPU9E7M8PR3/AP8Sf/Wf20/7Z9nzU+iMK5SlkdSwREbFjtpw5vgG8YJrmNcMwnIGthmGsNk3zBxs+U0TswPUb13n7y7f5eP3HlChQglW9V9G4cmOrY4mISA5gs3JsmqYJXLv5rfPNL9NWz3tgUVGQJw/oR7QiWcK6/evoMrcLxy4co2ednoxpMQYPVw+rY4mISA7hYMvBDcNwNAxjD/A38J1pmjts+bwH0r49lC8PgwfDzp1gZr3+LpITXI6+jP8cf+p/VB9nR2c2D9zM5LaTVYxFRCRT2bQcm6aZaJqmD1ACeNowDO//vscwjEDDMCINw4g8d+6cLeOkzM8PHnkEJkyAp5+GMmXgzTdh61ZISsr8PCI50Je7v8TzXU9CIkIY3HgwPw//mZrla1odS0REciDDzKSZUsMwhgPXTdOccLf3+Pr6mpGRkZmS5w4XL8LXX0NYGHz3HcTFQbFi0KIFtGwJtWqBkzb3EMlIf139i14LerF011KeKPEEszrNomrpqlbHEhERO2cYxi7TNH1Tes1mM8eGYRQxDCP/zX92A+oBv9nqeelWsCB06gTh4XDuHISGQvXq8MUXULcuFC8OAQGwenVycRaRB2aaJvMi5uH5ridf/fwVwX7B7By2U8VYREQsZ7OZY8MwHgdCAEeSS/gS0zTfu9c9ls4c383167BmDSxbllyco6IgXz5o2hRatYIGDcDNzeqUItnGiQsn6Dq/K2v2raH6o9WZ1XEWFYtXtDqWiIjkIPeaOc60ZRWpkSXL8b/FxiYvuVi2LHkJxqVLkDs3NGmSvPTixReTd74QkTskJSXx+abPGbxsMCYmY5qPoWednjg42PSjDyIiIndQObaF+HjYuDG5KK9YkbwUw9UVGjZMLspNm0L+/FanFMkSDp49SEBIAFuPbKW+Z32mt59OmcJlrI4lIiI5lMqxrSUmwpYtyUV5+XI4fRqcnaFeveSi3KwZFC5sdUqRTJeQmMCEbycw4usRuLm48dGrH9Gxekcd/SwiIpZSOc5MSUmwY0dyUV62DI4dA0dHeP755KLcvHnyh/tE7NyeE3vwD/HnpxM/0aJKC6a0nUKxfMWsjiUiIqJybBnThN27k0tyWBgcOpR8El/16skf5mvRAkqVsjqlSIaKjY9lVPgoxq0ZR+E8hZnSdgotq7a0OpaIiMgtKsdZgWnCr7/+M6P8yy/J1596KnlGuWVLeOwxazOKpNO2I9sICAngt7O/0al6Jz589UMK5i5odSwREZHbqBxnRYcP/1OU//d7fvzx5JLcqhV4elqbTyQNrsVeY+iKoUzeOJlSBUsxvf10Gng1sDqWiIhIilSOs7pjx5I/yLdsGWzfnnytYsV/ZpR9fJKXY4hkQd/++i2B8wI5cfEEQXWCGN18NHlctaWhiIhkXSrH2cnp08lbwy1bBps2JX/Ar2zZf46xfvpp0L6wkgVcvH6R/kv6M2f7HCoUq8CsjrN47rHnrI4lIiJyXyrH2dW5c/DVV8lFed06SEiAhx/+Z0b5ueeSd8IQyWTLdi2j54KenL92nkGNBvHOS+/g6uxqdSwREZFUUTm2B5cuwcqVyUV57Vq4cQOKFk3eGq5lS6hdO3lvZREbOnvlLEELglj20zKeLPUkszvOxqeUj9WxRERE0kTl2N5ERcGqVclF+ZtvIDoaChaEl19OLsr160OuXFanFDtimiYh20Pot6Qf0XHRjHx5JP0b9MfJ0cnqaCIiImmmcmzPYmKSZ5LDwpJnlq9ehbx54aWXkotyo0bg7m51SsnGjp0/RuC8QL7b/x01HqvBzI4zqVCsgtWxREREHpjKcU5x4wasX588o/zVV3DhQnIxbtw4uSg3aZJcnEVSISkpiSkbpzBkxRAMDMa1HEe357vhoA+EiohINqdynBMlJCTvdrFsWfI2cX/9lbzUokGD5KL88stQoIDVKSWLOnDmAAEhAWw/up1G3o2Y+vpUShcqbXUsERGRDKFynNMlJibvn/y/onzyJDg5wQsvJBdlP7/kD/dJjhefEM/4teMZGT6SPLnyMKn1JF5/9nUM7bMtIiJ2ROVY/mGasHPnP6fzHT2avG9yzZrJJ/M1b568XZzkOD8d/wn/EH/2nNzDq76v8kmbT3go70NWxxIREclwKseSMtOEvXuTP8y3bBkcOJB8vVq1f/ZSLlPG0ohiezFxMbwX/h7j146niEcRPm/3OX5P+lkdS0RExGZUjiV1Dhz4Z0Z5z57ka1Wq/FOUK2iHAnuz5dAWAuYGcOivQ/jX8Gd8q/EUyK216CIiYt9UjiXtjh5NXp+8bBns2JF8zdv7n6Ls7Q1ah5ptRcVGMXjZYD77/jPKFCrDjA4zqOdZz+pYIiIimULlWNLn5Ml/ivLWrcnLMcqV+6coV62qopyNrP5lNV3nd+XUpVP0qduH9/3eJ3eu3FbHEhERyTQqx5Jxzp6FL79MLsobNybvhFG6NLRokfyBvmefTf6An2Q5F65d4M3FbzLvh3l4FvdkZseZVHu0mtWxRERE/r+9+w6PskrbAH6fJJPeE3pISEhfFhFdVAQBkd4zrKuLwq4BXBULu/bysbrqZ8WKBYFrgbWtDL2KomJB5QNFIAVCIAQIENJ7JjPn++MMM4EZMMEk75T7d11cMZMT55nX1+TmzDnP6XAMx9Q+SkqAtWvVhr6tWwGjEejeXXW80OtVBwwfHi+sNSklVuxagTkfzEFpbSkeGfMIHhv7GPx0PGKciIg8E8Mxtb+KCmD9ejWjvGkTUF8PREerHsp6veqp7OurdZUe50T5Cdz1/l1Y/fNqXBF3BZb8ZQn6xvTVuiwiIiJNMRxTx6qpUQHZYFCBuboaCA9Xp/Lp9eqUPn9/rat0a1JKLPlmCf7xyT/Q0NSApyY9hbk3zIWPN2fyiYiIGI5JO/X1wKefqqC8di1QXg4EBwPjxqmgPHYsEMTNYG0pvzgfs5fPxufZn+O65OuwaPoilMC2MQAAIABJREFUJHVJ0rosIiIip8FwTM6hsVFt4jMY1Ka+4mI1gzx6tNrMN348EBamdZUuy2Q24Y1tb+CxVY/B28sbL+hfwOzrZsOLGySJiIjOwXBMzqepCfj6axWUV64EiooAnQ4YMULNKE+aBERFaV2ly8g6kYXMpZn4Pv97jPv9OLx9y9voGdlT67KIiIicEsMxOTezGfj+e9vpfAUFgLc3MHSoCspTpgBdu2pdpVNqbGrE85ufx9MbnkaIfwhev+l13DzgZgj2nSYiIroghmNyHVICu3erkLxiBXDwoDpgZNAgFZQzMoCenBEFgJ2HdyJzaSb2Ht+Lm/5wE16/+XV0CumkdVlEREROj+GYXJOUwL59thnlffvU4wMG2E7n691b2xo1UNtQi3+u+yde/vRldA3rirenvY2J/SZqXRYREZHLYDgm93DggC0o79qlHrvsMrWZT68H0tK0ra8DfJX7FWYum4m803mYNXgWXpz6IsICuYmRiIioNRiOyf0cPqw28hkMwI4d6rG0NNuM8mWXqeUYbqKyrhIPGR7CO1+9g4ROCXjv1vdwfdr1WpdFRETkkhiOyb0dPw6sWqWC8vbtaoNfQoItKA8Y4NJBecMvG/C3//wNJ8pPYO6IuXhq4lMI9AvUuiwiIiKXxXBMnuP0aWDNGhWUP/9ctYzr2VNt5NPrgYEDVScMF1BcVYz7ProPH/z4Afr06IPFMxZjQPwArcsiIiJyeQzH5JnKytSpfAaDOqWvoQHo0kW1htPrVas4H+c7TllKiY93foy7P7wbFXUVeGzsY3hk7CPw9fHVujQiIiK3wHBMVFUFbNiggvLGjUBtLRAZqQ4bmToVGD4c8PPTukocLzuOO96/A+v2rMOA+AFYPGMx+vToo3VZREREboXhmKi52lpgyxbVR3ndOhWcQ0OBCRPUjPLo0UBAQIeWJKXEoq8X4f4V98NoMuLpSU/j3hvuhbeXaywBISIiciUMx0QX0tAAfPaZmlFeswYoLQUCA4GxY1VQHjcOCAlp1xIOnT6EWctm4YvcLzAsZRjem/4eenf2vP7NREREHeVi4dj5FlwSdSQ/PxWAx40DjEbgq69UUF61Ss0s+/kBo0apoDxhAhAR0WZPbTKb8Npnr+HxNY9D563DwlsXYubgmTz6mYiISEOcOSZyxGQCvv1WBeWVK4Fjx9TmveHDVVCePBnodOlHNe87vg+ZSzPx4+EfMeGyCXh72tvoEdGjDV8AERERXYgmyyqEED0BLAPQFYAZwEIp5WsX+x6GY3JKZjOwc6ftdL78fMDLCxgyRAXlKVOA7t1b9K9qbGrEsxufxbMbn0V4YDjeuPkN3HjljZwtJiIi6kBaheNuALpJKXcLIUIA7AIwWUqZdaHvYTgmpyclsGePCskrVgA5OerxgQNth47ExQEA3t/7Ph77/DEcrTiK2LBY3Pb72/Df7f/F/hP7Me2qaXj1T68iOiRawxdDRETkmZxiQ54QYg2AN6WUWy80huGYXE5Wlm1Gec8e9dgVV+CnQb0xw2st9obV28aagYjGCCy/eTnG9R2nTb1ERESkfTgWQvQCsB1AHyll5YXGMRyTS8vLU+uTDQbgxx8BAL90BgzpgCEN2N8Z6BnWE0fnHtW4UCIiIm1UVzciJ6cUWVklyM4uwfDhcbjhhrgOr0PTbhVCiGAABgD3OQrGQojZAGYDQGxsbHuXQ9R+EhOBBx8EHnwQcXMFpmQD+mxg3pfAk18CuVHAyrRC4LpdQP/+ANcZExGRmyorq7cG4OYfjx6tso7R6bwQHu6nSTi+mHadORZC6ACsB7BFSjn/18Zz5pjcwdqf1yJjVQZMXiYAQNcqYHIOoM8ChhYAPmYAvXrZ1ihfdZXa4EdERORCpJQ4darWLgBnZ5fi5Mka67iAAB+kpkYiPT0KaWlRlo+R6N07HDqdNoddabUhTwBYCqBUSnlfS76H4Zhc2enK07jno3vw8c6PEdMjBmd8zqDeZFtzHKgLxNJBL2Nqnq/azPfZZ6q3co8equOFXg8MHgx481Q8IiJyHlJKFBZW2QXgrKwSlJXZfs+Fhvpaw+/ZAJyeHoW4uDB4eTnXu6VaheNBAL4GsBeqlRsAPCql3Hih72E4JlckpcQHP3yAez++F1X1VXhi3BN4cPSD+CT7k3O6VTwz/BlM+/002zeWlwPr16s1yps3A/X1qnfy5MnA1KnAsGGATqfdCyMiIo9iMpmRn1/hcCa4psZoHdepU0CzGWBbEO7ePdhlWpNqviGvpRiOydUUlhbijv/cgQ17N+DqhKuxeMZipHdPb/2/qLoa2LRJBeUNG9TnERHAxIlqRnnECMDfv+1fABEReZzGRhMOHixDVta5IfjAgTI0NJis43r0CLYLwGlpUejUKVDD6tsGwzFRGzObzVi4fSEeNDwIk9mEZ6c8iznXz4G3VxssiairAz79VAXltWuBigogJEQdca3XA2PGAEFBv/15iIjIrdXWGs/pDHF2Fjgvrwwmk8p/QgDx8WHnBOD09GikpkYiLMxP41fQfhiOidrQwVMHMXPZTGw/sB03pN2AhbcuRHyn+PZ5ssZGYNs2FZRXrwbOnAECAlRA1uuB8eOB0ND2eW4iInIJ5eX1yM4utVsOUVBQibMxz8fHC4mJ4efMBKenRyE5OQKBgZ63hI/hmKgNNJmaMH/rfMxbOw9+Pn6Yf+N8/PXav3bc+qqmJmD7dhWUV60CiooAX1+15EKvByZNAiIjO6YWIiLqUFJKFBfXWjfCNZ8JPnGi2jrOz8/bYWeIxMQI+Ppyw/dZDMdEv9Gewj3IXJqJXQW7MLnfZCyYtgDdw7trV5DZDOzYYTud7+hR1eVi2DC1mW/yZKBLF+3qIyKiSyKlxPHj1Q47Q5SU1FnHBQfrzgvA6mOvXqHw9mZ70F/DcEx0iRqMDXh6w9N4bvNziAyMxJt/fhNTr5jqXLtxpQR27VIhecUKdVKfEKotnF4PZGQAMTFaV0lERM2YTGYcOVLpsDNEVVWjdVxkpH+z1mi2meCYmBDn+l3kYhiOiS7BjkM7kLk0E9lF2Zh+zXTMv3E+ooKjtC7r4qQE9u61zSjv368ev+oq26EjCQna1khE5EGMRhPy8srtOkPk5pahvr7JOq5btyC7AJyerjpDMAS3PYZjolaoaajBY6sew+vbXkfPiJ5499Z3MbrPaK3LujS5ubagvHu3euzyy21BOTVV2/qIiNxEXZ0RublldkcmHzxYjqYms3Vcr16h5x2UoYJweDjbdXYkhmOiFvos6zPMWjYLR0qO4K5hd+F/M/4XIf4hWpfVNg4fBlauVEsvvv9ePZaebgvKffuq5RhERHRBlZUNDjtDHD5cYe0M4e0t0Lt3uN1yiJSUCAQF+Wr7AggAwzHRryqrKcP9n9yPJd8uQXKXZCyavgiDkwdrXVb7OXZMdbwwGICvv1Yb/BITbUH5yisZlInIo505Y98ZIiurBMeP2zpD+Pp6IyUlwm45RFJSBPz8fDSsnn4NwzHRRazavQp3fnAniquK8cCoBzBvwjz46zzo7a3Tp1UPZYNB9VRuagJiY9VGPr0eGDgQ8OLOZyJyP1JKFBXV2AXg7OwSFBfbOkMEBemsp8M1D8Lx8WHw8eHPR1fEcEzkwKnKU7j7g7vxya5P0K9nPyyesRj94/prXZa2SkvVqXwGgzqlr7ER6NoVmDJFBeUhQwAfzoYQkWsxmyUKCiqsM8HNw3Blpa0zRHi4n8POED17hsLLi++muROGY6JmpJRYvmM57vv4PtQ01mDe+Hl4YNQD0Pl43glBF1VZCWzYoILyxo3qWOuoKNVDWa8Hhg9Xh5AQETkJo9GEQ4fK7ZZD5OSUoq7O1hmiS5dAu5Pi0tKi0KULO0N4CoZjIoujJUdx+39ux+Z9mzGw90AsnrEYqd3YseFX1dYCmzerzXzr1wNVVUBYGDBhggrKo0apY62JiDpAfX0TDhyw7wxx4EAZjEZbZ4jY2BCHnSEiI/nzytMxHJPHM5vNePurt/Gw4WFISDyX8RzuHHonvLiWtvXq64HPPlMzymvWAGVlQFAQMHasOp1v7FggOFjrKonIDVRXNyI7u8RuJjg/vwJms8ovXl4CCQlhdsshUlMjERzMd7fIMYZj8mi5J3Mxc+lMfJP3DUamj8S7t76LXtG9tC7LPRiNwJdfqqC8apXa3Ofvr2aS9Xo1sxwernWVROTkSkvrHHaGKCysso7R6byQnHx+Z4goJCdHwN+feyGodRiOySMZm4x4eevL+OfafyLQNxCv/OkVTL9mOteTtReTCfjmGxWUV64Ejh8HdDq1NlmvV2uVo6O1rpKINCKlxKlTtQ47Q5w6VWsdFxDgY13+0DwIJySEQafz1vAVkDthOCaP89PRn5C5NBM/Hf0J+v56vPnnN9E1rKvWZXkOsxn48Ufb6XyHD6t2cEOHqqA8ZQrQrZvWVRJROzCbJQoLKx12higvb7COCwvzswvAaWmRiIsLY2cIancMx+Qx6o31+Nf6f+H5zc8jOjgaC/68APor9FqX5dmkBH7+WYXkFSvUkdZCqP7JZw8diY3VukoiaqWmJjPy8x13hqipMVrHdeoU4LAzRLduQXwnjzTDcEwe4du8b5G5NBO5J3Px12v/ipf++BIigyK1LouakxLIyrLNKP/yi3r8D3+wBeXERG1rJKJzNDQ04eDBs50hbEE4N7cMjY0m67iYmBDrTHDzzhDR0YEaVk/kGMMxubXq+mo8uupRvPnFm4iNjMXCWxdi5O9Gal0WtUReni0o79ypHuvb1xaU09N5jDVRB6mpaUROTqndTPChQ+UwmVRWEAJISAi3Ww6RmhqJ0FA/jV8BUcsxHJPb2rJvC2Yvn43CskLcff3deGbyMwj2Zxsxl1RQoDbyGQzAd9+pWeaUFNUeTq8H+vVjUCZqA+Xl9Q47QxQUVFrH+Ph4ISkp3K4zREpKBAICeGASuT6GY3I7pTWl+PvHf8fSHUuR2jUVi2YswrWJ12pdFrWVoiLVGs5gUK3izGYgPt42ozxggNrgR0QOSSlRXFx7zlKIs2G4qKjGOs7f3wepqZF2M8GJieHsDEFujeGY3IphlwF3fXAXzlSfwcOjH8bj4x+Hv85f67KovZw5ow4bMRjU4SNGI9CjB5CRoYLyoEGAN3+Jk2eSUuLYsSq79cBZWSUoLa23jgsJ8bULwOnpUYiLC4W3N/+iSZ6H4ZjcQlF5EeZ8OAcrd69E/9j+WDxjMfrF9tO6LOpI5eXAunUqKG/eDDQ0AJ07q9Zwer1qFafjW77kfkwmMw4frrAE4DPnBOHqaltniKiogGYb4mxhuEePYHaGIGqG4ZhcmpQSS79birn/nYu6xjo8OfFJ/GPkP+DjzRORPFp1NbBxowrKGzYANTVAZCQwcaIKyiNGAH7cIESupbHRhLw8+84QOTmlaGiwdYbo3j3YLgCnp0ehUyd2hiBqCYZjcllHzhzB7OWzsTVrKwYnDcai6YuQ3DVZ67LI2dTVAVu2qKC8bh1QUQGEhADjx6sNfaNHA4EMDeQ8amuNyM217wyRl1eOpiazdVx8fJjDzhDh4VxKRvRbMByTyzGZTVjwxQI8uupRCAg8r38efxvyN3hxExb9msZG4PPPVVBevRooKVHBeMwYNaM8bhwQGqp1leQhKisbrEshms8GHzlSgbO/fr29BRITI+xmglNSIhAU5KvtCyByUwzH5FKyi7Ixc+lMfHfoO4zuMxrv3vIuYqN4ghpdgqYmYPt2dTLfqlXAyZOAry8wcqQKyhMnqqUYRL/RmTO1zTpC2GaDjx+vto7x8/NGSop9Z4ikpAj4+nJTKVFHYjgml2BsMuKFLS/gqfVPIdgvGK/+6VXccvUt3ERCbcNsVv2Tzx46UlgI+PgA11+vgvLkyWpzH9EFSClx4kS1w84QZ87UWccFBekcdoaIjw9jZwgiJ8FwTE5vd8Fu3Pbv27Dn2B7ceOWNeP3m19EltIvWZZG7khL4v/+zBeW8PNU3efBgFZQzMlS7OPJIZrPEkSMV1uDbPAhXVjZax0VE+DvsDBETEwIvL/6lnsiZMRyT06prrMOT657ES5++hM4hnfHWtLcw+fLJWpdFnkRKYO9etfTCYACystTjV19tO3QkPl7bGqldGI0mHDpUbndQRm5uKerqmqzjunYNctgZonPnQL6zReSiGI7JKW0/sB0zl87EwdMHkTkoEy/98SWEB4ZrXRZ5upwc24zyTz+px/r3twXllBRt66NWq69vatYZwtYj+ODBMhiNts4QcXGhdgE4LS0KERHsDEHkbhiOyalU1lXikZWP4K0v30J8dDzem/4ehqcN17osInv5+cDKlWpW+Ycf1GO/+50KyVOnAn36AJw5dBpVVY3IySmx2xh3+HAFzGb1u87LS6B373C7meDU1EgEB7MzBJGnYDgmp7Fp7ybc/p/bcazsGO4dfi+envw0gvyCtC6L6NcdO6aCssEAfP21Wo6RlGSbUb7iCgblDlJSUudwPXBhYZV1jE7ndU5niLMhOCkpAv7+PECIyNMxHJPmSqpLMPfjuVj+/XKkd0vH4hmLcXXvq7Uui+jSnDqleigbDMC2bYDJBMTFqY18ej1wzTVqgx9dMiklTp6sOacjxNkgfPp0rXVcYKAPUlOj7DbGJSSEw8eH/w2IyDGGY9KMlBKf/N8nmPPhHJTVluHRMY/i0bGPwk/HY33JTZSUAGvXqqC8das6hKRbN1tQHjxYtYwjh8xmiaNHK+0CcFZWCSoqGqzjwsL8HHaGiI0NZWcIImo1hmPSxInyE7jz/Tux5uc1uDLuSiz+y2L0jemrdVlE7aeyEli/XgXlTZvUsdbR0aqHsl6veir7eua61qYmM/Lz7TtD5OSUoLbW1hmic+dAh50hunYNYmcIImozDMfUoaSUWPLNEvzjk3+goakB/5r0L9x3w33w8ebsGXmQmhpg82a1mW/9eqC6GggPByZMUEF55EggIEDrKttcQ0MTDhwos84Enw3DBw6UobHRZB3Xs2dIs44QtiAcFeV+14SInA/DMXWY/OJ8zFo2C9tytmFI8hAsmrEIiZ0TtS6LSFv19WrJhcGglmCUlQHBwcC4cSoojxmjPnch1dWNyMkptVsOcehQubUzhBBAQsK5nSHS06OQmhqFkBDPnEEnIufAcEztzmQ24Y1tb+CxVY/B28sbL059EbMGz4IXNyURnctoBL74QgXlVauA4mLA3x8YPVoF5QkTgLAwrau0KiurbxaAbUG4oKDSOsbHxwvJyRHNAnA00tIikZwcgYAAnYbVExE5xnBM7Wr/8f3IXJqJHw7/gHG/H4d3bnkHMZExWpdF5PxMJtUWzmBQbeJOnAB0OuCGG1Qf5UmTgKiodi9DSonTp2sddoY4ebLGOs7f3wepqZF2G+N69w6HTufd7nUSEbUVTcKxEGIJgPEATksp+7TkexiOXUtjUyOe3/w8/rX+XwgNCMXrN72OmwfczE0zRJfCbFYHjZw9ne/IEcDbGxg6VM0oT5kCdO36m55CSonCwiqHnSHKyuqt40JDfR2uB46LC4W3N98NIiLXp1U4vg5ANYBlDMfuZ+fhnchcmom9x/fi5gE347WbXkOnkE5al0XkHqRUR1cbDGpD34EDagHvtdeqoJyRgfe/jsVjjwFHjwKxscAzzwDTpqlvN5nMOHy44pyZYNUZohTV1Ubr00RHB5x3VLIKwt27B/MvuUTk1jRbViGE6AVgPcOx+6htqMW8tfMwf+t8dAvrhrdveRsTLpugdVlE7ktKYP9+24zy3r0AgJ1eA/CJWQ8D9MhHb+h0Jlxxxc+ord2H3NxSNDTYOkP06BHscCa4U6dArV4VEZGmGI6pTXyZ+yVmLZuFvNN5mH3dbLygfwFhgc6zcYjIXdXWGq2dIU598zPKFuVhYtM6/AHq5+XPuAwG6LHaazTixhSeszEuNTUSYWE8dIeIqDmnDsdCiNkAZgNAbGzsFQUFBe1WD12aitoKPGR4CO9ufxe9O/XGe9Pfw7DUYVqXReR2KioaHHaGOHKkAmd/VHt7C5hMfwcgEIcjyMBK6GHAtfhODUhLU0sv9HrgssvUcgwiIjqHU4fj5jhz7Hw2/LIBty+/HUUVRZg7Yi6emvgUAv34VizRb1Fc7LgzxIkT1dYxfn7eSE2NtFsOkZgYgeRkb5w/j9ANJ5AZuQr/uswAfPWV2uCXkGALygMGMCgTEVlcLBzzyDJyqLiqGPd9dB8++PED9OnRByvvXIkB8QO0LovIZUgpcfx4td1JcVlZJSgpqbOOCw7WIS0tCiNGxFkCsOoRHB8fdsHOEM88A8yeDdTW2h6rCOyO1NfvAqbdpXonr1mj1ii/8grw4otATAyQkaFaxA0cqDphEBGRnfbsVvEhgKEAogGcAjBPSrn4Yt/DmWPtSSnx0Y8f4Z6P7kFFXQUeH/c4Hh7zMHx9eJoVkSMmkxkFBZV2nSGys0tRVdVoHRcZ6e+wM0RMTMgldYZ4/31csFvFOcrKgHXrVFDesgVoaAC6dFGt4fR6YMgQ1VuZiMiD8BAQapFjpcdwx/t3YP0v6zEgfgAWz1iMPj1atCKGyO0ZjSbk5ZXbLYfIySlFfX2TdVy3bkF2ATg9XXWG0Lw9WlUVsHGjCsobNqip58hIddiIXq8OH/Hj5j0icn8Mx3RRZrMZi75ZhAdWPACjyYhnJj+De4bfA28vvu1Knqeuzojc3DK7jXEHD5ajqclsHRcXF9rspDhbGA4P99ew+laoq1MzyStWqJnlykogNFQdX63XA6NGAYHcX0BE7onhmC4o73QeZi2bhS9zv8SwlGF4b/p76N25t9ZlEbW7qqpGh50h8vPLrZ0hvLwEEhPD7QJwamokgoLcaKlRQwPw+edqRnn1aqC0VAXjsWNVUB43DggJ0bpKIqI2w3BMdkxmE1797FU8seYJ6Lx1ePmPLyNzUKb2b/sStbGSkjqHnSGOHauyjvH19UZKSoTdcoikpAj4+XnYvuWmJtXtYsUKYNUq4NQptdRi5Ei1mW/CBCAiQusqiYh+E4ZjOsfeY3uRuTQTO4/sxMTLJuKtaW+hR0QPrcsiumRSShQV1Zy3IU59LC62dYYICtIhNTXSbiY4ISEcPj6OO0N4NJMJ+O47NaO8ciVQWAj4+ADDh6sZ5cmTgU48Np6IXA/DMQEAGowNeHbjs3h207OICIzAGze/gRuvvJGzxeQyzGaJo0cr7QJwdnYpKioarOPCw/3OCcBnQ3DPnqHw8uL9fkmkBHbutB1jfegQ4OUFXHedCsoZGUD37lpXSUTUIgzHhB/yf0Dm0kzsP7Eft1x9C1658RVEh0RrXRaRQ01NZhw65KgzRAlqa22dIbp0CTxvKUQ00tOj0KWLE3SGcGdSAr/8opZeGAxAdrZ6fOBAW1Du1UvTEomILobh2IPVNNTgidVP4NXPX0WP8B5455Z3MK7vOK3LIgIA1Nc34cAB+84QBw6UwWi0dYbo2TPEYWeIyMgADasnq+xs24zyzz+rx664wnY6X3KytvUREZ2H4dhDbcvehlnLZyG/OB93DL0Dz2U8h9CAUK3LIg9UXd2InJxSu4My8vMrYDarn0FeXgIJCWEOOkNEISTEjTpDuLtDh9T6ZIMB+OEH9VifPiokT50K/O53PMaaiDTHcOxhymvL8cCKB7Do60VI6pyE96a/hyEpQ7QuizxAaWkdsrNL7TbGHT1q6wyh03khObl5Zwj1MTk5Av7+HtYZwt0VFtqC8jffqOUYycm2GeX+/RmUiUgTDMceZO3Pa3HH+3fgZMVJ3D/yfvxz4j8R4Mu3nqntSClx6lStw84Qp07VWscFBPjYdYZIT49CQkIYdDoeMONxTp5UPZQNBuCLL1QnjF691PrkqVOBq65SG/yIiDoAw7EHOF15Gvd8dA8+3vkx+sb0xeIZi3FlL4f/zYlaREqJwsIqh50hysrqreNCQ30ddoaIiwtjZwhyrKQEWLNGBeWtWwGjUXW6yMhQM8qDBwPe/AsUEbUfhmM3JqXE+z+8j3s/uhfVDdV4YtwTeGj0Q9D56LQujVyEyWRGfn7FeQFYheCaGqN1XKdOAecFYPWxW7cgdoagS1dRAaxfr4Lypk1Afb3qnTx5sgrK118P6PjzjIjaFsOxmyosLcTf/vM3bNy7EVcnXI3FMxYjvXu61mWRk2poaMLBg2XWE+LOBuHc3DI0Npqs42JiQqwnxDXfGBcdHahh9eQRampUQF6xAtiwAaiuBsLDgUmTVFAeMQLw99e6SiJyAwzHbsZsNuPd7e/iIcNDMJlNeHbKs5hz/Rx4e/FtSAJqahqRm1tm1xni0KFymEzq/3chgPj48ztDqBAcGuqn8SsggppB/vRTNaO8di1QXg4EBwPjx6ugPGYMEBSkdZVE5KIYjt3IwVMHMXPZTGw/sB03pN2AhbcuRHyneK3LIg2Ul9c77Axx5EildYyPjxeSksLtOkOkpEQgIIBvVZOLaGxUm/gMBrWpr7gYCAgARo9WQXn8eCAsTOsqiciFMBy7gSZTE+ZvnY95a+fBX+eP+TfOx18G/oVrPd2clBLFxbXWjXDNZ4OLimqs4/z9fZCSEmG3Ma5373D4+vIdBXIjTU3A11+roLxyJVBUBPj6qiUXej0wcSIQFaV1lUTk5BiOXdyewj247d+3YffR3Zhy+RQs+PMCdAvvpnVZ1IaklDh+vNphZ4iSkjrruOBgnYPOEFHo1SsU3t5sg0UexmwGvv/edjpfQYHqcjFsmArKU6YAXbpoXSUROSGGYxfVYGzA0xuexnObn0NkYCQWTFsAfX89Z4tdmMlkxpEjlQ47Q1RVNVrHRUb6W8Nv8zDco0cw//sTOSIlsHu3CskrVgAHD6rF9YMGqT7KGRlATIzWVRKRk2A4dkE7Du1A5tJMZBdlY/o10zH/xvmICuZbha6isdGEvDzHnSHq65vowS58AAAQDElEQVSs47p3D3bYGaJTp0CGYKJLJSWwb59tRnnfPvX4VVfZTudLSNC2RiLSFMOxC6mur8bjqx/H69teR8+Innj31ncxus9orcuiC6irMzrsDJGXV46mJrN1XK9eoQ47Q4SHsy0VUbs7cMAWlHftUo/162cLymlp2tZHRB2O4dhFbM3aitnLZuNIyRHMGTYHz2Y8ixD/EK3LIgCVlQ0OO0McPlyBs/8LeXsLJCZG2M0Ep6REICjIV9sXQETK4cNqI5/BAOzYoR5LT7cF5b591XIMInJrDMdOrqymDPd/cj+WfLsEKV1TsGj6IgxKGqR1WR7pzBnHnSGOH6+2jvH19XbYGSIxMRx+fj4aVk9ErXL8OLBqlQrK27erDX69e9uC8h/+wKBM5KYYjp3Yqt2rcOcHd6K4qhgPjnoQ/zPhf+Cv41vt7UlKiRMnqu0CcHZ2CYqLbZ0hgoJ0SEuLtOsMER8fBh8fdoYgciunTwNr1qig/PnnqmVcbKzayKfXAwMHAl78/57IXTAcO6GTFSdx94d3Y8WuFejXsx8Wz1iM/nH9tS7LrZjNEgUFFQ5ngisrbZ0hIiL8rRvhms8Gx8SEwMuLs0ZEHqesTJ3KZzCoU/oaGoCuXVVrOL0eGDIE8OG7RESujOHYiUgpsXzHctz38X2obazFvAnzcP/I+6Hz4Wlll8poNOHQoXJrAD4bgnNySlFXZ+sM0aVLoF1rtLS0KHTpws4QRHQBVVXAhg0qKG/cCNTWqkNGJk1SLeKGD1eHkBCRS2E4dhIFJQW4ffnt2LJ/C65NvBaLpi9CardUrctyGfX1TThwoMwSgM9Yw/DBg2UwGm2dIWJjQxx2hoiMDNCweiJyebW1wJYtqo/yunUqOIeFARMmqBnlUaPUsdZE5PQYjjVmNpvx9ldv42HDw5CQeC7jOdw59E54cf2aQ1VVjcjJKbFbDpGfXwGzWd2vXl4CvXuH2y2FSE2NRHAwZ3GIqJ01NACffaZmlNesAUpLgaAgYOxYFZTHjQOCg7WukogugOFYQ7knczFz6Ux8k/cNRv1uFN699V3ERcVpXZZTKC2tc7geuLCwyjpGp/NCcrJ9Z4ikpAj4+3PNHxE5AaMR+OorFZRXrlSb+/z8gNGjVVCeMAEID9e6SiJqhuFYA8YmI1769CU8ue5JBPoG4pU/vYLp10z3uLWtUkqcPFljCcBnzlkXfPp0rXVcQICPdflD83XBvXuHszMEEbkOkwn49ltbUD52DNDp1NpkvV6tVe7USesqiTwew3EH++noT8hcmomfjv6EqVdMxRs3v4GuYV21Lqtdmc0ShYWVDmeCy8sbrOPCwvwcdoaIjQ1lZwgici9mM7Bzp+10vvx81Q5uyBC1mW/KFKBbN62rJPJIDMcdpN5Yj6fWPYUXtryA6OBovDXtLWT0z9C6rDbV1GRGfr59Z4js7BLU1to6Q3TuHGgXgNPSotCtW5DHzZ4TEUFKYM8eFZJXrAByctQBIwMHqhnljAwgjkvuiDoKw3EH+DbvW2QuzUTuyVz89dq/4uU/voyIoAity7pkDQ1NOHiwrFkAVmH4wIEyNDaarONiYkLsZoLT0iIRHR2oYfVERE4uK8s2o7xnj3rsyittp/MlJWlbH5GbYzhuR1X1VXh05aNY8OUCxEXGYeH0hRiRPkLrslqspqYROTmldsshDh0qh8mk7g0hgISEcLv1wKmpkQgN9dP4FRARubi8PLU+2WAAfvxRPda3ry0op6fzGGuiNsZw3E627NuC2ctno7CsEHdffzeemfwMgv2ds3VPWVm9ZfnDuRvjCgoqrWN8fLyQlBRud1BGcnIEAgJ4SAkRUbs7etQWlL/9Vi3HSEmxBeXLL2dQJmoDDMdtrLSmFHM/notlO5YhtWsqFs9YjIGJA7UuC1JKnD5da90I13xd8MmTNdZx/v4+SE2NtFsOkZgYDp3OW8NXQEREVkVFwOrVKih/+aXqhBEfr9YnT50KDBigNvgRUasxHLehFbtW4K7370JpbSkeGvUQHh//OPx1/h1ag5QSx45V2QXg7OwSlJbWW8eFhPg67AwRFxcKb2/+QCUichlnzqjDRgwGdfiI0Qj06KGCsl4PDBoEeHNyg6ilGI7bQFF5EeZ8OAcrd69E/9j+WDxjMfrF9mvX5zSZzDh8uOKctmjZ2aXIzi5BdbXROi4qKqDZUgjbuuDu3YPZGYKIyN2UlwPr16ugvHkzUF8PdO4MTJ6sZpSHDlW9lYnoghiOfwMpJf793b/x9//+HfXGejw58Un8fcTf4ePddqezNTaakJdn3xkiN7cUDQ22zhDduwc7nAnu1ImdIYiIPFJ1NbBpk2oPt2EDUFMDRESow0b0emDECHVaHxGdg+H4Eh0uPozb/3M7tmZtxeCkwVg0fRGSuyZf8r+vttaI3Fz7zhAHD5ad0xmiV68wu84QaWlRCAvjDzgiIrqAujrg00/VjPLatUBFBRASAowfr4LymDFAICdTiACG41YzmU1Y8MUCPLLyEXgJL7ww9QXcft3t8GrhxoeKigbrwRjNg/CRIxU4e7m9vQUSEyPslkOkpEQiMJBvhxER0W/Q2Ahs26aC8urVas1yQAAwdqwKyuPGAaGhWldJpBnNwrEQYjSA1wB4A1gkpXzuYuOdIRxnF2Ujc2kmdhzagTF9xuCdW95BbFSsw7HFxY47Q5w4UW0d4+fnjZQU+84QSUkR8PXl5gkiImpnTU3A9u0qKK9apbpg+PoCI0eqoDxxIhAZqXWVRB1Kk3AshPAGcADACADHAOwEcLOUMutC36NlODY2GfHClhfw1PqnEOwXjNdueg3TrpoGADhxotphZ4gzZ+qs3x8UpHO4Hjg+PoydIYiIyDmYzcCOHbbT+Y4eBXx8gGHDVFCePBno0kXrKonanVbh+BoA/5RSjrJ8/ggASCn/90Lfo1U43lWwC7f9+zb8cuwXDOo+Ftf7342jB6S1M0RlZaN1bESEv8POEDExIewMQURErkNKYNcuFZJXrFAn9Xl5AYMHq6CckaHaxRG5Ia3C8VQAo6WUMy2f3wrgKinlnAt9T0eHYyklBt5zK76v+xCoCwa+nQIU9AEAdO0aZBeA09Ki0LlzIEMwERG5FymBvXttM8r796vHr74aePFF1UeZyI1cLBy3XT8yB8/r4DG7JC6EmA1gNgDExjpe29tehBCobqxETM0wTOx2L/o/0QtpaZFIS4tCRETHHuxBRESkGSGAvn3VnyefBHJzbUE5KEjr6og6lMcvqzCbzS3uQkFEREREru9iM8ftmQp3AkgSQsQLIXwB3ARgbTs+3yVhMCYiIiKis9ptWYWUskkIMQfAFqhWbkuklPvb6/mIiIiIiH6r9lxzDCnlRgAb2/M5iIiIiIjaCtcUEBERERFZMBwTEREREVkwHBMRERERWTAcExERERFZMBwTEREREVkwHBMRERERWTAcExERERFZMBwTEREREVkwHBMRERERWTAcExERERFZMBwTEREREVkwHBMRERERWTAcExERERFZMBwTEREREVkwHBMRERERWQgppdY1WAkhigEUaPDU0QDOaPC8rorXq3V4vVqH16t1eL1ah9er9XjNWofXq3W0ul5xUspOjr7gVOFYK0KI/5NSXql1Ha6C16t1eL1ah9erdXi9WofXq/V4zVqH16t1nPF6cVkFEREREZEFwzERERERkQXDsbJQ6wJcDK9X6/B6tQ6vV+vwerUOr1fr8Zq1Dq9X6zjd9eKaYyIiIiIiC84cExERERFZeFQ4FkKMFkLkCiHyhBAPO/i6EEK8bvn6L0KI/lrU6SxacL2GCiEqhBA/W/78jxZ1OgMhxBIhxGkhxL4LfJ33VjMtuF68t5oRQvQUQnwhhMgWQuwXQtzrYAzvMYsWXi/eYxZCCH8hxI9CiD2W6/WkgzG8vyxaeL14f51HCOEthPhJCLHewdec6v7y0fLJO5IQwhvAAgAjABwDsFMIsVZKmdVs2BgASZY/VwF42/LR47TwegHA11LK8R1eoPP5N4A3ASy7wNd5b53r37j49QJ4bzXXBOAfUsrdQogQALuEEFv58+uCWnK9AN5jZzUAuF5KWS2E0AH4RgixSUr5fbMxvL9sWnK9AN5f57sXQDaAUAdfc6r7y5NmjgcAyJNS5kspGwF8BGDSeWMmAVgmle8BhAshunV0oU6iJdeLLKSU2wGUXmQI761mWnC9qBkpZZGUcrfln6ugfsH0OG8Y7zGLFl4vsrDcM9WWT3WWP+dvSOL9ZdHC60XNCCFiAIwDsOgCQ5zq/vKkcNwDQGGzz4/B/odlS8Z4ipZei2ssby1tEkL8rmNKc0m8t1qP95YDQoheAC4H8MN5X+I95sBFrhfAe8zK8pb3zwBOA9gqpeT9dREtuF4A76/mXgXwIADzBb7uVPeXJ4Vj4eCx8/+m15IxnqIl12I31PGLlwF4A8Dqdq/KdfHeah3eWw4IIYIBGADcJ6WsPP/LDr7Fo++xX7levMeakVKapJT9AMQAGCCE6HPeEN5fzbTgevH+shBCjAdwWkq562LDHDym2f3lSeH4GICezT6PAXDiEsZ4il+9FlLKyrNvLUkpNwLQCSGiO65El8J7qxV4b9mzrG00AHhfSrnSwRDeY8382vXiPeaYlLIcwJcARp/3Jd5fDlzoevH+Ose1ACYKIY5ALdG8Xgjxn/PGONX95UnheCeAJCFEvBDCF8BNANaeN2YtgOmWXZNXA6iQUhZ1dKFO4levlxCiqxBCWP55ANT9VNLhlboG3lutwHvrXJZrsRhAtpRy/gWG8R6zaMn14j1mI4ToJIQIt/xzAIAbAOScN4z3l0VLrhfvLxsp5SNSyhgpZS+oLLFNSnnLecOc6v7ymG4VUsomIcQcAFsAeANYIqXcL4T4m+Xr7wDYCGAsgDwAtQD+qlW9Wmvh9ZoK4A4hRBOAOgA3SQ89VUYI8SGAoQCihRDHAMyD2qTBe8uBFlwv3lvnuhbArQD2WtY5AsCjAGIB3mMOtOR68R6z6QZgqaVLkReA/0op1/P34wW15Hrx/voVznx/8YQ8IiIiIiILT1pWQURERER0UQzHREREREQWDMdERERERBYMx0REREREFgzHREREREQWDMdERC5GCNFTCHFYCBFp+TzC8nmc1rUREbk6hmMiIhcjpSwE8DaA5ywPPQdgoZSyQLuqiIjcA/scExG5IMvxyLsALAEwC8DlUspGbasiInJ9HnNCHhGRO5FSGoUQDwDYDGAkgzERUdvgsgoiItc1BkARgD5aF0JE5C4YjomIXJAQoh+AEQCuBjBXCNFN45KIiNwCwzERkYsRQgioDXn3SSmPAngRwEvaVkVE5B4YjomIXM8sAEellFstn78FIFUIMUTDmoiI3AK7VRARERERWXDmmIiIiIjIguGYiIiIiMiC4ZiIiIiIyILhmIiIiIjIguGYiIiIiMiC4ZiIiIiIyILhmIiIiIjIguGYiIiIiMji/wHEHU7G7wh25AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.figure(figsize=(12, 8))\n",
    "plt.plot([0, 4], [0, 2], color='darkblue', label='beta_1 = 1/2')\n",
    "plt.plot(2, 1, 'go', color='blue', label='(X, Y) = (2, 1)')\n",
    "plt.plot([0, 4], [0, 8], color='darkgreen', label='beta_1 = 2')\n",
    "plt.plot(1, 2, 'go', color='green', label='(X, Y) = (1, 2)')\n",
    "plt.plot([0, 3], [3, 0], color='red', label='E[Y | X = x]')\n",
    "plt.legend()\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, what we want to do to make $r(x) = \\mathbb{E}(Y | X = x)$ is to ensure the distribution generates points that have a lower $X$ for a larger $Y$, by selecting a smaller value of $X$ when $Y$ is larger.\n",
    "\n",
    "Let's choose a joint probability distribution as $\\mathbb{P}(\\beta_1, x) = \\mathbb{P}(\\beta_1) \\mathbb{P}(x | \\beta_1)$, and let's choose distribution of $x | \\beta_1$ as a discrete probability mass function with all its mass on a single value:\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(x | \\beta_1) = \\begin{cases}\n",
    "1 & \\text{if } x + y = 1, \\text{ where } y = \\beta_1 x\\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Now, by construction all points $(X, Y)$ lie on the line $X + Y = 1$, so $\\mathbb{E}(X + Y) = 1$ and $\\mathbb{E}(Y | X = x) = 1 - x$, which is decreasing in $x$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}