Skip to content

Instantly share code, notes, and snippets.

@justinian336

justinian336/source-ddc-showcase.ipynb Secret

Created September 2, 2020 06:12
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save justinian336/778493f3d5314a346ac08a61df7a10ba to your computer and use it in GitHub Desktop.
Save justinian336/778493f3d5314a346ac08a61df7a10ba to your computer and use it in GitHub Desktop.
Download ZIP
Source DDC Showcase
Display the source blob
Display the rendered blob
Raw
source-ddc-showcase.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Source DDC Showcase",
"provenance": [],
"collapsed_sections": [],
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/justinian336/778493f3d5314a346ac08a61df7a10ba/source-ddc-showcase.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "rRO3awGLVpOS",
"colab_type": "text"
},
"source": [
"## Source DDC Showcase\n",
"\n",
"Here we show a basic example of how you can use Source DDC to simulate and perform structural estimation on a simple dynamic discrete choice model.\n",
"\n",
"The model is a very simplified version of the retirement model in Rust & Phelan (1997).\n",
"\n",
"First, install the library. You might be asked to restart your kernel afterwards to use the correct dependencies."
]
},
{
"cell_type": "code",
"metadata": {
"id": "6gsRD4v_VqaD",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 343
},
"outputId": "6d4e72f1-0118-4b31-9d9e-de516d30dc6b"
},
"source": [
"!pip install git+https://github.com/sansan-inc/econ-source.git#subdirectory=source_ddc"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Collecting git+https://github.com/sansan-inc/econ-source.git#subdirectory=source_ddc\n",
" Cloning https://github.com/sansan-inc/econ-source.git to /tmp/pip-req-build-b1wm_95s\n",
" Running command git clone -q https://github.com/sansan-inc/econ-source.git /tmp/pip-req-build-b1wm_95s\n",
"Requirement already satisfied (use --upgrade to upgrade): source-ddc==0.0.2.dev1 from git+https://github.com/sansan-inc/econ-source.git#subdirectory=source_ddc in /usr/local/lib/python3.6/dist-packages\n",
"Requirement already satisfied: numpy==1.18.4 in /usr/local/lib/python3.6/dist-packages (from source-ddc==0.0.2.dev1) (1.18.4)\n",
"Requirement already satisfied: pandas==1.0.3 in /usr/local/lib/python3.6/dist-packages (from source-ddc==0.0.2.dev1) (1.0.3)\n",
"Requirement already satisfied: scipy==1.4.1 in /usr/local/lib/python3.6/dist-packages (from source-ddc==0.0.2.dev1) (1.4.1)\n",
"Requirement already satisfied: statsmodels==0.12.0rc0 in /usr/local/lib/python3.6/dist-packages (from source-ddc==0.0.2.dev1) (0.12.0rc0)\n",
"Requirement already satisfied: bidict==0.20.0 in /usr/local/lib/python3.6/dist-packages (from source-ddc==0.0.2.dev1) (0.20.0)\n",
"Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas==1.0.3->source-ddc==0.0.2.dev1) (2018.9)\n",
"Requirement already satisfied: python-dateutil>=2.6.1 in /usr/local/lib/python3.6/dist-packages (from pandas==1.0.3->source-ddc==0.0.2.dev1) (2.8.1)\n",
"Requirement already satisfied: patsy>=0.5 in /usr/local/lib/python3.6/dist-packages (from statsmodels==0.12.0rc0->source-ddc==0.0.2.dev1) (0.5.1)\n",
"Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.6.1->pandas==1.0.3->source-ddc==0.0.2.dev1) (1.15.0)\n",
"Building wheels for collected packages: source-ddc\n",
" Building wheel for source-ddc (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
" Created wheel for source-ddc: filename=source_ddc-0.0.2.dev1-cp36-none-any.whl size=13685 sha256=00065f630718f8e541e3f73f0609156c382efa03518f34bf7d01cc1aa79eb658\n",
" Stored in directory: /tmp/pip-ephem-wheel-cache-ofoeiwb9/wheels/32/5f/1f/8908db9ccb4922cccb776a2fdab33c8a5c1e933f75d7c7991e\n",
"Successfully built source-ddc\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "oqMtIwgkVzwV",
"colab_type": "text"
},
"source": [
"## The Setting\n",
"\n",
"A person decides every year whether to retire or to continue working after considering the state of the environment defined by his age, health, the pension he would receive, his medical expenses and a disutility from beginning to work after having retired.\n",
"\n",
"The utility of retirement is given by: $ U_{retire} = pension_t - medical_expenses_t $\n",
"\n",
"The wage is given by: $ wage_t = \\exp(\\theta_0 + \\theta_1*health_t + \\theta_2*( \\frac{age_t}{1 + age_t})) $\n",
"\n",
"And his utility from working is given by: $ U_{work} = wage_t - medical_expenses_t + \\theta_3 is_retired_{t-1} $\n",
"\n",
"$\\theta_3$ is a disutility, so we expect it to be negative.\n",
"\n",
"Assume a discount factor of $0.9$\n",
"\n",
"First, let's define the imports:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "pADClAWIVt1Q",
"colab_type": "code",
"colab": {}
},
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from source_ddc.probability_tools import StateManager, random_ccp\n",
"from source_ddc.simulation_tools import simulate\n",
"from source_ddc.algorithms import NPL, NFXP, CCP\n",
"from source_ddc.fixed_point import get_discounted_value, phi_map, lambda_map"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "WcbCHSZhWYdj",
"colab_type": "text"
},
"source": [
"## Defining the state-space:\n",
"\n",
"The state-space is discrete, defined in the following way:\n",
"\n",
"- age: 30 categories\n",
"- health: 2 categories (bad, good)\n",
"- is retired: 2 categories\n",
"- pension: 3 categories (low, medium, high, in equal monetary intervals)\n",
"- medical_exp: 3 categories (low, medium, high, in equal monetary intervals)"
]
},
{
"cell_type": "code",
"metadata": {
"id": "KvoRY91dV6hs",
"colab_type": "code",
"colab": {}
},
"source": [
"# Number of states\n",
"\n",
"n_age_states = 30\n",
"n_health_states = 2\n",
"n_retirement_states = 2\n",
"n_pension_states = 3\n",
"n_medical_exp_states = 3"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "TBkk8zY_Weal",
"colab_type": "text"
},
"source": [
"## The discount factor\n",
"\n",
"Assume a discount factor of $0.9$"
]
},
{
"cell_type": "code",
"metadata": {
"id": "E347MykbWWlo",
"colab_type": "code",
"colab": {}
},
"source": [
"discount_factor = 0.9"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "IpbEcpkWWiM_",
"colab_type": "text"
},
"source": [
"## The Utility Function\n",
"\n",
"Source DDC takes a function that accepts three arguments: the paramaters (list of float), a numpy array of possible choices and a numpy array of states.\n",
"\n",
"The function should return a numpy array with a first dimension size equal to the number of actions, representing the utility of each of the actions respectively."
]
},
{
"cell_type": "code",
"metadata": {
"id": "Jr6nP9JIWgJL",
"colab_type": "code",
"colab": {}
},
"source": [
"# Utility function\n",
"\n",
"def utility_fn(theta, choices, states):\n",
" [age, health, retirement, pension, medical_exp] = states\n",
" \n",
" u_retires = pension - medical_exp\n",
" \n",
" wage = np.exp(theta[0] + theta[1]*health + theta[2]*(age/(1 + age)))\n",
" \n",
" u_works = wage - medical_exp\n",
" \n",
" u_works = np.where(retirement == 1, u_works + theta[3], u_works)\n",
" \n",
" u = np.array([u_retires, u_works])\n",
" \n",
" return u"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "7xx88TM0WlVF",
"colab_type": "text"
},
"source": [
"## The True Parameters\n",
"\n",
"Let's assume that the true parameters are the following (in order of $\\theta_0$ to $\\theta_3$):"
]
},
{
"cell_type": "code",
"metadata": {
"id": "LQrxqBLrWkkD",
"colab_type": "code",
"colab": {}
},
"source": [
"theta = [1, 0.8, -1, -2]"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "UXqEQg_RWrnh",
"colab_type": "text"
},
"source": [
"## The Transition Matrices\n",
"\n",
"Assume that they are defined as below. Note that if an agent will transition to retirement state with probability 1 if he chooses the retirement option. Conversely, if he chooses to work, he transitions to non-retirement with probability 1. All matrices are row-normalized to 1, since they represent transition probabilities. We need to define transition matrices for each choice: retirement or work."
]
},
{
"cell_type": "code",
"metadata": {
"id": "pc7xaTkfWo_I",
"colab_type": "code",
"colab": {}
},
"source": [
"# Transition functions retired:\n",
"health_transition_retired = np.array(\n",
" [\n",
" [0.2, 0.8], \n",
" [0.1, 0.9]\n",
" ]\n",
")\n",
"\n",
"is_retired_transition_retired = np.array(\n",
" [\n",
" [0.0, 1.0],\n",
" [0.0, 1.0]\n",
" ]\n",
")\n",
"\n",
"pension_transition_retired = np.eye(3)\n",
"\n",
"medical_exp_transition_retired = np.array(\n",
" [\n",
" [0.7, 0.2, 0.1],\n",
" [0.1, 0.7, 0.2],\n",
" [0.2, 0.1, 0.7]\n",
" ]\n",
")\n",
"\n",
"age_transition_retired = np.diag(np.ones(shape=(29)), k=1)\n",
"age_transition_retired[-1, -1] = 1\n",
"\n",
"\n",
"# Transition matrices working:\n",
"health_transition_working = np.array(\n",
" [\n",
" [0.4, 0.6], \n",
" [0.1, 0.9]\n",
" ]\n",
")\n",
"\n",
"is_retired_transition_working = np.array(\n",
" [\n",
" [1.0, 0.0],\n",
" [0.9, 0.1]\n",
" ]\n",
")\n",
"\n",
"pension_transition_working = np.array(\n",
" [\n",
" [0.7, 0.2, 0.1],\n",
" [0.1, 0.7, 0.2],\n",
" [0.1, 0.2, 0.7]\n",
" ]\n",
")\n",
"\n",
"medical_exp_transition_working = np.array(\n",
" [\n",
" [0.5, 0.4, 0.1],\n",
" [0.4, 0.5, 0.1],\n",
" [0.1, 0.5, 0.4]\n",
" ]\n",
")\n",
"\n",
"age_transition_working = np.diag(np.ones(shape=(29)), k=1)\n",
"age_transition_working[-1, -1] = 1"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "sDdefAS6WyAG",
"colab_type": "text"
},
"source": [
"## The State Manager\n",
"\n",
"`StateManager` is a class that helps with handling the different states. In the inside, it assigns a state ID to each possible combination of states. You can define it in the following way:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "oH9DsetiWrD9",
"colab_type": "code",
"colab": {}
},
"source": [
"state_manager = StateManager(\n",
" age=n_age_states,\n",
" health=n_health_states,\n",
" retirement=n_retirement_states,\n",
" pension=n_pension_states,\n",
" medical_exp=n_medical_exp_states,\n",
")"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "CSgaQ_OcW10Z",
"colab_type": "text"
},
"source": [
"You can also use StateManager's static function `merge_matrices` to apply the Kronecker product to all the transition matrices defined before, to turn them into a single transition matrix to use from now.\n",
"\n",
"Create one such combined matrix for the *retirement* and *working* choices, and stack them into a single transition matrix:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "ge7TKrr1WzqR",
"colab_type": "code",
"colab": {}
},
"source": [
"transition_matrix_working = StateManager.merge_matrices(age_transition_working, health_transition_working, is_retired_transition_working, pension_transition_working, medical_exp_transition_working)\n",
"transition_matrix_retired = StateManager.merge_matrices(age_transition_retired, health_transition_retired, is_retired_transition_retired, pension_transition_retired, medical_exp_transition_retired)\n",
"transition_matrix = np.array([transition_matrix_retired, transition_matrix_working])"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "ZwFnC24PW_hE",
"colab_type": "text"
},
"source": [
"## Simulating data\n",
"\n",
"Use the `simulate` function to simulate a panel dataset of *20* time periods and *200* agents. This function returns the dataset and the set of conditional choice probabilities:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "BgjaArKnW6kK",
"colab_type": "code",
"colab": {}
},
"source": [
"df, ccp = simulate(20, 200, 2, state_manager, theta, utility_fn, discount_factor, transition_matrix)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "sRSM627CW-Ps",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 297
},
"outputId": "d6776f0a-6db2-4afb-9fcb-847c23d534cf"
},
"source": [
"df.describe()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>agent_id</th>\n",
" <th>t</th>\n",
" <th>state</th>\n",
" <th>action</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>count</th>\n",
" <td>4000.000000</td>\n",
" <td>4000.000000</td>\n",
" <td>4000.000000</td>\n",
" <td>4000.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>99.500000</td>\n",
" <td>9.500000</td>\n",
" <td>808.020000</td>\n",
" <td>0.754500</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>57.741523</td>\n",
" <td>5.767002</td>\n",
" <td>279.044473</td>\n",
" <td>0.430437</td>\n",
" </tr>\n",
" <tr>\n",
" <th>min</th>\n",
" <td>0.000000</td>\n",
" <td>0.000000</td>\n",
" <td>9.000000</td>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25%</th>\n",
" <td>49.750000</td>\n",
" <td>4.750000</td>\n",
" <td>610.000000</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50%</th>\n",
" <td>99.500000</td>\n",
" <td>9.500000</td>\n",
" <td>890.000000</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>75%</th>\n",
" <td>149.250000</td>\n",
" <td>14.250000</td>\n",
" <td>1065.000000</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>max</th>\n",
" <td>199.000000</td>\n",
" <td>19.000000</td>\n",
" <td>1079.000000</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" agent_id t state action\n",
"count 4000.000000 4000.000000 4000.000000 4000.000000\n",
"mean 99.500000 9.500000 808.020000 0.754500\n",
"std 57.741523 5.767002 279.044473 0.430437\n",
"min 0.000000 0.000000 9.000000 0.000000\n",
"25% 49.750000 4.750000 610.000000 1.000000\n",
"50% 99.500000 9.500000 890.000000 1.000000\n",
"75% 149.250000 14.250000 1065.000000 1.000000\n",
"max 199.000000 19.000000 1079.000000 1.000000"
]
},
"metadata": {
"tags": []
},
"execution_count": 11
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "el3H9IJCXJUi",
"colab_type": "text"
},
"source": [
"You can recover the original state variables from their IDs:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "o_TsYtdTXGiz",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 419
},
"outputId": "c8981da7-161d-478c-9a0e-e064f5aaebc7"
},
"source": [
"pd.DataFrame(state_manager.state_variables_for(df['state'].values), columns=state_manager.state_names)"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>age</th>\n",
" <th>health</th>\n",
" <th>retirement</th>\n",
" <th>pension</th>\n",
" <th>medical_exp</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>7</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>8</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>9</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>10</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>2</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>11</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3995</th>\n",
" <td>24</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>2</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3996</th>\n",
" <td>25</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3997</th>\n",
" <td>26</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3998</th>\n",
" <td>27</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3999</th>\n",
" <td>28</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>4000 rows × 5 columns</p>\n",
"</div>"
],
"text/plain": [
" age health retirement pension medical_exp\n",
"0 7 0 0 0 0\n",
"1 8 0 0 1 2\n",
"2 9 0 0 1 1\n",
"3 10 1 0 2 2\n",
"4 11 1 0 2 1\n",
"... ... ... ... ... ...\n",
"3995 24 1 0 2 2\n",
"3996 25 1 0 1 2\n",
"3997 26 1 0 1 2\n",
"3998 27 1 0 1 1\n",
"3999 28 1 0 1 0\n",
"\n",
"[4000 rows x 5 columns]"
]
},
"metadata": {
"tags": []
},
"execution_count": 12
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Z5YKwstiXOsh",
"colab_type": "text"
},
"source": [
"## Define the algorithm\n",
"\n",
"Create an instance of the algorithm by passing all the required arguments. You can choose between NFXP, CCP and NPL, all of them included in the `algorithms` module. Choose the best one for the task according to their properties or try different ones to verify the stability of the results. Some may be more stable than others depending on the situation.\n",
"\n",
"Here we will assume that you have estimated consistently the transition probabilities and the choice probabilities and use Hotz & Miller's CCP algorithm:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "iE_oT8bOXIZ4",
"colab_type": "code",
"colab": {}
},
"source": [
"parameter_names=['const', 'health', 'age', 'work_disutility']\n",
"\n",
"ccp_algorithm = CCP(df['action'],\n",
" df['state'],\n",
" transition_matrix,\n",
" utility_fn,\n",
" discount_factor,\n",
" state_manager=state_manager,\n",
" parameter_names=parameter_names,\n",
" initial_p=ccp\n",
" )"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "wdIe2XlCXTPQ",
"colab_type": "text"
},
"source": [
"## Estimating\n",
"\n",
"Each algorithm implements the `estimate` method, to which you pass a set of initial parameter values and the optimization method. Here we'll pass a randomly generated set of initial values and use the BFGS method. The optimization logs can be a bit annoying, you can silence them by passing the `disp=0` option:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "0IB15_85XRA_",
"colab_type": "code",
"colab": {}
},
"source": [
"ccp_results = ccp_algorithm.estimate(start_params=np.random.uniform(size=len(parameter_names)), method='bfgs', disp=0)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "irmjXkziXczQ",
"colab_type": "text"
},
"source": [
"## The Results\n",
"\n",
"Source DDC models that are estimated by Maximum Likelihood extend Statsmodels' classes under the hood, and return a `GenericLikelihoodModelResults` object, where you'll find all the necessary information.\n",
"\n",
"You can verify that indeed we have obtained parameters that are very close to the true ones we defined at the beginning of this notebook."
]
},
{
"cell_type": "code",
"metadata": {
"id": "9bS9ZctXXZ80",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 311
},
"outputId": "5ba00879-c6ea-48c1-cf15-9dbc39e1aaa0"
},
"source": [
"ccp_results.summary()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>CCP Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>action</td> <th> Log-Likelihood: </th> <td> -1421.8</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>CCP</td> <th> AIC: </th> <td> 2846.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Maximum Likelihood</td> <th> BIC: </th> <td> 2852.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 02 Sep 2020</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>05:56:59</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 4000</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 3999</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 0</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>const</th> <td> 1.2097</td> <td> 0.285</td> <td> 4.242</td> <td> 0.000</td> <td> 0.651</td> <td> 1.769</td>\n",
"</tr>\n",
"<tr>\n",
" <th>health</th> <td> 0.7098</td> <td> 0.083</td> <td> 8.602</td> <td> 0.000</td> <td> 0.548</td> <td> 0.871</td>\n",
"</tr>\n",
"<tr>\n",
" <th>age</th> <td> -1.1443</td> <td> 0.303</td> <td> -3.776</td> <td> 0.000</td> <td> -1.738</td> <td> -0.550</td>\n",
"</tr>\n",
"<tr>\n",
" <th>work_disutility</th> <td> -1.9708</td> <td> 0.086</td> <td> -22.816</td> <td> 0.000</td> <td> -2.140</td> <td> -1.801</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" CCP Results \n",
"==============================================================================\n",
"Dep. Variable: action Log-Likelihood: -1421.8\n",
"Model: CCP AIC: 2846.\n",
"Method: Maximum Likelihood BIC: 2852.\n",
"Date: Wed, 02 Sep 2020 \n",
"Time: 05:56:59 \n",
"No. Observations: 4000 \n",
"Df Residuals: 3999 \n",
"Df Model: 0 \n",
"===================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"const 1.2097 0.285 4.242 0.000 0.651 1.769\n",
"health 0.7098 0.083 8.602 0.000 0.548 0.871\n",
"age -1.1443 0.303 -3.776 0.000 -1.738 -0.550\n",
"work_disutility -1.9708 0.086 -22.816 0.000 -2.140 -1.801\n",
"===================================================================================\n",
"\"\"\""
]
},
"metadata": {
"tags": []
},
"execution_count": 15
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "RGysZKdaXlUl",
"colab_type": "text"
},
"source": [
"The CCP method is quite fast and stable as long as you have good estimates of the conditional choice probabilities. If you don't have them, you can employ the NFXP algorithm to simultaneously estimate the utility function parameters and the conditional choice probabilities. It is a bit slower though."
]
},
{
"cell_type": "code",
"metadata": {
"id": "3IbQTsdLXk2G",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 396
},
"outputId": "3fe958d6-c017-4136-8ce5-8390b8f5fc47"
},
"source": [
"nfxp_algorithm = NFXP(df['action'],\n",
" df['state'],\n",
" transition_matrix,\n",
" utility_fn,\n",
" discount_factor,\n",
" state_manager=state_manager,\n",
" parameter_names=parameter_names\n",
" )\n",
"\n",
"nfxp_results = nfxp_algorithm.estimate(start_params=np.random.uniform(size=len(parameter_names)), method='bfgs')\n",
"\n",
"nfxp_results.summary()"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.355438\n",
" Iterations: 28\n",
" Function evaluations: 30\n",
" Gradient evaluations: 30\n"
],
"name": "stdout"
},
{
"output_type": "execute_result",
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>NFXP Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>action</td> <th> Log-Likelihood: </th> <td> -1421.8</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>NFXP</td> <th> AIC: </th> <td> 2846.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Maximum Likelihood</td> <th> BIC: </th> <td> 2852.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 02 Sep 2020</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>05:58:59</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 4000</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 3999</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 0</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>const</th> <td> 1.2113</td> <td> 0.288</td> <td> 4.208</td> <td> 0.000</td> <td> 0.647</td> <td> 1.775</td>\n",
"</tr>\n",
"<tr>\n",
" <th>health</th> <td> 0.7098</td> <td> 0.082</td> <td> 8.609</td> <td> 0.000</td> <td> 0.548</td> <td> 0.871</td>\n",
"</tr>\n",
"<tr>\n",
" <th>age</th> <td> -1.1460</td> <td> 0.306</td> <td> -3.746</td> <td> 0.000</td> <td> -1.746</td> <td> -0.546</td>\n",
"</tr>\n",
"<tr>\n",
" <th>work_disutility</th> <td> -1.9707</td> <td> 0.086</td> <td> -22.864</td> <td> 0.000</td> <td> -2.140</td> <td> -1.802</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" NFXP Results \n",
"==============================================================================\n",
"Dep. Variable: action Log-Likelihood: -1421.8\n",
"Model: NFXP AIC: 2846.\n",
"Method: Maximum Likelihood BIC: 2852.\n",
"Date: Wed, 02 Sep 2020 \n",
"Time: 05:58:59 \n",
"No. Observations: 4000 \n",
"Df Residuals: 3999 \n",
"Df Model: 0 \n",
"===================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"const 1.2113 0.288 4.208 0.000 0.647 1.775\n",
"health 0.7098 0.082 8.609 0.000 0.548 0.871\n",
"age -1.1460 0.306 -3.746 0.000 -1.746 -0.546\n",
"work_disutility -1.9707 0.086 -22.864 0.000 -2.140 -1.802\n",
"===================================================================================\n",
"\"\"\""
]
},
"metadata": {
"tags": []
},
"execution_count": 16
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "IZP2djtYXqwO",
"colab_type": "text"
},
"source": [
"The results should have converged, to something close to the true parameters. the probabilities should also be very close to the true ones given by the `ccp` object:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "cn-pPogkXocR",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "1c7d64dc-a5cc-463d-ab72-6d547a028c47"
},
"source": [
"np.corrcoef(nfxp_algorithm.p[0].ravel(), ccp[0].ravel())"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1. , 0.99803546],\n",
" [0.99803546, 1. ]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 17
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "ZH8ZUyh_Xumt",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "51289d97-626d-42e5-e1a6-ec83c06c10b7"
},
"source": [
"np.corrcoef(nfxp_algorithm.p[1].ravel(), ccp[1].ravel())"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1. , 0.99803546],\n",
" [0.99803546, 1. ]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 18
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "JY8zJegiXwAy",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"outputId": "13ec682b-1dc5-4c06-c44f-66bd2e56e2d8"
},
"source": [
"np.abs(nfxp_algorithm.p - ccp).mean()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.01243738407070375"
]
},
"metadata": {
"tags": []
},
"execution_count": 19
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "TRG981J5X0Ro",
"colab_type": "text"
},
"source": [
"The correlation between the probabilities between the estimates and the true ones at each choice is close enough to 1. The mean absolute error is also decently low, so we could say that the algorithm has converged.\n",
"\n",
"Finally, let's use the NPL algorithm, which also estimates the probabilities as well as the parameters and can be faster than the NFXP."
]
},
{
"cell_type": "code",
"metadata": {
"id": "GGAN45nhXxbM",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 396
},
"outputId": "f5d9ad11-55ab-4a1c-b5f9-7a8d9314d092"
},
"source": [
"npl_algorithm = NFXP(df['action'],\n",
" df['state'],\n",
" transition_matrix,\n",
" utility_fn,\n",
" discount_factor,\n",
" state_manager=state_manager,\n",
" parameter_names=parameter_names,\n",
" initial_p = random_ccp(state_manager.total_states, 2)\n",
" )\n",
"\n",
"npl_results = npl_algorithm.estimate(start_params=np.random.uniform(size=len(parameter_names)), method='bfgs')\n",
"\n",
"npl_results.summary()"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.355438\n",
" Iterations: 27\n",
" Function evaluations: 29\n",
" Gradient evaluations: 29\n"
],
"name": "stdout"
},
{
"output_type": "execute_result",
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>NFXP Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>action</td> <th> Log-Likelihood: </th> <td> -1421.8</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>NFXP</td> <th> AIC: </th> <td> 2846.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Maximum Likelihood</td> <th> BIC: </th> <td> 2852.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 02 Sep 2020</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>06:00:56</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 4000</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 3999</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 0</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>const</th> <td> 1.2118</td> <td> 0.288</td> <td> 4.210</td> <td> 0.000</td> <td> 0.648</td> <td> 1.776</td>\n",
"</tr>\n",
"<tr>\n",
" <th>health</th> <td> 0.7100</td> <td> 0.082</td> <td> 8.609</td> <td> 0.000</td> <td> 0.548</td> <td> 0.872</td>\n",
"</tr>\n",
"<tr>\n",
" <th>age</th> <td> -1.1467</td> <td> 0.306</td> <td> -3.748</td> <td> 0.000</td> <td> -1.746</td> <td> -0.547</td>\n",
"</tr>\n",
"<tr>\n",
" <th>work_disutility</th> <td> -1.9707</td> <td> 0.086</td> <td> -22.864</td> <td> 0.000</td> <td> -2.140</td> <td> -1.802</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" NFXP Results \n",
"==============================================================================\n",
"Dep. Variable: action Log-Likelihood: -1421.8\n",
"Model: NFXP AIC: 2846.\n",
"Method: Maximum Likelihood BIC: 2852.\n",
"Date: Wed, 02 Sep 2020 \n",
"Time: 06:00:56 \n",
"No. Observations: 4000 \n",
"Df Residuals: 3999 \n",
"Df Model: 0 \n",
"===================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"const 1.2118 0.288 4.210 0.000 0.648 1.776\n",
"health 0.7100 0.082 8.609 0.000 0.548 0.872\n",
"age -1.1467 0.306 -3.748 0.000 -1.746 -0.547\n",
"work_disutility -1.9707 0.086 -22.864 0.000 -2.140 -1.802\n",
"===================================================================================\n",
"\"\"\""
]
},
"metadata": {
"tags": []
},
"execution_count": 20
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "WzZZEXRoX7pc",
"colab_type": "text"
},
"source": [
"The NPL algorithm should also have also converged to the true parameter values or something close. The probabilities should have as well."
]
},
{
"cell_type": "code",
"metadata": {
"id": "zAKY90TJX5Ms",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "02167d75-2984-4160-8ad8-8710cbba3916"
},
"source": [
"np.corrcoef(npl_algorithm.p[0].ravel(), ccp[0].ravel())"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1. , 0.99803395],\n",
" [0.99803395, 1. ]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 21
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "V0IxMgjgX9go",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "18df60d1-591a-45d9-bbb6-882332255b11"
},
"source": [
"np.corrcoef(npl_algorithm.p[1].ravel(), ccp[1].ravel())"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1. , 0.99803395],\n",
" [0.99803395, 1. ]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 22
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "WL1jrH5BX-xu",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"outputId": "577bacbc-3e6a-432f-b085-88e2f5a27e3c"
},
"source": [
"np.abs(npl_algorithm.p - ccp).mean()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.012425489193289835"
]
},
"metadata": {
"tags": []
},
"execution_count": 23
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "HvaYX4W7YD3Q",
"colab_type": "text"
},
"source": [
"We have shown how to use the three algorithms and that all of them return good approximations of the true utility function parameters in this non-trivial toy model. Which one to use should be decided on the base of each algorithm's asymptotic properties, speed and the data available.\n",
"\n",
"We'll appreciate your bug reports and feature requests at https://github.com/sansan-inc/econ-source"
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment