{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "\n",
    "The unique feature of structural models compared to other econometric models is the ability to provide counterfactuals. There exist three different ways to simulate counterfactuals with **respy**. Ordered from no data to a panel data on individuals, there is:\n",
    "\n",
    "1. *n-step-ahead simulation with sampling*: The first observation of an individual is sampled from the distribution implied by the initial conditions, i.e., the distribution of observed variables or initial experiences, etc. in the first period. Then, the individuals are guided for $n$ periods by the decision rules from the solution of the model.\n",
    "\n",
    "2. *n-step-ahead simulation with data*: Take the first observation of each individual from the data and do as in 1..\n",
    "\n",
    "3. *one-step-ahead simulation*: Take the complete data and find for each observation the corresponding outcomes, e.g, choices and wages, using the decision rules from the model solution.\n",
    "\n",
    "You can find more information on initial conditions and how to express them with **respy** in the guide below."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<div class=\"d-flex flex-row gs-torefguide\">\n",
    "    <span class=\"badge badge-info\">How-to Guide</span>\n",
    "\n",
    "    Find out more about<a\n",
    "    href=\"how_to_initial_conditions.html\">Initial Conditions</a>.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, you will learn how to use all of the three methods and what their purpose is. We take the basic model from Keane and Wolpin (1997) and their original data.\n",
    "\n",
    "We further restrict the model to 12 periods to reduce the complexity and runtime. Thus, the following simulations will differ from the real data a lot because the parameters belong to a model with fifty periods which distorts choice patterns for shorter time frames. We start by importing the model and data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import respy as rp\n",
    "\n",
    "params, options, df = rp.get_example_model(\"kw_97_basic\")\n",
    "\n",
    "options[\"n_periods\"] = 11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. n-step-ahead simulation with sampling\n",
    "\n",
    "The first exercise is an *n-step-ahead simulation with sampling*. This method and the following *n-step-ahead simulation with data* are the cornerstone for providing predictions of behavior in counterfactural environments. Here, the first observation of each individual is sampled from the distribution implied by the initial conditions of the model. The initial conditions determine, for example, the distribution of observed variables, experiences, and previous choices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = rp.get_simulate_func(params, options, \"n_step_ahead_with_sampling\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_1 = simulate(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1e8856ad248>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5yOdf7H8dfHNIwUyWmJjEoHM46DMJRDIbRISdgckrK2tbW16LCd1qbdtl8/rdrVVrRboiQd6CSHiBiRUwdsYqKI2rDxY/r8/rgv0+AeM2Nm7pu53s/HYx73dX/u67o+32uYz/29v9d1fy9zd0REJBxKxbsBIiISOyr6IiIhoqIvIhIiKvoiIiGioi8iEiInxbsBealcubInJyfHuxkiIieUZcuWfePuVQ6PH/dFPzk5mYyMjHg3Q0TkhGJmX0SLa3hHRCREVPRFREJERV9EJESO+zF9Kbn2799PZmYme/fujXdTQicpKYmaNWuSmJgY76ZIjKnoS9xkZmZy6qmnkpycjJnFuzmh4e7s2LGDzMxM6tSpE+/mSIxpeEfiZu/evVSqVEkFP8bMjEqVKukTVkjlq+ib2UYzW2VmK8wsI4idbmZvm9m64LFijvVHm9l6M/vUzDrliKcF+1lvZuNMf+2hp/8C8aHfe3gVpKffzt0buXvT4PkoYLa71wVmB88xs3pAHyAF6Aw8ZmYJwTaPA0OBusFP58IfgoiI5Fdhhne6A5OC5UlAjxzx5919n7t/DqwHmptZdaC8uy/yyCT+z+TYRiQuvvrqK/r06cPZZ59NvXr16NKlCxMmTKBbt24F2s+QIUNYu3ZtMbVSpOjk90SuA2+ZmQN/d/cJQDV33wrg7lvNrGqw7hnA4hzbZgax/cHy4fEjmNlQIp8IOPPMM3Nt1Pgb381n86Mb/rf2x7xtYXKfiHnjmXvbF98f87ZVa5fP9TV3p2fPngwYMIDnn38egBUrVvDqq68WOO8f73r4iG2OljsvxXXMedHfVOzyxit3fnv66e7eBLgMGG5mFx1l3WiDhX6U+JFB9wnu3tTdm1apcsTUESJFYs6cOSQmJnLjjTdmxxo1akSbNm3YvXs31w37BentmzJsxBAO3mFu/sK5dOjSmos7tWTEbcPZt28fAD2v7sqKlR8C8O7cd7ikaxsaNmxIhw4dANizZw+DBw+mWbNmNG7cmBkzZsT4aEUi8lX03X1L8LgNmA40B74OhmwIHrcFq2cCtXJsXhPYEsRrRomLxMXq1atJS0uL+try5cu5//djee+dJXyxaSNLMhazd+9eRtz6Syb8dSLz3lxEVtYBJv7ryUO2+2bHN/x29K958m//5KOPPuKFF14AYMyYMbRv356lS5cyZ84cbrvtNvbs2VPsxyhyuDyLvpmVM7NTDy4DHYHVwCvAgGC1AcDBrssrQB8zK2NmdYicsF0SDAXtMrMWwVU71+bYRuS40rx5c2pUP4NSpUqRWq8+mzM3seHf6zizVm3OPuscAHr3uobFSxYest2y5Utp0bwVtWslA3D66acD8NZbbzF27FgaNWpE27Zt2bt3L5s2bYrpMYlA/sb0qwHTg0u8TgKec/c3zGwpMNXMrgM2AVcBuPsaM5sKrAUOAMPdPSvY1zBgIlAWmBX8iMRFSkoKL774YtTXypQpk72ckJDAgQMHsod4jsbdo14O6e5MmzaN884779gbLFIE8uzpu/u/3b1h8JPi7mOC+A537+DudYPHnTm2GePuZ7v7ee4+K0c8w91Tg9d+5fn5KxIpJu3bt2ffvn088cQT2bGlS5cyb968qOufc/a5bM7cxOcbNwDw4ktTaHlh60PWadqkOYs+WMgXmzcCsHNn5M+iU6dOPProo9lvHMuXLy/qwxHJF30jV0LLzJg+fTpvv/02Z599NikpKdxzzz3UqFEj6vpJSUk88ufxDPnlAC7u1BIrVYoB/QYfsk7lSpV56I//y+Ab+tOwYUOuvvpqAO666y72799PgwYNSE1N5a677ir24xOJRnPvSKjVqFGDqVOnHhG//vrrsy+bfOC+h7LjF6W3ZfbMBUesP33K69nLHdpdSod2lx5y6WTZsmX5+9//XpRNFzkm6umLiISIir6ISIio6IuIhIiKvohIiKjoi4iEiIq+iEiI6JJNOW4kj3o975UKYMmwNkW6P5GSQD19CbWNGzeSmpp6RLxt27bZs2bGQ9P0+uzYuQOAOvWif1lM5Fio6IuUIO7Ojz/+GO9myHFMRV9C78CBAwwYMIAGDRpw5ZVX8t///veQ13P2tF+d+TK//u0wIDKN8uAb+9Pp523p9PO2LMlYTG52797NoEGDqF+/Pg0aNGDatGkATJ48mfr165OamsrIkSOP2s49e3bTq+/lXNK1DRd3asmstyLDYZs2f0HrDs0YeectXNK1DZs3bz6m34OEg8b0JfQ+/fRTnnzySdLT0xk8eDCPPfZYvra7896R3HDdcC5s1pLMLzfT59orWDB7adR177//fipUqMCqVasA+Pbbb9myZQsjR45k2bJlVKxYkY4dO/Lyyy/To0f0u4iWKZPExL8/y6mnlmfHzh106dmBzpd2AWD9v9fxyJ8f48E/PFyoO2dJyaeiL6FXq1Yt0tPTAejfvz/jxo3L13bzF87ls3WfZj/ftXsXu3fv4pRTTj1i3XfeeSf7lowAFStWZP78+bRt25aDd4fr168f8+fPz7Xouzt//PN9LFryPqWsFF99tZXt2yP3Lqp5Ri2aNmmWvwOWUFPRl9A7fP77oz0/eHtEAP/xR16f/jZlk8rmmSPaPPsFnVl82stT+WbHN7z96jwSExNpml6fvfv2AnDyyeUKtC8JLxV9OW5sHNs1arwwNwnPj02bNrFo0SJatmzJ5MmTad26dfbN0QGqVK7CZ+s/5Zyz6jLzzVc5pVykJ39xm/Y8NWkCw28YAcDqNStJTWkQNUfHjh3561//yiOPPAJEhncuvPBCRowYwTfffEPFihWZPHkyN910U67t/H7X91SuXIXExEQWvD+fzV/qzltScDqRK6F3wQUXMGnSJBo0aMDOnTsZNmzYIa/fOfIe+g/uTa9rLqda1Z9lx8fc8ydWrFxO286taHNJcyY9+1SuOe68806+/fZbUlNTadiwIXPmzKF69eo88MADtGvXjoYNG9KkSRO6d++e6z569ejNRyuX0/Hyi5k2Yyp1zz638AcvoaOevoRacnIya9euPSI+d+7c7E8Yl3fpweVdjhxnr3R6JZ4YPzFfeU455RQmTZp0RLxv37707dv3iHjGwlXZy5+v3ZKdb+b0d6Luf/5buV85JJKTevoiIiGinr5IEZo89V888fTfADipdKRPlZ6ezvjx4+PZLJFsKvoiReia3v25pnd/AF0vL8clDe+IiISIir6ISIio6IuIhIjG9OX4cU+FqOGqx7i7bYPynnhs48aNdOvWjdWrVx9jFpETi3r6IiIhoqIvoZeVlcX1119PSkoKHTt25IcffmDFihVc1qMDbTu3YuDQfnz3n28B6Hl1V+66bzTde19G6w7NWP7RMgbd0I8WbRvzwEP3Z+/zxelTaN68OY0aNeKGG24gKysrXocncggVfQm9devWMXz4cNasWcNpp53GtGnTuPbaa7lr1L3MfeN9Lji/Hg898mD2+qUTE5kxdRbX9hvMgOv7Mvb+vzDvzcVMefE5dn67k8/Wf8rLr73EwoULWbFiBQkJCTz77LNxPEKRn2hMX0KvTp06NGrUCIC0tDQ2bNjAd999R6sWrQG4utc1DPnlgOz1OwVz2F9wfj3OO/f87Pl4ap+ZzJYtmXyQsZiVq1bQrFlkquMffviBqlWP9cyESNFS0ZfQK1OmTPZyQkIC33333VHXL126NAClrBSlS/+0bSkrxYGsA7g7vXtdw7jHHi6eBosUgoZ3RA5ToUIFKlasyOIl7wPwwktTaHlh63xv3yb9Yl6bNYNt2yI3ONm5cydffPFFsbRVpKDU05fjxz3/iRou7vn0o5k0aRLXDbqeH/b+QO1ayfzvQ/mfO+e8uucz6rd30rFjR3788UcSExMZP348tWvXLsYWi+RPvou+mSUAGcCX7t7NzE4HpgDJwEagt7t/G6w7GrgOyAJ+7e5vBvE0YCJQFpgJjPCC3j5IpAglJycfco3+rbfemr086+XZR6w/fcrr2cvpLduQ3rJN1Nd6XN6Lob8aVNTNFSm0ggzvjAA+zvF8FDDb3esCs4PnmFk9oA+QAnQGHgveMAAeB4YCdYOfzoVqvYiIFEi+ir6Z1QS6Av/IEe4OHLwrxCSgR4748+6+z90/B9YDzc2sOlDe3RcFvftncmwjIiIxkN+e/iPA74Afc8SquftWgODx4DVpZwA5v/+eGcTOCJYPjx/BzIaaWYaZZWzfvj2fTRQRkbzkWfTNrBuwzd2X5XOfFiXmR4kfGXSf4O5N3b1plSpV8plWRETykp8TuenAz82sC5AElDezfwFfm1l1d98aDN1sC9bPBGrl2L4msCWI14wSFxGRGMmzp+/uo929prsnEzlB+6679wdeAQ5+TXEAMCNYfgXoY2ZlzKwOkRO2S4IhoF1m1sLMDLg2xzYiIhIDhblOfyww1cyuAzYBVwG4+xozmwqsBQ4Aw9394GxTw/jpks1ZwY8IAPUn1S/S/c1uu7BI9vPG2zP5bN0n/PqXt/Dn/3mAcuXK8cuhv+bBh8fQonkrLm7djr8/+Ri/6DuQk8ueXCQ5RYpLgYq+u88F5gbLO4AOuaw3BhgTJZ4BpBa0kSLx1PnSLnQO5tvJaeQtd2QvP/HU41zZ8+oCFf2srCwSEhLyXlGkCGkaBgm1jRs3cv755zNkyBBSU1Pp168f77zzDunp6bRo25gPVyzj+ReeZfTvbz1i21//dhivznyZJ57+G19t28oV13SjZ59uAPzujptp2rQpKSkp3H333dnbJCcnc99999G6dWvGjh1LkyZNsl9bt24daWlpxX/QEmoq+hJ669evZ8SIEaxcuZJPPvmE5557jgULFnD37X/gf8f/Jc/trx90Iz+rWp2XJr/G9OdfA2D0bXeRkZHBypUrmTdvHitXrsxePykpiQULFnDHHXdQoUIFVqxYAcDTTz/NwIEDi+UYRQ5S0ZfQq1OnDvXr16dUqVKkpKTQoUMHzIwLzq/H5sxNx7TPV16bTpMmTWjcuDFr1qxh7dq12a9dffXV2ctDhgzh6aefJisriylTptC3b99CH4/I0ajoS+jlnFq5VKlS2c8PTpVcUF9s3shjTzzK7NmzWblyJV27dmXv3r3Zr5crVy57uVevXsyaNYvXXnuNtLQ0KlWqVIgjEcmbir5IESh3yins3rMLgN27dnFy2XJUqFCBr7/+mlmzcr9ILSkpiU6dOjFs2DAGDdIEbVL8NLWyHDdWDVgVNR6PqZUL6hfXDKTvgCupWvVnTH/+NeqnNCAlJYWzzjqL9PT0o27br18/XnrpJTp27Bij1kqYqehLqB0+tfLEiROzl8+sVZv5by0GoM9V/QC47ebR2a+P+8vj2ctDBt7AkIE3HPJa1drlj8i3cePGI2ILFixg8ODBunxTYkJFXySOevbsyYYNG3j33Xfj3RQJCRV9kTiaPn16vJsgIaMTuSIiIaKiLyISIir6IiIhoqIvIhIiOpErx42Pz7+gSPdX6c0Pimxfz7/wLB+tWs4D9z1U6H0lJyeTkZFB5cqVi6BlIgWjnr6ISIio6Euo7dmzh65du9KwYUNSU1OZMmUKS5cupVWrVrTrnE6n7u3YvTsyvcJXX2+lz7VX0KJtY+574K7sfbw040Uu7tSSizq24P4Hfp8dnzx5MvXr1yc1NZWRI0fG/NhEotHwjoTaG2+8QY0aNXj99dcB+M9//kPjxo2ZMmUKtauex65d35OUVBaA1WtXMfv19yhdugzpHZpy3YAbSEhI4A8P3s1br87jtAqn0fsXPZj55ms0aZTGyJEjWbZsGRUrVqRjx468/PLL9OjRI56HK3JiF/32c4cXcg8fF0k75MRVv359br31VkaOHEm3bt047bTTqF69Os2aNWPbF99z6qk/TaXQJv1iypevAMC555xH5peb2fntTlq1aE3lSpHx+V49erN4yfuYGW3btqVKlSpAZH6d+fPnq+hL3Gl4R0Lt3HPPZdmyZdSvX5/Ro0czffp0zCzqumVK55iCOSGBAwcO4O5R180tLhJvKvoSalu2bOHkk0+mf//+3HrrrSxevJgtW7awdOlSAHbv3sWBA7nPqZ/WuCmLPljIjp07yMrKYvorL9LywnTSGjdl3rx5fPPNN2RlZTF58mQuvvjiWB2WSK5O6OEdKVku+CT6cFtxTq28atUqbrvtNkqVKkViYiKPP/447s5NN93Erv/sISkpiRefnZHr9tWq/ozbf3c3V1zTDXfnknaXclnHrgA88MADtGvXDnenS5cudO/evdiOQyS/VPQl1Dp16kSnTp2OiC9evPiQN5s+V/XLnl4Z4NmnpmYv9+p+Fb26X3XEPvr27Rv19ofRplcWiRUN74iIhIiKvohIiKjoi4iEiIq+iEiIqOiLiISIir6ISIjokk05boy/sWhvDn7V6KZFuj+RkkA9fZGjqFOvRtT4pH89ydRpk4HIXPtffb212NuyafMXXNSxBQALF71Hv8G9iz2nlDzq6YscgwH9r8tenvLic5x/3gX8rFr1OLbo6A4cOMBJJ+nPXdTTl5D705/+xLhx4wC4+eabad++PQCzZ8/ml7+5HoA//vk+2nVO57IeHdi2fRsAf/6fB3hswjhenfkyK1Yt55e/uZ72l7Xmh70/8NGq5fTo3YW0tDQ6derE1q25fwpYv349l1xyCQ0bNqRJkyZs2LABd+feP97JRR1bcHGnlrz86rSjHsOHK5bR9YpL6dClNV2vuJRPP/0UgIkTJ3LVVVdx+eWX07Fjx0L/rqRkyLPom1mSmS0xs4/MbI2Z3RvETzezt81sXfBYMcc2o81svZl9amadcsTTzGxV8No4y206Q5EYueiii3jvvfcAyMjIYPfu3ezfv58FCxbQollL/vvfPaQ1bsacNxbS8sJW/Ov5SYdsf3mXHjSq35jHHnmCd2ct4KSEk7j97t/xj8efYdmyZQwePJg77rgj1/z9+vVj+PDhfPTRR7z//vtUr16dl156idVrVzFn1kJe+NcM7nvg93y97atc91H37LrMmDqL2TMX8Ltbbuf222/Pfm3RokVMmjSJd98t2vMlcuLKz+e9fUB7d99tZonAAjObBVwBzHb3sWY2ChgFjDSzekAfIAWoAbxjZue6exbwODAUWAzMBDoDs4r8qETyKS0tjWXLlrFr1y7KlClDkyZNyMjI4L333uPuUX+kdOnSdOzQGYAGqY2Yt2DOUfe3/t/r+OSzj+ndvwcnlS5FVlYW1atHH/bZtWsXX375JT179gQgKSkJgAULFtDz51eSkJBA1SpVaXlhOss/+pB656dE3c/3u77npt8O498bN2BmOFnZr1166aWcfvrpBf69SMmVZ9H3yMTgu4OnicGPA92BtkF8EjAXGBnEn3f3fcDnZrYeaG5mG4Hy7r4IwMyeAXqgoi9xlJiYSHJyMk8//TStWrWiQYMGzJkzhw0bNnDuOedx0kmJ2fPrJyQkkHUg66j7c3fOq3s+M6e/Q9Xa5fNctyDx3Dz4lzGkt2zDxAnPsmnzF1zZ//Ls18qVK1egfUnJl68zO2aWACwDzgHGu/sHZlbN3bcCuPtWM6sarH4GkZ78QZlBbH+wfHhcBIDhf2sfNV6cUytDZIjnoYce4qmnnqJ+/frccsstpKWl5XozlcOVO+UUdu+O9IvOOasuO3Z+w9JlS+ha+xL279/PZ599RkrKkb308uXLU7NmzezbKO7bt4+srCwuuugiHn1kPFf36su3333L4iXvc/fo+9m7b2/U/N/v+p7qP4t8mpjy4nPH+FuQsMjXiVx3z3L3RkBNIr321KOsHu0vxY8SP3IHZkPNLMPMMrZv356fJoocszZt2rB161ZatmxJtWrVSEpKok2bNvnevs+VffndnTfT/rLWZP2YxZOPPcMfHrybhg0b0qhRI95///1ct/3nP//JuHHjaNCgAa1ateKrr76iZ8+e1LsghXaXpdOr7+XcNepeqlatlus+ht8wgjF/upduvTqS9ePRP4mIFOgaLnf/zszmEhmL/9rMqge9/OrAtmC1TKBWjs1qAluCeM0o8Wh5JgATAJo2bar7zkmx6tChA/v3789+/tlnnwGRTxifr/3pv+jlXXpweZfIPW5vu3l0drzbZd3pdtlPN0hJTWnAjKmz8hzeAahbt27Uk6x33/4H7r79D4fEzqxVm/lvRT5Ep7dsQ3rLyBtTs7TmLJrzYfZ6D4/7EwADBw5k4MCBebZBwiXPom9mVYD9QcEvC1wCPAi8AgwAxgaPB28v9ArwnJk9TOREbl1gibtnmdkuM2sBfABcCzxa1AcUK4W7KbtuyC4i8ZGfnn51YFIwrl8KmOrur5nZImCqmV0HbAKuAnD3NWY2FVgLHACGB1fuAAwDJgJliZzA1UlcCYXhw4ezcOHCQ2IjRoxg0KBBcWqRhFV+rt5ZCTSOEt8BdMhlmzHAmCjxDOBo5wMkZNw93ydMT2Tjx4+PdxMOUdArhKTk0DdyJW6SkpLYsWOHClCMuTs7duzI/l6AhIsm45C4qVmzJpmZmeR1hdauHdEvVcyPHf899sJWmLzxzJ2fvElJSdSsWTPP9aTkUdGXuElMTKROnTp5rleYKZdzu/Y/Pwo71XO8chcmr5R8Gt4REQkRFX0RkRBR0RcRCREVfRGREFHRFxEJERV9EZEQUdEXEQkRFX0RkRBR0RcRCRF9I/cEoymdRaQw1NMXEQkRFX0RkRBR0RcRCREVfRGREFHRFxEJERV9EZEQUdEXEQkRFX0RkRBR0RcRCREVfRGREFHRFxEJERV9EZEQUdEXEQkRzbIp+aYZPkVOfOrpi4iEiIq+iEiIqOiLiISIir6ISIio6IuIhIiKvohIiORZ9M2slpnNMbOPzWyNmY0I4qeb2dtmti54rJhjm9Fmtt7MPjWzTjniaWa2KnhtnJlZ8RyWiIhEk5+e/gHgt+5+AdACGG5m9YBRwGx3rwvMDp4TvNYHSAE6A4+ZWUKwr8eBoUDd4KdzER6LiIjkIc+i7+5b3f3DYHkXkW/ZnAF0ByYFq00CegTL3YHn3X2fu38OrAeam1l1oLy7L3J3B57JsY2IiMRAgcb0zSwZaAx8AFRz960QeWMAqgarnQFszrFZZhA7I1g+PB4tz1AzyzCzjO3btxekiSIichT5LvpmdgowDfiNu39/tFWjxPwo8SOD7hPcvam7N61SpUp+mygiInnIV9E3s0QiBf9Zd38pCH8dDNkQPG4L4plArRyb1wS2BPGaUeIiIhIjeU64Flxh8yTwsbs/nOOlV4ABwNjgcUaO+HNm9jBQg8gJ2yXunmVmu8ysBZHhoWuBR4vsSKTE0kRvIkUnP7NspgO/AFaZ2YogdjuRYj/VzK4DNgFXAbj7GjObCqwlcuXPcHfPCrYbBkwEygKzgh8REYmRPIu+uy8g+ng8QIdcthkDjIkSzwBSC9JAEREpOvpGrohIiKjoi4iEiIq+iEiIqOiLiISIir6ISIio6IuIhIiKvohIiKjoi4iEiIq+iEiI5GcaBpFQKtycP6B5f+R4pJ6+iEiIqOiLiISIir6ISIio6IuIhIiKvohIiKjoi4iEiIq+iEiIqOiLiISIir6ISIio6IuIhIiKvohIiGjuHZHjUOHm/dGcP5I79fRFREJERV9EJERU9EVEQkRFX0QkRFT0RURCRFfviEg23S2s5FNPX0QkRE7onn7v0YVr/qoiaoeIyIlCPX0RkRA5oXv6qz7fFO8miIicUPIs+mb2FNAN2ObuqUHsdGAKkAxsBHq7+7fBa6OB64As4Nfu/mYQTwMmAmWBmcAId/eiPZzYKczQkoaVRCRe8lO5JgJ/BZ7JERsFzHb3sWY2Kng+0szqAX2AFKAG8I6ZnevuWcDjwFBgMZGi3xmYVVQHEhZ6sxGRwsizgrj7fDNLPizcHWgbLE8C5gIjg/jz7r4P+NzM1gPNzWwjUN7dFwGY2TNAD07goq+hJZGipUnmYuNYu43V3H0rgLtvNbOqQfwMIj35gzKD2P5g+fB4VGY2lMinAs4888xjbGLJpDcbESmMoj6Ra1FifpR4VO4+AZgA0LRp0xN23L+k0dCSyInvWP+Kvzaz6kEvvzqwLYhnArVyrFcT2BLEa0aJi+RJbzYiRedY/5peAQYAY4PHGTniz5nZw0RO5NYFlrh7lpntMrMWwAfAtcCjhWq5xFy8hpY0pCUlVTzOY+Tnks3JRE7aVjazTOBuIsV+qpldB2wCrgJw9zVmNhVYCxwAhgdX7gAM46dLNmdxAp/EFSmp9C33ki8/V+9ck8tLHXJZfwwwJko8A0gtUOtE4iiMBVCfqkq+E/obuSLFKZ4FMIznMeJ1zPG8VDQex6yiLyLHhTB+yojHMavoixyHwlgAJTY0y6aISIiopy8ioRa28ycq+iISamEbStPwjohIiKjoi4iEiIq+iEiIqOiLiISIir6ISIjo6h0RkThJ3vvcMW+78Ri3U9EXkWyFKUJw7IVIYkdFX0SOC/Ho9YaRir7IcUgFUIqLir5ILjTUEQ5he4NV0Zd8i9cfR9j+KEWKk4r+CUYFUEQKQ9fpi4iEiIq+iEiIaHjnGGmYRUROROrpi4iEyAnd09cldSIiBaOevohIiKjoi4iEiIq+iEiIqOiLiISIir6ISIio6IuIhIiKvohIiKjoi4iEiIq+iEiIqOiLiIRIzIu+mXU2s0/NbL2ZjYp1fhGRMItp0TezBGA8cBlQD7jGzOrFsg0iImEW655+c2C9u//b3f8PeB7oHuM2iIiElrl77JKZXQl0dvchwfNfABe6+68OW28oMDR4eh7w6TGmrAx8c4zbFka88sYzt445HLnDljeeuQubt7a7Vzk8GOuplS1K7Ih3HXefAEwodDKzDHdvWtj9nCh545lbxxyO3GHLG8/cxZU31sM7mUCtHM9rAlti3AYRkdCKddFfCtQ1szpmVhroA7wS4+T9KmcAAAXtSURBVDaIiIRWTId33P2Amf0KeBNIAJ5y9zXFmLLQQ0QnWN545tYxhyN32PLGM3ex5I3piVwREYkvfSNXRCREVPRFREKkRBb9eE31YGZPmdk2M1sdq5xB3lpmNsfMPjazNWY2Ioa5k8xsiZl9FOS+N1a5g/wJZrbczF6Lcd6NZrbKzFaYWUYM855mZi+a2SfBv3fLGOU9LzjWgz/fm9lvYpT75uD/1mozm2xmSTHKOyLIuaa4jzVa7TCz083sbTNbFzxWLJJk7l6ifoicIN4AnAWUBj4C6sUo90VAE2B1jI+5OtAkWD4V+CyGx2zAKcFyIvAB0CKGx34L8BzwWox/5xuByrHMGeSdBAwJlksDp8WhDQnAV0S+/FPcuc4APgfKBs+nAgNjkDcVWA2cTOSCl3eAusWY74jaAfwJGBUsjwIeLIpcJbGnH7epHtx9PrAzFrkOy7vV3T8MlncBHxP5Y4lFbnf33cHTxOAnJlcHmFlNoCvwj1jkizczK0+kODwJ4O7/5+7fxaEpHYAN7v5FjPKdBJQ1s5OIFOFYfLfnAmCxu//X3Q8A84CexZUsl9rRncibPMFjj6LIVRKL/hnA5hzPM4lRATwemFky0JhIjztWORPMbAWwDXjb3WOV+xHgd8CPMcqXkwNvmdmyYNqQWDgL2A48HQxp/cPMysUod059gMmxSOTuXwIPAZuArcB/3P2tGKReDVxkZpXM7GSgC4d+sTQWqrn7Voh07ICqRbHTklj08zXVQ0lkZqcA04DfuPv3scrr7lnu3ojIN6ybm1lqcec0s27ANndfVty5cpHu7k2IzBg73MwuikHOk4gMATzu7o2BPUQ+9sdM8KXKnwMvxChfRSI93jpADaCcmfUv7rzu/jHwIPA28AaRYeIDxZ03Fkpi0Q/lVA9mlkik4D/r7i/Fow3BUMNcoHMM0qUDPzezjUSG8Nqb2b9ikBcAd98SPG4DphMZVixumUBmjk9SLxJ5E4ily4AP3f3rGOW7BPjc3be7+37gJaBVLBK7+5Pu3sTdLyIy9LIuFnlz+NrMqgMEj9uKYqclseiHbqoHMzMi47wfu/vDMc5dxcxOC5bLEvkj/aS487r7aHev6e7JRP6N33X3Yu8BAphZOTM79eAy0JHIcECxcvevgM1mdl4Q6gCsLe68h7mGGA3tBDYBLczs5OD/eQci56yKnZlVDR7PBK4gtscNkbo1IFgeAMwoip3GepbNYuexn+ohm5lNBtoClc0sE7jb3Z+MQep04BfAqmBsHeB2d58Zg9zVgUnBDXJKAVPdPaaXT8ZBNWB6pAZxEvCcu78Ro9w3Ac8GHZp/A4NilJdgbPtS4IZY5XT3D8zsReBDIsMry4ndtAjTzKwSsB8Y7u7fFleiaLUDGAtMNbPriLz5XVUkuYLLgUREJARK4vCOiIjkQkVfRCREVPRFREJERV9EJERU9EVEQkRFX0LNzLKCWSNXm9kLwWWJBdn+H2ZWrwDrDzSzvxa8pSJFQ0Vfwu4Hd2/k7qnA/wE35ndDM0tw9yHuHusvSIkcMxV9kZ+8B5wDYGb9g/sErDCzvwdfPsPMdpvZfWb2AdDSzOaaWdPgtWuCOfZXm9mDB3dqZoPM7DMzm0fki3QicaOiLwIE0/ZeRuRbzRcAVxOZVK0RkAX0C1YtR2TO8wvdfUGO7WsQmaCrPdAIaGZmPYI5U+4lUuwvBfI9FCRSHErcNAwiBVQ2x9QV7xGZw2gokAYsDaZaKMtPk11lEZnY7nDNgLnuvh3AzJ4lMvc9h8WnAOcWw3GI5IuKvoTdD0FvPlswsdckdx8dZf297p4VJR5tSu+DNNeJHDc0vCNypNnAlTlmWTzdzGrnsc0HwMVmVjkY/7+GyN2WPgDaBjfjSKSIJs0SOVbq6Yscxt3XmtmdRO6MVYpglkUg19sDuvtWMxsNzCHS65/p7jMAzOweYBGROz99SGT2V5G40CybIiIhouEdEZEQUdEXEQkRFX0RkRBR0RcRCREVfRGREFHRFxEJERV9EZEQ+X/AMfjNlKRSmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_1.groupby(\"Period\").Choice.value_counts().unstack().plot.bar(stacked=True, rot=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. n-step-ahead simulation with data\n",
    "\n",
    "For the second variant of the n-step-ahead simulation, you do not need to specify, e.g., the distribution of initial experiences. Instead, the first observations in the data are taken as the starting points for the simulation. Thus, restrict the data to first period and pass it to the constructor for the simulation function.\n",
    "\n",
    "This routine is especially useful if the policy simulation results should be very similar to sample or the distribution of characteristics in the first period is complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = rp.get_simulate_func(params, options, \"n_step_ahead_with_data\", df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2 = simulate(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1e8848127c8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5yOdf7H8dfHNIxUkkMRGbU6mHE+hKFEIRRWyWkjSbV+ZWsr6bDaWpt2235tu7SrLfRbOUTSSSc5pIgROZWwiQkRHbCxTJ/fH/dldnAPM3PP3Deu9/PxmMd93d/re12f7zXG5/7e3+t7f29zd0REJBxKJLoBIiISP0r6IiIhoqQvIhIiSvoiIiGipC8iEiInJboBR1OhQgVPTU1NdDNERI4rixcv/sbdKx5afswn/dTUVDIzMxPdDBGR44qZfRmtXMM7IiIhctSkb2bPmdlWM1sRZd9dZuZmViFX2VAzW2tmq82sXa7yhma2PNj3lJlZ0V2GiIjkR356+mOB9ocWmlk14ApgQ66yWkAPIC04ZpSZJQW7nwYGAjWDn8POKSIixeuoY/ruPtfMUqPs+l/gHmB6rrLOwER33wt8YWZrgSZmth44zd3nA5jZ80AXYEZMrZfj2r59+8jKymLPnj2JbkropKSkULVqVZKTkxPdFImzQt3INbOrga/c/ZNDRmnOBhbkep4VlO0Ltg8tlxDLysri1FNPJTU1FY32xY+7s337drKysqhRo0aimyNxVuAbuWZ2MnA/8Jtou6OU+RHK84ox0MwyzSxz27ZtBW2iHCf27NlD+fLllfDjzMwoX7683mGFVGFm75wH1AA+CYZtqgIfm9lZRHrw1XLVrQpsCsqrRimPyt1Hu3sjd29UseJh00zlBKKEnxj6vYdXgZO+uy9390runuruqUQSegN33wK8AvQws1JmVoPIDduF7r4Z2GlmTYNZO9dz8L0AERGJg/xM2ZwAzAcuMLMsM7sxr7ruvhKYDKwC3gQGuXt2sPtW4B/AWmAduokrx4AtW7bQo0cPzjvvPGrVqkWHDh0YPXo0nTp1KtB5BgwYwKpVq4qplSJFJz+zd3oeZX/qIc+HA8Oj1MsE0gvYviMaect7MR0/6G+tExL7eIybyNhbv/yh0MdWqn5anvvcna5du9K3b18mTpwIwNKlS3n11VcLHPf3Dz5x2DFHin00xXXNR6P/U/GLm6jY+kSuhNasWbNITk7mlltuySmrV68eLVu2ZNeuXdx46y/IaN2IWwcP4MA3zM39YDZtOrTg0nbNGHz3IPbu3QtA1+s6snTZxwC8N/tdLu/Ykrp169KmTRsAdu/eTf/+/WncuDH169dn+nSNbkpiKOlLaK1YsYKGDRtG3bdkyRIe+c0I3n93IV9uWM/CzAXs2bOHwXf9ktF/Hcuct+aTnb2fsf989qDjvtn+Db8eejvP/u3/+OSTT3jxxRcBGD58OK1bt2bRokXMmjWLu+++m927dxf7NYocSklfJIomTZpQpfLZlChRgvRatdmYtYF1/1rDOdWqc965PwOge7eeLFj4wUHHLV6yiKZNmlO9WioAZ5xxBgBvv/02I0aMoF69erRq1Yo9e/awYcMGROLtmF9lU6S4pKWlMWXKlKj7SpUqlbOdlJTE/v37c4Z4jsTdo06HdHemTp3KBRdcUPgGixQB9fQltFq3bs3evXt55plncsoWLVrEnDlzotb/2XnnszFrA1+sXwfAlJcm0eziFgfVadSgCfM/+oAvN64HYMeOHQC0a9eOv/zlLzkvHEuWLCnqyxHJFyV9CS0zY9q0abzzzjucd955pKWl8dBDD1GlSpWo9VNSUnjyjyMZ8Mu+XNquGVaiBH179z+oToXyFXj893+m/819qFu3Ltdddx0ADz74IPv27aNOnTqkp6fz4IMPFvv1iUSj4R0JtSpVqjB58uTDym+66aacaZOPPvx4TvklGa2Y+ca8w+pPm/R6znaby66gzWVXHDR1snTp0vz9738vyqaLFIp6+iIiIaKkLyISIkr6IiIhoqQvIhIiSvoiIiGipC8iEiKasinHjNR7Xz96pQJYeGvLIj2fyIlAPX0JtfXr15OefviK361atcpZNTMRGmXUZvuO7QDUqBX9w2IihaGkL3ICcXd++umnRDdDjmFK+hJ6+/fvp2/fvtSpU4drrrmGf//73wftz93TfvWNl7n917cCkWWU+9/Sh3ZXt6Ld1a1YmLkgzxi7du3ihhtuoHbt2tSpU4epU6cCMGHCBGrXrk16ejpDhgw5Yjt3795Ft15XcXnHllzarhkz3o4Mh23Y+CUt2jRmyAN3cnnHlmzcuLFQvwcJB43pS+itXr2aZ599loyMDPr378+oUaPyddwDvx3CzTcO4uLGzcj6aiM9rv8582Yuilr3kUceoWzZsixfvhyAb7/9lk2bNjFkyBAWL15MuXLlaNu2LS+//DJdunSJeo5SpVIY+/fxnHrqaWzfsZ0OXdvQ/ooOAKz91xqe/OMoHvvdEzF9c5ac+JT0JfSqVatGRkYGAH369OGpp57K13FzP5jN52tW5zzfuWsnu3bt5JRTTj2s7rvvvpvzlYwA5cqVY+7cubRq1YqKFSsC0Lt3b+bOnZtn0nd3fv/Hh5m/8ENKWAm2bNnMtm1bAah6djUaNWicvwuWUFPSl9A7dP37Iz0/8PWIAP7TT7w+7R1Kp5Q+aoxo6+znZ33+3Ka+PJlvtn/DO6/OITk5mUYZtdmzdw8AJ59cpkDnkvBS0pdjxvoRHaOWx/Il4fmxYcMG5s+fT7NmzZgwYQItWrTI+XJ0gIoVKvL52tX87NyavPHWq5xSJtKTv7Rla54bN5pBNw8GYMXKZaSn1Ykao23btvz1r3/lySefBCLDOxdffDGDBw/mm2++oVy5ckyYMIHbbrstz3b+sPMHKlSoSHJyMvM+nMvGr/TNW1JwupEroXfRRRcxbtw46tSpw44dO7j11lsP2v/AkIfo07873XpexZmVzsopH/7QH1i6bAmt2jen5eVNGDf+uTxjPPDAA3z77bekp6dTt25dZs2aReXKlXn00Ue57LLLqFu3Lg0aNKBz5855nqNbl+58smwJba+6lKnTJ1PzvPNjv3gJnaP29M3sOaATsNXd04OyPwJXAf8B1gE3uPt3wb6hwI1ANnC7u78VlDcExgKlgTeAwV7Q97ciRSw1NZVVq1YdVj579uycdxhXdejCVR0OH2cvf0Z5nhk5Nl9xTjnlFMaNG3dYea9evejVq9dh5ZkfLM/Z/mLVppx4b0x7N+r5576d98whkdzy09MfC7Q/pOwdIN3d6wCfA0MBzKwW0ANIC44ZZWZJwTFPAwOBmsHPoecUEZFidtSevrvPNbPUQ8rezvV0AXBNsN0ZmOjue4EvzGwt0MTM1gOnuft8ADN7HugCzIj1AkSOJRMm/5NnxvwNgJNKRvpUGRkZjBw5MpHNEslRFDdy+wOTgu2zibwIHJAVlO0Ltg8tj8rMBhJ5V8A555xTBE0UiY+e3fvQs3sfAM2Xl2NSTDdyzex+YD8w/kBRlGp+hPKo3H20uzdy90YH5jCLiEjsCt3TN7O+RG7wtsl1QzYLqJarWlVgU1BeNUq5iIjEUaF6+mbWHhgCXO3uuRcqeQXoYWalzKwGkRu2C919M7DTzJpa5BMq1wPTY2y7iIgUUH6mbE4AWgEVzCwLGEZktk4p4J3gU4YL3P0Wd19pZpOBVUSGfQa5e3Zwqlv575TNGegmrhzqobJRiysV8nRbbzj6wmPr16+nU6dOrFixopBRRI4v+Zm90zNK8bNHqD8cGB6lPBM4fOFyERGJG30iV0IvOzubm266ibS0NNq2bcuPP/7I0qVLubJLG1q1b06/gb357vtvAeh6XUcefHgonbtfSYs2jVnyyWJuuLk3TVvV59HHH8k555Rpk2jSpAn16tXj5ptvJjs7O6/wInGlpC+ht2bNGgYNGsTKlSs5/fTTmTp1Ktdffz0P3vtbZr/5IRddWIvHn3wsp37J5GSmT57B9b370/emXox45E/MeWsBk6a8wI5vd/D52tW8/NpLfPDBByxdupSkpCTGjx9/hBaIxI8WXJPQq1GjBvXq1QOgYcOGrFu3ju+++47mTVsAcF23ngz4Zd+c+u2CNewvurAWF5x/Yc56PNXPSWXTpiw+ylzAsuVLadw4stTxjz/+SKVKhb0zIVK0lPQl9EqVKpWznZSUxHfffXfE+iVLlgSghJWgZMn/HlvCSrA/ez/uTvduPXlq1BPF02CRGGh4R+QQZcuWpVy5cixY+CEAL740iWYXt8j38S0zLuW1GdPZujXyBSc7duzgyy+/LJa2ihSUevpy7Hjo+6jFxb2efjTjxo3jxhtu4sc9P1K9Wip/fjz/a+dcUPNC7v31A7Rt25affvqJ5ORkRo4cSfXq1YuxxSL5o6QvoZaamnrQHP277rorZ3vGyzMPqz9t0us52xnNWpLRrGXUfV2u6sbA/7mhqJsrEjMN74iIhIiSvohIiCjpi4iEiJK+iEiIKOmLiISIkr6ISIhoyqYcM2qPq12k55vZ6oMiOc+b77zB52s+4/Zf3skf//dRypQpwy8H3s5jTwynaZPmXNriMv7+7Ch+0asfJ5c+uUhiihQXJX2Ro2h/RQfaB+vt5Dbkzvtztp957mmu6XpdgZJ+dnY2SUlJRdJGkfzS8I6E2vr167nwwgsZMGAA6enp9O7dm3fffZeMjAyatqrPx0sXM/HF8Qz9zV2HHXv7r2/l1Tde5pkxf2PL1s38vGcnuvboBMA9999Bo0aNSEtLY9iwYTnHpKam8vDDD9OiRQtGjBhBgwYNcvatWbOGhg0bFv9FS6gp6UvorV27lsGDB7Ns2TI+++wzXnjhBebNm8ew+37Hn0f+6ajH33TDLZxVqTIvTXiNaRNfA2Do3Q+SmZnJsmXLmDNnDsuWLcupn5KSwrx587j//vspW7YsS5cuBWDMmDH069evWK5R5AAlfQm9GjVqULt2bUqUKEFaWhpt2rTBzLjowlpszNpQqHO+8to0GjRoQP369Vm5ciWrVq3K2XfdddflbA8YMIAxY8aQnZ3NpEmT6NWrV8zXI3IkSvoSermXVi5RokTO8wNLJRfUlxvXM+qZvzBz5kyWLVtGx44d2bNnT87+MmXK5Gx369aNGTNm8Nprr9GwYUPKly8fw5WIHJ2SvkgRKHPKKezavROAXTt3cnLpMpQtW5avv/6aGTNm5HlcSkoK7dq149Zbb+WGG7RAmxQ/zd6RY8byvsujlidiaeWC+kXPfvTqew2VKp3FtImvUTutDmlpaZx77rlkZGQc8djevXvz0ksv0bZt2zi1VsLsqEnfzJ4DOgFb3T09KDsDmASkAuuB7u7+bbBvKHAjkA3c7u5vBeUNgbFAaeANYLC7e9FejkjBHLq08tixY3O2z6lWnblvLwCgx7W9Abj7jqE5+5/609M52wP63cyAfjcftK9S9dMOi7d+/frDyubNm0f//v01fVPiIj/DO2OB9oeU3QvMdPeawMzgOWZWC+gBpAXHjDKzA3/JTwMDgZrBz6HnFAmdrl278vzzzzN48OBEN0VC4qg9fXefa2aphxR3BloF2+OA2cCQoHyiu+8FvjCztUATM1sPnObu8wHM7HmgC5D3YKdICEybNi3RTZCQKeyN3DPdfTNA8FgpKD8b2JirXlZQdnawfWh5VGY20MwyzSxz27ZthWyiiIgcqqhn71iUMj9CeVTuPtrdG7l7o4oVKxZZ40REwq6wSf9rM6sMEDxuDcqzgGq56lUFNgXlVaOUi4hIHBU26b8C9A22+wLTc5X3MLNSZlaDyA3bhcEQ0E4za2pmBlyf6xgREYmT/EzZnEDkpm0FM8sChgEjgMlmdiOwAbgWwN1XmtlkYBWwHxjk7tnBqW7lv1M2Z6CbuHKITy+8qEjPV/6tj4rsXBNfHM8ny5fw6MOPx3yu1NRUMjMzqVChQhG0TKRg8jN7p2ceu9rkUX84MDxKeSaQXqDWiYhIkdIyDBJqu3fvpmPHjtStW5f09HQmTZrEokWLaN68OZe1z6Bd58vYtSuyvMKWrzfT4/qf07RVfR5+9MGcc7w0fQqXtmvGJW2b8sijv8kpnzBhArVr1yY9PZ0hQ4bE/dpEotEyDBJqb775JlWqVOH1118H4Pvvv6d+/fpMmjSJ6pUuYOfOH0hJKQ3AilXLmfn6+5QsWYqMNo24se/NJCUl8bvHhvH2q3M4vezpdP9FF9546zUa1GvIkCFDWLx4MeXKlaNt27a8/PLLdOnSJZGXK3J8J/3WswfFeIZPi6QdcvyqXbs2d911F0OGDKFTp06cfvrpVK5cmcaNG7P1yx849dT/LqXQMuNSTjutLADn/+wCsr7ayI5vd9C8aQsqlI+Mz3fr0p0FCz/EzGjVqhUHphz37t2buXPnKulLwml4R0Lt/PPPZ/HixdSuXZuhQ4cybdo0IhPMDleqZK4lmJOS2L9/P3ktH6VlpeRYpaQvobZp0yZOPvlk+vTpw1133cWCBQvYtGkTixYtAmDXrp3s35/3mvoN6zdi/kcfsH3HdrKzs5n2yhSaXZxBw/qNmDNnDt988w3Z2dlMmDCBSy+9NF6XJZKn43p4R04sF30WfbitOJdWXr58OXfffTclSpQgOTmZp59+GnfntttuY+f3u0lJSWHK+Lw/UnJmpbO4755h/LxnJ9ydyy+7givbdgTg0Ucf5bLLLsPd6dChA507dy626xDJLyV9CbV27drRrl27w8oXLFhw0ItNj2t75yyvDDD+uck52906X0u3ztcedo5evXpF/frDaMsri8SLhndEREJESV9EJESU9EVEQkRJX0QkRJT0RURCRElfRCRENGVTjhkjb3mvSM937dBGRXo+kROBevoiR1CjVpWo5eP++SyTp04AImvtb/l6c7G3ZcPGL7mkbVMAPpj/Pr37dy/2mHLiUU9fpBD69rkxZ3vSlBe48IKLOOvMygls0ZHt37+fk07Sf3dRT19C7g9/+ANPPfUUAHfccQetW7cGYObMmfzyVzcB8Ps/Psxl7TO4sksbtm6LfB30H//3UUaNfopX33iZpcuX8Mtf3UTrK1vw454f+WT5Erp070DDhg1p164dmzfn/S5g7dq1XH755dStW5cGDRqwbt063J3f/v4BLmnblEvbNePlV6ce8Ro+XrqYjj+/gjYdWtDx51ewevVqAMaOHcu1117LVVddRdu2bWP+XcmJQUlfQu2SSy7h/fffByAzM5Ndu3axb98+5s2bR9PGzfj3v3fTsH5jZr35Ac0ubs4/J4476PirOnShXu36jHryGd6bMY+Tkk7ivmH38I+nn2fx4sX079+f+++/P8/4vXv3ZtCgQXzyySd8+OGHVK5cmZdeeokVq5Yza8YHvPjP6Tz86G/4euuWPM9R87yaTJ88g5lvzOOeO+/jvvvuy9k3f/58xo0bx3vvFe39Ejl+6f2ehFrDhg1ZvHgxO3fupFSpUjRo0IDMzEzef/99ht37e0qWLEnbNu0BqJNejznzZh3xfGv/tYbPPv+U7n26cFLJEmRnZ1O5cvRhn507d/LVV1/RtWtXAFJSUgCYN28eXa++hqSkJCpVrESzizNY8snH1LowLep5ftj5A7f9+lb+tX4dZoaTnbPviiuu4Iwzzijw70VOXEr6EmrJycmkpqYyZswYmjdvTp06dZg1axbr1q3j/J9dwEknJeesr5+UlET2/uwjns/duaDmhbwx7V0qVT/tqHULUp6Xx/40nIxmLRk7ejwbNn7JNX2uytlXpkyZAp1LTnxK+nLMGPS31lHLi3NpZYgM8Tz++OM899xz1K5dmzvvvJOGDRvm+WUqhypzyins2rULgJ+dW5PtO75h0eKFdKx+Ofv27ePzzz8nLe3wXvppp51G1apVc75Gce/evWRnZ3PJJZfwlydHcl23Xnz73bcsWPghw4Y+wp69e6LG/2HnD1Q+K/JuYtKUFwr5W5CwiGlM38zuMLOVZrbCzCaYWYqZnWFm75jZmuCxXK76Q81srZmtNrPD17MVSYCWLVuyefNmmjVrxplnnklKSgotW7bM9/E9runFPQ/cQesrW5D9UzbPjnqe3z02jLp161KvXj0+/PDDPI/9v//7P5566inq1KlD8+bN2bJlC127dqXWRWlcdmUG3XpdxYP3/pZKlc7M8xyDbh7M8D/8lk7d2pL905HfiYgUuqdvZmcDtwO13P1HM5sM9ABqATPdfYSZ3QvcCwwxs1rB/jSgCvCumZ3v7vorlYRq06YN+/bty3n++eefA5F3GF+s2pRTflWHLlzVIfIdt3ffMTSnvNOVnel05X+/ICU9rQ7TJ8846vAOQM2aNaPeZB123+8Ydt/vDio7p1p15r69AICMZi3JaBZ5YWrcsAnzZ32cU++Jp/4AQL9+/ejXr99R2yDhEuvwzklAaTPbB5wMbAKGAq2C/eOA2cAQoDMw0d33Al+Y2VqgCTA/xjYkRGxfyq4vZBeRxCj08I67fwU8DmwANgPfu/vbwJnuvjmosxmoFBxyNrAx1ymygrLDmNlAM8s0s8xt27YVtokix4xBgwZRr169g37GjBmT6GZJCMUyvFOOSO+9BvAd8KKZ9TnSIVHKok5TcPfRwGiARo0aFWwqgxxX3D3fN0yPZyNHjkx0Ew5S0BlCcuKI5Ubu5cAX7r7N3fcBLwHNga/NrDJA8Lg1qJ8FVMt1fFUiw0ESUikpKWzfvl0JKM7cne3bt+d8LkDCJZYx/Q1AUzM7GfgRaANkAruBvsCI4HF6UP8V4AUze4LIjdyawMIY4stxrmrVqmRlZXG0Ibyd26NPVcyP7f8ufGKLJW4iY+cnbkpKClWrVi10DDl+FTrpu/tHZjYF+BjYDywhMiRzCjDZzG4k8sJwbVB/ZTDDZ1VQf5Bm7oRbcnIyNWrUOGq9WJZczmvuf37EutRzomLHEldOfDHN3nH3YcCwQ4r3Eun1R6s/HBgeS0wRESk8LbgmIhIiSvoiIiGipC8iEiJK+iIiIaKkLyISIkr6IiIhoqQvIhIiSvoiIiGib846zmhJZxGJhXr6IiIhoqQvIhIiSvoiIiGipC8iEiJK+iIiIaKkLyISIkr6IiIhoqQvIhIiSvoiIiGipC8iEiJK+iIiIaKkLyISIjElfTM73cymmNlnZvapmTUzszPM7B0zWxM8lstVf6iZrTWz1WbWLvbmi4hIQcTa0/8z8Ka7XwjUJbKM473ATHevCcwMnmNmtYAeQBrQHhhlZkkxxhcRkQIodNI3s9OAS4BnAdz9P+7+HdAZGBdUGwd0CbY7AxPdfa+7fwGsBZoUNr6IiBRcLOvpnwtsA8aYWV1gMTAYONPdNwO4+2YzqxTUPxtYkOv4rKDsMGY2EBgIcM4558TQRClKWstf5PgXy/DOSUAD4Gl3rw/sJhjKyYNFKfNoFd19tLs3cvdGFStWjKGJIiKSWyxJPwvIcvePgudTiLwIfG1mlQGCx6256lfLdXxVYFMM8UVEpIAKnfTdfQuw0cwuCIraAKuAV4C+QVlfYHqw/QrQw8xKmVkNoCawsLDxRUSk4GL9jtzbgPFmVhL4F3ADkReSyWZ2I7ABuBbA3Vea2WQiLwz7gUHunh1jfBERKYCYkr67LwUaRdnVJo/6w4HhscQUEZHC0ydyRURCRElfRCRElPRFREJESV9EJESU9EVEQiTWKZsixU7LP4gUHfX0RURCRElfRCRElPRFREJESV9EJESU9EVEQkRJX0QkRJT0RURCRElfRCRElPRFREJEn8gVyUNsnwQGfRpYjkXq6YuIhIiSvohIiCjpi4iEiJK+iEiIKOmLiIRIzEnfzJLMbImZvRY8P8PM3jGzNcFjuVx1h5rZWjNbbWbtYo0tIiIFUxQ9/cEcPDftXmCmu9cEZgbPMbNaQA8gDWgPjDKzpCKILyIi+RRT0jezqkBH4B+5ijsD44LtcUCXXOUT3X2vu38BrAWaxBJfREQKJtae/pPAPcBPucrOdPfNAMFjpaD8bGBjrnpZQdlhzGygmWWaWea2bdtibKKIiBxQ6KRvZp2Are6+OL+HRCnzaBXdfbS7N3L3RhUrVixsE0VE5BCxLMOQAVxtZh2AFOA0M/sn8LWZVXb3zWZWGdga1M8CquU6viqwKYb4IiJSQIXu6bv7UHev6u6pRG7QvufufYBXgL5Btb7A9GD7FaCHmZUysxpATWBhoVsuIiIFVhwLro0AJpvZjcAG4FoAd19pZpOBVcB+YJC7ZxdDfBERyUORJH13nw3MDra3A23yqDccGF4UMUWk6Gll0ROfPpErIhIiWk9f5BgUW49bvW3Jm3r6IiIhoqQvIhIiSvoiIiGipC8iEiJK+iIiIXJcz97pPjS25i8vonaIiBwv1NMXEQmR47qnv/yLDQmLHcu7DL3DEJFEOa6TfiIl8gVHRKSwlPSPM3qHISKxUNI/zugdhojEQklfRHIkckac1huKDyV9OeZpSCt+9E7yxKekL/mWqOSbqESkz4HIiUjz9EVEQkQ9fcm3sL31D9v1Svwl4j6Gkr6ISIIkYshUSV9EJEES8W5SSV9EQi1sU0ULnfTNrBrwPHAW8BMw2t3/bGZnAJOAVGA90N3dvw2OGQrcCGQDt7v7WzG1XkROGJqaGx+x9PT3A79294/N7FRgsZm9A/QDZrr7CDO7F7gXGGJmtYAeQBpQBXjXzM539+zYLkHkxBPGBKgb5/FR6L8sd98MbA62d5rZp8DZQGegVVBtHDAbGBKUT3T3vcAXZrYWaALML2wbRE5USoDxE7YX2CKZp29mqUB94CPgzOAF4cALQ6Wg2tnAxlyHZQVlIiISJzHfyDWzU4CpwK/c/Qczy7NqlDLP45wDgYEA55xzTqxNFBHJU9jeVcXU0zezZCIJf7y7vxQUf21mlYP9lYGtQXkWUC3X4VWBTdHO6+6j3b2RuzeqWLFiLE0UEZFcYpm9Y8CzwKfu/kSuXa8AfYERweP0XOUvmNkTRG7k1gQWFja+iBS91D0vxHT8+qJphhSjWIZ3MoBfAMvNbGlQdh+RZD/ZzG4ENgDXArj7SjObDKwiMvNnkGbuiIjEVyyzd+YRfZweoHXzOPcAAAdoSURBVE0exwwHhhc2pojIiSSWd1brC3mcVtkUEQkRJX0RkRBR0hcRCRElfRGRENEqmyJ5SOT0xUTc4Eu0MF5zIijpyzFPyUCKU9j+vjS8IyISIurpS76FrUckciJS0i+kRCVAJV4RiYWGd0REQkRJX0QkRJT0RURC5Lge09cysCIiBaOevohIiCjpi4iEiJK+iEiIKOmLiISIkr6ISIgo6YuIhIiSvohIiCjpi4iEiJK+iEiIxD3pm1l7M1ttZmvN7N54xxcRCbO4Jn0zSwJGAlcCtYCeZlYrnm0QEQmzePf0mwBr3f1f7v4fYCLQOc5tEBEJLXP3+AUzuwZo7+4Dgue/AC529/85pN5AYGDw9AJgdSFDVgC+KeSxsUhU3ETG1jWHI3bY4iYydqxxq7t7xUML473KpkUpO+xVx91HA6NjDmaW6e6NYj3P8RI3kbF1zeGIHba4iYxdXHHjPbyTBVTL9bwqsCnObRARCa14J/1FQE0zq2FmJYEewCtxboOISGjFdXjH3feb2f8AbwFJwHPuvrIYQ8Y8RHScxU1kbF1zOGKHLW4iYxdL3LjeyBURkcTSJ3JFREJESV9EJEROyKSfqKUezOw5M9tqZiviFTOIW83MZpnZp2a20swGxzF2ipktNLNPgti/jVfsIH6SmS0xs9fiHHe9mS03s6VmlhnHuKeb2RQz+yz4924Wp7gXBNd64OcHM/tVnGLfEfxtrTCzCWaWEqe4g4OYK4v7WqPlDjM7w8zeMbM1wWO5Ignm7ifUD5EbxOuAc4GSwCdArTjFvgRoAKyI8zVXBhoE26cCn8fxmg04JdhOBj4Cmsbx2u8EXgBei/PvfD1QIZ4xg7jjgAHBdkng9AS0IQnYQuTDP8Ud62zgC6B08Hwy0C8OcdOBFcDJRCa8vAvULMZ4h+UO4A/AvcH2vcBjRRHrROzpJ2ypB3efC+yIR6xD4m5294+D7Z3Ap0T+s8Qjtrv7ruBpcvATl9kBZlYV6Aj8Ix7xEs3MTiOSHJ4FcPf/uPt3CWhKG2Cdu38Zp3gnAaXN7CQiSTgen+25CFjg7v929/3AHKBrcQXLI3d0JvIiT/DYpShinYhJ/2xgY67nWcQpAR4LzCwVqE+kxx2vmElmthTYCrzj7vGK/SRwD/BTnOLl5sDbZrY4WDYkHs4FtgFjgiGtf5hZmTjFzq0HMCEegdz9K+BxYAOwGfje3d+OQ+gVwCVmVt7MTgY6cPAHS+PhTHffDJGOHVCpKE56Iib9fC31cCIys1OAqcCv3P2HeMV192x3r0fkE9ZNzCy9uGOaWSdgq7svLu5Yechw9wZEVowdZGaXxCHmSUSGAJ529/rAbiJv++Mm+FDl1cCLcYpXjkiPtwZQBShjZn2KO667fwo8BrwDvElkmHh/cceNhxMx6YdyqQczSyaS8Me7+0uJaEMw1DAbaB+HcBnA1Wa2nsgQXmsz+2cc4gLg7puCx63ANCLDisUtC8jK9U5qCpEXgXi6EvjY3b+OU7zLgS/cfZu77wNeAprHI7C7P+vuDdz9EiJDL2viETeXr82sMkDwuLUoTnoiJv3QLfVgZkZknPdTd38izrErmtnpwXZpIv9JPyvuuO4+1N2runsqkX/j99y92HuAAGZWxsxOPbANtCUyHFCs3H0LsNHMLgiK2gCrijvuIXoSp6GdwAagqZmdHPydtyFyz6rYmVml4PEc4OfE97ohkrf6Btt9gelFcdJ4r7JZ7Dz+Sz3kMLMJQCuggpllAcPc/dk4hM4AfgEsD8bWAe5z9zfiELsyMC74gpwSwGR3j+v0yQQ4E5gWyUGcBLzg7m/GKfZtwPigQ/Mv4IY4xSUY274CuDleMd39IzObAnxMZHhlCfFbFmGqmZUH9gGD3P3b4goULXcAI4DJZnYjkRe/a4skVjAdSEREQuBEHN4REZE8KOmLiISIkr6ISIgo6YuIhIiSvohIiCjpS6iZWXawauQKM3sxmJZYkOP/YWa1ClC/n5n9teAtFSkaSvoSdj+6ez13Twf+A9yS3wPNLMndB7h7vD8gJVJoSvoi//U+8DMAM+sTfE/AUjP7e/DhM8xsl5k9bGYfAc3MbLaZNQr29QzW2F9hZo8dOKmZ3WBmn5vZHCIfpBNJGCV9ESBYtvdKIp9qvgi4jsiiavWAbKB3ULUMkTXPL3b3ebmOr0Jkga7WQD2gsZl1CdZM+S2RZH8FkO+hIJHicMItwyBSQKVzLV3xPpE1jAYCDYFFwVILpfnvYlfZRBa2O1RjYLa7bwMws/FE1r7nkPJJwPnFcB0i+aKkL2H3Y9CbzxEs7DXO3YdGqb/H3bOjlEdb0vsArXUixwwN74gcbiZwTa5VFs8ws+pHOeYj4FIzqxCM//ck8m1LHwGtgi/jSKaIFs0SKSz19EUO4e6rzOwBIt+MVYJglUUgz68HdPfNZjYUmEWk1/+Gu08HMLOHgPlEvvnpYyKrv4okhFbZFBEJEQ3viIiEiJK+iEiIKOmLiISIkr6ISIgo6YuIhIiSvohIiCjpi4iEyP8Dhkk2l+ZT52UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_2.groupby(\"Period\").Choice.value_counts().unstack().plot.bar(stacked=True, rot=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. one-step-ahead simulation\n",
    "\n",
    "The one-step-ahead simulation takes the full dataset and simulates the outcomes for each observation. Thus, the original sample size is preserved. The results are used to measure the within-sample fit of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = rp.get_simulate_func(params, options, \"one_step_ahead\", df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_3 = simulate(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEeCAYAAACXAbTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXyU5bn/8e+VhC0syhIWEzZLAiQoRRB3UKiKx55itSpuoD+XQrHUpUdF/aHtKb7qOT3+FPcVca9iW6w9WldwRwOKgSwsgoAEiIIsskiS6/fHPNgYJ4RknsxMks/79corM/c8z/1cM4G55nru+7nH3F0AAAAAgNilJDoAAAAAAGgqKLAAAAAAICQUWAAAAAAQEgosAAAAAAgJBRYAAAAAhIQCCwAAAABCQoEFJICZXW9mD4W97X705WbWL4y+AACNg5mdZ2avNFDfj5rZHxqi7yjHutDM3onHsYBYUGABIQje9AvMbIeZrTeze83swJq2d/db3P2S/em7LtsCAJovMzvWzN4zsy1mtsnM3jWzw939SXc/KQnim2tm5DM0eRRYQIzM7GpJt0r6D0kHSDpSUm9Jr5pZyyjbp8U3QgBAU2dmHSS9KOlOSZ0kZUr6naTdiYwLaI4osIAYBAntd5J+7e4vu/sed18l6SxFiqzzzexmM5ttZk+Y2VZJFwZtT1TpZ7yZfW5mX5nZ/zWzVWb2k+Cx77Y1sz7BNL8JZrbazL40sxuq9DPczN43s6/NrNTM7opW5AEAmpwcSXL3p929wt13uvsr7v5p9al1QR75lZktM7NtZvafZvajIH9sNbNn9+aOaNPyappubmYdzexFMyszs83B7azgsemSjpN0l5ltN7O7gvYBZvZqMOJWYmZnVemvs5m9EMT0oaQfNcDrBoSOAguIzdGSWkv6S9VGd98u6SVJJwZNYyXNlnSgpCerbmtmuZLukXSepB6KjIJl1nLcYyX1lzRa0jQzGxi0V0i6UlIXSUcFj/+qHs8LANC4LJVUYWazzOwUM+tYy/ZjJA1VZNbFNZIeUCQP9ZQ0SNI59YghRdJMRU4w9pK0U9JdkuTuN0h6W9Ll7t7O3S83s7aSXpX0lKSuwTHvMbO8oL+7Je1SJDf+n+AHSHoUWEBsukj60t3LozxWGjwuSe+7+9/cvdLdd1bb7heS/u7u77j7t5KmSfJajvu74OzkIkmLJA2WJHdf4O4fuHt5MJJ2v6SR9XtqAIDGwt23KnLyzSU9KKksGP3pVsMut7r7VndfImmxpFfc/TN336LICcIh9YjhK3d/3t13uPs2SdO17xz0U0mr3H1mkLcWSnpe0i/MLFXSGZKmufs37r5Y0qy6xgQkAgUWEJsvJXWp4bqqHsHjkrRmH30cVPVxd98h6atajru+yu0dktpJkpnlBFMy1gfTEW/Rv4o8AEAT5u5F7n6hu2cpMgp1kKTba9h8Q5XbO6Pcb1fX45tZupndH0x53yrpLUkHBsVSNL0lHRFMa//azL5WZBStu6QMSWn6fv78vK4xAYlAgQXE5n1FLiA+vWpjMO3hFEmvB037GpEqlZRVZd82kjrXM557JRVLynb3DpKul2T17AsA0Ei5e7GkRxUptGLxjaT0vXfMrPs+tr1akenrRwQ5aMTe3faGVW37NZLmufuBVX7aufskSWWSyhWZsrhXrxieBxA3FFhADIKpFL+TdKeZjTGzFmbWR9JzktZKenw/upkt6d/N7OjgouLfqf5FUXtJWyVtN7MBkibVsx8AQCMSLBZxdZVFJXoqck3TBzF2vUhSnpn92MxaS7p5H9u2V2T062sz6yTppmqPb5B0cJX7L0rKMbMLgvzZwswON7OB7l6hyPXNNwcjY7mSJsT4XIC4oMACYuTu/6XISNGfFClu5ityVm60u9e6PG4w//3Xkp5RZDRrm6SNqt/Sur+VdG7Qx4OS/lyPPgAAjc82SUdImm9m3yhSWC1WZFSp3tx9qaTfS3pN0jJJ+/qi39sltVFkevwHkl6u9vgdilxftdnMZgTXaZ0kaZykdYpMf79VUqtg+8sVmaq4XpHRuJmxPBcgXsy9tmvpAcSTmbWT9LUi0/xWJjoeAAAA7D9GsIAkYGb/HkyBaKvISFiBpFWJjQoAAAB1RYEFJIexikyPWCcpW9I4Z3gZAACg0WGKIAAAAACEhBEsAAAAAAgJBRYAAAAAhCQt0QHUpkuXLt6nT59EhwEAiKMFCxZ86e4ZiY5jf5CnAKB5qilXJX2B1adPH+Xn5yc6DABAHJnZ54mOYX+RpwCgeaopV9U6RdDMHjGzjWa2OMpjvzUzN7MuVdqmmtlyMysxs5OrtA81s4LgsRlmZvV9MgAAVEWuAgAki/25ButRSWOqN5pZT0knSlpdpS1XkW/jzgv2ucfMUoOH75V0mSJLUGdH6xMAgHp6VOQqAEASqLXAcve3JG2K8tD/k3SNpKrrvI+V9Iy773b3lZKWSxpuZj0kdXD394Pv9nlM0mkxRw8AgMhVAIDkUa9VBM3sZ5K+cPdF1R7KlLSmyv21QVtmcLt6OwAADYJcBQBIhDovcmFm6ZJukHRStIejtPk+2ms6xmWKTNFQr1696hoiAKCZa+hcRZ4CANSkPiNYP5LUV9IiM1slKUvSQjPrrsjZvp5Vts2StC5oz4rSHpW7P+Duw9x9WEZGo1ilFwCQXBo0V5GnAAA1qXOB5e4F7t7V3fu4ex9FEtJh7r5e0guSxplZKzPrq8gFwh+6e6mkbWZ2ZLAi03hJc8J7GgAA/Au5CgCQKLVOETSzpyUdL6mLma2VdJO7PxxtW3dfYmbPSiqUVC5psrtXBA9PUmSVpzaSXgp+Ynb3xDfC6EaT7xsVSj8AgPhrDrmKPAUAjUOtBZa7n1PL432q3Z8uaXqU7fIlDapjfLUaNXdySD0VhdQPACDekj1XAQCajzovcoGaMZoGAEh2jKYBQMOq1zLtAAAAAIAfYgQrRMk0XZHRNABAsmM0DUBTRIGFBkcCBQAkO3IVgLBQYKFZIYECAJIZM1CAxo9rsAAAAAAgJBRYAAAAABASCiwAAAAACAnXYAEJwvVgAIBkxvVgQP0wggUAAAAAIaHAAgAAAICQMEUQaOaYAgIASHbkKjQmjGABAAAAQEgosAAAAAAgJBRYAAAAABASrsFqokbNnRxST0Uh9QMAAAA0fYxgAQAAAEBIGMECkDRYJQoAkOzCyFXkqaaNESwAAAAACAkjWGhw4VwPxrVgAAAASH61jmCZ2SNmttHMFldp+28zKzazT83sr2Z2YJXHpprZcjMrMbOTq7QPNbOC4LEZZmbhPx0AQHNErgIAJIv9mSL4qKQx1dpelTTI3Q+VtFTSVEkys1xJ4yTlBfvcY2apwT73SrpMUnbwU71PAADq61GRqwAASaDWAsvd35K0qVrbK+5eHtz9QFJWcHuspGfcfbe7r5S0XNJwM+shqYO7v+/uLukxSaeF9SQAAM0buQoAkCzCWOTi/0h6KbidKWlNlcfWBm2Zwe3q7VGZ2WVmlm9m+WVlZSGECABo5kLNVeQpAEBNYiqwzOwGSeWSntzbFGUz30d7VO7+gLsPc/dhGRkZsYQIAGjmGiJXkacAADWp9yqCZjZB0k8ljQ6mUkiRs309q2yWJWld0J4VpR0AgAZDrgIAxFu9CiwzGyPpWkkj3X1HlYdekPSUmd0m6SBFLhD+0N0rzGybmR0pab6k8ZLujC10oO5YMh5oPshVAIBEqLXAMrOnJR0vqYuZrZV0kyIrMbWS9Gqwgu0H7j7R3ZeY2bOSChWZjjHZ3SuCriYpsspTG0Xmwb8kAABCQK4CACSLWgssdz8nSvPD+9h+uqTpUdrzJQ2qU3QAAOwHchUAIFmEsYogAAAAAEAUWAAAAAAQGgosAAAAAAhJvZdpBxAbVjQEAABoehjBAgAAAICQUGABAAAAQEgosAAAAAAgJBRYAAAAABASCiwAAAAACAkFFgAAAACEhAILAAAAAEJCgQUAAAAAIaHAAgAAAICQUGABAAAAQEgosAAAAAAgJBRYAAAAABASCiwAAAAACAkFFgAAAACEJC3RAQBIrFFzJ4fUU1FI/QAAADRetY5gmdkjZrbRzBZXaetkZq+a2bLgd8cqj001s+VmVmJmJ1dpH2pmBcFjM8zMwn86AIDmiFwFAEgW+zNF8FFJY6q1XSfpdXfPlvR6cF9mlitpnKS8YJ97zCw12OdeSZdJyg5+qvcJAEB9PSpyFQAgCdRaYLn7W5I2VWseK2lWcHuWpNOqtD/j7rvdfaWk5ZKGm1kPSR3c/X13d0mPVdkHAICYkKsAAMmivotcdHP3UkkKfncN2jMlramy3dqgLTO4Xb09KjO7zMzyzSy/rKysniECAJq5BstV5CkAQE3CXkUw2lx130d7VO7+gLsPc/dhGRkZoQUHAIBCyFXkKQBATepbYG0IplIo+L0xaF8rqWeV7bIkrQvas6K0AwDQUMhVAIC4q2+B9YKkCcHtCZLmVGkfZ2atzKyvIhcIfxhMzdhmZkcGKzKNr7IPAAANgVwFAIi7Wr8Hy8yelnS8pC5mtlbSTZL+KOlZM7tY0mpJZ0qSuy8xs2clFUoqlzTZ3SuCriYpsspTG0kvBT8AAMSMXAUASBa1Fljufk4ND42uYfvpkqZHac+XNKhO0QEAsB/IVQCAZBH2IhcAAAAA0GxRYAEAAABASGqdIggAAGIzau7kEHopCqEPAEBDYwQLAAAAAEJCgQUAAAAAIaHAAgAAAICQcA0WAERx98Q3Quln8n2jQukHAIDqwshV5KnwMYIFAAAAACFhBAtA0ghnpTWJ1dYAAECiMIIFAAAAACGhwAIAAACAkFBgAQAAAEBIKLAAAAAAICQUWAAAAAAQElYRBBLkrKmx//crSJI4pHBiAQAAaOwYwQIAAACAkFBgAQAAAEBImCKIZiVZpuUBAACgaWIECwAAAABCElOBZWZXmtkSM1tsZk+bWWsz62Rmr5rZsuB3xyrbTzWz5WZWYmYnxx4+AAD7Rq4CAMRTvQssM8uUNEXSMHcfJClV0jhJ10l63d2zJb0e3JeZ5QaP50kaI+keM0uNLXwAAGpGrgIAxFusUwTTJLUxszRJ6ZLWSRoraVbw+CxJpwW3x0p6xt13u/tKScslDY/x+AAA1IZcBQCIm3pf8e/uX5jZnyStlrRT0ivu/oqZdXP30mCbUjPrGuySKemDKl2sDdrQxLGwBIBEIVcBAOKt3p98g/nqYyX1lfS1pOfM7Px97RKlzWvo+zJJl0lSr1696hsigEaGLz1G2BoqV5GnAAA1ieXTzE8krXT3Mkkys79IOlrSBjPrEZwR7CFpY7D9Wkk9q+yfpcg0jR9w9wckPSBJw4YNi1qEAQCwHxokVzXmPDVq7uQQeikKoQ8AaJpiuQZrtaQjzSzdzEzSaEXecV+QNCHYZoKkOcHtFySNM7NWZtZXUrakD2M4PgAAtSFXAQDiKpZrsOab2WxJCyWVS/pYkbN57SQ9a2YXK5LYzgy2X2Jmz0oqDLaf7O4VMcYPAECNyFUAgHiL6YIHd79J0k3VmncrcoYw2vbTJU2P5ZgAwlWwcnWiQwAaFLkKABBPsS7TDgAAAAAIhLNkF4A6Y+QIAACg6WEECwAAAABCQoEFAAAAACGhwAIAAACAkHANFpoVrntKbvx9AABAY8cIFgAAAACEhAILAAAAAEJCgQUAAAAAIeEarCbqrKnh/GkLwuiD62oAAADQTDCCBQAAAAAhocACAAAAgJBQYAEAAABASCiwAAAAACAkFFgAAAAAEBIKLAAAAAAICQUWAAAAAISEAgsAAAAAQsIXDQNAFKPmTg6pp6KQ+gEAAI0BBVYTVbBydaJDAAAAAJqdmKYImtmBZjbbzIrNrMjMjjKzTmb2qpktC353rLL9VDNbbmYlZnZy7OEDALBv5CoAQDzFOoJ1h6SX3f0XZtZSUrqk6yW97u5/NLPrJF0n6Vozy5U0TlKepIMkvWZmOe5eEWMMAADsC7kqSYUzFZdpuACSS71HsMysg6QRkh6WJHf/1t2/ljRW0qxgs1mSTgtuj5X0jLvvdveVkpZLGl7f4wMAUBtyFQAg3mKZIniwpDJJM83sYzN7yMzaSurm7qWSFPzuGmyfKWlNlf3XBm0/YGaXmVm+meWXlZXFECIAoJlrkFxFngIA1CSWAitN0mGS7nX3IZK+UWSKRU0sSptH29DdH3D3Ye4+LCMjI4YQAQDNXIPkKvIUAKAmsVyDtVbSWnefH9yfrUjS2mBmPdy91Mx6SNpYZfueVfbPkrQuhuMDQLNw98Q3Yu5j8n2jQoikUSJXAUAckKv+pd4jWO6+XtIaM+sfNI2WVCjpBUkTgrYJkuYEt1+QNM7MWplZX0nZkj6s7/EBAKgNuQoAEG+xriL4a0lPBqsyfSbpIkWKtmfN7GJJqyWdKUnuvsTMnlUksZVLmsyqTACS1VlTw/mawIJQekGMyFUAgLiJ6ROEu38iaViUh0bXsP10SdNjOSYAAHVBrgIAxFM4p2gBNFp9dj0VSj+rQukFAACgcaPAAgAAzR5fegwgLBRYAJIGo2kAmrtwCj2JYg9InFi+BwsAAAAAUAUjWAAAAPgBRtOA+qHAAoAoClauTnQIAIAAxR4aE6YIAgAAAEBIKLAAAAAAICQUWAAAAAAQEq7BAhIkjCXJV8UeBgAAqAO+Mw21YQQLAAAAAEJCgQUAAAAAIaHAAgAAAICQUGABAAAAQEgosAAAAAAgJKwiiGaFlfsAAADQkBjBAgAAAICQUGABAAAAQEgosAAAAAAgJDEXWGaWamYfm9mLwf1OZvaqmS0Lfnessu1UM1tuZiVmdnKsxwYAYH+QqwAA8RLGIhe/kVQkqUNw/zpJr7v7H83suuD+tWaWK2mcpDxJB0l6zcxy3L0ihBgAIFRhLIgisShKEiFXAQDiIqYRLDPLknSqpIeqNI+VNCu4PUvSaVXan3H33e6+UtJyScNjOT4AALUhVwEA4inWEazbJV0jqX2Vtm7uXipJ7l5qZl2D9kxJH1TZbm3QhiaOpdEBJBi5CgAQN/UewTKzn0ra6O4L9neXKG1eQ9+XmVm+meWXlZXVN0QAQDPXULmKPAUAqEksUwSPkfQzM1sl6RlJo8zsCUkbzKyHJAW/Nwbbr5XUs8r+WZLWRevY3R9w92HuPiwjIyOGEAEAzVyD5CryFACgJvUusNx9qrtnuXsfRS4IfsPdz5f0gqQJwWYTJM0Jbr8gaZyZtTKzvpKyJX1Y78gBAKgFuQoAEG9hrCJY3R8lPWtmF0taLelMSXL3JWb2rKRCSeWSJrMqEwDUbtTcySH0UhRCH00KuQoA0CBCKbDcfa6kucHtrySNrmG76ZKmh3FMAADqglwFAIiHhhjBAgAAANDAmOGQnGL6HiwAAAAAwL8wggUASe6sqbG/VReEEAcAAKgdBRYAAACwnzjphdowRRAAAAAAQsIIVojCOKMhhXNWo8+up0LoRVoVSi8AYlGwcnWiQwCAhEqmz1hAbSiwAABoYEwpAoDmgwILAIAGxigkADQfFFghIoECAICmgml5QP1QYAEAAAD7iRPqqA0FFgAAQJJg1Aho/CiwAABoRlhwAwAaFgUWAADNSDJNb6LYA9AUUWABAIBmL1mKvWQqgJMJrwsaEwosAEhyYXxx+KrYw5Ak3T3xjZj7mHzfqBAiQVPAh+bkxt8n+SXLiYFkEkaekmLLVRRYAACg2aOYABCWlEQHAAAAAABNBQUWAAAAAISEAgsAAAAAQsI1WAAAAEAjxLWDyaneBZaZ9ZT0mKTukiolPeDud5hZJ0l/ltRHkYWrznL3zcE+UyVdLKlC0hR3/2dM0QMAsA/kKgCIj1FzJ4fQS1EIfSReLFMEyyVd7e4DJR0pabKZ5Uq6TtLr7p4t6fXgvoLHxknKkzRG0j1mlhpL8AAA1IJcBQCIq3oXWO5e6u4Lg9vbFCk5MyWNlTQr2GyWpNOC22MlPePuu919paTlkobX9/gAANSGXAUAiLdQFrkwsz6ShkiaL6mbu5dKkcQmqWuwWaakNVV2Wxu0RevvMjPLN7P8srKyMEIEADRzYeYq8hQAoCYxL3JhZu0kPS/pCnffamY1bhqlzaNt6O4PSHpAkoYNGxZ1m73C+AZrqel9izUA4F/CzlV1yVMAmpY+u56KuY9VsYeBJBZTdWJmLRRJWE+6+1+C5g1m1sPdS82sh6SNQftaST2r7J4laV0sx5dYPQUAsG/JkKsAAM1HvacIWuT038OSitz9tioPvSBpQnB7gqQ5VdrHmVkrM+srKVvSh/U9PgAAtSFXAQDiLZYRrGMkXSCpwMw+Cdqul/RHSc+a2cWSVks6U5LcfYmZPSupUJFVnSa7e0UMxwcAxFkjXIaXXAUAiKt6F1ju/o6iz1WXpNE17DNd0vT6HhMAgLogVwFNQxjXPUlc+4T4CGeFCABAszBwHJcjAUCySKYFN8JYeK6pLDoXyjLtAAAAAABGsAAAABAF0/JQF6zs/S+MYAEAAABASBjBAgAASBKMGgGNHwUWAGC/JdMF1UCY+LcNICwUWAAAICEoagA0RRRYIWJYHwAAAGjeKLAAAAAAxCRZRqRHzZ0cQi+SVFTvPVlFEAAAAABCwggWAADNSLKcZQaApooCCwCABkZRAwDNBwUWAAAAgCbhrKnhlDcFMexLgQUAAACgSShYuTrRIbDIBQAAAACEhREsAAAAAE1CMnwvbaMvsJLhRQQAAAAAiSmCAAAAABAaCiwAAAAACEncCywzG2NmJWa23Myui/fxAQCoDbkKAFBfcS2wzCxV0t2STpGUK+kcM8uNZwwAAOwLuQoAEIt4j2ANl7Tc3T9z928lPSNpbJxjAABgX8hVAIB6i3eBlSlpTZX7a4M2AACSBbkKAFBv5u7xO5jZmZJOdvdLgvsXSBru7r+utt1lki4L7vaXVBLjobtI+jLGPsJCLNERyw8lSxwSsdSEWKILI5be7p4RRjB1tT+5qgHylNT0/oZhSJY4JGKpCbFElyyxJEscUtOMJWquivf3YK2V1LPK/SxJ66pv5O4PSHogrIOaWb67Dwurv1gQS3TEkrxxSMRSE2KJLpliqadac1XYeUpKrtctWWJJljgkYqkJsUSXLLEkSxxS84ol3lMEP5KUbWZ9zaylpHGSXohzDAAA7Au5CgBQb3EdwXL3cjO7XNI/JaVKesTdl8QzBgAA9oVcBQCIRbynCMrd/1fS/8b5sKFO44gRsURHLD+ULHFIxFITYokumWKpF3JV0sSSLHFIxFITYokuWWJJljikZhRLXBe5AAAAAICmLN7XYAEAAABAk9XkCywzG2NmJWa23MyuS2Acj5jZRjNbnKgYgjh6mtmbZlZkZkvM7DcJjKW1mX1oZouCWH6XqFiqxJRqZh+b2YsJjmOVmRWY2Sdmlp/gWA40s9lmVhz8uzkqQXH0D16PvT9bzeyKBMVyZfBvdrGZPW1mrRMRRxDLb4I4lsT79Yj2vmZmnczsVTNbFvzuGM+YGqNkyVNBLOSqH8aSVLkqWfJUEEtS5CryVI3xkKuUmFzVpAssM0uVdLekUyTlSjrHzHITFM6jksYk6NhVlUu62t0HSjpS0uQEvia7JY1y98GSfixpjJkdmaBY9vqNpKIEx7DXCe7+4yRY0vQOSS+7+wBJg5Wg18fdS4LX48eShkraIemv8Y7DzDIlTZE0zN0HKbIIwrh4xxHEMkjSpZKGK/K3+amZZccxhEf1w/e16yS97u7Zkl4P7qMGSZanJHJVNMmWq5IpT0nJkavIU9WQq77nUcU5VzXpAkuRP+Ryd//M3b+V9IyksYkIxN3fkrQpEceuFkepuy8Mbm9T5E0oM0GxuLtvD+62CH4SdlGgmWVJOlXSQ4mKIdmYWQdJIyQ9LEnu/q27f53YqCRJoyWtcPfPE3T8NEltzCxNUrqifJ9fnAyU9IG773D3cknzJP08Xgev4X1trKRZwe1Zkk6LVzyNVNLkKYlcVUMsSZOryFM/RJ7aJ3KVEpOrmnqBlSlpTZX7a5WgN+hkZGZ9JA2RND+BMaSa2SeSNkp61d0TFouk2yVdI6kygTHs5ZJeMbMFZnZZAuM4WFKZpJnBlJSHzKxtAuPZa5ykpxNxYHf/QtKfJK2WVCppi7u/kohYJC2WNMLMOptZuqR/0/e/IDcRurl7qRT5kCypa4LjSXbkqVqQq74nmfKUlBy5ijwVBbmqVg2aq5p6gWVR2lg2UZKZtZP0vKQr3H1rouJw94pgKD1L0vBgGDnuzOynkja6+4JEHD+KY9z9MEWmDU02sxEJiiNN0mGS7nX3IZK+UYKnfFnki19/Jum5BB2/oyJnvvpKOkhSWzM7PxGxuHuRpFslvSrpZUmLFJlahcaDPLUP5Kp/ScI8JSVHriJPRY+BXJVATb3AWqvvV8hZStzwaNIwsxaKJKwn3f0viY5HkoLh/LlK3Nz/YyT9zMxWKTJFZ5SZPZGgWOTu64LfGxWZvz08QaGslbS2ytna2YokskQ6RdJCd9+QoOP/RNJKdy9z9z2S/iLp6ATFInd/2N0Pc/cRikyBWJaoWAIbzKyHJAW/NyY4nmRHnqoBueoHkipPSUmTq8hT0ZGr9q1Bc1VTL7A+kpRtZn2DswnjJL2Q4JgSysxMkXnKRe5+W4JjyTCzA4PbbRR5MyhORCzuPtXds9y9jyL/Tt5w94Sc6TGztmbWfu9tSScpMrwed+6+XtIaM+sfNI2WVJiIWKo4RwmcdqHIdIsjzSw9+P80Wgm84NzMuga/e0k6XYl9baTIe+yE4PYESXMSGEtjQJ6Kglz1Q8mUp6TkyVXkqRqRq/atQXNVWpidJRt3LzezyyX9U5HVUx5x9yWJiMXMnpZ0vKQuZrZW0k3u/nACQjlG0gWSCoL55JJ0vbv/bwJi6SFpVrCKVoqkZ9094cvOJoFukv4aeT9UmqSn3P3lBMbza0lPBh/+PpN0UaICCeZunyjpl4mKwd3nm9lsSQsVmeLwsRL77fTPm1lnSXskTXb3zfE6cLT3NUl/lPSsmV2sSII/M17xNEbJlKckclUNyFXRJVOuIk9VQ676l0TkKnNnqjcAAAAAhKGpT7kJa+AAABKOSURBVBEEAAAAgLihwAIAAACAkFBgAQAAAEBIKLAAAAAAICQUWAAAAAAQEgosoIGZWYWZfWJmi83suWAJ17rs/5CZ5dZh+wvN7K66RwoAaI7IU0C4KLCAhrfT3X/s7oMkfStp4v7uaGap7n6Juyf6SxMBAE0XeQoIEQUWEF9vS+onSWZ2vpl9GJw1vD/4EkuZ2XYz+72ZzZd0lJnNNbNhwWPnmFlBcJbx1r2dmtlFZrbUzOYp8gWdAADUB3kKiBEFFhAnZpYm6RRJBWY2UNLZko5x9x9LqpB0XrBpW0mL3f0Id3+nyv4HSbpV0ihJP5Z0uJmdZmY9JP1OkYR1oqT9nqYBAMBe5CkgHGmJDgBoBtqY2SfB7bclPSzpMklDJX1kZpLURtLGYJsKSc9H6edwSXPdvUySzOxJSSOCx6q2/1lSTgM8DwBA00SeAkJEgQU0vJ3B2b/vWCRbzXL3qVG23+XuFVHabR/H8FgCBAA0a+QpIERMEQQS43VJvzCzrpJkZp3MrHct+8yXNNLMugTz4M+RNC9oP97MOptZC0lnNmTgAIBmgTwF1BMjWEACuHuhmd0o6RUzS5G0R9JkSZ/vY59SM5sq6U1FzhL+r7vPkSQzu1nS+5JKJS2UlNqwzwAA0JSRp4D6M3dGbAEAAAAgDEwRBAAAAICQUGABAAAAQEgosAAAAAAgJBRYAAAAABASVhEEQrRgwYKuaWlpD0kaJE5gAAAav0pJi8vLyy8ZOnToxlq3BkCBBYQpLS3toe7duw/MyMjYnJKSwhKdAIBGrbKy0srKynLXr1//kKSfJToeoDHgDDsQrkEZGRlbKa4AAE1BSkqKZ2RkbFFkZgaA/UCBBYQrheIKANCUBHmNz4zAfuI/CwAAAACEhGuwgAbU57p/DA2zv1V/PHVBmP01WzcfEOrfRTdvqfXvUlJS0vKnP/1p9rJly5aEeuw4OGTWIaG+XgUTCkL5d/zkk08esGTJkja33HLL+quuuuqgdu3aVfz+97/fcMUVVxx0/PHHbzvttNO2/f73v+965ZVXftm+ffvKMI4ZpqIBA0N9XQcWF4X2/jBjxozO+fn5bR977LHVsfaVmZl5SH5+flGPHj3Kw4htr7snvhHq6zf5vlG8vwIIBSNYQBNTUlLSMjs7O696+/Dhw/u/9dZb6YmISYp8yCotLU2TpPT09CGJigNNx3nnnbfllltuWV+9/fbbb1932mmnbZOk+++/v9v27dvrlOvKy0OtA9CM1PTe9l//9V8Zd911V2cpUryuWrWqRUPHUjUXvPjii+1POOGEfg19TAARFFgAklplZaUqKioSHUYoKioqNG7cuN79+vXLO+aYY7K3b99u7733XpvBgwcPyMnJyT3xxBN/VFZWlipFCuKLL76457Bhw/offPDBefPmzUs/6aSTftS7d+9BU6ZMOWhvn/fcc0+nQw45ZOCAAQNyzz333N5NpTgoKSlp2bdv37yzzz67d3Z2dt7Pfvazvn/729/aH3bYYQN69+496M0330yfMWNG5/Hjx/eqvu8ZZ5zRZ+bMmR3/8Ic/dN24cWOLkSNH5hxxxBE5knTeeef1GjRo0MB+/frlXXnlld+9jpmZmYf89re/7TF06ND+N9xwQ/fc3NyBex8rKCholZeXN7D6cRqjrVu3phx//PH9+vfvn5udnZ334IMPdpw3b176kCFDBvTv3z/3kEMOGbh58+YUSVq/fn2L4447Lrt3796DJk6cmLW3j/vvv79TTk5ObnZ2dt6kSZMya2uHdM0115RdfvnlX0nSE0880WX16tUNXmDFYs+ePYkOAWjUKLCAJqi8vFynn356n5ycnNwxY8YcvG3btu/9X696lnXmzJkdzzjjjD6StG7durSTTz75R4MGDRo4aNCgga+88krbmo6xZcuWlF/84hd9cnJycnNycnIfffTRA6W6fcjasmVLylFHHZWTm5s7MCcnJ/eJJ544UIp8uD744IPzzj///F55eXm5K1asaBnDy5E0Vq9e3XrKlCkbly9fvuSAAw6oeOyxxzpeeOGFfW+55Za1S5cuLczLy9t57bXXfvehv2XLlpX5+fklF110UdmZZ57Z78EHH1xdXFy85M9//nOX9evXpy5cuLD17NmzO+Xn5xcXFxcXpqSk+H333dc5kc8xTGvWrGl99dVXbywuLl6yYsWK1k8++WTn/Pz84unTp6+dPn16j9r2v/HGGzd27dp1z7x585bOnz9/qSTddtttXyxevLiouLh4ybvvvtt+/vz5bfZu37p168oFCxaU3Hrrrevbt29f8d5777WRpPvvv7/Lueee+1XDPdP4+ctf/tKhe/fue0pKSgqXLVu25PTTT9963nnn/ej2229fXVJSUjhv3rySdu3aVUpSYWFh+t/+9rfPioqKlrzwwgsdly9f3mLVqlUtbr755sy5c+cuLSwsXPLxxx+3ffzxxw+sqT3RzzdsN954Y7c//OEPXSXp4osv7nnkkUfmSNKcOXPajx07tq8k/frXv87s379/7uDBgwesWbMmTZKuuuqqg6ZNm9Zt5syZHRcvXpw+fvz4gwcMGJC7fft2e/vtt9MPP/zw/nl5eQOPPfbY7M8//7zG4mvx4sWtjj766Jz+/fvn5ubmDlyyZEmryspK/fKXv8zKzs7Oy8nJyX3wwQc77us5vPnmm+lDhgwZMHDgwNwhQ4YMWLRoUSspMrJ2yimnHDxq1Kh+xx13XE5YrxnQHFFgAU3QqlWrWk+cOLFs6dKlhe3bt6/87//+74z92e+Xv/xlz6uuumrD4sWLi/7617+umDhxYp+atr3uuut6dOjQoWLp0qWFS5cuLTz11FO31fVDVnp6euU//vGP5YWFhUXz5s1bev3112dVVlZ+9xwuuuiir4qKigpzcnK+retrkIwyMzN3H3300TslaciQITtWrFjRatu2bamnnnrqdkm69NJLv/rggw/a7d3+5z//+deSNHjw4J39+vXb2bt37z1t2rTxnj177v7ss89avvzyy+0XL16cPnjw4IEDBgzIfeeddzp89tlnrRLz7MKXmZm5e/jw4TtTU1OVk5Ozc9SoUVtTUlJ02GGH7Vi7dm29nuesWbM65ebmDszNzc1dtmxZ60WLFrXe+9j48eM377194YUXfvnggw92KS8v15w5czpefPHFTaLAOuyww3a+/fbbHSZNmpT58ssvt1uxYkXLrl277hk5cuQOSerUqVNlixaRz/fHHnvs1s6dO1ekp6d7v379dq1YsaLVO++80/bII4/cdtBBB5W3aNFCZ5999qZ58+a1q6k9oU+2AZxwwgnb33333XaS9Mknn6R/8803qbt377a33nqr3bHHHrtt586dKUcdddT2kpKSwqOOOmr7nXfe+b333osuumjzoEGDdjz22GOfFRcXF7Zo0UJTpkzpNWfOnBVLliwpmjBhwpe//e1vazwxde655/adOHHixpKSksL8/PziXr167XnssccOLCgoaFNUVLTk9ddfXzpt2rSsfRVpgwcP3vXhhx8WFxUVFd50001fXHPNNd+NTi5cuLDd008/vfKDDz5YGsbrBTRXLHIBNEHdu3f/9qSTTvpGki644IKvZsyY0XV/9nv33Xc7LFu27Lsz+tu3b0/dvHlzSseOHX+wQMBbb73V4Zlnnvls7/2MjIyKf/7zn+33fsiS9N2HrAsuuODraMerrKy0K664IuuDDz5ol5KSoo0bN7Zcu3ZtmiT16NHj29GjR39Tt2ee3Fq2bPndEv6pqan+9ddf73OaUOvWrV2SUlJS1KpVq+/2TUlJUXl5ubm7nXnmmV/dfffdXzRc1IlT9fVKSUn57vVITU1VRUWF1bW/4uLilnfddVe3BQsWFGVkZFScccYZfXbt2vXdicaqC2FMmDBh86233nrQM888s+2QQw7Z0b179yYxT/XQQw/dvXDhwsLnn3/+gBtuuCHzhBNO2GpmUb9aovq/1z179ph79G+hqKm9qTn22GN3TJgwoe3mzZtTWrVq5Yceeuj2t99+O/39999vf+edd65u0aKFjxs3boskDR069JvXXnutw776+/TTT1stW7aszahRo3KkyJTojIyMqPPzNm/enLJhw4aW48eP/1qS0tPTXZK//fbb7c8666xNaWlp6tmzZ/kRRxyx/Z133kkfNmzYzmj9bNq0KfXss8/uu2rVqtZm5nv27Pnu/9Jxxx23tVu3bk3i3zqQSIxgAU2Qme33/Z07d353x92Vn59fVFxcXFhcXFy4cePGT6MVV3u3rd5vXT9k3X///Z2++uqrtIKCgqLi4uLCzp0779m5c2eKFBndqlNnjdABBxxQ0aFDh4qXX365nSQ9/PDDnY866qjt+7v/mDFjtr744osdv/jiizRJ2rBhQ+rSpUubxHTKsLRt27Ziy5YtKZK0efPm1DZt2lR26tSpYs2aNWlz5849oKb90tPTfeTIkVuuuuqqXhdeeOGX8Yu4Ya1atapF+/btK3/1q19tuuKKKzZ89NFHbTds2NBy3rx56VLkQ/y+rr8ZMWLEN/Pnz29fWlqaVl5erueee67T8ccfv72m9rg9sThp1aqVZ2Vl7b777ru7DB8+fPuIESO2v/baa+0///zzVkOGDNmVlpbmKSmRj1ZpaWkqLy/f54kAd7d+/frt3Pueu3Tp0sJ33313WQ3b1tRHnZ7Dtddemzly5Mhty5YtW/L3v/99+bfffvvdZ8Hm8L4LxAMjWEADStSy6qWlpS1fe+21tj/5yU++eeqppzodffTR21966aXvpup17tx5z8KFC1sPHjx415w5czq2a9euQopMCbr11lu7/ud//ucGSXrvvffa7J3SVt3xxx+/9bbbbuv6yCOPrJGksrKy1BEjRnxz7bXX9iwtLU3LyMgof+655zr96le/2lhTnFu2bEnt0qXLnlatWvnf//739uvWrYtPcbAfy6rHy8yZM1dOmjSp95QpU1J69eq1++mnn161v/sOHTp014033vjF6NGjcyorK9WiRQufMWPG6rCnVIa1rHoiTJgw4ctTTjklu2vXrnvmz5+/dNCgQTuys7PzevXqtXvo0KH7LADGjx+/6aWXXup4+umnb22I2MJcVn1/LViwoM3UqVOzUlJSlJaW5vfcc8/n7q4pU6b02rVrV0rr1q0r33rrrRqnh/Xu3XvPtGnTvhg5cmSOu9vo0aO3nH/++V9LUk3tDSVRy6offfTR2+++++5u995776qhQ4fuvP7667MGDRq0Y29hVZt27dpVbNmyJVWSDj300F2bNm1K2/t+vXv3bisoKGg1bNiwXdX369SpU2X37t2/ffzxxw+84IILvt65c6eVl5fbyJEjtz344IMZl19++VcbN25M+/DDD9vNmDFjzd6TVdVt3bo1NSsr61spcn1hDC8FgBrUONwPoO4WLVq0avDgwQk9211SUtLy3/7t37KPOOKIbfn5+e369u27e/bs2StHjx6d/ac//WnNiBEjdsycObPjtGnTMnv06LFnwIABO7/55puU559/flVpaWnaJZdc0mvZsmWtKyoq7Igjjtj21FNPRf0enC1btqRcdNFFvQoKCtqmpKT49ddfv27ChAlf33fffZ1uu+227ns/ZN13331rpe9/F056evqQHTt2fFxaWpp2yimn9CsvL7e8vLwdH330UbuXXnppmSQ11u+MQtMxbdq0blu2bEm944471iU6FiSPOXPmtD/jjDOyN23a9EmHDh0q+/TpM+jCCy8su/nmmzfsfW+TIgsIvfjiiwc8//zzq6p+T9ujjz564M0335zVunXryvz8/KJPP/209ZQpU3pt27YttaKiwiZNmrTh6quvjppHCgoKWl166aW9N23alNaiRQt/7rnnVgwYMODbSZMmZb3xxhsHmJn/x3/8R+mll166uep377344ovt/+d//qfbm2++ufy1115re8kll/Tt1KlT+XHHHbd19uzZnb/44ouC2r77bNGiRV0GDx7cpwFfWqDJoMACQpQMBRaA2J144ok/+vzzz1vNmzdvadhfkAs0RhRYwP5jiiAAANW8+uqrKxIdAwCgcaLAArBPd9xxR+d77723W9W2ww8/fPvjjz8edRoJACA2F1xwQa+PPvroe8vcT5o0acNvfvObJvF1AUBTxxRBIESLFi367JBDDtmckpLCfywAQJNQWVlpBQUFHQcPHnxwomMBGgOWaQfCtbisrOyAysrKOn9HDwAAyaaystLKysoOkLQ40bEAjQVTBIEQlZeXX7J+/fqH1q9fP0icwAAANH6VkhaXl5dfkuhAgMaCKYIAAAAAEBLOsAMAAABASCiwAAAAACAkFFgAAAAAEJL/D+M/C4A+qM1YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "df.groupby(\"Period\").Choice.value_counts().unstack().plot.bar(ax=axs[0], stacked=True, rot=0, legend=False, title=\"Original\")\n",
    "df_3.groupby(\"Period\").Choice.value_counts().unstack().plot.bar(ax=axs[1], stacked=True, rot=0, title=\"Simulated\")\n",
    "\n",
    "handles, _ = axs[1].get_legend_handles_labels()\n",
    "axs[1].get_legend().remove()\n",
    "\n",
    "fig.legend(handles=handles, loc='lower center', bbox_to_anchor=(0.5, 0), ncol=5)\n",
    "\n",
    "plt.tight_layout(rect=[0,0.05,1,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "> Keane, M. P. and Wolpin, K. I. (1997). [The Career Decisions of Young Men](https://doi.org/10.1086/262080). *Journal of Political Economy*, 105(3): 473-522."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}