{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2956736f",
   "metadata": {},
   "source": [
    "\n",
    "File: 02-dynamic_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Thu Aug 07 2025, 17:43:46\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97fd6b01",
   "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, loglogit\n",
    "from biogeme.results_processing import (\n",
    "    EstimationResults,\n",
    "    get_pandas_estimated_parameters,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83331f18",
   "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": "ba3a4248",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12cffa2e",
   "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": "99df3f84",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_draws = 10\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d9b5a38",
   "metadata": {},
   "source": [
    "# Dynamic Choice Models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5428970",
   "metadata": {},
   "source": [
    "We analyze the smoking behavior of individuals, as a function of their age and the price of tobacco using synthetic\n",
    "data. We develop a model that predicts, for every year, the probability to smoke or not."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "320376d9",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8de6f3dd",
   "metadata": {},
   "source": [
    "We use synthetic data that has been generated as follows. We postulate a true model for the data generation process.\n",
    "It is a mixture of logit models with\n",
    "two alternatives: ``smoking`` or ``not smoking``.\n",
    "The utility for individual $n$ associated with \"not smoking\" in year $t$ 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 cigarettes 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)$, and constant over $t$.\n",
    "\n",
    "We generate a sample of 1000 individuals, and we simulate their smoking behavior between the age of 16 until the age\n",
    "of 100.\n",
    "\n",
    "The date of birth of each individual is uniformly distributed between 2000 and 2020.\n",
    "The price of cigarettes in 2000 is supposed to be 10. The price of cigarettes in year $t$ is\n",
    "$$P_t = 10 \\cdot 1.02^{t-2000},$$\n",
    "which represents a price increase of 2% per year."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e670b26",
   "metadata": {},
   "source": [
    "## True value of the parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b919858d",
   "metadata": {},
   "source": [
    "We store the true value of the parameters for future comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76894987",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "aaba528d",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be94adfe",
   "metadata": {},
   "source": [
    "The observations are available in the following data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea3fdaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_table('smoking.dat', sep=',')\n",
    "display(df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b6a7b7",
   "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": "6297d8ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "database = Database('smoking', df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a27483e",
   "metadata": {},
   "source": [
    "Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90e8fef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "Price = Variable('Price')\n",
    "Smoking = Variable('Smoking')\n",
    "LastYear = Variable('LastYear')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58d60c4c",
   "metadata": {},
   "source": [
    "## Estimation procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b6ce11",
   "metadata": {},
   "source": [
    "The following procedure estimates the choice model (or read the estimation results from file), and returns the\n",
    "estimated parameters in a Pandas format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b67839b",
   "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",
    "    print(f'Sample size: {biogeme.sample_size}')\n",
    "    print(f'Number of observations: {biogeme.number_of_observations}')\n",
    "\n",
    "    results: EstimationResults = 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": "daae420e",
   "metadata": {},
   "source": [
    "## Static model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dee090fd",
   "metadata": {},
   "source": [
    "The static model considers the data as cross-sectional. No state dependence, and no serial correlation is captured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75ab690",
   "metadata": {},
   "outputs": [],
   "source": [
    "cte_mean = Beta('cte_mean', 0, None, None, 0)\n",
    "coef_price = Beta('coef_price', 0, None, None, 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c8dc5af",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff13d9b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_smoking = coef_price * Price + cte_mean\n",
    "v_non_smoking = 0\n",
    "v = {0: v_non_smoking, 1: v_smoking}\n",
    "log_probability = loglogit(v, None, Smoking)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a31b13",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c73cf8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_static = estimate(log_probability, 'static_model', database)\n",
    "display(r_static)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3b7cbaa",
   "metadata": {},
   "source": [
    "## Comparison of the estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "038c4456",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = pd.concat([true_parameters['Value'], r_static['Value']], axis='columns')\n",
    "summary.columns = ['True', 'Static']\n",
    "summary.fillna('')\n",
    "display(summary)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "443a711c",
   "metadata": {},
   "source": [
    "The estimated price coefficient is not significant. Indeed, price is the only variable that the model considers.\n",
    "Ignoring state dependence generates endogeneity. The model \"thinks\" that individuals are insensitive to price,\n",
    "as they choose an alternative that is expensive."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81368b19",
   "metadata": {},
   "source": [
    "## Dynamic model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a675207",
   "metadata": {},
   "source": [
    "The dynamic model adds the choice of last year as an explanatory variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16e53196",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_last_year = Beta('beta_last_year', 0, None, None, 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "582ff6a4",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa807899",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_smoking = beta_last_year * LastYear + coef_price * Price + cte_mean\n",
    "v_non_smoking = 0\n",
    "v = {0: v_non_smoking, 1: v_smoking}\n",
    "log_probability = loglogit(v, None, Smoking)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7fda90c",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ed31bd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_dynamic = estimate(log_probability, 'dynamic_model', database)\n",
    "display(r_dynamic)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb41e989",
   "metadata": {},
   "source": [
    "### Comparison of the estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72674173",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = pd.concat(\n",
    "    [\n",
    "        true_parameters['Value'],\n",
    "        r_static['Value'],\n",
    "        r_dynamic['Value'],\n",
    "    ],\n",
    "    axis='columns',\n",
    ")\n",
    "summary.columns = ['True', 'Static', 'Dynamic']\n",
    "summary.fillna('')\n",
    "display(summary)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71927e45",
   "metadata": {},
   "source": [
    "The introduction of the lag variable has increased a lot the final log likelihood from `-57529.39` to `-2226.1`.\n",
    "Note that the error term in the model is not the same as in the true model. Indeed, serial correlation has been\n",
    "ignored. Therefore, the coefficients cannot be directly compared. But their ratio can be compared, as it cancels the\n",
    "scale. And it can be seen that the ratios are almost the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0dabc9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio_true = summary.loc['coef_price', 'True'] / summary.loc['beta_last_year', 'True']\n",
    "print(f'Ratio coef_price / beta_last_year (True): {ratio_true:.3g}')\n",
    "ratio_dynamic = (\n",
    "    summary.loc['coef_price', 'Dynamic'] / summary.loc['beta_last_year', 'Dynamic']\n",
    ")\n",
    "print(f'Ratio coef_price / beta_last_year (Dynamic): {ratio_dynamic:.3g}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42b7e4cc",
   "metadata": {},
   "source": [
    "## Static model with serial correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57dd8910",
   "metadata": {},
   "source": [
    "We now introduce the agent effect to capture serial correlation. First, we tell Biogeme that the data is organized\n",
    "as a panel, meaning that there are several observations corresponding to the same individuals.\n",
    "Therefore, instead of considering that there is a sample of 11000 independent observations, Biogeme knows that there\n",
    "is actually a sample of 1000 individuals, for which a trajectory is observed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ffb56d7",
   "metadata": {},
   "source": [
    "Declaring the panel nature of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f02ff2fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "database.panel('Id')\n",
    "print(\n",
    "    f'Sample size, accounting for the panel nature of the data: {database.num_rows()}'\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d88cd6b",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bbb2cc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cte_std = Beta('cte_std', 10, None, None, 0)\n",
    "cte = cte_mean + cte_std * Draws('agent', 'NORMAL_ANTI')\n",
    "v_smoking = 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": "109ef8b8",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e244205",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_serial_static = estimate(log_probability, 'static_model_serial', database)\n",
    "display(r_serial_static)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a1ea122",
   "metadata": {},
   "source": [
    "## Dynamic model with serial correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07e999b2",
   "metadata": {},
   "source": [
    "We now introduce the state dependence in the model, to make it dynamic."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "456300e5",
   "metadata": {},
   "source": [
    "Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "688aaa14",
   "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": "7e8c2d01",
   "metadata": {},
   "source": [
    "Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67beaf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_serial_dynamic = estimate(log_probability, 'dynamic_model_serial', database)\n",
    "display(r_serial_dynamic)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1247fb5",
   "metadata": {},
   "source": [
    "### Comparison of the estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3530db92",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = pd.concat(\n",
    "    [\n",
    "        true_parameters['Value'],\n",
    "        r_static['Value'],\n",
    "        r_dynamic['Value'],\n",
    "        r_serial_static['Value'],\n",
    "        r_serial_dynamic['Value'],\n",
    "    ],\n",
    "    axis='columns',\n",
    ")\n",
    "summary.columns = ['True', 'Static', 'Dynamic', 'Static + serial', 'Dynamic + serial']\n",
    "summary.fillna('')\n",
    "display(summary)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e16ff57",
   "metadata": {},
   "source": [
    "If the number of draws is sufficiently high, we observe that the parameters are quite well recovered. We can\n",
    "actually test the hypothesis that they are equal to their true value, using a $t$-test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b27d4444",
   "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_serial_dynamic.loc[param, 'Value'] - true_value) / r_serial_dynamic.loc[\n",
    "        param, 'Robust std err.'\n",
    "    ]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e915250",
   "metadata": {},
   "source": [
    "`coef_price`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6523a894",
   "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": "1e991700",
   "metadata": {},
   "source": [
    "`beta_last_year`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7909768b",
   "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": "bb731b57",
   "metadata": {},
   "source": [
    "`cte_mean`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd32dc3a",
   "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": "033c48de",
   "metadata": {},
   "source": [
    "`cte_std`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6fefbd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "the_test = t_test('cte_std', 50)\n",
    "print(f'Test for cte_std: {the_test:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4458cdc",
   "metadata": {},
   "source": [
    "Each of this $t$-test is sufficiently low so that the hypothesis cannot be rejected."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}