section_7-1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas_datareader\\compat\\__init__.py:7: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  from pandas.util.testing import assert_frame_equal\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "pd.core.common.is_list_like = pd.api.types.is_list_like\n",
    "\n",
    "from pandas_datareader import data, wb # This will import the data reader\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### [Equilibrium Technology Diffusion, Trade, and Growth](https://christophertonetti.com/files/papers/PerlaTonettiWaugh_DiffusionTradeAndGrowth.pdf) by Perla, Tonetti, and Waugh (AER 2020)\n",
    "---\n",
    "### 7.1. Calibration\n",
    "\n",
    "Table of Contents\n",
    "\n",
    "- [Firm Moments from SLBD](#firm)\n",
    "\n",
    "\n",
    "- [Firm moments from BEJK](#bejk)\n",
    "\n",
    "\n",
    "- [Trade moments](#trade)\n",
    "\n",
    "\n",
    "- [Productivity growth](#productivity)\n",
    "\n",
    "\n",
    "- [Real interest rate](#real_rate)\n",
    "\n",
    "\n",
    "- [Output of computations for calibration code](#output) and used in Section 7.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a name=\"firm\"></a> Firm/Establishment Dynamic Moments\n",
    "\n",
    "The micro-level data that we use is from the Synthetic Longitudinal Business Database (SLBD) (U.S.Census Bureau (2011)). The SLBD is a public access database that synthesizes information on establishments’ employment, payroll, industry, and birth and death year. It is synthesized from the confidential, restricted-access Longitudinal Business Database that is maintained by the U.S. Census Bureau. \n",
    "\n",
    "We construct SLBD moments using only establishments with at least 20 employees, which is the cutoff value used in Hurst and Pugsley (2011). Following Hurst and Pugsley (2011), the motivation for this sample restriction is that many of these small establishments appear to have no intention to grow or innovate and, thus, their motives do not correspond with the motives of firms in our model.\n",
    "\n",
    "#### How to access the SLBD\n",
    "\n",
    "Instructions on how to access the SLBD and access our code is at [``/src/calibration/SLBD_instructions/SLBD_instructions.md``](https://github.com/jlperla/PerlaTonettiWaugh.jl/blob/master/src/calibration/SLBD_instructions/SLBD_instructions.md).\n",
    "\n",
    "#### Entry\n",
    "This moment is measured as the share of employment (out of total employment) that is employed in new establishments each year. Averaging across all time periods, this value is two percent.\n",
    "\n",
    "- BGP entry (and exit rate): 0.02\n",
    "\n",
    "On the balanced growth path, the entry rate exactly corresponds with the exit rate and, hence, this moment directly pins down the exit shock faced by firms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transition Moments\n",
    "\n",
    "We construct these moments in the following way. For each year, we sort establishments into quartiles by employment size. We then compute the probabilities that establishments transition across size-quartiles at a five-year horizon, conditional on continuing to operate. We then average these transition probabilities across all years in the SLBD, i.e., from 1977–2000. The focus on the five-year horizons is to sweep out any transient fluctuations that may occur at higher frequencies due to shocks that are outside of our model. \n",
    "\n",
    "To give some perspective on the overall size of establishments, establishments larger than about 40 employees are in the top two quartiles; establishments larger than about 100 employees are in the\n",
    "top quartile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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>to_one</th>\n",
       "      <th>to_two</th>\n",
       "      <th>to_three</th>\n",
       "      <th>to_four</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>quartile</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.260</td>\n",
       "      <td>0.290</td>\n",
       "      <td>0.285</td>\n",
       "      <td>0.161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.094</td>\n",
       "      <td>0.136</td>\n",
       "      <td>0.270</td>\n",
       "      <td>0.501</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          to_one  to_two  to_three  to_four\n",
       "quartile                                   \n",
       "3          0.260   0.290     0.285    0.161\n",
       "4          0.094   0.136     0.270    0.501"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trans_moments = {\"to_one\": [0.26, 0.094], \"to_two\": [0.29, 0.136], \"to_three\": [0.285, 0.27], \"to_four\": [0.161, 0.501], \n",
    "                \"quartile\": [3, 4]}\n",
    "\n",
    "trans_moments = pd.DataFrame(trans_moments)\n",
    "trans_moments.set_index(\"quartile\", inplace = True)\n",
    "\n",
    "trans_moments.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These transition probabilities are quite similar when restricting the sample to only manufacturing establishments, with\n",
    "the only difference being a bit more persistence. Moreover, they look similar to the productivity transitions reported in Baily, Hulten, and Campbell (1992)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "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>to_one</th>\n",
       "      <th>to_two</th>\n",
       "      <th>to_three</th>\n",
       "      <th>to_four</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>quartile</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.18</td>\n",
       "      <td>0.288</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.05</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.30</td>\n",
       "      <td>0.49</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          to_one  to_two  to_three  to_four\n",
       "quartile                                   \n",
       "3           0.18   0.288      0.32     0.21\n",
       "4           0.05   0.150      0.30     0.49"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "manuf_trans_moments = {\"to_one\": [0.18, 0.05], \"to_two\": [0.288, 0.15], \"to_three\": [0.32, 0.30], \"to_four\": [0.21, 0.49], \n",
    "                \"quartile\": [3, 4]}\n",
    "\n",
    "manuf_trans_moments = pd.DataFrame(manuf_trans_moments)\n",
    "manuf_trans_moments.set_index(\"quartile\", inplace = True)\n",
    "\n",
    "manuf_trans_moments.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### <a name=\"firm\"></a> Exporters from BEJK \n",
    "\n",
    "This moment is constructed from the finding of Bernard, Eaton, Jensen, and Kortum (2003), who use the 1992 U.S. Census of Manufactures to report that 20 percent of manufacturing establishments export. Because we associate each variety in our model with the output of an establishment within the entire U.S. aggregate economy, we adjust the Bernard et al. (2003) number.\n",
    "\n",
    "In the SLBD data, when dropping establishments with < 20 employess, only 16.5 percent of establishments are in manufacturing and we assume, as a lower bound, that only manufactures are exported. This implies that the fraction of all establishments that export is 0.2*0.165 = 3.3 percent.\n",
    "\n",
    "We also target the finding of Bernard et al. (2003) that exporters’ domestic shipments are 4.8 times larger than non-exporters’ shipments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bejk = pd.DataFrame([0.033, 4.8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a name=\"firm\"></a> US Imports\n",
    "\n",
    "We construct the aggregate import share as the value of imports of goods and services divided by GDP. Over the 1977–2000 time period, the average import share was 10.63 percent and quite stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "codes = [\"GDP\", \"IMPGS\"] \n",
    "\n",
    "trade = data.DataReader(codes,\"fred\",1977)\n",
    "        \n",
    "trade = trade.resample(\"Y\", convention= 'end').sum()\n",
    "\n",
    "trade[\"imports_gdp\"] =  100* (trade.IMPGS / trade.GDP )        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over the 1977–2000 time period, the average import share was: 10.627\n"
     ]
    }
   ],
   "source": [
    "avg_trade = trade.loc[trade.index < '2001'].imports_gdp.mean()\n",
    "\n",
    "print(\"Over the 1977–2000 time period, the average import share was:\", round(avg_trade,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAGbCAYAAAARGU4hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5hU1f3H8c+BpQuKYkGQiAZBYoXV2GMBRUVRQY3EiCUWLFHUn7FFxRI19oIIYkcxFgwoFqyxmywoilKlSVFABEUQlt3z++M7+7Ds3Nk6587OzPv1PPPs7D135n6XC8xnzz33HOe9FwAAAMJpkOkCAAAAch2BCwAAIDACFwAAQGAELgAAgMAIXAAAAIEVxHmwXr16+ddeey3OQwIAANSWS9cbxdrDtXTp0jgPBwAAUC9wSREAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgsCoDl3PuEefcYufc5Ii2S51z3jnXJkx5AAAA2a86PVyPSepVcaNzbhtJPSXNS3NNAAAAOaXKwOW9f0/SsoimuyRdJsmnuygAAIBcUqsxXM65oyUt8N5Pqsa+ZznnipxzRUuWLKnN4QAAALJajQOXc665pKskXVOd/b33w733hd77ws0337ymhwMAAMh6tenh2l5SR0mTnHNzJLWXNNE5t1U6CwMAAMgVBTV9gff+S0lblH2fCF2F3vulaawLAAAgZ1RnWohRkj6W1Nk5N985d0b4sgAAAHJHlT1c3vuTqmjfNm3VAAAA5CBmmgcAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBFy5en9e0IXAAAABUNG5bWtyNwAQAAlLdmjXT33Wl9SwIXAABAeSNHSt99l9a3JHABAACUKS2Vbrst7W9L4AIAAChz993StGlpf1sCFwAAgGQD5S+5JMhbE7gAAAAefVQ655xgb0/gAgAA+W3kSOmMM4IegsAFAADy17PPSgMGSN4HPQyBCwAA5Kd//1vq39/uTKyob9+0HorABQAA8s+4cdIJJ0glJcltRx0lPf10Wg9H4AIAAPll/HjrwSouTm7r1Ut67jmpceO0HpLABQAA8sc770h9+tjyPRX16CGNHi01aZL2wxK4AABAfvjgA6l3b+nXX5PbDjhAGjNGatYsyKEJXAAAIPd98ol0xBHSqlXJbXvvLb38stS8ebDDE7gAAEBumzDBxmb9/HNy2x57SK++KrVsGbQEAhcAAMhdkyZJPXtKK1Ykt+2+u/T669LGGwcvg8AFAABy01df2UD4H39Mbtt5Z7tbsXXrWEohcAEAgNwzbZp0yCHS0qXJbV26SG++KbVpE1s5BC4AAJBbZs6UDj5Y+v775LZOnaS335a22CLWkghcAAAgd8yZY2Fr4cLkto4dLWy1bRt7WQQuAACQG7791sLWt98mt3XoYGGrffv46xKBCwAA5IKFCy1szZ6d3Lb11ha2tt029rLKELgAAEB2+/57GyA/c2Zy25ZbWtjafvv46yqHwAUAALLX0qU29cPUqcltbdpIb70lde4cf10VELgAAEB2WrbMJjWdPDm5rXVrm/rhd7+Lv64IBC4AAJB9VqyQDjtM+vzz5LaNN5beeEPaddf460qBwAUAALLLzz9Lhx8uFRUlt7Vsacv1dO8ef12VKMh0AQAAIIfMmiXdc4/08cfS2rWS97bd++RH1Pbq7LtypV1OrKhFC+mVV6Tf/z78z1lDBC4AAFB3P/8s3XyzdMcdFrTi1rSp9NJL0n77xX/saiBwAQCA2istlUaOlC6/XFq0KDM1NGkijRkjHXRQZo5fDQQuAABQO598Il14ofTf/2auhkaNpNGjpUMPzVwN1cCgeQAAUDMLFkh//rO0996ZDVtbbSWNHSsdcUTmaqgmergAAED1rF4t3Xmn9I9/SKtWpd6vSRPp0kulvn0l5zZ8SMnbUm2vbFtBgS1C3SA7+o4IXAAAoHLe22W7Sy+V5sypfN++faXbbpM6doyltGxB4AIAAKlNmiRddJH07ruV77fLLtLdd9frgeuZlB39cAAAIF5LlkgDB0rdulUetjbbTBo6VJowgbBVCXq4AADAesXF0gMPSNddJy1fnnq/ggLp/POla66xdQtRKQIXAAAwr70mDRokTZ1a+X6HHSbddZe0447x1JUDqryk6Jx7xDm32Dk3udy225xzU51zXzjnXnTObRK2TAAAEMz06VLv3rY+YWVhq1Mn6eWXpVdfJWzVUHXGcD0mqVeFbW9I2sl7v4uk6ZKuSHNdAAAgtBUr7M7DnXaSxo1LvV+rVtLtt0uTJ0tHHrl+egZUW5WBy3v/nqRlFbaN996vS3z7iaT2AWoDAAAhlJRII0ZYj9Udd9i4rSjOSWeeKc2YIV1yidS4cbx15pB0jOE6XdK/UjU6586SdJYkdejQIQ2HAwAANeK9zZ81caLdTThunPTFF5W/Zv/9pXvukXbfPZYSc12dApdz7ipJ6yQ9lWof7/1wScMlqbCw0NfleAAAoAreS998sz5cTZxoj2XLqn6tJG2zjV0+PP54Lh2mUa0Dl3NugKTekg7x3hOkAACIW2mpXe4rC1cTJkiffWZjs2qqWTPp8sttTFfz5umvNc/VKnA553pJ+pukP3jvK1lMCQAApEVJiTRt2vpeq7JwtXJl3d/7pJOkW2+13i0EUWXgcs6NknSgpDbOufmSrpXdldhE0hvOuhs/8d6fE7BOAADyx7p10pQpG4arzz+vfMHo2igstPm09tsvve+LJFUGLu/9SRGbHw5QCwAA+a24WLrlFhtD9dNP6X3vZs2k3XaTune35Xr23FPq2pVxWjFhpnkAAOqDdeuk/v2l55+v+3u1aGF3F5aFq+7dpc6dbTkeZAR/8gAAZFpJiXTKKbULW61arQ9XZQGrUyepYcP014laI3ABAJBJJSXSaadJo0ZVve8mm2zYa9Wtm7T99lKD6iwcg0wicAEAkCmlpdLZZ0tPPhnd3qOHDWwvC1cdOzLmKksRuAAAyATvpfPOkx5OcR/a0KHSOUwAkCvogwQAIG7eSxdeKD34YHT7vfcStnIMgQsAgDh5b7O533dfdPsdd0gXXBBvTQiOwAUAQFy8l668Urrzzuj2m2+WLr443poQCwIXAABxGTzYJjZN1Xb55fHWg9gQuAAAiMNNN1moinL11dI118RbD2JF4AIAILTbbrNQFeWyy6Trr4+3HsSOwAUAQEh3322hKspFF9klRubWynkELgAAQhkyRBo0KLrt3HNt8DxhKy8QuAAACOGhh6Tzz49uO/NMmxaCsJU3CFwAAKTbY4/Zkj1RTj3VJjxl/cO8wtkGACCdRo6UTj/d5tyq6E9/kkaMIGzlIc44AADp8uyz0oAB0WHrhBOs56thw9jLQuYRuAAASIcXX5T695dKS5PbjjvOer4KCuKvC/UCgQsAgLp66SXpxBOlkpLktqOOkkaNkho1ir8u1BsELgAA6uLVV6V+/aTi4uS2Xr2k556TGjeOvy7UKwQuAABq6803pWOPldauTW7r0UMaPVpq0iT+ulDvELgAAKiNd9+Vjj5aWrMmue3AA6UxY6RmzeKuCvUUgQsAgJr64AOpd29p9erktv32szFdzZvHXxfqLQIXAAA18fHH0uGHS7/8kty2117SK69IG20Uf12o1whcAABUV1GRDYRfuTK5rbDQBtC3bBl/Xaj3CFwAAFTHl19Khx0m/fRTcttuu0mvvy5tskn8dSErELgAAKjKtGl21+GyZcltO+8svfGGtOmm8deFrEHgAgCgMrNnS4ccIi1enNy24442NUSbNvHXhaxC4AIAIJX586WDD5YWLEhu++1vpbfekrbYIv66kHUIXAAARPn+e+vZmjMnua1DBwtbbdvGXhayE4ELAICKfvjBxmxNn57c1rat9PbbFrqAaiJwAQBQ3ooVdjfi5MnJbW3a2Jit7bePvy5kNQIXAABlVq6UjjhCmjAhuW2TTexuxK5d468LWY/ABQCAZMv09OkjffRRcttGG0mvvWbzbQG1QOACAGDtWqlfPxubVVGzZtK4cdLvfx9/XcgZBC4AQH5bt07q39/WQKyocWNpzBjpgAPirws5hcAFAMhfpaXSaadJL7yQ3FZQID33nNSzZ/x1IecQuAAA+cl76ZxzpJEjk9saNLDtRx8df13ISQQuAED+8V4aNEh66KHo9kcekU48Md6akNMIXACA/HP11dI990S3DRkiDRgQbz3IeQQuAEB+uekm6R//iG67/Xbp3HPjrQd5gcAFAMgfd99tvVtRBg+WLrkk3nqQNwhcAID8MGyYjduKctll0t//Hm89yCsELgBA7nvySWngwOi288+XbrlFci7empBXCFwAgNz2/PPSqafanYkVnX66DZ4nbCE0731sD0m+No9u3br5VLp161ar97QfPdqZZ55Z6/csKiqKfM9hw4bV+j2HDRsW+Z5FRUW1fs8zzzwz5c/PeeI8cZ44Tzl5niTvyz9OOsn7des4T/XtPNWjf08+jRmoQACQTy64wMbxdOokde0q7bijPaJ6P+qbX36RZs+WVq7MdCXZ75hjpMcflxo2zHQlyBMELgD5Ze1ae0yaZI/66ocfbBD3rFn2mD1b+v77TFeVG3r1kp55RmrUKNOVII8QuADUf95b2JgyxR7PP5/pisKbO1e68cZMV5F7DjzQ1k1s0iTTlSDPOB9jN3phYaEvKiqK7XgAstDcudKXX0pTp64PWFOmSMuXZ7qyaG3brr80WfZ1q62kb79d3ztV1kM1a5b044/pr6FZM+mf/7QJOxvk8b1Qa9ZYmDr//Og/5732ksaPl1q2jL82ZKu03U1B4AJQP8yfbwsJjxuX6Uqy1yGH2BqAHTpkupJ4TZ8uDR8uPfaYXYqNsvvu0ttvS5tsEmtpyHppC1xcUgSQedOnSz17SvPmpfd927ZdPyi+S5f1zxs0WN9z9vXX679+9116j59uDRpI7dtb71nUL8tvvSXtvLNNczBgQG5PdbBmjfTiixa03nmn8n27drWeLcIWMogeLgCZNXGiDWJesqR2r2/QQNpuu+RQ1aVLzT9gf/wxOojNnVu72mqjdWv7ecoeHTuuf96hgw30/uADC1SzZqV+n6OPtjCy5Zbx1R6HmTPt53r0UWnp0qr3/+1vpffes/AN1ByXFAHkgP/8RzrqKOnnn6vet0kTqXPn9YGq7NGpk9S0adg6f/ll/Ziy8kHsm2+kkpKavVejRtK22yaHqbLvqxsSV6605WiGDk29T5s20oMPSn371qzG+mbtWunf/7aled5+u/qv23dfadQoaZttwtWGXEfgApDlxo6VTjjBLg1V1KaN9dCUD1a/+U39mzNpzRppxozkILZ8udSu3YZhqixQtWuX3p9j/HibLX3BgtT79O8v3X+/9Z5lk2++kR56yHqzFi+u3muaNJH69ZPOOkvaf//cvqyKOBC4AGSxJ56wkBDVO/Sb30hvvGE9V6ieH3+ULrzQ1gtMZeutpYcftsu39VlxsTRmjPVmvflm9V/XubOFrAEDpM02C1cf8k3aAlce3z8MICPuuss+FKPC1u9+J334IWGrplq3thD7wgvS5ptH77NwoXT44XYnaH2cqX7WLOnKK+3y3/HHVy9sNW4snXSS9O671rt48cWELdRbBK588dFH0p/+JJ12mvTxx5muBvnIe+nqq+1DMcrvf2+Dm9u1i7euXHLccdLkydKxx6beZ9gwaZdd7M8604qLpdGjpcMOk7bfXrr55urNpt+pk3TbbXYZ9emnpT/8gUuHqPe4pJgPnn3WfgssLV2/7eSTbaJE7txBHEpKbDLKBx+Mbu/Z0z54N9oo3rpylffSyJG2buSKFdH7OGdrSt50U/ibDiqaNcsubz7ySPWn4mjUyALl2WfbbPEELMSDMVyoprFj7Q6ldeuS2zbaSLrmGhv70bhx/LUhP6xdK51yivSvf0W3H3+8jT1iqZX0mz/fxsq98UbqfXbc0S5HFham//hr19qlvkmTpC++WP+1ugPgJZvWoWxs1hZbpL9GoHIELlTD+PF2y/3atZXvt8MONlFifR9Mi+zzyy/WKzF+fHT72WdLQ4bUv7sPc4n31rN46aXSqlXR+zRsKF11lV3yre2Czt99t2GomjTJwlbUL3tVKSiwy6Jnny0ddFB+L1eETCNwoQrvvWcBavXq6r/m6KOlO++0sRRAXS1bJh15pPTJJ9HtV15pizNzaSgeM2dKp55qNyWk0q2b9Xb97nep90lHr1Uq221nvVmnnpp7E7YiWxG4UIlPP5V69Ii+E8m56CVByjRpYr8JX3GF1KJFuBqR2xYutIHQkydHt99xR+rB8winpMR+qbr66tQ9302aWBAeNMhm/09Xr1UqBQXSMcdYb9bBB9ObhfqGwIUUPv/cuuCXL09uK5t1+dZbbXbq8oPoK9pmG+n22218DT0QqImZM20Q/Jw5yW0NG0ojRlgPBjJn8mQbV/fZZ6n32WijcNNHNGpkaz7262d3Tm+1VZjjAHVH4EKEr7+226Oj1hcrLLR5bTbe2L6fNMnuYHr//crf88ADpXvvtf8cgapMmmQ9W1G39jdpYgPn+/SJvy4kW7vW7lC86aaaL09UE1ttJe26q01FUfa1S5fajxUD4kXgQgUzZ0oHHCAtWpTctvPONjHgpptuuN176ZlnpP/7v8qXBWnYUDr3XGnw4OxbGgTxef99u0kjahqCli3tjtkDD4y9LFThf/+z3q6pU+v2Po0bS127bhisdtmFOwuR7QhcKGfuXAtb8+Ylt3XubAsEVzYAdeVK6R//sHE1ld3R2KaN7Xf66dxVVluLFtml3Fyb3HPcOLs89OuvyW1t2kivvSZ17x5/Xaie1attXNddd1U+xrNM27YbBqtdd7X/a+i1Qu6JL3A55x6R1FvSYu/9Toltm0r6l6RtJc2RdIL3/seqDkbgCmDhQgtb33yT3Naxo/U6VPfDfeZMGyj78suV79e9uy2Eu9deNa83XxUXW1AdOdK+P/xw6YEHpG23zWhZafHUU6mX6unQwaaE6Nw5/rpQc//5j/SXv9j/BdL6XquKlwRTLR8E5J5YA9cBklZKeqJc4PqnpGXe+1ucc5dLau29/1tVByNwpdmSJTZma8qU5Lb27W1qiI4da/6+r7xik6GW/aebyoAB0i23MOC1OgYPlq67bsNtLVrYDQwDB2bvnVn33Sf99a/RbV26WNjaZpt4a0LdlJTYzTdNm9ocffRaIb/Fe0nRObetpJfLBa5pkg703i9yzrWV9K73vspfYQlcafTjj3Y34qRJyW1bbmlha4cdav/+a9bY5YUbb7TJK1Np2VK69lobgM9s9dG++krafXfr5Yqy//62zEk2LdjsvYXIwYOj2/fYw4J7mzbx1gUA6ZW2wFXbX6u39N4vkqTEV0ZFxumnn2xS06iwtdlmdjdiXcKWZHeUXX65NG2a1L9/6v1+/tnm7dp118qXD8lXJSV2iSZV2JLssu8uu9g0HCHvFkuX0lLr1UoVtg45RHrrLcIWAJQT/DqGc+4s51yRc65oyZIloQ+X+1atknr3lv773+S2Vq2k11+Xdtopfcdr187G6Lz3noWqVKZOlQ491JZxiZp/KV8NGZJ6pvXyfv3V7hbdZx/rEauviott4fP7749uP+44G0DfsmW8dQFAPVfbwPV94lKiEl9TrungvR/uvS/03hduzkDLuvn1V5uROWrurBYtpFdfDXcn2P77SxMmWICobGqIF1+0GqJ63/LNnDk2Y39Fm2yS+jX//a9dfrzhhsp7xTJh3jybQ2vUqOj2M86Qnn2WRagBIEJtA9dYSQMSzwdIGpOecpBScbF0wgnRl+2aNpVeesl6R0Iqm49rxgzpnHNSz0C/bJnUt2/0fEz5wntbqiRqseAxY+ySW6obGoqLpWuusXFQEyeGrbMq3tuda/36Wb2vvhq932WXSQ89xHQhAJBClYHLOTdK0seSOjvn5jvnzpB0i6SezrkZknomvkcoJSV2Geell5LbGjWSRo+2AfRx2WwzWxpowgRbLijKN99Yj0eM87zVK08+aXfoVTRwoE3jcfDB0pdf2t2gqYLrpEnSnnvaIs9R81uFtHq1DeTfbTebrPSFF1IvBXXrrfZgCSgASImJT+u70lKbv+nxx5PbGjaUnntOOvbY+Osq47309NM2/ihqlvt77kk9bUCuWrxY2nFH6+krr317G5/VqtWG2z/6yM7xtGmp37NLF+mRR6S9905/veXNm2fzgz30UHL9FTVoIA0fbsEaAHJTxu9SRBy8l847LzpsOWe9KJkMW2V1/OlPdtkpaqD0pZdKn34af12Z9Ne/RoeVoUOTw5Zkl4I//9zuCk11SW7qVOtNvPji6MuUdVF22bBvX7tseOutVYetli0t7BO2AKBaCFz1lfcWVh58MLp9xAjppJPirakynTrZJaiKysaeVfUBnivGjrUFmiv64x/t7tJUmjaVbr7Zwukuu0Tv473NjbbzztI779S91lWr7O9R2WXD0aNTXzYs07atDeifOdPuSAQAVAuBq7669lrpzjuj2+67zy5B1TfHHy+df37y9nnzbFb6qj7Ms92KFTZGq6JNN7VLq9XRvbstJjx4cOoZvmfNsjFgAwfanGw1NW+e9aZts4105pnSF19U/Zp99rGFzufOtTX3WJAYAGqEwFUf3XKL9SJEufXW6FBTX9x+u1RYmLz95ZetLZf97W+2tmVF99xTs4DSuLHdpThhQvSfZZkHH7Q51157rer39F56992aXTZs3NiCclGR9OGH0oknsswLANQSg+brm3vvtTvXolxzTerZveuT2bOlbt2k5cs33N6woX3o77dfRsoK6j//sctyFR1+uE0EWts7+Nats57Oa66x5ZZSGTDA9tt00w23r1plNzXce6/dFVkdW29tvWdnnUVPFoB8F+9aiulC4KrCiBF2iSfKpZdK//xn9tx6P3asTZJZ0dZb2wDxXJoEd/Vqm4V/xowNt7doYXcl/uY3dT/GtGk2QP3DD1Pvs9VWdofhscfapb8HHrC/U9UdP7fPPjbg/7jj6MkCAMNdijnnqaesRyHKuedmV9iSpKOPtpBY0cKFNqdYNqwZWF2DByeHLckGwacjbElS5862vNK990rNm0fv8913Fpa6dZO2287+znDZEADqBXq46oNRo6Q//zk6hJx6qt391yALs3FxsV1m++ij5Lbrr5f+/vfYS0q7iRNtctKK527vvW0JphAzr8+ebT2hb71V+/fYemsL8meeyWVDAEiNHq6c8cADNo9VVNg68US7JJSNYUuynpJ//ctmpq/o2mvrFhjqg3XrpL/8JfncNW5s5y3UMjcdO9oSTw89FD2vV2X23dfOyZw50lVXEbYAICZZ+kmeA7yXbrzRJjaN6mXs08cmNs32tenat5dGjky+HOq91L9/9Oz02eKOO6TPPkvefvXVUteuYY/tnIW9r76Sjjyy8n0bN7ae0gkTpA8+sHnRuGwIALEicGVCaal0ySWpL6kdeqj1QuTKh2KvXtabUtHixTZ567p18ddUV9OnS9ddl7x9p51seoi4tG9va2yOHJl8h2K7dhbqv/1WevRRG9sFAMgIxnDFbd06Gzfz2GPR7X362ASTTZvGWlZwJSVSz57RM6RfeaV0003x11RbpaW2WPh77224vUED6eOPbUxXJixebOstLlwo7b+/dMwxuRPaASAzmBYiK/36q/Xo/Pvf0e0DBtjYn4KCeOuKy3ff2TIy33+f3PbKKzZnVTYYNkw655zk7YMGpV4dAACQjRg0n3V+/lk64ojUYeuii6x3IlfDlmTzRI0aFX0TwMkn26Wv+m7BAumyy5K3d+yYenUAAEDeI3DFYelSW/su1YLDN9xgPSPZejdiTRx0kE0JUdGyZTaYe+3a+GuqLu9tKoWo9QuHD7eJTgEAiJAHn/AZNn++dMABNrFkRc5JQ4bYXW3ZNKlpXV1xhXTYYcnbP/nE2uqr556zGfQrOu00qUeP+OsBAGQNxnCFNH26DRSfNy+5raBAevxxmxohHy1dauO5FixIbhs92panqU9++EHacUdpyZINt2+5pTRlitS6dWbqAgCExBiueu+zz2yR5qiw1bSpjeXK17AlSW3a2NQXUfOMnXaaNGtW/DVV5uKLk8OWZD2UhC0AQBUIXCG8/74taRP1Ad2qlTR+fNWTVeaDffeVbrklefuKFTaea82a+GuK8vrr0hNPJG8/9lipb9/46wEAZB0CV7qNG2cTl0YNrN5iC+ndd22OJJhLLpGOOip5+4QJ1pZpK1dKZ5+dvH3jjaX774+/HgBAViJwpdPTT9tkk7/+mtzWoYP1fO2+e/x11WfO2Vi2bbdNbhsyxC47ZtLVV0tz5yZvv/12WwAaAIBqIHCly5AhNpdU1DI1XbpIH34o7bBD/HVlg9atpWefjZ4V/S9/sZsPMuHjj6V7703eftBB0hlnxF8PACBrEbjqqmwR6vPPj16EurDQerbat4+/tmyyxx7Rs7SvXCn16yetXh1vPWvWWNireE6bNbM5t/JpGg8AQJ0RuOqitNTuXku1CPXBB0tvv2135KFq550nHX988vYvv5QuuCDeWm6+Wfr66+Tt118v/fa38dYCAMh6zMNVW+vWWQ/I449Htx9zjC1jk2uLUIf2009S9+7SzJnJbY89ZutNhvbVVzbWrrh4w+3du9vkrEieQ04AAA5sSURBVLm8/BIAoDzm4cqoX3+1y1ypwtapp9qs5IStmmvVSnr+ealJk+S2gQOlyZPDHr+kxMZnVQxbBQXSww8TtgAAtULgqqmffrJFqMeMiW4fNIgP5rraddfoKRdWr7ZLjitXhjv2ffdJn36avP2yy6wuAABqgUuKNbF0qXT44dHrIko2eP7KKxlQnQ7e2+XDJ59Mbuvd2y7ZNmsmNW+e+mvZ8+qG39mzpZ12klat2nB7587S55/TYwkA+SdtH+gErur69lub0HTq1OS2skWoBw6Mv65c9ssv0p57Rg9er4lGjaoOZ82a2eD8L79Mfv3779syTQCAfJO2wMV1r6oUF0tjx9rdiKkWoX7iCemkk+KvLde1aGFj4fbYI7nXqSaKi+0RNft/Vc49l7AFAKgzxnClsmCBdN11NgN6v36pF6EeM4awFVLXrtKwYZk5dvv2Nj0EAAB1RA9XeaWlNm/W0KEWpEpKUu+78cbSyy/T+xGHk0+2cXP33BPvcR980O6aBACgjghckrRsmU3xMHSoNGNG1ftvsYX0+uvSbruFrw3mrrukHj1sPNXKlXbH4qpVyV+jttVmnOJll0lHHpn+nwMAkJfyO3D973/SAw9IzzwTveB0lO7dbULTTp3C1oYNOWd3J/buXbPXeS+tXRsdxKK+FhdL3bpJe+8d5ucAAOSl/Atcq1ZZYBo6VJowoXqvadjQpiEYONCW62Hah+zhnE2i2qSJLZINAEAG5E/gmjrVxuQ8/ri0fHn1XrP11tJZZ9kSPu3aha0PAADkrNwOXMXFNvj9gQekd96p/ut69LDerKOOsjmcAAAA6iA3A9f8+dLw4dKIEdKiRdV7TevWtgbiOedIO+wQtDwAAJBfcidwlZZKb75pY7NeeqnyKR3K22MP68068USbcRwAACDNsj9w/fCD9NhjNj5r5szqvaZZM6l/fwta3bsHLQ8AACC7A9eaNXb5b9my6u3fubOFrFNO4Y41AAAQm+xe2qdJE5uuoTIFBbY0z1tvSVOmSBdeSNgCAACxyu7AJdniwlHatZMGD5bmzrUFkJk/CwAAZEh2X1KUbAzWHnvYrPGS1LOnhbDeva13CwAAIMNyI5Fceqn06ac2pQNL7gAAgHrG+dos7FtLhYWFvqioKLbjAQAA1EHaxiJl/xguAACAeo7ABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgdUpcDnnBjnnvnLOTXbOjXLONU1XYQAAALmi1oHLOddO0l8lFXrvd5LUUNIf01UYAABArqjrJcUCSc2ccwWSmktaWPeSAAAAckutA5f3foGk2yXNk7RI0grv/fiK+znnznLOFTnnipYsWVL7SgEAALJUXS4ptpbUR1JHSVtLauGcO7nift774d77Qu994eabb177SgEAALJUXS4p9pA023u/xHtfLGm0pH3SUxYAAEDuqEvgmidpL+dcc+eck3SIpCnpKQsAACB31GUM16eSnpc0UdKXifcanqa6AAAAckZBXV7svb9W0rVpqgUAACAnMdM8AABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACB1SlwOec2cc4975yb6pyb4pzbO12FAQAA5IqCOr7+Hkmvee/7OecaS2qehpoAAABySq0Dl3OulaQDJJ0qSd77tZLWpqcsAACA3FGXS4rbSVoi6VHn3GfOuRHOuRYVd3LOneWcK3LOFS1ZsqQOhwMAAMhOdQlcBZK6SRrqvd9d0i+SLq+4k/d+uPe+0HtfuPnmm9fhcAAAANmpLoFrvqT53vtPE98/LwtgAAAAKKfWgct7/52kb51znRObDpH0dVqqAgAAyCF1vUvxAklPJe5QnCXptLqXBAAAkFvqFLi8959LKkxTLQAAADmJmeYBAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgRG4AAAAAiNwAQAABEbgAgAACIzABQAAEBiBCwAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAiMwAUAABAYgQsAACAwAhcAAEBgBC4AAIDACFwAAACBEbgAAAACI3ABAAAERuACAAAIjMAFAAAQGIELAAAgMAIXAABAYAQuAACAwAhcAAAAgTnvfXwHc26JpLmxHTA7tZG0NNNF5DnOQeZxDjKPc5B5nIPMa+q93ykdb1SQjjepLu/95nEeLxs554q894WZriOfcQ4yj3OQeZyDzOMcZJ5zrihd78UlRQAAgMAIXAAAAIERuOqf4ZkuAJyDeoBzkHmcg8zjHGRe2s5BrIPmAQAA8hE9XAAAAIERuAAAAAIjcAXmnHvEObfYOTe53LZdnXMfO+e+dM695JxrldjeyDn3eGL7FOfcFeVe08s5N805N9M5d3kmfpZsla5zkGhv6Jz7zDn3ctw/RzZL47+DQc65r5xzk51zo5xzTTPx82SjGp6Dxs65RxPbJznnDkxsb+6cG+ecm5o4D7dk6MfJSuk4B+XahjvnpifORd8M/DhZyTm3jXPuncT/LV855y5MbN/UOfeGc25G4mvrxHbnnLs38dn7hXOuW7n3GpDYf4ZzbkCVB/fe8wj4kHSApG6SJpfb9j9Jf0g8P13SDYnn/SU9k3jeXNIcSdtKaijpG0nbSWosaZKkrpn+2bLlkY5zUO51F0t6WtLLmf65sumRpn8H7STNltQs0faspFMz/bNly6OG5+A8SY8mnm8haYLsF/Tmkg5KbG8s6X1Jh2f6Z8uWRzrOQeL7wZJuTDxvIKlNpn+2bHlIaiupW+J5S0nTJXWV9E9Jlye2Xy7p1sTzIyS9KslJ2kvSp4ntm0qalfjaOvG8dWXHpocrMO/9e5KWVdjcWdJ7iedvSCr77cRLauGcK5DUTNJaST9J2lPSTO/9LO/9WknPSOoTuvZckaZzIOdce0lHShoRuuZck65zIJusuVmirbmkhSHrziU1PAddJb2VeN1iScslFXrvV3nv30lsXytpoqT2gUvPGek4B4m20yXdnGgr9d4zG301ee8Xee8nJp7/LGmK7Je5PpIeT+z2uKRjEs/7SHrCm08kbeKcayvpMElveO+Xee9/lJ27XpUdm8CVGZMlHZ14frykbRLPn5f0i6RFkuZJut17v0z2l+Hbcq+fn9iG2qvpOZCkuyVdJqk0xjpzWY3Ogfd+gaTbE9sWSVrhvR8fb8k5J9U5mCSpj3OuwDnXUVL3cm2SJOfcJpKOUiIUoNZqdA4Sf+6SdINzbqJz7jnn3JbxlpwbnHPbStpd0qeStvTeL5IslMl6FaXUn781/lwmcGXG6ZLOc85NkHVprk1s31NSiaStJXWUdIlzbjtZV2ZFzOdRNzU6B8653pIWe+8nZKTa3FTTc9Ba9ttmx0RbC+fcyfGXnVNSnYNHZB8gRbJfND6StK7sRYkexlGS7vXez4q14txT03NQIOtV/NB7303Sx7JfRFADzrmNJL0g6SLv/U+V7RqxzVeyPaVY11KE8d5PlXSoJDnndpBdppJs7Mpr3vtiSYudcx/KupC/1Ya/XbYXl1LqpBbnYHdJRzvnjpDUVFIr59xI7z0f+LVUi3PgJc323i9JvGa0pH0kjYy79lyR6hx479dJGlS2n3PuI0kzyr10uKQZ3vu746s2N9XiHPwgaZWkFxNNz0k6I8aSs55zrpEsbD3lvR+d2Py9c66t935R4pLh4sT2+Yr+/J0v6cAK29+t7Lj0cGWAc26LxNcGkq6W9GCiaZ6kgxN3RbSQDdCbKhtU2ck519E511jSHyWNjb/y3FHTc+C9v8J73957v63sz/9twlbd1OLfwTxJeyXulHOSDpGNv0AtpToHiT/jFonnPSWt895/nfj+RkkbS7ooI0XnmJqeA28jtl/S+g/7QyR9HXfd2Srxf8fDkqZ47+8s1zRWUtmdhgMkjSm3/ZTE/0d7yYYyLJL0uqRDnXOtE73vhya2pUQPV2DOuVGyfxhtnHPzJV0raSPn3HmJXUZLejTxfEji+WRZd+Wj3vsvEu9zvuxkNpT0iPf+q9h+iCyXrnOA2kvjv4PnZQO110n6TCx9Um01PAdbSHrdOVcqaYGkPyfeo72kq2QBeKJ9dul+7z03klRDOs5Bwt8kPemcu1vSEkmnxVB+rthX9mf5pXPu88S2KyXdIulZ59wZsl/ujk+0vSK7U3GmrGfxNEny3i9zzt0g6xCRpOvLjfeNxNI+AAAAgXFJEQAAIDACFwAAQGAELgAAgMAIXAAAAIERuAAAAAIjcAEAAARG4AIAAAjs/wHl8LStgZph8AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (10,7))\n",
    "\n",
    "ax.plot(trade.index, trade[\"imports_gdp\"], color = \"red\",\n",
    "        linewidth = 5, label = \"Imports\")\n",
    "\n",
    "ax.axhline(y= avg_trade,          # Set the value equall to the average\n",
    "           color='k',       # make the color red\n",
    "           linestyle='--',   # the line style\n",
    "           linewidth=4)     # thickness of the line\n",
    "\n",
    "ax.spines[\"top\"].set_visible(False)\n",
    "ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "ax.set_xlim(dt.datetime(1977,1,1), dt.datetime(2000,1,1))\n",
    "ax.set_ylim(5,15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a name=\"productivity\"></a> TFP Growth\n",
    "\n",
    "To measure U.S. productivity growth, we use the [Bureau of Labor and Statistics Multifactor Productivity Database](https://www.bls.gov/mfp/home.htm). We focus on the private non-farm business sector and measure labor-augmenting productivity growth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_path = \"src/calibration/data\"\n",
    "\n",
    "input_path = input_path + \"/LATFP.xlsx\"\n",
    "## This is the historical multifactor dataset. \n",
    "\n",
    "bls = pd.read_excel(input_path, sheet_name = \"XG\",skiprows=6,nrows=72, headers = True, usecols = \"A:G\")\n",
    "# Use the non-farm, private sector, sheet XG\n",
    "\n",
    "# Note the orginal source of the data is\n",
    "# https://www.bls.gov/mfp/special_requests/mfptablehis.xlsx\n",
    "# As of 4/06, it appears that there was a revison of the historical tables.\n",
    "# The TFP growth rate is 0.80 vs. 0.79 used in the paper. \n",
    "\n",
    "bls.dropna(inplace = True)\n",
    "\n",
    "bls.set_index(\"Year\", inplace = True)\n",
    "\n",
    "bls.columns = [\"RVA\", \"L\", \"K\", \"NVA\", \"NLC\", \"NKC\"]\n",
    "\n",
    "bls.drop(2018, axis = 0, inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "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>RVA</th>\n",
       "      <th>L</th>\n",
       "      <th>K</th>\n",
       "      <th>NVA</th>\n",
       "      <th>NLC</th>\n",
       "      <th>NKC</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year</th>\n",
       "      <th></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>2013.0</th>\n",
       "      <td>12138.731787</td>\n",
       "      <td>7024.670066</td>\n",
       "      <td>4376.594278</td>\n",
       "      <td>12312.7</td>\n",
       "      <td>7072.07</td>\n",
       "      <td>4481.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014.0</th>\n",
       "      <td>12520.291128</td>\n",
       "      <td>7201.584199</td>\n",
       "      <td>4500.352539</td>\n",
       "      <td>12922.7</td>\n",
       "      <td>7429.04</td>\n",
       "      <td>4697.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2015.0</th>\n",
       "      <td>12964.111194</td>\n",
       "      <td>7375.836859</td>\n",
       "      <td>4637.590224</td>\n",
       "      <td>13476.3</td>\n",
       "      <td>7825.63</td>\n",
       "      <td>4835.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016.0</th>\n",
       "      <td>13173.415497</td>\n",
       "      <td>7505.505874</td>\n",
       "      <td>4777.426761</td>\n",
       "      <td>13839.6</td>\n",
       "      <td>8047.16</td>\n",
       "      <td>4962.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017.0</th>\n",
       "      <td>13528.844631</td>\n",
       "      <td>7649.337123</td>\n",
       "      <td>4916.218202</td>\n",
       "      <td>14451.1</td>\n",
       "      <td>8435.82</td>\n",
       "      <td>5162.35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 RVA            L            K      NVA      NLC      NKC\n",
       "Year                                                                     \n",
       "2013.0  12138.731787  7024.670066  4376.594278  12312.7  7072.07  4481.46\n",
       "2014.0  12520.291128  7201.584199  4500.352539  12922.7  7429.04  4697.01\n",
       "2015.0  12964.111194  7375.836859  4637.590224  13476.3  7825.63  4835.53\n",
       "2016.0  13173.415497  7505.505874  4777.426761  13839.6  8047.16  4962.26\n",
       "2017.0  13528.844631  7649.337123  4916.218202  14451.1  8435.82  5162.35"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bls.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "bls[\"labor_share\"] = bls[\"NLC\"] / bls[\"NVA\"]\n",
    "\n",
    "reg_labor_share = 0.66\n",
    "\n",
    "# Labor's share is nominal compensation vs. nomional value added.\n",
    "\n",
    "bls[\"TFP\"] = (bls[\"RVA\"]/(bls[\"L\"]**(bls[\"labor_share\"])*bls[\"K\"]**(1-bls[\"labor_share\"]))).astype(\"float\")\n",
    "\n",
    "# TFP measure\n",
    "\n",
    "bls[\"LA_TFP\"] = (bls[\"TFP\"]**(1/(bls[\"labor_share\"]))).astype(\"float\")\n",
    "\n",
    "# Labor Augmenting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_year = 1976\n",
    "\n",
    "mean_tfp = 100.*np.log(bls.loc[start_year:2000].LA_TFP).diff().mean()\n",
    "\n",
    "#100.*np.log(bls.loc[start_year:2000].LA_TFP).diff()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Labor Augmenting TFP growth\n",
      "0.79\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean Labor Augmenting TFP growth\")\n",
    "print(round(mean_tfp,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Standard Deviation Labor Augmenting TFP growth\n",
      "3.3\n"
     ]
    }
   ],
   "source": [
    "print(\"Standard Deviation Labor Augmenting TFP growth\")\n",
    "\n",
    "std_tfp = 100.*np.log(bls.loc[start_year:2000].LA_TFP).diff().std()\n",
    "print(round(std_tfp,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standard deviation of these decadal averages\n",
      "0.78\n"
     ]
    }
   ],
   "source": [
    "mean_50 = 100.*np.log(bls.loc[1949:1959].LA_TFP).diff().mean()\n",
    "\n",
    "mean_60 = 100.*np.log(bls.loc[1959:1969].LA_TFP).diff().mean()\n",
    "\n",
    "mean_70 = 100.*np.log(bls.loc[1969:1979].LA_TFP).diff().mean()\n",
    "\n",
    "mean_80 = 100.*np.log(bls.loc[1979:1989].LA_TFP).diff().mean()\n",
    "\n",
    "mean_90 = 100.*np.log(bls.loc[1989:1999].LA_TFP).diff().mean()\n",
    "\n",
    "mean_00 = 100.*np.log(bls.loc[1999:2009].LA_TFP).diff().mean()\n",
    "\n",
    "mean_10 = 100.*np.log(bls.loc[2009:2017].LA_TFP).diff().mean()\n",
    "\n",
    "decadal_avg = [mean_50, mean_60, mean_70, mean_80, mean_90, mean_00, mean_10]\n",
    "\n",
    "df_dec_avg = pd.DataFrame(decadal_avg, columns=['DateValue'])\n",
    "\n",
    "std_dec_avg = df_dec_avg.std()[0]\n",
    "\n",
    "print(\"standard deviation of these decadal averages\")\n",
    "print(round(std_dec_avg,2))\n",
    "\n",
    "outdf = pd.DataFrame([std_dec_avg])\n",
    "outdf.to_csv(\"./src/calibration/data/dec_std_avg.csv\", index = False, header = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAGbCAYAAAAV7J4cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5jUxf0H8Pdc4eDggKN5ICCgAiIiICrGho3YsfxQo0YUFXsh9hJjYmKMLagxUSIQu9gi9oYRrOjRexEQONrR23F1fn/MruzuzN5t+Zb57r5fz8Mjzu3tDiy7996Zz/czQkoJIiIiInJfjt8TICIiIsoWDF5EREREHmHwIiIiIvIIgxcRERGRRxi8iIiIiDyS58eDnnzyyfLjjz/246GJiIiIkiWcuiNfVrw2bNjgx8MSERER+YpbjUREREQeYfAiIiIi8giDFxEREZFHGLyIiIiIPMLgRUREROQRBi8iIiIijzB4EREREXmEwYuIiIjIIwxeRERERB5h8CIiIiLyCIMXERERkUcYvIiIiIg8wuBFRERE5BEGLwqWbduAYcOA7t2BSy8Fysv9nhEREVHC8vyeAFFSRowAxo9Xv1+8GNi4EXjvPX/nRERElCCueFFwSAlMmBA99v77wNy5/syHiIgoSQxeFBzbtgG7d+vjY8Z4PxciIqIUMHhRcGzaZB5/4QWgstLbuRAREaWAwYuCI17w2rhR34IkIiKykCPBSwgxUggxVwgxRwjxqhCisRP3SxQlXvACuN1IRESBkHbwEkLsDeBGAAOklL0B5AK4IN37JdLUF7w++wxYvtyzqRAREaXCqa3GPABNhBB5AAoBrHbofon2qC94SQmMG+fdXIiIiFKQdvCSUpYBeBTACgBrAGyVUn4aezshxAghRKkQorScTS8pFfUFLwAYOxaorfVmLkRERClwYquxGMAQAF0BdADQVAhxceztpJSjpZQDpJQD2rZtm+7DUjZqKHitWqW2HImIiCzlxFbjiQCWSSnLpZTVAN4G8CsH7pcoWkPBCwCee879eRAREaXIieC1AsBAIUShEEIAOAHAfAfulyhaIsFrwgRg/Xr350JERJQCJ2q8pgB4E8A0ALND9zk63fsl0iQSvGpqgBdfdH8uREREKXDkqkYp5R+klD2llL2llL+VUrKNODnPFLxOOUUfe+45dZUjERGRZdi5noLDFLxuvVUfW7AA+PZb9+dDRESUJAYvCgYpzcHr8MOBQYP0cXayJyIiCzF4UTDs2gVUVUWPNWoEFBYCV1yh3378eGDbNm/mRkRElCAGLwoG02pXq1aAEMA55wAtWkR/bdcu4LXXvJkbERFRghi8KBg2b9bHWrVS/23SBLhY69nLnl5ERGQdBi8KhngrXmGm7cYffwRmzXJvTkREREli8KJgaCh49e0L9O+v34ZF9kREZBEGLwoGU/AqLo7+f9Oq14svArt3uzMnIiKiJDF4UTA0tOIFAL/5jar3irR5M/DOO+7Ni4iIKAkMXhQMiQSvli2B//s//XYssiciIksweFEwJBK8APN248SJwNKlzs+JiIgoSQxeFAyJBq+jjwb2318fHzfO+TkRERElicGLgiHR4CUEcPnl+vi4cUBNjfPzIiIiSgKDFwVDosELAIYNA3Jzo8fKyoBPPnF+XkRERElg8KJgSCZ4lZQAZ5yhj7PInoiIfMbgRcGQTPACzEX2778PrF3r3JyIiIiSxOBF9qusBHbujB7LyQGaN4//Pb/+NdChQ/RYTQ3wwgvOz4+IiChBDF5kP9MB2cXFKnzFk5cHXHaZPj5mDCClc3MjIiJKAoMX2S/Zbcaw4cP1sUWLgK+/Tn9OREREKWDwIvulGry6dQOOP14fZ5E9ERH5hMGL7Jdq8ALMRfZvvAFs3ZrenIiIiFLA4EX2Syd4nX22qgeLVFEBvPpq+vMiIiJKEoMX2S+d4NW4MXDxxfo4txuJiMgHDF5kv3SCF2Debpw6FZg+PfU5ERERpYDBi+yXbvDq0wc49FB9fMyY1OdERESUAgYvsl+6wQswr3q9/LKq9yIiIvIIgxfZz9RANdngdcEFQGFh9NiWLcDbb6c+LyIioiQxeJH9nFjxat4cOO88fZxF9kRE5CEGL7KfKXjFtohIhGm78csvgSVLkr8vIiKiFDB4kf2cWPECgF/9CujRQx8fOzb5+yIiIkoBgxfZrbZW1WLFSmXFSwjzqtd//gPU1CR/f0REREli8CK7mUJX8+ZAXl5q93fJJfr3rlkDfPRRavdHRESUBAYvsptT24xh7doBZ56pj7PInoiIPMDgRXZzOngB5u3GDz4AVq9O736JiIgawOBFdnMjeA0eDHTsGD1WWws8/3x690tERNQABi+ymxvBKzcXuOwyfXzMGEDK9O6biIioHgxeZDc3ghcADB+urnKM9NNPwKRJ6d83ERFRHAxeZDe3gleXLsCJJ+rjLLInIiIXORK8hBAthRBvCiEWCCHmCyGOcOJ+iVwLXoC5yP6tt8xnQxIRETnAqRWvJwB8LKXsCeBgAPMdul/Kdm4GryFDgNato8d27wZeecWZ+yciIoqRdvASQjQHcAyAMQAgpaySUhq6XhKlwM3gVVAA/Pa3+ji3G4mIyCVOrHh1A1AOYJwQYroQ4jkhRNPYGwkhRgghSoUQpeXl5Q48LGUFN4MXAFx+uT42YwYwbZpzj0FERBTiRPDKA9AfwL+klP0A7ARwZ+yNpJSjpZQDpJQD2rZt68DDUlZwO3j17g0cfrg+zlUvIiJygRPBaxWAVVLKKaH/fxMqiBGlz+3gBZiL7F9+Gdi1y9nHISKirJd28JJSrgWwUgjRIzR0AoB56d4vEerqzMGruNjZxzn/fKBpzO74tm3Am286+zhERJT1nLqq8QYALwshZgHoC+BBh+6Xstn27Sp8RSosBBo3dvZxiopU+Io1Zoyzj0NERFnPkeAlpZwRqt/qI6U8S0rJRkiUPi+2GcNM242TJwOLFrnzeERElJXYuZ7sZWpk6lbwGjgQ6NVLH+eqFxEROYjBi+zl5YqXEOZVr+efB6qr3XlMIiLKOgxeZC8vgxegmqnm50ePrVsHfPede49JRERZhcGL7OXFFY2R2rQBTj1VH58+3b3HJCKirMLgRfbyesULMDdTnTnT3cckIqKsweBF9vIjeB18sD7G4EVERA5h8CJ72RK85s4FamrcfVwiIsoKDF5kLz+CV4cOQOvW0WOVlcDChe4+LhERZQUGL7KXH8FLCG43EhGRaxi8yF5+BC+AwYuIiFzD4EX2YvAiIqIMw+BFdpKSwYuIiDIOgxfZqaJCFbVHys8HmjZ1/7EPOADIy4seW7sWWL/e/ccmIqKMxuBFdoq32iWE+49dUKDCVyyuehERUZoYvMhOfm0zhnG7kYiIXMDgRXZi8CIiogzE4EV2YvAiIqIMxOBFdrIxeM2frxf8ExERJYHBi+zkd/Bq1w4oKYkeq6lR4YuIiChFDF5kJ7+DF8DtRiIichyDF9mJwYuIiDIQgxfZicGLiIgyEIMX2WnzZn3MluAlpbfzICKijMHgRXayYcWrRw/VxT7Sxo3A6tXezoOIiDIGgxfZyRS8iou9nUNeHnDggfo4txuJiChFDF5kJxtWvADWeRERkaMYvMg+VVXAjh3RY0IALVp4P5e+ffUxBi8iIkoRgxfZx1RYX1wM5Pjwz5UrXkRE5CAGL7KPLduMANCnjz62aBFQUeH9XIiIKPAYvMg+NgWv4mKgc+fosbo6YM4cf+ZDRESBxuBF9rEpeAHcbiQiIscweJF9GLyIiChDMXiRfRi8iIgoQzF4kX2CELxmzeLRQURElDQGL7KPbcFr332Bpk2jx7ZuBX7+2Z/5EBFRYDF4kX1sC145OcBBB+nj3G4kIqIkMXiRfWwLXgDrvIiIyBGOBS8hRK4QYroQ4n2n7pOyFIMXERFlKCdXvG4CMN/B+6NsxeBFREQZypHgJYToCOA0AM8l/E3l5U48NGUiG4OXqcbrp5+A7du9nwsREQWWUyteowDcDqAu3g2EECOEEKVCiNK6GTOAoUMdemjKKLW1wJYt+nhxsfdziVRUpK5ujDV7tvdzISKiwEo7eAkhTgewXko5tb7bSSlHSykHSCkH5NTWApMnA2Vl6T48ZRpT6GreHMjL834usbjdSEREaXJixetIAGcKIZYDeA3A8UKIlxr8LimBN95w4OEpo9i4zRjG4EVERGlKO3hJKe+SUnaUUnYBcAGAL6SUFyf0za++mu7DU6bZvFkfY/AiIqIM4W8frx9+AJYu9XUKZJmgrXjNng3UxS1tJCIiiuJo8JJSfimlPD2pbxo/3skpUNCZgpffhfVh++wDtGgRPbZzp7q6kYiIKAH+d67ndiNFsnnFSwigTx99nNuNRESUIP+D1+zZwNy5fs+CbGFz8AJY50VERGnxP3gB3G6kPRi8iIgog9kRvF57TbWXIGLwIiKiDGZH8Fq8GJg2ze9ZkA1sD169ewM5MS+bFSvMbTCIiIhi+BO8mjXTx157zft5kH1sD15NmgDdu+vjs2Z5PxciIgocf4KX6Qfp+PHsh0T2By+A241ERJQyf4JXcTGQmxs9tnIl8O23vkyHLMLgRUREidi8GTjrLNXqp2tX4Isv/J5RQvwJXnl5wEkn6ePcbsxuUtrdQDWMwYuIyH9//CMwYYL6/fLlwEUXqabWlvOvuP6CC/SxN94Aamq8nwvZYft2oLY2eqxJE/XLJqbgNWcO/+0SEXlFSuCVV6LH1q4NxKqXf8HrrLOAgoLosfXrgS+/9GU6ZIEgbDMCQIcOQOvW0WOVlcCiRf7Mh4go28ydC5SX6+Off+79XJLkX/Bq0QI49VR9nEcIZa+gBC8huN1IROSneCtbn33m7TxS4G8fL9N249tvq9UDyj5BCV4AgxcRkZ/+9z/z+Pz5QFmZt3NJkr/B6/TTgaZNo8e2bAE+/dSf+ZC/GLyIiKghtbXApEnxvz5xondzSYG/wauwEBgyRB/ndmN2YvAiL82YAUyZol/QQUR2mzmz/tNCLN9u9P/IINN244QJgbgkNJDefBPo2xfo1w/473/9nk20IAWvAw5QbVEirVljLvYku0ipLjvv1w8YOBA45hhg1y6/Z0VEiYq3zRj2+edWn//sf/AaPBho2TJ6bNcu4IMP/JlPJps3Dxg6VH1amDEDOOccYOFCv2e1R5CCV0GBCl+xuOplv/fei74M/dtvgTFj/JsPuWP7dmDUKGD4cO6iZJqGWkasXauuerSU/8GroAA491x9nM1Unffyy/rY+PHezyOeIAUvgNuNQWV6b7F8a4KSNHcuMGAAMHIkMG4ccOGFwMMP+z0rckJ1NTB5csO3s7ithP/BCzBvN374IbB1q/dzyWSLF+tjK1Z4P494GLzIbVVV5tX0qVO9nwu54/XXgcMP1/vq3XuvXSv8lJqpU4EdOxq+HYNXAwYNAtq1ix6rrATeeceX6WSsJUv0sVWrvJ9HPKZiSQYvctIXXwDbtunjq1erXxRcNTXAbbcB559vrhGurgZuuMHq2h9KgKm+a5999LEvv1QftCxkR/DKywPOO08f53ajc6QEfvpJH7cpeGXCitf8+da+2An1X1DCVa/gWr9e1Qs/+mj9t/vsM+Ctt7yZE7nDFLxuvFE/03fnTnXVsoXsCF6Aebvxs894lZhTNm40f9K3PXjZdkB2pHbtgJKS6LHqahW+yD61tXsO1DVh8AqmKVOAQw5p+Eq3sJEjE9uqIvtUVgJff62Pn3gicPzx+riltZv2BK8jjgA6dYoeq63lpxOnmLYZAVVHt327t3OJJ2grXgC3G4NkyhRg3br4Xy8t9W4u5Ix//1u1A4n3AfLww/WxVauAP//Z3XmRO374AaioiB5r0wbo3VuFr1iW1nnZE7xycsyrXtxudIZpmzHMhuMVKiqA3bujx/LygGbN/JlPohi8gqOhvnWlpaz/CYrdu4ErrgBGjDBv7efkAH/7G/Ddd+YzgR97jCvTQWRqIzFokHq+TzpJ/9oPP1h5kZ49wQswB6/Jk+0IBkEXb8ULsOPvN95qlxDezyUZDF7BIGXDwWvdOhbYB8GKFcDRR8fvvdamjTp27vbb1fvHk0+qtkWRamqA669n0A4a03byccep/3brBnTpEv212lpVZG8Zu4JXv37A/vtHj0mpLg+m9NS34mVDnVcQtxmB+MGLb+h2mTOn/tdAGLcb7TZxoqrnivc8DRigavVOOGHP2L77Anfeqd/2iy/4syVIKirUCmascG2XEOZVLwu3G+0KXkJwu9EtDF7u6NFD/zS9YYM6PojskejxWCywt5OUautw8GD1+jK54grgq6+Azp31r91xB9C1qz7+u9/ZU+NK9fv2W31buaREvQeHmeq8LCywtyt4Aebg9cMPwNKl3s8lk9S31cjglbq8PODAA/VxbjfaxRS8evfWx7jiZZ/t29VRZ3feCdTV6V9v1AgYPVoV2jdubL6PJk3UlmOs1auBP/3J2fmSO0zbjMcfH12OEvv/gGqau3Klu3NLkn3Bq1cvoE8ffZyrXqnbvl31uYmHwSs9rPOy27Jl6mzSWKYfuCywt8uCBcBhh8W/ur1jR7XKdeWVDd/X6acDZ5yhj48aZfW5fhRiKqwP13eFtWmjSpZiWbbdaF/wArjd6LSGalsYvNLD4GU3U++ugw9WP4hjt4nLy+14Pdhg0yZg7FjgxRfV8TteB9K331aha8EC89ePO05tDR92WOL3+cQT+qpYTQ1w3XUM3Dbbvh348Ud93NS7KwBtJYITvGbP5qeSVDF4uYvBy26mbcazzwby883PHbcb1ZWDPXsCl18OXHKJqqNp316dMPKPfwCzZpm3/ZxQWwvcdRdw7rnx669uu01duRh71FxDunYF7r5bH580CXj11eTnSt74+msVkCN17myu24tXYO/Wv9cU2Bm8unY1N77jqldq6qvvAlSxamwPLa9lWvBauFBv9EfeKy83d7o++2z13wED9K8xeKli9NhTQ9atA954Q513ePDBQOvWauvukUdUc9rq6vQfd8MG4OSTgYceMn+9aVN1JeLDD6v6ylTcdpu60jHWLbeYT/cg/8VrI2FqN3TkkfpK9vr16spmS9gZvID4241cDk5eIpfR+92/KMjBq7hYP3Whro4rtDZ49139k263bsBBB6nfH3KI/j3ZfmXj4sWJtVnYsgV4/33VL2vgQKBlS7XN86c/qd5JyX7wKC1Vz0e8baHu3dWFVkOHJne/sRo3Bp56Sh9fuxa4//707pvcYarvMm0zAupCiqOP1sct2m60N3idd56eZpcsAaZN82c+QZZI8PJ7uzHIwQvgdqOt4m0zht9b4q14ZfMHvIcfTm1bZtcu1WfrD39QqxEtWqjVh7vuAj76qP4O4mPHAkcdpbY4TYYMUaGrV6/k52VyyinAWWfp408+qcpayB5btgDTp+vjsYX1kSxvK2Fv8OrQATj2WH2c243Ja2irEWDwSheDl322bze/2Ya3GQH1gzy22HrjxvgBINOVlQHPP6+Pp3KCRHW16r300EPq2J5WrdSK1s03q8L58nJ16PFVV6lasspK8+P+5S/q9i1aJD+H+owapVZHItXWstDeNpMn6x8E9ttP32WIZApekyeb/435wN7gBcTfbrSoSM56lZWJ9TBh8EoPg5d9PvpIb7i4117AEUfs+f+8PKBvX/17s7XO67HH9Fqtli1VfdekScADD6ji5aZNk7/vujq1Y/HEE6pwvl079QF79Gjz7Vu1Aj7+WBXD57jwo2qffYB779XHv/oKeOkl5x+PUpNIG4lY/frpPz927TJ3vveB3cHr3HP1AspVq9SnKErMsmWJfXpj8EoPjw6yzzvv6GNDhug/xFlgr2zcCDz7rD5+ww1A27bAMceooPLpp8DmzWrr79FHgTPPTP21anrdA+oH59SpqlO9m265RT+mDgBuvVVtcZH/4jVOrU9OTvSxUWGW1HnZHbzatDFfGsrtxsQlUt8F+Bu8qqqAHTuix4RwfmvBTfvuCxQWRo9t3Zq9W1Z+q6oCPvhAHzfV9bDAXnnySbUqEKmwELjxRv22+fnAoYeq4DJhgto2nD0bePpp4PzzVeuJVF16KfDNN/qBx24oKDAX2q9fr2rVyF/l5ap1SaxBgxr+XlN2sKTOK+3gJYToJIT4nxBivhBirhDiJicm9gvTduMbb+g9PcjMFLy6ddPHysrcn0s8mzfrY8XF7mwvuCU3d8+VcpG43eiPL77QWwMUFZk/KbPAXtXDmY7UGTFCfQBuSE6OOoLp2mvVB+OyMlVbOnasClKm95xY+fnAv/6lvie29spNv/612l2J9Y9/8PXrt0mT9LEDDlBnNDbEVOdVWmr+eeMxJ36y1QC4RUp5AICBAK4TQjh06QnUJ1RTTw7T8iPpTIX1posW/FzxCvo2YxjrvOxhuprxtNP09xJANQqNXa3cvFlt02eLZ5/Vt9by89WKViqEUKvAl10GjBunPgCuWqWalF5zjX6+aceO6ofs1VenVsifrr//Xf83UFengiRriv2TyjZjWNeuer+2ujrV6sRnaQcvKeUaKeW00O+3A5gPYO907/cXzZurN8xY3G5MjGnFyxS81qxxpgFiKkyfQBi8KFW1teZjgiKvZowUr8A+W7Ybd+9WRfWxLrlEBSKn7L232sH45z9VM8sNG9R28Lvvqt5hkRc9eK1TJ+C++/Txb78FXnjB+/mQkkphfSRL20o4upcjhOgCoB+AKYavjRBClAohSstjOyI3xLTd+Pbb1lwaajXTiteBB6qu05GkVA0E/cAVL3LS99+rq/AiFRSo3k3xZHOB/X/+o7/2c3JU93o3tW6t2kyccYbe0sMPI0eqo5Fi3X67FdtTWWfNGvM5nYnUd4VZem6jY8FLCNEMwFsAbpZSaucuSClHSykHSCkHtG3bNrk7P+00oFmz6LEtW4BPPkl9wtmgtta8XbLvvuZPsn5tN5qCV3Gx9/NIV58++thPP+kXDpC7TNuMJ56oarziMQWvbFjxqqlRDVNjDR1qvtovkzVqpOq6YpWXA7//vffzyXambcbwUVWJOv54fet68WLg55/Tm1uaHAleQoh8qND1spTybSfuM0phoboMPBa3G+u3apW+fdiqlQo1tgevIK54FRXpRcRSshO2l6Q0t5GIt80YZrqyMRsK7MePN384u/NO7+digxNPVKemxPrXv3hqitfSqe8KCzftjeXzqpcTVzUKAGMAzJdSPp7+lOIwbTdOmADs3OnaQwaeaZsxXGzI4OUObjf6a84cva4xJ0f1mqpPjx56U9CtWxNvxxJEdXXAX/+qj596qrnmLVs8/rj+b4GF9t5Lt74rzNRWIujBC8CRAH4L4HghxIzQr1MduN9ogwerDsqRdu1Sh7SSmemHxn77qf+agpdfLSUYvMgppm3GI49UDUDrk5urmnbGyuTtxvffNx/kfvfd3s/FJnvvbT4se8oUdYUmuW/FCmDp0uixnBzVxDdZ8eq8fAzRTlzV+LWUUkgp+0gp+4Z+fejE5KI0amTutcLtxvi44uU9Bi9/xTsUOxHZVGAvJfDgg/r40UeroJrtbrrJfCD3HXfE77ZPzjFtMx5ySGpNtX/1K70v3IYN5sasHglQh0oAv/mNPvbhh/Wfep/Nkl3xYvBKnyl4zZrFLQovLFsGzJihjzN46b78Uq3gxMr21a6w/Hxzof3GjcA993g/n2zj1DYjoK6YPfpofdzHthLBCl6DBqlDbiNVVZmLaan+Fa+9Da3WGLzS16WL6j0XaedOfdmcnGd6H+jbN/GjZ0xFuNOmZWZoNq129eunuriTctxx5g/7zz6buYHcBlI6U1gfybK2EsEKXrm56jLnWK++6v1cbCelecWrvq3GsjJ/fshkUvASwtxWgtuN7kvlasZI3bvrbWu2bTN/gAmyH380/9C56y5/usbb7NFH9TYkUqpC+9paf+aU6X76CVi5MnosLy+9LXBTgf3kyap5sA+CFbwA8yeQzz9XvVZoj/Xr9Ss+Cwv3nHFVVKSvzNTUqO/zWiYFL4B1Xn4oLwe+/lofTyZ45eRkx4HZpisZu3cHzjnH+7nYrkMH4I9/1Md//BEYM8b7+WQD02rXYYfpH4qS0aePfoHN7t3qZAIfBC94DRwIdO4cPVZbC7z1lj/zsVW8bcbIT7Q21HnV1upnxAHBbKAalg3Ba+ZM4IknzGHHD+++q6/W7ruvOrg5GfH6eWWKefPMFyDccYfaUSDd9deb/x3ddZcq0iZnmeq70tlmBNSHqhNO0Md92m4MXvDKyQHOP18f53ZjtPoK68NsaCmxdavepLKoSBW3BlWmB6/x41Xt1M03q6JVU+dzr8W7mjHZrbNML7B/6CF9rGNH4OKLvZ9LUOTnA08/rY9v2sSLEZwWr74r1cL6SBad2xi84AWYtxu/+sq/4nAb1VdYH2bDilembTMC6tNxTsxL6+efzSt7QbN9O3DNNdFjd90FTJ/uz3wANSfTG+hZZyV/X5lcYL98OfDKK/r4bbepdj0U3zHHmMPpc8+Zrw6l1CxYYD5n1YkD1E11XlOn+tIeJJjBq29fVZMQSUrgjTf8mY+N6iusD2PwckdhofmcOx/7xjjmX//SDwz2u6v3Rx+pq5sj7bVXam/W++2n1z7u2AEsWpT6/GzxyCN6QXibNsAVV/gzn6B55BH934aUwHXXsdDeKaZtxiOO0PtwpaJzZ/19WUrzY7osmMFLCPMRQtxu3CPVrUYGL2dk4nZjRQXw2GPmr33/PTB2rLfzCTNtMw4Zoq86JiJegX3QtxvXrjUXg998s/qgQA0rKQEeeEAfnzoVGD3a+/lkIqfbSMSypK1EMIMXYA5eP/6Y2WerJSORrUYbenkxeAXHc8/Vf9XrHXd4X2xcWQl88IE+nszVjLEy8crGUaPU31WkoiK1WkOJu/Zac7uYu+/mlfXpqqtzr74rzLTd6EOdV3CD1wEHmH+4jR/v/Vxss2WL6rAcKT8f6NQpeowrXu7JtOBVVdVwEf2mTarey0v/+5+q8YrUvHl6n5IzrcB+yxbgn//Ux6+9Vj//luqXl2cutN+yRTWgvftuYOFC7+eVCWbP1n8eFBaqVvFkzOcAACAASURBVBJOOe44fSV86VLPG1wHN3gB5lUvnt1oXvXr0kW9aUSKd1Vj7FWGbsqm4DVnjuqVFkQvvJBYKH/uObXt6BXTNuNpp6VXLG5a8Zo+Pbh1PE8/rYfTxo2BkSP9mU/QHXUUMGyYPl5Wpnqk9eypzgf89795nF0yTLVWRx3l7IUfLVsChx6qj0+c6NxjJCDYwcvUVmL2bGDuXO/nYpNECusB1SsrtmixokIvnnZTpgavvffW/xy7dwOLF/szn3TU1JjbEBx3nF5sDKiVFC8CZm0tMGGCPp7ONiOgXiuxh/Hu3BnMlYydO9U2Y6zLL9ePX6PEPfxw/Qc2f/cdMGKEqgu76CK1nRXU4O4Vt7cZwyxoKxHs4NW1q2qoGivbV70SKawH1EUKfm83ZmrwEiJzthvHjzf/m3rkEeDPf9bHp09XVz+67fvvzZeen3xyevcrROZsNz73nF53l5sL3HqrP/PJFO3aAU8+2fDtdu9WLTwGD1a7DvfcE8wPX26rqQEmTdLHnSysDzPVeU2c6OlV2cEOXkD87UYvt8tsk0hhfRiDl3syIXjV1ZkPVD75ZLUld801qr1LrHvvVVfSucm0zXjiifrZeqnIhAL7qip11mCsiy5K/OBwiu+SS9QVcWeemVjX/1Wr1Gupe3e1hfbcc+osUFIf1mL/LoqKgP79nX+sgQP1K3k3bfK0F2Hwg9fQoXp36iVLgvcm6aREV7wA/4OXaVuTwcse77yjjpmJde+96r95eebVrW3bVGNOt0gZv1u9EzJhxeull8yv5Tvu8H4umeqEE9R2d1mZarWS6BFV33wDXHml2oq8+GIV4DKhSW+qTNuMxx6r1yU7oaBANcSN5WFbieAHrw4dgEGD9PFs3m5MZsXL75YSXPGyl5TmrcRjjwWOPHLP/w8caG7C+dJLwJdfujO3OXP0K5FyctTqgxNMwWv69OBcHFFba67LO/tsoFcv7+eT6fbaC/jd71ST5KlTgRtuSOx9rKICePlltf3VpQvw+9+b378znamw3o36rjDTdiODV5JM243jx2fnJ4iKCv28RSFUPZyJ3ytepuAV5AOyI/XqpX9iW706OAfrfvyxefk9vNoV6a9/Nf+gue46vau8E0yrXUcdBbRt68z9d+mi/zusqFBHmgTB22+ba4m8bveRbYRQ22NPPqle62++CZx+emJbkStXqg86+++vzkAdM0a/GjUTVVcDX3+tj7sZvEwF9l99pV7jHsiM4HXOOfoPuFWrgG+/9Wc+fjL1I+nYUV0+buLnQdlSZvaKV0GBurQ8VhBWveKtdh12mNpeidWmDfC3v+nj8+aZr6pLl5vbjECwC+ylNNflnXii+VJ6ckdBAXDuucB776mfR48+Chx4YGLf+/XXahW5pETVkn3ySea2pvjxR3X1baTiYvOOgVMOOkhdIBGpslJtAXsgM4JXmzbmpcPXX/d+Ln5LtJVEmJ8rXjt26Fs3jRs7cy6XLYK63ThpkvmDy7336jWVYcOHm68y/uMf1ad5pyxbBsyYoY+ncih2fYJaYP/xx+a/n7vv9n4upJSUALfcotodlZYC11+f2AfMXbuAF19UF7O0bAl066YWGv70J+Ddd4EVK4J/IZlpm3HQoNSO/EqUEP62lZBSev4LgEzlV//+/WU8/bt0Sek+1V+B2ZVXXpnyfZaWlhrv89lnn035Pp999lnjfZaWlqZ8n1defLGU6qW751fz5lKqJ8r556l//+x6nh5++Je/19IU7w+AvPLKK+P++V15noqKUn+ecnL0f1Pnnuvc8/T447/c77Np/J268nry+nnKttdTNj5PcX5dmeL98XkKPU/jxml/p/2bNIn/PDiYgTJjxQsAmjXzewbB1LixOk4o0rZtvMzZKW4ul7vl++/Tqy0xnf/31lvOrXqZthmJiJJhKplgjVeSnDxWIJsIYb6y0as6r0wXxOD1l7+k9/0PPGDuiu5E/cT69eZCXCJK3OjRqg3Gzz8j8FuVqerUCejRw5eHzpzgRanzu6VEJttrr2AdzTJjBvD+++ndR4sWqqdRLCeu0Hrvvez9QUHZ6Te/URfpxKutTMXo0aomsksXVWt24onZuZJsqg33gpP7lon+OuSQQ+Lu26Zl+HB9L/xvf3PnsWxUXS1lXp7+d7B1a/3fd/75+veMHev+fJ95Rn/c4cPdf1yvDR6s/zlfeMHvWZkNHarPtUsXKauqkrufujopBw3S76ugQMolS1Kf32mn6fd5662p319D7rlHf7xrrnHv8dJx3nn6XLt1U+8LFHw7dkj5/ffqffPqq6UcOFDKwsK4NWAp/fL6femoo/Q5PPOMd4//zjvm97u6OtOtWeNlZOoaPGeO9/Pwy4oV+lWCbduaDzKO5FdLiUxuJREpKFc2zp+v+g7FuuMOvQ6wIUIATz+tt3mprARuvDG1Vavt281XHTnZRiJWUK5sXLQIeOMNffz2293p/k3ea9oUOPxw4Kqr1GkR332nanEXLFANw++8U139mM4K+1VXqSsvvbBzJzBlij7uxvmM8QwapPdYW77c3JbJQZkVvEw9UubO9X4efkmmY30kv1pKMHjZ5aGH9EDUoQNw6aWp3V+vXqqbd6wPP1RHESXro4/0ZqwlJeYWFk4x9fKaOVM1fbTJww/rz11JCTBsmD/zIW/k5qo6pfPPV02MP/pInZG6Zo36/V//qr7Wo0diW5UVFar3mBcXV33zjf462nvv+MfbuaFFC9WbMJbLbSUyK3iZVrzmzVPHZ2SDZM5ojMTg5a54wcumWqWlS9XRJbFuvTV+891E/P735n9fN92kN01siKkGZcgQd/v9dOyod8OvrLTrA92qVcALL+jjt9yS3nNHwVVSola/7rxTrYYtWKBWjL/7Tq2WXXWVWj0z9UxcvFj15HP7/cl0PuPxxztby5YIH44Pyqzg1b69fszH7t2q4WI24IqXnXr00K+6LS9Xn0xt8fDD+geUNm2AESPSu99mzYAnntDHw8ejJKqyEvjgA33czW1GIBgd7B97TF85KC5WP1yJwpo2VavDV18NPPOMahtTXm7eKXrrLfPr1klen88Yj6mR6hdfuLpgk1nBSwjzP6JsqfNKtmt9GIOXu/Lzzf8ubdluLCsDxo3Tx0eOVG/W6Tr7bPXpO9ajj6q6skR88YV+VWTz5t68UdscvMrL1RVqsW64ASgq8n4+FCxNm6qQZeqDedtt7h2hs3Wr+TXkZX1X2OGH6+9zmzcD06a59pCZFbwA83ajTdsCbkp1q7GkRN+u2bjR/WZy2RK8ALvrvB59VK+datHC3Ag1FUIATz2lzq6LVFOjHiORLQ1TTdhpp3nTv8/mAvsnn1THykRq2lRdwECUiB491IHcsWpqVH3Y+vXOP+ZXXwF1ddFjXbsC++zj/GM1pFEjVWQfy8XtxswLXtm64iVl6iteeXkqfMVy+8pGBi/v5xFr/Xrg2Wf18RtuUOHLKfvtp+pNYv3vf8Crr9b/vbW1qtljLLe3GcPiFdhXVnrz+PFs2wb84x/6+FVXAa1bez8fCq7zzjOH9bIy4MILnd92M9V3+bHNGObxuY2ZF7yytaXEmjX6ClVRkV4YHI8fLSUYvLyfR6xRo/R/N02bquJ3p91xhzrkN9Ytt6ith3i+/x5Yty56rKAAOOUUZ+cXT4cO+geT6mr/31eeeQbYsiV6LD/ffCUpUUMeeQQ44gh9fOJE4P77nX2seIX1fjEV2H/zjb6a7JDMC16mFa+FC+27/Ntp8QrrE71CxOs6r4oKdeFDpLy8zD1z0xS8Fi7U/w68tHmzecXk6qtVYb3TmjQxP97atcAf/hD/+0xXM550knf/VoSwb7uxogJ4/HF9/NJLzSdREDWkUSNg/Hjza//Pf1btKZywaZM6ISOWnytevXqpi/MiVVWpLVEXZF7watsWaNcueqy6Wl0im8lS3WYM8zp4xVvt8vpSYq+0aqX/HdfW+lt/+NRTesF6QYFagXLLKaeYtwifesr8ZiylOXh5tc0YZluB/X/+o68C5uSohqlEqerUCXjlFfP78MUXq7Md0zVpkl7X2aOHWln2ixDm7UaX6rwyL3gB2dlI1bTilUwjOluCVyazabtx+3bz5eKXX65/8nPaqFFAYWH0WF0dcO21esHt7Nl6F+mcHOCMM9ydYyybVrxqa9UFEbHOO8/b5pOUmU46yby1uGkTMHRo+rWNtrSRiOVhnVdmBq9srPMK2orX5s36GIOXd555Rg+/eXnerJh07gzcd58+/t13aiUnkmm16+ijE69ddIopeM2e7c9W8YQJ5iNN7rrL+7lQZrr3XuDXv9bHf/wx/RpC2+q7wkzBa+ZMV67qdCR4CSFOFkIsFEIsEUIYLl3yWDaueKXaSiLMVBfCFS9n2RK8KipU081Yv/2td5dzjxwJHHCAPn777aqVSZipjYTX24yA2gaJ3QqprvbuXLtIptquk08G+vTxfi6UmXJygJdeUluPsf75T7UdmYp168w/i03tHLzWoYOq9YplWqFLU9rBSwiRC+BpAKcA6AXgN0IIw+w9lI0rXql2rQ+zYasx9tSBTGPL0UFjxpjrg0ztHtzSqJE6RDvWxo3A3Xer3y9bZq77Oussd+cWjw3bjVOmmJtaulmXR9mpTRt18Hp+vv61K69Ux/El68sv9bHevb1fwY7Ho+1GJ1a8DgOwREq5VEpZBeA1AEMcuN/UmVa8lizx9woyN23apF9WXlBgDlPxmAob161z72rQbFzx2m8//Wy0LVvU8TleqapSxwPFOu88oHt37+YBqLqOCy/Ux//9bxUwTKtd/fr502QRsKPA/u9/18cOOgg44QRv50HZ4fDDzSusu3apw7RjL85piK3bjGGmthKffeb4h2MngtfeACJ/cqwKjfmnZUt966y2Vl2+n4lMq11duyZ3eHDjxvqnDilVfzA3ZGPwys1VPyRjPf+8XlTulhdfNAe98CqT1x59VB39E0lK4JprgDff1G/vxzZjmGnFy8vg9fPP5r+T3/0uc68GJv9ddx1wwQX6+IIF6izXZEKJrYX1Ycceq2pdI61c6XhXBCeCl+kVrz0TQogRQohSIURpeXm5Aw/bgGw6OijdwvowL7cbszF4AebtxvvuUy/4RM8tTFVNDfDQQ/r4WWeZA6EX2rcHHnhAH58+Hfj2W33ctuA1d677R2uFPfWU3kG8pAT4zW+8eXzKTkKoVeiePfWvvfaauWTAZNUqPcAIod77bFFUpA4Sj+VwWwkngtcqAJEVeB0BrI69kZRytJRygJRyQFsv9nOz6eigdAvrwxi83Gd6UQPA11+rUHbffe5tib/+unl19J573Hm8RF17LdC3b8O3228/8+vaKyUl+mukpgaYNcv9x962Tf3wi3X99foZmEROa9ZMHaYd2wYGUCuuU6Y0fB+mbcZ+/eyr7fWgn5cTwetHAPsLIboKIRoBuADAuw7cb3qyacUr3cL6MAYv9w0dav7kCKh6ugceUAHM9CaVjro64MEH9fFf/9pcu+SlvDx1pVRDzj7b/y01vwrsx4xR4StSkybqXEYiL/TqZQ7/1dXqfW3Dhvq/3/b6rjBTnZfDVzamHbyklDUArgfwCYD5AF6XUvqfcLjilfz9MHi5r6hIrW4NGxb/NosWqTekyy6Lbq2QjgkTzB88/F7tCjviCNW8tT5+bjOG+VFgX1NjbnY7bJg7RzsRxXPhhWqFOtbKlaqzfX2Hadt2MHY8hx6q3qcj1XeWbAoc6eMlpfxQStldSrmvlPIvTtxn2kz9OJYuBXbu9H4ubnNqxcvLXl7ZGrwAoHVr1Sh04sT6A/J//qNWx158Mb2raqRUZ63FOuYY1YzUFg89FP/fQPv26gorv/kRvP77X/NRLTff7O7jEpk8/rgKJ7E++QT4S5wf/8uWAcuXR4/l5tr1/hOWn+96X7HM7FwPqD3pLl30cbcLmL22c6c6ZDhSTo75z94Q04pXWVlK02pQNgevsOOPVw04773X3CsHUMv3l1wCDB5sDtiJ+OQTYNo0fdyW1a6wNm3Mxf8AMGRIclfpusW01Thvnrq83i2my/nPOEOdb0fktYIC1d/L9H59//3Ap5/q46bVLtPKki1M240OsuCdzEXZ0EjVtM3YubNqUJksr7Yaq6v1/i9CAC1aOP9YtmvcWNV1TZ8OHHlk/Nt9/rm68vDBB1UvrkTFW+069FDX31xScvnl5pWtc8/1fi4mbduq11ek2lr3TiD47jvg++/18XSPbSFKxz77qM72saRU25GxLWtsbyMRy1Rg76DMDl7ZcHSQU60kAPNW4+rV9e/bp8J0TmPLlmrpOVsdeCAwebI6QzFeAN29W61SHXKIudWCyeTJ5k7n99zjf6G6SU6O2mKNPKj7nHPsahDqZT8v09FO/frZdQk+ZadTTlGr9bE2blQNmcMfEKUMTmF9WM+e5p+HDsns4JWtK16pFNYDanu2ZcvosZoa5w8J5TajWU6Oukpt/nz1xhXPnDnAUUepJqOxJxbEMq12HXSQ2qqyVc+e6krBl14C3n9fbWvYFBJNdV5uXNm4dKn5kHA2TCVb3H+/eXXo+++B225Tv1+0SH2Aj5SfD/zqV65PL2VCuLrqldnBKxtWvJwqrA/zYruRwat+7dsD48cDH3wQ/3gcKdXq2AEHqGBiKr6fMsXcf+buu+2ol6pP+/bARRcBp51m31y9KrB/8kn9RIMOHeoP5UReys1VB2abVoeefFL1DjStdh1xhLknmE1cLMWw7B3NYT176m/aK1c6fmmor5zcagQYvGxy6qnqg8Itt8QPH2vXqh/EZ5yhX/lmusJo//1Vzx1KnWmrcf58Z6+Y3rJF9e6KdcMNqdVvErmlbVsVsGKP2gFUzebzz+vjNtd3hblY3pDZwatJE3MISeVUdVuZVrxS3WoEvGkpweCVuKZN1XmGpaXmH/hhH3ygVngff1xtD8+cCbz3nn67u+7K7lo6J7RurV81XFcHzJjh3GM89xywY0f0WGEhG6aSnX71K+CRR/TxHTvMF4cEIXiVlJjLlRyQ2cELyOw6r6oqYMUKfbxbt9Tv04uWEgxeyevXT20djhqlwpjJzp1qdezww4GRI/Wvd+6smhxS+tzcbqyuNjdMHT7cvuNViMJuugn4v/9r+HaNG8c/Os02gwe7crfZGbwypc7r55/1GpCSElUknypuNdorN1e9uc2bV39x/LRp5rqKO++M3y+MkuPm0UFvvqm/5oRQzz2RrYRQ2+Pdu9d/uyOPDM75oiNHxv+gm4bMD16ZfHSQ04X1AINXEHTurI4AevPN6LYL9WnfXh1BRM5wa8VLSnPD1CFD0ishIPJC8+bqfalJk/i3sbmNRKyOHdXr+oYbHL3bzA9embzi5XRhPcDgFRRCqKai8+ers9Maai9w661qiZ+c0b+/PrZggd4YOFlff20OcGyYSkFx0EHqiut4glDfFalnT3WFpoMyP3jtv79+tcXatc4dPuwnpwvrgfjBK52zAmMxeDmnRQvg6adVk9R4haCtW7Mo22mtWum1lFKmX2BvWu069FDVt40oKC65BBgxQh9v1sy8WpxlMj94NWpkPtMsE1a93FjxatFC39OurHQ2qJo61zN4peeII1Rt11//qq9sPfigK3UKWc/p7cYlS9QWciw2TKUgeuIJ/TVy/fWsM0U2BC8gc+u8nOxaHyaE+y0luOLljvx8VUA/Z47aWjzzTHX8jumTJ6XP6QL7UaP0leVOnew5p5IoGY0bqzMab7gBGDRInUl7331+z8oKho5nGah3b9XgLVLQg1ddnTpSJFa6K16A2m5ctCh6rKwM6Ns3/fsGzMGLl8k7Z999zT11yFlOrnht2gSMG6eP33gjVwgouIqKHK+PygTZu+IV9K3GsjK1BRipZUtnVo7cLLCvqzNvNTJ4UdCYCuwXLgS2bUv+vkaPBnbtih5r1gy44orU5kZE1sqO4BWviaqTBeNei9dKwolaEDeD19at+t97s2Y8BoWCp2VL89b+9OnJ3U9VFfDUU/r4FVfoh9YTUeBlR/Dad1+9YdumTcC6df7MxwluFNaHuRm8WN9FmcSJ7cbXXwdWr44ey8lR24xElHGyI3jl5gIHHKCPB7nOy41WEmEMXkSJMRXYJxO84jVMPeccoGvX1OdFRNbKjuAFZF6dF1e8iPxnWvFK5srGL780b02yYSpRxsqe4JVph2W70UoizNROwqmDshm8KJOYCuwXL1a1jIkwrXYNHKj6shFRRsqe4JVJK15SunNOY1ibNnqx+/btqV2tFYvBizJJ8+bmQ4GnTWv4excuBN5/Xx/nahdRRsue4BXvzMYgXtm4YYN+JlzjxokfmNyQnBz3mqgyeFGmSbXAftQofaxLF+Dss9OeEhHZK3uC1z776MembNvm/AHQXoi32pXj4NPpVp0XgxdlmlQK7DdsAJ5/Xh+/6Sb9bFkiyijZE7xycoBevfTxINZ5uVlYH8bgRZSYVArsn3kGqKiIHmveHBg+3Ll5EZGVsid4AfG3G4PGzcL6MAYvosT066c3Lv7pJ/MJDYA6ceIf/9DHr7xShS8iymjZFbwy5bBsNwvrwxi8iBJTVAT06KGPx1v1evVVvXlzbi4bphJliewKXpnSUsKLrUa3WkoweFEmSnS7MV7D1KFDgc6dnZ8XEVknu4KXacVr3jx1cHOQuNm1PowrXkSJS/TKxokTgdmz9fGRI52fExFZKbuC1957Ay1aRI9VVADLlvkzn1Rs3w6Ul0eP5eY6/2nZjeAlJYMXZSbTlY2mFS/TatdRRwGHHeb8nIjIStkVvIQIfiNV0zZjly5Afr6zj1NSogJdpE2bgF27Ur/PHTuAmproscaNgSZNUr9PIhv07au3c1m2DNi4cc//z5sHfPSR/r1smEqUVbIreAHBr/PyorAeUKHL1JA1nTovrnZRpmrWDOjZUx+PXPX6+9/1r++7L3Dmme7Ni4isk33BKxNXvNwIXoDz240MXpTJ6iuwX78eePFF/es336yvLBNRRsu+4JWJK15OF9aHOR28TH2NGLwoU9RXYP+vf6n+XZFatgQuvdT1aRGRXbIveJlWvBYs0GuPbOX3ihe3GonM4h0dtHs38PTT+teuukptURJRVsm+4NWuHdCmTfRYVZV5JclGXq54OX1Qtil4FRenfn9ENjEV2K9YoQ7Djr0SOS8PuP567+ZGRNbIvuAlRHCPDqqsNAefbt3ceTzWeBElrrDQvKJ+//362Pnnm19fRJTxsi94AcE9OmjZMtULK9Lee7vXjoHBiyg5pu3G2NougC0kiLJYWsFLCPGIEGKBEGKWEOK/QoiWTk3MVUFd8fKqlUQYgxdRckwF9rEGDQL693d9KkRkp3RXvD4D0FtK2QfAIgB3pT8lDwR1xcvLwnoA6NBBH1u3TtXEpYLBizKdacUrFle7iLJaWsFLSvmplDJ8OeD3AIJRtGAKXosWmbcEbOJlYT0ANGoE7LWXPr56dWr3x+BFme7gg+vvy7X//sBpp3k3HyKyjpM1XsMBGM7DUIQQI4QQpUKI0vLYK3y81qqV3pW9tlaFL5t5veIFONtSgsGLMl2TJuZShrCRI/UrH4koqzT4DiCE+FwIMcfwa0jEbe4BUAPg5Xj3I6UcLaUcIKUc0LZtW2dmn44gNlI1BS83V7wAZ1tKMHhRNoi33diqFTBsmLdzISLr5DV0AynlifV9XQgxDMDpAE6QMvaSO4sdeCDw2WfRYzYX2NfWqqsaY/mx4sXgRRTfgAHA2LH6+DXXqJYTRJTVGgxe9RFCnAzgDgDHSil3OTMljwRtxWvlSqC6OnqsVSt17IibnApeFRXqV6TcXKCoKLV5EdnKdGVjfj5w3XXez4WIrJNuscE/ABQB+EwIMUMI8YwDc/JG0A7L9rqwPsyp4BXvnEYhkr8vIpv176+K6CNdcYVeV0pEWSmtFS8ppQc/+V3Sq5c+9tNPwK5ddm4H+FFYDzgXvLjNSNkiNxd46y3g6qvVBTtnnAE89pjfsyIiS6QVvAKteXOgc2d1llqYlOrAbBubG/pRWA8weBGl4qCDgG++8XsWRGSh7L6uOUh1Xl53rQ8zXdW4Zo0q9k8GgxcRERGDl8bWOi+/thoLC4Hi4uix2lrVwT4ZDF5ERERZHryCcnSQlP5tNQLObDcyeBEREWV58ArKite6dcDOndFjTZuaj/NxA4MXERGRI7I7ePXsqbcz+PlnYPt2f+YTT7z6Lq9aMTB4EREROSK7g1dhoblOat487+dSH7/qu8IYvIiIiByR3cELCEadl5/1XQCDFxERkUMYvILQUsKvVhJhpuBVVpbcfcTrXE9ERJRFGLyCcHSQ31uNpl5eTqx4xbapICIiynAMXkFd8bJhq1HKxO+DW41EREQMXujeXZ2tFmnNGnNQ8MOWLfpc8vOBTp28m0Pz5kCzZtFjVVXAhg2JfX91NbBtmz7esmX6cyMiIgoQBq+CAhW+Ytmy3WjaZuzaVQ+LbhIivQL7LVv0sZYtvf0zEBERWYDBC7C7zsvvwvqwdIIXtxmJiIgAMHgpNtd5+V1YH8bgRURElDYGL8Duo4P8LqwPS6elBIMXERERAAYvxeYmqraseKXTUoLBi4iICACDl7LffkCjRtFjGzYA69f7M59IfnetD+NWIxERUdoYvAAgL08dmB3L71Wvigp9O08IdVWj1xi8iIiI0sbgFWZjndfSpfpYx46qBYbX0mmiyuBFREQEgMFrDxvrvGwprAeA1q31wLdzJ7B1a8Pfy+BFREQEgMFrDxtbSthSWA+k10SVwYuIiAgAg9ce8ZqoJnMeodNsWvECUm8pweBFREQEgMFrj65dgSZNose2bk28V5UbbFrxArjiRURElCYGr7CcHKBXL33czwJ724JXqr28GLyIiIgAMHhFs6nOq7oaWL5cHw/ailddHbB5sz5eXOzMnIiIiAKEwSuSTYdlr1gB1NZGNWG0QwAADzxJREFUj7VtCzRv7s98gNSC19atep1cs2Z6w1oiIqIswOAVyaYVL9sK64HUghe3GYmIiH7B4BXJtOI1b57aLvOabfVdQGpXNTJ4ERER/YLBK1KnTkBRUfTYzp3Azz97PxdbzmiM1K6dOl4p0ubN6u8oHgYvIiKiXzB4RRLCnqODTFuNfq945eYCHTro4/WtepkK6xm8iIgoSzF4xbLl6CAbtxqB5FtKmFa8eEUjERFlKQavWDaseNXV2bnVCCRfYM+tRiIiol8weMWyYcVrzRpg9+7osaIioE0bb+dhwuBFRESUMgavWKYVr/nz9Z5aborXSkII7+YQD4MXERFRyhi8Yu21lx4MKivNW39usbW+C0i+pQSDFxER0S8YvGLFu7LRy+3GoAUvrngRERElhMHLxO+jg2zsWh/G4EVERJQyR4KXEOJWIYQUQlhQ/e0ArnjF1769Xmu2fr3ajjVh8CIiIvpF2sFLCNEJwEkAVqQ/HUv4ueIlpd0rXvn5qg4u1urV+piUDF5EREQRnFjx+juA2wFIB+7LDqbgtXAhUFXl/mNv2gRs3Ro9VlBgblzql0S3G3fuBKqro8cKCoAmTdyZFxERkeXSCl5CiDMBlEkpZyZw2xFCiFIhRGl5eXk6D+u+Nm30VZ2aGmDxYvcf27TN2LUrkGNROV6iwSveapcNbTGIiIh8kNfQDYQQnwMoMXzpHgB3AxicyANJKUcDGA0AAwYMsH91rHdvYN266LE5c8yrYU6yeZsxLNGWEtxmJCIiitJg8JJSnmgaF0IcBKArgJlCrWB0BDBNCHGYlHKto7P0Q+/ewMSJ0WNe1HnZXFgflu6KFxERUZZqMHjFI6WcDaBd+P+FEMsBDJBSbnBgXv7z6+igoK54MXgRERE1yKLCIcv4dVh2EFa8TIX+DF5EREQNSnnFK5aUsotT92WFXr30sSVL1OHVjRu797hBCF5c8SIiIkoJV7ziadEC6NQpeqyuDliwwL3H3LEDWBtTHpeTA3Tp4t5jpsK04rVmjbryMxKDFxERURQGr/p4Xee1dKk+1rkz0KiRe4+ZiiZNgNato8fq6vTQyOBFREQUhcGrPl4fHRSEwvqwRFpKMHgRERFFYfCqj9dHBwWhvisskTovBi8iIqIoDF718XLFq6IC+PxzfTxIK16xwWvzZv02DF5ERJTFGLzqc8AB+tjy5aoI3kkzZwKHHgp8+qn+NVtXvBJpKWFa8Soudmc+REREAcDgVZ+mTYFu3fTxefOcuf/aWuCRR1ToireF2a+fM4/lNG41EhERJY3BqyFu1XmtWAGccAJw++1AdbX5NldeaV8ribCGgtfu3cCuXdFfz80Fmjd3d15EREQWY/BqiBt1Xq+8AvTpA0yaZP56o0bA448DzzyT3uO4qaHgZarvKi4G1LmeREREWcmxzvUZy8mjgzZvBq69Fnjttfi3Oegg4OWX1X9tZgpeq1erfl45OdxmJCIiMuCKV0OcaqL6xRdqlau+0HXLLcAPP9gfugCgqEjfNqyqAjaEzkhn8CIiItIweDWkRw9VmxSprAzYsiWx76+sBG69VdVzmc4zBNTq0cSJwKOPunsOpNPq225k8CIiItIweDWkcWNzL61Ethtnz1ZXLD72WPzbXHABMGsWcPzxqc/RLwxeRERESWHwSkSydV51dao4fsAAFb5MWrRQtVyvvhrc3lb19fJi8CIiItIweCUimTqvVauAk05S9VpVVebbHHusWuW68ELn5ugHrngRERElhcErEYm2lBg/XhXGf/GF+X7y84GHH1b1XJ07OztHPzB4ERERJYXtJBLRUBPVrVuB668HXnop/n306qW2Fvv2dX5+fjEFr7Iy9V8GLyIiIg1XvBKx//5qtSrS+vVAeblqgtqnT/2h66abgNLSzApdAFe8iIiIksTglYj8fNVWItbllwPHHaeO/zHp0EEdfD1qFNCkibtz9EO84CUlgxcREZEBg1eiTHVe772nQobJ//2fKqA/6SR35+Wn4mI9UO7apXqcMXgRERFpGLwSZarzMikqAp5/Hnj9daB1a3fn5Dch4reUYPAiIiLSMHglyrTiFevII4GZM4FLLsmew6BN243LlgHbtunjLVu6Px8iIiKLMXglqr7glZcHPPigKrTv2tW7OdnAFLxMrTZattSPXiIiIsoybCeRqK5dgb32Atatix7v2VNd0XjIIf7My2+m4GXq1s9tRiIiIq54JSw3F/jzn6O3EK+7Dpg6NXtDF2AOXrNm6WMMXkRERFzxSsoVV6jzF6dNU/VcphYT2cYUvBYu1McYvIiIiBi8kta3b+Y1Qk2HKXjV1upjDF5ERETcaqQ0mdpJmBQXuzsPIiKiAGDwovS0a6eu6mwIV7yIiIgYvChNOTmJrXoxeBERETF4kQNMdV6xGLyIiIgYvMgBDF5EREQJYfCi9DF4ERERJYTBi9LH4EVERJQQBi9KH4vriYiIEsLgRelLZMWLfbyIiIgYvMgBDQWvpk2BggJv5kJERGSxtIOXEOIGIcRCIcRcIcTDTkyKAqakRPXziofbjERERADSPKtRCHEcgCEA+kgpK4UQ7ZyZFgVKfr4KX6tXm7/O4EVERAQg/RWvawA8JKWsBAAp5fr0p0SBVN92I4MXERERgPSDV3cARwshpgghJgkhDo13QyHECCFEqRCitLy8PM2HJesweBERETWowa1GIcTnAEoMX7on9P3FAAYCOBTA60KIblJKGXtjKeVoAKMBYMCAAdrXKeDqaynB4EVERAQggeAlpTwx3teEENcAeDsUtH4QQtQBaAOAS1rZhiteREREDUp3q/EdAMcDgBCiO4BGADakOykKIAYvIiKiBqV1VSOAsQDGCiHmAKgCMMy0zUhZgMGLiIioQWkFLyllFYCLHZoLBRmDFxERUYPYuZ6c0aFD/K8xeBEREQFg8CKnNG4MtG1r/hqDFxEREQAGL3JSvO1GHpBNREQEgMGLnBSvlxdXvIiIiAAweJGTTCtejRoBhYXez4WIiMhCDF7kHFPwatUKEML7uRAREVmIwYucEy94EREREQAGL3JS5876WLt23s+DiIjIUgxe5JwjjwRKYs5TP+ssf+ZCRERkIQYvck6jRsB//wv066eucPzd74CrrvJ7VkRERNZI96xGomgDBwLTpvk9CyIiIitxxYuIiIjIIwxeRERERB5h8CIiIiLyCIMXERERkUcYvIiIiIg8wuBFRERE5BEGLyIiIiKPMHgREREReYTBi4iIiMgjDF5EREREHmHwIiIiIvIIgxcRERGRRxi8iIiIiDzC4EVERETkESGl9P5BhSgH8LPnDxwsbQBs8HsS5Co+x5mPz3F24POc+RpLKXs7cUd5TtxJsqSUbf143CARQpRKKQf4PQ9yD5/jzMfnODvwec58QohSp+6LW41EREREHmHwIiIiIvIIg5e9Rvs9AXIdn+PMx+c4O/B5znyOPce+FNcTERERZSOueBERERF5hMGLiIiIyCMMXh4RQowVQqwXQsyJGDtYCPGdEGK2EOI9IUTz0Hi+EOL50Ph8IcRdEd9zshBioRBiiRDiTj/+LGTm4HO8PDQ+w8lLmMkZST7PjYQQ40LjM4UQgyK+55DQ+BIhxJNCCOHDH4cMHHyOvwy9X88I/Wrnwx+HDIQQnYQQ/wu9/84VQtwUGm8lhPhMCLE49N/i0LgIvU6XCCFmCSH6R9zXsNDtFwshhjX44FJK/vLgF4BjAPQHMCdi7EcAx4Z+PxzAA6HfXwjgtdDvCwEsB9AFQC6AnwB0A9AIwEwAvfz+s/GXc89x6P+XA2jj95+Hvxx5nq8DMC70+3YApgLICf3/DwCOACAAfATgFL//bPzl+HP8JYABfv95+Mv4HLcH0D/0+yIAiwD0AvAwgDtD43cC+Fvo96eGXqcCwEAAU0LjrQAsDf23OPT74voemyteHpFSTgawKWa4B4DJod9/BuDc8M0BNBVC5AFoAqAKwDYAhwFYIqVcKqWsAvAagCFuz50S49BzTJZL8nnuBWBi6PvWA9gCYIAQoj2A5lLK76R6934BwFluz50S48Rz7ME0KQ1SyjVSymmh328HMB/A3lA/U58P3ex57HldDgHwglS+B9Ay9Dr+NYDPpJSbpJSbof5tnFzfYzN4+WsOgDNDvx8KoFPo928C2AlgDYAVAB6VUm6C+kexMuL7V4XGyF7JPseACmWfCiGmCiFGeDlZSlm853kmgCFCiDwhRFcAh4S+tjfU6zeMr2X7Jfsch40LbTP+ntvJdhJCdAHQD8AUAHtJKdcAKpxBrWIC8X/+Jv1zmcHLX8MBXCeEmAq11FkVGj8MQC2ADgC6ArhFCNENaokzFvuB2C3Z5xgAjpRS9gdwSuh7j/F4zpS8eM/zWKg34lIAowB8C6AGfC0HUbLPMQBcJKU8CMDRoV+/9XTG1CAhRDMAbwG4WUpZ365DvNds0q9lX85qJEVKuQDAYAAQQnQHcFroSxcC+FhKWQ1gvRDiG6il65WI/iTVEcBq72ZMyUrhOV4qpVwd+t71Qoj/QoW0ydqdkzXiPc9SyhoAI8O3E0J8C2AxgM1Qr98wvpYtl8JzDCllWei/24UQr0C9ll/wduYUjxAiHyp0vSylfDs0vE4I0V5KuSa0lbg+NL4K5p+/qwAMihn/sr7H5YqXj8JXuAghcgDcC+CZ0JdWADg+dBVFU6hCvgVQxZ37CyG6CiEaAbgAwLvez5wSlexzLIRoKoQoCn1PU6g3+jn6PZNN4j3PQojC0PMIIcRJAGqklPNCWxjbhRADQ9tPlwCY4M/sKRHJPsehrcc2ofF8AKeDr2VrhF53YwDMl1I+HvGldwGEr0wchj2vy3cBXBJ6zx4IYGvodfwJgMFCiOLQFZCDQ2NxccXLI0KIV6FScRshxCoAfwDQTAhxXegmbwMYF/r906Hfz4FaxhwnpZwVup/roZ7UXABjpZRzPftDUL2ceI5D243/DZWC5AF4RUr5sXd/CmpIks9zOwCfCCHqAJQheqvpGgD/gbq44qPQL7KAQ89xQWg8H+r9+nMA//bmT0AJOBLquZothJgRGrsbwEMAXhdCXA71AXlo6GsfQl3ZuATALgCXAYCUcpMQ4gGohREA+FNEva4RjwwiIiIi8gi3GomIiIg8wuBFRERE5BEGLyIiIiKPMHgREREReYTBi4iIiMgjDF5EREREHmHwIiIiIvLI/wMPhP7u5jU8gwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (10,7))\n",
    "\n",
    "ax.plot(100.*np.log(bls.loc[start_year:2007].LA_TFP).diff(), color = \"red\",\n",
    "        linewidth = 5, label = \"Imports\")\n",
    "\n",
    "ax.axhline(y= mean_tfp,   \n",
    "           # Set the value equall to the average\n",
    "           color='k',       # make the color red\n",
    "           linestyle='--',   # the line style\n",
    "           linewidth=4)     # thickness of the line\n",
    "\n",
    "ax.spines[\"top\"].set_visible(False)\n",
    "ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "ax.set_xlim(1977,2000)\n",
    "#ax.set_ylim(5,15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a name=\"real_rate\"></a> Real Interest Rates \n",
    "\n",
    "We construct the real interest rate as the difference between the rate of return on a one year U.S. treasury constant maturity nominal bond and the realized inflation rate. We use the U.S consumer price index, all items as our measure of realized inflation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "codes = [\"GS1\", \"CPALTT01USM659N\"] \n",
    "\n",
    "real_rates = data.DataReader(codes,\"fred\",1977)  # Read data from FRED and specify start date\n",
    "real_rates = real_rates.resample(\"Y\", convention= 'start').mean()\n",
    "\n",
    "real_rates = real_rates.dropna()\n",
    "\n",
    "real_rates[\"real_rate\"] = real_rates[\"GS1\"] - real_rates[\"CPALTT01USM659N\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "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>GS1</th>\n",
       "      <th>CPALTT01USM659N</th>\n",
       "      <th>real_rate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DATE</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1977-12-31</th>\n",
       "      <td>6.088333</td>\n",
       "      <td>6.498055</td>\n",
       "      <td>-0.409722</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1978-12-31</th>\n",
       "      <td>8.340000</td>\n",
       "      <td>7.615543</td>\n",
       "      <td>0.724457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1979-12-31</th>\n",
       "      <td>10.665000</td>\n",
       "      <td>11.224806</td>\n",
       "      <td>-0.559806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-12-31</th>\n",
       "      <td>12.048333</td>\n",
       "      <td>13.575444</td>\n",
       "      <td>-1.527111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1981-12-31</th>\n",
       "      <td>14.778333</td>\n",
       "      <td>10.352041</td>\n",
       "      <td>4.426293</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  GS1  CPALTT01USM659N  real_rate\n",
       "DATE                                             \n",
       "1977-12-31   6.088333         6.498055  -0.409722\n",
       "1978-12-31   8.340000         7.615543   0.724457\n",
       "1979-12-31  10.665000        11.224806  -0.559806\n",
       "1980-12-31  12.048333        13.575444  -1.527111\n",
       "1981-12-31  14.778333        10.352041   4.426293"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "real_rates.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average Real Rate 2.825\n"
     ]
    }
   ],
   "source": [
    "avg_real = real_rates[\"real_rate\"].loc[(real_rates.index < \"2001\") & (real_rates.index > \"1977\")].mean()\n",
    "# This has the real rate starting in 1977 (so inflation is from 76 to 77), then ending in 2000\n",
    "# so inflation is going from (99 to 00)\n",
    "\n",
    "print(\"Average Real Rate\", round(avg_real,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAGbCAYAAAAV7J4cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUVf7H8c8hoYoIIipFRcQCi4WyYpe1C6hYECwoouiuXbGjYq+sBRcVRFGBFRQLVtC1dwRkFUUQy0qTKk0QCDm/P074mcy9E5LMnXOnvF/Pk8d47uTOF4Yknzn33O8x1loBAAAg/arFXQAAAEC+IHgBAAB4QvACAADwhOAFAADgCcELAADAk8I4nvSoo46y48ePj+OpAQAAKstEdaJYZrwWL14cx9MCAADEikuNAAAAnhC8AAAAPCF4AQAAeELwAgAA8ITgBQAA4AnBCwAAwBOCFwAAgCcELwAAAE8IXgAAAJ4QvAAAADwheAEAAHhC8AIAAPCE4AUAAOAJwQvwxVrps8+kyy6T9tpLatdOGjs27qoAAB4Vxl0AkNOslaZOlUaPlp59Vvr557LHe/aUPv9cat8+lvIAAH4RvIB0+PZbF7bGjJFmzkz+uA0bpPvuk0aN8lcbACA2xlrr/Uk7dOhgJ02a5P15gbSaNcsFrTFjpK+/rvjX1awpzZsnbbll+moDAKTCRHUiZryAVPzyi7uEOHq0NHly1c6xdq00cqR08cXR1gYAyDgEL6Cyfv1Veu45N7P18ceV+9qCAqlJE2n27LLjjz0mXXSRZCJ7UwUAyEAEL6AiFi+WXnjBzWy9/75UXFzxrzVGOvBAt5D+xBOlJUuk1q3LPmbaNGniRKljx2jrBgBkFIIXkMzy5dJLL7mw9Z//SEVFlfv6jh1d2OreXWra9M/xrbeW9t8/OFv22GMELwDIcQQvoLQ1a6SXX5aeeUZ64w1p3brKfX3btlKPHtLJJ0s77pj8ceecEwxeo0dL998vbb555esGAGQFghewYYO7fDhypGtounJl5b6+VSs3s9Wjh7TrrhX7mu7dpUsukVas+HPs999d+Orbt3LPDwDIGgQv5K+vv3Zha9Qoae7cyn3tTju5oNWzp9SmTeUXxW+2mXTqqdKjj5YdHzaM4AUAOYzghfwyd667jDhihPTVV5X72u22c5cQe/Z0neZTvQOxb99g8Jo40dW1xx6pnRsAkJEIXsh9K1e6OxJHjJDeecdt41NR22zjLgv27Cntu69ULcLtTdu1cx9TppQdHzZMGjQouucBAGQMOtcjN61fL731lgtb48a5RfMVtfnm0kknSaedJnXq5Hpvpcsjj0jnn192rH5918m+du30PS8AoDIia7JI8ELusFaaNMmFrdGjpUWLKv61hYXSUUdJp58uHXOMVKdO+uosbflyqXHjYDAcMcLVAgDIBGwZBPy/n35yi+RHjix/Q+owHTu6gNOjh9SoUXrqK88WW7h1Y089VXZ82DCCFwDkIIIXstPSpW6PxJEjK79tT4sWLtScfrq0887pqa8y+vYNBq/333chcpdd4qkJAJAWBC9klwULpH79XOhav77iX7fllm5Wq1cvaZ99MmtPxP32c73Apk8vO/7449Ldd8dTEwAgLSK8RQtIM2ulU05xfbcqErpq1nSL5MeNk+bPlx5+2N2ZmEmhS3L1nHNOcPzJJyvfOR8AkNEIXsgekyZJ77676ccdfLBbI/Xrr9Jzz0nHHivVqJH++lLRq5dUvXrZsYULpVdfjaceAEBaELyQPYYOTX6sdWvpzjul//1Peu896eyzXVuGbNGokXT88cHxxx7zXwsAIG0iCV7GmPrGmLHGmO+MMdONMftGcV7g/61c6TrOJ+rVyzUgnTZNuuYaafvt/dcWlbDLjRMmuDAJAMgJUc14PShpvLV2N0l7Spq+iccDlfPMM24T6dLq1XMNSNu2zbx1W1Vx6KHSjjuWHbNWGj48nnoAAJFLOXgZY+pJOkjS45JkrV1nrV2W6nmBMsIuM552mttsOldUq+YukSZ64glpwwb/9QAAIhfFjFcLSYskDTfGfGmMGWaMCfw2NMaca4yZZIyZtKgyHcWBKVOkyZOD4+ee67+WdOvdO7gf5OzZ0ptvxlIOACBaUQSvQkntJD1irW0r6XdJ1yQ+yFo71FrbwVrboVEcHcKRvcIWmHfoIO21l/9a0q1pU6lLl+A4i+wBICdEEbzmSJpjrf285P/HygUxIHWrVrm+XYlycbZro7BF9q+84tpjAACyWsrBy1r7q6TZxphdS4YOlfRtqucFJLkO9StXlh2rW1fq2TOeenzo3NltnF1aUVFwWyEAQNaJ6q7GiySNMsZ8JWkvSXdEdF7ku7BF9aecIm2+uf9afCkslM46Kzg+bJi7yxEAkLUiCV7W2qkl67f2sNZ2s9b+FsV5kee++kr6/PPgeC5fZtwo7O7GWbPc5tkAgKxF53pkrrAF5W3bSu3b+6/FtxYtXF+vRCyyB4CsRvBCZlq9WhoxIjjet29uNEutiL59g2PPPy8tXeq/FgBAJAheyEzPPSctX152rE4d6dRT46knDt26SQ0blh1bu1YaOTKeegAAKSN4ITOFXVLr2VPaYgv/tcSlZk3pjDOC4489xiJ7AMhSBC9knm++kT7+ODgedukt14X19Jo2TZo40X8tAICUEbyQecJmu3bfXerY0X8tcWvdWtpvv+D4sGH+awEApIzghczyxx/S008Hx889N38W1ScKm+l75plgY1kAQMYjeCGzPP+89FtCG7hataTTTounnkzQvbtUr17Zsd9/l8aMiaceAECVEbyQWcI61Z98stSggf9aMsVmm4XfzUlPLwDIOgQvZI4ZM6QPPgiO50On+k0JW2Q/caLr7g8AyBoEL2SOsBmcVq3CF5fnm/btXdf+RCyyB4CsQvBCZli7VnryyeB4Pi+qTxS2yH7ECGnNGv+1AACqhOCFzPDSS9KSJWXHataUevWKp55MdOqpUu3aZceWLZNeeCGeegAAlUbwQmYIW1R/4onBLXPy2RZbuBsNErHIHgCyBsEL8Zs1S3rnneA4i+qDwhbZv/++NHOm/1oAAJVG8EL8whaI77KLdNBB/mvJdPvvL+22W3D88cf91wIAqDSCF+K1bp00fHhwvG9fFtWHMSZ81uvJJ93fJQAgoxG8EK+XX5YWLiw7Vr26dOaZ8dSTDc44w/0dlbZwofTqq/HUAwCoMIIX4hW2MPyEE6RGjfzXki0aNZK6dQuOs8geADIewQvx+ekn6c03g+Nh/apQVtjf0YQJ0i+/+K8FAFBhBC/EJ2xR/U47SX/7m/9ass2hh0rNm5cds1Z64olYygEAVAzBC/FYvz75ovpq/LPcpGrVpLPPDo4/8YS0YYP/egAAFcJvOMTjtdek+fPLjhUWSr17x1JOVjrrrGBInT07/PItACAjELwQj7BO9ccdJ22zjf9aslXTplLnzsFxNs4GgIxF8IJ/v/wijR8fHKdTfeWFLbJ/+WVpwQL/tQAANongBf8ef9wtBC+teXPpsMNiKSerde4sNW5cdqyoyDVUBQBkHIIX/CoqCt/e5pxzWFRfFYWFbq1XomHDguEWABA7ftPBr/Hjpblzy44VFISHB1RMnz7BsVmz3ObZAICMQvCCX2GL6rt2lZo08V9LrthpJ9fXKxGL7AEg4xC84M+cOa6NRCIW1acubJH92LHS0qX+awEAJEXwgj/Dh0vFxWXHtttOOvLIeOrJJd26SQ0blh1bu1YaNSqeegAAoQhe8GPDhvBLX2ef7dZ4ITU1a0pnnBEcf+wxFtkDQAYheMGPN98MbuBcrVr4wnBUzTnnBMe+/lp66CFp3Tr/9QAAAghe8OOxx4JjnTu7S42IRuvW0n77BccvuURq2dIFsDVr/NcFAPh/BC+k3/z5rpt6orAF4UhN2KyX5PZwvPhi16j27rulFSu8lgUAcAheSL/hw90ar9KaNAnfZxCp6dFD2mGH5McXLpSuucY9ZsAAackSf7UBAAheSLPi4uSL6gsL/deT6+rUkd5+e9N3ii5bJt1yiwtgV14p/fqrn/oAIM8RvJBeb78t/fRT2TFjXPBCeuy0k9shYOJE12aiPL//Lg0c6C5BXnih9L//eSkRAPIVwQvpFdap/sgjy78chmj89a/Siy+6OxtPPbX8vTDXrpUGD3aL8Pv0kWbO9FcnAOQRghfSZ8EC6aWXguN0qverTRvXSHXGDLf4vnr15I8tKnJr8lq1knr2lL76yl+dAJAHCF5In6eecr/IS9tmG7c3I/xr2dK19fjhB3eHY61ayR9bXCyNGSPtuad07LHS55/7qxMAchjBC+lhbXjvrj59yp9xQfptt5304INuPdc110ibb17+4195RdpnH+mww6R336UTPgCkgOCF9HjvPWnWrOB4sj5T8G/rraU773QB7JZbpC23LP/xb78tHXKItP/+brNzAhgAVBrBC+kRtqj+sMOkFi3814LyNWgg3XCDC2ADB0rbblv+4z/91F0ubtdOGjfOT40AkCMIXoje4sXSCy8Ex1lUn9nq1pX69XPtPx5+eNN3nk6d6tpVDBzopz4AyAEEL0Tv6aeDmzI3aiQdd1w89aByatWS/vEP6fvv3R2Ou+xS/uNvuklautRLaQCQ7QheiJa14ZcZe/eWatTwXg5SUL26e92+/VZ69ll3h2OY33+XHnnEa2kAkK0IXojWhx+6flGJWFSfvQoKpO7dpS+/lF59VWrbNviYQYOkNWv81wYAWYbghWiFtZD42982fbkKmc8YqUsX1w2/oKDssYUL3SVmAEC5CF6Izm+/Sc89Fxzv29d/LUifHXZwXe0TDRwobdjgvx4AyCIEL0Tno4/cnn+lNWwoHX98PPUgfa68Mjg2axbtJQBgEwheiM7PPwfHjjmm/K1pkJ323NNtdp7o7rtprAoA5SB4ITqzZwfHmjf3XgY8CZv1mjjR3WABAAhF8EJ0woLXdtv5rwN+HHKI616f6J57/NcCAFmC4IXoELzyizHSVVcFx197TZo2zX89AJAFCF6Izpw5wTGCV2478URpxx2D42wjBAChCF6IxoYN0ty5wfFmzfzXAn8KC93+jolGjQoP4gCQ5wheiMaCBVJRUdmx+vXdxsvIbWed5dqGlFZUJD34YDz1AEAGI3ghGqzvyl916kgXXRQcHzJEWrbMfz0AkMEIXogGwSu/XXCBVLt22bGVK134AgD8P4IXosHC+vy21VZSnz7B8QceCO5mAAB5jOCFaDDjhcsvl6ol/Ej59Vdp5Mh46gGADETwQjTCghd3NOaXFi2k7t2D4/feKxUX+68HADJQZMHLGFNgjPnSGPNqVOdEFmHGC1L4NkIzZkivvOK/FgDIQFHOeF0iaXqE50M2IXhBktq3lw49NDjONkIAICmi4GWMaSapi6RhUZwPWaaoSJo/PzjOpcb8FLaN0CefSB9/7L8WAMgwUc14PSDpKklJF3IYY841xkwyxkxatGhRRE+LjDB/fnANz1ZbBdsLID8cfri0557B8Xvv9V8LAGSYlIOXMaarpIXW2snlPc5aO9Ra28Fa26FRo0apPi0yCZcZUZox4Wu9xo2TvvvOfz0AkEGimPHaX9KxxpifJY2WdIgxhvvH8wl3NCLRySdL228fHGfzbAB5LuXgZa291lrbzFrbXFJPSe9Ya09PuTJkD2a8kKh6ddfXK9GIEdK8ef7rAYAMQR8vpI6u9Qhz9tlSgwZlx9atkwYNiqceAMgAkQYva+171tquUZ4TWYAZL4SpW9ft4ZjokUekFSv81wMAGYAZL6SO4IVkLrpIqlmz7NiKFdLQofHUAwAxI3ghdSyuRzJbby317h0cf+ABd9kRAPIMwQupWbdOWrAgON60qf9akJn69XMtJkqbO1d65pl46kmnFSukDz5wwbJ3b+mkk6RRo+KuCkAGKYy7AGS5uXMla8uObbNN8PIS8tfOO0snnCA9/3zZ8XvvlXr1kqpl6fu/RYukKVOkL7/887+zZgUf9/zz0urVUt++/msEkHEIXkgNdzSiIq68Mhi8vvlGeuMNqUuXeGqqKGvdv/PSIWvKFPemo6KuvtrNfiXe5Qkg7xC8kBoW1qMiOnaUDj5Yev/9suP33JNZwau42M1alZ7FmjJFWrIktfP+9pt0xx1smwSA4IUUsbAeFXXVVcHg9cEH0mefSfvs47+e9eul6dPLBqypU6VVq9LzfIMGufYazZun5/wAskKWLq5AxmDGCxV19NFSmzbBcd+zQOvWSTff7O643HNP6ayzXCj66KNoQtfOO7stkxI3iV+3TurfP/XzA8hqBC+khuCFijJGuuKK4PiLL0rff++nhhkzpP32k266SVq2LLVzFRRIe+whnXmm9OCDbvZu+XJp5kxpzJjwLZP+/W9p0qTUnhdAVjM28Y40Dzp06GAn8cMnN7Rv7y7RlPbRR9L++8dTDzLbunVSixbBhennnSc9+mj6ntda17T18svdHYaVVauWC1nt2klt27r/tmnjxpNZsUJq2dLd/Vhap07SO+8EW2wAqLyVK6Wff5aaNJEaNkznM0X2DcuMF1LDjBcqo0YN6bLLguNPPhneDy4KixZJxx0n/f3vFQtd9eq5GwEuvVR6+mnp66/dD/fPP3fbHZ17rtShQ/mha+N5bropOP7ee9Jrr1XlTwJAkn791b2R6txZ2mor96Zohx2khx6Ku7IKYcYLVffHH8F1LMZIa9dK1avHUxMy34oV0vbbu8typfXvL912W7TP9frrUp8+yUNdQYF06KFuBmvjbFaLFtH1Flu/3s2MzZxZdrxVK+mrr6RC7m8CKmTWLLcs4aWXpE8/DfaP3Oiqq6S77krHjHJkJyR4oepmzXILiUtr0qRy/Y2Qn6691v1wLK1+fTeDWrdu6udfvdr9AB48OPljWrZ0XeX33jv15yvPSy9Jxx8fHB8yxM2eAQiy1i1jeeklF7i++abiX3vmmdJjj0U9AcClRmQALjOiqi6+2F12LG3ZMmnYsNTP/eWX7lJgeaGrb1/3uHSHLsld5jzggOD4jTemr3UFkI2Kitz6x4svdpcOO3Rws+CVCV2S9NRT7s1OVdZzekDwQtXRtR5V1bixdMYZwfH77nOX56piwwbXkLVjR9efK0zDhu7d89Ch0cysVYQx0sCBwfEFC8LHgXyyerWb1TrzTLfd3KGHurVaYW/sK+O116TDDpOWLo2mzggRvFB1zHghFf36Bcdmz3atGCpr9mz3Q/bqq5MHtyOPdAvlu3Wr/PlT1bGj6+2V6N57pfnz/dcDxGnJEjcr1a2bWxx//PHuRpbKhKRddnHf75995vryhfn0UzfbnGqIixjBC1VH13qkYrfd3GW4RPfem3zhbJjRo91dTe+9F368Zk3XIPX1191MW1zuuCO45mT1amnAgHjqAXz65Rf3fXjIIW5mq3dvadw4ac2aip9j773d99G337qefHfd5d7U3Hija0cTdlPM9Omud1+yWfAYELxQdcx4IVVXXRUc++or6c03N/21y5dLvXpJp5ySvBnqHntIkydLF10U3Z2KVbXTTtL55wfHH3+88mtYgGxQVCT961+u3+MOO0iXXCK9+65bFlARhYVuJnvwYPf75vPP3Y05rVoFH3veedJzzwXXjkpuWcwBB7gZsAxA8ELVEbyQqv32C2+2e8895X/dRx+57X5Gjkz+mCuukCZOlP7yl9RqjNINN0hbbFF2rLjYXTIBcskff0hHHeXe9CQ22S5PnTrSiSe67+2FC6W33nJvWCpyNeWEE6QJE1wPvURLl7r1YxnQQ4/ghapjcT2iEDbr9c474VvrrF/v+n0dfLD0v/+Fn69pU+ntt90ly5o1o601VQ0bStddFxx/7TU3EwDkgvXrpe7d3fdhRWy1leu39/LL0uLF0tix0mmnSQ0aVP65O3Vy23dtu23w2Jo1bnnD009X/rwRoo8Xqmb1ammzzcqOFRS45qkFBfHUhOxUXOxmpb77ruz4ySeXXWg/c6b7YVzez47u3d1ajy23TE+tUfjjD2nXXd2al9LatZO++CL+S6JAKjZscN+nm7pJpnlzt6i+Wzc36x31740ff3Q31MyaFX78nnukK6+szBnp44WYhV1mbNKE0IXKq1YtfPPssWOlH374c5/Ftm2Th666dd22Q2PGZHboktxWQ7ffHhyfMkV65hn/9QBRKS52a62Sha4993Q3k0yd6oLRffdJBx2Unt8bLVpIH3/s1peFueoq93OnuDj6594EZrxQNf/5j3T44WXH9tvP/UMHKmvtWvcO+Ndfy46fdpprMjpuXPKv3W8/acQI94M2WxQXS3/9a3Dty/bbu7u1NrUPJPLD+vXSJ5+4NxO77x53NeWz1u3D+uCDwWPGSP/+t9Szp/+6Vq50M2vJLnuefrr0xBMV6XLPjBdixsJ6RKlmTbcpdaJRo5KHroIC6ZZbpPffz67QJblZvnvvDY7/8kvWbPSLNHvjDXcnbKdO7u7cnj0ze6eDG28MD12S25EijtAlSZtv7tZQ9ugRfnzkSOnYY6Xff/dWEsELVUPwQtTOO8/9kKyIli3d7OoNN2TvRtOHHCJ17hwcv/1212AS+WnpUtfFvXPnsj9nx4yRDjww45qBSnLrpZJtcP/AA27hfJxq1nQzbhdeGH58/Hj3/bh4sZdyCF6oGu5oRNTq13fha1POPtvts9ixY/prSrd77gkupl++PPkvMeS2F16QWrdOftfd1KnuEvVnn/mtqzyPPJK8Hcqtt7reXZmgWjXXwDXZ99bEia7XV7K7paMsJe3PgNxE13qkwyWXJF9r0bCh+8U0bJi/fRbT7S9/CZ8NGDzY3ViA/LBggbuL98QT3eebemynTm4GJ24jRoQ3BZbc4vX+/f3WsynGuJqGDg2/e3jGDHeH5bRpaS2D4IWq4VIj0qFZM7eVSKIjjnAd7Y8/3ntJaXfLLa5pZGkb+5Uht1nr1jG2bu26rlfU2rXuxpP+/WO5K0+SexMU9r0qSf/4h9vOx0S2Hj1afftKzz8f3udv7lx3Sfejj9L29AQvVA3BC+nywANSly5u8XyzZm7B7htvuHYluahx4/B2GmPGuMsfyE1z57pF3aefnnxz6C22cGv+kv3bv+MO6aST/C+6nzDBLZYPC329erltgjI1dG3UrZvbmixxJwnJbUF2+OHSK6+k5alpJ4HKW7Ei+I+1enXXGJLmj4jKqlVuJigf/k2tWuVuGEi8zHTgge6uzUz/JYaKs9btz9mvn/tZmswxx7hmwE2aSPPmuaDwxRfhj91rL9f13ceb3w8+cFsBhW1ufcIJ7g1DNt3w8tVX7s8zf37wWEGBuyzplgPQTgIxCltY36xZfvyChD916+bPv6m6daWbbw6Of/ih+4WK3PDTT24mpW/f5KGrYUO3fmvcuD9nupo0cQE8WUsEX4vuv/hC6to1PHQdeaSrO5tCl+RadXzyibTLLsFjGza4m3nuvDPSp8yTn2qIFJcZgeidfbbUqlVw/Oqr3ZovZK/iYtefrU2b8vcv7NFD+vZb6ZRTgrOctWu7nQ3CArr056L7UaMiK7uMadPczNDKlcFjBx7o1nxl2t6oFdW8uVvT1aFD+PGw/VVTQPBC5XFHIxC9wkLp7ruD4zNmuDs5kZ1mzHDb4lx8sdvjNsy227rgMnq0tPXWyc9ljGtU+uyzLoglWrvWrRmLetH999+7mbqwtWgdOkivvhq8QSTbNGrkNqo/4oi0PxXBC5XHjBeQHl27SgcfHBy/6abwmQZkrqIiF6T33LP8rdR693azXJW5Y7d7d3cZumnT8ONRLrr/5RfpsMOC23lJrh3K+PFSvXqpP08mqFvXLag/5ZS0Pg3BC5VH8ALSwxhp4MDg+MKFrtkqssPXX0v77CNdc42bhQqz3Xbubt3hw6UGDSr/HO3bu7te//rX8OMvvugagv7yS+XPvdGvv7rQFXaOli2lt95ya9JySY0abhuhdDZ+tdZ6/5Bkq/LRrl07m0y7du2qdE73VxCub9++VT7npEmTQs85ZMiQKp9zyJAhoeecNGlSlc/Zt2/fpH9+XideJ16nDHidJGsla2vXtnbOHF6nTH2d2rWzdu1aawcMsLZ6dfealXy0q+I5eZ0y5+eejTADMeMFANlgzRq3vgeZafVqt97p5pu5GQLlIngBQLYYPtxdxkLm+e47XhtUCMELALKFtW4PPGSPVq2kXXeNu4qykt1dWVTkt448Red6VM5XX7m7dEpr2dLdbgwgWvfd5zqcJ3rrLbfoGfF58UV352Cytg0FBW5x/Q03+O9vNX++dNxxyTvd77mna8y7/fZ/jhUVubv5xo4NPr5WLXf3Ytgdt/mDzvWICXc0Av5ccIG0447B8SuvjG9zZLiZx6uvTv4a7LWXCz233RZPU9HGjV2n+549w4//97/S3ntLn37q/r+4WDrnnPDQVb266zGW36ErUgQvVE7YdkEELyA9atZ0PZkSTZ3qbnlHPKZPD5/lr1HDha2JE6W2bf3XVVrt2m4Ln1tuCT++sdP9iBGuuetTTwUfU62aO8fRR6e11HxD8ELl0LUe8KtHj/BeTddfH75nHtLvlVeCYzvtJE2Z4rrGV6/uv6YwxrhLnWPHhne6X7dOOuMMafDg8K9/4gl3ORWRInihcrjUCPiVrKnq7NnSgw/6rwfhweu881wn90x04oluL8Jkne7DPPSQdOaZ6aspjxG8UDkEL8C/gw6Sjj02OH7nndKiRf7ryWeLF/+5Nqq0Y47xX0tltGvn1p0l63Rf2h13SBdemP6a8hTBC5VD8ALicffd7k650laskG69NZ568tXrrwcX1bdsmXktI8JsatG9JF17rftA2hC8UHHWsrgeiMtuu0l9+wbHH3mEdi4+vfxycOyYY9wl4WywcdF9WGC/6CLp9tv915RnCF6ouCVLpD/+KDtWp45Uv3489QD55qabpLp1y44VFUmXXebeGCG91q6VJkwIjmf6ZcZExribM15/3W3m3aqVdP/90gMPZE+AzGIEL1RcssuMfKMCfmyzTXjn+tdec784kV7vvy+tWlV2bIstpAMOiKeeVB19tFuv9u230qWXuvYRSDv+llFxrO8C4nf55dCncywAAB+iSURBVG6tTqKrr5Y++cR/Pfkk7G7Gzp0zp30EsgLBCxVH8ALit9lm0sMPB8eLilzPr8WL/deUD6wND17ZdpkRsSN4oeJYWA9khm7d3LquRHPmSKefznZC6TBtmvS//5UdKyiQjjoqnnqQtQheqDi61gOZ46673MLoRBMmhG8zhNSE3c144IFSgwb+a0FWI3ih4rjUCGSOGjWkZ5+VGjYMHhswQHrnHf815TIuMyIiBC9UHMELyCzbbec2OU5UXCydeqo0f77/mnLRggVu4+tEBC9UAcELFVNczBovIBMdfbTbmDnRggXSKae4RfdIzWuvBfuk7babtPPO8dSDrEbwQsUsWiStX192rF499wEgXjfdJHXqFBx//33pxht9V5N7uMyICBG8UDEsrAcyV2Gh9MwzrsFqojvvdB3KUTV//CG9+WZwnOCFKiJ4oWJY3wVktm23lUaPDu8+3quX9Msv/mvKBe+8I61eXXZsyy2lffeNpx5kPYIXKobgBWS+Tp3CNz9eulQ6+WRp3TrvJWW9ZN3qCwv914KcQPBCxRC8gOxwzTVuwX2izz8P3+cRyVkrvfpqcJzLjEgBwQsVQ/ACskO1aq7FRNj354MPSs8/77+mbDV1avBu7urVpSOPjKce5ASCFyqGVhJA9mjYUBozJvxyWJ8+0qxZ/mvKRmGXGQ8+WNpiC/+1IGekHLyMMdsZY941xkw3xnxjjLkkisKQYbirEcgu++4r3XtvcHzFCql7d2nNGv81ZZuwbYK4zIgURTHjVSSpn7W2laR9JF1gjGkdwXmRKTZskObODY4z4wVktksukU44ITg+dao7huTmzZMmTw6OE7yQopSDl7V2vrV2SsnnKyVNl9Q01fMig/z6qwtfpTVoIG22WTz1AKgYY6QnnpBatAgee+yx8O2G4IQtqv/LX6Qdd/RfC3JKpGu8jDHNJbWV9HmU50XMWFgPZK8ttpDGjpVq1gwe+/vfpW++8V9TNqBbPdIksuBljKkr6XlJl1prV4QcP9cYM8kYM2nRokVRPS18YGE9kN3atpUGDQqOr17t1nutWuW/pky2erX0n/8Ex4891n8tyDmRBC9jTHW50DXKWvtC2GOstUOttR2stR0aNWoUxdPCFxbWA9mvb1/ptNOC49Onu5mvxE2g89nbb7utgkpr1Ejae+946kFOieKuRiPpcUnTrbX3pV4SMg6XGoHsZ4z06KNSq1bBY6NGuTVfcMLuZuzSRSoo8F8Lck4UM177S+ol6RBjzNSSj84RnBeZguAF5Ia6dd16rzp1gscuvliaMsV/TZmmuJhu9UirKO5q/Mhaa6y1e1hr9yr5eD2K4pAhCF5A7mjdWhoyJDi+dq1b77V8uf+aMsnkye5O7tJq1JCOOCKeepBz6FyPTWNxPZBbTj/drflK9OOP0lln5fd6r7C7Gf/2NzdbCESA4IXyFRVJ8+cHx5vSqg3IaoMGSXvtFRx/8UW3p2O+Cgte3M2ICBG8UL5589yah9K22kqqXTueegBEo1Yt6bnnpHr1gseuvFL69FP/NcVt9mzX1T9R167+a0HOInihfKzvAnJXy5aus32ioiLp5JOlxYv91xSnsNmuPfeUtt/efy3IWQQvlI/gBeS2E08M37dxzhzpjDOCM965jG718IDghfIRvIDcd889UseOwfE33pDuust/PXFYtUp6553gOMELESN4oXxhdzTStR7ILTVqSM8+K225ZfDYDTdI777rvybf3npLWreu7Ni220odOsRTD3IWwQvlY8YLyA/bby+NGBEcLy6WevSQZs3yX5NPYZcZu3aVqvFrEtHiXxTKR/AC8kfnztK11wbHFy2SjjxSWrDAf00+bNhAt3p4Q/BC+QheQH655Rbp4IOD4z/+6ILZypX+a0q3iRNduCytVi3psMPiqQc5jeCF5NauDX+HS/NUIHcVFkqjR0s77BA8NmWKdMIJwbVQ2S7sMuOhh4bvaQmkiOCF5ObNC45ts41biAsgd227rTRhgtSwYfDYf/7jthXKpTYTtJGARwQvJMdlRiB/7bqrW/cUtkvFv/8tXXWV/5rS4eefpWnTguN0q0eaELyQHMELyG/77OPaTBQUBI/985/uI9uFzXa1b8+SCqQNwQvJEbwAdO0qPfZY+LErrpBGjfJbT9S4zAjPCF5IjuAFQHJrum6/PfxY797Sm296LScyK1ZI770XHCd4IY0IXkiOrvUANrr2WumCC4LjRUVuv8fJk/3XlKoJE6T168uONW0qtW0bTz3ICwQvJMeMF4CNjJEefFA66aTgsVWrpKOPzr7u9sm61RvjvxbkDYIXkiN4ASitoMBtKxTWYHXRIumoo7Knu/2GDdLrrwfHjz3Wfy3IKwQvhFuzRlq8uOxYtWpSkybx1AMgM9SqJb30krT77sFjP/yQPd3tP/1UWrKk7FidOtIhh8RTD/IGwQvhwtZ3NW7suloDyG/160tvvOE21k40ZYpb85Xp3e1ffjk4dvjhLlgCaUTwQjgW1gMoT9OmbnH6llsGj731ltSnT2Z3t6eNBGJC8EI41ncB2JTddkve3X7UqMztbj9rlvTdd8HxLl3814K8Q/BCOIIXgIrYd9/s624fNtu1995uj0ogzQheCEfwAlBRXbtKQ4aEH7viCre3YyYJC17czQhPCF4IR/ACUBlnny3ddlv4sd693bqvTPDbb9IHHwTHWd8FTwheCEfwAlBZ110nnX9+cHz9eumEEzKju/348a6HV2nbbx/eHgNIA4IXwnFXI4DKMkYaNMiFrESrVrkeXz/84L+u0pLdzUi3enhC8ELQ77+76fjSCgtZeApg0woK3B2NBx0UPLZwoXTkkfF1t1+/3vUfS8RlRnhE8EJQ2GXGJk3C71oCgES1aknjxklt2gSP/fCDa9sQR3f7jz+Wli0rO1a3rtSpk/9akLcIXghifReAVNWv79ZThf3smDzZbbbtu7t92GXGI4+Uatb0WwfyGsELQQQvAFEor7v9m2/67W5vbfg2QVxmhGcELwSxsB5AVFq1Kr+7/dVX+6ljxgzXsb40Y9yCf8AjgheCmPECEKV995XGjJGqhfzKGThQuusuNyOVTmGXGffdV2rUKL3PCyQgeCGI4AUgascck7y7/bXXSgccIE2cmL7nZ1NsZAiCF4IIXgDS4ZxzpFtuCT/2ySdSx47S6aeH/wxKxZIl7o7GRGwThBgQvBBE8AKQLtdfL/3jH8mPjxol7bqrNGCA6ykYhTfeCC7ib9HCrT8DPCN4oawVK4L9dapXZx0EgGgYIz30kLujMZk1a9zM2C67SE89lfqdj3SrRwYheFXGb7+5d2oHHijdcYe/26B9CpvtatYsfFEsAFRFQYH0+OOuz1fr1skfN2+e22B7772lDz+s2nOtW0e3emQUfptWVFGRa7T36KPSRx9J/fu7BaG5hsuMAHw58kjpv/+VHn5Y2mqr5I+bPNltQdS9u/TTT5V7jg8+CM7i16vn3kADMSB4VdTDD0tffFF27L77pB9/jKeedCF4AfCpsNBdSfj+e+mKK9zShmTGjpV22831/lqxomLnD7vMeNRRUo0aVasXSBHBqyIWLpRuvDE4XlQk3Xyz/3rSieAFIA7160v33it9+610/PHJH7dunXTPPdLOO0tDh0obNiR/rLXhwYu7GREjgldFXHuttHx5+LERI9wPilxB13oAcWrZUnrhBendd6W2bZM/buFC6bzz3GPefjv8Md9+G7w0WVAgHX10dPUClUTw2pSJE6Unnkh+3Nrw2bBsxYwXgEzQqZNb3vHEE9K22yZ/3NdfS4cd5maxZs4seyxstmv//cP3jgQ8IXiVp7hYuvDCTT/u+efd4s9cQPACkCkKCqSzznKBqn9/qVat5I995RXpL3+RLr1UWrrUjbEpNjIQwas8Tz4ZXFAvuR8Gia6/Pu3lpJ21BC8AmWfzzaXbbpO++07q2TP544qKpAcfdOu/7rxT+uyz4GMIXogZwSuZZcuka64Jjrdp476xE40fX/U+M5nit9+k1avLjtWqJTVsGE89AFDaDjtIzzzz5/ZCySxdKl13XXDj7Z13dl3xgRgRvJIZMEBatCg4/tBD0rnnuu0mEoV9o2eTZM1T6e4MIJPsu68LXyNHVu7mH+5mRAYgeIWZNk0aPDg43qOHW/BZvXp4G4mPPpImTEh7eWkTdkcjlxkBZKJq1aTTTpNmzHDbC9Wps+mv4TIjMgDBK5G10kUXBXvD1KkjDRz45/+fckr4Vhf9+2fvrBfruwBkmzp1pBtucA1Ye/dO/rgGDdwdjUDMCF6Jnn1Weu+94Hj//mWntAsK3GLPRFOmuB402YjgBSBbNWkiDR8uTZoUvh3Quee6LvlAzAhepf3+u9uyIlHLllK/fsHxbt2kDh2C4zfcUH435UxF8AKQ7dq3l95/320vdNBBUqNGLnTddFPclQGSCF5l3XFH+DqnBx6QatYMjhsj3X57cHz6dGnUqOjrS7dki+sBIJsYI514ogtgCxdKQ4aU3wMM8IjgtdGsWWXXcG3UpYv7SObww927qkQDBrg9xbIJi+sBAEgrgtdGl10WDEo1arjZrvIkm/X6+Wfp8ccjKy/trCV4AQCQZgQvSXrtNenVV4PjV1zh1ndtygEHhG+6euutwYakmWrxYumPP8qObbaZVL9+PPUAAJCDCF5r17q9vRI1a+YaolZU2B2O8+dLDz9c9dp8SrawnuapAABEhuB1331ufVeigQPdjE9FtWsnde8eHL/zTmnFiqrX5wt3NAIAkHb5HbzmzAmfqerUSTr55Mqf75ZbXDfl0pYule6/v0rleRW2vos7GgEAiFR+B68rrwyuwSookAYNqtoltt12k844Izj+z39KS5ZUrUZfmPECACDt8jd4vfeeNHp0cPyCC6Tdd6/6eQcMcHs5lrZypXT33VU/pw8ELwAA0i4/g1dRkXTxxcHxRo3CN7+ujObNXZfkRA89JM2bl9q504ngBQBA2uVn8HrkEenrr4Pjd94ZTfuE/v2l2rXLjv3xR3i/r0xB8AIAIO3yL3gtWiTdeGNw/K9/lc46K5rnaNxYuuii4PjQodKPP0bzHFEqLpbmzg2Os7geAIBIRRK8jDFHGWNmGGNmGWOuieKcaXPdddKyZcHxf/0reEdiKq66SqpXr+xYUVHqlzLTYeFCaf36smP16gXrBwAAKUk5aRhjCiQNlnS0pNaSTjHGtE71vGnxxRfh2/j06SPtvXe0z9WwodSvX3B85Ejp22+jfa5UcZkRAAAvopji2VvSLGvtj9badZJGSzougvNGq7hYuvBCtydhaVts4dZ2pcOll7oAllhH2KXOOBG8AADwIorg1VRS6d/cc0rGyjDGnGuMmWSMmbRo0aIInraSnnpKmjgxOH7zzdLWW6fnOevVk669Njj+/PPS5Mnpec6qIHgBAOBFFMErrNOoDQxYO9Ra28Fa26FRo0YRPG0lLF8uXROy9KxNG9e3K53OP19q0iQ4fv316X3eyggLXiysBwAgclEErzmSSk+PNJOUWQ2rbrrJLSBP9NBDUmFhep+7dm3phhuC4+PHSx9+mN7nrqiw7YKY8QIAIHJRBK8vJO1sjNnRGFNDUk9JL0dw3mhMm+YCVqKTT3Z7MvrQp4/UokVwvH//4JqzOHCpEQAAL1IOXtbaIkkXSpogabqkZ62136R63khY6zrUb9hQdrxOHWngQH911KjhZt0Sffih9Oab/upIhuAFAIAXkTSusta+bq3dxVq7k7U2c9qzjx0rvftucPy66/wHi1NPlVqHdNmIe9Zrw4bwrYxY4wUAQORyt3P977+H99Haaafw8XQrKJBuvTU4Pnmy9OKL/uvZaP784IxggwbSZpvFUw8AADksd4PXXXeFX0J74AGpVi3/9UjS8cdL7dsHx6+/Phh+fGFhPQAA3uRm8PrhB+mee4LjnTtLXbv6r2cjY8I3yp4+Xfr3v/3XI7G+CwAAj3IzeF12mbRuXdmxGjXcbFfcjjhCOuig4PiAAcGafSB4AQDgTe4FrzfekF55JTjer5+0887+60mUbNbrp5/C95FMN4IXAADe5FbwWrtWuuSS4HjTpu5OxkxxwAHS0UcHx2+9VVqzxm8tdK0HAMCb3Ape998vff99cHzgQKluXf/1lOe224Jj8+dLgwf7rYPF9QAAeJM7wWvu3PAwc/DBUo8e/uvZlHbtpJNOCo7fdZe0YoW/OrjUCACAN7kTvK680vXuKq2gQBo0yK2rykS33CJVS3gJlixxM3c+rF/vZtkScakRAIC0yI3gNW+e9HLI9pDnny/tsYf/eiqqVSupV6/g+D//6QJYus2bF+ya36hRfH3OAADIcbkRvJo0cb2wSl9S3Gor6eab46upogYMkKpXLzu2cqV0993pf24W1gMA4FVuBC/JrUsaPVp65x2pTRvpzjvd1jeZbscdpb59g+P/+lf4HopRYn0XAABe5U7w2uhvf5O+/FLq0yfuSiru+uul2rXLjq1ZE36zQJS4oxEAAK9yL3hJUmFhcNF6JmvcWLrwwuD4kCHSxx+n73mZ8QIAwKssSic57uqrpc03LztWXCydcYZb85UOBC8AALwieGWKhg2l/v2D4z/+KF1+eXqek8X1AAB4RfDKJP36SfvsExwfNix8/8lUMeMFAIBXBK9MUlgojRgh1akTPHbOOdLChdE919q1wfMZ4/a1BAAAaUHwyjQtW4Z3rl+4UDr33GDD06qaOzc4ts02Uo0a0ZwfAAAEELwyUd++UpcuwfFx46Thw6N5Di4zAgDgHcErExnj1nU1bBg8dsklbsF9qgheAAB4R/DKVNtuKw0dGhxftUo680xpw4bUzs8djQAAeEfwymQnnCD17h0c/+gjaeDA1M5N13oAALwjeGW6Bx+UdtghOH7DDdLUqVU/L5caAQDwjuCV6erVk55+2q37Km39eun006U//qjaeQleAAB4R/DKBgcdJF1xRXD8m2/Cu91XBMELAADvCF7Z4tZbpd13D47ff7/07ruVO9fq1dKSJWXHqlVzm3UDAIC0IXhli5o1pZEjgw1OrXV3OS5fXvFzhS2sb9zYdc4HAABpQ/DKJnvsId12W3B89mzpoosqfh7uaAQAIBYEr2xz+eVuzVeiESOk556r2DlY3wUAQCwIXtmmoEB66ilp882Dx/7+d2n+/E2fg+AFAEAsCF7ZqHlz6aGHguNLl0p9+mx6I2261gMAEAuCV7Y64wzX2T7R+PHSo4+W/7XMeAEAEAuCV7YyRhoyRNpmm+Cxfv2kmTOTfy2L6wEAiAXBK5tttZX0+OPB8TVrXFf79evDv44ZLwAAYkHwynZdukjnnRcc/+IL6Y47guOrVknLlpUdKywMnzkDAACRInjlgoEDpZYtg+O33uoCWGlhs11Nmri7JQEAQFoRvHJB3bquj1e1hJdzwwZ3yXH16j/HuMwIAEBsCF65Yp99pOuuC47PnCldddWf/8/CegAAYkPwyiU33ii1bx8cHzzYtZmQmPECACBGBK9cUr26u+RYq1bwWJ8+0pIlBC8AAGJE8Mo1rVpJ99wTHJ8/X/rHP+haDwBAjAheueiCC6TDDguOP/ec9N57wXFmvAAA8ILglYuqVZOGD5fq1w8eW7cuOEbwAgDAC4JXrmrWTHrkkU0/rkYNqVGj9NcDAAAIXjmtZ0/plFPKf0yzZsH+XwAAIC34jZvrBg+WmjZNfpyF9QAAeEPwynUNGkhPPpn8OOu7AADwhuCVDw47TLr44vBjBC8AALwheOWLu+6SdtstON66tf9aAADIUwSvfFG7tjRqVNkWE40bSyedFF9NAADkmcK4C4BH7dpJkydLjz8uFRRIl1ziAhkAAPCC4JVvWrSQbr897ioAAMhLXGoEAADwhOAFAADgCcELAADAE4IXAACAJwQvAAAATwheAAAAnhC8AAAAPCF4AQAAeELwAgAA8ITgBQAA4ElKwcsYc68x5jtjzFfGmBeNMfU3/VUAAAD5KdUZr7cktbHW7iFppqRrUy8JAAAgN6UUvKy1b1pri0r+9zNJzVIvCQAAIDdFucarj6Q3kh00xpxrjJlkjJm0aNGiCJ8WAAAgOxhrbfkPMOY/krYNOdTfWjuu5DH9JXWQdILd1AkldejQwU6aNKkK5QIAAHhnojpR4aYeYK09rNxKjDlTUldJh1YkdAEAAOSrTQav8hhjjpJ0taSDrbWroykJAAAgN6W6xutfkjaX9JYxZqox5tEIagIAAMhJKc14WWtbRlUIAABArqNzPQAAgCcELwAAAE8IXgAAAJ4QvAAAADwheAEAAHhC8AIAAPCE4AUAAOAJwQsAAMATghcAAIAnBC8AAABPCF4AAACeELwAAAA8IXgBAAB4QvACAADwhOAFAADgCcELAADAE4IXAACAJwQvAAAATwheAAAAnhC8AAAAPCF4AQAAeELwAgAA8ITgBQAA4AnBCwAAwBOCFwAAgCcELwAAAE8IXgAAAJ4QvAAAADwheAEAAHhC8AIAAPCE4AUAAOAJwQsAAMATghcAAIAnBC8AAABPCF4AAACeELwAAAA8IXgBAAB4QvACAADwhOAFAADgCcELAADAE4IXAACAJwQvAAAATwheAAAAnhC8AAAAPCF4AQAAeELwAgAA8ITgBQAA4Imx1vp/UmMWSfqf9yfOLltJWhx3EXmO1yB+vAbx4zWIH69B/GpZa9tEcaLCKE5SWdbaRnE8bzYxxkyy1naIu458xmsQP16D+PEaxI/XIH7GmElRnYtLjQAAAJ4QvAAAADwheGWuoXEXAF6DDMBrED9eg/jxGsQvstcglsX1AAAA+YgZLwAAAE8IXgAAAJ4QvDwxxjxhjFlojJlWamxPY8ynxpivjTGvGGPqlYxXN8Y8VTI+3RhzbamvOcoYM8MYM8sYc00cf5ZsFdVrUHK8wBjzpTHmVd9/jmwW4ffBZcaYb4wx04wxzxhjasXx58lGlXwNahhjhpeM/9cY06lkvI4x5jVjzHclr8NdMf1xslIUr0GpY0ONMTNLXosTY/jjZCVjzHbGmHdLfrZ8Y4y5pGR8S2PMW8aY70v+26Bk3BhjBpX87v3KGNOu1LnOLHn898aYMzf55NZaPjx8SDpIUjtJ00qNfSHp4JLP+0i6teTzUyWNLvm8jqSfJTWXVCDpB0ktJNWQ9F9JreP+s2XLRxSvQamvu1zSvyW9GvefK5s+Ivo+aCrpJ0m1S449K6l33H+2bPmo5GtwgaThJZ9vLWmy3Bv2OpL+VjJeQ9KHko6O+8+WLR9RvAYl/3+zpNtKPq8maau4/2zZ8iGpsaR2JZ9vLmmmpNaS7pF0Tcn4NZLuLvm8s6Q3JBlJ+0j6vGR8S0k/lvy3QcnnDcp7bma8PLHWfiBpacLwrpI+KPn8LUkb361YSZsZYwol1Za0TtIKSXtLmmWt/dFau07SaEnHpbv2XBHRayBjTDNJXSQNS3fNuSaq10Cu+XPtkmN1JM1LZ925pJKvQWtJb5d83UJJyyR1sNautta+WzK+TtIUSc3SXHrOiOI1KDnWR9KdJceKrbV0t68ga+18a+2Uks9XSpou96buOElPlTzsKUndSj4/TtLT1vlMUn1jTGNJR0p6y1q71Fr7m9xrd1R5z03witc0SceWfN5d0nYln4+V9Luk+ZJ+kTTQWrtU7h/F7FJfP6dkDFVX2ddAkh6QdJWkYo915rJKvQbW2rmSBpaMzZe03Fr7pt+Sc06y1+C/ko4zxhQaY3aU1L7UMUmSMaa+pGNUEg5QZZV6DUr+3iXpVmPMFGPMc8aYbfyWnBuMMc0ltZX0uaRtrLXzJRfO5GYZpeS/fyv9e5ngFa8+ki4wxkyWm+pcVzK+t6QNkppI2lFSP2NMC7kpzkT0A0lNpV4DY0xXSQuttZNjqTY3VfY1aCD37nPHkmObGWNO9192Tkn2Gjwh94tkktwbjk8kFW38opIZx2ckDbLW/ui14txT2degUG6W8WNrbTtJn8q9IUElGGPqSnpe0qXW2hXlPTRkzJYznlQsezXCsdZ+J+kISTLG7CJ3+Upya1vGW2vXS1pojPlYbmp5tsq+22wmLrGkpAqvQVtJxxpjOkuqJameMWaktZZf/FVUhdfASvrJWruo5GtekLSfpJG+a88VyV4Da22RpMs2Ps4Y84mk70t96VBJ31trH/BXbW6qwmuwRNJqSS+WHHpO0tkeS856xpjqcqFrlLX2hZLhBcaYxtba+SWXEheWjM9R+O/fOZI6JYy/V97zMuMVI2PM1iX/rSbpekmPlhz6RdIhJXdRbCa3kO87ucWXOxtjdjTG1JDUU9LL/ivPHZV9Day111prm1lrm8v9/b9D6EpNFb4PfpG0T8mddUbSoXLrM1BFyV6Dkr/jzUo+P1xSkbX225L/v03SFpIujaXoHFPZ18C6ld2v6M9f+odK+tZ33dmq5GfH45KmW2vvK3XoZUkb70w8U9K4UuNnlPw82kduicN8SRMkHWGMaVAyG39EyVhSzHh5Yox5Ru4bZCtjzBxJAyTVNcZcUPKQFyQNL/l8cMnn0+SmMYdba78qOc+Fci9qgaQnrLXfePtDZLmoXgNUXYTfB2PlFnQXSfpSbKlSYZV8DbaWNMEYUyxprqReJedoJqm/XBCe4n6H6V/WWm44qYAoXoMSV0saYYx5QNIiSWd5KD9X7C/3d/m1MWZqydh1ku6S9Kwx5my5N3ndS469Lndn4yy5mcazJMlau9QYc6vcxIgk3VJqPXAotgwCAADwhEuNAAAAnhC8AAAAPCF4AQAAeELwAgAA8ITgBQAA4AnBCwAAwBOCFwAAgCf/BzaETR5O/e7hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (10,7))\n",
    "\n",
    "ax.plot(real_rates[\"real_rate\"].dropna(), color = \"red\",\n",
    "        linewidth = 5, label = \"Imports\")\n",
    "\n",
    "ax.axhline(y=avg_real,   \n",
    "           # Set the value equall to the average\n",
    "           color='k',       # make the color red\n",
    "           linestyle='--',   # the line style\n",
    "           linewidth=4)     # thickness of the line\n",
    "\n",
    "ax.spines[\"top\"].set_visible(False)\n",
    "ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "ax.set_xlim(dt.datetime(1977,1,1), dt.datetime(2000,1,1))\n",
    "#ax.set_ylim(5,15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a name=\"output\"></a> Save the Empirical Moments for Use as Input into Calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "outpath = \"src/calibration/data\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are the firm transition dyanics moments\n",
    "trans_moments.to_csv(outpath + \"/firm_moments.csv\", index = False, header = False)\n",
    "\n",
    "# Exit rate entry rate\n",
    "entry_moment = pd.DataFrame([0.02])\n",
    "entry_moment.to_csv(outpath + \"/entry_moment.csv\", index = False, header = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The bejk moments\n",
    "bejk.to_csv(outpath + \"/bejk_moments.csv\", index = False, header = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Trade moments\n",
    "trade_moment = pd.DataFrame([1-.01*avg_trade])\n",
    "trade_moment.to_csv(outpath + \"/trade_moments.csv\", index = False, header = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Growth and r moments\n",
    "\n",
    "growth_and_r_moments = pd.DataFrame([.01*avg_real, .01*mean_tfp])\n",
    "growth_and_r_moments.to_csv(outpath + \"/growth_and_r_moments.csv\", index = False, header = False)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}