Problem Set 1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preliminary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "MISSING_INT = -99\n",
    "INVALID_FLOAT = -99.00\n",
    "HOURS = [8, 12, 16]\n",
    "T = 24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_LABLES_SIM = [\n",
    "    \"Period\", \n",
    "    \"Experience\", \n",
    "    \"Choice_Lagged\",\n",
    "    \"Choice\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_spec = dict()\n",
    "model_spec[\"num_periods\"] = 3\n",
    "model_spec[\"num_choices\"] = 3\n",
    "model_spec[\"experience_factor_0\"] = 0.25\n",
    "model_spec[\"experience_factor_1\"] = 0.15\n",
    "model_spec[\"delta\"] = 1\n",
    "\n",
    "model_spec = collections.namedtuple(\"model_specification\", model_spec.keys())(\n",
    "        **model_spec\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params = dict()\n",
    "model_params[\"alpha\"] = 0.5\n",
    "model_params[\"beta\"] = 0.5\n",
    "model_params = collections.namedtuple(\"model_specification\", model_params.keys())(\n",
    "        **model_params\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "model_specification(num_periods=3, num_choices=3, experience_factor_0=0.25, experience_factor_1=0.15, delta=1)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "model_specification(alpha=0.5, beta=0.5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create states and indexer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_state_space(model_spec):  \n",
    "    \"\"\"This function creates the state space matrix\n",
    "    and assignes an index value to every state\"\"\"\n",
    "    \n",
    "    data = []\n",
    "    \n",
    "    shape = (\n",
    "        model_spec.num_periods,\n",
    "        model_spec.num_periods*2,\n",
    "        model_spec.num_choices,\n",
    "    )\n",
    "    \n",
    "    indexer = np.full(shape, MISSING_INT)\n",
    "    \n",
    "    # Initialize index counter\n",
    "    i = 0\n",
    "    \n",
    "    # Loop over all periods\n",
    "    for period in range(model_spec.num_periods):\n",
    "        \n",
    "        # Define entry state in period zero:\n",
    "        if period == 0:\n",
    "\n",
    "            indexer[period, 0, 0] = i\n",
    "            row = [period, 0, 0]\n",
    "\n",
    "            # Update counter\n",
    "            i += 1\n",
    "\n",
    "            # Add to state space matrix\n",
    "            data.append(row)\n",
    "        \n",
    "        else:\n",
    "        \n",
    "            # Loop over all admissible values of experience\n",
    "            for exp in range((period)*2 + 1):\n",
    "\n",
    "                # Loop over all choices:\n",
    "                for choice_lagged in range(model_spec.num_choices):\n",
    "                    \n",
    "                    # Exclude impossible states:\n",
    "                    if (exp < choice_lagged):\n",
    "                        continue\n",
    "                        \n",
    "                    if (choice_lagged == 0) and (exp > 2*(period - 1)):\n",
    "                        continue\n",
    "                    \n",
    "                    if (choice_lagged == 1) and (exp > 2*(period - 1) + 1):\n",
    "                        continue\n",
    "                    \n",
    "                    if (choice_lagged != 0) and (exp == 0):\n",
    "                        continue\n",
    "\n",
    "                    # Check for duplicate states:\n",
    "                    if (\n",
    "                        indexer[period, exp, choice_lagged] != MISSING_INT\n",
    "                    ):\n",
    "                        continue\n",
    "\n",
    "                    # Record index of currently reached admissible state\n",
    "                    # space point\n",
    "                    indexer[\n",
    "                        period, exp, choice_lagged\n",
    "                    ] = i\n",
    "\n",
    "                    # Update count\n",
    "                    i += 1\n",
    "\n",
    "                    # Add to state space matrix\n",
    "                    row = [period, exp, choice_lagged]\n",
    "                    data.append(row)\n",
    "    \n",
    "        states = np.array(data)\n",
    "    \n",
    "    return states, indexer\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "states, indexer = create_state_space(model_spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0],\n",
       "       [1, 0, 0],\n",
       "       [1, 1, 1],\n",
       "       [1, 2, 2],\n",
       "       [2, 0, 0],\n",
       "       [2, 1, 0],\n",
       "       [2, 1, 1],\n",
       "       [2, 2, 0],\n",
       "       [2, 2, 1],\n",
       "       [2, 2, 2],\n",
       "       [2, 3, 1],\n",
       "       [2, 3, 2],\n",
       "       [2, 4, 2]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indexer[2,0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 0])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "states[4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### State space derived covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "produce = np.full(states.shape[0], MISSING_INT)\n",
    "produce = np.where(states[:, 1] == 0, 1, produce)\n",
    "produce = np.where((states[:, 1] != 0) & (states[:, 0] == 1), (1 + states[:, 1]*model_spec.experience_factor_0), produce)\n",
    "produce = np.where((states[:, 1] != 0) & (states[:, 0] == 2), (1 + (states[:, 1] - states[:, 2])*model_spec.experience_factor_0 + states[:, 2]*model_spec.experience_factor_1), produce)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backward induction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_continuation_values(states_subset, indexer, emaxs):\n",
    "    \"\"\"Obtain continuation values for each of the choices at each state\n",
    "    of the period currently reached by the parent loop.\"\"\"\n",
    "    \n",
    "    for i in range(states_subset.shape[0]):\n",
    "    \n",
    "        # Unpack parent state and get index\n",
    "        (\n",
    "            period,\n",
    "            exp,\n",
    "            choice_lagged,\n",
    "        ) = states_subset[i]\n",
    "\n",
    "        k_parent = indexer[period, exp, choice_lagged]\n",
    "\n",
    "        # Choice 8 hours\n",
    "        k_0 = indexer[period + 1, exp, 0]\n",
    "\n",
    "        # Choice 12 hours\n",
    "        k_1 = indexer[period + 1, exp + 1, 1]\n",
    "\n",
    "        # Choice 16 hours\n",
    "        k_2 = indexer[period + 1, exp + 2, 2]\n",
    "\n",
    "        # Get emax\n",
    "        emaxs[k_parent, 0] = emaxs[k_0, 3]\n",
    "        emaxs[k_parent, 1] = emaxs[k_1, 3]\n",
    "        emaxs[k_parent, 2] = emaxs[k_2, 3]\n",
    "    \n",
    "    return emaxs\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_emax(model_spec, model_params, produce_period, flow_utilities_period, value_functions_period, emaxs):\n",
    "    \"\"\"Calculates the maximum value function as the flow ulity plus continuatin value\n",
    "    of the most profitable choice\"\"\"\n",
    "    emax_period = np.full(emaxs.shape[0], INVALID_FLOAT) \n",
    "    \n",
    "    for i in range(emaxs.shape[0]):\n",
    "        current_max_value_function = INVALID_FLOAT\n",
    "\n",
    "        for j in range(model_spec.num_choices):\n",
    "\n",
    "            #flow_utility = ((produce[i]*HOURS[j])**model_params.alpha)*((T - HOURS[j])**model_params.beta)\n",
    "            flow_utility = ((produce[i])**model_params.alpha)*((T - HOURS[j])**model_params.beta)\n",
    "            flow_utilities_period[i, j] = flow_utility\n",
    "\n",
    "            value_function_choice = flow_utility + model_spec.delta * emaxs[i, j]\n",
    "            value_functions_period[i,j] = value_function_choice\n",
    "\n",
    "            if value_function_choice > current_max_value_function:\n",
    "                current_max_value_function = value_function_choice\n",
    "            \n",
    "            emax_period[i] = current_max_value_function\n",
    "\n",
    "    return flow_utilities_period, value_functions_period, emax_period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_induction(model_spec, model_params, states, indexer, produce):\n",
    "    \"\"\"This function performs the backward induction and returns\n",
    "    the model solution\"\"\"\n",
    "    \n",
    "    emaxs = np.zeros((states.shape[0], model_spec.num_choices + 1))\n",
    "    flow_utilities = np.zeros((states.shape[0], model_spec.num_choices))\n",
    "    value_functions = np.zeros((states.shape[0], model_spec.num_choices))\n",
    "    \n",
    "    \n",
    "    # Loop backwards over all periods:\n",
    "    for period in reversed(range(model_spec.num_periods)):\n",
    "        \n",
    "        # Extract period information\n",
    "        # States\n",
    "        states_period = states[np.where(states[:, 0] == period)]\n",
    "        # Produce\n",
    "        produce_period = produce[states[:, 0] == period]\n",
    "        \n",
    "        # Continuation value calculation not performed for last period\n",
    "        # since continuation values are known to be zero\n",
    "        if period == model_spec.num_periods - 1:\n",
    "            pass\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            # Fill first block of elements in emaxs for the current period\n",
    "            # corresponding to the continuation values\n",
    "            emaxs = get_continuation_values(\n",
    "                states_period,\n",
    "                indexer,\n",
    "                emaxs,\n",
    "            )\n",
    "            \n",
    "        # Extract current period information for current loop calculation\n",
    "        emaxs_period = emaxs[np.where(states[:, 0] == period)]\n",
    "        flow_utilities_period = flow_utilities[np.where(states[:, 0] == period)]\n",
    "        value_functions_period = flow_utilities[np.where(states[:, 0] == period)]\n",
    "            \n",
    "        # Calculate emax for current period reached by the loop\n",
    "        flow_utilities_period, value_functions_period, emax_period = construct_emax(\n",
    "            model_spec,\n",
    "            model_params,\n",
    "            produce_period,\n",
    "            flow_utilities_period,\n",
    "            value_functions_period,\n",
    "            emaxs_period[:, :3],\n",
    "        )\n",
    "        \n",
    "        emaxs_period[:, 3] = emax_period\n",
    "        emaxs[np.where(states[:, 0] == period)] = emaxs_period\n",
    "        flow_utilities[np.where(states[:, 0] == period)] = flow_utilities_period\n",
    "        value_functions[np.where(states[:, 0] == period)] = value_functions_period\n",
    "        \n",
    "    return flow_utilities, value_functions, emaxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow_utilities, value_functions, emaxs = backward_induction(model_spec, model_params, states, indexer, produce)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4.        , 3.46410162, 2.82842712],\n",
       "       [4.        , 3.46410162, 2.82842712],\n",
       "       [4.        , 3.46410162, 2.82842712],\n",
       "       [4.47213595, 3.87298335, 3.16227766],\n",
       "       [4.        , 3.46410162, 2.82842712],\n",
       "       [4.        , 3.46410162, 2.82842712],\n",
       "       [4.47213595, 3.87298335, 3.16227766],\n",
       "       [4.89897949, 4.24264069, 3.46410162],\n",
       "       [4.        , 3.46410162, 2.82842712],\n",
       "       [4.47213595, 3.87298335, 3.16227766],\n",
       "       [4.28952212, 3.71483512, 3.03315018],\n",
       "       [4.89897949, 4.24264069, 3.46410162],\n",
       "       [4.73286383, 4.09878031, 3.34664011]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flow_utilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[12.        , 11.46410162, 12.19954257],\n",
       "       [ 8.        ,  7.93623757,  7.30056308],\n",
       "       [ 8.        ,  7.46410162,  7.72740661],\n",
       "       [ 9.37111544,  8.16250546,  7.89514149],\n",
       "       [ 4.        ,  3.46410162,  2.82842712],\n",
       "       [ 4.        ,  3.46410162,  2.82842712],\n",
       "       [ 4.47213595,  3.87298335,  3.16227766],\n",
       "       [ 4.89897949,  4.24264069,  3.46410162],\n",
       "       [ 4.        ,  3.46410162,  2.82842712],\n",
       "       [ 4.47213595,  3.87298335,  3.16227766],\n",
       "       [ 4.28952212,  3.71483512,  3.03315018],\n",
       "       [ 4.89897949,  4.24264069,  3.46410162],\n",
       "       [ 4.73286383,  4.09878031,  3.34664011]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value_functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.        ,  8.        ,  9.37111544, 12.19954257],\n",
       "       [ 4.        ,  4.47213595,  4.47213595,  8.        ],\n",
       "       [ 4.        ,  4.        ,  4.89897949,  8.        ],\n",
       "       [ 4.89897949,  4.28952212,  4.73286383,  9.37111544],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.47213595],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.89897949],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.47213595],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.28952212],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.89897949],\n",
       "       [ 0.        ,  0.        ,  0.        ,  4.73286383]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emaxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0],\n",
       "       [1, 0, 0],\n",
       "       [1, 1, 1],\n",
       "       [1, 2, 2],\n",
       "       [2, 0, 0],\n",
       "       [2, 1, 0],\n",
       "       [2, 1, 1],\n",
       "       [2, 2, 0],\n",
       "       [2, 2, 1],\n",
       "       [2, 2, 2],\n",
       "       [2, 3, 1],\n",
       "       [2, 3, 2],\n",
       "       [2, 4, 2]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate(model_spec, model_params, produce, emaxs):\n",
    "    \n",
    "    data = []\n",
    "    hours = np.asarray(HOURS)\n",
    "    \n",
    "    for period in range(model_spec.num_periods):\n",
    "        \n",
    "        if period == 0:    \n",
    "            current_states = states[states[:, 0] == period]\n",
    "        \n",
    "        else:\n",
    "            current_states = current_states\n",
    "\n",
    "        idx = indexer[\n",
    "            current_states[:, 0],\n",
    "            current_states[:, 1],\n",
    "            current_states[:, 2]\n",
    "        ]\n",
    "        \n",
    "        #flow_utilities = ((produce[idx]*hours)**model_params.alpha)*((T - hours)**model_params.beta)\n",
    "        flow_utilities = ((produce[idx])**model_params.alpha)*((T - hours)**model_params.beta)\n",
    "\n",
    "        value_functions = flow_utilities + model_spec.delta * emaxs[idx, :3]\n",
    "        \n",
    "        # Determine choice as option with highest choice specific value function\n",
    "        choice = np.argmax(value_functions, axis=1)\n",
    "        \n",
    "        # Record period experiences\n",
    "        rows = np.column_stack(\n",
    "            (\n",
    "                current_states.copy(),\n",
    "                choice,\n",
    "                #flow_utilities,\n",
    "                #continuation_values,\n",
    "                #value_functions,\n",
    "            )\n",
    "        )\n",
    "\n",
    "        data.append(rows)\n",
    "        \n",
    "        # Update current states according to choice\n",
    "        current_states[:, 0] += 1\n",
    "        current_states[:, 2] = choice\n",
    "        current_states[:, 1] = np.where(\n",
    "            choice == 1, current_states[:, 1] + 1, current_states[:, 1]\n",
    "        )\n",
    "        current_states[:, 1] = np.where(\n",
    "            choice == 2, current_states[:, 1] + 2, current_states[:, 1]\n",
    "        )\n",
    "\n",
    "        dataset = pd.DataFrame(np.vstack(data), columns=DATA_LABLES_SIM)\n",
    "    \n",
    "    return dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = simulate(model_spec, model_params, produce, emaxs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Period</th>\n",
       "      <th>Experience</th>\n",
       "      <th>Choice_Lagged</th>\n",
       "      <th>Choice</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Period  Experience  Choice_Lagged  Choice\n",
       "0       0           0              0       2\n",
       "1       1           2              2       0\n",
       "2       2           2              0       0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}