Risk 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3e879630",
   "metadata": {},
   "source": [
    "# Tropical Cyclones\n",
    "## Introduction\n",
    " \n",
    "In this session we will consider the threat posed by *tropical cyclones* -- the biggest beast of them all when it comes to weather hazards. We will be assessing the risk they currently pose (or, rather, have historically posed). \n",
    "\n",
    "The use of the term *risk* there was not accidental. We will be going beyond simpy assessing the hazard; we will be estimating impact as measured through total finanical losses (normalised to 2018 USD [$] to account for inflation). \n",
    "\n",
    "To work through these important concepts, we will be following a hypothetical scenario, outlined below, in which the US Treasury have hired you and your team as experts in tropical cyclone risk. \n",
    "\n",
    "Note that, as last week, some questions are included in the text; the answers to these are at the end of the notebook (and may also be disclosed in the text). \n",
    "\n",
    "### Scenario\n",
    "\n",
    "Your task is to summarise the potential financial impact of tropical cyclones over the next 100 years. They have asked you to assess: \n",
    "\n",
    "[1] The expected annual cost of tropical cyclones to the USA, assuming no further climate change. \n",
    "\n",
    "[2] The cost of a very expensive (i.e., rare) year -- again in a stationary climate. The Treasury define this as a year with costs exceeded in no more than 1 % of years: the one-in-100-year cost. \n",
    "\n",
    "### Data\n",
    "\n",
    "You will be working estimates of tropical cyclone damage from [ICAT](http://www.icatdamageestimator.com/viewdata). All costs are provided in US Dollars, adjusted for inflation to 2018 equaivalent.   \n",
    "\n",
    "### Analysis\n",
    "\n",
    "Follow the workbook to address the requests of the Treasury. Remember to run code cells when you get to them by either pressing \"run\" at the of the browser window, or by hitting \"ctrl\" + \"enter\" on your keyboard. \n",
    "\n",
    "The first code cell (below) sets the notebook up for analysis. It also summarises the economic impacts from past hurricanes, printing the mean losses per year and per hurricane. The top ten years (by cost) are printed, and the most costly year is identified. Note that the values here are simple descriptive statistics. \n",
    "\n",
    "Run this first code cell now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "144f00fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing and preliminaries\n",
    "import pandas as pd, numpy as np, matplotlib.pyplot as plt\n",
    "from random import randint\n",
    "from scipy import stats\n",
    "from IPython.display import YouTubeVideo\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "_nsim=10000\n",
    "data=pd.read_csv(\"Data/HurricaneDamage.csv\",index_col=\"lf_ISO_TIME\",parse_dates=True)\n",
    "data=data.loc[data[\"ND\"]>0] # throw away TCs that caused no damage\n",
    "data[\"ND\"]/=1e9 # Transform cost to billions\n",
    "anncost=data[\"ND\"].resample(\"Y\").sum()\n",
    "annidx=pd.date_range(start=\"%.0f/01/01\"%anncost.index.year[0],end=\"%.0f/01/01\"%anncost.index.year[-1],freq=\"1y\")\n",
    "anncost=anncost.reindex(annidx)\n",
    "ntc=len(data)\n",
    "ny=len(anncost)\n",
    "maxcost_storm=data.loc[data[\"ND\"]==data[\"ND\"].max()]\n",
    "maxcost_year=anncost.loc[anncost==anncost.max()]\n",
    "rate=ntc/ny # this is the rate for our TC occurrence model \n",
    "\n",
    "print(\"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\")\n",
    "print(\"%.0f damaging hurricanes in %.0f years\"%(ntc,ny))\n",
    "print(\"Rate = %.1f storms/year\"%rate)\n",
    "print(\"*** Cost summaries (amounts in billion, 2018 $ equivalent)***\")\n",
    "print(\"Mean annual cost is $%.1f\"%anncost.mean())\n",
    "print(\"Top ten years: \")\n",
    "for i in range(10):\n",
    "    print(\"\\t[%.0f] %.0f ($%.2f)\"%(i+1,(anncost.sort_values()[::-1]).index.year[i],\n",
    "          (anncost.sort_values()[::-1])[i]))\n",
    "print(\"Mean cost/storm = $%.2f\"%data[\"ND\"].mean())\n",
    "print(\"Most costly storm = %.0f ($%.2f)\"%(maxcost_storm.index.year[0],maxcost_storm[\"ND\"][0]))\n",
    "print(\"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9adeecff",
   "metadata": {},
   "source": [
    "Take a moment to appreciate the scale of the costs we're talking about here. For context, the cost of the Furlough scheme in the UK during the Covid-19 pandemic cost ~£70 billion (i.e., just under $100 billion). \n",
    "\n",
    "You may recall the devastating 2005 Hurricane season (in which Hurricane Katrina devestated New Orleans); and maybe 2012 (when Hurricane Sandy struck New York). But perhaps 1926 -- and the Great Miami Hurricane -- is an event you have not yet heard of. Run the code below to find out more. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b160b6c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A video explaining the most impactful hurricane you've never heard of. \n",
    "YouTubeVideo('LRAJ7Bc0O5E',width=800, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2879440a",
   "metadata": {},
   "source": [
    "## The expected value\n",
    "\n",
    "- [**Q1**] One way of estimating the expected cost is to compute the *sample* mean -- the 'average' of the annual costs. Using that approach, we already have an estimate of the expectation (in 2018 $) from the summary printed above. What is it?\n",
    "\n",
    "That was easy! But is it a useful answer? Possibly *not*, because as shown by the histogram below, *very* large values (i.e., expensive years) are *not* that unlikley, and that can make the sample mean -- of a relatively short dataset -- unrepresentative of the *expected* annual cost (i.e., the mean value we would observe if we could observe all *possible* years). \n",
    "\n",
    "Run the code to see the distribution of annual costs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9cc0ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate histogram of the annual costs\n",
    "fig,ax=plt.subplots(1,1)\n",
    "h=ax.hist(anncost,bins=15,facecolor='green',edgecolor='k',linewidth=1,alpha=0.25)\n",
    "ax.set_ylabel(\"Count\")\n",
    "ax.set_xlabel(\"Cost ($2018 Billion)\")\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "655b7f1b",
   "metadata": {},
   "source": [
    "To put the magnitude of variability into perspective, the maximum annual bill (\\\\$265 billion) is not too far short of 100 times larger than the median value ($3 billion). If such variability existed in human heights, we should expect to encounter people well over 200 m tall! The possibility of such large \"events\" can make the mean *very* sensitive to sampling variability. To appreciate this, imagine calculating the mean height of all people in this room -- hoping it's a representative sample of all students on campus -- and then picture how that would change if a ~200 m tall student accidentally (tried) to pop in; your mean would go up -- a lot! \n",
    "\n",
    "We would be on a much firmer footing to compute the expected value if we knew the statistical distribution of the *population* (of heights, in this example). Knowing the population distribution is equivalent to being able to measure the heights of *all* students on campus. However, seldom can we measure all \"members\" of a population, so we instead make an assumption about the *probability distribution* (e.g., the Normal, GEV, etc.,) that represents the population, and use the properties of *that* to compute things like the expected value and, for eample, the probability of values beyond given thresholds (via the *CDF* -- as discussed last session).\n",
    "\n",
    "It would be very helpful, then, to compute the expected value from the probability distribution that represents *annual costs*. We'll get into that soon; first, let's deal with a slightly simpler question: what does the distribution of *costs per cyclone* look like?  \n",
    "\n",
    "\n",
    "## The distribution of cyclone impacts\n",
    "\n",
    "The reason why we're breaking things down like this is because, when you think on the processes for a little bit, it makes sense to deal with the damage per cyclone *seperately* from the number of cylones per year. The conceptual basis is explained in [Katz (2002)](https://journals.ametsoc.org/view/journals/apme/41/7/1520-0450_2002_041_0754_smohd_2.0.co_2.xml), but one intuitive reason to seperate them is because the annual cost can equal **exactly** zero (a year in which no damaging hurricanes make landfall), *or **any** other amount*, so our final statistical model must be capable of predicting exactly zero damages, *and* any other amount -- both with the right frequency. For this reason (and some others) it turns out separating our statistical approach into models (distributions) that deal with cyclone occurrence (number of storms per year), from those which consider impact (the cost per storm), is a good idea. So, with that explainer out the way, let's get back to assessing the distribution of costs per hurricane. Run the code to take a look. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1e8e9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking the fit for a log-normal distribution \n",
    "fig,ax=plt.subplots(1,2)\n",
    "fig.set_size_inches(8,3)\n",
    "h1=ax.flat[0].hist(anncost,bins=15,facecolor='green',edgecolor='k',linewidth=1.,density=True,alpha=0.25)\n",
    "# Taking logs of the costs/storm, and then computing the mean and stdv of the logs\n",
    "lncost=np.log(data[\"ND\"]) # log transform\n",
    "lnmean=np.mean(lncost) # mean of log-transformed data\n",
    "lnstd=np.std(lncost) # stdv of log-transformed data\n",
    "x=np.linspace(np.min(lncost)-1,np.max(lncost)+2,100)\n",
    "ax.flat[1].plot(x,stats.norm.pdf(x, loc=lnmean,scale=lnstd),color='red')\n",
    "h2=ax.flat[1].hist(lncost,bins=15,density=True,facecolor='green',edgecolor='k',linewidth=1.,alpha=0.25)\n",
    "ax.flat[0].grid(); ax.flat[1].grid()\n",
    "ax.flat[0].set_ylabel(\"Density\")\n",
    "ax.flat[0].set_xlabel(\"Cost ($2018 Billion)\")\n",
    "ax.flat[1].set_xlabel(\"ln{Cost ($2018 Billion)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e86406b",
   "metadata": {},
   "source": [
    "You will see in the plots that the distribution of costs/hurricane (left) looks very similar to the cost/year distribution: it has very strong positive skew, meaning values *much* greater than the median are not that unlikely. \n",
    "\n",
    "The ditribution on the right, however, looks very familiar. It seems to be very well described by the Normal distribution, as shown by the way the red line (the Normal pdf) matches the histogram well. The quantity plotted in this right-hand panel is the cost/hurricane series *after it has been log transformed*. \n",
    "\n",
    "*[Note that a log transform just means \"taking the logarithm\" of all values in the dataset; and the logarithm of a number is the exponent that a base must be raised to in order to equal that number. For example, if the base is 10 and the number is 100, the logarithm would be 2 (as 10$^{2}$=100). In our log transform we use the base e, which has a value of ~2.71]*\n",
    "\n",
    "The log transform has the effect of reducing larger numbers by greater amounts, and it results in our transformed data following a Normal distrubtion. Technically, this means that the underlying data follow what's called a *Log-Normal distribution*. \n",
    "\n",
    "The formula for the expected value (E[X]) from the Log-Normal distribution (given on pg. 88 of the [Wilks textbook](https://ebookcentral.proquest.com/lib/kcl/reader.action?docID=689817)) is: \n",
    "\n",
    "\\begin{aligned}\n",
    "E[X]=e^{\\mu + \\frac{\\sigma^{2}}{2}}\n",
    "\\end{aligned}\n",
    "\n",
    "where $\\sigma$ and $\\mu$ are the mean and standard deviation, respectively, of the *log-transformed* variable.\n",
    "\n",
    "The equation therefore tells us that to compute the expected value (i.e., the mean) for the cost/hurricane, we simply (a) log transform the data; (b) compute the mean and standard deviation of the transformed data; and (c) plug those values into the equation above. The code below does this. Run it to find out the expected cost per hurricane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8383f63",
   "metadata": {},
   "outputs": [],
   "source": [
    "lncost=np.log(data[\"ND\"]) # Log transform\n",
    "lnmean=np.mean(lncost) # Computing the mean of the log-transformed data\n",
    "lnstd=np.std(lncost) # Computing the standard deviation of the log-transformed data\n",
    "E=np.exp(lnmean+lnstd**2/2.)\n",
    "print(\"The expected value is $%.2f billion per hurricane\"%E)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "986a895a",
   "metadata": {},
   "source": [
    "- [**Q2**] How does this compare with the sample average computed at the beginning of this notebook?\n",
    "\n",
    "## The number of events per year\n",
    "\n",
    "We now have a statistical model the cost per hurricane, but in order to estimate the expected cost per year, and the cost of a very expensive year (i.e., the one-in-200-year event), we need to also consider a statistical model for the *number of hurricanes per year*. \n",
    "\n",
    "What does this distribution look like? Run the code below to take a look. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "993abe28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking the Poisson distribution\n",
    "# Compute the frequency with which ncylones/year are observed\n",
    "x=np.arange(0,10)\n",
    "of=np.array([np.sum(data.index.year==ii) for ii in anncost.index.year])\n",
    "off=np.array([np.sum(of==ii) for ii in x])\n",
    "gen=stats.poisson.rvs(rate,size=_nsim)\n",
    "sff=np.array([np.sum(gen==ii) for ii in x])/(_nsim/ny)\n",
    "\n",
    "fig,ax=plt.subplots(1,1)\n",
    "ax.bar(x,off,color='blue',alpha=0.8,label=\"Observed\")\n",
    "ax.bar(x,sff,color='red',alpha=0.25,label=\"Poisson\")\n",
    "ax.set_xticks(x)\n",
    "ax.set_xlabel(\"Hurricanes/year\")\n",
    "ax.set_ylabel(\"Frequency\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efd7527a",
   "metadata": {},
   "source": [
    "You will again observe a distribution with positive skew, albeit much less pronounced than in the costs dataset. Although it's not a perfect match, you should additionally notice that the *Poisson* distribution is a reasonably good fit to our data. It succeeds in reproducing the basic shape of our histogram. \n",
    "\n",
    "As explained by [Katz (2002)](https://journals.ametsoc.org/view/journals/apme/41/7/1520-0450_2002_041_0754_smohd_2.0.co_2.xml), the Poisson distribution is also theoretically suited to describing the annual frequency of damaging hurricanes, so we will persist with its use here. \n",
    "\n",
    "Conveniently, the expected value (E[N]) for the Poisson distribution (i.e., the average number of damaging hurricanes per year) is simply the mean *observed* frequency, known as the \"rate\". We'll denote this quantity $\\lambda$. \n",
    "\n",
    "- [**Q3**] What is the value for $\\lambda$ here? (Hint look at the stats printed out at the top of this notebook)**.\n",
    "\n",
    "Once we have determined $\\lambda$, we can compute the *expected cost per year* (E[Cy]) by combining the expected value for cost/hurricane E[Ch]) with the expected value of the number of hurricanes per year (E[N]):\n",
    "\n",
    "\\begin{aligned}\n",
    "E[Cy]=\\lambda e^{\\mu + \\frac{\\sigma^{2}}{2}}\n",
    "\\end{aligned}\n",
    "\n",
    "In other words, we just multiply the two expected values together! \n",
    "\n",
    "Run the code below to evaluate the expected cost below and print the number to screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca361dac",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ecy=np.exp(np.mean(lncost)+lnstd**2/2)*rate\n",
    "print(\"The expected cost is $%.2f billion/year\"%Ecy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e58d829a",
   "metadata": {},
   "source": [
    "- [**Q4**] How does this compare to the sample average computed at the beginning of this notebook?\n",
    "\n",
    "## The magnitude of the one-in-100-year event\n",
    "\n",
    "The answer above addresses the first request of the Treasury. We must now turn our attention to their second question. How severe should we expect the costs of hurricane damage to be in a very *bad* year? \n",
    "\n",
    "This is harder to address with simple formuale, so we will approach the problem via a statistical *simulation*, where we generate a very long sequence of \"possible years\" that are statistically consistent with our observations. We then summarise the plausible scenario via simple descriptive statistics. \n",
    "\n",
    "The algorithm to implement this should be quite intuitive. For each year, we will:\n",
    "\n",
    "[1] Select a random number from the Poisson distribution to represent the number of damaging hurricanes in that year \n",
    "\n",
    "[2] Draw a random number from the Log-Normal distribution to represent its damage\n",
    "\n",
    "[3] Add up the damage from all hurricanes in the year to obtain the total damage. \n",
    "\n",
    "\n",
    "If we do this for many, many years (say 100,000), we should get a very good idea of what's likely for the 1-in-100 year event (because we will experience 1000 events as 'extreme' in the 100,000 year series). Computing the 99th percentile from the simulated timeseries provides us with the estimate of the 1-in-100 year event. \n",
    "\n",
    "When you have read (and understood) the above, run the code below to generate some plausible scenarios! \n",
    "\n",
    "Once the simulations are complete a timeseries plot showing all the plausible years is plotted, and the 99th percentile is printed to screen. For reference, we also highlight (with a red line) the most costly year in the observations (1926). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c018a08f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stochastic simulation\n",
    "nsim=100000# This sets the number of years in our simulation\n",
    "nstorms=np.zeros(nsim,dtype=np.int)\n",
    "cost=np.zeros(nsim)\n",
    "for i in range(nsim):\n",
    "    nstorms[i]=int(stats.poisson.rvs(rate))\n",
    "    cost[i]=np.sum(np.exp(stats.norm.rvs(loc=lnmean,scale=lnstd,\n",
    "                                         size=nstorms[i])))\n",
    "fig,ax=plt.subplots(1,1)\n",
    "fig.set_size_inches(7,3)\n",
    "fig.set_dpi(200)\n",
    "ax.plot(cost)\n",
    "ax.set_xlabel(\"Simulated year number\")\n",
    "ax.set_ylabel(\"Cost ($ Billion)\")\n",
    "ax.axhline(anncost.max(),color='red',label=\"1926 cost\")\n",
    "ax.set_ylim(0,np.max(cost))\n",
    "ax.legend()\n",
    "print(\" 99th percentile/1-in-100-year event would be $%.2f billion\"%np.percentile(cost,99))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee8c2ee0",
   "metadata": {},
   "source": [
    "You will all get slightly different answers for the above due the randomness of the simulation, but with a large enough value for nsim, I expect you are unlikely to differ by more than ~20 %. \n",
    "\n",
    "- [**Q5**] How reaslistic, do you think the largest events in the simulated series are? If you judge them to be unreaslitic, why?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82e69502",
   "metadata": {},
   "source": [
    "## Wrap up\n",
    "\n",
    "We've covered a lot of ground here. It's now time to take stock. \n",
    "\n",
    "#### What did we do? \n",
    "Under a scenario of responding to a request from the US Treasury, we summarised the *observed* economic impacts from hurricanes making landfall in the USA. We then assessed the statistical distributions of the impact per hurricane (the Log-Normal), and the number of hurricanes per year (Poisson). We used a simple formula to compute the expected hurricane cost per year using paramaters from the the Log-Normal and Poisson distributions. To evaluate the expected magnitude of the 1-in-100-year hurricane season, we then used a \"*stochastic simulation*\" (i.e., the random sampling from the Poisson and Log-Normal distributions). This is a flexible, widely-used approach to explore risk resulting from potentially complex interacting processes.  \n",
    "\n",
    "#### Answers \n",
    "\n",
    "- [Q1] \\$18.9 Billion (2018 adjusted)\n",
    "\n",
    "- [Q2] It is substantially higher (23.06 vs. 9.8 billion \\$ per hurricane)\n",
    "\n",
    "- [Q3] The value for $\\lambda$ is 2.1 storms/year\n",
    "\n",
    "- [Q4] It is substantially higher (48.29 vs. 18.9 billion \\$ per year)\n",
    "\n",
    "- [Q5] Some points to consider: (i) the total value of all assets may be below some of the highest simulated losses; (ii) there may be a physical limit to the number of tropical cyclones per year; (iii) perhaps the damage/hurricane is *not* independent from the total number of hurricanes \n",
    "\n",
    "#### Additional questions\n",
    "[1] Upon reading the report you prepared for the Treasury, the *US Secretary of Homeland Security* concludes: \n",
    "\n",
    "\"*Hurricanes are bad news. We know that and we've dealt with them just fine for over 100 years. We remember the awful 2005 season and already spend plenty on reducing vulnerability, so I see no reason for us to rethink anything now.\"*\n",
    "\n",
    "How would you respond?\n",
    "\n",
    "[2] List the factors that you think would *increase* hurricane risk (as measured by expected financial losses) over the rest of the century\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4f1b9b9",
   "metadata": {},
   "source": [
    "#### Further reading\n",
    "\n",
    "[Blackwell 2015](https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29NH.1527-6996.0000162) \n",
    "\n",
    "[Conte and Kelly 2018](https://miami.pure.elsevier.com/en/publications/an-imperfect-storm-fat-tailed-tropical-cyclone-damages-insurance-) \n",
    "\n",
    "[Flyvbjerg (2020)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7533687/)\n",
    "\n",
    "[Katz 2002](https://journals.ametsoc.org/view/journals/apme/41/7/1520-0450_2002_041_0754_smohd_2.0.co_2.xml) \n",
    "\n",
    "[Malamud (2004)](https://iopscience.iop.org/article/10.1088/2058-7058/17/8/35/meta) \n",
    "\n",
    "[Cirillo and Taleb (2020)](https://www.nature.com/articles/s41567-020-0921-x#Tab1)\n"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}