standard error of the trend question

notebooks/analysis/error-sources/trend-std-err.ipynb

@@ -0,0 +1,141 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "ce60cc80",
+   "metadata": {},
+   "source": [
+    "# What is the confidence interval of the trend of the broken linear model?\n",
+    "This notebook answers the question on what the trend is of the sea-level rise trend of the broken linear model. This question is asked in the context of a question from [parliament](https://www.tweedekamer.nl/kamerstukken/kamervragen/detail?id=2023Z07751&did=2023D21773).\n",
+    "\n",
+    "The trend estimate is reported in the main sealevel document. There the following numbers are reported:\n",
+    "\n",
+    "| Parameter     | estimate | stderr| ci l.   | ci u.   |\n",
+    "|---------------|----------|-------|---------|---------|\n",
+    "| Constant      | -39.2860 | 3.427 | -46.004 | -32.568 |\n",
+    "| Trend         |   1.7827 | 0.085 | 1.615   | 1.950   |\n",
+    "| +trend (1993) |\t1.0779 | 0.342 | 0.407   | 1.748   |\n",
+    "\n",
+    "We can use [error propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty) to combine both parameters. We'll use the [uncertainties](https://pythonhosted.org/uncertainties/) package for this."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d378ef55",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import uncertainties"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "800af798",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(1.7827+/-0.085, 1.0779+/-0.342)"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# standard error of the trend\n",
+    "trend_se = uncertainties.ufloat(1.7827, 0.085)\n",
+    "\n",
+    "# standard error of the added trend after 1993\n",
+    "extra_trend_se = uncertainties.ufloat(1.0779, 0.342)\n",
+    "trend_se, extra_trend_se"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "df7b1751",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2.8606+/-0.3524045970188244"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# combine the 2 trends\n",
+    "combined_se = trend_se + extra_trend_se\n",
+    "combined_se"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "5d15a56e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute the confidence interval assuming normal distribution\n",
+    "ci_lower = combined_se.nominal_value - combined_se.std_dev * 1.96\n",
+    "ci_upper = combined_se.nominal_value + combined_se.std_dev * 1.96"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "b5f6311f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The broken linear is 2.9mm/yr with a standard error of 0.4mm/yr. This implies a 95% confidence interval of (2.2, 3.6) in mm/yr.\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(f\"The broken linear is {combined_se.nominal_value:.1f}mm/yr with a standard error of {combined_se.std_dev:.1f}mm/yr. This implies a 95% confidence interval of ({ci_lower:.1f}, {ci_upper:.1f}) in mm/yr.\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5a6b14c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.16"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}