causal-inference-for-the-brave-and-true/14-Panel-Data-and-Fixed-Effects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 14 - Panel Data and Fixed Effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import style\n",
    "from matplotlib import pyplot as plt\n",
    "import statsmodels.formula.api as smf\n",
    "import graphviz as gr\n",
    "from linearmodels.datasets import wage_panel\n",
    "\n",
    "%matplotlib inline\n",
    "pd.set_option(\"display.max_columns\", 6)\n",
    "style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이전 챕터에서, 실험군과 대조군(각각 POA와 FLN)과 2개의 기간(Intervention 전후)이 있는 기본적인 DID에 필요한 가정과 적용방법에 대해 배웠습니다. 만약에 더 많은 기간과 3개 이상의 그룹인 경우에 대해 인과효과를 추정하기 위해서는 어떻게 해야 할까요? DID의 설정은 인과추론을 적용하는데 있어 매우 일반적이며, 강력하다보니 이에 대한 자체적인 용어도 있습니다. 바로 패널 데이터(Panel data)인데요, 패널은 여러 기간에 걸쳐 동일한 단위(unit)를 반복적으로 관찰한 경우입니다. 패널 데이터는 여러 해에 걸쳐 여러 도시 또는 주에 대한 데이터를 추적할 수 있는 정부 정책 평가에서 많이 발생합니다. 패널 데이터는 정부 관련 뿐만 아니라, 기업들이 몇 주 또는 몇 달에 걸쳐 유저 데이터를 추적하는 산업에서도 매우 일반적으로 볼 수 있어요.\n",
    "\n",
    "패널 데이터 구조를 활용하는 방법을 이해하기 위해, 먼저 POA(Porto Alegre)시에 광고판을 배치할 때(Treatment)의 영향을 추정하고자 했던 DID 예제를 계속 진행해보도록 할게요. 이 예시에서는 오프라인 마케팅 전략이 투자 상품의 활용도를 높일 수 있는지 파악하고자 했습니다. 특히, 광고판을 설치하면 투자 계좌의 예금이 얼마나 증가하는지 알고 싶었던 것이죠.\n",
    "\n",
    "이전 챕터에서, 광고판을 배치하지 않았다면 Porto Alegre에 어떤 일이 일어났을지에 대한 대체 전략(Imputation strategy)으로 DID 추정량을 활용했어요. Porto Alegre **개입 후**(광고판 배치)에 대한 반사실 결과 $Y_0$는 **개입 전** Porto Alegre의 예금수에 성장 요인(growth factor)을 더한 값으로 볼 수 있다고 했었습니다. 참고로 이 성장 요인은 대조군인 Florianopolis(FLN) 시에서 추정되었습니다. 몇 가지 표기법 요약을 통해, 어떻게 반사실 결과를 추정할 수 있는지 다음과 같이 나타내보겠습니다.\n",
    "\n",
    "$$\n",
    "\\underbrace{\\widehat{Y_0(1)|D=1}}_{\\substack{\\text{POA outcome after intervention}} \\\\ \\substack{\\text{if no intervention had taken place}}} \n",
    "= \\underbrace{Y_0(0)|D=1}_{\\substack{\\text{POA outcome}} \\\\ \\substack{\\text{before intervention}}} \n",
    "+ \\big( \\underbrace{Y_0(1)|D=0}_{\\substack{\\text{FLN outcome after}} \\\\ \\substack{\\text{intervention in POA}}} \n",
    "- \\underbrace{Y_0(0)|D=0}_{\\substack{\\text{FLN outcome before}} \\\\ \\substack{\\text{intervention in POA}}} \\big)\n",
    "$$\n",
    "\n",
    "여기서 $t$는 시간을 나타내고, $D$는 처치($t$를 사용하기 때문에)를 나타내며, $Y_D(t)$는 $t$ 기간의 처치 $D$에 대한 잠재적 결과(potential outcome)를 나타냅니다(예시 : $Y_0(1)$는 기간이 1 일 때, 대조군에 대한 결과). 이제 counterfactual에 대한 대체적인 잠재적 결과를 취하면 다음과 같이 POA에 대한 처치 효과(ATT)를 표현할 수 있어요.\n",
    " \n",
    "$$\n",
    "\\widehat{ATT} = \\underbrace{Y_1(1)|D=1}_{\\substack{\\text{POA outcome} \\\\ \\text{after intervention}}} - \\widehat{Y_0(1)|D=1}\n",
    "$$\n",
    "\n",
    "즉, 광고판을 POA에 배치하는 효과는 광고판을 배치한 후, POA에서 측정한 결과에서 광고판을 배치하지 않았다면 일어날 일에 대한 추정치를 뺀 것에 해당해요. 또한 DID의 핵심은 위에서 언급된 반사실을 추정하기 위해, POA에서의 예금의 증가가 FLW의 예금 증가와 일치해야 한다는 사실에서 비롯된다는 사실을 명심해야 해요. 이 부분이 바로 DID에 핵심인 평행 추세 가정(parallel trends assumption)입니다. 평행 추세 가정은 나중에 매우 중요한 부분을 차지하기 때문에 이번에 확실하게 배워봅시다.\n",
    "\n",
    "## Parallel Trends\n",
    " \n",
    "평생 추세 가정을 보는 한 가지 방법은 독립성 가정(Independence assumption)입니다. 초반에 학습한 챕터를 떠올려보면, 독립성 가정은 처치 할당(treatment assign)이 잠재적 결과로부터 독립적이어야 합니다:\n",
    "\n",
    "$$\n",
    "Y_d \\perp  D\n",
    "$$\n",
    "\n",
    "해당 가정은 더 높은 결과(효과 추정에서 상승 편향을 일으킬 수 있음) 또는 더 낮은 결과(하향 편향을 유발함)를 가진 대상에 더 많은 treatment를 주지 않는다는 것을 의미해요. 조금 구체적으로 설명하자면, 마케팅 매니저가 이미 은행 예금이 매우 높은 도시에만 광고판을 추가하기로 결정했다고 가정해 볼게요. 그렇게 하면 나중에 광고판이 있는 도시가 더 많은 예금을 만들었다고 자랑할 수 있게 됩니다. 물론 마케팅 캠페인은 성공적이었겠죠. 도덕적 논의를 제쳐두고, 이 부분이 독립성 가정을 위반한다는 것을 알 수 있을 거에요: $Y_0$가 높은 도시에 treatment를 주는 것이죠. 또한 독립성 가정의 자연스러운 확장은 조건부 독립에 대한 가정으로 이어집니다. 잠재적 결과가 처음에는 treatment에 종속되지만 교란변수 $X$를 condition하게 된다면 독립적 조건부 독립 가정이 됩니다.\n",
    "\n",
    "$\n",
    "Y_d \\perp D | X\n",
    "$\n",
    " \n",
    "그런데, 독립성 가정이 DID와 평행 추세 가정과 정확히 어떻게 연결될까요? 클래식한 독립성 가정에서 treatment assignment 잠재적 결과의 **수준(levels)**과 관련될 수 없다고 되어 있는 경우, 평행 추세 가정은 treatment assignment이 잠재적 결과의 시간에 따른 **성장(growth)**과 관련될 수 없다고 말하고 있어요. 평행 추세 가정을 작성하는 한 가지 방법은 다음과 같습니다.\n",
    "\n",
    "$\n",
    "\\big(Y_d(t) - Y_d(t-1) \\big)  \\perp D\n",
    "$\n",
    " \n",
    "조금 덜 수학적으로는, 이 조건부 독립 가정은 결과 수준이 높거나 낮은 대상에 treatment를 할당하는 것이 괜찮다고 이야기하고 있습니다. 이 가정에서는 결과가 어떻게 증가하는지에 따라 treatment를 실험 대상에 할당할 수 없게 됩니다. 가정은 앞서 살펴본 광고판 예시에서 원래 예금 수준이 높은 도시에만 광고판을 배치하는 것이 좋다는 것을 의미합니다. 반면, 예금이 가장 많이 증가하는 도시에만 광고판을 배치할 수 없다고 말하고 있어요. DID가 대조군의 추세와 함께 실험군의 counterfactual 추세를 input으로 한다는 것을 기억한다면 해당 가정은 많은 의미가 있어요. 통제하에 있는 실험군의 추세가 대조군의 추세와 다른 경우, 앞서 말한 가정에서 문제가 생기게 됩니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Controlling What you Cannot See\n",
    " \n",
    "성향 점수, 선형 회귀 및 매칭과 같은 방법은 랜덤하지 않은 데이터(non-random data)의 교란변수를 통제하는 ​​매우 좋은 방법입니다. 하지만, 해당 방법들은 중요한 가정인 conditional unconfoundedness에 의존합니다.\n",
    " \n",
    "$\n",
    "(Y_0, Y_1) \\perp T | X\n",
    "$\n",
    "\n",
    "한마디로 말하면, 해당 방법들은 모든 교란변수를 알고있고 측정이 되어야, 교란변수를 조건부로 두고 random assignment와 같이 treatment를 구성 할 수 있습니다. 그런데, 이 부분에 대한 한가지 중요한 문제는 교란변수를 항상 측정(관측) 할 수 없는 경우가 존재합니다. 예를 들어, 결혼이 남성의 수입에 미치는 영향을 알아내는 고전적인 노동 경제학 문제를 생각해 보시죠. 기혼 남성이 미혼 남성보다 더 많이 번다는 것은 경제학에서 잘 알려진 사실입니다. 그러나 이 관계가 인과관계인지 아닌지는 분명하지 않아요. 교육 수준이 높은 남성이 결혼할 가능성 및 높은 수입을 가질 가능성이 둘다 높을 수 있어요. 즉, 이 부분은 교육이 수입에 대한 결혼의 영향을 교란효과임을 의미합니다. 해당 경우에는 교란변수에 대해, 실험에 참여한 사람의 교육 수준을 측정하고 이를 통제하는 ​​회귀 분석을 실행해볼 수 있습니다. 그러나, 교육 이외의 또 다른 교란변수는 beauty (외모)가 될 수 있습니다. 더 잘생긴 남자는 결혼할 가능성 및 높은 보수를 주는 직업을 가질 가능성이 더 높을 수 있어요. 아쉽게도, 교란변수로 볼 수 있는 외모는 지능과 같이 정량화하기 어려운 특성 중 하나입니다. \n",
    "\n",
    "측정되지 않은 교란 변수가 있으면 bias가 존재하기 때문에, 인과효과 추정에 어려운 상황에 놓이게 됩니다. 교란변수를 효과적으로 처리하는 한 가지 방법은 이전에 다룬 것처럼 도구변수를 사용하는 것입니다. 하지만, 좋은 도구변수를 만드는 것은 쉬운일이 아니며, 많은 창의성을 요구해요. 그래서 해당 챕터에서는 도구변수 대신, 패널 데이터 구조를 활용해 보겠습니다.\n",
    "\n",
    "앞서 어떻게 패널 데이터가 unconfoundedness 가정을 평행 추세 가정으로 대체 가능하게 해주는 부분을 살펴 보았습니다. 그러나 패널 데이터를 이용한 방법이 측정되지 않은 교란변수를 다루는데 어떻게 도움이 될까요? 먼저, 시간이 지남에 따라 반복된 관찰이 있는 설정을 나타내는 인과 그래프를 보도록 할게요. 이 부분에서 4개의 기간에 걸쳐 동일한 관찰을 추적합니다. 결혼(처치)과 소득(결과)은 시간이 지남에 따라 변합니다. 구체적으로, 결혼은 기간 3과 4에 나타나고 (0에서 1로) 같은 기간에 소득이 증가합니다. 측정할 수 없는 교란변수 beauty는 모든 기간에 걸쳐 동일합니다 (과감한 표현이지만 기간이 불과 몇 년정도면 사람의 외모는 일반적으로 크게 변하지 않으니 합리적이라고 볼 수 있어요). 그렇다면 소득이 늘어나는 이유가 단순히 교란변수인 beauty의 증가 때문이 아니라 결혼 때문이라는 것을 어떻게 알 수 있을까요? 그리고 더 중요한 건 볼 수 없는 교란변수에 대해 어떻게 통제할 수 있을까요?\n",
    " \n",
    "![img](./data/img/fixed-effects/fe-graph.png)\n",
    " \n",
    "여기서 사용해볼 요령은 실험 대상을 확대하고 시간이 지남에 따라 대상이 어떻게 발전하는지 추적해서, 시간에 따라 고정된 모든 것을 이미 통제하고 있음을 확인하는 것입니다. 이 부분에서 시간에 대해 고정되어있고 측정되지 않은 교란 변수가 포함됩니다. 예를 들어 위의 그래프에서 시간에 따라 소득이 증가하는 것이 단순히 외모가 동일하게 유지되기 때문에(결국 시간이 고정되어 있기 때문에) 외모가 좋아진 것 때문이 아님을 이미 알 수 있습니다. 결론은 외모라는 변수를 제어할 수는 없고 측정할 수 없으므로, 패널 구조를 계속 사용할 수 있게되어 더 이상 문제가 되지 않는다는 것입니다.\n",
    " \n",
    "패널 구조에서 교란변수를 바라보는 또 다른 방법은 시간에 고정된 (time fixed) 교란변수를 각 실험대상에 고유한 속성으로 생각하는 것입니다. 이 부분은 인과 그래프에 중간 단위 노드를 추가하는 것과 같아요. 이제 실험 대상에 대한 통제가 어떻게 결과와 관찰되지 않은 시간 고정된 교란변수 사이의 backdoor path를 차단하는지 주목해주세요.\n",
    " \n",
    "![img](./data/img/fixed-effects/control-unit.png)\n",
    " \n",
    "잠깐 생각을 해보았을 때, 사람의 외모나 지성과 같은 속성을 측정할 수 없지만 그 특징을 가진 사람은 시간이 흘러도 같은 개인이라는 것을 알고 있습니다. 실제로 이러한 통제를 하는 메커니즘은 매우 간단해요. 여기서 여러분이 해야할 일은 특정 사람을 나타내는 더미 변수를 만들고 선형 모델에 추가하는 것입니다. 즉, 해당 더미변수를 추가함으로, 사람 자체를 제어할 수 있다는 것을 의미하고 있어요. 위에서 설명한 모델에서 사람에 대한 더미변수와 함께 결혼이 소득에 미치는 영향을 추정할 때, 회귀분석은 **사람에 대한 변수를 고정된 상태로 유지하면서** 결혼의 효과를 찾습니다. 그래서, 실험대상에 대한 더미변수를 추가하는 것을 고정효과 모델(fixed effect model)이라고 합니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed Effects\n",
    "\n",
    "문제를 좀 더 공식화하기 위해, 먼저 가지고 있는 데이터를 한 번 살펴보겠습니다. 앞서 설명드린 예시에 따라 결혼이 소득에 미치는 인과적 영향을 추정하려고 노력할 것입니다. 아래 데이터에는 여러 해 동안 여러 개인(`nr`)에 대한 `married (기혼여부)` 및 `lwage (소득의 로그변환)` 이 2가지 변수가 포함되어 있습니다. 소득은 로그로 변환한 형식입니다. 이 외에도 해당 연도의 근무 시간, 교육 년수 등과 같은 다른 제어할 수 있는 변수가 있습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nr</th>\n",
       "      <th>year</th>\n",
       "      <th>black</th>\n",
       "      <th>...</th>\n",
       "      <th>lwage</th>\n",
       "      <th>expersq</th>\n",
       "      <th>occupation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>1980</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.197540</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13</td>\n",
       "      <td>1981</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.853060</td>\n",
       "      <td>4</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>13</td>\n",
       "      <td>1982</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.344462</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>13</td>\n",
       "      <td>1983</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.433213</td>\n",
       "      <td>16</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>13</td>\n",
       "      <td>1984</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.568125</td>\n",
       "      <td>25</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   nr  year  black  ...     lwage  expersq  occupation\n",
       "0  13  1980      0  ...  1.197540        1           9\n",
       "1  13  1981      0  ...  1.853060        4           9\n",
       "2  13  1982      0  ...  1.344462        9           9\n",
       "3  13  1983      0  ...  1.433213       16           9\n",
       "4  13  1984      0  ...  1.568125       25           5\n",
       "\n",
       "[5 rows x 12 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = wage_panel.load()\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "일반적으로, 고정효과 모델은 아래와 같이 정의됩니다.\n",
    "\n",
    "$\n",
    "y_{it} = \\beta X_{it} + \\gamma U_i + e_{it}\n",
    "$\n",
    "\n",
    "여기서 $y_{it}$는 시간 $t$에서 각 실험대상 $i$의 결과이고, $X_{it}$는 시간 $t$에서 실험대상 $i$에 대한 변수들의 벡터입니다. $U_i$는 실험대상 $i$에 대해 관측 불가능한 변수의 집합입니다. 이러한 관측할 수 없는 변수들은 시간이 지남에 따라 변하지 않으므로(가정), 시간 첨자가 없음을 주목해주세요. 마지막으로 $e_{it}$는 오차항입니다. 위 예시에서 $y_{it}$는 임금에 대한 로그변환, $X_{it}$는 결혼 및 경력과 같이 시간에 따라 변하는 관찰 가능한 변수이고 $U_i$는 관찰되지 않지만 각 개인에 대해 일정한 변수입니다 (ex. 외모와 지성).\n",
    "\n",
    "고정효과 모델과 함께 패널 데이터를 사용하는 것이 실험대상에 대한 더미(dummy)를 추가하는 것만큼 간단하다고 설명드렸던 것을 기억하시나요? 맞는 말이지만, 실제로는 그렇게 하지 않아요. 예를 들어, 백만명의 고객이 있는 데이터 셋을 생각해봅시다. 각 고객에 대해 하나의 더미를 추가하면 백만 개의 열이 생성되고, 이런 결과가 나오는 것은 좋은 방법이 아닐거에요. 대신, 선형 회귀를 2개의 개별 모델로 분할하는 트릭을 사용할 거에요. 이전 챕터에서도 확인하셨겠지만, 지금이 회귀 모델에 대해 요약하기에 좋은 타임인 것 같아요. 변수들의 집한인 $X_1$과 $X_2$가 있는 선형 회귀 모델이 있다고 가정해봅시다.\n",
    "\n",
    "$\n",
    "\\hat{Y} = \\hat{\\beta_1} X_1 + \\hat{\\beta_2} X_2\n",
    "$\n",
    "\n",
    "여기서 $X_1$과 $X_2$는 변수 행렬(피쳐당 1개의 행, 관측값당 1개의 열)이고 $\\hat{\\beta_1}$ 및 $\\hat{\\beta_2}$는 행 벡터입니다. 아래와 같은 방법을 통해, 똑같은 $\\hat{\\beta_1}$ 매개변수 값을 얻을 수 있어요.\n",
    "\n",
    "1. 두 번째 변수 집합에서 $y$ 결과를 회귀 $\\hat{y^*} = \\hat{\\gamma_1} X_2$\n",
    "2. 두 번째 $\\hat{X_1} = \\hat{\\gamma_2} X_2$의 첫 번째 변수 집합을 회귀\n",
    "3. 잔차 $\\tilde{X}_1 = X_1 - \\hat{X_1}$ 및 $\\tilde{y}_1 = y_1 - \\hat{y^*}$를 구하기\n",
    "4. 변수 $\\hat{y} = \\hat{\\beta_1} \\tilde{X}_1$의 잔차에 대한 결과의 잔차 회귀\n",
    "\n",
    "마지막 회귀 모델의 매개변수는 모든 변수를 사용해 회귀를 하는 것 결과와 똑같습니다. 그런데, 이 부분이 어떤 도움이 될까요? 실험 대상에 대한 더미변수가 있는 모델의 추정을 2개로 나눌 수 있습니다. 먼저 더미를 사용하여 변수와 결과를 예측합니다. 위에서 설명한 1,2단계 입니다.\n",
    "\n",
    "이제 더미 변수를 가지고 회귀를 모델링을 하는 것이 해당 더미 변수의 평균을 추정하는 것만큼 간단하다는 것을 기억하시나요? 데이터를 가지고 이 부분이 맞는지 확인해 봅시다. 연도(year) 더미변수에 대한 함수로 임금을 예측하는 모델을 돌려 보겠습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    1.3935</td> <td>    0.022</td> <td>   63.462</td> <td> 0.000</td> <td>    1.350</td> <td>    1.437</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1981]</th> <td>    0.1194</td> <td>    0.031</td> <td>    3.845</td> <td> 0.000</td> <td>    0.059</td> <td>    0.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1982]</th> <td>    0.1782</td> <td>    0.031</td> <td>    5.738</td> <td> 0.000</td> <td>    0.117</td> <td>    0.239</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1983]</th> <td>    0.2258</td> <td>    0.031</td> <td>    7.271</td> <td> 0.000</td> <td>    0.165</td> <td>    0.287</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1984]</th> <td>    0.2968</td> <td>    0.031</td> <td>    9.558</td> <td> 0.000</td> <td>    0.236</td> <td>    0.358</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1985]</th> <td>    0.3459</td> <td>    0.031</td> <td>   11.140</td> <td> 0.000</td> <td>    0.285</td> <td>    0.407</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1986]</th> <td>    0.4062</td> <td>    0.031</td> <td>   13.082</td> <td> 0.000</td> <td>    0.345</td> <td>    0.467</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1987]</th> <td>    0.4730</td> <td>    0.031</td> <td>   15.232</td> <td> 0.000</td> <td>    0.412</td> <td>    0.534</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ C(year)\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "해당 회귀 모델이 1980년의 평균 소득을 1.3935로, 1981년에는 1.5129(1.3935+0.1194)과 같이 어떻게 예측하는지 봐주세요. 만약 연도별로 평균을 계산한다면, 회귀 모델과 똑같은 결과를 얻게됩니다. (기준 연도인 1980이 y절편임을 기억해주세요. 따라서 해당 연도의 평균 '임금'을 구하려면, 다른 연도의 매개변수에 y절편값을 추가해야 해요.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year\n",
       "1980    1.393477\n",
       "1981    1.512867\n",
       "1982    1.571667\n",
       "1983    1.619263\n",
       "1984    1.690295\n",
       "1985    1.739410\n",
       "1986    1.799719\n",
       "1987    1.866479\n",
       "Name: lwage, dtype: float64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"year\")[\"lwage\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "즉, 패널 데이터의 모든 사람에 대한 평균을 얻으면 본질적으로 다른 변수에 대한 개별 더미에 회귀를 해야해요. 이 부분은 다음과 같은 추정 절차에 대한 의미를 부여합니다.\n",
    "\n",
    "1. 개별 대상에 대한 평균을 빼는 방식으로 time-demeaned(시간에 대한 차분) 변수 생성:   \n",
    "$\\ddot{Y}_{it} = Y_{it} -  \\bar{Y}_i$  \n",
    "$\\ddot{X}_{it} = X_{it} -  \\bar{X}_i$\n",
    "\n",
    "2. $\\ddot{Y}_{it}$를 $\\ddot{X}_{it}$에 회귀\n",
    "\n",
    "위와 같은 단계를 거치면, 관측되지 않은 $U_i$ 변수가 사라집니다. $U_i$는 시간에 따라 일정하므로 $\\bar{U_i}=U_i$가 됩니다. 아래와 같이 2개의 식이 있다고 보면, \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_{it} & = \\beta X_{it} + \\gamma U_i + e_{it} \\\\\n",
    "\\bar{Y}_{i} & = \\beta \\bar{X}_{it} + \\gamma \\bar{U}_i + \\bar{e}_{it} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "그리고 우리는 1번식에서 2번식을 빼주게 된다면,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "(Y_{it} - \\bar{Y}_{i}) & = (\\beta X_{it} - \\beta \\bar{X}_{it}) + (\\gamma U_i - \\gamma U_i) + (e_{it}-\\bar{e}_{it}) \\\\\n",
    "(Y_{it} - \\bar{Y}_{i}) & = \\beta(X_{it} - \\bar{X}_{it}) + (e_{it}-\\bar{e}_{it}) \\\\\n",
    "\\ddot{Y}_{it} & = \\beta \\ddot{X}_{it} + \\ddot{e}_{it} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "시간이 지남에 따라 일정한 관측되지 않은 모든 변수들를 제거하는 효과를 가집니다. 정확히 말하면, 관측되지 않은 변수만 사라지는 것은 아니에요. 시간에 일정한 모든 변수가 해당이 됩니다. 이러한 이유로, 시간에 일정한 변수가 더미변수의 선형 조합이 되고 모델이 실행되지 않게되므로, 시간에 따라 일정한 변수를 포함할 수 없게됩니다.\n",
    "\n",
    "![img](./data/img/fixed-effects/demeaned.png)\n",
    "\n",
    "어떤 변수가 그런지 확인하기 위해, 개인별로 데이터를 그룹화 해서 표준 편차의 합계를 얻을 수 있어요. 만약 합계가 0이면 변수가 개별 대상에 대해 시간이 지남에 따라 변하지 않음을 의미합니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year            1334.971910\n",
       "black              0.000000\n",
       "exper           1334.971910\n",
       "hisp               0.000000\n",
       "hours         203098.215649\n",
       "married          140.372801\n",
       "educ               0.000000\n",
       "union            106.512445\n",
       "lwage            173.929670\n",
       "expersq        17608.242825\n",
       "occupation       739.222281\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"nr\").std().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "데이터를 보면, 개별 실험대상에 대해 일정하므로 인종과 관련된 더미변수인 `black`와 `hisp` 및 교육 변수를 제거할 필요가 있어요. 그리고 또한 직업 변수를 사용하지 않을 계획입니다. 해당 변수를 사용하지 않는 이유는 결혼이 임금에 미치는 영향을 중재할 것이기 때문입니다(미혼 남성이 더 많은 시간을 요구하는 포지션을 맡을 수 있기 때문입니다). 사용할 변수를 골랐다면, 이제는 해당 변수를 이용해 모델을 추정할 차례입니다.\n",
    "\n",
    "고정 효과 모델을 돌리기 위해, 먼저 평균 데이터를 가져와 볼게요. 실험 대상별로 모든 것을 그룹화하고 평균을 취해서, 평균 데이터를 만들 수 있어요."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.000</td>\n",
       "      <td>25.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2807.625</td>\n",
       "      <td>1.255652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2504.125</td>\n",
       "      <td>1.637786</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2350.500</td>\n",
       "      <td>2.034387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>0.125</td>\n",
       "      <td>35.5</td>\n",
       "      <td>0.250</td>\n",
       "      <td>2225.875</td>\n",
       "      <td>1.773664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>0.500</td>\n",
       "      <td>77.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2108.000</td>\n",
       "      <td>2.055129</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     married  expersq  union     hours     lwage\n",
       "nr                                              \n",
       "13     0.000     25.5  0.125  2807.625  1.255652\n",
       "17     0.000     61.5  0.000  2504.125  1.637786\n",
       "18     1.000     61.5  0.000  2350.500  2.034387\n",
       "45     0.125     35.5  0.250  2225.875  1.773664\n",
       "110    0.500     77.5  0.125  2108.000  2.055129"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = \"lwage\"\n",
    "T = \"married\"\n",
    "X = [T, \"expersq\", \"union\", \"hours\"]\n",
    "\n",
    "mean_data = data.groupby(\"nr\")[X+[Y]].mean()\n",
    "mean_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "데이터의 평균을 차분하기 위해서는 원본 데이터의 인덱스를 개별 식별자 'nr'로 설정해야해요. 그런 다음에 평균 데이터에서 기존 데이터를 단순히 빼면 됩니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-24.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>-135.625</td>\n",
       "      <td>-0.058112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-21.5</td>\n",
       "      <td>0.875</td>\n",
       "      <td>-487.625</td>\n",
       "      <td>0.597408</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-16.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>132.375</td>\n",
       "      <td>0.088810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-9.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>152.375</td>\n",
       "      <td>0.177561</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>263.375</td>\n",
       "      <td>0.312473</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    married  expersq  union    hours     lwage\n",
       "nr                                            \n",
       "13      0.0    -24.5 -0.125 -135.625 -0.058112\n",
       "13      0.0    -21.5  0.875 -487.625  0.597408\n",
       "13      0.0    -16.5 -0.125  132.375  0.088810\n",
       "13      0.0     -9.5 -0.125  152.375  0.177561\n",
       "13      0.0     -0.5 -0.125  263.375  0.312473"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demeaned_data = (data\n",
    "               .set_index(\"nr\") # set the index as the person indicator\n",
    "               [X+[Y]]\n",
    "               - mean_data) # subtract the mean data\n",
    "\n",
    "demeaned_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "마지막으로 평균을 뺀 데이터에 대해 고정 효과 모델을 돌릴 수 있게 됩니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>-6.679e-17</td> <td>    0.005</td> <td>-1.32e-14</td> <td> 1.000</td> <td>   -0.010</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1147</td> <td>    0.017</td> <td>    6.756</td> <td> 0.000</td> <td>    0.081</td> <td>    0.148</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0040</td> <td>    0.000</td> <td>   21.958</td> <td> 0.000</td> <td>    0.004</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.0784</td> <td>    0.018</td> <td>    4.261</td> <td> 0.000</td> <td>    0.042</td> <td>    0.115</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -8.46e-05</td> <td> 1.25e-05</td> <td>   -6.744</td> <td> 0.000</td> <td>   -0.000</td> <td>   -6e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(f\"{Y} ~ {'+'.join(X)}\", data=demeaned_data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "고정 효과가 생략된 모든 변수들의 편향을 제거한다고 생각한다면, 위 모델은 결혼이 남성의 임금을 11% 증가시킨다는 것을 의미하고 있어요. 해당 결과는 매우 중요한데요, 여기서 한 가지 고려할 세부 사항은 고정 효과 모델의 경우 표준 오차가 클러스터되어야(clustered) 한다는 것입니다. 따라서 모든 추정을 직접 수행하는 대신(교육학적인 이유로만 유용해요!) 관련 패키지인, `linearmodels`를 사용 및 해당 `cluster_entity` 인자를 True로 설정할 수 있어요."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>0.0040</td>    <td>0.0002</td>   <td>16.552</td>  <td>0.0000</td>   <td>0.0035</td>    <td>0.0044</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0784</td>    <td>0.0236</td>   <td>3.3225</td>  <td>0.0009</td>   <td>0.0322</td>    <td>0.1247</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.1147</td>    <td>0.0220</td>   <td>5.2213</td>  <td>0.0000</td>   <td>0.0716</td>    <td>0.1577</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>   <td>-8.46e-05</td> <td>2.22e-05</td>  <td>-3.8105</td> <td>0.0001</td>   <td>-0.0001</td> <td>-4.107e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels.panel import PanelOLS\n",
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "파라미터 추정치가 시간 차분 데이터에서 얻은 것과 어떻게 동일한지 확인해주세요. 유일한 차이는 표준 오류가 약간 더 크다는 것입니다. 이제, 데이터의 시간 부분을 고려하지 않은 간단한 OLS 모델과 비교해보겠습니다. 해당 OLS 모델의 경우, 시간에 일정한 변수들을 다시 추가해야 합니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    0.2654</td> <td>    0.065</td> <td>    4.103</td> <td> 0.000</td> <td>    0.139</td> <td>    0.392</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0032</td> <td>    0.000</td> <td>   15.750</td> <td> 0.000</td> <td>    0.003</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.1829</td> <td>    0.017</td> <td>   10.598</td> <td> 0.000</td> <td>    0.149</td> <td>    0.217</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1410</td> <td>    0.016</td> <td>    8.931</td> <td> 0.000</td> <td>    0.110</td> <td>    0.172</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -5.32e-05</td> <td> 1.34e-05</td> <td>   -3.978</td> <td> 0.000</td> <td>-7.94e-05</td> <td> -2.7e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>black</th>     <td>   -0.1347</td> <td>    0.024</td> <td>   -5.679</td> <td> 0.000</td> <td>   -0.181</td> <td>   -0.088</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hisp</th>      <td>    0.0132</td> <td>    0.021</td> <td>    0.632</td> <td> 0.528</td> <td>   -0.028</td> <td>    0.054</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>educ</th>      <td>    0.1057</td> <td>    0.005</td> <td>   22.550</td> <td> 0.000</td> <td>    0.097</td> <td>    0.115</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ expersq+union+married+hours+black+hisp+educ\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "위 모델은 결혼을 하면, 남성의 임금이 14% 증가한다고 말하고 있네요. 해당 효과는 고정 효과 모델에서 찾은 것보다 다소 더 큰 효과입니다. 지능과 아름다움 변수와 같은 고정된 개별 요인으로, 일부 생략된 변수의 편향이 모델에 반영되지 않음을 의미해요.\n",
    "\n",
    "## Visualizing Fixed Effects\n",
    "\n",
    "고정 효과 모델이 어떻게 동작하는지에 대한 직관을 키우고자 다른 예시를 조금 살펴볼게요. 만약, 여러분이 대기업 IT 회사에서 일하고 광고판 마케팅 캠페인이 인앱 구매에 미치는 영향을 추정한다고 가정해봅시다. 과거의 데이터를 보면 마케팅 부서가 구매 수준이 낮은 도시에 광고판을 배치하는 데 더 많은 비용을 지출하는 경향이 있음을 알 수 있어요. 그런데 여러분 이런 부분이 합리적이라고 생각하시나요? 매출이 급격히 치솟는다면 광고를 많이 할 필요가 없을 거에요. 해당 데이터를 가지고 회귀 모델을 돌린다면, 마케팅 비용이 높을수록 인앱 구매가 감소하는 것처럼 보이지만, 이는 마케팅 투자가 지출이 적은 지역에 편중되어 있기 때문입니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAE0CAYAAACirQ3aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABeVklEQVR4nO3dd1hT59vA8W+YLhRECYgKAu6Bs25x7wHu8VPrrNvWOuuso1bFXfdetU6cda+60CquaqsoihNFEWSpjLx/8JJKWQkGk+j9uS4vkzOec+eQ5M455z7PowgNDVUhhBBCGCETfQcghBBCZJQkMSGEEEZLkpgQQgijJUlMCCGE0ZIkJoQQwmhJEhNCCGG0JIkJvTt9+jTW1tZMnz5d36EkERgYiLW1Nf3799d3KEIDmzZtwtramk2bNhlEO+LTkCQmdC4+Pp7169fTvHlzChUqRJ48eXB1daVq1ar079+fHTt26DtEgxMWFsbs2bOpX78+zs7O2NnZUbx4cbp06cLevXtTXa9Zs2ZYW1tz+vRpjbZz7949vv32W8qXL4+9vT358uWjdOnStGzZkmnTpvHo0SON2kn84WFtbU3p0qWJj49PcbmoqCgKFiyoXvbOnTsatS+Epsz0HYD4vMTHx9OpUycOHTpEzpw5ady4Mfny5SMmJobbt2+zZ88eLl26RJs2bdTrVKhQgYsXL2Jra6vHyPXH19eXrl27EhwcjJubG23atCFXrlwEBgZy5MgR9u/fT/369Vm9ejU5c+bM8HZOnz5N+/btiY6OplKlStSrV4+cOXMSFBTEpUuXmDVrFsWKFaNAgQIat2lmZsajR484ceIE9erVSzbfx8eHN2/eYGZmRmxsbIZjFyI1ksSETu3YsYNDhw5RqlQp9u/fT65cuZLMf/fuHefPn08yLVu2bBQpUuRThmkw7ty5Q7t27YiIiGDq1KkMGDAAE5N/T5C8evWKXr16cfToUXr06MG2bduSzNfGt99+S3R0NIsXL6Zz587J5vv7+2NqaqpVmx4eHpw7d45169almMTWr19P3rx5cXFx4cKFCxmKW4i0yOlEoVOJX1SdO3dOlsAALC0tqV27dpJpqV0TSzxV9uDBA5YvX07lypVRKpWULl2a2bNno1Il9Ji2a9cu6tWrR758+XBzc2PEiBFER0cn23biqa/Q0FC+//57ihUrhlKppEqVKqxYsULdnibevn3LwoUL8fDwwNHRkXz58lG7dm1Wr16tVTsjR44kPDycwYMHM2jQoGQJytbWlo0bN1KwYEGOHTvGzp07NW77Qy9fvuTevXvkzJkzxQQGULhwYVxcXLRqN1euXLRs2ZIDBw4QHBycZN4///zDhQsX6Ny5M2Zmqf9eTjxCLFSoEHZ2dri7uzN69GhevnyZ4vIBAQF0794dJycn8uXLR8OGDTl06FCacT5//pzRo0dTvnx5lEolTk5OeHl5cerUKa1erzA8ksSETllbWwMJ1150Zfz48cycOZOKFSvSrVs33r9/z5QpU5g5cyaLFy9m0KBBuLq60r17d3LmzMmKFSsYO3Zsim3FxMTg5eXFyZMnadu2Ld26dSMkJIQRI0YwevRojeIJDw+nWbNmjB8/HpVKRefOnenSpQtv3rxh2LBhDBgwQKN2Hjx4wMmTJ7G0tOS7775LdbkcOXIwePBgANasWaNR2/9lZWWFmZkZkZGRBAUFZaiN1HTv3p2YmBh+/fXXJNPXrVsHQLdu3VJdd/369bRs2ZIzZ87QuHFjBg4ciKOjI0uXLqV27do8efIkyfL37t2jfv367N69m0qVKtGvXz8cHR3TvHZ48+ZNatasybJlyyhUqBB9+vShWbNm+Pn54enpyYYNGz5yDwh9ktOJQqdatGjB/PnzWb16tfrL3t3dHScnJxQKRYba/Ouvvzh37hx2dnYADBo0iIoVK7JgwQKyZcvGqVOncHV1BVD/2t64cSNjxowhb968SdoKCgrC2dmZc+fOYWlpqV6nTp06LFu2jNatW1O5cuU04/nhhx+4fPkykyZN4ttvv1VPf/fuHV27dmXz5s20bNmSJk2apNmOr68vAGXLlsXGxibNZevUqQPApUuXiIuL0/q0n6WlJU2bNmXPnj00atSInj17UrVqVUqWLEn27Nm1auu/qlWrRpEiRVi/fj1Dhw4FEvbFli1bqFGjhvpv81+PHz9m+PDhZMuWjaNHj1K8eHH1vKlTp+Lt7c2wYcPYsmWLevrw4cMJCQlhypQp6sQOcPDgQTp27JhsG3FxcXTv3p2wsDD27t1LjRo11POCgoKoV68eI0eOpHHjxsneK8I4yJGY0Cl3d3eWL1+OnZ0dW7dupXv37pQtWxYnJyfat2/Pjh07Uq1kS83w4cPVCQzAycmJqlWrEhkZSc+ePZN8SebKlYvGjRvz/v17bt++nWJ7EyZMUCcwSDhlN2zYMIB0y6pfv37N5s2bKVOmTJIEBgmJYsKECQBJvnhT8/z5cwAcHR3TXTZxmXfv3hESEpLu8ilZsGABLVq04OHDh0ycOJGGDRuSP39+qlSpwrhx43j48GGG2gXo2rUr9+7dU1dJ7t27l5CQELp3757qOr/99hvv37+nV69eSRIYJPzNHRwcOHToEM+ePQPgyZMnnDhxgvz58ye77aFx48bUrFkz2TYOHz7M3bt36dWrV5IEBmBvb8/gwYOJjo5m9+7dGXrdQv/kSEzonJeXF82bN+f06dOcP3+e69ev4+vry+HDhzl8+DC//vormzdvxsLCQqP2ypQpk2yavb19qvMcHBwAePr0abJ5ZmZmKR5pVa9eHYDr16+nGcvly5eJjY3FxMQkxfvaEivwDLGU3Nramg0bNhAYGMjx48e5cuUKfn5+3Lx5k3/++YfVq1ezbt06GjRooHXbnTt3ZsqUKaxbt46aNWuydu1abGxsaNmyZarrXLt2DYBatWolm5clSxaqVKmCj48P169fx8HBQf23qVKlSorX2KpXr57sVoPEa7SPHz9O8e8VEBAAkOoPHmH4JImJTGFubk7dunWpW7cukHBaZ8+ePQwaNIhjx46xatUqjW8iTqmsPPF0WlrzYmJiks2ztbVN8VRc4qmkN2/epBlL4lHQ1atXuXr1aqrLRUREpNkOoD66/O91n5QkLmNhYUHu3LnTXT4tTk5O9OjRgx49eqjb/v777zl48CADBgzg1q1bmJuba9Wmra0tzZo1Y+/evfz555+cOXOGfv36JTni/a/Eff3hUfaHlEolkHAP3YfLp3baL6V2Ev9ee/bsYc+ePanGEhkZmeo8YdjkdKL4JExNTfHy8lInLn1Vhb169Yq4uLhk0xMr69K7Dytxft++fQkNDU31X3pHdABVq1YFEhJiaGhomsuePHkSgEqVKml9PSw9jo6OrF69GgsLC4KDg7l582aG2vn666959+4dX3/9tfp5WhL35YsXL1Kcn3i6NXG5xP//WwWZKKV2EtdZv359mn+vxYsXp/3ihMGSJCY+KSsrKwCtytB1KTY2NsX7lc6ePQukfHryQxUrVsTExCTZvW4Z4ezsTK1atXj37h3z5s1LdbnIyEh++eUXAPXRk65ZWlpqfHo3NbVq1cLZ2ZknT55QuXJlihUrluby7u7uACn2NvLu3Tv13ylxucS/zYULF1K8cTrxb/ihSpUqAejk7yUMkyQxoVPbt2/nxIkTKRZvPH/+XF12nXgNSh+mTJnCu3fv1M9fvXrFnDlzAOjSpUua6+bJk4cOHTpw48YNpk+fnuKX6ZMnTzS+JjZjxgyyZ8/OggULWLJkSbLkHhISQteuXXnw4AH16tWjdevWGrX7X5GRkcycOTPVo57FixcTERGBtbV1siILTSkUCtavX8/GjRtZsGBBusu3b98eCwsLVq1alWx/zZkzh6dPn9KwYUP1NU5HR0fq1KnDo0ePWLJkSZLlDx48mGIybNq0KS4uLqxZs4bff/89xTiuXbuW4WIZoX9yTUzo1KVLl1i6dKn6JmInJycgoTPdw4cPEx0dzVdffUWfPn30Ep+9vT3R0dFUq1aNJk2a8O7dO/bs2cPz58/55ptv0i2vB5g5cyYBAQHMmDGDLVu2UK1aNZRKJc+fP+fu3bv8+eefTJs2TaNeSIoXL87WrVvp1q0bY8aMYfXq1dSuXRsrKysePXrEoUOHePPmjbrbqdR665g3b16y+7QSdevWjeLFi/PTTz8xY8YMKlSoQOnSpbG2tub169f4+vpy69YtzMzMmD9/fprXsdJTpkyZdI9mExUsWJAZM2YwbNgw6tSpg6enJ0qlkgsXLnD27FkcHR2ZPXt2knW8vb1p0KAB48eP59SpU5QpU4b79++zd+9eGjduzMGDB5Msb25uzsaNG2ndujWdO3emYsWKuLu7kz17dp48ecL169fx9/fnjz/++OhrjUI/JIkJnRo8eDCFCxfmxIkT3Lp1ixMnThAVFYWNjQ1fffUVnp6e/O9//9O6cEBXzM3N2bVrF5MnT2b79u2EhIRQqFAhvv/+e40Tq5WVFfv27WPDhg1s27aNffv28fbtW/LmzYuTkxMTJ07Ey8tL45iqV6/O5cuXWbFiBQcPHmTLli1ER0eTO3duqlevTqdOnWjRokWa99kdO3Ys1Xk1atSgcuXK7NixgxMnTuDr66vuYcPc3JwCBQrQo0cP+vbtm+GjsIzq0aMHLi4uLFy4kP379xMZGYmDgwN9+/ZNdmsFgKurK0ePHmXSpEmcPHmSc+fOUbJkSTZt2sTLly+TJTGAEiVKcPbsWZYsWcLvv//O5s2bUalUKJVKihUrpn7PCuOkCA0N1c/FCSE+MWtrawoUKMCNGzf0HYoQQkfkmpgQQgijJUlMCCGE0ZIkJoQQwmhJYYf4YqR3Q7EQwvjIkZgQQgijJUlMCCGE0ZIkJoQQwmhJEjMC/v7++g7BqMj+0o7sL+3I/tJOZu8vSWJCCCGMliQxIYQQRkuSmBBCCKMlSUwIIYTRkpudhRBAwoChkZGR+g7D4GXJkoWwsDB9h2E0NNlf2bNnx8wsY+lIkpj4bASGxzDVL5yAl5a4PA1hXHkrnKz0M+SLsYmNjSU8PBxra+s0h3wRCaNgZ8mSRd9hGI309pdKpSI0NBQrK6sMJTJJYuKzEBgeg+ehV9wPjwNMufwmmkvB79nVyFYSmQYiIyMlgQm9UCgUWFtb8+bNG3LlyqX1+nJNTHwWpvqF/38C+9f98Dim+oXrKSLjIwlM6MvHvPc0PhJTqVTcvHmT27dv8+rVKxQKBba2thQpUoSSJUvKB0Do1bOouBSnB6UyXQjxeUg3if3xxx9s2rSJAwcOEBERgUqVdCBohUJBjhw5aNy4MV26dMHDwyPTghUiNQ7ZTFOcbp/KdCHE5yHVJHb06FGmTZvG1atXKV68OF27dqVs2bI4OztjbW2tvhgXGBjI1atXOXHiBJ6enri7uzN+/Hjq1av3KV+H+MKNK2/FpeD3SU4pFrIyZVx5Kz1GJUTGBQYG4u7uzokTJyhXrpze4pg+fTp79uzh/PnzeoshLYrQ0FBVSjPy5ctH165d6dmzJ0WLFtWosdu3b7Nq1So2bdrEkydPdBrol8zf35/ChQvrOwyDp65OfBWBi20OqU7UkL+/P3Z2dhm6qK5v/fv3Z/PmzQCYmpri4OBAw4YNmTBhAtbW1pmyzbdv336S6sS4uDhevnyJra1thsvPNXH69GlatGjBvXv3sLW1TTY/IiKC9+/fkzt37gy1r+n+CgsLy9B7MNU9c+PGjRRfUFqKFi3KzJkzGTVqlNaBCPGxnKzMWeGRG3//VxQuXFDf4YhPpHbt2ixbtozY2Fhu377NoEGDCAsLY9WqVZm63ZiYGMzNM+9HkqmpKUqlMtPa11SOHDn0HUKaUq1O1DaB6WpdIYTQhqWlJUqlEkdHR+rWrYuXlxfHjx9PsszGjRupXLkySqWSChUqsGjRIuLj49Xz7969S9OmTVEqlVSsWJHDhw/j6OjIpk2bgIRTe9bW1mzfvp02bdpgb2/PmjVrNGp7zZo1VKhQAaVSiYuLC61btyY2NhaAmzdv0rJlSwoUKICjoyPVq1fnjz/+SLLNK1euqNs6e/Ys9erVQ6lUUrhwYcaMGcP79+/V85s1a8b333/P5MmTcXFxwc3NjXHjxiWJR1vTp0+natWq6uf9+/enQ4cOLFmyhOLFi+Pk5MSAAQOIiopSL6NSqZg/f776ElS1atXYsmVLhmNIi1bHqGFhYQQFBREdHU3WrFmxt7c3ylMQQgjNzHfy/qTbGxo4/KPWf/DgAceOHUtyhLRu3Tp++uknZs6cibu7O3///TdDhw7F3Nycvn37Eh8fz//+9z/s7Ow4cuQIb9++ZcyYMbx79y5Z+z/++CMTJkxg0aJFmJubp9v2lStXGD58OEuWLKFKlSqEhYWpkxRAnz59KFWqFMeOHcPMzIybN2+meurt6dOntGvXjg4dOrB48WLu37/PkCFDMDExYdq0aerltm3bxjfffMPhw4e5ceMGvXv3pmzZsrRt2/aj9u2Hzp8/j1KpZNeuXTx58oSvv/4aNzc3hg0bBsDUqVPZvXs33t7eFChQgOvXrzN06FCsra1p1KiRzuIADZJYSEgICxYswMfHh0ePHiWbnz9/fry8vBg8eDB58uTRaXDCMCVee3oWFYdDNlO59iT06ujRozg6OhIXF8fbt28Bknypz5o1ix9//JFWrVoB4OzszP3791m1ahV9+/blxIkT+Pv7s3PnTvLlywfATz/9lOKXbd++fWnRooU60aTX9qNHj8iePTtNmjTByiqhyKh06dLq9h49esSgQYMoUqQIAC4uLqm+zlWrVmFvb8/s2bMxMTGhaNGiTJw4ke+++46xY8eSLVs2IOGyztixYwFwc3Nj3bp1nDp1SqdJzMrKirlz52JqakrRokXx9PTk1KlTDBs2jMjISBYtWsTOnTupVq0ab9++pWjRoly+fJmVK1d+2iR2//59mjVrxosXL/Dw8KB9+/bY29tjaWnJu3fvCAoK4vLlyyxatIht27axb9++NP8Iwvgl7RkjgfSMIfSpWrVqzJ8/n+joaNatW8eDBw/o168fAC9fvuTx48d89913fP/99+p1YmNj1bcL3blzBwcHB3UCAyhfvjwmJsmvtnxYJahJ23Xq1CF//vy4u7tTr1496tSpQ4sWLdQJbcCAAQwZMoTNmzfj4eFBy5Yt1Qntv27fvk3FihWTxFW1alXev39PQEAApUqVAqBkyZJJ1rO3tyc4OFiDPam5okWLYmr67+0r9vb2XLp0SR3n27dvadu2LQqFApVKhUKhICYmhoIFdX+tOs0kNmbMGExMTDh//nya1XH+/v54eXkxduxYdaWQ+Dyl1TPGCo+MVS8J8TGyZcum/vE8c+ZMmjdvzsyZMxkzZoz6WtCcOXOoXLnyR28re/bs6seatG1lZcUff/zB2bNnOXnyJHPnzmXKlCkcP34cBwcHxowZQ/v27Tly5AjHjx9nxowZzJkzh65du2oV14edTfy32CQxkehSWttI3C+bN2+mQIECvHv3DktLS4BMqbJMs8WzZ8/yww8/pFveXbhwYQYMGMDPP/+s0+CE4ZGeMb4sH3uNSh9GjRpFu3bt+Prrr3FwcMDBwYH79+/TqVOnFJcvUqQIz54949mzZzg4OABw5cqVdIsh7Ozs0m0bEr64PTw88PDwYMyYMbi5uXHo0CG+/vprAFxdXXF1daVfv34MGzaMDRs2pJjEihYtio+PD/Hx8eqjsfPnz2NhYUGhQoU02TWfRNGiRbG0tOTRo0d4eHhk+i0JaSYxU1NTdRVNeuLi4qTrqS+A9IwhDF3NmjUpWrQo3t7ezJ49mzFjxjBy5Ehy5cpFw4YNiYmJ4dq1azx79oxhw4ZRp04dChcuTP/+/ZkyZQpv375l7NixmJmZpfudll7bBw8e5P79+1SrVg0bGxtOnz5NREQERYoUITo6mvHjx9OqVSsKFixIcHAwvr6+VKhQIcVt9erViyVLlvD999/Tr18/Hjx4wI8//kifPn3U18M+xq1bt5IV6iWeotSGlZUVgwcPZvz48ahUKipWrEhMTAyXLl3CxMREnbx1Jc0kVqdOHX755ReqV69O+fLlU13Oz8+PhQsXUrduXZ0GJwyP9IwhjMGgQYMYOHAgQ4cOpVu3bmTLlo0FCxYwefJksmTJQvHixenTpw8AJiYmbNy4kcGDB1OvXj0KFizI1KlT6dq1a7pHEOm1nStXLvbv38/MmTOJjo6mUKFCLFiwgGrVqvH+/XtCQ0MZMGAAz58/J3fu3DRq1IgpU6akuK18+fKxbds2JkyYQM2aNcmVKxdt27ZlwoQJOtlnLVq0SDbt8ePHGWpr7Nix5M2bl19++YX79+9jZWVF6dKlGTp06MeGmUyqPXYAPHv2jObNm3P//n3KlClDuXLlUCqV6sKO58+fc+XKFa5fv46zszP79+9XH44L3TG0HjsSqxODouKwN8DqREPbX4bOmHvsyCw3btygZs2anDx5krJlyyaZ96l67Phc6K3HDgAHBwf++OMPVq5cyZ49e9i8eXOSeycsLS0pUaIEEydOpFevXgZ/Z7fQjcSeMYT4XOzdu5fs2bPj4uLCw4cPGTt2LKVKlcLd3V3foYl0pFsqkj17doYOHcrQoUNRqVS8fv1afbOzjY2NXAcTQhi9iIgIJk2axJMnT7C2tqZGjRr89NNP8v1mBLSqd1QoFBnuBFIIIQxVp06d0qwwFIZLoyQWHx/P+fPnuX79erJup0qXLk3VqlWT3PimqbNnz7Jw4UJ1Nc+iRYvo0qVList+++23rF27lilTpjB48GCttyWEEOLzk24S2759OxMmTCAoKCjFG+YUCgVKpZIff/yR9u3ba7XxyMhISpQoQadOndR32Kdk9+7dXL58+YsrGlEPLfLSEpenIQZXQCGEEPqWai/2AFu3bqVPnz44OzuzYsUKrly5wtOnTwkJCeHp06dcuXKFZcuW4ezsTL9+/di6datWG08c96dVq1YpdvEC8PDhQ0aPHs3KlSszdUwdQ5PYvdO2gGguvzFlW0A0nodeERgeo+/QxGdK1706CKGpj3nvpZnE5s2bR6NGjfj9999p06YNzs7OZM2aFYVCQdasWXF2dqZdu3YcOHCABg0aMG/evAwHkpLY2Fh69+7N8OHDNR6Y83ORVvdOQuha9uzZCQ0NlUQmPjmVSkVoaGiSLr20keahTUBAQJqn+T7UrFkzRo4cmaEgUjN9+nRy585Nr169NF7H399fpzHoS8BLSyD5dcaAVxH4+7/69AEZmc/lffCp3L9/H4A3b96kelZEiMwQHx/P+/fvefnyZYrz07vnM80kZm9vz5UrV+jWrVu6gfj5+WFvb5/ucpo6ffo0v/76K6dPn9Zqvc/lJleXpyFcfhOdfLptDhm1OB1ys7N2ZH9pR/aXdjJ7f6X5k6tHjx6sXbuWESNGcOvWrRSXuXXrFsOHD2f9+vX06NFDZ4GdOXOGoKAgihYtiq2tLba2tjx69IiJEydSokQJnW3HUI0rb0Uhq6RHYtK9kxBCJJXmkdiQIUMIDQ1l0aJFrFq1CktLS/Lmzavudio4OJh3795hZmbG4MGDddovVu/evdUDzSVq06YNbdq0oXv37jrbjqFysjJnVyPbhOrEVxG42OaQ6kQhhPiPNJOYQqFg4sSJ9O3bl/3793Pjxg2CgoLUfWHVqVOH0qVL07Rp0yQDymkqIiKCgIAAIOG86OPHj7l+/To2NjYUKFCAvHnzJg3WzAylUvnFHMondu/k7/9KTiEKIUQKNKpZd3BwoHfv3jrf+JUrV5L0nDx9+nSmT59Op06dWLJkic63J4QQ4vOi1xuvatasSWhoqMbL37hxI/OCEUIIYXQ0SmJHjx7Fx8dH3e1U4unExG6nvLy8aNCgQWbHKoQQQiSRZhKLjIyka9eunDx5kuzZs1OqVClq1KihLuwICgpi7969/Pbbb3h4eLBhwwYZjkUIIcQnk2YSmzRpEufPn2fevHl06tQJc/PklXExMTFs3ryZUaNG8eOPPzJr1qxMC1YIIYT4UJr3ie3evZtBgwbRrVu3FBMYgLm5Od26dWPgwIHs2rUrM2IUQgghUpRmEouIiNC4dN7R0ZHIyEidBCWEEEJoIs0k5u7uztq1a4mIiEizkYiICNasWSNDeQshhPik0r0m5unpScWKFWnXrh1ly5ZFqVSqCzueP3/OlStX2L59O2FhYXI6UQghxCeVZhKrXLkyhw4dYsqUKSxdupTY2FgUCoV6vkqlwszMjNq1azNu3Dg5EhNCCPFJpXufWJkyZdi2bRvh4eHcunUr2X1ixYsXJ2fOnJ8iViGEECIJjXvssLKyonLlypkZixBCCKEVGf1OCCGE0dJZEvP29sbW1lZXzX1ygeEx9DkVQvMDwfQ5FUJgeIy+QxJCCJEOnXYArFKpdNncJxMYHoPnoVfcD49TT7sU/J5djWxl/C4hhDBgaSYxHx8fjRtKbeRnYzDVLzxJAgO4Hx7HVL9wVnjk1lNUQggh0pNmEuvZsycKhULjI6wPy++NybOouBSnB6UyXQghhGFIM4lZWVlRpkwZRo8enW5Dv/32G7/++qvOAvuUHLKZpjjdPpXpQgghDEOaSaxcuXI8fvyYGjVqpNuQr6+vzoL61MaVt+JS8PskpxQLWZkyrryVHqMSQgiRnjSrE8uVK0dAQIBGoy+rVCqjLexwsjJnVyNb2rlkpaa9Be1cskpRhxBCGIE0j8T69+9P/fr1Ux2G5UMjRoxgxIgROgvsU3OyMpciDiGEMDJpJjGlUolSqfxUsQghhBBakR47hBBCGC1JYkIIIYyWJDEhhBBGS5KYEEIIoyVJTAghhNGSJCaEEMJoZagX+4iICEJDQ1O8ublAgQIfHZQQQgihCY2T2Nu3b5kxYwYbNmwgJCQk1eXSmieEEELoksZJ7Pvvv2fz5s00a9aMqlWrYm1tnYlhCSGEEOnTOInt3buXbt26MW/evEwMRwghhNCcxoUdCoUCd3f3zIxFCCGE0IrGSaxp06acPHlSpxs/e/YsHTt2pHjx4lhbW7Np0yb1vJiYGCZOnEi1atXIly8fRYsWpXfv3jx69EinMQghhDBeGiex77//nvv37zNkyBAuXbpEUFAQwcHByf5pIzIykhIlSvDzzz+TNWvWJPOioqK4du0aw4cP59SpU/z66688efKEtm3bEhsbq9V2NPHs8lPOzTrDVq9feXHjuc7bF1+uwPAY+pwKofmBYPqcCiEwPEbfIQnx2VCEhoZqNAiYjY3NvyspFKkul9HqREdHR2bOnEmXLl1SXeaff/6hSpUqnD17lpIlS2ZoOyl5eTuYX5usRxWXdFdY5spCw9mNcWngprNtZYS/vz+FCxfWawzGxJD2V2B4DJ6HXiUbcNWQxqszpP1lDGR/aSez95fGhR0jR45MM3l9CuHh4QA6r4x8fjUoWQIDeBf2lr29d6mf15pQB/fu5TAxk3vEhWam+oUnSWAA98PjmOoXLuPXCaEDGh+JZbb0jsTev39PixYtsLGx4bfffku1HX9/f623/fZlNGf6neJ96DuN13Fu7UKRr4thls0wfk0Lw9TvuiWX35gmm14hVxxLS2v+fhPiS5XeUVyGeuz41GJjY+nbty9hYWFs3rw5zWUzdNhaGFyOuLLOYyUxUZpdr3iwM4AHOwMAKFTfldqT65HTMaf229aAnL7QjiHtL5enIVx+E518um0OChcuqIeIkjOk/WUMZH9pR2+nExOTRceOHVEoFOkmj0SdOnXSTWT/LzY2ll69enHr1i327dtH7tyZcwomu112Bvw9FIBnfk85NuoQr+680mjd+0fvcf/oPQBsi+Wh/sxG2Ls7ZEqcwriMK2/FpeD3ya6JjStvpceohPh8pJrEBgwYgEKhoE2bNlhYWDBgwIB0G1MoFDpNYjExMfTs2ZO///6bffv2oVQqddZ2WhzK5+N/R3oA8OZRGCcmHOPB8QCN1n31z0u2tEy4VcA8mzkNvBvj1rSI3q8nCv1wsjJnVyNbpvqFExQVh322hARmKEUdQhi7VJPYtWvXALCwsEjyXJciIiIICEhIDvHx8Tx+/Jjr169jY2ODg4MD3bt358qVK2zevBmFQsHz5wml7zlz5kxWkp9ZchbIRas1rQF4F/6O895nuLb2ikbrxkTF8PuAvernNcbUomyvCpiaJ79GIj5fTlbmUsQhRCbRa2HH6dOnadGiRbLpnTp1YvTo0an2ELJo0aI0S/E/hfjYeK6u8eP01JMZWr/0/9ypPrImlrmypLusnIPXjuwv7cj+0o7sL+0YTIl9ZqhZsyahoaGpzk9rnr6ZmJlQvk9FyvepCMC9Q/4cHn6Q9280qzi7sfEaNzYmHN06eThTZ2p9chW0zqxwhRDis5TqDU/9+vXjzp07Wjd4584d+vXr91FBGSPXRoXpf2MwQwOH03Hv/8hTIq/G6waeesDamiuZ7+TN+rqreXrpSSZGKoQQn49Uj8Rev35N1apVqVq1Kq1bt6Z27dq4uLikuGxAQAAnTpzAx8cHX19fGjRokGkBGwNlGXu6HOgOQPjTN5yadIJ7hzS7f+31vRC2tUmoBDW1NKXBrMYoiss1NCGESEma18QuXLjAggULOHz4MHFxcVhZWeHk5IS1tTUqlYrQ0FAePnxIeHg4ZmZmNG7cmMGDB1OpUqVP+RqMxvuI9/jOOcuVVZcztH7VETWo0LcSphaS1NIi1yy0I/tLO7K/tJPZ+0ujwo6XL19y8OBB/vzzT+7cucPr168ByJ07N0WLFqVy5co0aNAAW1vbTAv0cxMfF8+NDVc5OfF4htYv2bE0NcbUIov1p6nSNCbyJaMd2V/akf2lHYNIYiLzBRy7x9HhB4kOSd67Q3oKVC9InWkNsClkk/7CXwD5ktGO7C/tyP7SzmddnSj+5VLPlb5XBgIQfPMFx8Yc5vm1II3WfXT2IetrrwIgl5M1DWY1wrFygUyLVQghDIUkMQOUt6QdHff8D4CI5xH8PmIPz0491WjdsMBQtrffkvBEAQ1nN6FY6xLSY4gQ4rMkSczA5VDmoNz4SrRfX5iYqPdcWODL5SUXNVtZBYeHHeDwsAMAVBlWjQr9vsLMUv7sQojPg3ybGRHzbBbUGF2LGqNroYpX8dfm6xz/4YjG6/vOOYfvnHMAFG9bkppjPciaO1tmhSuEEJlOkpiRUpgoKN3FndJdErrmenDyPkdHHCTyRaRG6/+9/SZ/b78JgGOV/NT7qSE2rtK/nxDCuEgS+0w41y5E7z/7A/DydjDHxxzh2WXNrqM98X3M+rqrAbBytKLBrCYUqG4YY10JIURatE5ib9684fLlywQHB1O7dm3s7OwyIy7xEfIUzUv7nZ0BiAyO5PTUk9ze9bdG64Y/CWdn563q5/VnNqJE+1JSGCK+eIHhMUz1CyfgpSUuT0NkSB0DodV9YrNnz2bOnDlERUWhUCjw8fHBw8ODV69eUapUKaZNm0bPnj0zM94vkq7us4h9G8Ofv1zg4kLfDK3/1eAqVBpUBbMshn0AL/fxaEf2V/oCw2PwPPQq2eCmuxrZSiJLR2a/v1LtAPi/Vq9ezdSpU2nbti1r1qxBpfo399na2tK0aVN27dqVGTEKHTHLYk7V4TUYGjicIfe/p/7MRlqtf3GhL4uKzmO+kzcHh+4n6qVm19+EMHZT/cKTJDCA++FxTPUL11NEIpHGP6mXLVuGp6cn8+fPJyQkJNn8MmXKsGTJEp0GJzKPwkRByQ6lKdmhNAAPzwRydORBwp9o9qG8vevvJKcomy5uQeFmRTMlViH07VlUXIrTg1KZLj4djZPYgwcP6N+/f6rzra2t1X0qCuNTsIYTPc99A0DI3VccH3uEJ76PNV4/YQTrhFGsy/etSM2xtTMhSiH0wyFbyp1u26cyXXw6Gicxa2trgoODU53/999/o1QqdRKU0K/cbra03dIRgOiQKE7/dIq/t93UeH2/5ZfwW34JALtSSjrs7oKJmcZnroUwOOPKW3Ep+H2ya2LjylvpMSoBWlwTa9iwIevWrUvxaOuvv/5i/fr1NG3aVKfBCf3LmjsbDb2bMDRwOAPvfEvl76pptf6Lv56z0HUO8528me/kTXRIVCZFKkTmcbIyZ1cjW9q5ZKVCrjjauWSVog4DoXF14vPnz6lXrx6xsbE0bNiQjRs30rZtW2JjY9m/fz/58uXj2LFj5M4tN8zqmiFWj6niVezqvoOHfzzIcBud9nfFrpTuj94NcX8ZMtlf2pH9pR2D6cVeqVRy8uRJpkyZwp49e1CpVGzbtg0rKyvatWvHpEmTJIF9QRQmCrw2tFU/v7jgPOdnn9Wqjc3NNqgfN5rblGKtS+gsPiHElyHD44m9fPmS+Ph48uTJg4mJXO/ITMb2y+/Byfvs7r4jw+uX6VqWOlPrZ3h9Y9tf+ib7Szuyv7RjMEdi/5UnTx4AgoKCCAsLo2hRKa8WCZxrF2Jo4HAAQgNDWVdrpVbrX99wlesbrgJg45qbLoe6Y2ouVWBCiOQ0TmJr167l4sWLLF68WD1txIgRrFqVMBhjqVKl8PHxwdbWVvdRCqNl7WStTmgxUe9ZWWkp7yPea7z+63sh/OI2V/2896X+ZM+bXedxCiGMk8ZJbNWqVVSuXFn9/PTp06xcuZJ27dpRokQJvL298fb2Zvr06ZkSqDB+5tks6H9zCAAqlYp9fXcTcPiuVm2srPjvDfUddnXBvpyDTmMUQhgXjZNYYGAg3bt3Vz/38fHB0dGRpUuXYmJiQlhYGD4+PpLEhEYUCgUtVniqn/st/5PT005p1cYWz03qx/V+bkipTmV0FZ4QwkhonMTi4uIwN//3nogTJ05Qv359dVGHi4sLQUFBuo9QfBHK961E+b6VgIQusHy6bNNq/WOjD3Ns9GEAHBsVoPByufAuxJdA47JCJycnTp1K+KV85coVHjx4QN26ddXzX7x4gZWV3L0uPl7BGk4MDRzO0MDh9DjbR+v1nxx6pL65er6TN7FvYzMhSiGEIdD4SKxnz56MGDGCf/75h6dPn+Lo6EjDhg3V8319fSlWrFimBCm+XDnz51IXhsS+jWF19RVEv9Su149FReepH399uje5ClrrMEIhhD5pnMR69+6NhYUFhw8fpmzZsnz77bdkyZIFgNevXxMcHCxjiYlMZZbFnL6XBwAJhSGHhv7O7d2aDfaZaG3Nf8v9W672olA9V53GKIT4tDJ8s7P4dOTmyvRdW3eFkxOOZXj9igO+ovqoWjqMyHjI+0s7sr+0Y7A3OwthSNy7l8O9ezkAzm48w6Wx2o1efWnxRS4tvgiAbRFbuhz+GoVCofM4hRC6pVUSe/HiBRs2bODq1au8efOG+Pj4JPMVCgV79uzRuL2zZ8+ycOFCrl27xrNnz1i0aBFdunRRz1epVPz888+sW7eO0NBQKlSogLe3N8WLF9cmbPGFsausVF9HC3sYmuQUoiZe3XnFAufZ6ucD/h6CeTYLncYohNANjZPYrVu3aN68OVFRUbi5uXHr1i2KFStGaGgoz549o1ChQjg6Omq18cjISEqUKEGnTp3o169fsvnz589n0aJFLFq0iMKFCzNz5ky8vLz4888/pRJSaCRXQesPCkNikxR5aGpx8QXqx92O98TGVTq6FsJQaJzEfvzxR7JkycKJEyfIkSMHbm5uTJ8+HQ8PD7Zv387IkSNZvXq1Vhtv2LChusJxwIABSeapVCqWLFnCt99+S6tWrQBYsmQJhQsXZvv27fTo0UOrbQlhlsVMndAAtrb+lWeXn2rVxvq6/77Hmy5uQeFm0meoEPqkcRLz9fVl4MCBODk5qQfGVKkSakLatm2Lr68v48ePZ+/evToJLDAwkOfPnye5Fy1r1qxUq1aNCxcuSBITH639zs7qxxcX+nLe+4xW6/8+YC+Q8H4v26M8HpPqpr2CEELnNE5iMTEx2NvbA6hL68PCwtTzS5cuzW+//aazwJ4/fw5A3rx5k0zPmzcvz54909l2hAD4anAVvhpcBchYjyFX1/hxdY0fAFb5c9LjTB8pDBHiE9A4iRUoUIDHjx8DCUdE9vb2XLx4UX2q79atW2TPrv/exf39/fUdQqb4XF9XZvmo/aWEpkcT3tfRwdGc6HRYq9XDH79JUhjSYHdTzLMb9jD28v7Sjuwv7XzM/kqvPF/jJFazZk3279/PDz/8AEC7du1YvHixukpxy5YtdO3aNcOB/pdSmTBsfXBwMAUKFFBPDw4Oxs7OLtX1Psf7N+S+FO3odH8VhjKBCR0Lx72P45fCc9NZIbkjrX5XP+5+qhfWzja6iU1H5P2lHdlf2jGY+8SGDh1KzZo1effuHZaWlowdO5bQ0FB2796NqakpHTp0YMqUKToLzMnJCaVSyYkTJyhfvjwAb9++5fz580yePFln2xFCU6YWpkkKQ3z+t42HpwO1amOdxyr141ZrW+Ncx0Vn8QnxJdLqdOKHR0SWlpYsWLCABQsWpLFW2iIiIggICAAgPj6ex48fc/36dWxsbChQoAD9+/dnzpw5FC5cGDc3N7y9vcmePTtt27bN8DaF0BWvje3UjzMylMzur3eqH1f9vjpfDamqs9iE+FLotdup06dP06JFi2TTO3XqxJIlS9Q3O69duzbJzc4lSpTQQ7T6I6cvtKPv/fXk4mO2t8t4kVOB6gXx2tgOhcmnKQzR9/4yNrK/tJPZ+0urJBYaGsr27dt58OABoaGh6hJ7dWMKBb/88ovOg/zSyYdGO4a0vyJfRLKy0pL0F0xDv78GY2llqaOIkjOk/WUMZH9px2CuiR07dozu3bsTGRmJlZUV1tbWyZaRkmIhkspul119HS0+Lp6tXr/y/Jp2g8cuLbVQ/bjrsR7kdrPVaYxCGDONk9i4ceOws7Njw4YNlCxZMjNjEuKzZGJqQsc9/1M/Pz3tJH7LL2nVxoZ6a9SPmy9vhWujz++IIDA8hql+4TyLisMhmynjylvhZGXYtygI/dE4iQUEBDB58mRJYELoSM2xtak5tjYA/vtv/38PIJrb13e3+nGlgZWpNrKmLsPTi8DwGDwPveJ+eJx62qXg9+xqZCuJTKRI4yTm6upKREREZsYixBercLOiDA1M6Ifx5e1gNjVcp9X6fy66wJ+LLgBgX86Bdjs6YWJqovM4M9tUv/AkCQzgfngcU/3CWeEhHS+L5DROYmPHjmXkyJG0adMGZ2fnTAxJiC9bnqJ51dfR3oW9ZWkZ7Yqlgq48Y6HLHPXzb64PIkuuLDqNMbM8i4pLcXpQKtOFSDWJff/998mm2djYULlyZWrVqoWjoyOmpqZJ5isUCry9vXUfpRBfKMtcWdQJTRWvYkfHLTy58FirNpZ9kAS7HOpOnmJ501havxyymaY43T6V6UKkWmJvY6N91zgKhYKQkJCPDkokJSW92vlS9td57zNcXKjdCNYfavJLc4q0KGZQ+yula2KFrEwN6pqYIe0vY2BQ94kJ/ZAPjXa+xP0VcOQue3vvyvD65XpVoNaEOroL6CMkVicGRcVhb4DViV/i++tjSBIT8qHR0pe+v17fC0kyeKe28pTIS6e9XTExM77CkE/hS39/actgbnb29fXl3LlzDBs2LMX5c+fOpXr16nz11Vc6C04IoT0b19zq62jvI96zvNwi4t5rXhjx8lYwC13/LQzpe2UAWXNn03mcQuiCxklsxowZKfbSkeivv/7izJkz7NixQxdxCSF0wCKHBYP8vwMSRmLf3X0HgaceaNXG8nKL1Y877vsfytL2ugxRiI+icRK7fv06w4cPT3V+pUqVpDJRCAOmUCjwXP/vCBAXf/Hl/KwzWrXxW/ON6scN5zSheBvp/EDol8ZJLCoqKt2+EeVmaKFPiQUBAS8tcXkaYnAFAYbmq0FV+GpQFfz9/bF4asaubtqdRTk87ACHhx0AoFQXd+r91CAzwhRG6lN9HjUu7KhZsyYODg5s3bo1xflt27bl6dOnnDt3TqcBCrmQrAljKM02VP99f4U9DGVtzZUZbi+XszVdj/TA1OLzvLdLPo/p+5SfR43Lj7p168aRI0cYOXIkr1+/Vk8PCQlhxIgRHD9+nK5du+o0OCE0lVZ3RUI7uQpaMzRwOEMDhzPg7yFY5NRuGJiwB6H8Ungu8528me/kTeSLyEyKVBiqT/l51Ph0Yp8+fbhx4wYrVqxg5cqV2NnZAfDixQtUKhWdO3emf//+Og9QCE1Id0WZwzybBf1vDAYSCkP2f7OHe4f8tWrjw/HU2vt0xqF8Pp3GKAzPp/w8apzEABYsWEC7du3Ys2cPDx48AMDZ2ZlWrVpRo0YNnQcnhKaku6LMp1AoaL68lfq538pLnJ5yUqs2tnr9qn5cd3oDSnd211F0wpB8ys+jRtfEoqOjGT58OA0bNqRVq1bpLS50TM7Bp0+uiWWcLt5fj849ZGenlK+Xa6JY6xI0nNPEKAbWlc9j+j7l51Hjwg5HR0emT59Ot27ddBqASJ98aDSjroZ6FYGLbQ6pTtSQrt9fb568YU215RleP4d9Drqf6o1ZFq1OFH0y8nnUzKf6PGr8LilXrhw3btzQeQBC6IqTlTkrPHLj7/+KwoUL6jucL1ZOx5zqHkNi38aytuYKrYo7IoIiWFR0nvp5T99vsHKw0nWYIpN9qs+jxkls+vTptG3bliJFitC9e3csLCwyLSghxOfBLIsZvf9MKPhSqVQc+vZ3bu/6W6s2VldZpn7cdmsHHCsX0GmMwrhpfDqxcuXKhIaGEhwcjJmZGfb29mTNmjVpYwoFvr4ZHxpCpExOX2hH9pd29LW/rm+4yolxRzO8vsekupTtUV6HEWlG3l/aMZgOgPPkyUPevHnljyeE0IkyXctSpmtZAJ7++YRtbTdrtf6pScc5Nek4AG5NCtN0SUujKAwRuqVxEtu/f39mxiGE+ILlq+Sovo4W8TyCVV8t1Wr9uwf8WeA8G0gYDbvXhW8wzypFPV8Cwyz/EUJ8sXIoc6gTWtz7ODbUX0NYYKjG678Le8viYvPVz3uc6UPOArl0HaYwEBonsbNnz2q0XPXq1TMcjBBCfMjUwpSv/+itfn505CFubtGuSnpNjRXqx54b2+JU01lX4QkDoHESa968uUbnm0NCQj4qICGESE39mY2oP7MRADe33ODoyENarb/rf9vVj6uPqUXFfjKIr7HTOInt3bs32bS4uDgePnzIunXriI+PZ+LEiToNTgghUlOyQ2lKdigNQNC1Z2xpuUmr9c9O/4Oz0/8AwKl2IVqtbS2FIUZI4ySWVt+IXbp0oUmTJpw5cwYPDw+dBCaEEJqyd3dQX0eLehnJigpL0lkjqcCT99WFIaaWpvS9MhCL7HIvrDHQSWGHiYkJrVu3Zu7cuYwdO1YXTQohRIZky5P938KQmDg2N9vAq9svNV4/7l0cS0osUD/vdrIXNoVsdB6n0A2dVSe+fv2asLAwXTUnhBAfzdTclP8d/lr9/OTEY1xbe0WrNtbXXqV+3HK1F0iPZgZF4yT26NGjFKeHhYVx7tw5Fi5cSNWqVXUWmBBC6FrtH+tR+8d6APzjc4tD3/6u1fp7evqoH7f5rT35q0pG0zeNk1iZMmVSveipUqmoVKkSc+fO1VlgkFA4Mn36dLZu3crz589RKpW0b9+e0aNHY2Ymt7gJITKumFcJinmVACD45gt+bbpeq/V3dPx36JkG3o0p3rakFIbogcaZ4Jdffkn2B1IoFFhbW1OoUCGKFSum8+DmzZvHypUrWbJkCSVKlODmzZv0798fCwsLRo4cqfPtCfElUQ+V8dISl6chX/TQNXlL2qmvo0W/jmZ52UVarX9k+EGODD8IwFdDq1JpQGWDHUrmc6NxB8D60KFDB2xsbFi69N8uaPr168fr16/ZsmWLHiP7tKTDUe3I/kqfDCKqmfi4eLa12UzQlWcZWr94mxLUGFubbLbZdByZ8cjsz6NJegscO3aMtm3bUqlSJRo1asSSJdqVrn6MKlWqcObMGe7cuQPAP//8w+nTp2nQoMEni0GIz9FUv/AkCQzgfngcU/3C9RSRYTIxNaHDri4MDRzO0MDh1PhBu1uI/t5xixXlFzPfyZttbTcT4v8qkyL9cqV5JHb48GE6duyISqXC1taW0NBQ4uPjGTlyJKNHj8704FQqFVOnTmXOnDmYmpoSGxvL8OHDGTduXKrr+Pv7Z3pcQhi7ftctufzGNNn0CrniWFr6nR4iMj7hgW/4a951Xt/QPjFlyZuFMsPLkaeCXSZE9nlJ7yguzSTWvHlz7t69y86dOylRogShoaF0794dPz8/7t69i6Wlpc4D/tCOHTuYMGECkydPplixYty4cYPRo0czefJkunXrlqnbNiRyekw7sr/S1+dUCNsCopNNb+eSlRUeufUQkfFI6f0V9SqKM9NO8veOWxlqs97PDSnZoTQKk8+vMCSzP49pJrGCBQsyePBgRowYoZ7m5+dH/fr18fX1pUiRIpkWGEDJkiUZNGgQ/fv3V0+bNWsWv/76K1euaHevhzGTL2XtyP5Kn1wTy7j03l+xb2O5tOQCF+adz1D7lQZW5qshVTDL8nn8HfQ6KGZ4eDhOTk5Jpjk5OaFSqXj+/HmmJ7GoqChMTZOe8jA1NSU+Pj5TtyvE587JypxdjWwTqhNfReBim+OLrk7UJbMsZlT5rjpVvquOSqXi7+031ZWLmvhz0QX+XHQBgCIti1FrfB2y22XPrHCNXro1oCmV1QOfJJE0btyYefPm4eTkRLFixbh+/TqLFi2iY8eOmb5tIT53TlbmrPDIjb//KwoXlpt2M4NCoaBEu1KUaFcKgMfnH3JkxCHePNKsd6M7e/7hzp5/ALAv50C9nxuSp1jeTIvXGKV5OtHGxgYPDw9cXV3V096+fcvmzZtp1KgR+fLlS9qYQoG3t7fOggsPD2fatGns27ePly9folQqadOmDSNHjiRLliw6246hk9Nj2pH9pR3ZX9rR1f56HRDC8bFHeXzuodbrZs2TjYbejXGu4/LRcWQ2vV4Ts7HRrtNLhUIh44llAvmS0Y7sL+3I/tJOZuyv6NfRnJ3+h9YDfiaqM7U+pbu4G2RhiF6vib1+/TrTNiyEECJBVpus6gE/Y9/F4rfsT87PPqvx+ifGHeXEuKMAlO9bkSrfVcM825cxlIz0iyL0KrHro2dRcThkM5XiAvHFM7M046shVflqSFVUKhW3d//DkeEHiI/RrA7Bb/kl/JZfAsCtaRE8JtUlhzJHZoasV5LEhN6kVOZ9Kfi9lHkL8f8UCgXFPItTzLM4AE/+fMzRkYcIDdDsLNnd3+9w9/eEHo/sSimpN6MhdqWUmRavPkgSE3qTVtdHcsOtEMk5VspP9xO9AAgNDOXk+KMEnnqg0bov/nrO5mYbAMhinYUGs5vgUt81nbUMnyQxoTfPouJSnB6UynQhxL+snazxXN8WgLdhbzk74zR/bbqm0bpvQ9+yt9e/Y6N5TKpLmW5lMTFNtztdgyNJTOiNQ7bkffcB2KcyXQiRsiy5slDvpwbU+6kBcTFxXFl1mbPT/9B4/VOTjnNq0nEAyvaqQNVh1bHIYRyFIZLEhN6MK2/FpeD3ybo+GlfeSo9RCWHcTM1NqdjvKyr2+wqVSoX//tsc+f4gsW9jNVr/6qrLXF11GQCXhm7U/rEuVvlyZmbIHyXDSezNmzcA5MxpuC9OGLYPuz4KiorDXqoThdAphUJBkebFKNI8YdDiZ35POTrykMZDwgQcvkvA4bsA5Cmel/ozG6EsY59p8WaEVknswYMHTJ8+ncOHDxMWltBtSq5cuWjYsCGjRo3CxcXw7x4XhiWx6yMhROZzKJ+Prkd7APDmURgnJhzjwfEAjdZ9+Xcwv7XYCIBFDgsaeDfGtXHhZF0TfmoaJ7ErV67QqlUroqOjqVu3Lm5ubgDqoVoOHDjA7t27KVeuXKYFK4QQQjdyFshFqzWtAXj35h3nvc9wbZ1mo4O8j3jP/n571M9rjvWgbM8KmJh9+sKQNLud+lCjRo149OgRe/fuTdKXIiQkshYtWuDk5MTBg5r31iw0I90CaUf2l3Zkf2nnc99f8bHxXF3jx+mpJzO0fpluZak2oiaWORPGm9Rrt1Mfun79OiNGjEiWwADc3Nzo3bu3Tjv/FYZLetkQ4vOk/mznd8JhfU/GlstB7LlAjgw/wPvw9xq1cX39Va6vvwqAc10XnHu6QibmfI2TmJ2dHSYmqR8qmpqaolR+XneCi+Sklw0hPk+pf7ad6f/XEACeXw/i6KhDvLwVrFGbD44H8OB4AKpJcZTtUT5T4tb4BObAgQNZuXIljx8/Tjbv4cOHrFy5koEDB+o0OGF40uplQwhhvDT5bCvL2NPlQHeGBg6n5/m+uDZy06jtP6acIDokSqfxJtL4SCw6OppcuXJRqVIlmjRpoq5EvHfvHgcPHsTV1ZWoqCjmz5+vXkehUDBkyBDdRy30RnrZEOLzpO1n2ypfTpov9wQSCj1855zlyv/fX/ZfChMFJmaZ04mBxkls0qRJ6sc+Pj7J5t+8eZObN28mmSZJ7PMjvWwI8Xn6mM+2RQ4Lak2oQ60JdYiPi+fGhqucnJjQA0g2x+zUGF5LXeihaxonsWvXNOuTS3zepJcNIT5Puvpsm5ia4P51edy/TrgGZjDViQULFsy0IITxkF42hPg8GetnW+tup8LDwzl9+jQPHz4EEpJbzZo1sbKSX+JfCullQ4jPkzF+trVKYgsWLGDmzJlERUWhUv17j3S2bNkYOXIkQ4cO1XmAQgghRGo0TmILFy5k4sSJ1KhRg969eyfpdmrFihX8+OOPmJiYMHjw4EwLVgghhPiQxkls2bJl1KlThx07diTp8LFkyZK0bNkSLy8vli1bJklMCCHEJ6Pxzc4hISE0bdo0xR6LFQoFzZs3JyQkRKfBCSGESFtgeAx9ToXQ/EAwfU6FEBgeo++QPimNj8Tc3d35559/Up3/999/4+7urpOghBBCpE+6gdPiSGzWrFns3buX+fPnExERoZ4eERHBvHnz2L9/v3QALIQQn5B0A6fFkVivXr1QKBT8+OOPTJkyBTs7OwBevHhBfHw8SqWSnj17JllHoVDg6+ur24iFEEIA0g0caJHE8uTJQ968edVViYlkNGchhNAP6QZOiyS2f//+zIxDCCGElqQbuAz02CGEEMIwGGtXUbqkdRKLiYnhzp07vHnzhvj4+GTzq1evrpPAhBBCpM8Yu4rSJY2TmEqlYsqUKaxYsYLIyMhUl5N7xYQQQnwqGpfYz5s3j7lz59KmTRuWLl2KSqVi0qRJzJ07l+LFi1O6dOkUxxn7WEFBQfTr1w9XV1eUSiWVK1fmzJkzOt+OEEII46PxkdjGjRtp2bIl8+bNUx9tubu74+HhQceOHalXrx5nzpzBw8NDZ8GFhobSqFEjqlSpwtatW7G1tSUwMJC8efPqbBtCCCGMl8ZHYo8fP1YnKBOThNXevXsHgKWlJR06dGDz5s06DW7BggXY29uzbNkyKlSogLOzMx4eHhQtWlSn2xFCCGGcNE5i1tbWvH37FoCcOXNiYWHBkydP1PMtLS11fj1s//79VKhQgR49euDm5kaNGjVYvnx5kmFghBBCfLkUoaGhGmUET09PHBwcWLJkCQBNmjQhPDycLVu2EB8fT4cOHTA3N+fUqVM6C06pVAIwYMAAPD09uXHjBqNGjWLixIn07ds3xXX8/f11tn0hhBD6Vbhw4TTna5zENm3axKpVq/j999/JkiUL58+fx8vLi/fv3wNgbm7Or7/+Sr169T4+6v+XN29eypUrx+HDh9XTJk+ezL59+7h48aLOtmPo/P390/1Din/J/tKO7C/tyP7STmbvL40LO7p06UKXLl3Uz6tWrYqvry8HDhzA1NSUevXq4erqqtPglEplsutfRYoU4fHjxzrdjhBCCOP0UT12ODs7079/f13FkkyVKlW4e/dukml3796lQIECmbZNIYQQxkPjwo7/evPmDQMHDuTOnTu6jCeJAQMG8Oeff+Lt7U1AQAC7du1i+fLl9O7dO9O2KYQQwnhkOIm9ffuWzZs38+zZM13Gk0T58uXZtGkTPj4+VK1alSlTpvDDDz9IEhNCCAF85OnET1Hq3qhRIxo1apTp2xFCCGF8MnwkJoQQQuhbhpOYhYUF1atXx9raWofhCCGEEJrL8OlEa2tr9u3bp8tYhBBCCK3I6UQhhBBGS6vxxNauXcuGDRt48OABoaGhyZZRKBS8evVKl/EJIYQQqdI4iU2YMIFFixZRunRp2rdvL9fChBBC6J3GSWzz5s20bNmStWvXZmI4QgghhOY0vib29u1bateunYmhCCGEENrROInVqlULPz+/zIxFCCGE0IrGSWz27NlcunQJb29vXrx4kZkxCSGEEBrR+JpYuXLlUKlU/PTTT/z000+Ym5tjYpI0ByoUCp4+farzIIUQQoiUaJzEvLy8UCgUmRmLEEIIoRWNk9iSJUsyMw4hPlpgeAxT/cIJeGmJy9MQxpW3wsnKXN9hic+EvL8M00f1Yi+EoQgMj8Hz0Cvuh8cBplx+E82l4PfsamQrXzTio8n7y3ClmcQuX76sdYMVKlTIcDBCZNRUv/D//4L51/3wOKb6hbPCI7eeohKfC3l/Ga40k1j9+vU1vg6mUqlQKBSEhIToJDAhtPEsKi7F6UGpTBdCG/L+MlxpJrFFixZ9qjiE+CgO2UxTnG6fynQhtCHvL8OVZhLr3Lnzp4pDiI8yrrwVl4LfJznlU8jKlHHlrfQYlfhcyPvLcMlQLOKz4GRlzq5GtrRzyUqFXHG0c8kqF92Fzsj7y3BJdaL4bDhZmbPCIzf+/q8oXLigvsMRnxl5fxkmORITQghhtCSJCSGEMFqSxIQQQhgtSWJCCCGMliI0NFSl7yCEEEKIjJAjMSGEEEZLkpgQQgijJUlMCCGE0ZIkJoQQwmhJEhNCCGG0JIkZqDlz5lCnTh0KFCiAq6srHTp04NatW/oOy2jMmTMHa2trRowYoe9QDFZQUBD9+vXD1dUVpVJJ5cqVOXPmjL7DMkhxcXFMnTqVMmXKoFQqKVOmDFOnTiU2NlbfoRmEs2fP0rFjR4oXL461tTWbNm1KMl+lUjF9+nSKFSuGvb09zZo14++//9bJtiWJGagzZ87Qq1cvDh06xJ49ezAzM8PT05PXr1/rOzSD9+eff7J27VpKliyp71AMVmhoKI0aNUKlUrF161YuXLjAzJkzyZs3r75DM0jz5s1j5cqVzJgxg4sXL/Lzzz+zYsUK5syZo+/QDEJkZCQlSpTg559/JmvWrMnmz58/n0WLFjFjxgyOHz9O3rx58fLyIjw8/KO3LfeJGYmIiAgKFizIpk2baNKkib7DMVhhYWF4eHiwYMECZsyYQYkSJZg1a5a+wzI4kydP5uzZsxw6dEjfoRiFDh06YGNjw9KlS9XT+vXrx+vXr9myZYseIzM8jo6OzJw5ky5dugAJR2HFihWjT58+DB8+HIDo6GgKFy7MlClT6NGjx0dtT47EjERERATx8fFYW1vrOxSD9u2339KqVStq1aql71AM2v79+6lQoQI9evTAzc2NGjVqsHz5clQq+U2bkipVqnDmzBnu3LkDwD///MPp06dp0KCBniMzfIGBgTx//py6deuqp2XNmpVq1apx4cKFj25fhmIxEqNHj6Z06dJ89dVX+g7FYK1bt46AgACWL1+u71AM3oMHD1i1ahUDBgzg22+/5caNG4waNQqAvn376jk6w/Ptt98SERFB5cqVMTU1JTY2luHDh9O7d299h2bwnj9/DpDsVHXevHl59uzZR7cvScwI/PDDD/j6+nLw4EFMTWU49JT4+/szefJkDh48iLm5DFSYnvj4eMqVK8fEiRMBcHd3JyAggJUrV0oSS8HOnTv57bffWLlyJcWKFePGjRuMHj2aggUL0q1bN32H90WTJGbgxowZw86dO9m7dy/Ozs76DsdgXbx4kVevXlGlShX1tLi4OM6dO8fq1at5+vQplpaWeozQsCiVSooWLZpkWpEiRXj8+LGeIjJsEyZMYNCgQbRp0waAkiVL8ujRI+bOnStJLB1KpRKA4OBgChQooJ4eHByMnZ3dR7cvScyAjRo1Ch8fH/bu3UuRIkX0HY5Ba9asGeXKlUsybeDAgbi6ujJs2DAsLCz0FJlhqlKlCnfv3k0y7e7du0m+ZMS/oqKikp0FMTU1JT4+Xk8RGQ8nJyeUSiUnTpygfPnyALx9+5bz588zefLkj25fkpiBGj58OFu2bGHjxo1YW1urzytnz56dHDly6Dk6w2NtbZ2s6CVbtmzY2NhQokQJ/QRlwAYMGEDDhg3x9vamdevWXL9+neXLlzN+/Hh9h2aQGjduzLx583BycqJYsWJcv36dRYsW0bFjR32HZhAiIiIICAgAEk5VP378mOvXr2NjY0OBAgXo378/c+bMoXDhwri5ueHt7U327Nlp27btR29bSuwNVGpViKNGjWLMmDGfNhgj1axZMymxT8OhQ4eYPHkyd+/eJX/+/PTp04dvvvkGhUKh79AMTnh4ONOmTWPfvn28fPkSpVJJmzZtGDlyJFmyZNF3eHp3+vRpWrRokWx6p06dWLJkCSqVip9//pm1a9cSGhpKhQoV8Pb21skPTEliQgghjJbcJyaEEMJoSRITQghhtCSJCSGEMFqSxIQQQhgtSWJCCCGMliQxIYQQRkuSmDAYpUuXVnfrk5maNWtGs2bNMn07X4orV65ga2vLvXv31NP69+9P6dKl9RiV/owfP5569erpO4wvhiQxobZp0yZ1zxfnzp1LcZly5cphbW1t8EngwoULTJ8+ndDQUH2HkqJbt27Rv39/SpUqhZ2dHQULFqRJkyYsXbqUt2/f6nx7UVFRTJ8+ndOnT+u87cmTJ9OyZUtcXV113vaHdu7cyTfffEOlSpWwsbGhUqVKqS777t07Jk2aRPHixbG3t6du3bocO3YsxWVv375N27ZtyZ8/P87OzvTp04cXL16kuOyGDRuoXLkySqWScuXKsXjx4mTD1wwcOJC//vqL33//PeMvVmhMkphIJkuWLGzfvj3Z9D///JP79+8bRQ8FFy9eZMaMGYSFhSWb5+Pjg4+Pjx6iSrBp0yZq1arFsWPH8PLywtvbm7Fjx2Jvb8/48ePVQ6LoUnR0NDNmzODMmTM6bffGjRucOHEi2cCGCxYs4NKlSzrd1qpVq9i3bx958+ZNt+PYAQMG8Msvv9CmTRt+/vlnzM3Nad++fbLX/+TJE5o2bcrdu3cZN24cQ4YM4ciRI3h6eib7MbFmzRoGDx5MkSJFmDVrFpUrV+aHH35g9uzZSZazt7enSZMmLFy4UDcvXKRJ+k4UyTRo0IBdu3YxY8aMJMOabNu2jSJFiuh8OJioqCiyZcum0zbTos/OgC9fvszQoUMpX74827ZtI1euXOp533zzDffu3TOq0ZY3btyInZ0dNWrUSDI9M4bDWbZsGQ4ODpiamtKsWbNUj5YuX77Mjh07mDRpEt9++y0AHTt2pGrVqkyYMIHjx4+rl50zZw4RERGcOHGCggULAlC+fHk8PT3ZuHGjeryw6OhopkyZQv369dmwYQMA3bp1Iz4+ntmzZ9OjRw9sbW3V7Xp5edG9e3fu3buX6UeoXzo5EhPJtG3bltevXyc5/RIXF4ePj0+qHXYuXLiQxo0b4+LiglKppFq1aqxfvz7ZconXvU6dOkW9evVQKpXMnz8/1Vh8fHzImzcv3333nfq0jZ+fH+3ataNgwYLY29vTuHFj/vjjD/U606dPV3dk6+7urj5Fmngq7b/XxAIDA7G2tmbu3LmsW7eOsmXLYmdnR506dfDz80sW065du9SnlKpWrcqePXs0vgb0888/o1KpWLFiRZIElsjV1ZUBAwaon8fFxeHt7U25cuWws7OjVKlSTJgwgejo6CTrXb16lXbt2uHq6opSqaRUqVL07duX6OhoAgMD1V+kM2bMUO+P/v37Awmdt44bN44yZcpgZ2eHq6srzZo14+zZs+m+nv3791OrVi1MTJJ+laS0P6ytrfnuu+/Yt28fVatWxc7OjipVqnD06NF0twOQP39+jX5A7d69GxMTE77++mv1tCxZstC1a1f8/PwIDAxUT9+zZw8NGjRQJzCA2rVr4+bmluRo/fTp04SEhNCzZ88k2+rduzfR0dEcPHgwyfTatWujUCjYv3+/Rq9NZJwciYlk8uXLR9WqVdm+fTuNGzcG4OTJkwQHB9OuXbsUT8UtWbKERo0a4enpiUKh4Pfff2fIkCHExsYm++AHBATQrVs3unfvTteuXcmfP3+KcWzevJlBgwbRt29fpk+fDsCZM2do06YNpUuXZsSIEZibm7NlyxZat26Nj48PNWvWpEWLFty7d4/t27fz008/qX8h/3f8rP/auXMnkZGR9OjRA4VCwfz58+natStXr15VH1kcOnSIHj16ULx4ccaPH09YWBhDhgwhX7586e7X6OhoTp48SbVq1XByckp3eUgYUXjDhg20aNGCgQMHcuXKFRYsWMDff//N1q1bUSgUvHz5Ei8vL2xtbRk6dCjW1tY8fvyYAwcOEBUVRZ48eZgzZw7Dhg2jefPm6o5aCxUqBMCwYcPYtWsXvXv3plixYoSFhXHp0iX++usvqlevnmpsT58+5fHjx5QtW1aj1wIJp3kPHjxIz549yZEjB8uWLaNbt2789ddf5M6dW+N20nL9+nUKFSqUrBPtChUqqOc7OTnx9OlTgoODkw3hAwlHYx8mpuvXrwMkW7Zs2bKYmJhw/fp1unTpop6eK1cuChUqhK+vL0OGDNHJ6xIpkyQmUtSuXTvGjh1LZGQk2bNnZ+vWrVSsWFH9xfdfly5dSnJK8JtvvsHLy4uFCxcmS2L379/n119/pWnTpqluf82aNQwbNoxhw4apj6pUKhXfffcdVapUYdeuXere1nv27EmtWrWYMmUKhw8fplSpUri7u7N9+3aaNWumccJ48uQJfn5+6i8/Nzc3OnfuzLFjx9TJ/Mcff8Te3p5Dhw5hZWUFgIeHB82bN093LK579+4RExNDyZIlNYrnr7/+YsOGDXTu3JnFixerp+fPn58ZM2Zw6NAhGjduzIULF3j9+jU7d+5M8iX7ww8/qB+3atWKYcOGUbJkSTp06JBkO4cOHaJ79+789NNPGsWVyN/fH0CrwVrv3LnDhQsXcHFxAaBmzZrUqFGD7du362xE6efPn6sHYvxQ4rSgoCD1ch9O/5C9vT1v3rxRv/8Tl7W3t0+ynIWFBblz51a3+SFnZ2du3779cS9GpEtOJ4oUeXp6EhMTw/79+4mOjub333+nXbt2qS6fmMBiYmJ4/fo1r169ombNmty/fz9ZcYWjo2OaCWzx4sUMGzaMH374Icn4Vjdu3MDf35+2bdsSEhLCq1evePXqFeHh4dSuXZtLly4RFRWV4dfcsmXLJL/eq1WrBsCDBw8AePbsGbdu3aJ9+/bqBAZQo0YNjYaUCA8PB9B4PLjDhw8DCdVuHxowYACmpqbq+Tlz5gTg4MGDxMTEaNT2h3LmzMmlS5d4+vSpVuuFhIQAqQ8blJKaNWuqExhAqVKlyJkzp3of60J0dHSKo3gnFiQlFmwknpJNadnEaR8ua25unuIwNVmyZEmxotTa2ppXr15l8FUITcmRmEiRjY0NdevWZdu2bZiZmREVFUXr1q1TXX7//v3MmjWLGzduEBcXl2Temzdvklz/SevI6MKFCxw7doxBgwYxYsSIJPMS70MaPHgwgwcPTnH9kJCQDBeJ/Pe0ZuKXc2KZ/qNHjwCSfAkncnFx4dq1a2m2n5j4IiIiNIrn0aNHKBQK3NzckkzPlSsX9vb2PHz4EEhIoq1atWLGjBksXryYatWq0bRpU9q2bUv27NnT3c6UKVPU5f5lypShfv36dOjQgcKFC2sU539LzNOS0qnjXLly6fRWiKxZs/Lu3btk0xMTTWIyy5o1K0CKyyZO+3DZmJgY4uPjk13/e/v2bYoVuyqVSsZm+wQkiYlUtWvXjn79+qmPdPLmzZvicufPn+d///sfVapUYe7cudjb22NhYcHhw4dZvHhxsiHcE788UlKkSBGioqLYtm0bPXr0SFLZldjOpEmTUr0OkydPHi1f5b9SKxrQ5ks6LS4uLpibm3Pr1i2dtJdIoVCwbt06Ll++zMGDBzl58iRDhw5lzpw5HD16NNW/WyJPT0+qVq3K77//zvHjx1m2bBnz589n8eLFaR59J17D0iYBZfY+hoTTg4k/OD7031OCiacRE6d/KCgoiJw5c6p/BHy4rIODg3q59+/fExISkuw0IyTslw8rFkXmkNOJIlVNmzbF0tISX1/fNIcR3717N1myZMHHx4du3brRsGFDateunWaySo2NjQ27du0iW7ZstGrVKsmXUeL1uBw5clC7du0U/2XmPWyJ17wSh2H/UErT/itbtmx4eHhw9uxZ9VFUettTqVTcvXs3yfQ3b94QFBSUpKIOEgoXxo4dy5EjR9i2bRuBgYGsW7cOIN0jAqVSSY8ePdiwYQM3btzA2dlZXUyTmsRCmQ+r/QxB6dKluX//frLkmnjfWmLVZL58+ciTJw9XrlxJ1oafnx+lSpVK0iaQbNkrV64QHx+fYmVqYGAgRYoU+ajXItInSUykKlu2bMyePZtRo0bRvHnzVJczNTVFoVAkOeIKDQ1l48aNGdquvb09u3btQqVS4enpqf6lXLZsWVxcXFi0aJH6+tKHXr58qX6c+Atal6epHBwcKFGiBFu3bk2y/TNnzmh8dDVy5EggofDlzZs3yebfv3+fJUuWANCwYUMA9fNES5cuJS4ujkaNGgEJr/G/RzLu7u4A6uuRiT8o/rs/4uLikl2ztLa2xsnJKcUbxT9kb2+Pk5MTV69eTXO5T61Vq1bEx8ezdu1a9bR3796xadMmypYtm6QQpWXLlhw5ciTJj4pTp05x9+5dPD091dNq1aqFjY0Nq1evTrKtVatWkSVLFnXhT6KwsDDu379P5cqVdfraRHJyOlGkqWPHjuku07hxYxYtWoSXlxcdOnTg9evXrFu3Djs7uxRP1WjCycmJXbt20bRpU7y8vNi/fz82NjYsXLiQtm3bUqVKFbp06YKjoyPPnj3j7NmzqFQq9u3bB/xbCj158mTatm2LhYUFtWrVSvfUWnrGjx9P586dadSoEV26dCEsLIwVK1ZQokQJja51ffXVV8ydO5fvvvuOSpUq0aFDB9zc3Hj79i0XL15kz549dO7cGUgoeujatSsbNmzgzZs31KpVi2vXrrFx40bq16+vTnK//vorK1eupHnz5hQqVIjo6Gg2bdqEqakprVq1AhKSWPHixdm5cydubm7kzp0bJycn3NzcKFGiBC1atFAXWfj6+nL06FH69OmT7utp2rQpO3bsSPFaka6dPXtW3R3ao0ePiIyMZNasWUBCEU7i7QAVK1bE09OTqVOn8urVK1xdXfntt98IDAxMdntI4u0FLVu2pF+/fkRHR7NgwQKKFStGt27d1MtlzZqVsWPHMnz4cLp27UqDBg04d+4cW7duZcyYMclOY588eRKVSpVmAZPQDUli4qPVrFmTJUuWMHfuXMaMGUO+fPno27cv1tbWDBo0KMPtFi5cmJ07d9KiRQtat27N7t27qV69OkeOHGHWrFmsWrWK8PBw7OzsKF++fJIvnXLlyjFx4kRWrVrFwIEDiY+PZ+/evR+dxJo0acKqVav4+eefmTx5Mq6urixatIjffvuNf/75R6M2unXrRvny5fnll1/YuXMnL168IEuWLJQsWZJp06YleR3z5s3DycmJjRs3cuDAAezs7Bg8eDBjxoxRnyKsXr06V65cwcfHhxcvXmBlZUWZMmWYOXMmFStWVLe1cOFCRo0axbhx43j37h2dOnVi/vz59O7dmxMnTnDgwAFiY2NxcnJSF3uk53//+x9LlizhzJkz1KpVS8u9qZ0//viDGTNmJJk2bdo0AEaNGpXknralS5fy008/sXXrVl6/fk3x4sX57bffksWYP39+9u/fz7hx45gyZQpmZmY0aNCAadOmJTs13bt3bywsLPjll184dOgQ+fLlY+rUqcmqR+HfG+I1LY4RGacIDQ3V3RVVIb5QNWrUIE+ePOzatUvfoXxyrVu3JleuXKxZs0bfoRiEoKAg3N3dWbVqVZqn4YVuyDUxIbQQExNDbGxskmmnT5/mr7/+StZ/4Jdi/Pjx7NmzJ8lQLF+yRYsWUbJkSUlgn4gciQmhhcDAQDw9PWnfvj0ODg7cuXOHNWvWkDNnTs6fP6+zrpOEEJqRa2JCaMHa2pqyZcuyYcMGXr58SbZs2WjYsCGTJk2SBCaEHsiRmBBCCKMl18SEEEIYLUliQgghjJYkMSGEEEZLkpgQQgijJUlMCCGE0ZIkJoQQwmj9H7jwMoKqDHN/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "toy_panel = pd.DataFrame({\n",
    "    \"mkt_costs\":[5,4,3.5,3, 10,9.5,9,8, 4,3,2,1, 8,7,6,4],\n",
    "    \"purchase\":[12,9,7.5,7, 9,7,6.5,5, 15,14.5,14,13, 11,9.5,8,5],\n",
    "    \"city\":[\"C0\",\"C0\",\"C0\",\"C0\", \"C2\",\"C2\",\"C2\",\"C2\", \"C1\",\"C1\",\"C1\",\"C1\", \"C3\",\"C3\",\"C3\",\"C3\"]\n",
    "})\n",
    "\n",
    "m = smf.ols(\"purchase ~ mkt_costs\", data=toy_panel).fit()\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase)\n",
    "plt.plot(toy_panel.mkt_costs, m.fittedvalues, c=\"C5\", label=\"Regression Line\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\")\n",
    "plt.title(\"Simple OLS Model\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이제 여러분은 인과추론에 대해 많이 계시므로, 고정 효과 모델을 실행하기로 결정하고 도시의 지표를 모델에 더미변수로 추가합니다. 고정 효과 모델은 시간에 따라 일정한 도시 고유의 특성을 통제하므로, 도시가 여러분의 상품에 대해 덜 개방적이면 이를 포착하게 됩니다. 해당  모델을 돌리면, 마침내 여러분은 더 많은 마케팅 비용이 더 많은 인앱 구매로 이어진다는 것을 알 수 있게 됩니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAE0CAYAAACirQ3aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABoqUlEQVR4nO3dd1gU19cH8O/Si+AiZbFQlKJgR42ggh0QbNhLxGiMNdaoqLGXIDYEY4sl9h47CnYFFAuioFjAgqKASO9t5/2DH/u6zgKL7rK7eD7Pk+cJc2dnzgy4Z+bOmXs56enpDAghhBAFpCTrAAghhJBvRUmMEEKIwqIkRgghRGFREiOEEKKwKIkRQghRWJTECCGEKCxKYkQmuFwu3N3dZR2GgLe3N7hcLoKDg6W6n+vXr8PV1RXm5uascxAREQEPDw9YWlqCy+WiefPmUo2lJmrevLlEzpuktkOkj5IYkSgul1vhf1u2bJF1iBJRlvQq+u/rJB0XF4fhw4cjOjoagwcPhpeXF0aMGAEAyMzMxJAhQ3D79m24ubnBy8sLkyZNkvpxBAcHg8vlftO+3N3dBce6f//+ctdbt26dYL2xY8d+T7iEsKjIOgBSM3l5eYlc3q5dOwDAvXv3oKmpWZ0hSUXHjh3RqVMnkW2mpqZCP9+8eRP5+fnw8vLCzJkzhdoePnyI5ORkjBkzBr6+vlKLVxpUVFSwb98+jBo1itXGMAz2798PFRUVFBcXyyA6UtNREiNSMX/+/Arbra2tqykS6erUqVOlx1omISEBAGBkZFSlNnnn4uKCgIAAREdHw9bWVqjtxo0biIuLQ+/evXH+/HkZRUhqMupOJDLxdXfb27dvYWZmBlNTU7x9+1Zo3YKCAjg6OoLL5eLixYtCbbdu3cKwYcNgYWEBQ0NDNGvWDH/88QeSkpJE7vfRo0cYOHAgGjRoABMTE/Tr1w/37t2T+PF9qazLztvbGwAwZcoUQffawYMHhbrzfHx8hNrK5OfnY9OmTejcuTPq16+PevXqoUuXLti9ezcYRvTIcQ8fPsTYsWNhY2MDQ0NDWFtbo0+fPjh06BCA0i7RPn36AAAOHz4s1BX65b4r4+npCQ6Hg71797La9u7dCw0NDQwdOrTcz2dlZWHlypVo164deDweTE1N0adPHwQEBIhcn2EY/PPPP7C3twePx4ONjQ3mzJmDjIyMCuM8c+YM+vXrB3NzcxgZGcHOzg5Lly5FZmam2MdK5A/diRG5YG5ujk2bNsHT0xNjxoxBUFAQ1NTUAAB//vknoqKiMHnyZPTq1UvwmY0bN2Lp0qXQ09ODs7MzeDwenj59il27duHixYu4fPky6tevL1j/7t276N+/PwoKCtCnTx9YWFjg6dOn6N27N5ycnKR2bKampvDy8kJISAhCQ0Ph5uYmKBpo3rw5vLy8EBUVhQsXLgh1T5atk5WVhf79+yM8PBwtWrQQPEe7evUqZs2ahfv372Pr1q1C+9y3bx9mzpwJJSUluLq6wsrKCikpKXj8+DG2bt2KESNGoFOnTnj37h0OHz6MZs2aCV1UVKWooWHDhnB0dMTRo0exbNkyaGhoAAA+f/6MCxcuoH///uByuSI/m5GRgV69eiE6OhotWrTAxIkTkZGRgdOnT2PkyJGYP38+q2t63rx52L59O3g8Hjw9PaGuro4LFy7gwYMHKCoqgqqqKms/f/zxB3bt2oX69eujd+/e4HK5ePDgATZu3IhLly4hKCgIOjo6Yh8zkR+UxIhUlN11fInH41X4YL9v374YP348/vnnHyxevBirV6/GmTNnsHPnTtjZ2WHZsmWCdUNDQ7Fs2TK0a9cOx48fF/qSPHLkCCZOnIh58+YJCg4YhsHvv/+OvLw87N27F/369ROsv2PHDsyZM+ebjjMkJETksQKlXY2Ojo4wMzPD/Pnz4e3tjdDQULi7u2PkyJGC9Vq0aIGDBw/iwoULIrsnFyxYgPDwcCxduhQzZswQLC8oKMCoUaNw+PBh9O3bV5Dgnz9/jlmzZkFbWxsXL15E06ZNhbYXHx8PAHB0dARQehfWvHlzsbtFRRk9ejR+/fVXnD17FkOGDBFst7CwEJ6enuXeLS5btgzR0dEYOXIk/v77b3A4HADAnDlz0K1bN6xevRo9e/aEnZ0dgNILke3bt8PU1BTXr1+Hvr4+AGDRokXo168fEhMTYWJiIrSPo0ePYteuXejduzd27Ngh9Cx27dq1WLVqFVavXo1Vq1Z98/ET2aEkRqTCx8eHtaxZs2aVVqetXLkS9+7dw7Zt22BmZobVq1dDV1cXu3fvFrrC3rZtGxiGga+vL+sqf9iwYdiyZQsuXLiArKws6Ojo4O7du4iJiUH79u2FEhgA/Prrr9i6dStev35d5eMMDQ1FaGhoue1lieJbpaWl4fDhw2jRooVQAgMAdXV1LF68GJcuXcLRo0cFSWzXrl0oLi7G7NmzWQkMABo0aPBdMYnSu3dv1KlTB3v37hUksX379sHS0hKdOnUS+epCUVERjhw5Ai0tLSxfvlyQwACgfv36mDVrFry8vLBv3z5BEivr5vzjjz8ECQwoPReLFi0S+drGli1boKysjE2bNrGKiWbNmoVt27bh2LFjlMQUFCUxIhXp6enf9Dk1NTX8+++/6Ny5s+DOYO/evTA3Nxda7+7du1BRUcG5c+dw7tw51nYKCwtRUlKCV69eoVWrVnj8+DGA0mrCrykpKcHe3v6bkpiXl9d33cFUJjw8HMXFxVBSUhJ5x1dW8ffy5UvBsgcPHgAAevToIbW4vqauri64eIiNjUViYiJiYmKwfPnycj/z8uVL5Obmom3btkIJqUyXLl0AQPC7+/L/Rf0e7e3toaIi/JWWm5uLyMhI6OnpYdu2bSLjUFNTQ0JCAlJTU1GnTp1Kj5XIF0piRO40bNgQLVu2RHBwMMzMzEReXaempqK4uFjkHd+XsrOzAUDw8N7Q0FDkevJaFZiamgqgtCDl0aNH5a5XdpwABAUO9erVk2psXxs9ejS2bNmCvXv3IikpCaqqqhg+fHi565f9Tso79zweDwCECjYq+j0qKyuzklB6ejoYhkFqaqpYfyuUxBQPJTEid/7++28EBwdDX18fcXFxWLlyJZYuXSq0jq6uLoqKivD+/XuxtqmrqwsASE5OFtn+6dOn74pZWsriHj9+PNasWSPWZ2rXrg0A+PjxY7kFFdLQuHFjODg44NChQ8jJyYG7u3u5Fw3A/x9beee+rMK0bL0v/z85OVlwnGVKSkqQmpqKunXrsta3tbXF7du3v+GoiLyjEnsiV8LDw7Fs2TKYmZnhzp07aNasGfz8/HD16lWh9dq1a4esrCxERUWJtd2WLVsCgMjnV3w+H2FhYd8fvBS0bdsWSkpKuHPnTpU+AwBXrlypdF1lZWUApQlAEjw9PZGSkoL8/HyMHj26wnWtra2hpaWF6OhopKSksNpv3rwJAGjVqpVgWUW/x7CwMNYL1bVq1YKtrS1iYmJE7oMoPkpiRG6kp6djzJgxAIDdu3fDyMgIe/bsgba2NiZMmIDExETBulOmTAEAzJgxAx8+fGBtKz8/X+iLv3379rCyssLdu3dx5swZoXV37dr1Tc/DqoOBgQGGDh2KqKgoeHt7ixz14sOHD0LPxH799VeoqKhg3bp1iI6OFrl+mbLus7KKxe/l4eGBAwcO4ODBg4JnWuVRVVXF0KFDkZubi2XLlglVMCYkJMDX1xccDgc///yzYHnZ6wXr168XdLUCpZWaK1asELmfKVOmoKioCJMnT0ZaWhqrPSsrS/AckSge6k4kcmPq1Kl49+4dVqxYgTZt2gAALC0tsX79ekyYMAG//fYbzpw5AyUlJTg5OWHFihVYsmQJ2rRpg549e8Lc3Bz5+fl4//49bt++DVNTU4SEhAAAOBwONm3aBA8PD4wZM0boPbHr16+jR48eYt25fK2iEnsNDQ3W8FLfYs2aNXj9+jV8fHxw9OhRdOjQATweD0lJSYiNjcX9+/exatUqwSgoTZo0wfr16zFz5kx06dJF8J5YWloaIiMjUVBQIKgWtLKyQoMGDXDnzh389ttvsLCwgLKyMnr16oVmzZpVOVYNDQ307t1b7PWXLFmCO3fuYN++fYiMjESXLl0E74mlpaVh7ty5gjtLoLR4o+w1DAcHB/Tt21fwnljt2rVhbGzM2sfIkSPx+PFj/PPPP2jVqhW6d+8OU1NTZGRk4N27d7h9+za6du0qeAmcKBZKYkQubN++HefOnYOzszN+//13obahQ4fi5s2bOHToEHx8fATVgFOnToW9vT22bduGO3fuIDAwELVq1ULdunUxZMgQeHh4CG3H3t4eFy9exIoVK3D16lVcvXoVbdq0wfnz53H16tVvSmIVldjr6upKJInp6Ojg/Pnz2L9/P44fP47z588jPz8fhoaGMDMzw5IlS1jHOnr0aNja2mLTpk0ICwvDxYsXUadOHTRu3Bjjxo0TrKesrIwDBw5g6dKlCAoKQlZWFhiGQb169b4piVUVl8tFUFAQ/Pz8cPbsWWzZsgXq6upo0aIFJkyYgL59+7I+4+PjA0tLS+zcuRN79+5FnTp10Lt3byxatKjccSzXrFkDZ2dn7Nq1CyEhIUhLS0Pt2rVRr149jBs3DoMGDZL2oRIp4aSnp4t+C5EQQgiRc/RMjBBCiMKiJEYIIURhURIjhBCisCiJEUIIUViUxAghhCgsSmKEEEIUFiUxQgghCouSmAKIiYmRdQgKhc5X1dD5qho6X1Uj7fNFSYwQQojCoiRGCCFEYVESI4QQorAoiRFCCFFYlMQIIYQoLEpipEb5cmJFQkjNR/OJkRqBX8LHw3/uI+NdBkzHNpR1OISQakJJjCi8rIQsXJpxAfFh7wEAqo3UYGVlJeOoCCHVQewkxjAMnj59ihcvXiAlJQUcDgf6+vqwtrZG06ZNweFwpBknISLFXHiJq/MuoSAjX7AsyvcRWrnZQbe+rgwjI4RUh0qT2K1bt3Dw4EFcvHgR2dnZrGcOHA4HtWrVgqurK0aOHInOnTuLvfPQ0FBs2rQJjx8/RkJCAjZv3oyRI0cKrRMbG4ulS5fi1q1bKCoqgpWVFXbs2IHGjRuLvR9S8xTmFOLWsut4ejSK1cZRVkLm+wxKYoT8AMpNYleuXMGqVavw6NEj2NjYYNSoUWjVqhXMzc3B5XLBMAzS09MRFxeHR48e4fr16+jfvz9atmyJRYsWoXv37pXuPCcnB7a2thg+fDgmTpzIan/79i1cXFwwbNgwnD17FlwuFy9fvoS2tvb3HTVRaImPExA0/QLS36Sx2kw6mcHq98ZoYG8ig8gIIdWNk56eLrKcq169ehg1ahTGjh0r9l3PixcvsGvXLhw8eBAfPnyoUiD169fHmjVrhO7Exo0bBw6Hgx07dlRpWzVNTEwMPeNBafFG+Lb7CNsQCn4xX6hNSVUJHec6ovW4toh9FUvnqwro76tq6HxVjbTPV7kl9lFRUfDx8alSt13jxo2xZs0aREZGfndgfD4fgYGBaNy4MQYOHAgLCwt07doVJ0+e/O5tE8WT9TETJ0ccw+01wawEpmdRB0NPj4Td+HbgKNGzWUJ+JOXeiVW3r+/EkpKS0LhxY2hpaWHBggVwcnLCrVu3sGTJEhw6dAguLi4it0MjTNc8CTc/IMr3MYqzi1htpr3NYTOxKZQ1qNCWkJqosru4Kv3Lz8jIQGJiIvLy8qCpqQljY2PUrl37uwIsD59ferXt5uaG33//HQDQokULPHr0CDt27Cg3idXE2/wftfuiMKcQN5dcQ/TxJ6w2DT1N9FjjAgtnS1bbj3q+vhWdr6qh81U10j5flSax1NRU+Pv749SpU3j//j2rvUGDBvDw8MDUqVNhYGAgscD09fWhoqLC6s60tramLsUfQOKjBAROD0DG23RWm6mjGZzX94I2r1b1B0YIkSsVJrE3b97A3d0dnz59QufOnTFkyBAYGxtDXV0dBQUFSExMRHh4ODZv3ozjx4/j/PnzaNSokUQCU1NTg52dHat7MDY2FiYmVHlWU/FL+Hiw9R7u+t5mPftSVlNGBy9HtB7bhp59EUIAVJLE5s+fDyUlJdy5c6fC28GYmBh4eHjgzz//xOHDh8XeeXZ2Nl6/fg2gtPswPj4ekZGR0NPTg4mJCaZNm4YxY8agQ4cOcHJyQnBwME6ePImDBw+KvQ+iODI/ZOLSzAv4cDee1VbHsg5c/XvDsKmRDCIjhMirCgcADg0NxZQpUyrtz7SyssLkyZMRGhpapZ1HRETAyckJTk5OyMvLg7e3N5ycnPDXX38BAHr37o2NGzdi06ZN6NChA7Zv345t27aV+zyMKK6X55/jkOtekQmsxahWGHZ+FCUwQghLhXdiysrKKC4uFmtDJSUlVR56ytHREenp6RWuM3LkSNYoHqTmKMwuxI0lV/HsxFNWm2YdTfRY64pGPSxkEBkhRBFUmMS6du2Kv//+Gx07doSdnV256z18+BCbNm1Ct27dJB4gqbkSI/5XvBGXzmozdTKH8zpXKt4ghFSowiT2119/oXfv3ujRowdatGiB1q1bg8fjCQo7kpKSEBERgcjISJibmwu6AQmpCL+Ejwdb7iHMNxRMifBrispqyug4zwmtxthR8QYhpFIVJrG6devi1q1b2LlzJ86ePYvDhw+joKBA0K6urg5bW1ssWbIEv/76K2rVoqtmUrHM+AwEzbyIj/dEFG9Y6cN1U28Y2hjKIDJCiCKq9D0xbW1tTJ8+HdOnTwfDMEhLSxO87Kynp0dTsBCxvTj7HNf+vIzCzAJWW8vRrdFpgRNUNFRlEBkhRFFVacQODoeDOnXqSCsWUkMVZBXg5pKrePZfNKtNU18TPde6omF3Kt4ghFSdWEmMz+fjzp07iIyMZA071bx5czg4OEBZWVnasRIFlPDwI4KmByDjXQarzaxLQ/Rc5wptQ5pahxDybSpNYidOnMDixYuRmJjImhATKL074/F4WLZsGYYMGSKVIIni4Zfwcf/vu7jrd5tdvKGujE7zO6PlL62pO5oQ8l0qTGLHjh3DhAkT4ODggBUrVqBNmzbg8XjQ0NBAfn4+kpKScP/+fezevVswqSUlMpL5PgNBMy/g4332nHL61qXFGwZNqHiDEPL9KpyKpUOHDjA1NcWRI0cq3dDQoUPx/v173L59W6IBEsUaNfvFmWelxRtZhay2lr+0Rqf5Uize4JdA+WEIMsLvQNfJBfwmrQC606uUIv19yQM6X1Uj01HsX79+LbjDqoy7uzvmzp0rkaCI4inIKsCNRVfx/JSI4g0DrdLijW6SGRxapLwcaK6dA6V3MdAsKgLz8CZKGrdE/vSVgDLNNUZITVXh2InGxsaIiIgQa0MPHz6EsbGxRIIiiuXjgw841GufyARm3rUhRgaOlm4CA6B2eCuUX0WDU1Q6cSYnPw/KUfehcvWMVPdLCJGtCpPYmDFjsGfPHsyZMwfR0ewvKACIjo7G7NmzsW/fPowZM0YqQRL5xC/mI2zjbZwYcgSZ74WrD5XVldFlWTf0/XdAtVQfKse/Zi3j8Eug8uS+1PdNCJGdCvtZpk2bhvT0dGzevBm7du2Curo6DA0NBcNOJScno6CgACoqKpg6dSqmT59eXXETGct8n4HA6QFICP/IatNvYgBXf3cYNK6+4g1GVU30cnWNaouBEFL9KkxiHA4HS5Yswfjx4xEQEICoqCgkJiYiPz8fGhoa6Nq1K5o3bw43NzfUq1evumImMvb8VDSuL7oisnij1Vg7dPRygopG9T6HKu7QE8pxL8HJyxUs4+vqoajXsGqNgxBSvcT6pqlbty7GjRsn7ViInCvILMD1RVfw4vQzVpuWoRZ6rusF8y4NZRAZUNzZHZzMdKiEXUFJRjqU9Q1R5DwI/EZNZBIPIaR6UNkWEcvH+x8QNCMAmfGZrDbzbo3Qc60LtAxkO/JGUZ+RKHIfjtjn0bC0aUrl9YT8AMRKYleuXMGpU6cEw06VdSeWDTvl4eGBnj17SjtWIgP8Yj7ubbqDe/5hYPhfj7yhAsc/O6OFZyv5GXlDSQmMqjolMEJ+EBUmsZycHIwaNQo3btyAtrY2mjVrhk6dOgkKOxITE3Hu3DkcOXIEnTt3xv79+2k6lhok4106gqZfQMJDdvGGQRNDuG7oCX3bupQwCCEyU2ESW7p0Ke7cuYONGzdi+PDhUFVlj7RQVFSEw4cPw8vLC8uWLcPatWulFiypHgzD4PmpaNxYdBWF2ezijTY9dNDVKBiqO8+D0eGiqLsHiju7ySBSQsiPrsIkdubMGfz+++/w9PQsdx1VVVV4enri3bt32Lt3LyUxBVeQkY9rC6/g5dnnrDYtQ224/lYXjSN2QOljTunC1GRwjm8H37g++I1bVnO0hJAfXYUvO2dnZ4tdOl+/fn3k5ORIJCgiGx/ux+Og2z6RCaxh90YYGTQaVhkhUMoT/j0rZWVALfB4dYVJCCECFSaxli1bYs+ePcjOzq5wI9nZ2fj333/RsiVdiSsifjEfd9aH4L8hR5H1VfWhsroKuq7sgT67PKClrwVOEXtWZgBAYTnLCSFEiip9Jta/f3+0bdsWgwcPRqtWrcDj8QSFHUlJSYiIiMCJEyeQkZGB06dPV1PYRFLS49IRND0AiREJrDYDW0O4+rlD39pAsKzEzArKL6OE1mPAQUkTuoAhhFS/CpNY+/btERQUhBUrVmDbtm0oLi4WKqVmGAYqKiro0qULFi5cSHdiCoRhGDw/WTryRlFOEau99bg26DDXESrqwn8ihYN/g/LbGCi9fQFOUSEYdQ2UWDZFUa+h1RU6IYQIVPqeWIsWLXD8+HFkZWUhOjqa9Z6YjY0NdHV1qyNWIiEFGfm49udlvDz3gtWmZagN5w29YOZkLvrD6prIW+AH5YjbUI59gmLbNuA3a0tl9oQQmRB7xA4dHR20b99emrGQavDh7nsEzbyArA9ZrLZGPS3QY40LNOtoVbwRJSWUtOmEkjadpBQlIYSIp8LCDmkLDQ3FsGHDYGNjAy6Xi4MHD5a77owZM8DlcrFp06ZqjLDmKCkqwe21Ifhv2DFWAlPRUEHXVT3Qe0f/yhMYIYTIEYklsXXr1kFfX79Kn8nJyYGtrS1Wr14NTU3Nctc7c+YMwsPDUbdu3e8NU7FkZ0Dt0GY0PLoJaoc2A9nscQvFkf42DccHHcb9v9lDRxnaGmH4+VFo8bMcDR1FCCFikugAwAzDVL7SF5ydneHs7AwAmDx5ssh13r17h3nz5uH06dMYNGjQd8eoKDgZqdBYPRPKH+OgBgAxkVCOuof8eb5gatcRaxsMw+DZiae4seSqyOINu9/awmFOJ1bxBiGEKIoKv71OnTol9obKm/n5exQXF2PcuHGYPXs2GjduLPHtyzO14zug/DFOaJnyxziontiJwl/nVvr5/Ix8XJt/CTEBL1lt2kbacPZ1g2knM4nFSwghslBhEhs7diw4HI7Yd1iS7o7y9vZGnTp18Ouvv0p0u4qAk8x+bwsAlMpZ/qX4sPcImnEB2QkiijecLdHDx5mefRFCaoQKk5iOjg5atGiBefPmVbqhI0eO4NChQxILLDg4GIcOHUJwcHCVPhcTEyOxGGTJHErQE7E8C0p4W84x8ov5iNn7HK+OxABfXXcoqSvDdlIzmLibIT7lA5Ai+ZjlSU35O6gudL6qhs5X1XzP+bKysqqwvcIk1rp1a8THx6NTp8pLqcPCwqoWWSVCQkKQmJgo1I1YUlKCJUuWYOvWreV2X1Z2wIqCM3IK+L7zoJSaLFjGr2MEtZ+nwKpBI9b6aW/SEPRHAJIeJ7LaDJsawdXfHXUsq1Z4o6hiYmJqzN9BdaDzVTV0vqpG2uer0iR269YtpKeng8vlVrghhmGqXNhRkXHjxqFfv35CywYOHIiBAwdi9OjREtuPvGJMLZD/+zKond6LguQkqBvyUOjxC5ivEhjDMIg+9gQ3l15DUS67eKPNhHaw/6MjFW8QQmqkCr/ZJk2ahB49eoicR+xrc+bMwZw5c6q08+zsbLx+/RoAwOfzER8fj8jISOjp6cHExASGhobCwaqogMfj/TBXQXwLW+T/4VPulUx+eh6uzr+M2Asiijd4teC8oRcVbxBCarQKkxiPxwOPx5PaziMiItCnTx/Bz97e3vD29sbw4cOxdetWqe23Joi/8w5BMy+KLN6wcLFCdx9naOqV/+4dIYTUBDLtY3J0dER6errY60dFRVW+Ug1XUliCOxtCEb7tHqt4Q0VTBZ2XdEPTYc3pxWVCyA+BHpQokLTXqQicFoBPUUmsNqNmPLj6u0PPQrwXoQkhpCagJKYAGIbBk8ORuLnsGorzioUbOaXFGw5/dIKymrJsAiSEEBmhJCbn8tLy8HDZfSSFsF9yrmVcC84b3GDS0VQGkRFCiOxREpNj70Pf4dKsC8hOzGa1WbhaoftqKt4ghPzYKInJoZLCEtxZH4Lw7fdFF28s7Y6mQ5tR8QYh5If3TUksOzsb6enpIl9uNjEx+e6gfmRpr/5XvPFERPFG8/8VbzSi4g1CCAGqkMTy8/Ph4+OD/fv3IzU1tdz1Kmoj5WMYBk+PRJVbvNF20k+wn9mRijcIIeQLYiexP/74A4cPH4a7uzscHBwqHYaKiC8vLQ9XvS7hVRB7kMxaxrXQdHYL2A/uIIPICCFEvomdxM6dOwdPT09s3LhRiuH8eN6FxOHSrIvISWIXb1i6WaO7d0+8T46XQWSEECL/xE5iHA4HLVu2lGYsMqP0+hnU/tsNTlYGmFo6KOw7CvwmraS6z+KCYtxZF4KH/zxgtalqqaLzsm6wHfy/4o1kERsghBAifhJzc3PDjRs3MGbMGGnGU+04H+Og8fcSKKV8EixT+hCH/KnLwbe0lco+U2NTEDgtAMlPP7HaeC2N4eLnDr2GomYTI4QQ8iUlcVf8448/8ObNG0ybNg0PHjxAYmIikpOTWf8pGrUze4USGAAopX+G2rkDEt8XwzCIOvgYh933sxMYB2g3pT0G/zecEhghhIhJ7Duxdu3aASgdhPfAgfK/4BWtOpGTnSl6eS77GdX3yEvNxRWvS3h9KZbVVqueDlx83dDAnl5PIISQqhA7ic2dO7dGvlzLNzYBnrCfS/ENJDcFTVzwW1yedRE5n3JYbVa9G6PbXz2hUVtDYvsjhJAfhdhJbP78+dKMQ2YKPcZA+UUklN+/EiwrqWeGwiETvnvbxQXFuL0mGBE7w1ltqlqq6LK8O2wGNa2RFweEEFIdaNipWrrI+9MfqgGHofzhLfi8+ijsPRKopftdm02N+V/xRrSI4o1WxnD1cwfXnJ59EULI9yg3iR0+fBgAMGzYMHA4HMHPlRk+fLhkIqtOmtooGjQORRLYFMMwiDrwGLdW3EBJAXvkjXZT7NF+hgOUVWnkDUII+V7lJrHJkyeDw+Fg4MCBUFNTw+TJkyvdGIfDUcwkJiG5Kbm46hWE15dfsdp06pcWb9RvT8UbhBAiKeUmscePHwMA1NTUhH4mosXdeotLsy4iN5ldvGHdpzG6reoJdSreIIQQiSo3iZmamlb4MylVXFCM2z7BiNglonhDWxVdV/RAkwG2VLxBCCFSQIUd3yHl5WcETgvA52fsl7yNW9eFi587uGbc6g+MEEJ+EOWO2DFx4kS8fPmyyht8+fIlJk6c+F1ByTuGYfB4XwQO9z7ASmAcJQ5+mmaPQceHUQIjhBApK/dOLC0tDQ4ODnBwcMCAAQPQpUsXNGrUSOS6r1+/xvXr13Hq1CmEhYWhZ8+eUgtY1nI/5+DK3CC8ufqa1abTQBcuG91Qv10DGURGCCE/nnKT2NGjR3H37l34+/vDy8sLJSUl0NHRgZmZGbhcLhiGQXp6Ot69e4esrCyoqKjA1dUVFy9eFAxRVdPE3XyDS39cRG5yLqvNum8TdFvZg4o3CCGkGlX4TKx9+/Y4ePAgPn/+jMDAQNy/fx8vX77Ep0+lL/DWqVMHAwYMQPv27dGzZ0/o6+tXS9DVrTi/GKFrgvFIRPGGWi01dF3RA409bKh4gxBCqplYhR0GBgb4+eef8fPPP0s7HrlTUfFGXbt6cPFzQ21TbvUHRgghhKoTy8MwDCL3RiD4r1uskTfKijd+muoAJRWxZ7MhhBAiYTL9Bg4NDcWwYcNgY2MDLpeLgwcPCtqKioqwZMkSdOjQAfXq1UPjxo0xbtw4vH//Xupx5X7OwdkxJ3FjyTVWAtNtoItBx4bBfmZHSmBEbMV8BhmFfDAMI+tQCKlRZPotnJOTA1tbW6xevRqamppCbbm5uXj8+DFmz56Nmzdv4tChQ/jw4QMGDRqE4uLicrb4/d5ef40DLnvx9vobVlvj/jYYcXE06rWrL7X9k5qFYRisCM9Ah9Of0P5kEjqfTcbxV+zCIELIt5Fpd6KzszOcnZ0BgDU2Y+3atXH69GmhZb6+vrC3t8eLFy/QtGlTicZSnF+MEO+beLwngtWmplNavNHEw1ai+yQ135bobGx9mo3cktKfE/P4+PNeBmz1VNC0jppsgyOkBlCo/rCsrCwAAJfLlfi2n59+JjKB1W1TDyMueFICI9/k7Nt8QQIr8ymfD/8oyc4cTsiPSmEKOwoLC7Fw4UK4urqifv3yu/NiYmK+aftqrdVh0MYQn8NLqxA5ShxYjmoMixFW+FSQjE8x7OrE6vStx/WjkpfzlZmrDoA97U5iehZiYlKqP6ByyMv5UhR0vqrme86XlZVVhe1VTmKZmZkIDw9HcnIyunTpAiMjo28OTlzFxcUYP348MjIyKp3XrLIDrki9rfVw0GUv1GqpwdXPHXXb1PvmbUlSTEzMdx3Xj0aezlerhFQ8y84TWqamBAyzNYCVlbaMohImT+dLEdD5qhppn68qdSeuX78eNjY2GDBgACZOnIhnz54BAFJSUlC3bl3s3r1b4gEWFxfj119/xdOnT3HmzBnUqVNH4vsoU4tXC/32DsCIi55yk8CIYvNuz0VbQ1Vo/O9fWm1VDnqZamCopZZsAyOkhhA7ie3evRsrV67EoEGD8O+//wqVCuvr68PNzY1ViPG9ioqKMGbMGDx9+hTnzp0Dj8eT6PZFMW5ZF+o66lLfD/kxcNWVcMndELu71sG8Vjo45WKAvV31oUSjuxAiEWJ3J27fvh39+/eHn58fUlNTWe0tWrTA1q1bq7Tz7OxsvH5dOpAun89HfHw8IiMjoaenh7p162L06NGIiIjA4cOHweFwkJSUBADQ1dVlleQTIq+UOBy4mWrCzZT+ZgmRNLHvxN6+fYvOnTuX287lcpGWllalnUdERMDJyQlOTk7Iy8uDt7c3nJyc8Ndff+HDhw+4cOECEhIS0KVLFzRu3Fjw38mTJ6u0H0IIITWT2HdiXC4XycnlV+g9e/asyt19jo6OSE9PL7e9ojZCCCFE7DsxZ2dn7N27V+Td1pMnT7Bv3z64ublJNDhCCCGkImInsYULFwIAHBwcsHTpUnA4HBw8eBBjx45F9+7dwePxMHfuXKkFSgghhHxN7CTG4/Fw48YNuLi44Ny5c2AYBsePH8eVK1cwePBgXL58Warl74QQQsjXqvSys4GBAfz8/ODn54fPnz+Dz+fDwMAASkoKNXoVIYSQGuKbh50yMDAAACQmJiIjIwONGzeWWFCEEEKIOMS+hdqzZw9rpPk5c+bA1tYWDg4OcHJyQkqK/IwFRwghpOYTO4nt2rULWlr/P1ROcHAwdu7ciUGDBmHx4sV4/fo11q1bJ5UgCSGEEFHE7k6Mi4vD6NGjBT+fOnUK9evXx7Zt26CkpISMjAycOnUK3t7eUgmUEEII+ZrYd2IlJSVQVVUV/Hz9+nX06NFDUNTRqFEjJCYmSj5CQgiRE0lp8fiQGov8QpqdW16IncTMzMxw8+ZNAKXDRb19+xbdunUTtH/69Ak6OjqSj5AQQmQsvzAXuwO9sfPCKlx9dgRbzi5GcFSArMMiqEJ34tixYzFnzhw8f/4cHz9+RP369eHs7CxoDwsLQ5MmTaQSJCGEyNKZ23vwJvG54Oe07GSEPL0Ia5NW4HHLn6SXSJ/Yd2Ljxo2Dn58fGjVqBDc3N5w8eRIaGhoAgLS0NCQnJ2Pw4MFSC5QQQmQlMe0da1lufhbuRl+WQTTkS1V6T8zT0xOenp6s5Xp6erhx44akYiKEELnCAQdqz3XAqDAosswWLKeBHmSPfgOEEFKBorwiaJ2rC93D5tA52QCc7NJr/1qaXHRs2kvG0ZEq3Yl9+vQJ+/fvx6NHj5CZmQk+ny/UzuFwcPbsWYkGSAghspL89BMCp51HTmwhAEApRxW6p02gOjkfHZu6Qk/HUMYRErGTWHR0NHr37o3c3FxYWloiOjoaTZo0QXp6OhISEtCwYUPUr08POAkhio/hM4jY+QCha4LBLxK+WFd7rYNBDX9Dvcb0fScPxO5OXLZsGTQ0NHD37l2cOXMGDMPA29sb0dHR2LFjB9LT07FixQppxkoIIVKXnZSNU6NOIHjVTVYCq23OhYOfI+q1pgQmL8ROYmFhYfjll19gZmYmeJjJMAwAYNCgQRgwYAAWLVoknSgJIaQavAqKwUGXvXgfEsdqsx3SDCMueILbRE8GkZHyiN2dWFRUBGNjYwAQlNZnZGQI2ps3b44jR45IODxCCJG+otxC3FpxA08ORbLa1HXV0X21M6zcaaYOeST2nZiJiQni4+MBAJqamjA2Nsa9e/cE7dHR0dDW1pZ8hIQoOIZhUJhTKOswSDk+RSXhcO8DIhNYffsGGBk0mhKYHBP7TszR0REBAQFYsGABAGDw4MHYsmWLoErx6NGjGDVqlNQCJUQR5abk4qpXEIrzi9F/3yBwlDiyDon8D8Nn8PCf+7i9LoT17EtJRQn2f3REmwntoKRMbyLJM7GT2PTp0+Ho6IiCggKoq6vjzz//RHp6Os6cOQNlZWUMHTqUCjsI+ULcrbe4NOsicpNzAAARu8Jh91tbGUdFACA7MQuXZl3E+1D2SBzchnpw8XODccu6MoiMVJXYSczExAQmJiaCn9XV1eHv7w9/f3+pBEaIoiouKMZtn2BE7AoXWn57TTBMOpjCsKmRjCIjABAbGIOrXkHIT89ntTUd2hxOS7pCTVtNBpGRb1Gll50JIRVLefkZgdMD8Dk6mdVmYGsI1Vr05SgrRbmFuLX8Bp4cFlG8UVujtHjDzVoGkZHvUaUklp6ejhMnTuDt27dIT08XlNiX4XA4+PvvvyUaICGKgGEYRO5/hOCVN1FSUCzUxlHioN2U9vhpugOUVZVlFOGPLSkqEYHTApD+Oo3V1sDBBM6+btCpS1NJKSKxk9jVq1cxevRo5OTkQEdHB1wul7UOh0MPrcmPJzclF1fmBuHNlVesNp36OnDZ6I76PzWQQWSE4TMI/+c+7pRTvOEwpxPsfmtLxRsKTOwktnDhQhgZGWH//v1o2rSpRHYeGhqKTZs24fHjx0hISMDmzZsxcuRIQTvDMFi9ejX27t2L9PR0tGnTBuvWrYONjY1E9k/I94q7+QaX/riI3GT2TL/WfZug28oeUK+tIYPISFZCFi7NvID4O+9ZbdxGenD1dwevubEMIiOSJPblx+vXrzFhwgSJJTAAyMnJga2tLVavXg1NTU1Wu5+fHzZv3gwfHx9cu3YNhoaG8PDwQFZWlsRiIORbFOcX4+by6zjt+R8rganVUoOzby+4+rtTApORmAsvcdBlr8gE1mx4C4wIGEUJrIYQ+07MwsIC2dnZla9YBc7OzoLZoSdPnizUxjAMtm7dihkzZqBfv34AgK1bt8LKygonTpzAmDFjJBoLIeJKefkZgdMC8PkZu3jDuHVduPq7o7Ypt/oDIyjMKcStZdfx9GgUq02DW1q8YdmLijdqErHvxP7880/s3r0bb9++lWI4/y8uLg5JSUno1q2bYJmmpiY6dOiAu3fvVksMhHyJYRg83vMQh3sfYCUwjhIHP013wOATwymByUhSZCIOu+8XmcAadDDFyKDRlMBqoHLvxP744w/WMj09PbRv3x5OTk6oX78+lJWFK604HA7WrVsnkcCSkpIAAIaGwvP1GBoaIiEhQSL7IERcuZ9zcHlOEN5ee81q02mgC5eNbqjfjoo3ZIFfwkf49vsIWx8KfvFXxRuqSugwuxPsxrej0VJqqHKT2O7du8v90JUrV0Qul2QS+1YxMTEy3b+01NTjkhZJnq/ke0l4vDYChWkFrLa6Xeuj2fSWyK2Vp9C/I0WNPe9THh77PETq48+sNm2TWmi1oA10rbiIfRUr0f0q6vmSle85X1ZWVhW2l5vE0tLY71NUJx6PBwBITk4WGikkOTkZRkblj3hQ2QEropiYmBp5XNIiqfNVnF+M0NW38Ojfh6w2tVpq6LqiB5oMsP3u/ciaov59xQS8wO35t1CQwR55o9nIlnBa2BmqWpJ/uVxRz5esSPt8ye2IHWZmZuDxeLh+/Trs7OwAAPn5+bhz5w6WL18u4+hITff5RTICpwYg5QX7Cr+uXT24+LnRsy8ZKcwpxM0l1xB9/AmrTUNPEz18nGHhQknmRyF2EgsLC8Pt27cxa9Yske2+vr7o2LEjfvrpJ7F3np2djdevS58x8Pl8xMfHIzIyEnp6ejAxMcGkSZOwYcMGWFlZwdLSEuvWrYO2tjYGDRok9j4IqYrS4o0IhHjfRElBiVBbWfHGT7/bQ0mFXo6VhcRHCQicHoCMt+msNpNOZnDe0Au1eLWqPzAiM2InMR8fH5GjdJR58uQJQkJC8N9//4m984iICPTp00fws7e3N7y9vTF8+HBs3boV06dPR15eHubMmSN42fnkyZPQ0aHhYYjk5STn4PLsQMTdeMNq022gCxc/d9RrS9PSywK/hI8HW+/hru9tkcUbHec6ovW4tlS88QMSO4lFRkZi9uzZ5ba3a9euykUdjo6OSE9PL7edw+Fg/vz5mD9/fpW2S0hVvbn6CpfnBCIvJY/V1ri/Dbqu6AF1XXUZRPZjKiouREFRHrQ1dJGdkIWgGRfw4W48az09izpw9XeHUTOeDKIk8kDsJJabm1vp2IiSfhmaEGkrzi9CyF+38HhvBKtNTed/xRseil+8oSj4DB/nw/bj1ccnKCwugPYzA3CO66I4u5i1bvOfW8JxYReoaqrKIFIiL8ROYpaWlrh27RomTpwosv3KlSto1KiRxAIjRNo+P09G4NTzSHmZwmqr26YeXDZS8UZ1uxx+HOEvb4LJZ6B9oS5KHmkBEE5gGnqa6LHGBRbOlrIJksgVsZ9Oe3p64vLly5g7d65Q+X1qairmzJmDa9euYdSoUVIJkhBJYhgGEbvDcaTvAVYC4yhzYD+rAwYdG0YJTAZefXwKpfdq4G6zhMajOqx2U0cz/Bw0mhIYERD7Tuy3335DVFQUduzYgZ07dwre1fr06RMYhsGIESMwadIkqQVKiCTkfMrB5dkXEXfzLatN16Q2XP3cUbdNveoPjIBfwkdxoApqX7QAh//VowsVwHF+F7Qe24aKN4iQKr0n5u/vj8GDB+Ps2bOCMRTNzc3Rr18/dOrUSRrxESIxr6++wpVyijeaDLBFl+Xdoa5DxRuykBmfgaCZF8Hc08LXKarEsABd1nWHXZe2MomNyDexklheXh5mz54NZ2dn9OvXD46OjtKOixCJKc4vQvCqm4jc94jVpqarjm4re6BxP5qjTlZenH2Oa39eRmEme1ivYvtMWExsCLsu9jKIjCgCsZKYpqYmTp8+jfbt20s7HkIkKjn6EwKnBSA1hl28Ua9dfbhsdINug9oyiIwUZhfixuIrePZfNKtNtbYqTKcbw35QZxjUpnm/SPnE7k5s3bo1oqLYUxwQIo8Yfmnxxm2fYJQUfjXyhjIH9jM6oO2U9jQtvYwkPPyIoOkByHiXwWozdTKH8zpXaNPIG0QMYicxb29vDBo0CNbW1hg9ejTU1CQ/sCYhkpCfko/Ty//Du1tvWW21TWvDxc8dde2oeEMW+CV83P/7Lu763QZTwgi1Kaspo+M8J7QaY0fFG0RsnPT0dKby1YD27dsjPT0dycnJUFFRgbGxMTQ1NYU3xuEgLCxMKoH+yGjUbPHEZRbBf/dT8P4JhnoO+/mKzaCm6Ly0GxVvfKW6/r4y32cgaOYFfLz/gdVWx0ofrpt6w9DGUMQn5Qv9exRPcdoTFL//DzlZadCp9xNUzQaDoyT5F9PFvhMzMDCAoaEh/fKIXHryMRd/z7oCmzsvWW1quuro/ldPWPdpIoPICAC8OPOstHgjq5DV1nJ0a3Ra4AQVDRp5o6Yo+ngJha92AkWZ0ABQ9OYFStKjoNHqr0pHfqoqsZNYQECARHdMiKQkP/2E87+egU0C+/lKoY0xxuzqC936ujKIjBRkFeDGoqt4fopdvKGpr4mea13RsLuFDCIj0sIwDIrizwJFmV8uBT8jGiWpD6Ci306i+5Pb+cQIqQzDZxCxKxy31wRD86vijRIOB9e7NofmsFaYQwlMJhLCPyJwegAy37MvLsy6NETPda7QNtSWQWREqvhFQBH7dw5+AUpSH8kuiYWGhoq1XseOHb85GELElZOUjUt/XMS74DhWW4peLfw30AEfGhhgrA51UVU3fjEf9zeH4a7fHXbxhroyOs3vjJa/tJZ4txKRE0qqgGotoCCZtVy5tuTfxxQ7ifXu3VusP7rU1NTvCoiQyry+HIvLc4KQn8YeeeNhq4a42KsNCtVVYamrDK9WdBdWnTLfZyBwxgUkPGAXb+hblxZvGDSR/+IN8u04HA5U6rqg6PV+oCRHsFxJxwrKhh0kvj+xk9i5c+dYy0pKSvDu3Tvs3bsXfD4fS5YskWhwhHypKK8IwStvIOrAY1abuq46Gi3ohjBDY9im58DCQANL2uiCp6Usg0h/TM9PReP6oiuiizd+aY1O86l440ehZtIfSur6KPpwAXk5mdA2agq1RmPA4Uj+vUyxk1hFYyOOHDkSvXr1QkhICDp37iyRwAj50qcnSQicFoC0V+w7/frtG8DZ1w269XXhDCAmJhVWVqbVH+QPqiCzANcXXcGL089YbZoGWqXFG91omqYfjYqRI1SMHBEfE4M6Uqxql0haVFJSwoABA7B//35JbI4QAYbPIHz7fRztf5CVwJRUlNBhriMGHB5C1Ycy8vH+BxzqtVdkAjPv2hAjA0dTAiNSJbHqxLS0NGRkiKhIIeQbZSdl49Ksi3gfwi7eqG3Ohau/O4xb1pVBZIRfzMe9TXdwzz8MDJ9dvOG4oDNajKbiDSJ9Yiex9+/fi1yekZGB27dvY9OmTXBwcJBYYOTH9iooBle8Loks3rAd0gydl3aDmjYNfSYLGe/SETT9AhIefmS16TcxgKu/OwwaU/EGqR5iJ7EWLVqUe1XFMAzatWsHX19fiQVGfkxFuYW4teIGnhyKZLWp19ZAd++esHJvXO1xkVLPT/6veCObXbzRaqwdOno5QUWDXj8l1Ufsv7a///6blcQ4HA64XC4aNmyIJk1oSB/yfT5FJSFwuujijQb2JnD27QWdevTsSxYKMgtwfeEVvDjDfvalZaiFnut6wbxLQxlERn50YiexkSNHSjMO8gNj+Awe/nMft9eFgF/EF2pTUlGCw+yOsBvfjqZNkZEP9+MRNOMCsuIzWW3m3Rqh51oXaBnQyBtENipNYlevXsXWrVsRFxeHOnXqoH///pg0aVJ1xEZ+ANmJWaXFG6HvWG3cRnpw9XcHrzlNiigL/GI+7vrdxv2/74oo3lCB45+d0cKzFRVvEJmqMIldunQJw4YNA8Mw0NfXx5s3b3D//n1kZGRg3rx51RUjqaFiA2Nw1SsI+en5rLamw5qj85KuUNWi4g1ZyHiXjsBpAUiMSGC1GdgYwtXfHfrWBjKIjBBhFfbP+Pv7w8jICKGhoYiNjUVsbCwcHR2xefNmFBSw52siRBxFuYW4Ou8SAiacYSUw9doacNvaFz18XCiByQDDMHh28ikO9donMoG1/rUNhp4eSQmMyI0K78QiIyMxdepU2NraAgC4XC6WLFmCHj16IC4uDtbW1tUSJKk5kqISETgtAOmv01htDRxM4OzrBp26OjKIjBRk5OPawit4efY5q03LUBvO611h1pmKN4h8qTCJZWVlwczMTGiZmZkZGIZBUlKS1JNYSUkJvL29cezYMSQlJYHH42HIkCGYN28eVFSojFeRMHwG4f/cxx1RxRuqSugwuxPsxrejaell5MO9eATNCEDWhyxWW8MeFuixxgVa+loyiIwoIn4JHw+23isd7Nms8vW/R6WZQFRZPQDw+XxRq0vUxo0bsXPnTmzduhW2trZ4+vQpJk2aBDU1NcydO1fq+yeSkZWQhUszLyD+DvuFeSrekJ3wlzfx4PlN5JwugfK1WgAj/G9dWV0FTou6oPnPLal4g4gt80MmLs28gA9346FZRxMOW50A6Q2dWHkSO3jwIMLCwgQ/5+fng8PhYNu2bTh79qzQuhwOB+vWrZNYcPfu3YOrqyt69eoFoPQusFevXggPD5fYPoh0xVx4iavzLqEgg1280WxECzgt6kLPvmTgYUwwgoJOQO2IIVTj2d23BraG6OXfG3Ws9GUQHVFUL88/x7X5l1GQWVozkZeah8g1D9HsRDOp9bJUmsRu3ryJmzdvspYHBgaylkk6idnb22PXrl14+fIlrK2t8fz5cwQHB2PmzJkS2weRjsKcQtxadh1Pj0ax2jS4Guju4wJLVylenpFyMQyDu/uCoXmkAZQK2VPVtB7XBh3mOkJFnbrsiXgKswtxY+lVPDv+lNWWEZOBjHfp4JrrSWXfnPT0dKby1WSDYRisXLkSGzZsgLKyMoqLizF79mwsXLiw3M/ExMRUY4RElPQXaXj0VzhyP+Sw2vRbG6Cllx00DDRlEBkpyirEE79IJNxgT1rJr1UEzdEcdPcYKIPIiKJKf5aKR94PkfuR/e/doI0hWsy1g4a+xjdv36qSaVzk+lLr5MmTOHLkCHbu3IkmTZogKioK8+bNg6mpKTw9PUV+prIDVkQxMTEKcVz8Ej7Ct99H2PpQ8ItFFG/McYTdb22lXryhKOerun24+x6BM64h+yO7eKOgcSZy+39Ex+4j6NxVgv6+SvFL+Hiw5R7CfEPBlHz1MryaMjp4OaL12DaIfRUr1fMl10ls8eLF+P333zFwYOmVYdOmTfH+/Xv4+vqWm8SIbGR9zMSlmRcRH8Yu3tCzqANXP3cYNefJIDJSUlSCu763cX/LXeCrfhdGlY8clwQUtsuAmbE12lg5ySZIolAy4zMQNPMiPt6LZ7XVsawDV//eMGxqVC2xyHUSy83NhbKycJ+9srJytVRGEvHFBLzA1fmXRRdvjGxZWryhSdPSy0L62zQETgtA0uNEVpu+jT5qT1JHkhKD1o37o5VlRygryfVXApEDL84+x7U/L6Mwkz3gRYtRrdDpz87V+u9drv9iXV1dsXHjRpiZmaFJkyaIjIzE5s2bMWzYMFmHRlBavHFzyTVEH3/CatPQ00SPNS6wcLaUQWSEYRhEH3+Cm0uuoSi3iNVuN74tHGZ3goq6CnWPEbEUZhfixuIrePZfNKtNs44meqxzRaPuFtUel1wnsTVr1mDVqlX4448/8PnzZ/B4PIwePZreEZMDiY8TEDgtABlv01ltpo5m6Lm+F2rxalV/YAT5Gfm4Nv8SYgJestq0jbTh7OsG005SfgOV1CgJDz8iaHoAMt5lsNrMOpuj57pe0DaSzUwGcp3EdHR0sHr1aqxevVrWoZCvZCdksxKYspoyOsx1ROtf24hdvJFfzGDvyxw8/FyILnXVMdhCCyo0asc3i7/zDkEzLyI7gV280cjZEj18nKFZh0beIOLhl/Bx/++7uOt3W2TxRsf5Tmj1i51MR9r55iSWmVk6t5CuLk1S+COydLVCs+Et8ORw6QzM3/IwNzW/BP2DPuNpajFKAJx6k4d9Mbk46WwATRVKZFVRUlSCsA2heLD1Hqt4Q0VDBU5LuqLZ8PJnZyfka5nxGQiacQEf77Nfx6hjpQ/XTb1haGMog8iEVSmJvX37Ft7e3rh06RIyMkpvK2vXrg1nZ2d4eXmhUaNGUgmSyCenxV0Qf/c9TDqYwnFh1Ys3Fj3IRGRqseDnQj5wJ6kQ/lFZ8GpNF0fiSnuThsBp5/EpMonVZtSMBxc/N9SxpJE3iPgqKt5oObo1Oi1wgoqGfBRriZ3EIiIi0K9fP+Tl5aFbt26wtCx9YB8bG4uTJ0/i4sWLOHPmDFq3bi21YIl8UdVSw/Bzo6BW69uGjYpJZxccAMD95MLvCeuHwTAMoo89wc2lIoo3OECbCe3g8EcnKKuxR+UgRJSCrALcWHwVz0+KKN7Q10TPta5oKIPijYqIncTmzZuHWrVq4fr167CwED6I2NhY9OnTB/Pnzxc5HBWpub41gQGAloro6ex01ajLqzL56Xm4Ou8SYi+yR6jR5tWCi68bTDqayiAyoqgSwj8icHoAMt+LKN7o0hA917lC21A2xRsVETuJRUZGYs6cOawEBgCWlpYYN26cRMdNJDXfeBttRKYWIrXg/x/iGGsqYVZzmk+sIu9vv8OlmReQnZjNarNwsUJ3H2do6tGwXkQ8/GI+7m8Ow12/O+ziDXVldJrfGS1/aS23z1PFTmJGRkZQUip/ImhlZWXweDQiAxGfm5kmVhTysft5DtIL+DDQVMKsFjpopk+j2otSUliCOxtCEb5NRPGGpgo6L+mGpsOay+2XDZE/me8zEDjjAhIesIs39BsbwNXfvXROMDkmdhKbMmUK/P39MWjQIDRo0ECo7d27d9i5cyemT58u8QBJzTbSShsjreSvi0LepL1KReD0AHyKElG80ZwHVz936FnUkUFkRFE9PxWN64uuoDCL/Qy61Rg7dJznKDfFGxURO4nl5eWhdu3aaNeuHXr16iWoRHz16hUCAwNhYWGB3Nxc+Pn5CT7D4XAwbdo0yUdNyA+CYRg8PRKFm8uuoTivWLiRA7SZ+BMcZnWk4g0itoLMAlxfdAUvTj9jtWkaaMF5nSvMuypOpbnYSWzp0qWC/z916hSr/enTp3j6VHguGUpiNVdBCYNPeSUw1lKGKr2cLBV5aXm46nUJr4LYxRu1jGvB2dcNJh2oeIOI7+ODDwiaHoDM+ExWm3nX0uINTS7Az08GR91AIbqmxU5ijx8/lmYcRIGsCM/Ambd5SCvgQ19DGSMstTCjBRVjSNL70He4NEt08YZlLyt0X+0MDS4VbxDx8Iv5uLfpDu75h4Hhs4s3HBd0RvOfm6Lw2TrkvXgOhl8IjroB1CzGQkW/jYyiFo/YSczUlK74CLD/ZQ62Rmcj9389WykFxfCNzIKtngqcTehL9XuVFJbgzroQhP9zn1W8oaqlis5Lu8F2SDOFuEIm8iHjXTqCpl9AwsOPrDb9Jv8r3mhsiILotShJDha0MUUZKHyxCcrtNoGjKr8XqVUediorKwvBwcF49+4dgNLk5ujoCB0d+T1IIjnHX+cKEliZjCIGu1/kUhL7TqmxKQiafgGfnogo3mjBg6t/b+g1lM4U76Rmen7yf8Ub2SKKN8baoaOXE1Q0VMAwDEoynrPWYfITUfQhAGrm8jtzSJWSmL+/P9asWYPc3FwwzP9fJmppaWHu3LlUnfgDKC5nKreir7ooiPgYhsGTw5G4tew6ivPZxRttJ7eH/cwOUFal4g0inoLMAlxfeAUvzrCLN7QMtdBzXS+Yd2n4xVIGDFPMWhcAUMweTFqeiJ3ENm3ahCVLlqBTp04YN26c0LBTO3bswLJly6CkpISpU6dKLVgiez8ZqeFOUqFQT5cKB+hZX0NmMSmyvNRcXJ13Ca+CYlltterqwMW3Fxo4UFc+Ed+H+/EImnEBWSKKNxp2b4Qea1ygZSD8WguHowQlzfrg53/VC6DKhXI9N2mG+93ETmLbt29H165d8d9//wn1xzdt2hR9+/aFh4cHtm/fTkmshlvQWhdPUotw71MhMosY6Klx4FhXHb/Z0LteVfUuJA6XZl5AzqccVpulmzW6e/ek4g0iNn4xH3f9buP+33dFFG+owHFhZ7QY1arc56lqTWagMGoZ+DnvAKYIHDUDqNR3g7JW/eoI/5uJncRSU1Ph5uYm8gRwOBz07t0bixcvlmhwRP6oKXNwwtkAD5MLEf65EB2N1WCrRyNsVEVxQTHurA3Bwx0PWG2qWqrovKwbbAdT8QYRX3pcOoKmByAxIoHVZmBjCFd/d+hbG1S4DWVNI2i080dJygMwBSlQMXQAR03+n8GKncRatmyJ58/ZD/7KPHv2DC1btpRIUET+2Rmqwc6QkldVpcamIHBqAJKjP7HaeC2N4ervDq65/H9xEPnAMAyen4zGjcVXRRZvtP61DTp4OUJFXbyveg5HGSoG7SUdplSJncTWrl2LgQMHwsTEBL/++itq1Sqdej47Oxs7d+5EQEAA/vvvP6kFSogiYxgGTw5F4tZy0cUb7abYo/0MByreIGIryMjHtYVX8PIs++ZCy1AbzutdYda5oYhP1ixiJ7Fff/0VHA4Hy5Ytw4oVK2BkVDqD76dPn8Dn88Hj8TB27Fihz3A4HISFhUk2YkIUTF5qLq54XcLrSyKKN+rpwHWjG+q3N5FBZERRfbgXj6AZAcj6wK4cbNjdHD3WukFLX0sGkVU/sZOYgYEBDA0NBVWJZWg2Z0LKFxf8FpdmXkRuMrt4w6p3Y3T7qyc0alNlJxFPSVEJ7vndwf3NIoo3VEvQfthLNO2fAM06HjKKsPqJncQCAgKkGQchNUpxQTFurwlGxM5wVpuqtiq6LO8Om4FNqXiDiK2i4o06JlnoOjEKevVywWR9QklSMFSMu1R/kDJQ5RE7CCEVS3n5GYHTA/A5OpnVxmtlDFc/Kt4g4mMYBs9OPMWNJVdRlFPEam/uEoe2A2OhrPq/OzOmCCXpkZTEylNUVISXL18iMzMTfD57+IaOHTtKJDBCFA3DMIjc/wjBK2+ipEC4eIOjxEG7Ke3x03Qq3iDiy8/Ix7UFlxFz/gWrTdtIG04TXqN+o69mOeCoQonbrJoilD2xkxjDMFixYgV27NiBnBx2/36Z1NRUiQRGiCLJTcnFlblBeHPlFatNp74OXDa6o/5PDUR8khDR4sPeI2jmBWR/ZBdvNOppgR5rXKCSH4LC2Dig+P9H51DStYSKUefqDFWmxE5iGzduhK+vL0aPHo0OHTpgwoQJWLZsGWrXro1//vkHKioqWL58uTRjJUQuxd16i0uzRBdvWPdtgm4re0CdijeImEqKSnDX9zbub7nLmslARUMFTou7otmIFv97nuoMjro+iuJPAyUFUNJpBLWGo8FR+nHu9sVOYgcOHEDfvn2xceNGwd1Wy5Yt0blzZwwbNgzdu3dHSEgIOnf+ca4AyI+tOL8YoWuC8WgXu3hDrZYauqzojiYetlS8QcSW/jYNgdMCkPQ4kdVmaGsEV3931LHSF1quot9G7uf8kiYlcVeMj48XJCglpdKPFRQUAADU1dUxdOhQHD58WAohEiJ/Ul5+xtH+B0UmMOPWdTHioidsBlD1IREPwzB4eiwKh3rtE5nA7Ma3xZDTI1gJjFThTozL5SI/Px8AoKurCzU1NXz48EHQrq6uLpXnYYmJiVi6dCkuX76M7OxsmJubY/369ejUqZPE90VIZRiGQeS+RwheVU7xxu/t8dM0Kt4g4svPyMe1+ZcQE/CS1abNqwXnDb1g2slMBpEpBrGTmI2NDaKiogCU3onZ2dlh165dcHZ2Bp/Px549e2BlZSXR4NLT0+Hi4gJ7e3scO3YM+vr6iIuLg6GhoUT3Q4g4cj/n4PKcILy99prVptNAFy4b3VC/HRVvEPHF33mHoJkXkZ3ALt6wcLFEdx8XaOrRTAYVETuJDR48GLt27UJ+fj40NDSwePFieHh4oHnz5gAAVVVVHDp0SKLB+fv7w9jYGNu3bxcsMzc3l+g+CBHH2xtvcHn2ReQm57LaqHiDVFVJYQnCfEPxYOs90cUbS7qi2fAW1B0tBrGT2MiRIzFy5EjBzw4ODggLC8PFixehrKyM7t27w8LCQqLBBQQEoHv37hgzZgyCg4NhbGwMT09P/Pbbb/TLJdWiOL8Yoatv4dG/D1ltarXU0HVFDzQZYCuDyIiiSnudisDpAfgUmcRqM2rGg4ufG+pY0rMvcXHS09Pldl55Ho8HAJg8eTL69++PqKgoeHl5YcmSJRg/frzIz8TExIhcTkhVZb3JxKO/wpH1hj1DLtdWD63mt4FWXZoMlIiHYRjEB75D9OYolOSXCDdygEZDLGH9iw2UVMWut/shVPaY6puTWGZmJubPn4/p06fD2tr6m4KrjKGhIVq3bo1Lly4Jli1fvhznz5/HvXv3pLJPeRQTEyPx54012feeL4ZhELk3AsF/3URJgfCXDUeJg5+mO+Cn3+2hpFIzvmzo76tqvuV85afn4eq8S4i9yL7I1ubVgouvG0w6mkoqRLki7b+vbx47MT8/H4cPH8aQIUOklsR4PB4aN24stMza2hrx8fFS2R8hOck5uDInEG+vv2G16TbQhYufO+q1le/p2ol8eX/7HS7NvIDsxGxWm4WLFbr7OFPxxnf4rgGAGUa6PZH29vaIjRWegyk2NhYmJjT3EpG8N9de4/KcQOR9ZhdvNO5vg64rekBdV10GkRFFVFJYgjsbQhG+TUTxhqYKOi/phqbDmtPz/e8k16PYT548Gc7Ozli3bh0GDBiAyMhI/PPPP1i0aJGsQyM1SHF+EUK8b+HxnghWm5rO/4o3PKh4g4gv7dX/ijeiRBRvNOfB1c8dehZ1ZBBZzfPNSUxNTQ0dO3YEl8uVYDjC7OzscPDgQSxfvhxr165FgwYNsGDBAowbN05q+yQ/ls/PkxE4LQApLz6z2uq2qQeXjW6obcqt/sCIQmIYBk+PROHmsmsozhN+GR4coM3En+AwqyOU1ehleEn55iTG5XJx/vx5ScYikouLC1xcXKS+H/JjYRgGj/59iNDVt9jFG8octJ/ugHZTak7xBpG+vLQ8XPW6hFdB7OKNWsa14OzrBpMONbN4Q5bkujuREGnI+ZSDy7MvIu7mW1abrkltuPq5o26betUfGFFY70LicGnWReQksYs3LHtZoftqZ2hwqXhDGqo0n9iePXuwf/9+vH37Funp6ax1OBwOUlJSJBkfIRL15uqr0uKNlDxWW5MBtuiyvDvUdah4g4inpLAEd9aFIPyf+6ziDVUtVXRe2g22Q5pR8YYUiZ3EFi9ejM2bN6N58+YYMmSIVJ+FESJpxflFCF51E5H7HrHa1HTU0G1VTzTuZ1P9gRGFlRqbgsBpAUh++onVZtSCB1f/3tBrqCeDyH4sYiexw4cPo2/fvtizZ48UwyFE8pKfJSNw6nmkxrB7Ceq1qw8XXzfomtSWQWREETEMg6hDj3Fr2XUU57OLN9pO+gn2M6l4o7qIncTy8/PRpUsXKYZCiGQx/C+KNwpFFW90QLsp7al4g4gtLzUXD5feQ1Ioe86vWnV14OLbCw0cqHijOomdxJycnPDw4UP88ssvUgyHEMnIScrG5TmBIos3apvWhoufO+raUfEGqZobS66JTGCWbtbo7t2TijdkQOxL0PXr1+PBgwdYt24dPn1i9wETIi+S7iTioOtekQnMZlBTDL/gSQmMfJNOCzpDVUdV8LOqlip6rHWB25Y+lMBkROwBgI2NjcEwDIqKigCUzh+mpCScAzkcDj5+/Cj5KH9wNECr+B7ufIDgFTdYy9V01dH9r56w7tOk2mOSd/T3VTXBu2/i4bL74LU0hqu/O7jmVLxREbkZANjDw4PKRIlcKywqQEqDd4AaAxT+/99qvZ8awMW3F3QbUPEG+XZMYQaK3v2HJjavwFvTHo08ekFFTbXyDxKpEjuJbd26VZpxEPJdsvMysOfSGiSlxUPdTQ86pxsASgx+mtEB7X93gJIyFW+Qb1eSFYuCqFVg8hOgBUCLF4Xi5xFQbr6YLu5ljEbsIDXCxXuHkJRWOkVPQas0KH9SR2HTDBR0saYERr5bYcwOMPkJ/7+AX4iS1IcoSbkPFYOfZBcYqTiJhYeHV3mDbdq0+eZgCPlWqdnJ//8DB8h1Ka0ge//plYwiIjUJU8AeIBr8AhR/ukVJTMYqTGI9evQQ+1aZYRhwOBykpqZKJDBCqkJNRfRQUZrq2tUcCamRlEVXHippUpWrrFWYxDZv3lxdcRDyXeyb9ERC6jvkFfz/AKy6Wnro3LKPDKMiNYWKkSOKcuMBfr5gGUezAVRN+skwKgJUksRGjBhRXXEQ8l1szOxQUJyPu8+vIDM7HXV0DdCtlQcMa9OVMvl+qmZDAHBQnByCwrwMqOs2gJrVBHBU6E5f1qiwg9QYrSw6oJVFB3rviUgch8OBmvkQqJkPwQf6+5IrVLZFCCFEYVESI4QQorAoiRFCCFFYlMQIIYQoLEpihBBCFJbYo9gTQggh8obuxAghhCgsSmKEEEIUFiUxQgghCouSGCGEEIVFSYwQQojCoiQmpzZs2ICuXbvCxMQEFhYWGDp0KKKjo2UdlsLYsGEDuFwu5syZI+tQ5FZiYiImTpwICwsL8Hg8tG/fHiEhIbIOSy6VlJRg5cqVaNGiBXg8Hlq0aIGVK1eiuLhY1qHJhdDQUAwbNgw2Njbgcrk4ePCgUDvDMPD29kaTJk1gbGwMd3d3PHv2TCL7piQmp0JCQvDrr78iKCgIZ8+ehYqKCvr374+0tDRZhyb37t+/jz179qBp06ayDkVupaenw8XFBQzD4NixY7h79y7WrFkDQ0NDWYcmlzZu3IidO3fCx8cH9+7dw+rVq7Fjxw5s2LBB1qHJhZycHNja2mL16tXQ1GTPvebn54fNmzfDx8cH165dg6GhITw8PJCVlfXd+6b3xBREdnY2TE1NcfDgQfTq1UvW4citjIwMdO7cGf7+/vDx8YGtrS3Wrl0r67DkzvLlyxEaGoqgoCBZh6IQhg4dCj09PWzbtk2wbOLEiUhLS8PRo0dlGJn8qV+/PtasWYORI0cCKL0La9KkCX777TfMnj0bAJCXlwcrKyusWLECY8aM+a790Z2YgsjOzgafzweXy5V1KHJtxowZ6NevH5ycnGQdilwLCAhAmzZtMGbMGFhaWqJTp074559/wDB0TSuKvb09QkJC8PLlSwDA8+fPERwcjJ49e8o4MvkXFxeHpKQkdOvWTbBMU1MTHTp0wN27d797+zSfmIKYN28emjdvjp9++knWocitvXv34vXr1/jnn39kHYrce/v2LXbt2oXJkydjxowZiIqKgpeXFwBg/PjxMo5O/syYMQPZ2dlo3749lJWVUVxcjNmzZ2PcuHGyDk3uJSUlAQCrq9rQ0BAJCQnfvX1KYgpgwYIFCAsLQ2BgIJSVlWUdjlyKiYnB8uXLERgYCFVVVVmHI/f4fD5at26NJUuWAABatmyJ169fY+fOnZTERDh58iSOHDmCnTt3okmTJoiKisK8efNgamoKT09PWYf3Q6MkJufmz5+PkydP4ty5czA3N5d1OHLr3r17SElJgb29vWBZSUkJbt++jd27d+Pjx49QV1eXYYTyhcfjoXHjxkLLrK2tER8fL6OI5NvixYvx+++/Y+DAgQCApk2b4v379/D19aUkVgkejwcASE5OhomJiWB5cnIyjIyMvnv7lMTkmJeXF06dOoVz587B2tpa1uHINXd3d7Ru3Vpo2ZQpU2BhYYFZs2ZBTU1NRpHJJ3t7e8TGxgoti42NFfqSIf8vNzeX1QuirKwMPp8vo4gUh5mZGXg8Hq5fvw47OzsAQH5+Pu7cuYPly5d/9/Ypicmp2bNn4+jRozhw4AC4XK6gX1lbWxu1atWScXTyh8vlsopetLS0oKenB1tbW9kEJccmT54MZ2dnrFu3DgMGDEBkZCT++ecfLFq0SNahySVXV1ds3LgRZmZmaNKkCSIjI7F582YMGzZM1qHJhezsbLx+/RpAaVd1fHw8IiMjoaenBxMTE0yaNAkbNmyAlZUVLC0tsW7dOmhra2PQoEHfvW8qsZdT5VUhenl5Yf78+dUbjIJyd3enEvsKBAUFYfny5YiNjUWDBg3w22+/YcKECeBwOLIOTe5kZWVh1apVOH/+PD5//gwej4eBAwdi7ty50NDQkHV4MhccHIw+ffqwlg8fPhxbt24FwzBYvXo19uzZg/T0dLRp0wbr1q2TyAUmJTFCCCEKi94TI4QQorAoiRFCCFFYlMQIIYQoLEpihBBCFBYlMUIIIQqLkhghhBCFRUmMyI3mzZsLhvWRJnd3d7i7u0t9Pz+KiIgI6Ovr49WrV4JlkyZNQvPmzWUYlewsWrQI3bt3l3UYPwxKYkTg4MGDgpEvbt++LXKd1q1bg8vlyn0SuHv3Lry9vZGeni7rUESKjo7GpEmT0KxZMxgZGcHU1BS9evXCtm3bkJ+fL/H95ebmwtvbG8HBwRLf9vLly9G3b19YWFhIfNtfOnnyJCZMmIB27dpBT08P7dq1K3fdgoICLF26FDY2NjA2Nka3bt1w9epVkeu+ePECgwYNQoMGDWBubo7ffvsNnz59Ernu/v370b59e/B4PLRu3RpbtmxhTV8zZcoUPHnyBBcuXPj2gyVioyRGWDQ0NHDixAnW8vv37+PNmzcKMULBvXv34OPjg4yMDFbbqVOncOrUKRlEVergwYNwcnLC1atX4eHhgXXr1uHPP/+EsbExFi1aJJgSRZLy8vLg4+ODkJAQiW43KioK169fZ01s6O/vjwcPHkh0X7t27cL58+dhaGhY6cCxkydPxt9//42BAwdi9erVUFVVxZAhQ1jH/+HDB7i5uSE2NhYLFy7EtGnTcPnyZfTv3591MfHvv/9i6tSpsLa2xtq1a9G+fXssWLAA69evF1rP2NgYvXr1wqZNmyRz4KRCNHYiYenZsydOnz4NHx8foWlNjh8/Dmtra4lPB5ObmwstLS2JbrMishwMODw8HNOnT4ednR2OHz+O2rVrC9omTJiAV69eKdRsywcOHICRkRE6deoktFwa0+Fs374ddevWhbKyMtzd3cu9WwoPD8d///2HpUuXYsaMGQCAYcOGwcHBAYsXL8a1a9cE627YsAHZ2dm4fv06TE1NAQB2dnbo378/Dhw4IJgvLC8vDytWrECPHj2wf/9+AICnpyf4fD7Wr1+PMWPGQF9fX7BdDw8PjB49Gq9evZL6HeqPju7ECMugQYOQlpYm1P1SUlKCU6dOlTtg56ZNm+Dq6opGjRqBx+OhQ4cO2LdvH2u9sudeN2/eRPfu3cHj8eDn51duLKdOnYKhoSFmzpwp6LZ5+PAhBg8eDFNTUxgbG8PV1RW3bt0SfMbb21swkG3Lli0FXaRlXWlfPxOLi4sDl8uFr68v9u7di1atWsHIyAhdu3bFw4cPWTGdPn1a0KXk4OCAs2fPiv0MaPXq1WAYBjt27BBKYGUsLCwwefJkwc8lJSVYt24dWrduDSMjIzRr1gyLFy9GXl6e0OcePXqEwYMHw8LCAjweD82aNcP48eORl5eHuLg4wRepj4+P4HxMmjQJQOngrQsXLkSLFi1gZGQECwsLuLu7IzQ0tNLjCQgIgJOTE5SUhL9KRJ0PLpeLmTNn4vz583BwcICRkRHs7e1x5cqVSvcDAA0aNBDrAurMmTNQUlLCL7/8IlimoaGBUaNG4eHDh4iLixMsP3v2LHr27ClIYADQpUsXWFpaCt2tBwcHIzU1FWPHjhXa17hx45CXl4fAwECh5V26dAGHw0FAQIBYx0a+Hd2JEZZ69erBwcEBJ06cgKurKwDgxo0bSE5OxuDBg0V2xW3duhUuLi7o378/OBwOLly4gGnTpqG4uJj1D//169fw9PTE6NGjMWrUKDRo0EBkHIcPH8bvv/+O8ePHw9vbGwAQEhKCgQMHonnz5pgzZw5UVVVx9OhRDBgwAKdOnYKjoyP69OmDV69e4cSJE/jrr78EV8hfz5/1tZMnTyInJwdjxowBh8OBn58fRo0ahUePHgnuLIKCgjBmzBjY2Nhg0aJFyMjIwLRp01CvXr1Kz2teXh5u3LiBDh06wMzMrNL1gdIZhffv348+ffpgypQpiIiIgL+/P549e4Zjx46Bw+Hg8+fP8PDwgL6+PqZPnw4ul4v4+HhcvHgRubm5MDAwwIYNGzBr1iz07t1bMFBrw4YNAQCzZs3C6dOnMW7cODRp0gQZGRl48OABnjx5go4dO5Yb28ePHxEfH49WrVqJdSxAaTdvYGAgxo4di1q1amH79u3w9PTEkydPUKdOHbG3U5HIyEg0bNiQNYh2mzZtBO1mZmb4+PEjkpOTWVP4AKV3Y18mpsjISABgrduqVSsoKSkhMjISI0eOFCyvXbs2GjZsiLCwMEybNk0ix0VEoyRGRBo8eDD+/PNP5OTkQFtbG8eOHUPbtm0FX3xfe/DggVCX4IQJE+Dh4YFNmzaxktibN29w6NAhuLm5lbv/f//9F7NmzcKsWbMEd1UMw2DmzJmwt7fH6dOnBaOtjx07Fk5OTlixYgUuXbqEZs2aoWXLljhx4gTc3d3FThgfPnzAw4cPBV9+lpaWGDFiBK5evSpI5suWLYOxsTGCgoKgo6MDAOjcuTN69+5d6Vxcr169QlFREZo2bSpWPE+ePMH+/fsxYsQIbNmyRbC8QYMG8PHxQVBQEFxdXXH37l2kpaXh5MmTQl+yCxYsEPx/v379MGvWLDRt2hRDhw4V2k9QUBBGjx6Nv/76S6y4ysTExABAlSZrffnyJe7evYtGjRoBABwdHdGpUyecOHFCYjNKJyUlCSZi/FLZssTERMF6Xy7/krGxMTIzMwV//2XrGhsbC62npqaGOnXqCLb5JXNzc7x48eL7DoZUiroTiUj9+/dHUVERAgICkJeXhwsXLmDw4MHlrl+WwIqKipCWloaUlBQ4OjrizZs3rOKK+vXrV5jAtmzZglmzZmHBggVC81tFRUUhJiYGgwYNQmpqKlJSUpCSkoKsrCx06dIFDx48QG5u7jcfc9++fYWu3jt06AAAePv2LQAgISEB0dHRGDJkiCCBAUCnTp3EmlIiKysLAMSeD+7SpUsASqvdvjR58mQoKysL2nV1dQEAgYGBKCoqEmvbX9LV1cWDBw/w8ePHKn0uNTUVQPnTBoni6OgoSGAA0KxZM+jq6grOsSTk5eWJnMW7rCCprGCjrEtW1Lply75cV1VVVeQ0NRoaGiIrSrlcLlJSUr7xKIi46E6MiKSnp4du3brh+PHjUFFRQW5uLgYMGFDu+gEBAVi7di2ioqJQUlIi1JaZmSn0/KeiO6O7d+/i6tWr+P333zFnzhyhtrL3kKZOnYqpU6eK/Hxqauo3F4l83a1Z9uVcVqb//v17ABD6Ei7TqFEjPH78uMLtlyW+7OxsseJ5//49OBwOLC0thZbXrl0bxsbGePfuHYDSJNqvXz/4+Phgy5Yt6NChA9zc3DBo0CBoa2tXup8VK1YIyv1btGiBHj16YOjQobCyshIrzq9LzCsiquu4du3aEn0VQlNTEwUFBazlZYmmLJlpamoCgMh1y5Z9uW5RURH4fD7r+V9+fr7Iil2GYWhutmpASYyUa/DgwZg4caLgTsfQ0FDkenfu3MHPP/8Me3t7+Pr6wtjYGGpqarh06RK2bNnCmsK97MtDFGtra+Tm5uL48eMYM2aMUGVX2XaWLl1a7nMYAwODKh7l/yuvaKAqX9IVadSoEVRVVREdHS2R7ZXhcDjYu3cvwsPDERgYiBs3bmD69OnYsGEDrly5Uu7vrUz//v3h4OCACxcu4Nq1a9i+fTv8/PywZcuWCu++y55hVSUBSfscA6Xdg2UXHF/6ukuwrBuxbPmXEhMToaurK7gI+HLdunXrCtYrLCxEamoqq5sRKD0vX1YsEumg7kRSLjc3N6irqyMsLKzCacTPnDkDDQ0NnDp1Cp6ennB2dkaXLl0qTFbl0dPTw+nTp6GlpYV+/foJfRmVPY+rVasWunTpIvI/ab7DVvbMq2wa9i+JWvY1LS0tdO7cGaGhoYK7qMr2xzAMYmNjhZZnZmYiMTFRqKIOKC1c+PPPP3H58mUcP34ccXFx2Lt3LwBUekfA4/EwZswY7N+/H1FRUTA3NxcU05SnrFDmy2o/edC8eXO8efOGlVzL3lsrq5qsV68eDAwMEBERwdrGw4cP0axZM6FtAmCtGxERAT6fL7IyNS4uDtbW1t91LKRylMRIubS0tLB+/Xp4eXmhd+/e5a6nrKwMDocjdMeVnp6OAwcOfNN+jY2Ncfr0aTAMg/79+wuulFu1aoVGjRph8+bNgudLX/r8+bPg/8uuoCXZTVW3bl3Y2tri2LFjQvsPCQkR++5q7ty5AEoLXzIzM1ntb968wdatWwEAzs7OACD4ucy2bdtQUlICFxcXAKXH+PWdTMuWLQFA8Dyy7ILi6/NRUlLCembJ5XJhZmYm8kXxLxkbG8PMzAyPHj2qcL3q1q9fP/D5fOzZs0ewrKCgAAcPHkSrVq2EClH69u2Ly5cvC11U3Lx5E7Gxsejfv79gmZOTE/T09LB7926hfe3atQsaGhqCwp8yGRkZePPmDdq3by/RYyNs1J1IKjRs2LBK13F1dcXmzZvh4eGBoUOHIi0tDXv37oWRkZHIrhpxmJmZ4fTp03Bzc4OHhwcCAgKgp6eHTZs2YdCgQbC3t8fIkSNRv359JCQkIDQ0FAzD4Pz58wD+vxR6+fLlGDRoENTU1ODk5FRp11plFi1ahBEjRsDFxQUjR45ERkYGduzYAVtbW7Gedf3000/w9fXFzJkz0a5dOwwdOhSWlpbIz8/HvXv3cPbsWYwYMQJAadHDqFGjsH//fmRmZsLJyQmPHz/GgQMH0KNHD0GSO3ToEHbu3InevXujYcOGyMvLw8GDB6GsrIx+/foBKE1iNjY2OHnyJCwtLVGnTh2YmZnB0tIStra26NOnj6DIIiwsDFeuXMFvv/1W6fG4ubnhv//+E/msSNJCQ0MFw6G9f/8eOTk5WLt2LYDSIpyy1wHatm2L/v37Y+XKlUhJSYGFhQWOHDmCuLg41ushZa8X9O3bFxMnTkReXh78/f3RpEkTeHp6CtbT1NTEn3/+idmzZ2PUqFHo2bMnbt++jWPHjmH+/PmsbuwbN26AYZgKC5iIZFASI9/N0dERW7duha+vL+bPn4969eph/Pjx4HK5+P333795u1ZWVjh58iT69OmDAQMG4MyZM+jYsSMuX76MtWvXYteuXcjKyoKRkRHs7OyEvnRat26NJUuWYNeuXZgyZQr4fD7OnTv33UmsV69e2LVrF1avXo3ly5fDwsICmzdvxpEjR/D8+XOxtuHp6Qk7Ozv8/fffOHnyJD59+gQNDQ00bdoUq1atEjqOjRs3wszMDAcOHMDFixdhZGSEqVOnYv78+YIuwo4dOyIiIgKnTp3Cp0+foKOjgxYtWmDNmjVo27atYFubNm2Cl5cXFi5ciIKCAgwfPhx+fn4YN24crl+/josXL6K4uBhmZmaCYo/K/Pzzz9i6dStCQkLg5ORUxbNZNbdu3YKPj4/QslWrVgEAvLy8hN5p27ZtG/766y8cO3YMaWlpsLGxwZEjR1gxNmjQAAEBAVi4cCFWrFgBFRUV9OzZE6tWrWJ1TY8bNw5qamr4+++/ERQUhHr16mHlypWs6lHg/1+IF7c4hnw7Tnp6uuSeqBLyg+rUqRMMDAxw+vRpWYdS7QYMGIDatWvj33//lXUociExMREtW7bErl27KuyGJ5JBz8QIqYKioiIUFxcLLQsODsaTJ09Y4wf+KBYtWoSzZ88KTcXyI9u8eTOaNm1KCaya0J0YIVUQFxeH/v37Y8iQIahbty5evnyJf//9F7q6urhz547Ehk4ihIiHnokRUgVcLhetWrXC/v378fnzZ2hpacHZ2RlLly6lBEaIDNCdGCGEEIVFz8QIIYQoLEpihBBCFBYlMUIIIQqLkhghhBCFRUmMEEKIwqIkRgghRGH9H/1tNoThjybcAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fe = smf.ols(\"purchase ~ mkt_costs + C(city)\", data=toy_panel).fit()\n",
    "\n",
    "fe_toy = toy_panel.assign(y_hat = fe.fittedvalues)\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase, c=toy_panel.city)\n",
    "for city in fe_toy[\"city\"].unique():\n",
    "    plot_df = fe_toy.query(f\"city=='{city}'\")\n",
    "    plt.plot(plot_df.mkt_costs, plot_df.y_hat, c=\"C5\")\n",
    "\n",
    "plt.title(\"Fixed Effect Model\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "잠시 시간을 내서보면, 위의 그래프가 고정 효과가 어떤 역할을 하고 있는지 알려주고 있어요. 고정 효과는 **도시당 하나의 회귀선**을 적합(fitting)한다는 점과 적합된 회귀선이 평행하다는 점을 주목해주세요. 회귀선의 기울기는 인앱 구매에 대한 마케팅 비용의 효과입니다. 따라서 **고정 효과는 인과 효과가 모든 개체(이 예시에서는 도시)**에 걸쳐 일정하다고 가정합니다. 이러한 가정은 관점에 따라 약점이 될 수도 있고 장점이 될 수도 있습니다. 만약, 여러분이 도시별 인과효과를 추정하는데 관심이 있다면 약점이 될 수 있어요. 고정 효과 모델은 해당 효과가 개체(도시)간에 일정하다고 가정하기 때문에, 도시에 따른 인과효과 차이를 찾지 못할 거에요. 그런데, 인앱 구매에 대한 마케팅의 전반적인 효과를 찾으려고 한다면, 해당 데이터의 패널 형태가 고정 효과를 탐색할 수 있게 도움을 줍니다.\n",
    "\n",
    "## Time Effects\n",
    "\n",
    "개별 수준에 대해 고정 효과를 수행한 것처럼, 시간 수준에 대해서도 고정 효과를 설계할 수 있어요. 만약 고정된 개별 특성에 대한 개별 더미변수를 추가한다면, 시간 더미를 추가하는 것은 각 기간 동안 고정되는 변수들에 대해 통제를 할 수 있지만, 시간이 지남에 따라 바뀔 수 있어요. 해당 변수의 한 예시는 인플레이션이에요. 물가와 급여는 시간이 지남에 따라 상승하는 경향이 있으나, 각 기간의 인플레이션은 모든 개체에서 동일합니다. 조금 더 구체적인 예를 들자면, 결혼이 시간이 지남에 따라 증가하고 있다고 가정해 봅시다. 임금과 결혼 비중도 시간이 지남에 따라 변한다면, 시간은 교란 변수로 볼 수 있게됩니다. 인플레이션은 또한 시간이 지남에 따라 급여를 증가시킬 수 있으므로, 결혼과 급여 사이에 긍정적인 연관성이 있는 것은 단순히 시간이 지남에 따라 두 가지 변수가 모두 증가하기 때문일 거에요. 해당 교란효과를 보정하기 위해 각 기간에 대한 더미 변수를 추가할 수 있습니다. 패키지 'linearmodels'에서는 수식에 'TimeEffects'를 추가하고 'cluster_time'을 true로 설정할 수 있어서, 매우 간단하게 해결할 수 있어요."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>-0.0062</td>   <td>0.0008</td>   <td>-8.1479</td> <td>0.0000</td>   <td>-0.0077</td>   <td>-0.0047</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0727</td>    <td>0.0228</td>   <td>3.1858</td>  <td>0.0015</td>   <td>0.0279</td>    <td>0.1174</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.0476</td>    <td>0.0177</td>   <td>2.6906</td>  <td>0.0072</td>   <td>0.0129</td>    <td>0.0823</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>    <td>-0.0001</td>  <td>3.546e-05</td> <td>-3.8258</td> <td>0.0001</td>   <td>-0.0002</td> <td>-6.614e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects+TimeEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True, cluster_time=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "위의 새로운 모델에서는 결혼이 임금에 미치는 영향이 '0.1147'에서 '0.0476'으로 크게 감소했어요. 그래도 해당 결과는 99% 신뢰 수준에서 유의미하므로 남자는 여전히 결혼으로 인한 수입 증가를 기대할 수 있습니다.\n",
    "\n",
    "## When Panel Data Won't Help You\n",
    "\n",
    "패널 데이터와 고정 효과 모델을 활용하는 것은 인과추론을 위한 매우 강력한 도구입니다. 랜덤 데이터나 인과추론을 위한 마땅한 방법이 없을 때, 고정 효과는 비실험 데이터(non-experimental)로 인과추론을 할 때만큼 설득력이 존재합니다. 그래도 언제나 적용될 수 있는 만병 통치약은 아닙니다. 패널 데이터조차 도움이 되지 않는 상황이 있기 마련이에요.\n",
    "\n",
    "가장 분명한 것은 시간이 지남에 따라 변하는 교란변수가 있을 때입니다. 고정 효과는 각 대상의 일정한 속성에 대해서만 편향을 제거할 수 있어요. 예를 들면, 책을 읽고 많고 좋은 지방을 섭취함으로써 지능을 높일 수 있다고 가설을 세워봅시다. 그러면, 이렇게 지능을 높이는 요소가 여러분이 고소득 직업과 아내를 얻을 수 있게해요. 고정 효과는 관측되지 않은 교란변수(지능)로 인해, 해당 편향을 제거할 수 없습니다. 왜냐하면, 해당 예시에서 지능은 시간에 따라 변하기 때문이죠.\n",
    "\n",
    "![img](./data/img/fixed-effects/time-travel.png)\n",
    "\n",
    "고정 효과가 실패하는 또 다른 경우는 **역인과관계**가 있는 경우입니다. 예를 들어, 더 많은 돈을 벌게 하는 것은 결혼이 아니라고 가정해봅시다. 더 많이 버는 것은 여러분이 결혼할 확률을 높이고 있어요. 이 경우 양의 상관관계가 있는 것처럼 보이지만 수입이 먼저임을 알 수 있습니다. 해당 변수들은 시간에 따라 같은 방향으로 변할 것이므로, 고정 효과로는 시간에 변화에 대한 통제를 할 수 없습니다.\n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "이번 챕터에서 우리는 여러 기간에 걸쳐 동일한 대상에 대해 여러 측정값이 있는 데이터인 패널 데이터를 사용하는 방법을 살펴보았습니다. 이 경우 개체(entity)를 제어하는 ​​고정 효과 모델을 사용해, 모든 대상과 시간이 변하지 않는 변수를 고정할 수 있습니다. 고정 효과 모델은 교란변수를 통제하는 ​​강력하고 매우 설득력 있는 방법이며 랜덤이 아닌 데이터를 사용하는 것만큼이나 효과가 좋습니다.\n",
    "\n",
    "마지막으로 우리는 고정 효과 모델이 인과추론의 만병통치약이 아님을 확인해 볼 수 있었습니다. 고정 효과가 잘 작동하지 않는 두 예시를 확인했어요: 역인과관계가 있을 때와 관측되지 않은 교란변수가 시간에 따라 변할 때입니다.\n",
    "\n",
    "\n",
    "## References\n",
    "\n",
    "이 책을 쓰기 위해 Joshua Angrist, Alberto Abadie, Christopher Walters의 대단한 계량 경제학 수업 자료를 많이 참고했습니다. 해당 자료에 있는 대부분의 아이디어는 전미경제학회(American Economic Association)의 수업에서 가져왔어요. 이렇게 좋은 참고자료를 보는 것이 2020년의 힘든 한 해를 지탱하게 만들어준 원동력이었다고 생각해요.\n",
    "\n",
    "* [Cross-Section Econometrics](https://www.aeaweb.org/conference/cont-ed/2017-webcasts)\n",
    "* [Mastering Mostly Harmless Econometrics](https://www.aeaweb.org/conference/cont-ed/2020-webcasts)\n",
    "\n",
    "또한, Angrist의 정말 좋은 책들을 참고자료에 넣고 싶어요. 해당 저자가 쓴 책들은 계량경제학(Econometrics) 또는 '메트릭스(Metrics, 계량적 분석)'가 매우 유용할 뿐만 아니라 매우 흥미롭다는 걸 알려주었어요.\n",
    "\n",
    "* [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com/)\n",
    "* [Mastering ‘Metrics](https://www.masteringmetrics.com/)\n",
    "\n",
    "마지막으로 참고한 자료는 Miguel Hernan과 Jamie Robins의 책입니다. 이 책들은 제가 대답해야 했던 까다로운 인과적인 질문에서 신뢰할 수 있는 동반자 같은 존재였어요.\n",
    "\n",
    "* [Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)\n",
    "\n",
    "마지막으로, Scott Cunningham과 인과 추론과 Mixtape 형태를을 결합한 훌륭한 자료에도 감사를 전합니다.\n",
    "\n",
    "* [Causal Inference: The Mixtape](https://www.scunning.com/mixtape.html)\n",
    "\n",
    "![img](./data/img/poetry.png)\n",
    "\n",
    "## Contribute\n",
    "\n",
    "Causal Inference for the Brave and True는 인과추론, 통계학에 대한 오픈소스 자료입니다. 이 자료는 금전적으로나 모든 분들이 쉽게 접근하실 수 있도록 하는 것이 목표입니다. 또한, 이 책은 Python 기반의 무료 소프트웨어만 사용해요.\n",
    "여러분들께서 이 자료가 가치 있다고 생각하시고, 금전적으로 지원을 원하신다면 [Patreon](https://www.patreon.com/causal_inference_for_the_brave_and_true)를 방문해주세요. \n",
    "만약 여러분이 금전적으로 기여하기가 쉽지 않으시다면, 오타 수정, 수정 제안, 이해하기 난해한 부분에 대한 피드백 제공 등을 통해 도움을 주실 수 있어요. 이 책의 Github 저장소 [이슈 페이지](https://github.com/CausalInferenceLab/Causal-Inference-with-Python/issues)를 방문해주세요. 마지막으로 이 자료가 여러분의 마음에 드셨다면 도움이 될 수 있는 다른 사람들과 공유해주시고, [한국어 번역 자료](https://github.com/CausalInferenceLab/Causal-Inference-with-Python/stargazers)와 [해당 번역본의 원서](https://github.com/matheusfacure/python-causality-handbook/stargazers)에 star 부탁드립니다!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3.9.7 ('Workspace_JS')",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "a812027cb36d52afc6bf34b1f842ed1960d479f5dea7272b4a52eb53764f32d8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}