{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77b94cf9",
   "metadata": {},
   "source": [
    "# Session 5 — Numerical Optimization Basics in DICE (Python)\n",
    "\n",
    "This notebook complements the slides of the **Lecture 2** on carbon pricing and optimal control. We move from scenario-building (*what could happen*) to the social planner's problem (*what should happen*), and show how the **optimal carbon price (SCC)** emerges and is computed numerically.\n",
    "\n",
    "**Roadmap**  \n",
    "1. Set up and initialize the DICE model.  \n",
    "2. Present the model and solve the planner's problem.  \n",
    "3. Exercise: compare the joint optimum with the no-mitigation optimum.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "session5-cell-1",
   "metadata": {},
   "source": [
    "## Session 5 — Objectives and prerequisites\n",
    "\n",
    "**Learning objectives**\n",
    "- Identify the original DICE planner objective, controls, and finite planning window.\n",
    "- Compare savings-only and joint savings–abatement optimization.\n",
    "- Read optimal controls and the implied carbon-price path without changing the scientific model.\n",
    "\n",
    "**Prerequisites.** Slides for Session 5; Session 3 DICE structure; constrained optimization and discounted utility.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "x3a-setup-5-markdown",
   "metadata": {
    "x3a_technical_cell": true
   },
   "source": [
    "## 1. Setup\n",
    "\n",
    "Run the notebook from this session directory. The original course module and data files are stored next to it.\n",
    "\n",
    "**Convention.** At the initial date, control variables are `NaN`: the initial state is given, but the control is chosen only for subsequent transitions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "x3a-setup-5-code",
   "metadata": {
    "x3a_technical_cell": true
   },
   "outputs": [],
   "source": [
    "# Install dependencies in the current notebook kernel before running imports.\n",
    "%pip install -r requirements.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02d86b13",
   "metadata": {},
   "source": [
    "### 1.A Setup\n",
    "We import the DICE model components and ensure Python can find the module.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46864897",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import os, sys, math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Make sure we can import the user's module\n",
    "from importlib import reload\n",
    "import DICE\n",
    "reload(DICE)\n",
    "\n",
    "from DICE import Params, init_states, update_path, obj_fun, run_optimal_policy, mat_to_df\n",
    "print(\"DICE loaded from:\", DICE.__file__)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90998489",
   "metadata": {},
   "source": [
    "## 2. Presentation\n",
    "\n",
    "### 2.A Model skeleton (states & controls)\n",
    "\n",
    "- **States (examples):** capital $K_t$, TFP $A_t$, population $L_t$, climate stocks $M_{\\mathrm{AT}}, M_{\\mathrm{UP}}, M_{\\mathrm{LO}}$, temperatures $T_{\\mathrm{AT}}, T_{\\mathrm{LO}}$.\n",
    "- **Controls:** saving rate $s_t$, abatement rate $\\mu_t$ (which, together with technology parameters, implies a **carbon price / tax** path).\n",
    "- **Units:** $C_t$ in trillion USD, $L_t$ in million → per-capita consumption in **thousand USD**: $c_t = 1000 \\times C_t/L_t$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a7b46f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Instantiate parameters and construct an initial path\n",
    "p = Params()\n",
    "sim = init_states(p)\n",
    "\n",
    "# Baseline controls for initialization (not optimal): small abatement, constant saving\n",
    "sim[1:, p.i_mu] = 0.03\n",
    "sim[1:, p.i_s]  = 0.20\n",
    "\n",
    "# Update path for all endogenous variables given the controls\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(\"Horizon length (periods):\", p.nT, \" — years:\", sim[-1, p.i_time])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3c6f853",
   "metadata": {},
   "source": [
    "### 2.B Planner's objective (CRRA utility)\n",
    "\n",
    "We maximize discounted social welfare\n",
    "$$\n",
    "W = \\sum_{i=0}^{I^\\star} \\beta_\\Delta^{\\,i}\\, L_{t+i\\Delta}\\, u\\!\\left( 1000\\,\\frac{C_{t+i\\Delta}}{L_{t+i\\Delta}} \\right), \n",
    "\\qquad \\beta_\\Delta=(1+\\rho)^{-\\Delta},\n",
    "$$\n",
    "with CRRA utility $u(c)=\\frac{c^{1-\\gamma}-1}{1-\\gamma}$. The notebook uses the `obj_fun` provided in the module.\n",
    "\n",
    "**Note on truncation:** We replace the infinite horizon with a finite $I^\\star$ chosen so the tail weight is below a tolerance.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "714a7150",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate discounted welfare with a constant saving rate of 20%.\n",
    "# obj_fun returns NEGATIVE welfare because the optimizer minimizes it.\n",
    "x_const = np.full(p.nT - 1, 0.20)\n",
    "W_val = -obj_fun(x_const, sim, timevec, p, [p.i_s])\n",
    "start_year = int(sim[timevec[0], p.i_time])\n",
    "end_year = int(sim[timevec[-1], p.i_time])\n",
    "print(f\"Discounted welfare evaluated at {start_year}, over {start_year}–{end_year}, is {W_val:.2f} (constant saving rate: 20%).\")\n",
    "# Total discounted welfare over the horizon, not utility in one period.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b23158dc",
   "metadata": {},
   "source": [
    "### 2.C Finite-horizon planning window $I^{\\star}$\n",
    "\n",
    "We select a planning window so that the discount weight of the last term is below the numerical tolerance `p.toly`. This implements the finite-horizon approximation discussed in the slides.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca736d64",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Compute planning window similar to 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(\"Chosen planning window (periods):\", Tplanner, \" (each period = Δ years =\", p.Delta, \")\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa566dc8",
   "metadata": {},
   "source": [
    "### 2.D Solve the full problem: **savings + abatement**\n",
    "\n",
    "Now we optimize the pair $(s_t, \\mu_t)$. The carbon tax path is then implied by the model and reported in the `Tax` column. The problem reads as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{\\{s_{t+i},\\,\\{\\mu_{t+i}\\}_{i=0}^{I^\\star}} \\quad \n",
    "& \\sum_{i=0}^{I^\\star} \\beta^{i}\\, L_{t+i\\Delta}\\;\n",
    "u\\!\\left( 1000\\,\\frac{C_{t+i\\Delta}}{L_{t+i\\Delta}} \\right) \\\\[0.4em]\n",
    "\\text{s.t.}\\quad \n",
    "& y_{t+i} \\;=\\; G_p\\!\\big(y_{t+i-1},\\, x_{t+i-1},\\, s_{t+i},\\, \\mu_{t+i};\\big), \n",
    "\\quad i=0,\\dots,I^\\star\\!-\\!1,\\\\\n",
    "& s_{\\min} \\;\\le\\; s_{t+i} \\;\\le\\; s_{\\max}\\ \\ (\\text{box constraints}),\\\\\n",
    "& 0 \\;\\le\\; \\mu_{t+i} \\;\\le\\; 1\\ \\ (\\text{box constraints}),\\\\\n",
    "& y_{t-1} \\;=\\; y^{\\mathrm{given}}_{t-1},\n",
    "\\end{aligned}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c9a760",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2.D — Joint optimum: the planner chooses both saving and abatement.\n",
    "sim_smu = sim.copy()\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_smu, timevec, p, bounds_smu, control_id)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91840d6f",
   "metadata": {},
   "source": [
    "### 2.E Visualize the optimal controls and the implied carbon price\n",
    "\n",
    "Below we plot (one chart per figure) the optimal saving \\(s_t\\), abatement \\(\\mu_t\\), and the implied carbon tax (USD/tCO\\(_2\\)) reported in `Tax`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1c5aef9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the joint optimum after you have completed Exercise 2.D.\n",
    "# Until then, this cell deliberately remains executable without revealing the solution.\n",
    "if \"path_opt_smu\" not in globals():\n",
    "    print(\"Complete Exercise 2.D and run its cell before producing these plots.\")\n",
    "else:\n",
    "    # Plot saving rate\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.xlabel(\"Year\"); plt.ylabel(\"Saving rate s_t\"); plt.title(\"Optimal saving rate\")\n",
    "    plt.legend(); plt.grid(True); plt.tight_layout()\n",
    "\n",
    "    # Plot abatement rate\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(\"Year\"); plt.ylabel(\"Abatement rate μ_t\"); plt.title(\"Optimal abatement rate\")\n",
    "    plt.legend(); plt.grid(True); plt.tight_layout()\n",
    "\n",
    "    # Plot implied carbon tax (USD per 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(\"Year\"); 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": "950c7fb6",
   "metadata": {},
   "source": [
    "## 3. Exercise — Optimal climate policy versus no-mitigation optimum\n",
    "\n",
    "### 3.A Joint optimal policy\n",
    "\n",
    "This exercise compares the **joint optimum**, in which the planner chooses saving and abatement, with a **constrained optimum without mitigation**, in which the planner chooses saving but abatement is fixed at zero.\n",
    "\n",
    "> This second regime is not strict laissez-faire: saving remains a planner's choice. It isolates the value of having the mitigation instrument available.\n",
    "\n",
    "Reuse the joint optimum computed above, in which the planner chooses both $s_t$ and $\\mu_t$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77f0979",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.A — Reuse the joint optimum computed in Section 6.\n",
    "# path_opt_full = ...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddfd70c7",
   "metadata": {},
   "source": [
    "### 3.B Optimal saving with abatement fixed at zero\n",
    "\n",
    "Construct the constrained optimum with no mitigation:\n",
    "\n",
    "1. start from a copy of the initial simulation;\n",
    "2. impose $\\mu_t=0$ for every transition period;\n",
    "3. remove p.i_mu from control_id, so that the optimizer chooses only saving;\n",
    "4. solve the problem and store the path as path_opt_no_mitigation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1faaa37e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.B — Constrained optimum without mitigation.\n",
    "# sim_no_mitigation = sim.copy()\n",
    "# sim_no_mitigation[1:, p.i_mu] = ...\n",
    "# path_opt_no_mitigation = run_optimal_policy(\n",
    "#     sim_no_mitigation, timevec, p, (p.s_lower, p.s_upper), [p.i_s]\n",
    "# )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2de81c9a",
   "metadata": {},
   "source": [
    "### 3.C Compare the two policy regimes\n",
    "\n",
    "Produce a three-panel figure comparing the joint optimum and the no-mitigation optimum for atmospheric temperature $T_{AT,t}$, industrial emissions $E_t$, and consumption $C_t$. Briefly explain why the two regimes differ.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5e90d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.C — Plot T_AT, E, and C for the two policy regimes.\n",
    "# Hint: use plt.subplots(1, 3, figsize=(15, 4)) and the columns\n",
    "# p.i_T_AT, p.i_E, and p.i_C.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6223c6bd",
   "metadata": {},
   "source": [
    "### 3.D Timing of mitigation costs and climate benefits\n",
    "\n",
    "Compute the consumption gap and temperature gap:\n",
    "\n",
    "$$100\\left(\\frac{C_t^{\\mathrm{joint}}}{C_t^{\\mathrm{no\\ mitigation}}}-1\\right),\n",
    "\\qquad T_t^{\\mathrm{joint}}-T_t^{\\mathrm{no\\ mitigation}}.$$\n",
    "\n",
    "Plot both gaps and comment: when does the consumption cost appear, and when does the climate benefit become visible?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3f0f771",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.D — Compute and plot the consumption and temperature gaps.\n",
    "# consumption_gap = 100 * (\n",
    "#     path_opt_full[:, p.i_C] / path_opt_no_mitigation[:, p.i_C] - 1\n",
    "# )\n",
    "# temperature_gap = ...\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv (3.12.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.12.3"
  },
  "x3a_merged_sources": [
   "Session2/DICE_Notebook2-1_optimization_basics.ipynb"
  ],
  "x3a_provenance": {
   "allowed_changes": [
    "file rename",
    "session title",
    "local-path normalization",
    "optional installation cell",
    "stored-output normalization"
   ],
   "session": 5,
   "source_files": [
    "Session2/DICE_Notebook2-1_optimization_basics.ipynb"
   ]
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}