{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "606982f9",
   "metadata": {},
   "source": [
    "\n",
    "File: 01-testing_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Sat Aug 02 2025, 17:42:44\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc1ad2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.expressions import Beta\n",
    "from biogeme.models import loglogit\n",
    "from biogeme.results_processing import (\n",
    "    EstimationResults,\n",
    "    calc_p_value,\n",
    "    compile_estimation_results,\n",
    "    get_pandas_correlation_results,\n",
    "    get_pandas_estimated_parameters,\n",
    ")\n",
    "from scipy.stats import chi2, norm\n",
    "\n",
    "from logit_airline_alt_spec import (\n",
    "    Opt1 as Opt1_alt_spec,\n",
    "    Opt2 as Opt2_alt_spec,\n",
    "    Opt3 as Opt3_alt_spec,\n",
    "    results as res_alt_spec,\n",
    ")\n",
    "from logit_airline_base import (\n",
    "    Fare_1,\n",
    "    Fare_2,\n",
    "    Fare_3,\n",
    "    chosenAlternative,\n",
    "    database,\n",
    "    results as res_base,\n",
    ")\n",
    "from logit_airline_boxcox import results as res_boxcox\n",
    "from logit_airline_income import (\n",
    "    Opt1 as Opt1_income,\n",
    "    Opt2 as Opt2_income,\n",
    "    Opt3 as Opt3_income,\n",
    "    results as res_income,\n",
    ")\n",
    "from logit_airline_piecewise import (\n",
    "    Opt1 as Opt1_piecewise,\n",
    "    Opt2 as Opt2_piecewise,\n",
    "    Opt3 as Opt3_piecewise,\n",
    "    results as res_piecewise,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0709c66e",
   "metadata": {},
   "source": [
    "The objective of this laboratory is to investigate various specification tests."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35cdd3a7",
   "metadata": {},
   "source": [
    "Base model:\n",
    "\\begin{align*}\n",
    "V_class_1 &= \\beta_\\text{fare}  \\text{fare}_1 + \\beta_\\text{legroom}  \\text{legroom}_1\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_1  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_1\n",
    "+ \\beta_\\text{time} \\text{elapsed\\_time}_1 \\\\\n",
    "V_class_2 &= \\text{cte}_2 + \\beta_\\text{fare}  \\text{fare}_2 + \\beta_\\text{legroom}  \\text{legroom}_2\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_2  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_2\n",
    "+ \\beta_\\text{time} \\text{elapsed\\_time}_2 \\\\\n",
    "V_3 &= \\text{cte}_3 + \\beta_\\text{fare}  \\text{fare}_3 + \\beta_\\text{legroom}  \\text{legroom}_3\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_3  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_3\n",
    "+ \\beta_\\text{time} \\text{elapsed\\_time}_3 \\\\\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "139a1050",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_base.short_summary())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a48620d",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(get_pandas_estimated_parameters(estimation_results=res_base))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "503eee1b",
   "metadata": {},
   "source": [
    "Alternative specific time coefficient:\n",
    "\\begin{align*}\n",
    "V_class_1 &= \\beta_\\text{fare}  \\text{fare}_1 + \\beta_\\text{legroom}  \\text{legroom}_1\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_1  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_1\n",
    "+ \\beta_\\text{time, 1} \\text{elapsed\\_time}_1 \\\\\n",
    "V_class_2 &= \\text{cte}_2 + \\beta_\\text{fare}  \\text{fare}_2 + \\beta_\\text{legroom}  \\text{legroom}_2\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_2  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_2\n",
    "+ \\beta_\\text{time, 2} \\text{elapsed\\_time}_2 \\\\\n",
    "V_3 &= \\text{cte}_3 + \\beta_\\text{fare}  \\text{fare}_3 + \\beta_\\text{legroom}  \\text{legroom}_3\n",
    "+ \\beta_\\text{sd\\_early} \\text{sched\\_delay\\_early}_3  + \\beta_\\text{sd\\_late} \\text{sched\\_delay\\_late}_3\n",
    "+ \\beta_\\text{time, 2} \\text{elapsed\\_time}_3 \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c43e2a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_alt_spec.short_summary())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb92b2e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_alt_spec_pandas = get_pandas_estimated_parameters(\n",
    "    estimation_results=res_alt_spec\n",
    ")\n",
    "display(parameters_alt_spec_pandas)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff66011",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_alt_spec = res_alt_spec.get_beta_values()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65377c73",
   "metadata": {},
   "source": [
    "# Question 1: test the null hypothesis that `beta_elapsed_time_2` is equal to `beta_elapsed_time_3`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aee48964",
   "metadata": {},
   "source": [
    "Under the null hypothesis, the statistic\n",
    "$$\n",
    "\\frac{\\widehat{\\beta}_2 - \\widehat{\\beta}_3}{\\sqrt{\\operatorname{Var}(\\widehat{\\beta}_2 - \\widehat{\\beta}_3)}},\n",
    "$$\n",
    "is distributed as a $N(0, 1)$, where\n",
    "$$\n",
    "\\operatorname{Var}(\\widehat{\\beta}_2 - \\widehat{\\beta}_3) = \\operatorname{Var}(\\widehat{\\beta}_2) +\n",
    "\\operatorname{Var}(\\widehat{\\beta}_3) - 2 \\operatorname{Covar}(\\widehat{\\beta}_2, \\widehat{\\beta}_3).\n",
    "$$\n",
    "Note that it is called a $t$-statistic, as in theory, it follows a $t$ distribution. But with a large sample size,\n",
    "it is very well approximated by the normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dde61b1",
   "metadata": {},
   "source": [
    "We first calculate the numerator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed16a214",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = (\n",
    "    parameters_alt_spec['beta_elapsed_time_2']\n",
    "    - parameters_alt_spec['beta_elapsed_time_3']\n",
    ")\n",
    "print(f'{numerator=:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea6241b",
   "metadata": {},
   "source": [
    "For the denominator, we need the covariance between the estimates. It is available in the following table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1759435b",
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation_alt_spec = get_pandas_correlation_results(estimation_results=res_alt_spec)\n",
    "display(correlation_alt_spec)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc9b6dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "covar_beta2_beta3 = correlation_alt_spec.loc[\n",
    "    (correlation_alt_spec['First parameter'] == 'beta_elapsed_time_3')\n",
    "    & (correlation_alt_spec['Second parameter'] == 'beta_elapsed_time_2'),\n",
    "    'Robust covariance',\n",
    "].values[0]\n",
    "print(f'{covar_beta2_beta3=:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14fde2bf",
   "metadata": {},
   "source": [
    "Therefore, the square of the denominator is calculated as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a30a370",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_err_beta_elapsed_time_2 = res_alt_spec.get_parameter_std_err(\n",
    "    parameter_name='beta_elapsed_time_2'\n",
    ")\n",
    "std_err_beta_elapsed_time_3 = res_alt_spec.get_parameter_std_err(\n",
    "    parameter_name='beta_elapsed_time_3'\n",
    ")\n",
    "squared_denominator = (\n",
    "    std_err_beta_elapsed_time_2**2\n",
    "    + std_err_beta_elapsed_time_3**2\n",
    "    - 2 * covar_beta2_beta3\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0319502a",
   "metadata": {},
   "source": [
    "And the $t$-statistics is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88fe534",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_stat = numerator / (squared_denominator**0.5)\n",
    "print(f't-statistic: {t_stat:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19392f99",
   "metadata": {},
   "source": [
    "Note that this value is immediately available in the table `correlation_alt_spec` generated by the estimation\n",
    "procedure. Note also that the sign of the value is irrelevant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80f8e240",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_test_beta_2_beta_3 = correlation_alt_spec.loc[\n",
    "    (correlation_alt_spec['First parameter'] == 'beta_elapsed_time_3')\n",
    "    & (correlation_alt_spec['Second parameter'] == 'beta_elapsed_time_2'),\n",
    "    'Robust t-stat.',\n",
    "].values[0]\n",
    "print(f'{covar_beta2_beta3=:.3g}')\n",
    "\n",
    "print(f't-test: {t_test_beta_2_beta_3:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4972f29",
   "metadata": {},
   "source": [
    "Threshold for the 10% confidence level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33aacebb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'{norm.ppf(0.9):.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31b5f672",
   "metadata": {},
   "source": [
    "Given the low value of the $t$-statistic, we cannot reject the null hypothesis that `beta_elapsed_time_2` is equal\n",
    "to `beta_elapsed_time_3` at 10% confidence. Actually, as the high $p$-value suggests, the null hypothesis cannot\n",
    "be rejected here for any reasonable confidence level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52bfd66e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'{100*calc_p_value(t_test_beta_2_beta_3):.3g}%')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1d45465",
   "metadata": {},
   "source": [
    "If we reject the null hypothesis, there is 94% probability to make an error."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa177c7",
   "metadata": {},
   "source": [
    "# Question 2: test the null hypothesis that, between the two models, the base model is the true model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1160aa2",
   "metadata": {},
   "source": [
    "The base model is the restricted model. It is obtained from the unrestricted model using the linear restrictions:\n",
    "$$\\beta_{t,1} = \\beta_{t, 2} = \\beta_{t, 3}.$$ Therefore, we can perform a likelihood ratio test.\n",
    "Under the null hypothesis, the statistic\n",
    "$$ -2(\\mathcal{L}(\\widehat{\\beta}_R) - \\mathcal{L}(\\widehat{\\beta}_U))$$\n",
    "is distributed as $\\chi^2$ with $K_U - K_R$ degrees of freedom."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36dfe1e",
   "metadata": {},
   "source": [
    "We collect the general statistics for the restricted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0563f738",
   "metadata": {},
   "outputs": [],
   "source": [
    "general_statistics_base = res_base.get_general_statistics()\n",
    "display(general_statistics_base)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91e8f3bc",
   "metadata": {},
   "source": [
    "We collect the general statistics for the unrestricted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "471f6b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "general_statistics_alt_spec = res_alt_spec.get_general_statistics()\n",
    "display(general_statistics_alt_spec)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62e63871",
   "metadata": {},
   "source": [
    "We calculate the statistic for the test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30c05a41",
   "metadata": {},
   "outputs": [],
   "source": [
    "LR = res_base.final_loglikelihood\n",
    "LU = res_alt_spec.final_loglikelihood\n",
    "test = -2 * (LR - LU)\n",
    "print(f'Likelihood ratio statistic: {test:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a18e1527",
   "metadata": {},
   "source": [
    "We calculate the number of degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e86ad506",
   "metadata": {},
   "outputs": [],
   "source": [
    "KR = res_base.number_of_free_parameters\n",
    "KU = res_alt_spec.number_of_free_parameters\n",
    "degrees_of_freedom = KU - KR\n",
    "print(f'Degrees of freedom: {degrees_of_freedom}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "558cf6b5",
   "metadata": {},
   "source": [
    "What is the threshold for the statistics at a 5% significance level?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "281ca47d",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = chi2.ppf(0.95, degrees_of_freedom)\n",
    "print(f'Threshold for the test if 5%: {threshold:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d2cfebb",
   "metadata": {},
   "source": [
    "As the value of the statistic is below the threshold, the null hypothesis cannot be rejected at the 5% level."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89847d3",
   "metadata": {},
   "source": [
    "There is a tool that allows to perform directly this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc43f09",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_result = res_alt_spec.likelihood_ratio_test(res_base, 0.05)\n",
    "print(f'{lr_result.statistic=:.3g}')\n",
    "print(f'{lr_result.threshold=:.3g}')\n",
    "print(lr_result.message)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "603edb64",
   "metadata": {},
   "source": [
    "# Question 3: Test the null hypothesis that the fare coefficient does not vary with income."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c791b6c7",
   "metadata": {},
   "source": [
    "Base model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22d57080",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_base.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66c3f83a",
   "metadata": {},
   "source": [
    "Model where the fare coefficient varies with income."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48379c82",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_income.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc9cd8e4",
   "metadata": {},
   "source": [
    "Under the null hypothesis, the true model is the base model. Note that we cannot perform a likelihood ratio test\n",
    "here, as the base model cannot be obtained from a linear restriction of the parameters of the unrestricted model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac305b78",
   "metadata": {},
   "source": [
    "If we compare the results, we note that the model where the fare coefficient varies with income has a poorer fit\n",
    "than the base model, although it has one more parameter. It seems to suggest that the base model should be preferred.\n",
    "In order to formally test it, we perform a Cox test."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a43d8ca",
   "metadata": {},
   "source": [
    "Composite model: we add the term involving cost from the base model to each utility function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba62d7ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_fare = Beta('beta_fare', 0, None, None, 0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58afec7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Opt1_coxtest = Opt1_income + beta_fare * Fare_1\n",
    "Opt2_coxtest = Opt2_income + beta_fare * Fare_2\n",
    "Opt3_coxtest = Opt3_income + beta_fare * Fare_3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42433ab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = {1: Opt1_coxtest, 2: Opt2_coxtest, 3: Opt3_coxtest}\n",
    "logprob_coxtest = loglogit(V, None, chosenAlternative)\n",
    "biogeme_coxtest = BIOGEME(database, logprob_coxtest)\n",
    "biogeme_coxtest.model_name = 'logit_airline_coxtest'\n",
    "results_coxtest: EstimationResults = biogeme_coxtest.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "132c11dd",
   "metadata": {},
   "source": [
    "General statistics of the composite model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3941820",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_coxtest.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5268108d",
   "metadata": {},
   "source": [
    "Estimated parameters of the composite model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82eb6b51",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(get_pandas_estimated_parameters(estimation_results=results_coxtest))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "739ca9d2",
   "metadata": {},
   "source": [
    "We perform an LR test to compare the composite model (unrestricted) with the base model (restricted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43961b1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_composite_base = results_coxtest.likelihood_ratio_test(res_base)\n",
    "print(f'{lr_composite_base.statistic=:.3g}')\n",
    "print(f'{lr_composite_base.threshold=:.3g}')\n",
    "print(lr_composite_base.message)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba3d7593",
   "metadata": {},
   "source": [
    "We perform an LR test to compare the composite model (unrestricted) with the \"income\" model (restricted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8aeefea",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_composite_income = results_coxtest.likelihood_ratio_test(res_income)\n",
    "print(f'{lr_composite_income.statistic=:.3g}')\n",
    "print(f'{lr_composite_income.threshold=:.3g}')\n",
    "print(lr_composite_income.message)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f259dcee",
   "metadata": {},
   "source": [
    "Both models are rejected against the composite model. It shows that none of the two models can be accepted and that\n",
    "better models should be investigated. Note that it does not mean that the composite model is the correct one.\n",
    "This model has been specified for the sake of the test, and does not necessarily correspond to a behaviorally\n",
    "meaningful specification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "558cdc95",
   "metadata": {},
   "source": [
    "# Question 4: test the linear specification with alternative specific time parameters, against the nonlinear specification with Box-Cox transform of the time variable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e56c303",
   "metadata": {},
   "source": [
    "Alternative specific time parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34816d24",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_alt_spec.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f06b0f8",
   "metadata": {},
   "source": [
    "Box-Cox model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a55544ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_boxcox.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0f12483",
   "metadata": {},
   "source": [
    "There are two ways to test the null hypothesis that the true model is the linear model:\n",
    "\n",
    "- using a $t$-test,\n",
    "- using a likelihood ratio test."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "757401e9",
   "metadata": {},
   "source": [
    "Under the null hypothesis, the true value of $\\lambda$ is 1. Therefore, we perform a $t$-test to test that\n",
    "hypothesis. Note that the $t$-test reported in the results is to test if the true value of the parameter is zero.\n",
    "Therefore, it cannot be used here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7724150",
   "metadata": {},
   "outputs": [],
   "source": [
    "boxcox_estimates = get_pandas_estimated_parameters(estimation_results=res_boxcox)\n",
    "display(boxcox_estimates)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e38be0b",
   "metadata": {},
   "source": [
    "Under the null hypothesis, the statistic $$\\frac{\\widehat{\\lambda}-1}{\\widehat{\\sigma}_\\lambda}$$ follows\n",
    "approximately a $N(0, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc6bf3c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_value = res_boxcox.get_parameter_value('lambda_boxcox')\n",
    "lambda_std_err = res_boxcox.get_parameter_std_err('lambda_boxcox')\n",
    "t_test_lambda = (lambda_value - 1) / lambda_std_err\n",
    "print(f't-test testing lambda=1: {t_test_lambda:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f13d479b",
   "metadata": {},
   "source": [
    "The threshold for the 10% confidence level is the 95% quantile of the normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1f5acaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'95% quantile of the normal distribution: {norm.ppf(.95):.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bff7a5c",
   "metadata": {},
   "source": [
    "Therefore, the null hypothesis cannot be rejected at the 10% level of confidence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd487b24",
   "metadata": {},
   "source": [
    "We can calculate the p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cece6eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'p-value: {calc_p_value(t_test_lambda):.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06c8bf34",
   "metadata": {},
   "source": [
    "If we reject the null hypothesis that the true model is the linear model, there is a 23.5% probability to make\n",
    "an error."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf51ac6a",
   "metadata": {},
   "source": [
    "We can also perform a likelihood ratio test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "256fda7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_boxbox_alt_spec = res_boxcox.likelihood_ratio_test(res_alt_spec, 0.10)\n",
    "print(f'{lr_boxbox_alt_spec.statistic=:.3g}')\n",
    "print(f'{lr_boxbox_alt_spec.threshold=:.3g}')\n",
    "print(lr_boxbox_alt_spec.message)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ec06815",
   "metadata": {},
   "source": [
    "We reach the same conclusion. The null hypothesis cannot be rejected at the 10% level."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c33af9a",
   "metadata": {},
   "source": [
    "# Test the piecewise linear specification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ca1c499",
   "metadata": {},
   "source": [
    "First, we test the piecewise linear specification against the base model. To do so, we can perform a likelihood\n",
    "ratio test. Indeed, the base model can be obtained from the piecewise linear specification by imposing the linear\n",
    "constraints that all the time coefficients `beta_time_0_2`, `beta_time_2_3` and `beta_time_3_more` are equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9cf09f",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_piecewise_base = res_piecewise.likelihood_ratio_test(res_base, 0.01)\n",
    "print(f'{lr_piecewise_base.statistic=:.3g}')\n",
    "print(f'{lr_piecewise_base.threshold=:.3g}')\n",
    "print(lr_piecewise_base.message)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8049cc8b",
   "metadata": {},
   "source": [
    "It seems that the null hypothesis that the base model is the true one can be rejected at 1% level of confidence.\n",
    "Therefore, out of the two, we prefer the piecewise linear specification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19912534",
   "metadata": {},
   "source": [
    "Second, we test the piecewise linear specification against the specification with the alternative specific time\n",
    "coefficient. Both models have 9 parameters, and the log likelihood for the piecewise linear model is higher."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5a4de5",
   "metadata": {},
   "source": [
    "We test the null hypothesis that the alternative specific model is the true model. To test this hypothesis, we need\n",
    "to perform a **Davidson McKinnon test**, where $M_1$ is the alternative specific model, and $M_2$ is the piecewise\n",
    "linear model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c98cb68",
   "metadata": {},
   "source": [
    "We write a composite specification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b38ab81f",
   "metadata": {},
   "source": [
    "We fix the values of the parameters of the piecewise linear specification to their estimated value. We also change\n",
    "the name of the beta to avoid duplicates in the composite model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b257dd17",
   "metadata": {},
   "outputs": [],
   "source": [
    "piecewise_estimates = res_piecewise.get_beta_values()\n",
    "print(piecewise_estimates)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ad6b26e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Opt1_piecewise.fix_betas(beta_values=piecewise_estimates, prefix='piecewise_')\n",
    "Opt2_piecewise.fix_betas(beta_values=piecewise_estimates, prefix='piecewise_')\n",
    "Opt3_piecewise.fix_betas(beta_values=piecewise_estimates, prefix='piecewise_')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56d37693",
   "metadata": {},
   "source": [
    "We create the composite specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "361d5573",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = Beta('alpha', 0, None, None, 0)\n",
    "Opt1 = (1 - alpha) * Opt1_alt_spec + alpha * Opt1_piecewise\n",
    "Opt2 = (1 - alpha) * Opt2_alt_spec + alpha * Opt2_piecewise\n",
    "Opt3 = (1 - alpha) * Opt3_alt_spec + alpha * Opt3_piecewise\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44c9d6bf",
   "metadata": {},
   "source": [
    "We estimate the composite model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e82d7b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = {1: Opt1, 2: Opt2, 3: Opt3}\n",
    "logprob = loglogit(V, None, chosenAlternative)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'logit_airline_piecewise'\n",
    "results_composite: EstimationResults = biogeme.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4a06e85",
   "metadata": {},
   "source": [
    "Estimated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f19b852",
   "metadata": {},
   "outputs": [],
   "source": [
    "composite_parameters = get_pandas_estimated_parameters(\n",
    "    estimation_results=results_composite\n",
    ")\n",
    "display(composite_parameters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c10c3fd",
   "metadata": {},
   "source": [
    "$t$-test for $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bf2e26e",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_test_alpha = results_composite.get_parameter_t_test('alpha')\n",
    "print(f't-test for alpha: {t_test_alpha:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e1b04f5",
   "metadata": {},
   "source": [
    "$p$-value for $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83372f0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value_alpha = results_composite.get_parameter_p_value('alpha')\n",
    "print(f'p-value for alpha: {p_value_alpha:.3g}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b69b768d",
   "metadata": {},
   "source": [
    "Under the null hypothesis, the true value of $\\alpha$ is 0. This hypothesis can be safely rejected here. Therefore,\n",
    "the piecewise linear specification is preferred."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69111d0c",
   "metadata": {},
   "source": [
    "Note that we reach the same conclusion when comparing the AIC and the BIC of both models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15da97b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison, _ = compile_estimation_results(\n",
    "    {'Alt. spec': res_alt_spec, 'Piecewise linear': res_piecewise}\n",
    ")\n",
    "display(comparison)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}