{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Robinson Crusoe or 10 minutes to respy\n", "\n", "This is a short introduction to respy for new users. As economists love Robinsonades [1](#fn1), we will showcase the implementation of a Robinson Crusoe economy as a discrete choice dynamic programming model. Throughout the notebook you find italic text which make you familiar with Robinson's story. We will first set the scene with a broad introduction and then turn to the precise model specification. We continue by simulating the model and analyze its comparative statics. We then extend the model and showcase the estimation of the model parameters. \n", "\n", "Just to be clear, don't misinterpret the fact that we explain **respy** using such a simplistic model. **respy** is not a toy and can just as well solve state-of-the-art structural models. It's just easier to explain **respy** in a situation where we don't have to explain a complicated model at the same time. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import io\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import respy as rp\n", "import yaml\n", "import seaborn as sns\n", "import numpy as np\n", "\n", "from pathlib import Path\n", "from time import time\n", "\n", "plt.style.use(\"../_static/respy.mplstyle\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction \n", "\n", "---\n", "*After setting sail against his parents' wishes, being captured by pirates, escaping from them, building a plantation, and setting sail again to capture slaves in Africa, [Robinson Crusoe](https://en.wikipedia.org/wiki/Robinson_Crusoe) stranded on a small island. He is alone with one dog, two cats, and only some supplies. He goes fishing to make ends meet and if he is too tired he will relax in his hammock. But, he cannot relax to often as storing food is a difficult task on a tropical island.*\n", "\n", "---\n", "\n", "In the discrete choice dynamic programming model, Robinson chooses every period $t = 0, \\dots, T - 1$ to either go fishing, $a = 0$, or spend the day in the hammock, $a = 1$, to maximize his expected sum of discounted lifetime utility. The utility of a choice, $U(s_t, a_t)$, depends on the state $s_t$, which contains information on the individual's characteristics, and the chosen alternative $a_t$. For working alternatives like fishing utility consists of two components, a wage and a non-pecuniary component.\n", "\n", "$$\n", " U(s_t, a_t) = W(s_t, a_t) + N(s_t, a_t)\n", "$$\n", "\n", "For non-working alternatives like the hammock, $W(s_t, a_t) = 0$. The wage is defined as\n", "\n", "$$\\begin{align}\n", " W(s_t, a_t) &= r_a \\exp\\{x^w_{at} \\beta^w_a + \\epsilon_{at}\\}\\\\\n", " \\ln(W(s_t, a_t)) &= \\ln(r_a) + x^w_{at} \\beta^w_a + \\epsilon_{at}\n", "\\end{align}$$\n", "\n", "\n", "where $r_a$ is normally a market rental price for the skill units generated in the exponential expression. Another interpretation is that $ln(r_a)$ is simply the constant in the skill units. The skill units are generated by two components. $x^w_{at}$ and $\\beta^w_a$ are the choice- and time-dependent covariates and parameters related to the wage signaled by superscript $w$. $\\epsilon_{at}$ is a choice-specific random shock from the shock vector $\\epsilon_t \\sim \\mathcal{N}(0, \\Sigma)$ for all choices. Shocks are usually correlated across choices in one period, but are independent across periods.\n", "\n", "The non-pecuniary rewards for working alternatives are a vector dot product of covariates $x_t^w$ and parameters $\\beta^w$. The superscript $w$ signals that the components belong to working alternatives.\n", "\n", "$$\n", " N^w(s_t, a_t) = x_t^w\\beta^w\n", "$$\n", "\n", "The non-pecuniary reward for non-working alternatives is very similar except that the shocks enter the equation additively. Superscript $n$ stands for non-pecuniary.\n", "\n", "$$\n", " N^n(s_t, a_t) = x_t^n\\beta^n + \\epsilon_{at}\n", "$$\n", "\n", "Along with the lower triangular elements of the shock variance-covariance matrix of $\\epsilon_t$, the utility parameters $\\beta_a^w$, $\\beta_a^n$ and $r_a$ form the main parameters of the model.\n", "\n", "If Robinson chooses to go fishing, he gains one additional unit of experience in the next period. Experience starts at zero and goes over 1, 2, 3 up to $T - 1$.\n", "\n", "The general assumption imposed on Robinson is that he is forward-looking and maximizes the expected present value of utility over the remaining lifetime. W.l.o.g. $t = 0$ and let $V(s_0)$ define the value of the maximization which is achieved by a sequence of choices, $\\{a_t\\}^T_{t = 0}$, such that every action is in the choice set, $a_t \\in C(s_t)$ and $s_{t + 1}$ adheres to the law of motion. Then, the expected present value of lifetime utility in state $s_0$ is\n", "\n", "$$\n", " V(s_0) = \\text{E} \\max_{\\{a_t\\}^T_{t = 0}} \\left[\n", " \\sum^T_{t = 0} \\delta^t U(s_t, a_t) \\, \\Big|\n", " \\, a_t \\in C(s_t), s_{t+1} = m(s_t, a_t)\n", " \\right]\n", "$$\n", "\n", "Note that the shocks in period $t = 0$ are not stochastic. Thus, one can extract the utility in the current period $U(s_0, a_0)$, also called the flow utility, and the discount factor $\\delta$ from the expectation. Then, the formula can be rewritten so that the second term becomes the maximum over alternative-specific value functions at time $t = 1$ which are also called continuation values.\n", "\n", "$$\\begin{align}\n", " V(s_0) &= \\max_{a_0} \\, U(s_0, a_0) + \\delta \\text{E} \\max_{\\{a_t\\}^T_{t = 1}}\n", " \\left[\n", " \\sum^T_{t = 1} \\delta^{t - 1} U(s_t, a_t) \\, \\Big|\n", " \\, a_t \\in C(s_t), s_{t + 1} = m(s_t, a_t)\n", " \\right] \\\\\n", " &= \\max_{a_0} \\, U(s_0, a_0)\n", " + \\delta \\text{E} \\max_{a_1} V_{a_1}(s_1)\n", "\\end{align}$$\n", "\n", "The maximum over alternative-specific value functions can be rewritten as the value function of state $s_1$ or $V(s_1) = \\max_{a_1} V_{a_1}(s_1)$ which yields the famous Bellman equation. Due to the recursive nature of the problem, the alternative-specific value functions are defined as\n", "\n", "$$\\begin{equation}\n", " V_a(s_t) = \\begin{cases}\n", " U(s_t, a_t) + \\delta \\text{E} V(s_{t+1}) & \\text{if } t < T \\\\\n", " U(s_t, a_t) & \\text{if } t = T\n", " \\end{cases}\n", "\\end{equation}$$\n", "\n", "The former equation shows that the shocks in period $t + 1$ are unknown to the individual in period $t$. Thus, utility must be maximized given the joint distribution of shocks in period $t + 1$ which is a maximization problem over a two-dimensional integral. Denote the non-stochastic part of a state as $s^-$. Then, Robinson maximizes\n", "\n", "$$\\begin{equation}\n", " V(s_t) = \\max_{a_t}\\{\n", " U(s_t, a_t) + \\delta \\int_{\\epsilon_{0, t + 1}} \\dots \\int_{\\epsilon_{2, t + 1}}\n", " \\max_{a_{t + 1}} V_{a_{t + 1}}(s^-_{t + 1}, \\epsilon_{t + 1})\n", " f_\\epsilon(\\epsilon_{t + 1})\n", " d_{\\epsilon_{0, t + 1}} \\dots d_{\\epsilon_{2, t + 1}}\n", " \\}\n", "\\end{equation}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Specification\n", "\n", "How can we express the equations and parameters with **respy**? The following cell contains the code to write a `.csv` file which is the cornerstone of a model as it contains all parameters and some other specifications. With `io.StringIO` we can pretend it is an actual file on the filesystem and easily loaded with `pandas`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "params = \"\"\"category,name,value\n", "delta,delta,0.95\n", "wage_fishing,exp_fishing,0.1\n", "nonpec_fishing,constant,-1\n", "nonpec_hammock,constant,2.5\n", "nonpec_hammock,not_fishing_last_period,-1\n", "shocks_sdcorr,sd_fishing,1\n", "shocks_sdcorr,sd_hammock,1\n", "shocks_sdcorr,corr_hammock_fishing,-0.2\n", "lagged_choice_1_hammock,constant,1\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ " value\n", "category name \n", "delta delta 0.95\n", "wage_fishing exp_fishing 0.10\n", "nonpec_fishing constant -1.00\n", "nonpec_hammock constant 2.50\n", " not_fishing_last_period -1.00\n", "shocks_sdcorr sd_fishing 1.00\n", " sd_hammock 1.00\n", " corr_hammock_fishing -0.20\n", "lagged_choice_1_hammock constant 1.00" ], "text/html": "
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
value
categoryname
deltadelta0.95
wage_fishingexp_fishing0.10
nonpec_fishingconstant-1.00
nonpec_hammockconstant2.50
not_fishing_last_period-1.00
shocks_sdcorrsd_fishing1.00
sd_hammock1.00
corr_hammock_fishing-0.20
lagged_choice_1_hammockconstant1.00
\n
" }, "metadata": {}, "execution_count": 3 } ], "source": [ "params = pd.read_csv(io.StringIO(params), index_col=[\"category\", \"name\"])\n", "params" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The DataFrame for the parameters contains a two-level MultiIndex to group parameters in categories. `name` should be uniquely assigned in each category or otherwise only the sum of identically named parameters is identified. `value` contains the value of the parameter. Note that we named Robinson's alternatives `\"fishing\"` and `\"hammock\"` and we have to use the names consistently. As long as you stick to lowercase letters separated by underscores, you can choose any name you want.\n", "\n", "The parameter specification contains following entries:\n", "\n", "- The first entry contains the discount factor of individuals.\n", "- The second category `\"wage_fishing\"` contains the parameters of the log wage equation for fishing. The group contains only one name called `\"exp_fishing\"` where `\"exp_*\"` is an identifier in the model for experience accumulated in a certain alternative. **respy** requires that you respect those identifiers of which there are not many and reference your alternatives consistently with the same name. If you stick to lowercase letters possibly separated by underscores, you are fine.\n", "- The third and fourth categories concern the non-pecuniary reward of fishing and relaxing in the hammock.\n", "- `\"shocks_sdcorr\"` groups the lower triangular of the variance-covariance matrix of shocks.\n", "- `\"lagged_choice_1_hammock\"` governs the distribution of previous choices at the begin of the model horizon.\n", "\n", "`params` is complemented with `options` which contains additional information. Here is short description:\n", "\n", "- `\"n_periods\"` defines the number of periods for which decision rules are computed.\n", "- `\"_seed\"`: Seeds are used in every model component to ensure reproducibility. You can use any seed you would like or even repeat the same seed number. Internally, we ensure that randomness is completely uncorrelated.\n", "- `\"estimation_draws\"` defines the number of draws used to simulate the choice probabilities with Monte Carlo simulation in the maximum likelihood estimation.\n", "- `\"estimation_tau\"` controls the temperature of the softmax function to avoid zero-valued choice probabilities.\n", "- `\"interpolation_points\"` controls how many states are used to approximate the value functions of others states in each period. `-1` turns the approximation off. The approximation is detailed in Keane and Wolpin (1994).\n", "- ``\"simulation_agents\"`` defines how many individuals are simulated.\n", "- ``\"solution_draws\"`` defines the number of draws used to simulate the expected value functions in the solution.\n", "- `\"covariates\"` is another dictionary where the key determines the covariate's name and the value is its definition. Here, we have to define what `\"constant\"` means." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "options = \"\"\"n_periods: 10\n", "estimation_draws: 200\n", "estimation_seed: 500\n", "estimation_tau: 0.001\n", "interpolation_points: -1\n", "simulation_agents: 1_000\n", "simulation_seed: 132\n", "solution_draws: 500\n", "solution_seed: 456\n", "covariates:\n", " constant: \"1\"\n", " not_fishing_last_period: \"lagged_choice_1 != 'fishing'\"\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "{'n_periods': 10,\n", " 'estimation_draws': 200,\n", " 'estimation_seed': 500,\n", " 'estimation_tau': 0.001,\n", " 'interpolation_points': -1,\n", " 'simulation_agents': 1000,\n", " 'simulation_seed': 132,\n", " 'solution_draws': 500,\n", " 'solution_seed': 456,\n", " 'covariates': {'constant': '1',\n", " 'not_fishing_last_period': \"lagged_choice_1 != 'fishing'\"}}" ] }, "metadata": {}, "execution_count": 5 } ], "source": [ "options = yaml.safe_load(options)\n", "options" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "\n", "We are now ready to simulate the model." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "simulate = rp.get_simulate_func(params, options)\n", "df = simulate(params)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ " Experience_Fishing Lagged_Choice_1 Shock_Reward_Fishing \\\n", "Identifier Period \n", "0 0 0 hammock -0.035035 \n", " 1 1 fishing 0.074254 \n", " 2 1 hammock -0.354560 \n", " 3 1 hammock -0.109397 \n", " 4 2 fishing -1.063705 \n", " 5 2 hammock -0.106235 \n", " 6 3 fishing -0.692603 \n", " 7 3 hammock -0.555217 \n", " 8 3 hammock -0.943424 \n", " 9 3 hammock -1.541738 \n", "1 0 0 hammock -0.713137 \n", " 1 0 hammock -0.464134 \n", " 2 0 hammock 0.066755 \n", " 3 1 fishing -0.748565 \n", " 4 1 hammock -0.130574 \n", "\n", " Meas_Error_Wage_Fishing Shock_Reward_Hammock \\\n", "Identifier Period \n", "0 0 1 0.040965 \n", " 1 1 1.506491 \n", " 2 1 1.185316 \n", " 3 1 -0.785877 \n", " 4 1 1.245234 \n", " 5 1 -1.854834 \n", " 6 1 -1.238533 \n", " 7 1 -0.542048 \n", " 8 1 -0.581919 \n", " 9 1 0.487632 \n", "1 0 1 1.418950 \n", " 1 1 -0.384774 \n", " 2 1 -0.844060 \n", " 3 1 -1.348821 \n", " 4 1 -0.539693 \n", "\n", " Meas_Error_Wage_Hammock Dense_Key Core_Index Choice \\\n", "Identifier Period \n", "0 0 1 0 1 fishing \n", " 1 1 1 0 hammock \n", " 2 1 2 0 hammock \n", " 3 1 3 1 fishing \n", " 4 1 4 6 hammock \n", " 5 1 5 1 fishing \n", " 6 1 6 10 hammock \n", " 7 1 7 13 hammock \n", " 8 1 8 1 hammock \n", " 9 1 9 3 hammock \n", "1 0 1 0 1 hammock \n", " 1 1 1 3 hammock \n", " 2 1 2 1 fishing \n", " 3 1 3 5 hammock \n", " 4 1 4 9 fishing \n", "\n", " Wage ... Nonpecuniary_Reward_Fishing Wage_Fishing \\\n", "Identifier Period ... \n", "0 0 0.965571 ... -1 0.965571 \n", " 1 NaN ... -1 1.190358 \n", " 2 NaN ... -1 0.775258 \n", " 3 0.990647 ... -1 0.990647 \n", " 4 NaN ... -1 0.421597 \n", " 5 1.098301 ... -1 1.098301 \n", " 6 NaN ... -1 0.675297 \n", " 7 NaN ... -1 0.774749 \n", " 8 NaN ... -1 0.525490 \n", " 9 NaN ... -1 0.288882 \n", "1 0 NaN ... -1 0.490104 \n", " 1 NaN ... -1 0.628679 \n", " 2 1.069034 ... -1 1.069034 \n", " 3 NaN ... -1 0.522795 \n", " 4 0.969889 ... -1 0.969889 \n", "\n", " Flow_Utility_Fishing Value_Function_Fishing \\\n", "Identifier Period \n", "0 0 -0.034429 18.430145 \n", " 1 0.190358 17.891912 \n", " 2 -0.224742 15.433115 \n", " 3 -0.009353 13.645001 \n", " 4 -0.578403 11.723370 \n", " 5 0.098301 10.090155 \n", " 6 -0.324703 7.806439 \n", " 7 -0.225251 5.408186 \n", " 8 -0.474510 2.664255 \n", " 9 -0.711118 -0.711118 \n", "1 0 -0.509896 17.954678 \n", " 1 -0.371321 16.268399 \n", " 2 0.069034 14.828094 \n", " 3 -0.477205 13.177150 \n", " 4 -0.030111 11.562961 \n", "\n", " Continuation_Value_Fishing Nonpecuniary_Reward_Hammock \\\n", "Identifier Period \n", "0 0 19.436393 1.5 \n", " 1 18.633215 2.5 \n", " 2 16.481955 1.5 \n", " 3 14.373005 1.5 \n", " 4 12.949235 2.5 \n", " 5 10.517741 1.5 \n", " 6 8.559098 2.5 \n", " 7 5.929934 1.5 \n", " 8 3.303963 1.5 \n", " 9 0.000000 1.5 \n", "1 0 19.436393 1.5 \n", " 1 17.515494 1.5 \n", " 2 15.535854 1.5 \n", " 3 14.373005 2.5 \n", " 4 12.203234 1.5 \n", "\n", " Wage_Hammock Flow_Utility_Hammock Value_Function_Hammock \\\n", "Identifier Period \n", "0 0 NaN 1.547145 18.414638 \n", " 1 NaN 3.961203 20.072216 \n", " 2 NaN 2.732280 16.936717 \n", " 3 NaN 0.751880 13.073149 \n", " 4 NaN 3.932816 14.967629 \n", " 5 NaN -0.296111 8.560856 \n", " 6 NaN 1.425011 8.519145 \n", " 7 NaN 1.079947 5.832661 \n", " 8 NaN 1.118523 3.427303 \n", " 9 NaN 2.286128 2.286128 \n", "1 0 NaN 3.032909 19.900402 \n", " 1 NaN 1.215827 16.377564 \n", " 2 NaN 0.659643 14.053384 \n", " 3 NaN 1.328144 13.649412 \n", " 4 NaN 0.997326 11.386079 \n", "\n", " Continuation_Value_Hammock \n", "Identifier Period \n", "0 0 17.755256 \n", " 1 16.958962 \n", " 2 14.952039 \n", " 3 12.969757 \n", " 4 11.615592 \n", " 5 9.323124 \n", " 6 7.467509 \n", " 7 5.002857 \n", " 8 2.430295 \n", " 9 0.000000 \n", "1 0 17.755256 \n", " 1 15.959723 \n", " 2 14.098675 \n", " 3 12.969757 \n", " 4 10.935529 \n", "\n", "[15 rows x 22 columns]" ], "text/html": "
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
Experience_FishingLagged_Choice_1Shock_Reward_FishingMeas_Error_Wage_FishingShock_Reward_HammockMeas_Error_Wage_HammockDense_KeyCore_IndexChoiceWage...Nonpecuniary_Reward_FishingWage_FishingFlow_Utility_FishingValue_Function_FishingContinuation_Value_FishingNonpecuniary_Reward_HammockWage_HammockFlow_Utility_HammockValue_Function_HammockContinuation_Value_Hammock
IdentifierPeriod
000hammock-0.03503510.040965101fishing0.965571...-10.965571-0.03442918.43014519.4363931.5NaN1.54714518.41463817.755256
11fishing0.07425411.506491110hammockNaN...-11.1903580.19035817.89191218.6332152.5NaN3.96120320.07221616.958962
21hammock-0.35456011.185316120hammockNaN...-10.775258-0.22474215.43311516.4819551.5NaN2.73228016.93671714.952039
31hammock-0.1093971-0.785877131fishing0.990647...-10.990647-0.00935313.64500114.3730051.5NaN0.75188013.07314912.969757
42fishing-1.06370511.245234146hammockNaN...-10.421597-0.57840311.72337012.9492352.5NaN3.93281614.96762911.615592
52hammock-0.1062351-1.854834151fishing1.098301...-11.0983010.09830110.09015510.5177411.5NaN-0.2961118.5608569.323124
63fishing-0.6926031-1.2385331610hammockNaN...-10.675297-0.3247037.8064398.5590982.5NaN1.4250118.5191457.467509
73hammock-0.5552171-0.5420481713hammockNaN...-10.774749-0.2252515.4081865.9299341.5NaN1.0799475.8326615.002857
83hammock-0.9434241-0.581919181hammockNaN...-10.525490-0.4745102.6642553.3039631.5NaN1.1185233.4273032.430295
93hammock-1.54173810.487632193hammockNaN...-10.288882-0.711118-0.7111180.0000001.5NaN2.2861282.2861280.000000
100hammock-0.71313711.418950101hammockNaN...-10.490104-0.50989617.95467819.4363931.5NaN3.03290919.90040217.755256
10hammock-0.4641341-0.384774113hammockNaN...-10.628679-0.37132116.26839917.5154941.5NaN1.21582716.37756415.959723
20hammock0.0667551-0.844060121fishing1.069034...-11.0690340.06903414.82809415.5358541.5NaN0.65964314.05338414.098675
31fishing-0.7485651-1.348821135hammockNaN...-10.522795-0.47720513.17715014.3730052.5NaN1.32814413.64941212.969757
41hammock-0.1305741-0.539693149fishing0.969889...-10.969889-0.03011111.56296112.2032341.5NaN0.99732611.38607910.935529
\n

15 rows × 22 columns

\n
" }, "metadata": {}, "execution_count": 7 } ], "source": [ "df.head(15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can inspect Robinson's decisions period by period." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2021-06-01T13:22:43.769935\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.2, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEVCAYAAACv2pHlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXhklEQVR4nO3de7SddX3n8fcHglwSglzSWJUQFXAQp9AxXTq1KN6qTrXSYXRZkaJTpQvUNVqto44MCFgvdeyMinSYEVERvIyoiGK9oFNwrBrbiTMZayoaLkIg3EJOuGn6nT+eJ7o5nOTsQ3bO/mXn/Vprr5z97N9+9vfsc04++/c8v+f3S1UhSdK47TbuAiRJAgNJktQIA0mS1AQDSZLUBANJktQEA0mS1IQF4y4A4KCDDqrly5ePuwxJ0g72/e9//5aqWjLTY00E0vLly1m5cuW4y5Ak7WBJrtnaYx6ykyQ1wUCSJDXBQJIkNcFAkiQ1YahASvLqJCuT3Jvkglnavi7JuiQbkpyfZM+RVCpJmmjD9pBuAM4Gzt9WoyTPBt4EPANYDjwaeNt21CdJ2kUMFUhVdUlVfQ64dZamJwEfqqrVVXU7cBbwsu2qUJK0Sxj1OaQjgVUD91cBS5McOOLXkSRNmFFfGLsI2DBwf8vX+zKtd5XkZOBkgGXLlg239zP22+4Cf7WvDbO3GXpf1jW3fVnX3PZlXXPbl3XNbV/t1DXqHtIUsHjg/pavN05vWFXnVdWKqlqxZMmMs0hIknYhow6k1cBRA/ePAm6qqtnOPUmSdnHDDvtekGQvYHdg9yR7JZnpcN9HgT9O8rgk+wNvBS4YWbWSpIk1bA/prcDddEO6X9p//dYky5JMJVkGUFVfBt4NfAO4pr+dPvKqJUkTZ6hBDVV1BnDGVh5eNK3te4H3bldVkqRdjlMHSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmjBUICU5IMlnk2xKck2Sl2ylXZKcneRnSTYk+WaSI0dbsiRpEg3bQzoHuA9YCpwAnLuVoHkh8G+BY4ADgG8DHxtBnZKkCTdrICVZCBwPnFZVU1V1FXApcOIMzR8FXFVVP6mqzcCFwONGWbAkaTIN00M6HNhcVWsGtq0CZuohfQI4NMnhSfYATgK+vP1lSpIm3YIh2iwCNkzbtgHYd4a2NwJXAj8CNgPXAU+faadJTgZOBli2bNmQ5UqSJtUwPaQpYPG0bYuBjTO0PR34LeBgYC/gbcAVSfaZ3rCqzquqFVW1YsmSJXOrWpI0cYYJpDXAgiSHDWw7Clg9Q9ujgE9W1fVV9YuqugDYH88jSZJmMWsgVdUm4BLgzCQLkzwZeAEzj577HvDCJEuT7JbkRGAP4MejLFqSNHmGOYcEcCpwPnAzcCtwSlWtTrIM+H/A46rqWuBdwK8B/xtYSBdEx1fVHSOuW5I0YYYKpKq6DThuhu3X0g162HL/HuBV/U2SpKE5dZAkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQlDBVKSA5J8NsmmJNckeck22j46yWVJNia5Jcm7R1euJGlSDdtDOge4D1gKnACcm+TI6Y2SPAT4KnAF8DDgkcCFoylVkjTJZg2kJAuB44HTqmqqqq4CLgVOnKH5y4Abquq9VbWpqu6pqh+MtGJJ0kQapod0OLC5qtYMbFsFPKCHBDwJWJvk8v5w3TeT/POZdprk5CQrk6xcv3793CuXJE2UYQJpEbBh2rYNwL4ztH0k8GLgfcDDgS8Cn+8P5d1PVZ1XVSuqasWSJUvmVrUkaeIME0hTwOJp2xYDG2doezdwVVVdXlX3Ae8BDgSO2K4qJUkTb5hAWgMsSHLYwLajgNUztP0BUKMoTJK0a5k1kKpqE3AJcGaShUmeDLwA+NgMzS8EnpTkmUl2B14L3AL8cHQlS5Im0bDDvk8F9gZuBi4GTqmq1UmWJZlKsgygqn4EvBT4K+B2uuD6/f7wnSRJW7VgmEZVdRtw3Azbr6Ub9DC47RK6HpUkSUMbKpBasfyei0a2r7Uj25MkaRScy06S1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1IQF4y5gEiy/56KR7WvtyPYkSTsXe0iSpCbYQ5pg9twk7UzsIUmSmmAgSZKa4CE7zbtWDyW2Wpe0qzCQpMYZlNpVGEiSNA/8YDE7zyFJkppgIEmSmjDUIbskBwAfAn4XuAV4c1Vts/+Z5ArgacAeVfWL7S1UUltaPQTVal2a3bDnkM4B7gOWAkcDX0yyqqpWz9Q4yQlz2LckSbMfskuyEDgeOK2qpqrqKuBS4MSttN8POB144ygLlSRNtmHOIR0ObK6qNQPbVgFHbqX9nwPnAuu2szZJ0i5kmEBaBGyYtm0DsO/0hklWAE8G3j/bTpOcnGRlkpXr168fplZJ0gQbJpCmgMXTti0GNg5uSLIb8EHg3w0ziKGqzquqFVW1YsmSJcPWK0maUMME0hpgQZLDBrYdBUwf0LAYWAF8Msk64Hv99uuTHLPdlUqSJtqsI+GqalOSS4Azk7yCbpTdC4DfntZ0A/DwgfsHA98FngB4TE6StE3DXhh7KrA3cDNwMXBKVa1OsizJVJJl1Vm35cavQuimqrpvB9QuSZogQ10rVFW3AcfNsP1aukEPMz1nLZDtqE2StAtx6iBJUhMMJElSEwwkSVITDCRJUhMMJElSEwwkSVITDCRJUhMMJElSEwwkSVITDCRJUhMMJElSEwwkSVITDCRJUhOGmu1bkjSZlt9z0cj2tXY7n28PSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUhKECKckBST6bZFOSa5K8ZCvtTkry/SR3Jrk+ybuTuAigJGlWw/aQzgHuA5YCJwDnJjlyhnb7AK8FDgKeCDwDeMP2lylJmnSz9l6SLASOBx5fVVPAVUkuBU4E3jTYtqrOHbj7syQfB542wnolSRNqmB7S4cDmqlozsG0VMFMPabqnAKsfTGGSpF3LMIG0CNgwbdsGYN9tPSnJy4EVwHu28vjJSVYmWbl+/fphapUkTbBhAmkKWDxt22Jg49aekOQ44J3Ac6vqlpnaVNV5VbWiqlYsWbJkyHIlSZNqmEBaAyxIctjAtqPYyqG4JM8B/hvw/Kr6P9tfoiRpVzBrIFXVJuAS4MwkC5M8GXgB8LHpbZM8Hfg4cHxVfXfUxUqSJteww75PBfYGbgYuBk6pqtVJliWZSrKsb3casB/wpX77VJLLR1+2JGnSDHXRalXdBhw3w/Zr6QY9bLnvEG9J0oPi1EGSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmGEiSpCYYSJKkJhhIkqQmDBVISQ5I8tkkm5Jck+Ql22j7uiTrkmxIcn6SPUdXriRpUg3bQzoHuA9YCpwAnJvkyOmNkjwbeBPwDGA58GjgbSOpVJI00WYNpCQLgeOB06pqqqquAi4FTpyh+UnAh6pqdVXdDpwFvGyE9UqSJtQwPaTDgc1VtWZg2yrgAT2kftuqae2WJjnwwZcoSdoVpKq23SA5Bvh0VT1sYNsrgROq6thpba8GXlVVX+7v70F3qO9RVbV2WtuTgZP7u48FfrRd38mvHATcMqJ9jZJ1zY11zY11zY11zc0o6zqkqpbM9MCCIZ48BSyetm0xsHGItlu+fkDbqjoPOG+I15+TJCurasWo97u9rGturGturGturGtu5quuYQ7ZrQEWJDlsYNtRwOoZ2q7uHxtsd1NV3frgS5Qk7QpmDaSq2gRcApyZZGGSJwMvAD42Q/OPAn+c5HFJ9gfeClwwwnolSRNq2GHfpwJ7AzcDFwOnVNXqJMuSTCVZBtCfO3o38A3gmv52+ujL3qaRHwYcEeuaG+uaG+uaG+uam3mpa9ZBDZIkzQenDpIkNcFAkiQ1YWICaS7z7c2nJK9OsjLJvUkuGHc9AEn2TPKh/n3amOTvkzx33HUBJLkwyY1J7kyyJskrxl3ToCSHJbknyYXjrgUgyTf7eqb626iu59tuSV6c5If93+TV/TWN46xnatptc5L3j7OmLZIsT/KlJLf3c4F+IMkwl+Xs6LqOSHJFPzfpj5P8wY58vYkJJIacb28MbgDOBs4fdyEDFgDXAU8F9gNOAz6VZPk4i+q9A1heVYuB3wfOTvKEMdc06Bzge+MuYppXV9Wi/vbYcRcDkORZwLuAlwP7Ak8BfjLOmgbeo0V0/0/cDXx6nDUN+CDdoLFfB46m+9s8dZwF9YH4eeAy4AC6iQwuTHL4jnrNiQikOc63N6+q6pKq+hzQzLVYVbWpqs6oqrVV9U9VdRnwU2Ds//H38yDeu+Vuf3vMGEv6pSQvBu4Avj7mUnYGbwPOrKq/7X/HflZVPxt3UQP+DV0AXDnuQnqPAj5VVfdU1Trgy8w8Pdt8+mfAw4G/rKrNVXUF8C124P+rExFIzG2+PU2TZCndezjTxc7zLskHk9wF/ANwI/ClMZdEksXAmcDrx13LDN6R5JYk30py7LiLSbI7sAJY0h/mub4/BLX3uGsbcBLw0WpnmPF/AV6cZJ8kjwCeSxdK45StbHv8jnrBSQmkRcCGads20B0q0Db08w1+HPhIVf3DuOsBqKpT6X52x9BdlH3vtp8xL86im8n+unEXMs2/p1vm5RF014p8Icm4e5RLgT3oeiHH0B2C+k26C+XHrr9u8qnAR8Zdy4D/SfcB+k7gemAl8LlxFkT3gfBm4M+S7JHkd+net3121AtOSiDNZb499ZLsRjfjxn3Aq8dczv30hwiuAh4JnDLOWpIcDTwT+Mtx1jGTqvpOVW2sqnur6iN0h1T+1ZjLurv/9/1VdWNV3QK8l/HXtcUfAVdV1U/HXQj88u/wr+k+fC2km8h0f7pzcGNTVT8HjgN+D1hHd3TgU3SBuUNMSiDNZb49AUkCfIju0+zx/S9fixYw/nNIx9ItOHltknXAG4Djk/zdOIvaimLmQy3zV0C3Ftr1fS0t+iPa6h0dABwMfKD/YHEr8GEaCPCq+kFVPbWqDqyqZ9P1xr+7o15vIgJpjvPtzaskC5LsBewO7J5krxaGcwLnAkcAz6+qu2drPB+S/Fo/VHhRkt37FYj/ELhizKWdRxeKR/e3vwK+CDx7fCVBkocmefaW36kkJ9CNZvvrcdbV+zDwmv5nuj/wWrrRWmOV5LfpDm+2MrqOvgf5U+CU/uf4ULpzXKu2+cR5kOQ3+t+vfZK8gW4U4AU76vUmIpB6M863N96SgO64+d10S7u/tP96rMfSkxwC/Andf67rBq7LOGGcddF9oj6F7tP17cB7gNdW1efHWlTVXVW1bsuN7hDxPVW1fpx10Z2nORtYT7dWzWuA46qqhWuRzqIbHr8G+CHw98Dbx1pR5yTgkqpq7XD+vwaeQ/ez/DHwC+B1Y62ocyLdwKKbgWcAzxoYBTtyzmUnSWrCJPWQJEk7MQNJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSRqjJJcnOelBPndtkmeOuiZpXFqYMUDa6SRZSzft0mZgE92M5K+pqqm57KeqmlgYUWqBPSTpwXt+v9jbvwB+iznMwJGOf3/SAP8gpO3ULzx3OfD4JE9K8r+S3JFk1eD6RP1y429P8i3gLuDR/bZX9I/vluSt/dLyNyf5aJL9Bp5/Yv/YrUn+w/x+l9KOZyBJ2ynJwXQzM99IN+nq2XQzOL8B+EySJQPNT6RbCnpf4Jppu3pZf3sa3azKi4AP9K/xOLoJcU+kW8XzQLqlOaSJYSBJD97nktwBXEW3wNr1wJeq6kv9st1fpVtobXAZgQv6Zdp/McOSHycA762qn/Tnot5Mt4roArrF7i6rqr/pJ7c8DfinHfvtSfPLQJIevOOq6qFVdUi/yu1S4IX94bo7+rD6Hbop+7fY1oqzD+f+vaZr6AYeLe0f++Vz+yVXbh3NtyG1wVF20uhcB3ysql65jTbbml7/BuCQgfvL6JYhuInucOARWx5Isg/dYTtpYthDkkbnQuD5/aJ5WxZjPDbJsOd6LgZel+RRSRYBfw58sqp+AfwP4HlJfifJQ4Az8e9XE8ZfaGlEquo6upWK30K30Np1wJ8x/N/Z+XSrHP8N3Qqi99Atuke/2OSrgIvoektblgmXJoYL9KkZSR4LfAI4FFgInF5VZ83ynLXAK6rqazM8dgzw36vqsTug3Im0rfdzZ5PkDODQqnrpuGvRcDyHpJa8EfhmVf3mKHZWVVcChpG0kzCQdnHL3/TFHdpFXvvO38scmh9C10Pa9Zyx3449VHHGhrn8HKSx8BySmpDkCroLQj+QZCrJRUnO7h87KMll/VDq25JcOW3anaOT/CDJhiSfTLJX/7xjk1w/8Bprk7xhprb9429McmOSG5K8IkklOXSe3oKWPOD9TLJ//zNYn+T2/utfDtboZ5w4u5+lYirJF5IcmOTjSe5M8r0kywfaV5JTk/xjko1JzkrymCTf7tt/qh+8saX9K5P8uP/5X5rk4QOPHZnkq/1jNyV5y/RvKMkeSS5O8pnB/aotBpKaUFVPB64EXt3PD3ffwMOvpzuBv4Tumpy3cP/h0y8CngM8CvgNutkOtmbGtkmeA/wp8Ey6c1hP3c5vaWc203u0G/Bhul7sMuBu+lkkBryYbiaJRwCPAb7dP+cA4IfA6dPaPwd4AvAkusO159FdHHww8HjgDwGSPB14R1/Xr9Ndn/WJ/rF9ga8BX6a7VutQ4OuDL5Jkb+BzwL3Ai6pq8HdLDTGQtDP4Od1/RIdU1c+r6sq6/2ic91XVDVV1G/AF4Oht7GtrbV8EfLifReEu4G0j/y52Hg94j6rq1qr6TFXdVVUbgbfzwND+cFVdXVUb6Ob2u7qqvtYPW/80MP3c4Luq6s5+BOH/Bb7Sz1Kx5flb2p8AnF9Vf9fPUvFm4F/2Pa7nAeuq6j9V1T1VtbGqvjPwGovpwupq4OVVtXkk75B2CANJO4O/AH4MfCXJT5K8adrj6wa+votuDrit2Vrb+82EwLZnVJh0D3iPkuyT5L/2k7veSTc0/aFJdh9oe9PA13fPcH/6z2XY9vebwaKfVulWup7YwXRhszVPouvlvXPahxg1yEBS8/pPva+vqkcDzwf+NMkzRvwyN3L/yUoPHvH+d3avpxux+MSqWgw8pd8+H4Ml7jeDRZKFdLNU/Izug8NjtvHcr9Ad7vt6kqU7skhtPwNJzUvyvCSHJglwJ92ieKM+9PIp4OVJjuin5fmPI97/zm5ful7LHUkO4IHng3aki+h+Nkcn2ZNuBovvVNVa4DLgYUlem2TPJPsmeeLgk6vq3f0+vp7koHmsW3NkIGlncBjdiespuhPlH6yqb47yBarqcuB9wDfoDg9+u3/o3lG+zk7sPwN7A7cAf0t3XmZeVNXX6WY3/wxdT/YxdAMo6M9nPYuu57wO+Ee60ZrT93EW3cCGr/WBqgY5U4M0gyRH0J1o37M/KS9pB7OHJPWS/EGShyTZH3gX8AXDSJo/BpL0K39CNynq1XTnqE4ZbznSrsVDdpKkJthDkiQ1wUCSJDXBQJIkNcFAkiQ1wUCSJDXBQJIkNcFAkiQ1wUCSJDXBQJIkNcFAkiQ1wUCSJDXh/wOIybW13LLfbwAAAABJRU5ErkJggg==\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "fig, ax = plt.subplots()\n", "\n", "df.groupby(\"Period\").Choice.value_counts(normalize=True).unstack().plot.bar(\n", " stacked=True, ax=ax\n", ")\n", "\n", "plt.xticks(rotation=\"horizontal\")\n", "\n", "plt.legend(loc=\"lower center\", bbox_to_anchor=(0.5,-0.275), ncol=2)\n", "\n", "plt.show()\n", "plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also analyze the persistence in decisions." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "" ] }, "metadata": {}, "execution_count": 9 }, { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2021-06-01T13:22:44.276396\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.2, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEUCAYAAAD5i0vIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtFUlEQVR4nO3dd5xcVf3/8dd7U0gHQglJSIOEjgSkhk5AQ49GqdK+AlLkBwoiAkGafAEbIEWDgIKooBQDoQqKIAgCEvmGQAgllRASSNlNTz6/P+7dMJndnZ1JZiazs+8nj/vIzLnnnnNmd9nPnnvOPUcRgZmZWSnUrO0GmJlZ9XKQMTOzknGQMTOzknGQMTOzknGQMTOzknGQMTOzkmm7thtQqEXL8JxrK4q+37p/bTfBqsjMO45SscrquOO38/49t/A/Nxet3lJocUHGzKzqqXpuMjnImJlVGlV056QgDjJmZpXGPRkzMysZ92TMzKxkatqs7RYUjYOMmVml8e0yMzMrGd8uMzOzknFPxszMSsY9GTMzKxn3ZMzMrGTckzEzs5JxT8bMzEqmjZ+TMTOzUnFPxszMSsZjMmZmVjLuyZiZWcm4J2NmZiXjnoyZmZWMezJmZlYyXurfzMxKxrfLzMysZKrodln1hEszs2qhmvyPfIqTukt6SFKdpEmSjmsi30mSXpM0T9JUSddLaptxvr+kxyR9JmmGpJszzzfGQcbMrNIUOcgAtwBLgB7A8cBtkrZtJF8n4DxgQ2A3YChwQcb5W4GZQE9gMLAvcFauin27zMys0hTxdpmkzsAIYLuIqAVekDQaOAG4KDNvRNyW8XaapHuB/TPSBgA3R8QiYIakJ4DGgtVK7smYmVWaAnoykk6X9GrGcXpWaVsAyyNiQkbaWJoJDql9gHEZ728EjpHUSVJv4GDgiVwFuCdjZlZpCpjCHBGjgFE5snQB5malzQW65ipX0inAzsCpGcnPAacB84A2wG+Bh3OV456MmVmlkfI/mlcLdMtK6wbMb7p6DQeuBQ6OiFlpWg3wJPAg0Jlk3GZ94LpclTvImJlVGEl5H3mYALSVNCgjbQdWvQ2WWfcw4Hbg8Ih4M+NUd6APyZjM4oiYDdwFHJKrcgcZM7MKU8wgExF1JL2PKyV1lrQncCRwTyP1HgDcC4yIiFeyypkFfACcKamtpPWAk0jGd5rkIGNmVmlUwJGfs4COJNOP/wCcGRHjJPWVVCupb5pvJLAu8FiaXivp8YxyvgoMAz4BJgLLgO/kqtgD/2ZmFSbP22B5i4hPgeGNpE8mmRhQ/37/7DxZ+d8A9iukbgcZM7MKU+wgszY5yJiZVRgHGTMzKxnVOMiYmVmJuCdjZmYl4yBjZmYl4yBjZmYl4yBjZmalUz0xxkHGzKzSuCdjZmYlU1NTPSt+OciYmVWa6unIOMiYmVUa3y4zM7OScZAxM7OScZAxM7OScZAxM7PSqZ4Y4yBjZlZpPIXZzMxKpppul1VPuKxyc+fM4bz/dza77TyYYQfuz2OPPtJovnffncAZp32TfffcjR223bLJ8iZN+pBddtyeH3z/glI12SrYep3b85uzh/DBrV/ltesP5au79W0039FD+vH0yAN57+av8MaPD+Oyr32BNhl7nXxwy1dWOT66/Wtcc9yO5foY1UsFHPkUJ3WX9JCkOkmTJB3XRL6TJL0maZ6kqZKul9Q2K88xksanZb0nae9cdbsn00Jcc/WVtGvXjr8990/efns855z1LbbYaisGDhy0Sr52bdvypWHDOPrYYznvnLNzlrftdtuXutlWoa49fieWLFvBdt8ZzXZ91uPec/di3JQ5vDN93ir5OrZvy6V/fIPX3/+UDbquwz3n7MlZX96SXzz+NgADzn5oZd5O7dsw7udHMPrVKWX9LNWoBD2ZW4AlQA9gMDBG0tiIGJeVrxNwHvAysBEwGrgAuDZt10HAdcDRwCtAz+Yqdk+mBViwYAF/ffopzj7nXDp17sxOX9yZffc/gEdH/6VB3v4DNuOrI77O5psPaqSkxOOPjaFb167stvsepWy2VahO7dtw2Bd7c+3D/0fd4mW8PHEWT46dztf36Ncg72/+/h4vvzuLpctXMGPOQh54eTK7Dtyw0XIP33lTZs1fzL8mzCr1R6h6kvI+8iirMzACGBkRtRHxAknwOCE7b0TcFhHPR8SSiJgG3AvsmZHlCuDKiPhXRKyIiGlpviaVLchI2qyJo7ckB7scJk36kDZtaujff8DKtC233Ir3Jk4suKza2lpuvfkmzv/eRcVsorUgm23SleUrgvc/rl2ZNm7KXLbstW6z1+6xxUa8M31uo+eOGtKf+1+cVLR2tmbFDDLAFsDyiJiQkTYW2DaPa/cBxqVtagPsDGwkaWJ6O+1mSR1zFVDOX+4TgXfTI/P1ZGCxpAck9Shje1qMhQsW0KVL11XSunTpyoIFdQWXdcsvbuArXx3BJj2b7eValeq8TlvmL1y6Stq8hUvp0iH33fNj9uzPDv3W59Yn32lwrnf3TgzZciPue/HDYja11SokyEg6XdKrGcfpWcV1AbL/MpgLdCUHSaeQBJWfpEk9gHbA14C9SW677QhcmquccgaZ00i6XlsAHYAtgd8BZwHbk4wP3dLYhZlfxDK1taJ07NSJurraVdJq62rp1KlzQeW8PX48/3rpJU448eQits5amrrFy+jSod0qaV07tKV20bImrzl4x16MHPEFjr3heT6tXdLg/FFD+vHyu7OYPKvwP3ysIdUo7yMiRkXEzhnHqKziaoFuWWndgPlN1i8NJxmHOTgi6u9/Lkz//UVEfJSm/ww4JNdnKefA/xXAwIhYlL6fKOlMYEJE/ErSySQ9mwbSL9oogEXLiHI0tpL069efZcuWM2nSh/Tr1x+ACe+8zeYDBxZUzqv/fpnp06fx5QP3B5KxnhUrlnP0e1/hvj8/1MzVVi3enzGftm3EgI278MHM5I+Xbfus1+RtsP2324SfnrQzx9/4AuOnNXGrbI/+3JROBrA1V+SB/wlAW0mDIqL+d+wOpLfBGql7GHA7cGhEvFmfHhGfSZoKhf0OLmdPpgbon5XWF2iTvq7Fs90a1alTJ4YedBC3/uImFixYwH9ef42/P/sMhx1xZIO8EcHixYtZujS5HbJ48WKWLEn+8hzx9aMZ8/jT3P/Aw9z/wMN8/ehj2Huf/bht1B1l/Ty2di1Yspwxr0/j+8O3o1P7Nuw6cAOGDe7Fn15qOJ6y11Ybc9tpu/E/t77Ifz74tNHydtl8AzZZv6NnlRWRlP/RnIioAx4ErpTUWdKewJHAPQ3r1QEkd5xGRMQrjRR3F3COpI0lrU8yE+3RXPWX85f6DcCzku4CpgCbAqek6QCHAi+VsT0tyiWX/pAfjryY/fcZwnrrrsclIy9n4MBBfDR9Ol854lAeGj2Gnr16MX36NA750tCV1+260xfo1as3jz/9LB07dqRjx8/H6Dp16kT7ddrTvXv3tfGRbC36/u9e58ZTdmHcDUfyWe1iLvzd67wzfR69u3fihau+zF4jn2Tapwv47uHb0K1jO/5w7uePQvzr3Vkce8PzK98fPaQ/j70+lboct9usMCWYwnwWcCcwE5gNnBkR4yT1Bd4CtomIycBIYF3gsYw2PB8RB6evrwI2JOkdLQLuB36Uq2JFlO/uU9oN+zrQC/gIuD8iniikjNZ4u8xKo++37l/bTbAqMvOOo4oWGba48Im8f89NuH5YRS8PUNbbU2lAKSiomJm1NtW0rEzZgoyk9sDJJNPeumSei4gTy9UOM7NKV0Uxpqw9md+SzGh4BPi4jPWambUoNTXVE2XWKMikT+p/IyLuziP7MGBARMxZkzrNzKpdNQWZNZ3C3I5kSls+JgPrrGF9ZmZVr5hTmNe2Znsyki7LcbpdjnPZ7gb+IulGsm6XRcSzBZRjZlbVWtvA/0iSh21qGzlXSE/o2+m/12SlB7BZAeWYmVW11hZkxgO/jIgns09I6gAcm09FETGg+VxmZlZFMSavIPMwsHET55aRzBozM7MiaVU9mYhockwmIpaRLA3TKEnjI2Lr9PUUmlhYLSIa3/vVzKwVqqIYU9znZCTNi4jMJaVPy3j9jWLWZWZWrappCnOxH8Zc5SuTbvNZ//q5ItdlZlaVWtXtsgI1uaibl5UxM8tPFcUYLytjZlZp3JNZPV5WxswsD1UUY0o7JpPFy8qYmeWhVfdkJPUBekfEvxo5fXBW3gMy3npZGTOzPFRRjMk/yKTbdP6BZOA+gC6SvgYMi4hTYdXZZKnGNo/3sjJmZjm01inMvwLGAHuT7BEN8DTw06Yu8FIyZmaFq6bbZYUscLkrcG1ErCCdqhwRc4F1V6diSftL2nt1rjUzq2aS8j7yLK+7pIck1UmaJOm4JvKdJOk1SfMkTZV0vaQGnRFJgyQtkvS75uouJMh8DAzMqmgbkgH9Zkl6TtKe6evvA38E/ijp4gLaYGZW9Uqwn8wtwBKgB3A8cJukbRvJ1wk4D9gQ2A0YClzQRHn/zqfiQoLMT4BHJZ0CtJV0LHAfcF2e128H1E8WOA3YD9gdOKOANpiZVb1i9mQkdQZGACMjojYdOx8NnJCdNyJui4jnI2JJREwD7gX2zCrvGGAO8Ew+nyXvMZmIuFPSp8DpwBTgxLTRD+dZRA0QkjYHFBHj0wavn28bzMxag0KGZCSdTvJ7ud6oiBiV8X4LYHlETMhIGwvsm0fx+wDjMurqBlxJ0sP5Zj7tK2gKcxpQHi7kmgwvADcDPYGHANKAM2s1yzMzq0qFDPynAWVUjixdgLlZaXOBrs204RRgZ+DUjOSrgDsiYkq+bcz7dpmkmyQNyUobIumGPIs4maSL9V/g8jRtK+DGfNtgZtYaFHlMphbolpXWDZjfdP0aDlwLHBwRs9K0wcCBwM8L+SyF9GSOpeEA0GskPZvzmrs4ImYDF2eljSmgfjOzVqFNcZ+TmUAyjj4oIt5N03Yg4zZYJknDgNuBQyPizYxT+wH9gclpL6YL0EbSNhGxU1OVFxJkgoY9nzaNpGU29pKI+FH6+somC86xMZqZWWtTzOdkIqJO0oPAlZJOJXmg/khgSHbedJWWe4GvRMQrWadHkcwKrncBSdA5M1f9hcwuex64WlJN2pgakttez+e45oqM15sDfZo4zMwsVaP8jzydBXQEZpKs3HJmRIyT1FdSbbqiC8BIkmcfH0vTayU9DhARCyJiRv1BchtuUUR8kqviQnoy5wKPAh9JmgT0BT4CDs9xzYKM14dn7ZppZmaNKPYT/xHxKTC8kfTJZOzvFRH7F1Dm5fnkK2QK81RJO5E8oLMpyTTmV9IVAJoyUdJPSe79tU1nKzT46kXEnfm2w8ys2lXRqjIFT2FeAbxUwCXHABeSTBpoT/JsTYNiAQcZM7OUcu6a0rLkDDKSxkfE1unrKTSxvXJE9G0ifQLpHGtJz0TE0DVrrplZ9auiRZib7cmclvH6G2tSkQOMmVl+Ws1S/5n7w0TEc6VvjpmZ1VTRoEwhT/y3k3SFpPfTJZ7fT9+3L2UDzcxamxKswrzWFDLwfz3JnjJnAJOAfiRzqrsB3yl+08zMWqdq2rSskCDzdWCHdHkYgHckvU6ymqeDjJlZkVRRjCkoyDT1savoy2Fmtva1yjEZ4E/AI5K+LGnrdBG1h4H7S9IyM7NWSgUcla6QnsyFwKUk2272AqaTrIFzdQnaZWbWahV5Fea1qpBlZZYAl6WHmZmVSGsd+EfSliT7EHTJTPfaY2ZmxVNFMSb/ICPpYpJezFhWXV3Za4+ZmRVRa+3JnAfsGhH/LVFbzMyM1rV2WaaFwNulaoiZmSWqqSeTcwqzpJr6g+Tp/l9I6pmZXr9TppmZFUdrmsK8jM+X96//PKdmnFd6vk2R22Vm1mpV0xTm5nohA4DN0mNA1vvNMt6bmVmRSMr7yLO87pIeklQnaZKk45rId5Kk1yTNkzRV0vWS2qbn1pF0R3r9fEn/kXRwc3U3t9T/JCWfomtEzGukQd2A+Xl9SjMzy0sJhmRuAZYAPYDBwBhJYyNiXFa+TiSTvF4GNgJGAxcA15LEiynAvsBk4BDgfknbR8SHTVWcz3jKecCtORp+Th5lmJlZnmqkvI/mSOoMjABGRkRtuk/YaOCE7LwRcVtEPB8RSyJiGnAvsGd6ri4iLo+IDyNiRUQ8CnwAfDHnZ8nj854EXNHEuSuAU/Iow8zM8lTIfjKSTpf0asZxelZxWwDLI2JCRtpYYNs8mrIPkN3bSduoHmnZjZ6vl88U5n4R8W5jJyJioqT+eZRRNOvv8u1yVmfVrGO3td0CqypHFa2kQqYwR8QoYFSOLF2AuVlpc4GuzbThFGBnVp3sVX+uHUkv57cRkfPRlnx6MsvSiNVYI3oAy/Mow8zM8lRTwJGHWpLNJTPlHE+XNJxkHObgiJiVda4GuIdkjKfZv/rzaePfSAZ+GvNd4Nk8yjAzszwVeXbZBKCtpEEZaTvQ9G2wYcDtwOER8WbWOQF3kEwgGBERS5urPJ/bZZcC/5K0FfBn4COgJ8lA0hBgjzzKMDOzPLUt4iPuEVEn6UHgSkmnkswuO5Lk9/cqJB1AchvsKxHxSiPF3QZsDRwYEQvzqb/Zj5IOFu0CzCHpPj2a/juXZC2zRsdrzMxs9RT7ORngLKAjMJNkH7AzI2KcpL6SaiX1TfONBNYFHkvTayU9nrapH/AtkiA1I+P88bkqzmvtsoh4j0amu2WTdFFEXJtPmWZm1rhiP/AfEZ8CwxtJn0zG1i0RsX+OMiaxGivZFHvdsYuLXJ6ZWatTyBTmSlfQpmV5aAEf2cyssuXzkGVLUewgE81nMTOzXKppaftiBxkzM1tDVdSR8e0yM7NKU01L/Rc7yDxf5PLMzFqdKooxuYNM+mBOsyLi2fTfQ4rRKDOz1qw1DfzfkfW+N8ng/mxgA5LbY1PxxmVmZkVTRTGm2U3LBtS/lnQxSWAZGRELJHUCriQJOGZmViSt5nZZlu8AveoXREsDzQ+A6cD/lqJxZmatkapoDlUh07HrgF2z0nYBFhSvOWZmVqP8j0pXSE9mJPCEpEdI9nnuAxwGnF2KhpmZtVbVNIU5755MRNwD7AaMJ9nw5m1g9zTdzMyKpLX2ZIiItyS9DfSIiI9K1CYzs1atmmaX5d2TkbSepN8Di4CJadoRkq4uVePMzFqjGinvo9IVMvD/S5KNyvqR7O0M8BJwdLEbZWbWmrXW22VDSacwSwqAiPhE0salaZqZWevUAjooeSskyMwFNgRWjsWkW3Z6bMbMrIhqWulzMr8GHpC0P1AjaQ/gtyS30czMrEja1OR/5ENSd0kPSaqTNEnScU3kO0nSa5LmSZoq6XpJbQstJ1MhPZnrSAb9bwHaAXcCvwJuLKAMMzNrRgkG9G8hGUvvAQwGxkgaGxHjsvJ1As4DXgY2AkYDFwDXFljOSnkHmYgI4Ib0MDOzEilmjJHUGRgBbBcRtcALkkYDJwAXZeaNiNsy3k6TdC+wf6HlZMo7yORY9n8xMDUiJuVblpmZNa3IPZktgOURMSEjbSywbx7X7gPU91JWq5xCbpfdAfRKX9cv9Q8wE9hE0n+BYyLi3QLKNDOzLIXEGEmnA6dnJI2KiFEZ77uQTNzKNBfo2ky5pwA7A6euSTmFBpl1gcsiYqGkjsAVaSU3AD8FbgUOKqBMMzPLUsiMrDSgjMqRpZZkKbBM3YD5TV0gaTjJOMyBETFrdcuBwoLMuUDPiFgGkAaaS4DpEfEjSeeTbGBmZmZrQMW9XTYBaCtpUMadph34/DZYdt3DgNuBQyPizdUtp16hS/3vkpX2RT5f6n9FAWWZmVkTVMDRnIioAx4ErpTUWdKewJFAg8WN07H3e4EREfHK6paTqZCezGXAU+lsginApsDhwDnp+aHAnwsoz8zMGtGm+FOYzyJ57GQmyZj6mRExLn2g/i1gm4iYTLKly7rAYxm9qecj4uBc5eSquJApzHdLepVkClsvkq7THhHxVnr+UeDRfMszM7PGFTvGRMSnwPBG0ieTDOjXv99/dcrJpeCl/kminpmZlUiRx2TWqkKek7kHiEZOLSYZ8H84IsYWq2FmZq1VIYPlla6QzzKXZJBHJEFFwBHAcmBr4CVJJxa9hWZmrYykvI9KV8jtsi2AQyLin/UJ6SKZV0bEQem0txuAu4vbRDOz1qXyQ0f+Cgkyu5EsmpbpVWDX9PWTJDPOzMxsDbSEHkq+Crld9gbwI0kdANJ/ryJZuwZgAPBpUVtnZtYKtZHyPipdIUHmJGBvYJ6kGcA8ksXTTkrPdyeZQ21mZmugmA9jrm2FPCfzITAkfXinJ/BROse6/vyrxW+emVnr0wI6KHkr6DkZSB7ekTQFkKSaNM1LypiZFUmr3H5ZUq90283ZwDJgacZhZmZFIuV/VLpCxmR+RbLt5lCSJZ93Itma84wStMvMrNVSAf9VukJulw0B+kZEnaSIiLGSvgm8SLIstJmZFUFL6KHkq5Ags5zkNhnAHEkbkcww6130VpmZtWItYWpyvgoJMi8DhwAPkTx4eR+wEPh3CdplZtZqVVGMKSjInMDnYzjnAecD65PsPWBmZkXSEsZa8lXIczJzMl4vBK5On/qvA04tftPMzFqnmuqJMYU/J5MlaBkPnZqZtRitsieTQ2N7zJiZ2WpqVWMykg7Icbp9EdtiBVi/Wyd++cPjGbrHVsyeU8dlN43mvicaruzz9S9/kUvPOIQeG3Rj8dJlPPXPt/judX9ift2itdBqqxTrd+3ILy/+KkN3HcTsuXVcdttT3Pd0wz0Hjz94R876+hAG9tmA+XWLue+psVz2q6dYvjxZ5GPLfhtxwwVHsOOWvZk1p46Lb36c0f/w5rlrqrX1ZO5o5vzkZs5bCdzwg6NYsnQZ/Yb+gB223JQHbzqT/06Yyvj3Z6yS76U33uOAU37G7Dl1dO7YnpsvPZbLzz6M86//81pquVWCGy44giVLl9PvsGvYYVBPHvzJSfx34keM/2DmKvk6dWjPhTeO4ZVxU9hovc786foT+M78vfjJPf+gTZsa/nTdCfz64Zc59Nw72XvHATxw/YnsfvIvmDhl9lr6ZNWh2GMykrqT/C7/EjAL+EFE/L6RfNsBPwW+CGwQEco63x+4FdiDZFfkPwPnRcQymtDsE/8RMaC5o4AP2mhQk7RevmVY8j/+8KGDueLWMdQtXMKLb7zPmOfe5LjDdm2Qd+rHc5g9p27l++UrVrBZn43K2VyrMJ06tGP4fttyxe1PJz8//53EmBfGc9ywHRvkvf2hl/nn2A9Zumw502fN476n3mD37fsBSS+m54ZduemP/2TFiuC5197npTcnNVqOFaZGyvvI0y0kK7b0AI4HbpO0bSP5lgL3A99sopxbgZkkiyQPBvalmdX3y72V9B+UtRuPpA2AZ8vcjhZtUL+NWb58BRMnf/5X55sTprH1Zj0bzT9k8GbM+MePmfXizxg+dDA33/u3cjXVKtCgvhuyfEWs0tt4892P2HrAxs1eu9fgASt7O439ehOw7WY9itTS1quYS/1L6gyMAEZGRG1EvECyJNgJ2Xkj4p2IuAMY10RxA4D7I2JRRMwAngAaC1YrlTvILAV+Xf9G0sbA34ExZW5Hi9al0zrMrV11TGVu7UK6du7QaP4X33ifTfb5Hpt/6RJ+/ttnmDTde8u1Zl06NvLzU7eIrp3WyXndCYfuxE5b9eaG3z8PwDuTPuGTz+r47vF707ZNDUN3HcjeOw6gY4d2JWt7a1FIT0bS6ZJezThOzypuC2B5REzISBtLM8GhCTcCx0jqJKk3cDBJoGn6s6xGJWviRKCHpJ9L2oQkwNwXESNzXZT5RSxHIytd7YLFdMsKKN26dGh2MH/6J3N5+sW3uPvaU0rZPKtwtQsX063zqgGlW+cOzF+wuMlrDt9na646cxhHfvc3zJ67AIBly1dw1EW/Y9iQrfjw0Ys599i9eODZN5k2c15J298aFNKTiYhREbFzxjEqq7guwNystLlA19Vo2nMkwWkeMBV4FXg41wVlDTLp4NDXgB1JumN3RsTVeVy38otY6ja2BO9OmknbtjVs3vfzsZXtt+jN+Pc/avbatm1q2GzTDUvZPKtw706eRds2NWy+6QYr07YfuEmDQf96B+02iFu+/1W+9r27Gff+x6uc+7/3ZvCls29n04Ov5ojv/IYBvbrz6ltTS9r+VqG4W2PWAt2y0roB8wtqUrJ/2JPAg0BnYEOSVV+uy3VdyYOMpHsk3V1/AKNIougy4AsZ6ZanBYuW8Jdnx3LZmYfSqUN79thhMw7b9wv8/tFXGuQ95uCd6bPJ+gD07bk+l3/7cP7+yjvlbrJVkAWLlvKX597istMOpFOHduyxfV8O23sbfv/Efxrk3feLm3HX5Udz7CX38ur4hsFju803YZ32bem4TjvOO3YvNtmgK/c89lo5PkZVK/JS/xOAtpIGZaTtQNPjLk3pDvQBbo6IxRExG7iLZE3LJhXjYczmTGwi/fUy1F21zr3mPn51+fFMfvZ/+XROHedecx/j359Bn03W5/UHLmWnEVczZcZnbLVZT64+90jW69aJOfMW8MQLb3HZL0av7ebbWnbuj//Cry4ZweQxl/Dp3AWc++O/MP6DmfTpsS6v33seOx1/A1M+nssPTj6AdTuvw8M/OWnltf8c+yHDz/8tAMcNG8zJh+9Cu7Y1/HPshxx67p0sWbp8bX2sqlHMhzHT7VkeBK6UdCrJrLAjSbZvyapXAtYhfQYyXTos0qAyS9IHwJmSfkJyG+4kkvGdpj9LRMt6YL/jjt9uWQ22ytUx+w6C2epb+OI1RQsN//5gbt6/53YZsG6z9abPydwJHATMBi6KiN9L6gu8BWwTEZPT52A+yLp8UkT0T8sZDNxA0hNaDvwNODsiGr/XSnl6MitJugh4JiL+nZG2K7BfRFxfzraYmVWqYj/xHxGfAsMbSZ9M0iOpf/8hOUZ6IuINYL9C6i737LJzSaJmprdItg4wMzOS22X5HpWurD0Zkvt8S7PSlgCNP+BhZtYKtYDYkbdy92Reo+ESBGfgSQBmZp8r7hTmtarcPZnvAE9LOgF4DxhIspbOQWVuh5lZxWptqzAXTUSMk7QFcBjJfOsHgUcjorac7TAzq2QtYawlX+XuyRARtZJeBHoD0xxgzMxWVU1BpqxjMpJ6SnoOeJekFzNR0j8k9SpnO8zMKlmRn/hfq8o98H8bydOh3SOiJ8m6N/8BflnmdpiZVSxPYV59ewE9I2IprFzu4EJgWpnbYWZWsVpA7MhbuXsynwHbZKVtCcwpczvMzCqXpzCvtuuBv0q6A5gE9ANOAXLuJ2Nm1pq0hLGWfJV7CvPtkt4DjgO+AEwHjo0Ib79sZpZqCWMt+VobU5ifBRxUzMya4CCzmiS1BY4l2RmzS+a5iMjel9rMrFXy7bLV9ztge+Bx4ONm8pqZtUruyay+YUCfiChob2kzs9akimJM2acwv0WyT7SZmTXFU5hX2zeAX0t6iqzbZRFxd5nbYmZWkappTKbcPZmTgb2Bo4HTMo5Ty9wOM7OKVexlZSR1l/SQpDpJkyQd10S+7SQ9KWmWpGgizzGSxqdlvSdp71x1l7sncy6wY0SML3O9ZmYtRgn6MbeQ7ELcAxgMjJE0NiLGZeVbCtwP3Ao83KBd0kHAdSQdhVeAns1VXO4g8zEwucx1mpm1KCri9DJJnYERwHbp1iovSBoNnABclJk3It4B3pE0sInirgCujIh/pe+bXXey3LfLfg7cK2l3SZtlHmVuh5lZxSrkdpmk0yW9mnFkP3O4BbA8IiZkpI0Fti2sTWoD7AxsJGmipKmSbpbUMdd15e7J3JL+e0RWegBtytwWM7OKVEg/JiJGAaNyZOkCzM1Kmwt0LbBZPYB2wNdIxtaXAn8BLgUuaeqisvZkIqKmicMBxsysXnGnMNcC3bLSugGFPq+4MP33FxHxUUTMAn4GHJLronLfLjMzs2YUeWfMCUBbSYMy0nYAsgf9c4qIz4CpJHee8lbutcv6Aj+k8bXLtihnW8zMKlUxl5VJN4d8ELhS0qkks8uOBIY0rFcC1gHap+87JEXE4jTLXcA5kp4guV12HvBorvrLPSbzJ+Bt4DI+73qZmVmGEkxhPgu4E5gJzAbOjIhx6R/+bwHbRMRkkj2+Psi4biHJ3l/90/dXARuS9I4WkUx3/lGuissdZLYC9oiIFWWu18ysxSjmFGaAiPgUGN5I+mQy7ipFxIfkiHERsZQkYJ2Vb93lHpN5BNi3zHWambUoxX7if20qd0/m/wEvprtjZq9d9j9lbouZWUVqAbEjb+UOMncBy4HxeEzGzKxRLaGHkq9yB5kDgF7eT8bMrGlehXn1/RfYoMx1mpm1LN5PZrU9Czwl6S4ajsncWea2mJlVpBYQO/JW7iCzF8mqnV/KSg+SOdxmZq1eTRUNypQ1yETE/uWsz8ysRaqeGFP2nsxK6fIFK7+UfkDTzCxRRTGmvAP/knqnW4DOBpaRrH1Tf5iZGdX1MGa5Z5f9kmQL0KEky0/vBIwGzihzO8zMKlaRV2Feq8p9u2wI0DddFTQiYqykbwIvAreXuS1mZhWpJfRQ8lXunsxykttkAHMkbQTUAb3L3A4zMyuDcvdkXibZRe0h4EngPpLlZV4tczvMzCpWNfVkyh1kTuDz3tN5wPkky0zfUOZ2mJlVLD8ns/oWACdLGsznexgIuB44scxtMTOrSNUTYsofZH5Lsrf0I2QtK2NmZqkqijLlDjLDgAERMafM9ZqZtRgtYWpyvso9u2wysE6Z6zQza1GK/TCmpO7pg/B1kiZJOq6JfNtJelLSLEmRo7xBkhZJ+l1zdZe8JyPpgIy3dwN/kXQjDVdhfrbUbTEzawlK0I+5heRB+B7AYGCMpLERMS4r31LgfuBW4OFmyvt3PhWX43bZHY2kXZP1PoDNytAWM7PKV8QoI6kzMALYLiJqgRckjSaZ7XtRZt6IeAd4R9LAHOUdA8wheYi+yXz1Sh5kImJAqeswM6smRZ7CvAWwPCImZKSNBfYttCBJ3YArSZYG+2Y+16y1VZhX18L/3Fw9I2IlJOn0iBi1ttth1cE/T+XVoW3+fRlJpwOnZySNyvpedQHmZl02F+i6Gk27CrgjIqYoz0DY4oKM5e10wL8UrFj881Sh0oCS63tTC3TLSusGzC+knvT5xgOBHQu5zkHGzKy6TQDaShoUEe+maTsA2YP+zdkP6A9MTnsxXYA2kraJiJ2auqjcU5jNzKyMIqIOeBC4UlJnSXsCRwL3ZOdVogPQPn3fQVL9YyejgM1JZqcNJtm6ZQzw5Vz1O8hUL9/asGLyz1PLdhbQEZgJ/AE4MyLGSeorqVZS3zRfP5JFi+t7OQuBdwAiYkFEzKg/SG7DLYqIT3JVrIgmn7cxMzNbI+7JmJlZyTjIVBBJW0r6j6T5klZIGpnHNR9KOrCJc3tLeqf4LbW1Kdf3vKWRdHk+S5NYy+XZZZXlQuDvEVHQFMGmRMTzwJbFKMvMbHW4J1NZ+lH4tEIzs4rlIFMhJD0L7A/cnM72+L2kq9NzG0p6VNIcSZ9Kel5S5vdusKT/Spor6b50CiKS9pM0NaOODyVd0Fje9PyFkj6SNF3SqZIi1xpGtlY1+J5LWj/9OflE0mfp603rL5D0d0lXS3ox/Rl7RNIGku6VNE/SvyX1z8gfks6S9G56C/cqSZtLeinNf7+k9hn5T5M0Mf0ZHS2pV8a5bSU9nZ77WNLF2R9IUjtJf5D0QGa51rI5yFSIiDgAeB74dkR0IVkxtd75wFRgI5JVVC8mWVS03lGke/UAXwBOzlFVo3klDQO+S/JE70BWY10jK6vGvo81wF0kPeK+JNNPb8667hiShRF7kzzz8FJ6TXdgPPDDrPzDgC8Cu5Pczh0FHA/0AbYDjoWVq63/b9qunsAk4I/pua7AX4EngF4kP1/PZFYiqSPJqr+LgaMiIvPn31owB5mWYSnJ/7j9ImJpRDwfq849vykipkfEpyS7jg7OUVZTeY8C7oqIcRGxALii6J/CiqnB9zEiZkfEA+nzDPOBH9Hwj4W7IuK9iJgLPA68FxF/jYhlwJ9ouGTIdRExL10S/v+ApyLi/Yzr6/MfD9wZEa9HxGLgB8Aeac/oMGBGRPw0IhZFxPyIeDmjjm4kAeg94JSIWF6Ur5BVBAeZluHHwETgKUnvS7oo6/yMjNcLSJZ7aEpTeXsBUzLOZb62ytPg+yipk6RfKdmUah7wD2A9SW0y8mbu47SwkffZPzv55u9F0nsBIF1SfjZJj6kPSQBpyu4kvbFrs/54sirgINMCpH/5nR8RmwGHA9+VNLTI1XwEbJrxvk+Ry7fSO59kNuFuEdEN2CdNL8fK5dNJbtMlFSZ7mGwATCP5g2XzHNc+RXKr7RlJPUrZSCs/B5kWQNJhkgYqWZVuHrA8PYrpfuAUSVtL6gRcVuTyrfS6kvQu5kjqTsPxlVL6PcnPz+B0ratrgJcj4kPgUWATSedJWkdSV0m7ZV4cEdenZTwjacMytttKzEGmZRhEMnBaSzJQe2tE/L2YFUTE48BNwN9Ibs29lJ5aXMx6rKRuIFmfahbwL5JxjrKIiGeAkcADJL3izUkmGZCODx1E0gufAbxLMpMyu4yrSAb//5oGSasCXrvMGiVpa5KB3nXSQWEzs4K5J2MrSfqKpPaS1geuAx5xgDGzNeEgY5m+BXxCMhNoOXDm2m2OmbV0vl1mZmYl456MmZmVjIOMmZmVjIOMmZmVjIOMVZU12QRL0i/z2SjOzPLnTcusRZJ0HMmq0VsB84E3SBaEXG0Rccaat8zMMjnIWIsj6bvARcAZwJMk2yIMA44E6tZi08wsi2+XWYsiaV3gSuDsiHgwIurS7Q8eiYjvpdnaS7o73WhrnKSdM67fOt28a0567oiMc7+p3ygufX+kpDfSDbreS/fcQdK6ku5IN3iblm4ElrnSsZmlHGSspdkD6AA8lCPPESQbZq0HjCbduEtSO5K9V54CNgbOAe6VtGV2AZJ2Be4GvpeWsw/wYXr6t8Ayks23dgS+BJy6Jh/KrFo5yFhLswEwq5nlbl6IiMfSza/uAXZI03cn2f/k2ohYEhHPkqwQfGwjZXyTZBOupyNiRURMi4i306XoDwbOS3tRM4Gfky4GaWar8piMtTSzgQ0ltc0RaLI39OogqS3pxmwRsSLj/CSSjbWy9QEeayS9H9AO+CjZeQFI/ljzJm9mjXBPxlqal4BFwPDVuHY60EdS5s99X5KNtbI1tdHWFJLtDzaMiPXSo1tEbLsa7TGreg4y1qKke8tfBtwiaXi65XA7SQdLur6Zy18mmX12YXrNfiR7nPyxkbx3kGzCNVRSjaTekraKiI9IxnR+Kqlbem5zSfsW7UOaVREHGWtxIuJnJM/IXEqyavQU4NskG17lum4JyaSAg0k29roVODEi3m4k7yvAKSTjLXOB5/h8e+ETgfbAW8BnwJ+Bnmv4scyqkldhNjOzknFPxszMSsZBxszMSsZBxszMSsZBxszMSsZBxszMSsZBxszMSsZBxszMSsZBxszMSsZBxszMSub/A/y9kHWKtmivAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "data = pd.crosstab(df.Lagged_Choice_1, df.Choice, normalize=True)\n", "sns.heatmap(data, cmap=\"Blues\", annot=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "We now study how Robinson's behavior changes as we increase the returns to experience. We do so by plotting the average level of final experience in the sample under the different parameterizations.\n", "\n", "This analysis of the comparative statics of the model is straightforward to implement. In models of educational choice, this type of analysis is often applied to evaluate the effect of alternative tuition policies on average educational attainment. See Keane & Wolpin (1997, 2001) for example. The basic structure of the analysis remains the same." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Specification of grid for evaluation\n", "num_points = 15 \n", "grid_start = 0.0\n", "grid_stop = 0.3\n", "\n", "grid_points = np.linspace(grid_start, grid_stop, num_points)\n", "\n", "rslts = list()\n", "for value in grid_points:\n", " \n", " params.loc[\"wage_fishing\", \"exp_fishing\"] = value\n", "\n", " df = simulate(params)\n", "\n", " stat = df.groupby(\"Identifier\")[\"Experience_Fishing\"].max().mean()\n", " rslts.append(stat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We collected all results in `rslts` and are ready to create a basic visualization." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2021-06-01T13:22:54.639519\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.2, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEUCAYAAABkhkJAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqLElEQVR4nO3deXxU9b3/8dcnAcKSQIAsIBBQFkGoQIta61J31NbW2+Wqta1bS6tVa+1ib7drXdrb5SpSbS22VrQubdXaulTbgstPvS6ggKIIiOxLwhoSCAnJ5/fHOROGmIQTMpPZ3s/HYx6ZOefMmc/Jgfnke76f8/2auyMiIpJqeakOQEREBJSQREQkTSghiYhIWlBCEhGRtKCEJCIiaUEJSURE0oISkoiIpIUuS0hmdrmZzTWz3WZ2V4t1J5vZYjPbaWZPm9nwropLRETSQ1e2kNYBNwB3xi80sxLgYeCHwABgLvCnLoxLRETSQLeu+iB3fxjAzKYAQ+NWfQpY5O5/CddfC2wys7Huvrir4hMRkdRKhz6k8cCC2At3rwXeDZeLiEiO6LIWUjsKgaoWy7YDRS03NLNpwDSAPn36fGjs2LHJj05ERBJm3rx5m9y9tLV16ZCQaoC+LZb1BXa03NDdZwIzAaZMmeJz585NfnQiIpIwZrayrXXpcMluETAx9sLM+gAjw+UiIpIjurLsu5uZ9QTygXwz62lm3YC/AhPM7NPh+h8BC1XQICKSW7qyhfQDYBfwXeDz4fMfuHsV8GngRmArcBRwbhfGJSIiaaAry76vBa5tY92/AVUoiIjksHToQxIREVFCEhGR9KCEJCIiaUEJSURE0oISkoiIpAUlJBERSQtKSCIikhaUkEREJC0oIYmISFpQQhIRkbTQoYRkZnlmNjhZwYiISO6KlJDMrNjM7gPqgGXhsk+Y2Q3JDE5ERHJH1BbS7QSzuA4H6sNl/weck4ygREQk90Qd7ftk4CB3bzAzB3D3KjMrS15oIiKSS6ImpO1ACbA+tsDMKuJfi4jkmqYmp25PI7vqG9nV0EhdQyO76pvY1RC83lUfLKuLvW5opGGP0+iOu9PYFHsOjU1OkztN4bImD/bf5E5jE8G68D2x7eP30+TBNu7gOADue2ONPW+5rnmT/az3uJ3d8cUpDCwsSPBvM3pC+h3wkJl9H8gzs6OBnxBcyhMRyTh7GpvYVFNP5Y46NlbvZmN1HZXVdWzb1bBvgmmIJZymfZNLfSO79zSl+jBSoqHR97/RAYiakH5GUNBwG9AduBP4LXBLUqISETlAjU3O5prdbKzevW+yCZ/Hfm6q2b1PC0JSL1JC8qCtNj18iIikxO49jSzdWLM3uVTvZuOOoGUTSzZVO3bT1IWJplf3fHr1yKdntzx69sgPXseWdd/3dUH3PHrk55FnRn6ekWeQl2fkm5FnRl64LFhn4XZgFmyTn2fY+9bv3U+wLIjLsOYYrXkZ+zyJbdNyvVlby4Of/ft0T+jvMCZSQjKz7wKz3f3VuGVHAie4+8+TEpmICFC7ew/PvFPFk4s28PTiSmp270nYvksKe1BW1JOyvgWUF/WkvG8BA/r0oHePbs3JpWf3vPBnkFT2STDd8pq/vKXzol6y+zrwqxbL3gIeAZSQRCShtu2s599vV/Lkmxt4bmkV9R3sqxnQpwdlRQWU9+3Z/LO8bwFlfXs2LystKqB7vgarSSdRE1IPoKHFsnqgZ2LDEZFcVVldx1NvbeSpNzfwf8s309jGdbeh/XtxcEmfNpNNaWEBPbop0WSiqAlpHnAZ+/YhfRV4LdEBiUjuWL1lJ0++uYEnF23gtVVb2ywyGDuoiNMnDOKMCYMZU16oy2RZKmpC+gbwLzP7AvAuMAooB05NVmAikn3cnWWVNfzjzQ08+eYG3lpf3ea2kyuKOX38IKaOH8SIkj5dGKWkStQqu0VmNgY4CxgKPAw85u41yQxORDKfu/PG2u3NLaHlVbWtbpefZxx18ABOnzCI0w4bxKB+6hHINVFbSITJ5/4kxiIiWaKxyZm7YgtPLtrAU29uYN32ula365Gfx7GjSzh9wiBOGVfOgD49ujhSSSdRy74PBm4EJgGF8evcvSLxYYlIJlq5uZbbn32Xfy7ayOba+la36d0jnxMPLWPqhEGceGgpRT2Tc0+LZJ6oLaT7CPqOvgnsTF44IpKp3ttUy3/8+gW27WxZkAv9enXnlHHlnD5hEMeNLqFn9/wURCjpLmpCGg8c4+65OXCTiLRra209F/3hlX2SUWlRAVPHl3P6+MEcdcgA3fMj+xU1IT0HTCYo/xYRabZ7TyNfuWceKzYHF08KuuVx++c/xEfHlJKXp/JsiS5qQloBPGVmDwMb4le4+48SHZSIZAZ355oHF/LKii1AMNbZ9HMmceJYTZUmHRc1IfUBHiUY6XtY8sIRkUwy/d9LeWT+uubX3z19LGd8YHAKI5JMFvU+pIuSHYiIZJa/vr6GW2YvbX593pEVTDv+kBRGJJku8n1IZjYO+AxQ7u6Xm9mhQIG7L0xadCKSll5evpnvPLj3v/5xo0u47pPjNaSPdEqkshcz+yxBYcMQ4Ivh4iLgpiTFJSJpanlVDV/547zmWUPHlBdy2/kfVBWddFrUf0HXAae6+1eBxnDZAmBiogIxsxFm9oSZbTWzDWZ2q5lFbsGJSPJtqa3n4rtebS7vLiks4M4Lj6Cvbm6VBIiakMoIEhCAx/1M5LyMvwYqgcEEI0J8lGCEcRFJA0F599zm8u6e3fP43QVTGNq/d4ojk2wRNSHNA77QYtm5wCsJjOVg4M/uXufuG4AnCW7IFZEUc3e+8+BCXl2xFdhb3j1pWHFqA5OsEvWS2JXAP83sEqCPmT0FjAFOS2AstwDnmtkzQH/gDOCHCdy/iByg6f9eyt/iyrv/64yxnD5B5d2SWFHLvheb2Vjg48BjwGoSP/3Es8CXgWogH5hFMEV6MzObBkwDqKjQmK4iXeHh1/Yt7/7cURV8+TiVd0viRS6Lcfed7v5nd/+Fuz+QyGRkZnnAUwTzLPUBSghaST9rEcNMd5/i7lNKS0sT9fEi0oaXl2/mmof2Le/+8SdU3i3J0WYLycyedPfTw+f/jzYKGNz9+ATEMYBgBIhb3X03sNvM/gDcAHwnAfsXkQ5aXlXDtHv2lncfWl6k8m5JqvYu2d0d9/x3yQzC3TeZ2XvApWb2S4I5ly5gb2WfiHShWHn39l17y7t/f+EUlXdLUrWZkNz9PgAzywdGAjeGrZdk+RQwHbiG4F6np4FvJPHzRKQVdQ2NTLt73/Lu36u8W7rAfosa3L3RzL4GXJvMQNx9PnBCMj9DRNrn7lzz0ELmrowv757MRJV3SxeIejF4FvDVZAYiIql3c4vy7u+dMY7TJwxKYUSSS6Leh3QkcIWZfYeg5Lu5wCFBRQ0ikmIPzVvDjLjy7vOPquBLxx2cwogk10RNSHeEDxHJQi8t38x3H95b3n38mFKVd0uXi3pj7KxkByIiqfFuVQ1fiSvvHjuoiNs+N5luKu+WLhZ1+gkzsy+b2RwzWxguO97M/jO54YlIMrUs7y4tKuD3Fx5Bkcq7JQU6Mv3EJcBMIDZmzxqCEm0RyUCx8u6VLcq7hxT3SnFkkquiJqQLgY+7+wPsLWh4D9CAViIZKDZ6d3x59y3nTubwocWpDUxyWtSElA/Exq6LJaTCuGUikkFu/tcS/r5gb3n3988cx9TxKu+W1IqakJ4AbjKzAgj6lIDrgUeTFZiIJMdD89YwY86y5tfnH1XBJceqvFtSL2pCuho4CNgO9CNoGQ0HvpukuEQkCVqWd39U5d2SRqKWfVcDZ5tZGUEiWh3O6ioiGaKyuo7L73ttn/LuW1XeLWkkatn378yst7tXuvur7r7BzAab2ZPJDlBEOm9PYxNXPvA6m2rqgdjo3SrvlvQS9U+jImChmR0NYGbnAguB15MVmIgkzozZS3lp+RYA8gxmnDdJ5d2SdqJesjvHzM4H/mZm7wCDgbPd/YWkRicinfb/llbxq6f3FjF8/eQxfGRkSQojEmldRy4erwXqCO49eg94NykRiUjCbKyu46oH5uPhzRrHjBrI5SeNSm1QIm2I2of0S+AB4EpgBDCf4BLeZ5MWmYh0yp7GJq68/3U21wb9RqVFBUw/ZzL5eaqok/QUdbTvccBEd98Yvv62mT1KME/SX5ISmYh0yozZS3n5vb39RrecO4nSooIURyXStqh9SB9rZdlzZnZ44kMSkc56bon6jSTzRL1kV2BmN5rZcjPbHi47DbggqdGJSIdtrK7jG3/a22907KgS9RtJRoha1DAdmACcz96x7BYBlyYhJhE5QK31G918ziT1G0lGiNqHdDYwyt1rzawJwN3XmtmQpEUmIh12i/qNJINFbSHV0yJ5mVkpsDnhEYnIAXluSRW3xvUbXXWK+o0ks0RNSH8BZpnZwQBmNhi4laAUXERSrLV+o6+dqH4jySxRE9L3gBXAG0AxsBRYB/w4KVGJSGR7Gpu4Qv1GkgWiln3XA1cBV4WX6ja5x/4WE5FUmv7vpbyifiPJAlGLGpq5e1UyAhGRjnt2SRW3PaN+I8kOmghFJENt2L5vv9Fxo9VvJJlNCUkkA8XmN9oS9huVqd9IskCbCcnMfhH3/KSuCUdEonh/v9FkSgrVbySZrb0W0rS4548kOQ4Riahlv9E3ThnD0SMHpjAikcRor6hhgZk9CLwFFJjZda1t5O4/SkpkIvI+rfUbXaZ+I8kS7SWkzxC0koYDBgxrZRuVfot0kdg4deo3kmzVZkJy90rgBgAz6+buF3VZVCLyPjf/ewmvrNjbbzTjPPUbSXaJemPsRWbWHzgLGEIwnflj7r4lmcGJSODZJVXc9vS7za+vPnUMHz5E/UaSXaLOh3Q08C7wVeBw4CvAsnB5wpjZuWb2tpnVmtm7ZnZcIvcvkoli/UYxx40u4bIT1G8k2SfqSA3TgcvcvXkwVTM7B5gBHJGIQMzsVOBnwDnAK8DgROxXJJO11W+Up34jyUJRb4wdA/y5xbIHgUT+mfZj4Dp3f8ndm9x9rbuvTeD+RTLOTf9Sv5HkjqgJaSlwbotlnyW4jNdpZpYPTAFKzWyZma0xs1vNrFci9i+SiZ55p5JfP6N+I8kdUS/ZXQU8ZmZXAiuBEcBo4OMJiqMc6E5Qan4c0AD8DfgB8P3YRmY2jfCG3YqKigR9tEj6Wb99F1f/eUHza/UbSS6I1EJy9xeBkQST8s0DfkUwpfmLCYpjV/jzV+6+3t03ATcBZ7aIY6a7T3H3KaWlpQn6aJH00rLfqLyv+o0kN0SefsLdtwJ/TEYQ7r7VzNagG21FuOlfS3h1xVYg7DfSOHWSI9JptO8/AFeYWVl4z9NVwGOpDUmka/19wbp9+o2+edqhHKV+I8kRHZ6gL4muB0qAJUAdQVXfjSmNSKSLuDu/fW45//OPxc3LjhtdwqUfHZnCqES6VtokJHdvAC4LHyI5Y09jEz/6+yLue3lV87JDSvqo30hyTtokJJFcVLN7D1+79zWeXVLVvOzIEQOY+cUPUdy7RwojE+l6bSYkM1tNhCIDd1f9tcgBWL99FxffNZe311c3L/vkpIP4+WcOp6BbfgojE0mN9lpIn++yKERyzKJ127n4rlfZWL27edmVJ43iG6eOwUyX6SQ3tTf9xLNdGYhIrnj6nUouv/c1ausbAeiWZ/zkUx/gP6e0NuWYSO6IOtp3gZndaGbLzWx7uOw0M7s8ueGJZJc/vrSSL82a25yMigq6MeviI5WMRIh+H9LNwATgfPb2Ky0CLk1GUCLZpqnJ+ekTb/ODR96ksSn4LzSkuBcPXfYRjhlVkuLoRNJD1Cq7/yAYKqjWzJoA3H2tmQ1JXmgi2aGuoZGr/zyfJ97Y0Lzs8KH9+N0FUygr6pnCyETSS9SEVN9yWzMrBTYnPCKRLLK5Zjdfunsur6/a1rzslHHlzDhvEr176K4LkXhRL9n9BZhlZgcDmNlggoFWH2j3XSI57N2qGv7j1y/uk4wuOmYEv/3Ch5SMRFoRNSF9D1gBvAEUE8yPtI5gUj0RaeGl5Zv51K9fZNWWnQCYwX+fdRj/fdZ48jX6gkirIv2Z5u71BIOdXhVeqtvk7hqZW6QVf319Dd95cCENjcF/kV7d85lx3mROPaw8xZGJpLdICcnMHgHuBf7u7lX72VwkJ7k7v5qzjJv+taR5WUlhAXdeOIXDhxanLjCRDBH1kt2zwLeBSjObZWZTzSydpq4QSan6PU18+8GF+ySjMeWFPPK1jygZiUQUdcbYm939SGAKsByYDqwzsxlJjE0kI2zf1cCFf3iFB+etaV52zKiB/OWrH2Fo/94pjEwks3SolePuS939x8C5wELga0mJSiRDrN6yk8/85kVefHfvHRCf/dBQ/nDhkfTr1T2FkYlknsi1p2Y2EjgvfJQADwLXJSkukbS3YPU2Lpk1l001ewdI/dZpY/jaiaM0QKrIAYha1PAqMAb4O/At4J/u3pjMwETS2T8XbeDKB16nrqEJgB75efzis4fzyUkavETkQEVtIf2SoMJuVzKDEUl37s6dL6zghsffInbjQ79e3Zn5hQ9x1CEDUxucSIaLeh/Sn8xsoJl9Bhjs7j83s4OAPHdfs7/3i2SDeSu3cMPjb+8z8kLFgN784aIjGFlamLrARLJE1Et2HwUeAuYCxwA/B0YTXL47K2nRiaSBlZtr+dmTi/cZHBVgckUxv/viFAYWFqQoMpHsEvWS3XTgHHefbWZbw2UvA0cmJSqRNLBtZz0zZi/jnpdWNI+6AEF/0YXHjODqU8fQs7umGhdJlKgJaYS7zw6fx/5nvm8EcJFssHtPI/f830pmzF5Kdd2efdadNfEgvjP1UIYN0P1FIokWNaG8ZWZT3f2puGWnEAy2KpIV3J3H31jPz55czOot+9bvTBnen+9/bByTK/qnKDqR7Bc1IX0TeMzMHgd6mdlvCfqOPpm0yES60NwVW7jxiX0LFgBGDOzNd88Yy9Txg3RvkUiSRa2ye8nMJhJMYX4nsBo4UhV2kulWbAoKFv7x5r4FC8W9u/P1k0dz/lHD6dFNwzaKdIXIfUDuvpaguk4k422tredXc1ovWLjomBFcduIoDf0j0sXaTEhmdg97Cxja5O5fTGhEIkm0e08jd7+4kl/NUcGCSLppr4W0rMuiEEkyd+exhev5+VPvL1g4YkR/vnemChZEUq3NhBSO6i2S8eauCEZYmL962z7Lg4KFcUwdX66CBZE0oPuIJGupYEEksyghSdbZWlvPjDlL+eNLK1WwIJJBlJAkK7g7r6/exv0vr+LRheuap4WI+cTEg/i2ChZE0poSkmS07bsa+Nv8tdz38ioWb9jxvvVHjOjP9z92GJOGFXd9cCLSIe2VfV8cZQfufmfiwhHZv/21hgDGDiriqlPGqGBBJIO010L6QoT3O8HIDQlhZqMJxsd70N0/n6j9SnbYX2uoV/d8zpo4mPOOrGDSsGIlIpEM017Z94ldGUjoNuDVFHyupKmoraHzj6rgk5OH0LenihVEMlWH+5As+LOz+U9Pd3//N8QBMLNzgW3Ai8CoROxTMld1XQOPvK7WkEguiTpj7BDgVuB4oLjF6k7PUGZmfYHrgJOBSzq7P8lM7s781du4T60hkZwUtYV0O7CTIGE8S5CYrgWeSFAc1wO/d/fV7f2la2bTgGkAFRUVCfpoSTW1hkQEoiekjwAV7l5rZu7uC8zsEoLLa3d0JgAzm0Qw2d/k/W3r7jOBmQBTpkzZ78Cvkr7UGhKRlqImpEYgNjTyNjMrBaqBIQmI4QRgBLAq/Mu3EMg3s8Pc/YMJ2L+kCXdnWWUNcxZX8tfX16o1JCL7iJqQXgbOBP4KPAX8CdgFzE1ADDOBB+Jef4sgQV2agH1LitU1NPLye1uY8/ZG5rxT+b6RtmPUGhKRqAnpC0BsFMqrCJJGITC9swG4+06C/ikAzKwGqHP3qs7uW1JjY3UdTy+uZPbiSl5Ytomd9Y2tbqfWkIjEizqF+ba457sIihCSwt2vTda+JTmampw31m5n9uJK5izeyJtrq9vctrCgG8ePKeHEQ8uYOmGQWkMi0ixq2XcP4EJgEkHLqJlmjM1NNbv38PzSKma/XcnT71SxqWZ3m9seXNKHk8aWcfLYMqaMGKApH0SkVVEv2c0CJgKPAhuTF46ksxWbapm9uJKnF1fy8nub95naIV63POPIgwdw0tgyThpbxiGlha1uJyISL2pCOh04OP7SnWS/hsYmXl2xhTlvVzLnnUqWV9W2ue3APj044dAyTh5XxrGjS3QpTkQ6LGpCWgUUJDMQSb36PU28vb6a+au38cp7W3huSRU7du9pc/vDBvfl5HFlnDi2jIlDi8nPU1GCiBy4qAnpbuBvZnYLLS7ZufuchEclSefurNm6i9dXb2P+qm3MX72VN9dVU7+n7aEJe3bP49hRJZw0tpwTx5YyuF+vLoxYRLJd1IR0efjzJy2WO3BI4sKRZNlR18DCNduZv3obr4cJaFNN/X7fN6S4V9AXNK6Mow8ZSM/unR66UESkVVHLvg9OdiCSOI1NzpKNO8Lks5X5q7extLIGjzDY0rABvZg8rD+ThhVzzKgSxpQX6v4gEekSmsI8C1RW1/F6XMvnjTXbqW3jZtR4RQXdmDismEnDiplcUczEYcWUFKqrUERSo70pzN9293Hh89UEl+fex9017HYX2r6zgUXrt7NobXVzC2jd9rr9vi/P4NBBfZlcESagYcWMLC0kT4UIIpIm2mshfTnuuaYT72LuTuWO3by5djuL1lWzaF3wc83W1seCa2lQ355MGlbMpIog+XxgaD9691CDWETSV3vfUL8EPhw+P8Hdf9wF8eSkpiZn1ZadvLkulnyqeWvd9khFBxBUvx0+pLi59TOpolgVcCKScdpLSGPMrKe71wHfBJSQEqChsYlllTXNLZ+31lXz1vpqatq53ydej/w8xgwqZPzgfnxgaD8mVxRzaHkR3fI1HI+IZLb2EtLfgCVmtgLoZWbPtbaRux+fjMCywc76Pby9fgdvxbV83tm4o917feL16ZHPYQf1ZfxB/Rgf/hxVVqix4EQkK7WZkNz9IjM7lmBuoiOA33dVUJlm955G3ttUy5KNNSzbuIMlG2tYUrmDFZtqaYo4r+3APj2ak8+EIcHP4QN6q+hARHJGu73c7v488LyZ9XD3WV0UU9ravaeR5VW1LK2sYenGHSwNE8/KzTtpjJp5CG42Hd8i+ZT3LdD9PiKS06LeGHtnsgNJJ3UNsRbPDpZV1rAkTD4rt3Qs8eQZHFJaGCafvZfeinv3SGL0IiKZKafrgOsaYi2eIOHEfq7YHP1SW8ywAb0YXVbE6PJCxoQ/R5UVqtRaRCSinPy2vP+VVcx8bjkrDzDxjCkrYlSYeMaUFzGyrI8Sj4hIJ+Xkt+ieJue9TW3P7WMGw/r3ZnRZIaPLixhTXsjoMiUeEZFkijqFuQFfAs4DStz9cDM7Hhjk7n9OZoDJMLosmME0lnjGlAeJZ3RZYdDiKS2kVw+Nai0i0pWi/rl/HXAqMB24PVy2BrgZyLiENHFoMY9dcawSj4hIGomakC4EJrv7JjP7TbjsPTJ0LqRePfKZMKRfqsMQEZE4UW/5zwdqwuexMoDCuGUiIiKdEjUhPQHcZGYF0NyndD3waLICExGR3BI1IV0NHARsB/oRtIyGA9ckKS4REckxUUdqqAbONrNyoAJY7e4bkhqZiIjklKhl37GWVFX4wMzy3D3asNUiIiL7EfWS3R6goeXDzHab2Xtm9r9mVpisIEVEJPtFTUhXAHOA04BxwFRgNvAd4FLgIwT3KImIiByQqPchXQ180N23h6+XmNlcYJ67jzSzN4B5SYlQRERyQtQWUl+gd4tlvQkq7gA2AL0SFZSIiOSeqC2ku4F/mdktwGpgKPB1IDZp32nAO4kPT0REckXUhPRtYClwLsH9SOuB24A7wvVPA88kOjgREckdUe9DaiIYVPX2NtbXJTIoERHJPZEn9wlvij0SKAEstjwR05uHQxL9GjgFGAAsA77n7v/o7L5FRCQzRL0x9mzgjwSX7cYDi4AJwPNApxNSGMdq4KPAKuBM4M9m9gF3X5GA/YuISJqLWmV3A3CRu08GasOf00hQqbe717r7te6+wt2b3P0xguktPpSI/YuISPqLmpAq3P0vLZbNAr6Y4HiA5suDYwhaYiIikgOiJqTKMEkArDCzo4GRBPMkJZSZdQfuBWa5++IW66aZ2Vwzm1tVVZXojxYRkRSKmpDuAI4Nn99MUOa9gKAQIWHCQVzvAeqBy1uud/eZ7j7F3aeUlpYm8qNFRCTFolbZ/SI2sre7321mzwB93P3tRAUSTvr3e6AcONPdGxK1bxERSX/7TUhmlg/UmFmxu+8GcPdVSYjlNwQDt57i7ruSsH8REUlj+71k5+6NwBJgYLKCMLPhwFeAScAGM6sJH+cn6zNFRCS9RL1kdy/wWDiW3RrAYyvcfU5ng3D3lcTdbCsiIrknakK6NPx5bYvlDhySsGhERCRnRR3L7uBkByIiIrktatk3ZtbdzI4zs3PC133MrE/yQhMRkVwSKSGZ2QcIChvuICjNhmDcuUSMYyciIhK5hfQb4EfuPhaI3R/0LHtvlhUREemUqAlpPMFo3xBW2Ll7LZq2XEREEiRqQlpBi5G3zexIgnmLREREOi1q2fcPgcfN7Hagh5n9F/BV4MtJi0xERHJKpBZSOD/RGUApQd/RcOBT7v7PJMYmIiI5JOqMsSXu/hpwWZLjERGRHBW1D2mVmT1hZufr3iMREUmGyDPGAo8RDCG0wczuN7OzzCxqH5SIiEi7ovYhbXL3X7v7sQQl4AuAG4H1yQxORERyR+Shg+KUh48SYFtCoxERkZwVdeigw8zsejN7F3gkXHy2u49OWmQiIpJTovYBvQA8BEwD5ri7A5hZXmxqcxERkc6ImpDK3b0+9iIcbPUC4HPAQckITEREckvUooZ6Mys1s6+b2WvAfGAK8PVkBiciIrmj3RaSmXUHPgFcCEwlGLvufoKRGv7T3SuTHaCIiOSG/bWQNgK/Bd4BPuzuh7n79UB9+28TERHpmP0lpIVAMXAUcISZ9U96RCIikpPaTUjufgIwEvgn8C2CURoeBfoA3ZMenYiI5Iz9FjW4+0p3vz685+hkgtEZmoAFZvbzZAcoIiK5oUMjNbj78+4+DRgEXAF8IClRiYhIzjmQoYNw9zp3v9/dz0h0QCIikpsOKCGJiIgkmhKSiIikBSUkERFJC0pIIiKSFpSQREQkLSghiYhIWlBCEhGRtKCEJCIiaUEJSURE0kLaJCQzG2BmfzWzWjNbaWafS3VMIiLSdaJOYd4VbiOYZ6kcmAQ8bmYL3H1RSqMSEZEukRYtJDPrA3wa+KG717j788DfgS+kNjIREekqaZGQgDFAo7sviVu2ABifonhERKSLpcslu0Jge4tl24Gi+AVmNg2YFr6sMbN3Ovm5JcCmTu4jnWXz8WXzsUF2H182Hxvo+PZneFsr0iUh1QB9WyzrC+yIX+DuM4GZifpQM5vr7lMStb90k83Hl83HBtl9fNl8bKDj64x0uWS3BOhmZqPjlk0EVNAgIpIj0iIhuXst8DBwnZn1MbNjgE8C96Q2MhER6SppkZBClwG9gErgfuDSLij5TtjlvzSVzceXzccG2X182XxsoOM7YObuydq3iIhIZOnUQhIRkRymhCQiImkh6xJSR8bEM7NvmNkGM9tuZneaWcGB7KerJPDYnjGzOjOrCR+dvZ8rIaIen5lNMLOnzGyTmb3vmnMmn7sIx5bp5+4CM5tnZtVmtsbMfm5m3Tq6n66UwGPL9HN3rpm9E36nVJrZLDPr29H9tMvds+pBUBDxJ4KbbY8luMF2fCvbTQU2EowG0R94Bvifju4nQ4/tGeBLqT5XnTi+Q4FLCCox/UD3k6HHlunn7lLgOKAHMASYB3w3S87d/o4t08/dMKAkfF4I3AvMSOS5S/kvI8G/2D4EA7SOiVt2T/yXcdzy+4CfxL0+GdjQ0f1k2rGFr9PuP8aB/M6BUS2/tDP93LV3bNl07uK2uxp4NJvOXWvHlm3nLkw6dwNPJPLcZdslu46MiTc+XBe/XbmZDezgfrpKoo4t5qfhZaEXzOyERAd7ABL1O8/0cxdFNp2749l7A3y2nbv4Y4vJ6HNnZsea2XaCUXQ+DUw/kP20JdsSUqQx8drYNva8qIP76SqJOjaAa4BDCC4rzAQeNbORiQv1gCTqd57p525/subcmdlFwBTgl53ZT5Il6tggC86duz/v7v2AocAvgBUHsp+2ZFtCijQmXhvbxp7v6OB+ukqijg13f9ndd7j7bnefBbwAnJngeDsqUb/zTD937cqWc2dmZwP/A5zh7rGBOrPi3LVxbFlz7gDcfS3wJPBAZ/bTUrYlpI6MibcoXBe/3UZ339zB/XSVRB1baxywhER54BL1O8/0c9dRGXfuzOx04A7gLHd/40D300USdWytybhz10I3INbCS8y5S3WnWhI66R4gqPboAxxD2xUjpwMbgMMIKtHmsG8lWqT9ZNqxAcUEVXg9w39Q5wO1wKEZdO4sjP8wgv/UPYGCLDl3bR5blpy7k4DNwPGd2U+mHVuWnLvzgYrw3+hw4Fng4USeu5T+IpL0yx0APBKe7FXA58LlFQTNyoq4ba8mKI+uBv7Q4kut1f1k+rEBpcCrBE3pbcBLwKmpPraOHB8wguDLOv6xIhvOXXvHliXn7mlgT7gs9vhHlpy7No8tS87djcCacLs1BP1gAxN57jSWnYiIpIVs60MSEZEMpYQkIiJpQQlJRETSghKSiIikBSUkERFJC0pIIiKSFpSQROR9wvl6Dkl1HJJblJAkrZnZCjPbFX5BbjCzu8ysMOJ7nzGzLyU7xrjPG2FmHj8pW6Zy90J3X57qOCS3KCFJJjjL3QuBScBk4L+64kOzIbF0VC4es6QPJSTJGO6+AXiKIDEBYGYfNrMXzWybmS2IzTFjZjcSzN55a9i6urW1Fkx8K8rMLgznqbnZzLYA14YtstvM7HEz22FmL7czZcBz4c9t4WcebWZ5ZvaDcErnSjO728z6tXWMZvZxM5sfHs+LZnZ4uPwcM1semzLazM4IW4yl4Ws3syvDbTaZ2S/MLC9uvxeb2dtmttWCKdKHx61zM/uamS0FlsYtGxU+LzCzX5rZKjPbaGa3m1mvcN0JFkzX/c3w+NaHUy/E9t3LzP43PP7tZvZ83HtbPXeSw1I9jpIeerT3IJhv5ZTw+VDgDeCW8PUQgsEszyT44+rU8HVpuP4Z4mboZO84cd3iljVvA1xIMBbZFQQDYPYC7gK2AEeGy+4FHmgj1tb2fzGwjGAenELgYeCeNt7/QaASOArIBy4Ijz82DuG9YTwDgXXAx+Pe6wRjqQ0gGINsSdxxnR3GMC48hh8AL7Z477/C9/aKWzYqfD4d+Hu4vgh4FPhpuO6E8Hd2HdA9PBc7gf7h+tvC3/GQ8Jg+AhTs79zpkZuPlAeghx7tPcIv5BqCQSkdmA0Uh+uuafnlTtCCuiB8fiAJaVWL/d0F/C7u9ZnA4jZibW3/s4HL4l4fCjTEbxO37jfA9S2WvQN8NHxeTDBo5RvAb1ts58Dpca8vA2aHz/8BXBK3Li9MGsPj3ntSK/sbRTCycy0wMm7d0cB74fMTgF0tjrkS+HD4ObuAia0ca7vnTo/cfOiSnWSCs929iODLbyxQEi4fDnw2vOSzzcy2AccCgzvxWatbWbYh7vlOgpZOVAcBK+NeryRopZS3su1w4JstjmdYuA/cfRvwF2AC8L/7iX1l7H3hfm+J2+cWgkQzpI33xisFegPz4t7/ZLg8ZrO774l7HfsdlRBMt/BuG8ea6HMnGU4dmJIx3P1ZM7uLYFroswm+RO9x9y+39ZYWr2vDn70JpuUAGLSf93QoxFaWrSP48o2pILjEtbGVbVcDN7r7ja3t3MwmEVwCvB+YQTDvVbxh7J0QrSL87Pj93tvB2AE2EbRyxnswS2hHbALqCCZxW9Bi3f7OneQgtZAk00wHTg2/nP8InGVmU80s38x6hp3sQ8NtNxL03QDg7lXAWuDz4fYXs3fGy0SoApriP5MgeXzDzA4Oy9V/AvypRYsi5g7gq2Z2lAX6mNnHzKzIzHqGx/s94CJgiJld1uL93zaz/mY2DPg68Kdw+e3Af5nZeAAz62dmn41yQO7eFMZ1s5mVhe8fYmZTI773TuAmMzso/J0fbWYF7P/cSQ5SQpKMEiaVu4Efuvtq4JMEX9JVBH91f5u9/65vAT4TVpbNCJd9OdxmMzAeeDGBse0kmMTshfAy1IcJvpDvIajAe4+gxXBFG++fG8Z3K7CVoBDhwnD1T4E17v4bd98NfB64wfadMvpvwDxgPvA48Ptwv38FfgY8YGbVwJvAGR04tGvCWF4K3/9vgr6wKL5F0Of1KsGlwp8BeRHOneQgTdAnkgXMzIHR7r4s1bGIHCj9NSIiImlBCUlERNKCLtmJiEhaUAtJRETSghKSiIikBSUkERFJC0pIIiKSFpSQREQkLSghiYhIWvj/exqYCOliCR4AAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "fig, ax = plt.subplots()\n", "\n", "ax.plot(grid_points, rslts)\n", "\n", "ax.set_ylim([0, 10])\n", "ax.set_xlabel(\"Return to experience\")\n", "ax.set_ylabel(\"Average final level of exerience\")\n", "\n", "plt.show()\n", "plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the absence of any returns to experience, Robinson still spends more than two periods fishing. This share then increases with the return. Starting at around 0.2, Robinson spends all his time fishing.\n", "\n", "## Extension\n", "\n", "Let us make the model more interesting!\n", "\n", "---\n", "*At some point Crusoe notices that a group of cannibals occasionally visits the island and celebrate one of their dark rituals. But then, a prisoner can escape and becomes Crusoe's new friend Friday whom he teaches English. In return Friday can share his knowledge once to help Robinson improve his fishing skills, but that is only possible after Robinson tried at least once to go fishing.*\n", "\n", "---\n", "\n", "A common extension to structural models is to increase the choice set. Here, we want to add another choice called `\"friday\"` which affects the utility of fishing. The choice should be available once, starting with the third period, and only after Robinson has been fishing before.\n", "\n", "Note that, we load the example models with the function, `rp.get_example_model`. The command for the former model is `params, options, df = rp.get_example_model(\"robinson_crusoe_basic\")`. You can use `with_data=False` to suppress the automatic simulation of a sample with this parameterization." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "params, options = rp.get_example_model(\"robinson_crusoe_extended\", with_data=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At first, take a look at the parameterization. There is a new positive parameter called `\"contemplation_with_friday\"` which enters the wage equation of fishing. The choice `\"friday\"` itself has a negative constant utility term which models the effort costs of learning and the food penalty. The variance-covariance matrix is also adjusted." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ " value\n", "category name \n", "delta delta 0.95\n", "wage_fishing exp_fishing 0.10\n", " contemplation_with_friday 0.40\n", "nonpec_fishing constant -1.00\n", "nonpec_friday constant -1.00\n", " not_fishing_last_period -1.00\n", "nonpec_hammock constant 2.50\n", " not_fishing_last_period -1.00\n", "shocks_sdcorr sd_fishing 1.00\n", " sd_friday 1.00\n", " sd_hammock 1.00\n", " corr_friday_fishing 0.00\n", " corr_hammock_fishing 0.00\n", " corr_hammock_friday 0.00\n", "lagged_choice_1_hammock constant 1.00" ], "text/html": "
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
value
categoryname
deltadelta0.95
wage_fishingexp_fishing0.10
contemplation_with_friday0.40
nonpec_fishingconstant-1.00
nonpec_fridayconstant-1.00
not_fishing_last_period-1.00
nonpec_hammockconstant2.50
not_fishing_last_period-1.00
shocks_sdcorrsd_fishing1.00
sd_friday1.00
sd_hammock1.00
corr_friday_fishing0.00
corr_hammock_fishing0.00
corr_hammock_friday0.00
lagged_choice_1_hammockconstant1.00
\n
" }, "metadata": {}, "execution_count": 13 } ], "source": [ "params" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Turning to the `options`, we can see that the new covariate `\"contemplation_with_friday\"` is only affecting utility if Robinson is experienced in fishing and only for one interaction with friday. This naturally limits the interaction with Friday. The key `\"negative_choice_set\"` can be used to restrict the choice Friday to the third and following periods. The first key matches a choice. The value of the key can be a list of strings. If the string evaluates to `True`, a utility penalty ensures that individuals will never choose the corresponding states. There exist some states in the state space which will never be reached because choices are mutually exclusive or are affected by other restrictions. Filters under `\"core_state_space_filters\"` can be used to purge those states from the state space, reducing runtime and memory consumption." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "{'n_periods': 10,\n", " 'estimation_draws': 200,\n", " 'estimation_seed': 500,\n", " 'estimation_tau': 0.001,\n", " 'interpolation_points': -1,\n", " 'simulation_agents': 1000,\n", " 'simulation_seed': 132,\n", " 'solution_draws': 500,\n", " 'solution_seed': 456,\n", " 'covariates': {'constant': '1',\n", " 'contemplation_with_friday': 'exp_friday == 1 and exp_fishing >= 1',\n", " 'not_fishing_last_period': \"lagged_choice_1 != 'fishing'\"},\n", " 'negative_choice_set': {'friday': ['period < 2', 'exp_fishing == 0']},\n", " 'core_state_space_filters': [\"period > 0 and exp_fishing + exp_friday == period and lagged_choice_1 == 'hammock'\",\n", " 'period <= 2 and exp_friday != 0',\n", " 'period >= 3 and period - exp_friday < 2',\n", " 'exp_friday > 0 and exp_fishing == 0',\n", " \"exp_friday > 0 and exp_fishing == 1 and lagged_choice_1 == 'fishing'\",\n", " \"period - exp_friday == 2 and lagged_choice_1 != 'friday' and period > 2\",\n", " \"exp_{choices_w_exp} == 0 and lagged_choice_1 == '{choices_w_exp}'\"]}" ] }, "metadata": {}, "execution_count": 14 } ], "source": [ "options" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let us simulate a sample of the new model." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "simulate = rp.get_simulate_func(params, options)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "df = simulate(params)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2021-06-01T13:22:59.789117\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.2, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEVCAYAAACv2pHlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZTklEQVR4nO3de7hddX3n8fcHgggJQZA0ViUGFRxFCx3T0ani3UFHrXQYfVREtFVaUJ/Rqh1qdUDAehlrHRVxmBFBERRHVMRLvSBTcKwa64DNWCNouAfCxZAEAorf+WOt4M7hJNmH7JP9y8779Tznydlr/fba3+zL+ezfWr/1W6kqJEkat53GXYAkSWAgSZIaYSBJkppgIEmSmmAgSZKaYCBJkpowZ9wFAOyzzz61ePHicZchSZplP/zhD2+qqgXTrWsikBYvXszSpUvHXYYkaZYluXJT69xlJ0lqgoEkSWqCgSRJaoKBJElqwlCBlOR1SZYmuTPJGVto+8YkK5OsTnJ6kl1HUqkkaaIN20O6DjgZOH1zjZIcChwHPBNYDDwceMdW1CdJ2kEMFUhVdV5VfQG4eQtNjwI+VlXLqupW4CTglVtVoSRphzDqY0gHApcO3L4UWJjkgSN+HEnShBn1ibHzgNUDtzf8vgdTeldJjgaOBli0aNFwWz9hz60u8LfbWr3lNkNvy7pmti3rmtm2rGtm27KumW2rnbpG3UNaC8wfuL3h9zVTG1bVaVW1pKqWLFgw7SwSkqQdyKgDaRlw0MDtg4AbqmpLx54kSTu4YYd9z0lyf2BnYOck908y3e6+TwB/muQxSfYC3gacMbJqJUkTa9ge0tuAO+iGdL+8//1tSRYlWZtkEUBVfQ14L/Bt4Mr+5/iRVy1JmjhDDWqoqhOAEzaxet6Utu8H3r9VVUmSdjhOHSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWqCgSRJaoKBJElqgoEkSWrCUIGUZO8kn0+yLsmVSV62iXZJcnKSa5OsTnJRkgNHW7IkaRIN20M6BbgLWAgcAZy6iaB5EfAnwCHA3sB3gU+OoE5J0oTbYiAlmQscDry9qtZW1SXA+cCR0zTfD7ikqn5eVXcDZwGPGWXBkqTJNEwP6QDg7qpaPrDsUmC6HtKngUcmOSDJLsBRwNe2vkxJ0qSbM0SbecDqKctWA3tM0/Z64GLgp8DdwNXAM6bbaJKjgaMBFi1aNGS5kqRJNUwPaS0wf8qy+cCaadoeD/wBsC9wf+AdwIVJdp/asKpOq6olVbVkwYIFM6takjRxhukhLQfmJNm/qn7WLzsIWDZN24OAz1TVNf3tM5J8gO440tKtLVYz87j9Rtfz/PHItiRJ09tiD6mq1gHnAScmmZvkScALmX703A+AFyVZmGSnJEcCuwCXj7JoSdLkGaaHBHAscDpwI3AzcExVLUuyCPh/wGOq6irgPcDvAP8XmEsXRIdX1S9HXLckacIMFUhVdQtw2DTLr6Ib9LDh9nrgtf2PJElDc+ogSVITht1lJ42Mgy0kTccekiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQnb1eSqi9efPbJtrRjZliRJo2APSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktSE7WouO2k2PW6/RSPb1o9HtiVpx2EPSZLUBANJktQEA0mS1ASPIUmN89iWdhRD9ZCS7J3k80nWJbkyycs20/bhSS5IsibJTUneO7pyJUmTathddqcAdwELgSOAU5McOLVRkvsB3wAuBB4EPBQ4azSlSpIm2RYDKclc4HDg7VW1tqouAc4Hjpym+SuB66rq/VW1rqrWV9VlI61YkjSRhukhHQDcXVXLB5ZdCtyrhwQ8EViR5Kv97rqLkjxuuo0mOTrJ0iRLV61aNfPKJUkTZZhAmgesnrJsNbDHNG0fCrwE+CDwYODLwBf7XXkbqarTqmpJVS1ZsGDBzKqWJE2cYUbZrQXmT1k2H1gzTds7gEuq6qsASd4HvA14NF2vaiItXn/2yLa1YmRbkqTtyzA9pOXAnCT7Dyw7CFg2TdvLgBpFYZKkHcsWA6mq1gHnAScmmZvkScALgU9O0/ws4IlJnpVkZ+ANwE3AT0ZXsiRpEg077PtYYDfgRuAc4JiqWpZkUZK1SRYBVNVPgZcDHwVupQuuP6qqu0ZfuiRpkgw1U0NV3QIcNs3yq+gGPQwuO4+uRyVpgjmDhEbNqYOkxq35ybvHXYK0TRhIkiaKPbftl4Ek6T6x56ZRM5Cknn9gJ4Ov4/bL6yFJkppgIEmSmmAgSZKa4DEkSdoGnPNyywwkbXMedJY0HXfZSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmjBUICXZO8nnk6xLcmWSlw1xnwuTVJI5W1+mJGnSDRsWpwB3AQuBg4EvJ7m0qpZN1zjJETPYtiRJWw6NJHOBw4HHVtVa4JIk5wNHAsdN035P4HjgFcB3R1uuJGmUFq8/e2TbWrGV9x9ml90BwN1VtXxg2aXAgZto/zfAqcDKraxNkrQDGSaQ5gGrpyxbDewxtWGSJcCTgA9taaNJjk6yNMnSVatWDVOrJGmCDRNIa4H5U5bNB9YMLkiyE/AR4D9V1a+3tNGqOq2qllTVkgULFgxbryRpQg0TSMuBOUn2H1h2EDB1QMN8YAnwmSQrgR/0y69JcshWVypJmmhbHNRQVeuSnAecmOTVdKPsXgj84ZSmq4EHD9zeF/g+8HjAfXKSpM0a9sTYY4HdgBuBc4BjqmpZkkVJ1iZZVJ2VG374bQjdUFV3zULtkqQJMtS5QlV1C3DYNMuvohv0MN19VgDZitokSTsQpw6SJDXBQJIkNcHpfSbYmp+8e9wlSNLQ7CFJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaYCBJkppgIEmSmmAgSZKaMFQgJdk7yeeTrEtyZZKXbaLdUUl+mOS2JNckeW+SOaMtWZI0iYbtIZ0C3AUsBI4ATk1y4DTtdgfeAOwDPAF4JvDmrS9TkjTptth7STIXOBx4bFWtBS5Jcj5wJHDcYNuqOnXg5rVJPgU8fYT1SpIm1DA9pAOAu6tq+cCyS4HpekhTPQVYdl8KkyTtWIYJpHnA6inLVgN7bO5OSV4FLAHet4n1RydZmmTpqlWrhqlVkjTBhgmktcD8KcvmA2s2dYckhwHvBp5bVTdN16aqTquqJVW1ZMGCBUOWK0maVMME0nJgTpL9B5YdxCZ2xSV5DvA/gBdU1Y+3vkRJ0o5gi4FUVeuA84ATk8xN8iTghcAnp7ZN8gzgU8DhVfX9URcrSZpcww77PhbYDbgROAc4pqqWJVmUZG2SRX27twN7Al/pl69N8tXRly1JmjRDnbRaVbcAh02z/Cq6QQ8bbjvEW5J0nzh1kCSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQkGkiSpCQaSJKkJBpIkqQlDBVKSvZN8Psm6JFcmedlm2r4xycokq5OcnmTX0ZUrSZpUw/aQTgHuAhYCRwCnJjlwaqMkhwLHAc8EFgMPB94xkkolSRNti4GUZC5wOPD2qlpbVZcA5wNHTtP8KOBjVbWsqm4FTgJeOcJ6JUkTapge0gHA3VW1fGDZpcC9ekj9skuntFuY5IH3vURJ0o4gVbX5BskhwGer6kEDy14DHFFVT5vS9grgtVX1tf72LnS7+varqhVT2h4NHN3ffBTw0636n/zWPsBNI9rWKFnXzFjXzFjXzFjXzIyyrodV1YLpVswZ4s5rgflTls0H1gzRdsPv92pbVacBpw3x+DOSZGlVLRn1dreWdc2Mdc2Mdc2Mdc3MtqprmF12y4E5SfYfWHYQsGyatsv6dYPtbqiqm+97iZKkHcEWA6mq1gHnAScmmZvkScALgU9O0/wTwJ8meUySvYC3AWeMsF5J0oQadtj3scBuwI3AOcAxVbUsyaIka5MsAuiPHb0X+DZwZf9z/OjL3qyR7wYcEeuaGeuaGeuaGeuamW1S1xYHNUiStC04dZAkqQkGkiSpCRMTSDOZb29bSvK6JEuT3JnkjHHXA5Bk1yQf65+nNUl+lOS5464LIMlZSa5PcluS5UlePe6aBiXZP8n6JGeNuxaAJBf19aztf0Z1Pt9WS/KSJD/pP5NX9Oc0jrOetVN+7k7yoXHWtEGSxUm+kuTWfi7QDycZ5rSc2a7r0Uku7OcmvTzJH8/m401MIDHkfHtjcB1wMnD6uAsZMAe4GngqsCfwduDcJIvHWVTvXcDiqpoP/BFwcpLHj7mmQacAPxh3EVO8rqrm9T+PGncxAEmeDbwHeBWwB/AU4OfjrGngOZpH93fiDuCz46xpwEfoBo39LnAw3Wfz2HEW1AfiF4ELgL3pJjI4K8kBs/WYExFIM5xvb5uqqvOq6gtAM+diVdW6qjqhqlZU1W+q6gLgF8DY//D38yDeueFm//OIMZZ0jyQvAX4JfGvMpWwP3gGcWFX/2L/Hrq2qa8dd1ID/SBcAF4+7kN5+wLlVtb6qVgJfY/rp2balfwU8GPi7qrq7qi4EvsMs/l2diEBiZvPtaYokC+mew+lOdt7mknwkye3AvwDXA18Zc0kkmQ+cCLxp3LVM411JbkrynSRPG3cxSXYGlgAL+t081/S7oHYbd20DjgI+Ue0MM/5vwEuS7J7kIcBz6UJpnLKJZY+drQeclECaB6yesmw13a4CbUY/3+CngDOr6l/GXQ9AVR1L99odQndS9p2bv8c2cRLdTPZXj7uQKf4z3WVeHkJ3rsiXkoy7R7kQ2IWuF3II3S6o36c7UX7s+vMmnwqcOe5aBvxvui/QtwHXAEuBL4yzILovhDcCb0myS5J/R/e87T5bDzgpgTST+fbUS7IT3YwbdwGvG3M5G+l3EVwCPBQ4Zpy1JDkYeBbwd+OsYzpV9b2qWlNVd1bVmXS7VP79mMu6o//3Q1V1fVXdBLyf8de1wSuAS6rqF+MuBO75HP493ZevuXQTme5FdwxubKrqV8BhwPOAlXR7B86lC8xZMSmBNJP59gQkCfAxum+zh/dvvhbNYfzHkJ5Gd8HJq5KsBN4MHJ7kn8ZZ1CYU0+9q2XYFdNdCu6avpUWvoK3e0d7AvsCH+y8WNwMfp4EAr6rLquqpVfXAqjqUrjf+/dl6vIkIpBnOt7dNJZmT5P7AzsDOSe7fwnBO4FTg0cALquqOLTXeFpL8Tj9UeF6SnfsrEL8UuHDMpZ1GF4oH9z8fBb4MHDq+kiDJA5IcuuE9leQIutFsfz/OunofB17fv6Z7AW+gG601Vkn+kG73Ziuj6+h7kL8AjulfxwfQHeO6dLN33AaS/F7//to9yZvpRgGeMVuPNxGB1Jt2vr3xlgR0+83voLu0+8v738e6Lz3Jw4A/o/vjunLgvIwjxlkX3TfqY+i+Xd8KvA94Q1V9caxFVd1eVSs3/NDtIl5fVavGWRfdcZqTgVV016p5PXBYVbVwLtJJdMPjlwM/AX4EvHOsFXWOAs6rqtZ25/8H4Dl0r+XlwK+BN461os6RdAOLbgSeCTx7YBTsyDmXnSSpCZPUQ5IkbccMJElSEwwkSVITDCRJUhMMJElSEwwkSVITDCRpjJJ8NclR9/G+K5I8a9Q1SePSwowB0nYnyQq6aZfuBtbRzUj++qpaO5PtVFUTF0aUWmAPSbrvXtBf7O1fA3/ADGbgSMfPnzTAD4S0lfoLz30VeGySJyb5P0l+meTSwesT9Zcbf2eS7wC3Aw/vl726X79Tkrf1l5a/Mcknkuw5cP8j+3U3J/nrbfu/lGafgSRtpST70s3MfD3dpKsn083g/Gbgc0kWDDQ/ku5S0HsAV07Z1Cv7n6fTzao8D/hw/xiPoZsQ90i6q3g+kO7SHNLEMJCk++4LSX4JXEJ3gbVrgK9U1Vf6y3Z/g+5Ca4OXETijv0z7r6e55McRwPur6uf9sai/oruK6By6i91dUFX/0E9u+XbgN7P735O2LQNJuu8Oq6oHVNXD+qvcLgRe1O+u+2UfVk+mm7J/g81dcfbBbNxrupJu4NHCft099+0vuXLzaP4bUhscZSeNztXAJ6vqNZtps7np9a8DHjZwexHdZQhuoNsd+OgNK5LsTrfbTpoY9pCk0TkLeEF/0bwNF2N8WpJhj/WcA7wxyX5J5gF/A3ymqn4N/C/g+UmenOR+wIn4+dWE8Q0tjUhVXU13peK30l1o7WrgLQz/OTud7irH/0B3BdH1dBfdo7/Y5GuBs+l6SxsuEy5NDC/QN2GSPAr4NPBIYC5wfFWdtIX7rABeXVXfnGbdIcD/rKpHzUK5E2vK6/DXVfXBKes/Cly7qdcmSQH7V9Xls17shNjc+3h7k+QE4JFV9fJx17IteQxp8vwlcFFV/f4oNlZVFwOG0cxt9nWoqj/fxvVIzTOQRmTxcV+e1a7minc/L0M2fRjdN/MdzuPOfNysvgY/PurHw74GsJnXIcnOVXX3aKoasxP2nN1dLCesnslzru2cx5AmSJIL6U6q/HCStUnOTnJyv26fJBf0w5FvSXLxlKlrDk5yWZLVST6T5P79/Z6W5JqBx1iR5M3Tte3X/2WS65Ncl+TVSSrJI7fRU9CETbwOpyb5SpJ1wNOTnLHhtenv85aB5+1PpmzveUl+lOS2JFf3u3M2rPtyktdPaX9ZksNm9T/Zrnu9j5Ps1b/3VyW5tf/9noEm/WwZJ/czbKxN8qUkD0zyqf45/0GSxQPtK8mxSX6WZE2Sk5I8Isl3+/bn9gNPNrR/TZLL+8/d+UkePLDuwCTf6NfdkOStU/9DSXZJck6Szw1udxIZSBOkqp4BXAy8rp9j7a6B1W+iOwi+gO68lrey8RDkFwPPAfYDfo9uxoBNmbZtkucAfwE8i+7YyVO38r+0XdrE6/Ay4J10MzRcMti+f97eDDwb2J/u+Ru0DngF8ADgecAxA4FzJnDPcYYkBwEPoZvsdUc03XtzJ+DjdL3WRcAd9DNgDHgJ3SwYDwEeAXy3v8/ewE+A46e0fw7weOCJdLtnT6M7sXlf4LHASwGSPAN4V1/X79KdW/bpft0ewDeBr9GdZ/ZI4FuDD5JkN+ALwJ3Ai6tq8DM9cQykHcev6D4QD6uqX1XVxbXxiJYPVtV1VXUL8CXg4M1sa1NtXwx8vJ+J4HbgHSP/X2y/vlhV3+lncFg/Zd2G5+2f+xNeTxhcWVUXVdWP+/teRjc8fEPYfxHYP8n+/e0j6YaKT/Qfrs2413uzqm6uqs9V1e1VtYbui8HUL0sfr6orqmo13byEV1TVN/sh958Fph4LfE9V3daPfvxn4Ov9DBsb7r+h/RHA6VX1T/0MG38F/Nu+x/V8YGVV/W1Vra+qNVX1vYHHmE8XVlcAr5qY3bybYSDtOP4rcDnw9SQ/T3LclPUrB36/nW4etU3ZVNuNZhNg87MS7Gi2NEPD4PqN5rhL8oQk3+53Oa0G/hzYB6D/I3cu8PJ+F+xL6YaO76ju9d5MsnuS/55uYtrb6IbVPyDJzgNtbxj4/Y5pbk/9PAzbfqPZN/opoW6m64ntSxc2m/JEul7eu6d8eZxYBtIOov/29aaqejjwAuAvkjxzxA9zPRtP+LnviLe/PdvcH5Tr2fi5WjRl/dnA+cC+VbUn8FFg8GD/mXTfxJ8J3F5V3936cifKm+hGij6hquYDT+mXb4sBExvNvpFkLt0MG9fSfQl5xGbu+3W63X3fSrJwNotshYG0g0jy/CSPTBLgNroLy416F8C5wKuSPDrd1Db/ZcTbn1TnAq9M8pj+eZt6vGIP4JaqWp/k39Adj7pHH0C/Af6WHbt3tCl70PVafplkb+79/M6ms+k+Ewcn2ZVu9o3vVdUK4ALgQUnekGTXJHskecLgnavqvf02vpVkn21Y91gYSDuO/ekOoK6lO2D7kaq6aJQPUFVfBT4IfJtu9+CGb+p3jvJxJk3/vH0AuJDuebtwSpNjgROTrKEL+XOn2cwngMfRTV+kjX0A2A24CfhHuuMy20RVfYtuZvbP0fWEH0E3gIL+eNaz6fZYrAR+Rjc6c+o2TqIb2PDNPlAnljM1aNYkeTTdAd9d+4PDmiVJXgEcXVVPHnct0n1lD0kjleSPk9wvyV7Ae4AvGUazq9/Ndyzd0GNpu2UgadT+jG5i0SvojlEdM95yJluSQ+me7xvojjVI2y132UmSmmAPSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktQEA0mS1AQDSZLUBANJktSE/w+12FLV9SxuLwAAAABJRU5ErkJggg==\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "fig, ax = plt.subplots()\n", "\n", "df.groupby(\"Period\").Choice.value_counts(normalize=True).unstack().plot.bar(\n", " stacked=True, ax=ax, color=[\"C0\", \"C2\", \"C1\"],\n", ")\n", "\n", "plt.xticks(rotation=\"horizontal\")\n", "\n", "plt.legend(loc=\"lower center\", bbox_to_anchor=(0.5,-0.275), ncol=3)\n", "\n", "plt.show()\n", "plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimation\n", "\n", "For model calibration, **respy** supports estimation via maximum likelihood and the method of simulated moments. An appropriate criterion function function for both methods can be constructed in just a few steps.\n", "\n", "For optimization of the criterion, the use of external optimization libraries is requried. We recommend [estimagic](https://github.com/OpenSourceEconomics/estimagic), an open-source tool to estimate structural models and more. **estimagic** can be used for the optimization and standard error calculation of a criterion function produced by **respy**. \n", "\n", "Unlike other optimization libraries, ``estimagic`` does not optimize over a simple vector of parameters, but instead stores parameters in a ``pd.DataFrame``, which makes it easier to parse them into the quantities we need, store lower and upper bounds together with parameters and express constraints on the parameters. \n", "\n", "For ``estimagic``, we need to pass constraints on the parameters in a list containing dictionaries. Each dictionary is a constraint. A constraint includes two components: First, we need to tell ``estimagic`` which parameters we want to constrain. This is achieved by specifying an index location which will be passed to `df.loc`. Then, define the type of the constraint. Here, we only impose the constraint that the shock parameters have to be valid variances and correlations.\n", "\n", "*Note*: It is possible to utilize other optimization libraries but we recommend **estimagic** due to the reasons stated above." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "from estimagic import maximize" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "crit_func = rp.get_log_like_func(params, options, df)\n", "crit_func(params)\n", "\n", "constr = rp.get_parameter_constraints(\"robinson_crusoe\")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "results = maximize(\n", " criterion=crit_func, \n", " params=params, \n", " algorithm=\"scipy_lbfgsb\", \n", " algo_options={\"stopping_max_criterion_evaluations\": 3}, \n", " constraints=constr,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Running the minimal working example shown above will start an optimization that is limited to three criterion evaluations (for the sake of this tutorial). **estimagic** will also produce a logging database called `logging.db` that stores information about the optimization. The package offers many useful options to set up a tractable estimation for your model.\n", "\n", "More information can be found in the **estimagic** documentation: https://estimagic.readthedocs.io/." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Footnotes\n", "\n", "1\n", " One of the earliest references of Robinsonades in economics can be found in Marx (1867). In the 37th footnote, he mentions that even Ricardo used the theme before him.\n", "\n", "\n", "## References\n", "\n", "> Bellman, R. (1957). Dynamic Programming. *Princeton University Press*, Princeton, NJ.\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.\n", "\n", "> Marx, K. (1867). Das Kapital, Bd. 1. *MEW*, Bd, 23, 405" ] } ], "metadata": { "kernelspec": { "name": "python3710jvsc74a57bd0b549dadab0c2edb1c58f223f7584f57e60f2cc8e65ef8392efb9b23cb30dad20", "display_name": "Python 3.7.10 64-bit ('respy': conda)" }, "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.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }