Time_Series_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "pacific-discount",
   "metadata": {},
   "source": [
    "# Time series analyzer\n",
    "## A Climate Asset Pricing model -- AR6 analysis notebook\n",
    "\n",
    "In this notebook, we provide a smooth framework for users to analyze the output of CAP6. This is meant to be generic analysis, but can be used to make nicer figures; we use it here to make the preferred and featured runs figure in the paper. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alien-pledge",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import datetime\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "\n",
    "from src.tree import TreeModel\n",
    "from src.climate import BPWClimate\n",
    "from src.emit_baseline import BPWEmissionBaseline\n",
    "from src.damage import BPWDamage\n",
    "from src.tools import import_csv\n",
    "from src.analysis.output_unpacker import OutputUnpacker\n",
    "from src.analysis.tree_diagrams import TreeDiagram\n",
    "\n",
    "# plotting parameters\n",
    "color_list = ['#000000', '#E69F00', '#56B4E9', '#009E73', '#F0E442', '#0072B2', '#D55E00', '#CC79A7'] * 2\n",
    "marker_list = ['o', 's', 'P', '+', 'D', 'v', '3'] * 2\n",
    "markersize = 6\n",
    "linestyle_list = ['solid', 'dashed', 'dashdot', 'dotted'] * 4\n",
    "linewidth = 2\n",
    "node_times = [2020, 2030, 2060, 2100, 2150, 2200]\n",
    "x_label = 'Year'\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina' \n",
    "cdds_params={'axes.linewidth': 3,\n",
    " 'axes.axisbelow': False,\n",
    " 'axes.edgecolor': 'black',\n",
    " 'axes.facecolor': 'None',\n",
    " 'axes.grid': False,\n",
    " 'axes.labelcolor': 'black',\n",
    " 'axes.spines.right': False,\n",
    " 'axes.spines.top': False,\n",
    " 'axes.titlesize': 20,\n",
    " 'axes.labelsize': 20,\n",
    " 'axes.titlelocation': 'left',\n",
    " 'figure.facecolor': 'white',\n",
    " 'figure.figsize': (18, 10),\n",
    " 'lines.solid_capstyle': 'round',\n",
    " 'lines.linewidth': 2.5,\n",
    " 'patch.edgecolor': 'w',\n",
    " 'patch.force_edgecolor': True,\n",
    " 'text.color': 'black',\n",
    " 'legend.frameon': False,\n",
    " 'xtick.bottom': True,\n",
    " 'xtick.major.width': 3,\n",
    " 'xtick.major.size': 6,\n",
    " 'xtick.color': 'black',\n",
    " 'xtick.direction': 'out',\n",
    " 'xtick.top': False,\n",
    " 'ytick.color': 'black',\n",
    " 'ytick.direction': 'out',\n",
    " 'ytick.left': True,\n",
    " 'ytick.right': False,\n",
    " 'ytick.color' : 'black',\n",
    " 'ytick.major.width': 3, \n",
    " 'ytick.major.size': 6,\n",
    " 'axes.prop_cycle': plt.cycler(linestyle=linestyle_list),\n",
    " 'font.size': 16,\n",
    " 'font.family': 'serif'}\n",
    "plt.rcParams.update(cdds_params)\n",
    "\n",
    "# make base filename\n",
    "today = datetime.datetime.now()\n",
    "year = str(today.year)\n",
    "day = str(today.day)\n",
    "month = str(today.month)\n",
    "\n",
    "basefile = ''.join([month, '-', day, '-', year, '-'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "legislative-fortune",
   "metadata": {},
   "source": [
    "## Make lists of necessary model classes\n",
    "\n",
    "Make sure this matches the main file! To do this properly, the right classes with the right model parameters need to be initiated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "lonely-holocaust",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Import header, indices, and data from reference .csv file.\n",
    "\"\"\"\n",
    "\n",
    "data_csv_file = 'research_runs'\n",
    "header, indices, data = import_csv(data_csv_file, delimiter=',', indices=2)\n",
    "\n",
    "\"\"\"Run TCREZClimate for runs of interest. To change which runs the code does,\n",
    "change desired_runs.\n",
    "\"\"\"\n",
    "\n",
    "# runs needing analysis\n",
    "done_runs = [0,1,2,3] # main runs\n",
    "\n",
    "tree_list = []\n",
    "damage_list = []\n",
    "climate_list = []\n",
    "emit_baseline_list = []\n",
    "\n",
    "for i in done_runs:\n",
    "    ra, eis, pref, growth, tech_chg, tech_scale, dam_func,\\\n",
    "            baseline_num, tip_on, bs_premium, d_unc, t_unc, nfl = data[i]\n",
    "\n",
    "    baseline_num = int(baseline_num)\n",
    "    dam_func = int(dam_func)\n",
    "    tip_on = int(tip_on)\n",
    "    d_unc = int(d_unc)\n",
    "    t_unc = int(t_unc)\n",
    "    nfl = int(nfl)\n",
    "\n",
    "    \"\"\"Initialize model classes. First is the tree model.\n",
    "    \"\"\"\n",
    "\n",
    "    tmp_tree = TreeModel(decision_times=[0, 10, 40, 80, 130, 180, 230],\n",
    "                  prob_scale=1.0)\n",
    "\n",
    "    \"\"\"Emission baseline model. We also run its setup function.\n",
    "    \"\"\"\n",
    "\n",
    "    tmp_baseline_emission_model = BPWEmissionBaseline(tree=tmp_tree,\n",
    "                                                  baseline_num=baseline_num)\n",
    "    tmp_baseline_emission_model.baseline_emission_setup()\n",
    "\n",
    "    \"\"\"Climate class. We set draws to the number of Monte Carlo samples to take\n",
    "    from damage distributions, such as TCRE, and if the floor is on.\n",
    "    \"\"\"\n",
    "\n",
    "    draws = 3 * 10**6\n",
    "    tmp_climate = BPWClimate(tmp_tree, tmp_baseline_emission_model, draws=draws)\n",
    "\n",
    "    \"\"\"Damage class. We set draws to the number of Monte Carlo samples to take\n",
    "    from damage distributions, such as TCRE. We also pass a list of constant\n",
    "    values of mitigation for the damage simulation.\n",
    "    \"\"\"\n",
    "\n",
    "    d_m = 0.01\n",
    "    mitigation_constants = np.arange(0, 1 + d_m, d_m)[::-1]\n",
    "    tmp_df = BPWDamage(tree=tmp_tree, emit_baseline=tmp_baseline_emission_model,\n",
    "                   climate=tmp_climate, mitigation_constants=mitigation_constants,\n",
    "                   draws=draws)\n",
    "\n",
    "    \"\"\"Run the damage simulation. If a damage simulation has already been run,\n",
    "    then comment the line below and uncomment the df.import_damages() line.\n",
    "\n",
    "    NOTE: import_damages has the optional argument for a filename. This is\n",
    "    important if you're varying the damage function across simulations.\n",
    "\n",
    "    The default is currently \"simulated_damages_TCRE.csv\".\n",
    "    \"\"\"\n",
    "\n",
    "    damsim_filename = ''.join([\"simulated_damages_df\", str(dam_func),\n",
    "                               \"_TP\", str(tip_on), \"_SSP\", str(baseline_num),\n",
    "                               \"_dunc\", str(d_unc), \"_tunc\", str(t_unc)])\n",
    "\n",
    "    tmp_df.import_damages(file_name=damsim_filename)\n",
    "    \n",
    "    tree_list.append(tmp_tree)\n",
    "    emit_baseline_list.append(tmp_baseline_emission_model)\n",
    "    climate_list.append(tmp_climate)\n",
    "    damage_list.append(tmp_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hazardous-circumstances",
   "metadata": {},
   "source": [
    "## Make list of files of interest\n",
    "\n",
    "Here we make lists of file prefixes (i.e., whatever the `name` parameter was when the simulation was executed) and paths. We take the list of paths and prefixes and make lists of filenames, which will be passed to `OutputUnpacker` later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "yellow-homeless",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of prefixes to output files of EZClimate. i.e., the string before \"_node_period_output.csv\".\n",
    "# in EZClimate main file, this is the parameter \"name\" located at about line 42 \n",
    "\n",
    "# prefixes (\"name\" variable in main file)\n",
    "prefix_list = ['pref-run', 'feat-15', 'feat-25', 'feat-3']\n",
    "\n",
    "# path to data, in order of the prefixes listed above\n",
    "path_list = [os.getcwd() + '/data/'] * len(prefix_list) \n",
    "\n",
    "# file description - this will be used in plot titles and helps keep track of which file is which\n",
    "# in the prefix_list\n",
    "descriptions = ['2% discount rate', '1.5% discount rate', '2.5% discount rate', '3% discount rate']\n",
    "\n",
    "# make list of file names to be unpacked\n",
    "output_list = []\n",
    "picklefile_list = []\n",
    "\n",
    "# total files\n",
    "N_files = len(prefix_list)\n",
    "\n",
    "# populate output_list and picklefile_list\n",
    "for file in range(0, N_files):\n",
    "    tmp_filename = path_list[file] + prefix_list[file] + \"_node_period_output.csv\"\n",
    "    tmp_picklefile_name = path_list[file] + prefix_list[file] + \"_log.pickle\"\n",
    "    output_list.append(tmp_filename)\n",
    "    picklefile_list.append(tmp_picklefile_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subjective-greene",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use OutputUnpacker to \"unpack\" the files above\n",
    "\n",
    "# make one array of objects which correspond to the 'output' files (i.e., the files that end with\n",
    "# _node_period_output.csv) and another for the pickle files\n",
    "output_obj_list = []\n",
    "pickle_obj_list = []\n",
    "\n",
    "# populate these arrays\n",
    "for file in range(0, N_files):\n",
    "    tmp_output_obj = OutputUnpacker(output_list[file], descriptions[file], 'output', \n",
    "                                    tree=tree_list[file], emit_baseline=emit_baseline_list[file],\n",
    "                                    climate=climate_list[file], damage=damage_list[file])\n",
    "    output_obj_list.append(tmp_output_obj)\n",
    "    tmp_pickle_obj = OutputUnpacker(picklefile_list[file], descriptions[file], 'pickle',\n",
    "                                    tree=tree_list[file], emit_baseline=emit_baseline_list[file],\n",
    "                                    climate=climate_list[file], damage=damage_list[file])\n",
    "    pickle_obj_list.append(tmp_pickle_obj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cheap-concert",
   "metadata": {},
   "source": [
    "## SAVE FIGS?!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adjustable-gentleman",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_fig = False\n",
    "fig_filename_prefix = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aquatic-apache",
   "metadata": {},
   "source": [
    "## Time series plots\n",
    "\n",
    "Here we plot a number of time series of key model variables, such as:\n",
    "- average mitigation\n",
    "- expected mitigation\n",
    "- average temperature\n",
    "- average emissions\n",
    "- cost of carbon\n",
    "- CO$_2$ concentrations\n",
    "- economic damages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pointed-culture",
   "metadata": {},
   "source": [
    "### Cost of carbon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "driving-velvet",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(9.5,7))\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(node_times, np.mean(output_obj_list[file].price_path, axis=0), \n",
    "                  color=color_list[file], marker=None, \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "    print(np.mean(output_obj_list[file].price_path, axis=0)[0])\n",
    "\n",
    "\n",
    "ax.set_ylabel(\"Carbon price ($ USD 2020)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.set_ylim((0,250))\n",
    "ax.set_yticks([0,50, 100,150, 200, 250])\n",
    "ax.set_xticks([2020, 2100, 2200])\n",
    "sns.despine(trim=True, offset=10)\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, fig_filename_prefix, '-cost.png']), dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "invalid-crack",
   "metadata": {},
   "source": [
    "### Average mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "regular-peninsula",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(node_times, np.mean(pickle_obj_list[file].m_opt_path, axis=0), \n",
    "            color=color_list[file], marker=marker_list[file], \n",
    "            label=pickle_obj_list[file].description, markersize=markersize,\n",
    "            linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "ax.set_ylabel(\"Average mitigation\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "ax.tick_params(axis='both')\n",
    "#ax.axhline(1, color='g', linestyle='dashed')\n",
    "ax.set_ylim((0,1.6))\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, fig_filename_prefix, '-mit.png']), dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "persistent-audit",
   "metadata": {},
   "source": [
    "### Expected mitigation\n",
    "\n",
    "Difference from average, as this is the weighted average of mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressive-baking",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(node_times, np.mean(output_obj_list[file].exp_mit_path, axis=0), \n",
    "            color=color_list[file], marker=marker_list[file], \n",
    "            label=pickle_obj_list[file].description, markersize=markersize,\n",
    "            linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "ax.set_ylabel(\"Expected mitigation\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "#ax.axhline(1, color='g', linestyle='dashed')\n",
    "ax.set_ylim((0,1.3))\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, fig_filename_prefix, '-exp-mit.png']), dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interested-convert",
   "metadata": {},
   "source": [
    "### (Mitigated) emissions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "latest-species",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(tmp_baseline_emission_model.times[:-5], tmp_baseline_emission_model.baseline_gtco2[:-5], \n",
    "        color='slategrey', label=\"SSP2\")\n",
    "ax.plot(tmp_baseline_emission_model.times[:-5], (1 - 0.6) * tmp_baseline_emission_model.baseline_gtco2[:-5], \n",
    "        color='slategrey', label=\"60% Mitigation\", linestyle=\"dashed\")\n",
    "ax.plot(tmp_baseline_emission_model.times[:-5], (1 - 0.9) * tmp_baseline_emission_model.baseline_gtco2[:-5], \n",
    "        color='slategrey', label=\"90% Mitigation\", linestyle=\"dashdot\")\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].emis_gtco2_full, axis=0), \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "    \n",
    "ax.set_ylabel(\"Average emissions (GtCO$_2$ yr$^{-1}$)\")\n",
    "ax.set_xlabel(x_label)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polyphonic-manufacturer",
   "metadata": {},
   "source": [
    "### Average temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cooperative-excitement",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].temp_full, axis=0), \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "    \n",
    "ax.set_ylabel(\"Average temperature anomaly (K)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "#ax.set_ylim((0,5))\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, fig_filename_prefix, '-temp.png']), dpi=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "overhead-convergence",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].temp_full, axis=0) - np.mean(pickle_obj_list[file].temp_full, axis=0)[0], \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "    \n",
    "ax.set_ylabel(\"Average temperature anomaly from 2020 (K)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sacred-chaos",
   "metadata": {},
   "source": [
    "### CO$_2$ concentrations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "competent-investigation",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(node_times, np.mean(output_obj_list[file].ghg_lvl_path, axis=0), \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "\n",
    "ax.set_ylabel(\"Average CO$_2$ concentrations (ppm)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, fig_filename_prefix, '-co2.png']), dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "purple-width",
   "metadata": {},
   "source": [
    "### Economic damages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interior-health",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(node_times, np.mean(pickle_obj_list[file].econ_dam_path, axis=0), \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "ax.set_ylabel(\"Average economic damages (% GDP)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(x_label)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "champion-fellowship",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for file in range(0, N_files):\n",
    "    ax.plot(np.mean(pickle_obj_list[file].temp_path, axis=0), np.mean(pickle_obj_list[file].econ_dam_path, axis=0), \n",
    "                  color=color_list[file], marker=marker_list[file], \n",
    "                  label=pickle_obj_list[file].description, markersize=markersize,\n",
    "                  linestyle=linestyle_list[file], linewidth=linewidth)\n",
    "\n",
    "ax.set_ylabel(\"Average economic damages (% GDP)\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"Temperature (K)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "blind-sphere",
   "metadata": {},
   "source": [
    "## 6 Panels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interior-antique",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.transforms as mtransforms\n",
    "\n",
    "fig, ax = plt.subplot_mosaic([['A', 'B', 'C'], ['D', 'E', 'F']], sharex=True,\n",
    "                             gridspec_kw={'height_ratios': [1, 1], 'width_ratios': [1, 1, 1]},\n",
    "                            figsize=(22,10))\n",
    "\n",
    "pref_lw = 3.75\n",
    "feat_lw = 2\n",
    "pref, = ax['A'].plot(node_times, np.mean(output_obj_list[0].price_path, axis=0), \n",
    "              color=color_list[0], \n",
    "              label=pickle_obj_list[0].description,\n",
    "              linestyle=linestyle_list[0], linewidth=pref_lw)\n",
    "feat15, = ax['A'].plot(node_times, np.mean(output_obj_list[1].price_path, axis=0), \n",
    "              color=color_list[1], \n",
    "              label=pickle_obj_list[1].description,\n",
    "              linestyle=linestyle_list[1], linewidth=linewidth)\n",
    "feat25, = ax['A'].plot(node_times, np.mean(output_obj_list[2].price_path, axis=0), \n",
    "              color=color_list[2], \n",
    "              label=pickle_obj_list[2].description,\n",
    "              linestyle=linestyle_list[2], linewidth=linewidth)\n",
    "feat3, = ax['A'].plot(node_times, np.mean(output_obj_list[3].price_path, axis=0), \n",
    "              color=color_list[3], \n",
    "              label=pickle_obj_list[3].description,\n",
    "              linestyle=linestyle_list[3], linewidth=linewidth)\n",
    "\n",
    "#ax['B'].hlines(1.5, 2020, 2200, alpha=0.5, color='grey', linestyle='dashed', linewidth=1.75)\n",
    "#ax['B'].vlines(2100, 1, 2, alpha=0.5, color='grey', linestyle='dashed', linewidth=1.75)\n",
    "\n",
    "# other plots\n",
    "for file in range(0, N_files):\n",
    "    if file == 0:\n",
    "        ax['B'].plot(node_times, np.mean(output_obj_list[file].exp_mit_path, axis=0) * 100, \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=pref_lw)\n",
    "        \n",
    "        ax['C'].plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].emis_gtco2_full, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      linestyle=linestyle_list[file], linewidth=pref_lw)\n",
    "        \n",
    "        ax['D'].plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].temp_full, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=pref_lw)\n",
    "\n",
    "        ax['E'].plot(node_times, np.mean(output_obj_list[file].ghg_lvl_path, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=pref_lw)\n",
    "\n",
    "        ax['F'].plot(node_times, np.mean(pickle_obj_list[file].econ_dam_path * 100, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=pref_lw)\n",
    "        \n",
    "    else:\n",
    "        ax['B'].plot(node_times, np.mean(output_obj_list[file].exp_mit_path, axis=0) * 100, \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=feat_lw)\n",
    "        \n",
    "        ax['C'].plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].emis_gtco2_full, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      linestyle=linestyle_list[file], linewidth=feat_lw)\n",
    "        \n",
    "        ax['D'].plot(tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[file].temp_full, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=feat_lw)\n",
    "\n",
    "        ax['E'].plot(node_times, np.mean(output_obj_list[file].ghg_lvl_path, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=feat_lw)\n",
    "\n",
    "        ax['F'].plot(node_times, np.mean(pickle_obj_list[file].econ_dam_path * 100, axis=0), \n",
    "                      color=color_list[file], \n",
    "                      label=pickle_obj_list[file].description,\n",
    "                      linestyle=linestyle_list[file], linewidth=feat_lw)\n",
    "        \n",
    "        \n",
    "# plot SSP2 for clarity\n",
    "ax['C'].plot(tmp_baseline_emission_model.times[:-5], tmp_baseline_emission_model.baseline_gtco2[:-5], \n",
    "        color='slategrey', label=\"SSP2\", linewidth=feat_lw)\n",
    "ax['C'].legend(loc='right', fontsize=19)\n",
    "\n",
    "panels = ['A', 'B', 'C', 'D', 'E', 'F']\n",
    "for panel in panels:\n",
    "    ax[panel].set_xticks([2020, 2100, 2200])\n",
    "\n",
    "# set ticks\n",
    "ax['A'].set_yticks([0, 200, 400, 600])\n",
    "ax['B'].set_yticks([0, 20, 40, 60, 80, 100])\n",
    "ax['C'].set_yticks([0, 10, 20, 30, 40, 50])\n",
    "ax['D'].set_yticks([1, 1.5, 2, 2.5])\n",
    "ax['E'].set_yticks([350, 400, 450, 500, 550])\n",
    "ax['F'].set_yticks([0, 2, 4, 6])\n",
    "\n",
    "# set labels\n",
    "ax['A'].set_ylabel(\"Cost (2020 $USD)\", fontsize=19)\n",
    "ax['B'].set_ylabel(\"Expected abatement (%)\", fontsize=19)\n",
    "ax['C'].set_ylabel(\"Emissions (GtCO$_2$ yr$^{-1}$)\", fontsize=19)\n",
    "ax['D'].set_ylabel(\"Temperature ($^{\\circ}$C)\", fontsize=19)\n",
    "ax['E'].set_ylabel(\"CO$_2$ concentration (ppm)\", fontsize=19)\n",
    "ax['F'].set_ylabel(\"Econ. Damage (% GDP)\", fontsize=19)\n",
    "\n",
    "ax['D'].set_xlabel(\"Year\", fontsize=19)\n",
    "ax['E'].set_xlabel(\"Year\", fontsize=19)\n",
    "ax['F'].set_xlabel(\"Year\", fontsize=19)\n",
    "#ax['A'].set_ylim((0,310))\n",
    "\n",
    "# bottom legend\n",
    "fig.legend([feat15, pref, feat25, feat3], [\"1.5% discount rate\", \"2% discount rate\", \"2.5% discount rate\", \"3% discount rate\"], \n",
    "           bbox_to_anchor=(0.48, -0.045), \n",
    "           loc='lower center', ncol=4, fancybox=True, shadow=True, fontsize=18)\n",
    "\n",
    "# add labels\n",
    "right = ['A', 'B', 'C', 'D', 'E', 'F']\n",
    "for label in right:\n",
    "    # label physical distance in and down:\n",
    "    trans = mtransforms.ScaledTranslation(0, 0.0, fig.dpi_scale_trans)\n",
    "    ax[label].text(0.9, 1.0, label, transform=ax[label].transAxes + trans, fontsize=22, fontweight='bold',\n",
    "            verticalalignment='top', bbox=dict(facecolor='none', edgecolor='none', pad=1))\n",
    "    \n",
    "sns.despine(trim=True, offset=10)\n",
    "\n",
    "if save_fig:\n",
    "    fig.savefig(''.join([basefile, 'pref-runs-NFL.png']),bbox_inches='tight', dpi=400)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0c21eef7",
   "metadata": {},
   "source": [
    "### Print some relevant values (temperature and so on)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6f261e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(pickle_obj_list)):\n",
    "    print(np.mean(pickle_obj_list[i].econ_dam_path, axis=0) * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb7dc1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(pickle_obj_list)):\n",
    "    print(np.mean(pickle_obj_list[i].temp_full, axis=0)[0], \n",
    "         np.mean(pickle_obj_list[i].temp_full, axis=0)[1],\n",
    "         np.mean(pickle_obj_list[i].temp_full, axis=0)[4],\n",
    "         np.mean(pickle_obj_list[i].temp_full, axis=0)[8],\n",
    "         np.mean(pickle_obj_list[i].temp_full, axis=0)[13],\n",
    "         np.mean(pickle_obj_list[i].temp_full, axis=0)[18])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db026a21",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_baseline_emission_model.times[:-5], np.mean(pickle_obj_list[0].emis_gtco2_full, axis=0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}