{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77b94cf9",
   "metadata": {},
   "source": [
    "# TP5 — Optimisation numérique dans DICE\n",
    "\n",
    "**Séance du 16 septembre 2026 — M2 Actuariat, Sorbonne Université**\n",
    "\n",
    "Ce notebook traite de la tarification du carbone et du contrôle optimal. Nous passons des scénarios — ce qui pourrait arriver — au problème du planificateur social — ce qui devrait arriver — puis calculons numériquement le coût social du carbone (CSC).\n",
    "\n",
    "**Fichier de modèle :** le notebook charge `toolbox/DICE.py`, fourni avec la séance.\n",
    "\n",
    "**Programme**\n",
    "1. Charger le modèle et examiner les états et les contrôles.\n",
    "2. Définir l’objectif du planificateur et l’approximation à horizon fini.\n",
    "3. Résoudre le modèle avec une fenêtre glissante et lire le CSC.\n",
    "4. Étudier la sensibilité aux dommages et exporter les résultats.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02d86b13",
   "metadata": {},
   "source": [
    "## 1) Préparation : installer les dépendances et charger DICE\n",
    "\n",
    "Exécutez la cellule de préparation, puis les cellules suivantes. Les paquets pip manquants sont installés ; `DICE.py` est fourni dans `toolbox/`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "athena-dependency-setup",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Préparation commune des TP Python\n",
    "# Le notebook utilise uniquement les ressources fournies dans le dossier codes/.\n",
    "from pathlib import Path\n",
    "import importlib.util\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "REQUIRED_PACKAGES = {'numpy': 'numpy', 'pandas': 'pandas', 'matplotlib': 'matplotlib', 'scipy': 'scipy', 'numba': 'numba', 'tqdm': 'tqdm'}\n",
    "missing = [package for module, package in REQUIRED_PACKAGES.items()\n",
    "           if importlib.util.find_spec(module) is None]\n",
    "if missing:\n",
    "    print(\"Installation des paquets manquants :\", \", \".join(missing))\n",
    "    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", *missing])\n",
    "\n",
    "WORKDIR = Path.cwd().resolve()\n",
    "CODES_DIR = WORKDIR.parent if WORKDIR.name == \"correction\" else WORKDIR\n",
    "TOOLBOX_DIR = CODES_DIR / \"toolbox\"\n",
    "DATA_DIR = CODES_DIR / \"data\"\n",
    "if not CODES_DIR.is_dir():\n",
    "    raise FileNotFoundError(\"Exécutez le notebook depuis le dossier codes/.\" )\n",
    "if TOOLBOX_DIR.is_dir():\n",
    "    sys.path.insert(0, str(TOOLBOX_DIR))\n",
    "\n",
    "DICE_FILE = TOOLBOX_DIR / \"DICE.py\"\n",
    "if not DICE_FILE.is_file():\n",
    "    raise FileNotFoundError(\n",
    "        f\"Module du cours introuvable : {DICE_FILE}. \"\n",
    "        \"Téléchargez le dossier codes complet, avec toolbox/.\"\n",
    "    )\n",
    "\n",
    "print(\"Environnement prêt :\", CODES_DIR)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46864897",
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib.util\n",
    "\n",
    "if importlib.util.find_spec(\"DICE\") is None:\n",
    "    raise ModuleNotFoundError(\n",
    "        \"DICE.py est introuvable. Réexécutez la cellule de préparation.\"\n",
    "    )\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from DICE import (\n",
    "    Params,\n",
    "    init_states,\n",
    "    update_path,\n",
    "    mat_to_df,\n",
    "    obj_fun,\n",
    "    run_optimal_policy,\n",
    ")\n",
    "\n",
    "p = Params()\n",
    "print(\"Paramètres DICE chargés ; nombre de périodes :\", p.nT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90998489",
   "metadata": {},
   "source": [
    "## 2) Structure du modèle\n",
    "\n",
    "- **États :** capital $K_t$, productivité $A_t$, population $L_t$, stocks de carbone $M_t^{AT},M_t^{UP},M_t^{LO}$ et températures $T_t^{AT},T_t^{LO}$.\n",
    "- **Contrôles :** taux d’épargne $s_t$ et taux de réduction des émissions $\\mu_t$ ; ce dernier détermine implicitement le prix du carbone.\n",
    "- **Unités :** $C_t$ est en milliers de milliards de dollars et $L_t$ en millions d’habitants, donc $c_t=1000C_t/L_t$ est en milliers de dollars par habitant.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a7b46f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Paramètres instantanés et construire un chemin initial\n",
    "p = Params()\n",
    "sim = init_states(p)\n",
    "\n",
    "# Contrôles de référence pour l'initialisation (pas optimaux): petite réduction, économie constante\n",
    "sim[:, p.i_mu] = 0.03\n",
    "sim[:, p.i_s]  = 0.20\n",
    "\n",
    "# Mise à jour du chemin pour toutes les variables endogènes compte tenu des contrôles\n",
    "timevec = range(1, p.nT)\n",
    "sim = update_path(sim, timevec, p)\n",
    "\n",
    "# Peek at available columns\n",
    "print(\"Columns:\", p.col)\n",
    "print(\"Longueur de l’horizon (périodes) :\", p.nT, \" — années:\", sim[-1, p.i_time])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3c6f853",
   "metadata": {},
   "source": [
    "## 3) Objectif du planificateur\n",
    "\n",
    "Nous maximisons le bien-être social actualisé\n",
    "$$\n",
    "W=\\sum_{i=0}^{I^\\star}\\beta_\\Delta^iL_{t+i\\Delta}u\\!\\left(1000\\frac{C_{t+i\\Delta}}{L_{t+i\\Delta}}\\right),\\qquad \\beta_\\Delta=(1+\\rho)^{-\\Delta},\n",
    "$$\n",
    "avec une utilité CRRA $u(c)=(c^{1-\\gamma}-1)/(1-\\gamma)$. Le module fournit la fonction `obj_fun`.\n",
    "\n",
    "L’horizon infini est remplacé par un horizon fini $I^\\star$ dont la contribution résiduelle est inférieure à une tolérance.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "714a7150",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Examiner l'objectif sur un chemin de contrôle constant (contrôle de la sanité)\n",
    "# Ici, nous 'flattons' un contrôle (par exemple, l'épargne seulement) et évaluons le bien-être.\n",
    "x_const = np.full(p.nT-1, 0.20)  # constant saving\n",
    "W_val = obj_fun(x_const, sim, timevec, p, [p.i_s])\n",
    "print(\"Objective value (negative welfare) on constant s=0.20 path:\", W_val)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b23158dc",
   "metadata": {},
   "source": [
    "## 4) Fenêtre de planification à horizon fini $I^\\star$\n",
    "\n",
    "La fenêtre est choisie de façon que le poids d’actualisation du dernier terme soit inférieur à `p.toly`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca736d64",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Calculer une fenêtre de planification analogue à celle de DICE.run_optimal_policy\n",
    "disc, Tplanner = 1.0, 1\n",
    "while disc > p.toly:\n",
    "    Tplanner += 1\n",
    "    disc *= (1.0/(1.0 + p.rho))**p.Delta\n",
    "\n",
    "print(\"Fenêtre de planification retenue (périodes) :\", Tplanner, \" (each period = Δ années =\", p.Delta, \")\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "281626cc",
   "metadata": {},
   "source": [
    "## 5) Premier problème : optimiser seulement l’épargne\n",
    "\n",
    "Nous optimisons $s_t$ dans ses bornes et maintenons $\\mu_t$ constant. Cette solution constitue la référence sans transition.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ad344e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "print(\"Résolution : trajectoire optimale de l’épargne...\")\n",
    "bounds_s   = (p.s_lower, p.s_upper)\n",
    "control_id = [p.i_s]\n",
    "path_opt_s = run_optimal_policy(sim.copy(), timevec, p, bounds_s, control_id)\n",
    "\n",
    "print(\"Done. Preview of s_t path (first 5):\", path_opt_s[:5, p.i_s])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa566dc8",
   "metadata": {},
   "source": [
    "## 6) Optimiser l’épargne et la réduction des émissions\n",
    "\n",
    "Nous optimisons maintenant $(s_t,\\mu_t)$. Le prix implicite du carbone est stocké dans la colonne `Tax`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c9a760",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "print(\"Résolution : trajectoire optimale jointe de (s_t, μ_t)...\")\n",
    "bounds_smu   = [(p.s_lower, p.s_upper), (0.0, 1.0)]\n",
    "control_id   = [p.i_s, p.i_mu]\n",
    "path_opt_smu = run_optimal_policy(sim.copy(), timevec, p, bounds_smu, control_id)\n",
    "\n",
    "print(\"Done. Preview of μ_t path (first 5):\", path_opt_smu[:5, p.i_mu])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91840d6f",
   "metadata": {},
   "source": [
    "## 7) Visualiser les contrôles et le prix du carbone\n",
    "\n",
    "Tracez séparément l’épargne optimale $s_t$, le taux de réduction $\\mu_t$ et le prix implicite du carbone `Tax`, en USD/tCO₂.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1c5aef9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taux d'épargne du lot\n",
    "plt.figure(figsize=(6,3.8))\n",
    "plt.plot(path_opt_smu[:, p.i_time], path_opt_smu[:, p.i_s], label=\"Optimal transition (s_t)\")\n",
    "plt.plot(path_opt_s[:,  p.i_time], path_opt_s[:,  p.i_s],  linestyle=\"--\", label=\"No transition (s_t)\")\n",
    "plt.xlabel(\"Année\"); plt.ylabel(\"Saving rate s_t\"); plt.title(\"Optimal saving rate\")\n",
    "plt.legend(); plt.grid(True); plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c921bb0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taux de réduction des parcelles\n",
    "plt.figure(figsize=(6,3.8))\n",
    "plt.plot(path_opt_smu[:, p.i_time], path_opt_smu[:, p.i_mu], label=\"Optimal transition (μ_t)\")\n",
    "plt.xlabel(\"Année\"); plt.ylabel(\"Abatement rate μ_t\"); plt.title(\"Optimal abatement rate\")\n",
    "plt.legend(); plt.grid(True); plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9245f84b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taxe sur le carbone implicite (USD par tCO2)\n",
    "plt.figure(figsize=(6,3.8))\n",
    "plt.plot(path_opt_smu[:, p.i_time], path_opt_smu[:, p.i_Tax], label=\"Implied carbon tax (USD/tCO₂)\")\n",
    "plt.xlabel(\"Année\"); plt.ylabel(\"Carbon tax (USD/tCO₂)\"); plt.title(\"Optimal carbon price / SCC (model-implied)\")\n",
    "plt.legend(); plt.grid(True); plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1efad3ca",
   "metadata": {},
   "source": [
    "## 8) Lire le coût social du carbone\n",
    "\n",
    "- La colonne `Tax` contient le prix optimal du carbone, en USD/tCO₂.\n",
    "- Le CSC est le prix d’ombre des émissions dans le problème du planificateur. Il peut aussi être estimé par une petite perturbation des émissions puis conversion de la perte marginale de bien-être en dollars.\n",
    "\n",
    "Le squelette ci-dessous illustre cette seconde méthode ; il reste commenté pour réduire le temps d’exécution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5170e76",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# -- Optional: perturbation-based SCC at a chosen date (slow; left as a template) --\n",
    "# t_idx = 10         # pick a period index\n",
    "# eps = 1e-3 # ajouter 0,001 GtC (ajuster les unités au besoin) aux émissions à t idx\n",
    "# \n",
    "# base = path_opt_smu.copy()\n",
    "# pert = path_opt_smu.copy()\n",
    "# \n",
    "# # Injecter une petite augmentation des émissions à t idx en réduisant la réduction d'un cran\n",
    "# # (Ici comme une illustration; dans une expérience complète vous modifieriez E t directement ou ajouteriez au cycle du carbone)\n",
    "# pert[t_idx, p.i_mu] = max(0.0, pert[t_idx, p.i_mu] - 1e-5)\n",
    "# \n",
    "# base_W = obj_fun(np.hstack([base[1:, p.i_s], base[1:, p.i_mu]]), base, range(1, p.nT), p, [p.i_s, p.i_mu])\n",
    "# pert_W = obj_fun(np.hstack([pert[1:, p.i_s], pert[1:, p.i_mu]]), pert, range(1, p.nT), p, [p.i_s, p.i_mu])\n",
    "# dW = pert_W - base_W  # (negative welfare) difference\n",
    "# print(\"ΔW de ε-perturbation:\", dW)\n",
    "# # Convertir en USD/tCO2 en utilisant l'utilitaire marginal à t idx si nécessaire.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac5e356",
   "metadata": {},
   "source": [
    "## 9) Sensibilité à des dommages plus élevés\n",
    "\n",
    "Modifiez la fonction de dommages avec `p.user_damage_fn`, puis relancez l’optimisation pour observer la réaction de $\\mu_t$ et du prix du carbone.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8ac967e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def higher_damage(T, p_local: Params):\n",
    "    # Slightly steeper quadratic damages as an illustration\n",
    "    return 0.0030 * (np.asarray(T) ** 2)\n",
    "\n",
    "p_sens = Params()\n",
    "sim_sens = init_states(p_sens)\n",
    "sim_sens[:, p_sens.i_mu] = 0.03\n",
    "sim_sens[:, p_sens.i_s]  = 0.20\n",
    "sim_sens = update_path(sim_sens, range(1, p_sens.nT), p_sens)\n",
    "\n",
    "p_sens.user_damage_fn = higher_damage\n",
    "path_opt_smu_highDam = run_optimal_policy(sim_sens.copy(), range(1, p_sens.nT), p_sens, [(p_sens.s_lower,p_sens.s_upper),(0.0,1.0)], [p_sens.i_s, p_sens.i_mu])\n",
    "\n",
    "# Comparaison par parcelle pour μ t\n",
    "plt.figure(figsize=(6,3.8))\n",
    "plt.plot(path_opt_smu[:, p.i_time], path_opt_smu[:, p.i_mu], linestyle=\"--\", label=\"Baseline damages\")\n",
    "plt.plot(path_opt_smu_highDam[:, p.i_time], path_opt_smu_highDam[:, p.i_mu], label=\"Higher damages\")\n",
    "plt.xlabel(\"Année\"); plt.ylabel(\"Abatement rate μ_t\"); plt.title(\"Damage sensitivity: optimal abatement\")\n",
    "plt.legend(); plt.grid(True); plt.tight_layout()\n",
    "\n",
    "# Comparaison par parcelle pour la taxe implicite sur le carbone\n",
    "plt.figure(figsize=(6,3.8))\n",
    "plt.plot(path_opt_smu[:, p.i_time], path_opt_smu[:, p.i_Tax], linestyle=\"--\", label=\"Baseline damages\")\n",
    "plt.plot(path_opt_smu_highDam[:, p.i_time], path_opt_smu_highDam[:, p.i_Tax], label=\"Higher damages\")\n",
    "plt.xlabel(\"Année\"); plt.ylabel(\"Carbon tax (USD/tCO₂)\"); plt.title(\"Damage sensitivity: optimal carbon price\")\n",
    "plt.legend(); plt.grid(True); plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cef164da",
   "metadata": {},
   "source": [
    "## 10) Exporter les résultats\n",
    "\n",
    "Convertissez les tableaux en `DataFrame` pour poursuivre l’analyse ou produire d’autres graphiques.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98fa4a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "df_baseline = pd.DataFrame({\n",
    "    \"year\": path_opt_smu[:, p.i_time],\n",
    "    \"s\":    path_opt_smu[:, p.i_s],\n",
    "    \"mu\":   path_opt_smu[:, p.i_mu],\n",
    "    \"tax\":  path_opt_smu[:, p.i_Tax],\n",
    "    \"E\":    path_opt_smu[:, p.i_E],\n",
    "    \"T_AT\": path_opt_smu[:, p.i_T_AT],\n",
    "    \"C\":    path_opt_smu[:, p.i_C],\n",
    "})\n",
    "df_baseline.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dfc29bf",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "\n",
    "- Le prix optimal `Tax` constitue la recommandation normative du modèle.\n",
    "- La fenêtre finie met en œuvre la troncation temporelle.\n",
    "- La résolution par fenêtre glissante produit les trajectoires complètes des états et des contrôles.\n",
    "- Les hypothèses de dommages modifient directement $\\mu_t$ et le prix implicite du carbone.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}