{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "96656f8d",
   "metadata": {},
   "source": [
    "\n",
    "File: 04-wooldridge_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Sun Aug 10 2025, 17:31:54\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84553ae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import biogeme.biogeme_logging as blog\n",
    "import pandas as pd\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.database import Database\n",
    "from biogeme.expressions import (\n",
    "    Beta,\n",
    "    Draws,\n",
    "    Expression,\n",
    "    MonteCarlo,\n",
    "    PanelLikelihoodTrajectory,\n",
    "    Variable,\n",
    "    log,\n",
    ")\n",
    "from biogeme.models import logit\n",
    "from biogeme.results_processing import get_pandas_estimated_parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73bf54cd",
   "metadata": {},
   "source": [
    "As the estimation time may be long, we ask Biogeme to report the details of the iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2aa2eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4723e8b9",
   "metadata": {},
   "source": [
    "**Tip:**<div class=\"alert alert-block alert-info\">It is advised to start working with a low number of draws, until\n",
    "the script is working well. Then, increase the number of draws to 10000, say. Then, execute the script overnight.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25c128d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_draws = 10\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a7a0c31",
   "metadata": {},
   "source": [
    "# Dynamic Choice Models with Panel Effects"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a02793",
   "metadata": {},
   "source": [
    "We analyze again the smoking behavior of individuals, as a function of their age and the price of tobacco using\n",
    "synthetic data. We develop a model that predicts, for every year, the probability to smoke or not."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb88fc86",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9398707",
   "metadata": {},
   "source": [
    "We postulate a true model for the data generation process. It is a mixture of logit models where the utility\n",
    "associated with \"not smoking\" is\n",
    "\\begin{equation}\n",
    "U_{0nt}= \\varepsilon_{0nt}\n",
    "\\end{equation}\n",
    "and the utility associated with \"smoking\" is\n",
    "\\begin{equation}\n",
    "U_{1nt}= \\beta_{nt} y_{n,t-1} + \\beta^p_{nt} P_{t} + c_n + \\varepsilon_{1nt},\n",
    "\\end{equation}\n",
    "where\n",
    "\n",
    "- $\\beta_{nt} = 10$,\n",
    "\n",
    "- $y_{n,t-1}=1$ if $n$ is smoking at time $t-1$, $0$ otherwise,\n",
    "\n",
    "- $\\beta^p_{nt} = -0.1$,\n",
    "\n",
    "- $P_t$ is the price of cigarets at time $t$,\n",
    "\n",
    "- $c_n$ is an individual specific constant that captures the a priori, intrinsic attraction of each individual\n",
    "towards smoking. It is assumed to be normally distributed in the population, with zero mean and standard deviation\n",
    "50: $N(0, 50^2)$,"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84796d61",
   "metadata": {},
   "source": [
    "## True value of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d6fc4c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "true_parameters = pd.DataFrame(\n",
    "    {'Value': [-0.1, 10, 0, 50]},\n",
    "    index=['coef_price', 'beta_last_year', 'cte_mean', 'cte_std'],\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a664a44c",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23be8cb4",
   "metadata": {},
   "source": [
    "We observe every individual only from the age of 45 and the age of 55."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be01b16d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_table('smoking55.dat', sep=',')\n",
    "display(df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a7d571",
   "metadata": {},
   "source": [
    "Different values for age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57340d0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df['Age'].unique())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "937d4861",
   "metadata": {},
   "source": [
    "The data contains the following columns:\n",
    "\n",
    "- the age of the individual,\n",
    "- the price of the cigarettes,\n",
    "- a variable that is 1 if the individual is smoking, 0 otherwise,\n",
    "- a variable that is 1 if the individual was smoking last year, 0 otherwise,\n",
    "- a unique id for each individual,\n",
    "- a variable that is 1 if the individual was smoking at the age of 45, in the beginning of the observation period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f8df571",
   "metadata": {},
   "outputs": [],
   "source": [
    "database = Database('smoking55', df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25d82734",
   "metadata": {},
   "source": [
    "Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1022c6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Price = Variable('Price')\n",
    "Smoking = Variable('Smoking')\n",
    "LastYear = Variable('LastYear')\n",
    "Smoking45 = Variable('Smoking45')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dafda56",
   "metadata": {},
   "source": [
    "We declare that the data set contains panel data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97406bd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "database.panel('Id')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1631eb77",
   "metadata": {},
   "source": [
    "## Estimation procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe03d53d",
   "metadata": {},
   "source": [
    "The following procedure estimates the choice model, and returns the estimated parameters in a Pandas format. If the\n",
    "model happens to have been already estimated, the estimation results are read from the pickle file and reported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "342a9c32",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate(\n",
    "    the_logprob: Expression, the_name: str, the_database: Database\n",
    ") -> pd.DataFrame:\n",
    "    \"\"\"Estimates the choice model (or read the estimation results from\n",
    "        file if recycle is True), and returns the estimated parameters\n",
    "        in a Pandas format.\n",
    "\n",
    "    :param the_logprob: formula for the log likelihood function\n",
    "    :param the_name: name of the model (important for the output files)\n",
    "    :param the_database: database\n",
    "    :return: estimated values of the parameters and statistics.\n",
    "    \"\"\"\n",
    "    biogeme = BIOGEME(\n",
    "        the_database,\n",
    "        the_logprob,\n",
    "        number_of_draws=number_of_draws,\n",
    "    )\n",
    "    biogeme.model_name = the_name\n",
    "    results = biogeme.estimate()\n",
    "    print(results.short_summary())\n",
    "    pandas_results = get_pandas_estimated_parameters(estimation_results=results)\n",
    "    return pandas_results.set_index('Name')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b56e06cb",
   "metadata": {},
   "source": [
    "## Dynamic model with serial correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f83b75fd",
   "metadata": {},
   "source": [
    "In the previous quiz, we have estimated a dynamic model with panel effects to account for serial correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02febaf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cte_mean = Beta('cte_mean', 0, None, None, 0)\n",
    "cte_std = Beta('cte_std', 1, None, None, 0)\n",
    "cte = cte_mean + cte_std * Draws('agent', 'NORMAL_ANTI')\n",
    "coef_price = Beta('coef_price', 0, None, None, 0)\n",
    "beta_last_year = Beta('beta_last_year', 0, None, None, 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3498bc0",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df6b8b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_smoking = beta_last_year * LastYear + coef_price * Price + cte\n",
    "v_non_smoking = 0\n",
    "v = {0: v_non_smoking, 1: v_smoking}\n",
    "conditional_observation_probability = logit(v, None, Smoking)\n",
    "conditional_individual_probability = PanelLikelihoodTrajectory(\n",
    "    conditional_observation_probability\n",
    ")\n",
    "log_probability = log(MonteCarlo(conditional_individual_probability))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10bf339b",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b11bf80",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_serial_dynamic = estimate(log_probability, 'dynamic_model_serial', database)\n",
    "display(r_serial_dynamic)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54ddc0c6",
   "metadata": {},
   "source": [
    "### Comparison of the estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "942ae1bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = pd.concat(\n",
    "    [\n",
    "        true_parameters['Value'],\n",
    "        r_serial_dynamic['Value'],\n",
    "    ],\n",
    "    axis=1,\n",
    ")\n",
    "summary.columns = ['True', 'Dynamic + serial']\n",
    "summary.fillna('')\n",
    "display(summary)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "636e3dfa",
   "metadata": {},
   "source": [
    "We observe here the issue of the \"initial condition problem\". Although the model specification is correct (it is the\n",
    "same model as the data generation process), the values of the parameters are not correctly recovered. It is because\n",
    "the first observed choice, that is, the fact that an individual is smoking at the age of 45, is strongly correlated\n",
    "with the agent effect. This creates endogeneity. One visible consequence is the positive price coefficient. We are\n",
    "now using Wooldridge method to address it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b9360e",
   "metadata": {},
   "source": [
    "## Dynamic model with serial correlation and Wooldridge term"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20313cfc",
   "metadata": {},
   "source": [
    "We introduce in the specification of the constant a term that captures the fact that somebody is smoking during the\n",
    "first observation period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae2eaec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "coef_first_year = Beta('coef_first_year', 0, None, None, 0)\n",
    "cte = cte_mean + coef_first_year * Smoking45 + cte_std * Draws('agent', 'NORMAL_ANTI')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26fc36c7",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4786d47",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_smoking = beta_last_year * LastYear + coef_price * Price + cte\n",
    "v_non_smoking = 0\n",
    "v = {0: v_non_smoking, 1: v_smoking}\n",
    "conditional_observation_probability = logit(v, None, Smoking)\n",
    "conditional_individual_probability = PanelLikelihoodTrajectory(\n",
    "    conditional_observation_probability\n",
    ")\n",
    "log_probability = log(MonteCarlo(conditional_individual_probability))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "763e755b",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b993d8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_wooldridge = estimate(\n",
    "    log_probability, 'dynamic_model_serial_wooldridge_truncated_t_55', database\n",
    ")\n",
    "display(r_wooldridge)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e933fe1",
   "metadata": {},
   "source": [
    "### Comparison of the estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d541a72f",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = pd.concat(\n",
    "    [true_parameters['Value'], r_serial_dynamic['Value'], r_wooldridge['Value']], axis=1\n",
    ")\n",
    "summary.columns = [\n",
    "    'True',\n",
    "    'Dynamic + serial',\n",
    "    'Wooldridge',\n",
    "]\n",
    "summary.fillna('')\n",
    "display(summary)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fefbad0",
   "metadata": {},
   "source": [
    "The estimates of the coefficients `coef_price` and `beta_last_year` are now closer to their true value. We perform\n",
    "a $t$-test analysis, to test the hypothesis that the value of the parameter is equal to its true value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edbee180",
   "metadata": {},
   "outputs": [],
   "source": [
    "def t_test(param: str, true_value: float) -> float:\n",
    "    \"\"\"Calculates the t-statistic to test the hypothesis that the true\n",
    "        value of the parameter is indeed the true one.\n",
    "\n",
    "    :param param: name of the parameter\n",
    "    :param true_value: true value of the parameter\n",
    "    :return: t-statistic\n",
    "    \"\"\"\n",
    "    return (r_wooldridge.loc[param, 'Value'] - true_value) / r_wooldridge.loc[\n",
    "        param, 'Robust std err.'\n",
    "    ]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61d2ab9b",
   "metadata": {},
   "source": [
    "`coef_price`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfd77b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "the_test = t_test('coef_price', -0.1)\n",
    "print(f'Test for coef_price: {the_test:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37ae1834",
   "metadata": {},
   "source": [
    "`beta_last_year`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c603c2fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "the_test = t_test('beta_last_year', 10)\n",
    "print(f'Test for beta_last_year: {the_test:.3g}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1466bb32",
   "metadata": {},
   "source": [
    "`cte_mean`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8b9cdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "the_test = t_test('cte_mean', 0)\n",
    "print(f'Test for cte_mean: {the_test:.3g}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "677b8a80",
   "metadata": {},
   "source": [
    "`cte_std`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac6eacb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "the_test = t_test('cte_std', 50)\n",
    "print(f'Test for cte_std: {the_test:.3g}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d363bb24",
   "metadata": {},
   "source": [
    "Except for `cte_std`, the $t$-test are, in absolute value, below 1.96. It means that we cannot reject the null\n",
    "hypothesis that the value of the parameter is equal to the true value, at the 95% level of confidence. It is\n",
    "important to realize from the relatively large value of the standard errors that the precision of the estimates is\n",
    "not high.   This is due to a lack of observations. Indeed, we observe each individual only during 10 years.\n",
    "But the Wooldridge correction has allowed to address the endogeneity issue."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}